text:
module SVHN2_Tests
using Test
using ColorTypes
using ImageCore
using FixedPointNumbers
using MLDatasets
using DataDeps
using MAT

@testset "Constants" begin
    @test SVHN2.classnames() isa Vector{Int}
    @test SVHN2.classnames() == [1,2,3,4,5,6,7,8,9,0]
    @test length(SVHN2.classnames()) == 10
    @test length(unique(SVHN2.classnames())) == 10

    @test DataDeps.registry["SVHN2"] isa DataDeps.DataDep
end

@testset "convert2images" begin
    @test_throws DimensionMismatch SVHN2.convert2image(rand(100))
    @test_throws DimensionMismatch SVHN2.convert2image(rand(228,1))
    @test_throws DimensionMismatch SVHN2.convert2image(rand(32,32,4))

    data = rand(N0f8,32,32,3)
    A = SVHN2.convert2image(data)
    @test size(A) == (32,32)
    @test eltype(A) == RGB{N0f8}
    @test SVHN2.convert2image(vec(data)) == A
    @test permutedims(channelview(A), (3, 2, 1)) == data
    @test SVHN2.convert2image(reinterpret(UInt8, data)) == A

    data = rand(N0f8,32,32,3,2)
    A = SVHN2.convert2image(data)
    @test size(A) == (32,32,2)
    @test eltype(A) == RGB{N0f8}
    @test SVHN2.convert2image(vec(data)) == A
    @test SVHN2.convert2image(reinterpret(UInt8, data)) == A
end

# NOT executed on CI. only executed locally.
# This involves dataset download etc.
if parse(Bool, get(ENV, "CI", "false"))
    @info "CI detected: skipping dataset download"
else
    data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
        datadep"SVHN2"
    end

    @testset "Images" begin
        X_train = @inferred SVHN2.traintensor()
        X_test  = @inferred SVHN2.testtensor()
        X_extra = @inferred SVHN2.extratensor()
        @test size(X_train, 4) == 73_257
        @test size(X_test,  4) == 26_032
        @test size(X_extra, 4) == 531_131

        # Sanity check that the first trainimage is not the
        # first testimage nor extra image
        @test X_train[:,:,:,1] != X_test[:,:,:,1]
        @test X_train[:,:,:,1] != X_extra[:,:,:,1]
        @test X_test[:,:,:,1]  != X_extra[:,:,:,1]
        # Make sure other integer types work as indicies
        @test SVHN2.testtensor(0xBAE) == SVHN2.testtensor(2990)

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        for (image_fun, T, nimages) in (
                (SVHN2.testtensor, UInt8,   26_032),
                (SVHN2.testtensor, Int,     26_032),
                (SVHN2.testtensor, Float64, 26_032),
                (SVHN2.testtensor, Float32, 26_032),
                (SVHN2.testtensor, N0f8,    26_032),
            )
            @testset "$image_fun with T=$T" begin
                # whole image set
                A = @inferred image_fun(T)
                @test typeof(A) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                @test size(A) == (32,32,3,nimages)

                @test_throws BoundsError image_fun(T,-1)
                @test_throws BoundsError image_fun(T,0)
                @test_throws BoundsError image_fun(T,nimages+1)

                @testset "load single images" begin
                    # Sample a few random images to compare
                    for i = rand(1:nimages, 3)
                        A_i = @inferred image_fun(T,i)
                        @test typeof(A_i) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                        @test size(A_i) == (32,32,3)
                        @test A_i == A[:,:,:,i]
                    end
                end

                @testset "load multiple images" begin
                    A_5_10 = @inferred image_fun(T,5:10)
                    @test typeof(A_5_10) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test size(A_5_10) == (32,32,3,6)
                    for i = 1:6
                        @test A_5_10[:,:,:,i] == A[:,:,:,i+4]
                    end

                    # also test edge cases `1`, `nimages`
                    indices = [10,3,9,1,nimages]
                    A_vec   = image_fun(T,indices)
                    A_vec_f = image_fun(T,Vector{Int32}(indices))
                    @test typeof(A_vec) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test typeof(A_vec_f) <: Union{Array{T,4},Base.ReinterpretArray{T,4}}
                    @test size(A_vec)   == (32,32,3,5)
                    @test size(A_vec_f) == (32,32,3,5)
                    for i in 1:5
                        @test A_vec[:,:,:,i] == A[:,:,:,indices[i]]
                        @test A_vec[:,:,:,i] == A_vec_f[:,:,:,i]
                    end
                end
            end
        end
    end

    @testset "Labels" begin
        # Sanity check that the first trainlabel is not also
        # the first testlabel
        @test SVHN2.trainlabels(1) != SVHN2.testlabels(1)
        @test SVHN2.trainlabels(1) != SVHN2.extralabels(1)
        @test SVHN2.testlabels(1) != SVHN2.extralabels(1)

        # Check a few hand picked examples. I looked at both the
        # pictures and the native output to make sure these
        # values are correspond to the image at the same index.
        @test SVHN2.trainlabels(1) === 1
        @test SVHN2.trainlabels(2) === 9
        @test SVHN2.trainlabels(1337) === 2
        @test SVHN2.trainlabels(0xCAF) === 3
        @test SVHN2.trainlabels(73_257) === 9
        @test SVHN2.testlabels(1) === 5
        @test SVHN2.testlabels(4) === 10
        @test SVHN2.testlabels(0xDAD) === 4
        @test SVHN2.testlabels(26_032) === 7
        @test SVHN2.extralabels(1) === 4
        @test SVHN2.extralabels(3) === 8
        @test SVHN2.extralabels(531_131) === 4

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        # -- However, technically these tests do not check if
        #    these are the actual SVHN labels of that index!
        for (label_fun, nlabels) in
                    ((SVHN2.trainlabels, 73_257),
                     (SVHN2.testlabels,  26_032),
                     (SVHN2.extralabels, 531_131))
            @testset "$label_fun" begin
                # whole label set
                A = @inferred label_fun()
                @test typeof(A) <: Vector{Int64}
                @test size(A) == (nlabels,)

                @testset "load single label" begin
                    # Sample a few random labels to compare
                    for i = rand(1:nlabels, 10)
                        A_i = @inferred label_fun(i)
                        @test typeof(A_i) <: Int64
                        @test A_i == A[i]
                    end
                end

                @testset "load multiple labels" begin
                    A_5_10 = @inferred label_fun(5:10)
                    @test typeof(A_5_10) <: Vector{Int64}
                    @test size(A_5_10) == (6,)
                    for i = 1:6
                        @test A_5_10[i] == A[i+4]
                    end

                    # also test edge cases `1`, `nlabels`
                    indices = [10,3,9,1,nlabels]
                    A_vec   = @inferred label_fun(indices)
                    A_vec_f = @inferred label_fun(Vector{Int32}(indices))
                    @test typeof(A_vec)   <: Vector{Int64}
                    @test typeof(A_vec_f) <: Vector{Int64}
                    @test size(A_vec)   == (5,)
                    @test size(A_vec_f) == (5,)
                    for i in 1:5
                        @test A_vec[i] == A[indices[i]]
                        @test A_vec[i] == A_vec_f[i]
                    end
                end
            end
        end
    end

    # Check against the already tested tensor and labels functions
    @testset "Data" begin
        for (data_fun, feature_fun, label_fun, nobs) in
                 ((SVHN2.testdata,  SVHN2.testtensor,  SVHN2.testlabels,  26_032),)
            @testset "check $data_fun against $feature_fun and $label_fun" begin
                data, labels = @inferred data_fun()
                @test data == @inferred feature_fun()
                @test labels == @inferred label_fun()

                for i = rand(1:nobs, 10)
                    d_i, l_i = @inferred data_fun(i)
                    @test d_i == @inferred feature_fun(i)
                    @test l_i == @inferred label_fun(i)
                end

                data, labels = @inferred data_fun(5:10)
                @test data == @inferred feature_fun(5:10)
                @test labels == @inferred label_fun(5:10)

                data, labels = @inferred data_fun(Int, 5:10)
                @test data == @inferred feature_fun(Int, 5:10)
                @test labels == @inferred label_fun(5:10)

                indices = [10,3,9,1,nobs]
                data, labels = @inferred data_fun(indices)
                @test data == @inferred feature_fun(indices)
                @test labels == @inferred label_fun(indices)
            end
        end
    end
end

end


text:
struct LennardJones <: Continuous
    energy::Function
    force::Function

    function LennardJones()
        pot(x, y, i, j, d2, u) = u += 4.0 * ((1.0 / d2^6) - (1.0 / d2^3))

        function forces(x, y, i, j, d2, f)
            r = x - y
            dudr = -24.0 * r * ((1.0 / d2^7) - (1.0 / d2^4))
            @inbounds f[i] = f[i] + dudr
            @inbounds f[j] = f[j] - dudr

            return f
        end

        return new(pot, forces)
    end
end

ext: jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    T = Float64
    scalar_sets = (
        MOI.LessThan{T},
        MOI.GreaterThan{T},
        MOI.EqualTo{T},
        MOI.Interval{T},
        MOI.Integer,
        MOI.ZeroOne,
        MOI.Semiinteger{T},
        MOI.Semicontinuous{T},
    )
    scalar_functions = (
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{T},
        MOI.ScalarQuadraticFunction{T},
    )
    vector_sets = (
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.Zeros,
        MOI.Reals,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.PositiveSemidefiniteConeSquare,
        MOI.PositiveSemidefiniteConeTriangle,
    )
    vector_functions = (
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T},
        MOI.VectorQuadraticFunction{T},
    )
    constraints = vcat(
        [(F, S) for F in scalar_functions, S in scalar_sets],
        [(F, S) for F in vector_functions, S in vector_sets],
    )
    MOI.precompile_model(UniversalFallback{Model{T}}, constraints)
    MOI.precompile_model(
        CachingOptimizer{MOI.AbstractOptimizer,UniversalFallback{Model{T}}},
        constraints,
    )
    return
end


text:
"""
Analytic, quadratic Solov'ev equilibrium.
"""
module SolovevSymmetricField

    import ElectromagneticFields.SolovevSymmetric

    export initial_conditions_barely_passing, initial_conditions_barely_trapped,
           initial_conditions_deeply_passing, initial_conditions_deeply_trapped

    export hamiltonian, toroidal_momentum

    SolovevSymmetric.@code(2., 5., 1., 1.) # inject magnetic field code

    initial_conditions_barely_passing() = ([2.5, 0., 0., 3.425E-1], (μ = 1E-2,)) # Δt=2.5, nt=50
    initial_conditions_barely_trapped() = ([2.5, 0., 0., 3.375E-1], (μ = 1E-2,)) # Δt=3.0, nt=100
    initial_conditions_deeply_passing() = ([2.5, 0., 0., 5E-1], (μ = 1E-2,))     # Δt=2.5, nt=25
    initial_conditions_deeply_trapped() = ([2.5, 0., 0., 1E-1], (μ = 1E-2,))     # Δt=5.0, nt=50

    include("guiding_center_4d_common.jl")
    include("guiding_center_4d_equations.jl")

    function toroidal_momentum(t,q)
        R(t,q) * ϑ3(t,q)
    end

    include("guiding_center_4d_diagnostics.jl")

end


text:
#isdefined(:Tcl) || include("../src/Tcl.jl")

module TclBaseTests

using Tcl
using Test

# We need some irrational constants.
const π = MathConstants.π
const φ = MathConstants.φ

@testset "Basic Interface" begin

    @testset "Variables" begin
        i = Tcl.getinterp()
        for (name, value) in (("a", 42), ("1", 1), ("", "empty"),
                              ("π", π), ("w\0rld is beautiful!", true))
            # Check methods.
            Tcl.exec(TclStatus, "array", "unset", name)
            Tcl.setvar(name, value)
            if typeof(value) <: Union{String,Integer}
                @test Tcl.getvar(name) == value
            elseif typeof(value) <: AbstractFloat
                @test Tcl.getvar(name) ≈ value
            end
            @test Tcl.exists(name)
            Tcl.unsetvar(name)
            @test !Tcl.exists(name)

            # Check indexable interface.
            i[name] = value
            if typeof(value) <: Union{String,Integer}
                @test i[name] == value
            elseif typeof(value) <: AbstractFloat
                @test i[name] ≈ value
            end
            @test Tcl.exists(name)
            i[name] = nothing
            @test !Tcl.exists(name)
        end

        for (name1, name2, value) in (("a", "i", 42),
                                      ("1", "2", 12),
                                      ("", "", "really empty"),
                                      ("π", "φ", π),
                                      ("w\0rld is", "beautiful!", true))
            # Check methods.
            Tcl.unsetvar(name1, nocomplain=true)
            Tcl.setvar(name1, name2, value)
            if typeof(value) <: Union{String,Integer}
                @test Tcl.getvar(name1, name2) == value
            elseif typeof(value) <: AbstractFloat
                @test Tcl.getvar(name1, name2) ≈ value
            end
            @test Tcl.exists(name1, name2)
            Tcl.unsetvar(name1, name2)
            @test !Tcl.exists(name1, name2)

            # Check indexable interface.
            i[name1, name2] = value
            if typeof(value) <: Union{String,Integer}
                @test i[name1, name2] == value
            elseif typeof(value) <: AbstractFloat
                @test i[name1, name2] ≈ value
            end
            @test Tcl.exists(name1, name2)
            i[name1, name2] = nothing
            @test !Tcl.exists(name1, name2)
        end
    end

    @testset "Scalars" begin
        var = "x"
        @testset "Integers" begin
            for v in (1, -1, 0, 250, 1<<40)
                x = Tcl.exec(:set, var, v)
                @test x == Tcl.getvar(var)
                @test typeof(x) <: Integer
            end
        end
        @testset "Floats" begin
            for v in (1.0, 0.0, -0.5, π, φ, sqrt(2))
                x = Tcl.exec(:set, var, v)
                @test x == Tcl.getvar(var)
                @test typeof(x) <: Cdouble
            end
        end
        @testset "Strings" begin
            for v in ("", "hellow world!", "1", "true", " ", "\n", "\t", "\r",
                      "\a", "caleçon espiègle")
                x = Tcl.exec(:set, var, v)
                @test x == Tcl.getvar(var)
                @test typeof(x) <: String
            end
            for v in ("\u0", "\u2200", "\u2200 x \u2203 y")
                x = Tcl.exec(:set, var, TclObj(v))
                @test x == Tcl.getvar(var)
                @test typeof(x) <: String
            end
        end
        @testset "Booleans" begin
            for v in (true, false)
                x = Tcl.exec(:set, var, v)
                @test x == Tcl.getvar(var)
                @test typeof(x) <: Integer
            end
        end
    end

    @testset "Lists" begin
        wa = ["", "false", "hello world!", "caleçon espiègle"]
        wb = Tcl.exec(TclObj, "list", wa...)
        wc = Tcl.getvalue(wb)
        @test all(wc .== wa)
        @test eltype(wc) <: String

        xa = [false,true,0,1,2,3]
        xb = Tcl.exec(TclObj, "list", xa...)
        xc = Tcl.getvalue(xb)
        @test all(xc .== xa)
        @test eltype(xc) <: Integer

        ya = [-4,6,-7]
        yb = Tcl.exec(TclObj, "list", ya...)
        yc = Tcl.getvalue(yb)
        @test all(yc .== ya)
        @test eltype(yc) <: Integer

        za = [-1.0, 0.0, 1.0, -4.2, π, sqrt(2), φ]
        zb = Tcl.exec(TclObj, "list", za...)
        zc = Tcl.getvalue(zb)
        @test all(zc .≈ za)
        @test all(zc[1:3] .== za[1:3]) # no loss of precision?
        @test eltype(zc) <: AbstractFloat

        lst1 = Tcl.exec("list",xb,yb)
        @test eltype(lst1) <: Vector{<:Integer}
        @test all(lst1[1] == xa)
        @test all(lst1[2] == ya)

        lst2 = Tcl.exec("list",xb,"hello")
        @test typeof(lst2) == Vector{Any}
        @test all(lst2[1] .== xc)
        @test all(lst2[2] == "hello")

        lst3 = Tcl.exec("list",xb,yb,zb)
        @test typeof(lst3) == Vector{Any}
        @test all(lst3[1] .== xc)
        @test all(lst3[2] .== yc)
        @test all(lst3[3] .== zc)

        lst2b = Tcl.exec(TclObj,"list",xb,"hello")
        lst3b = Tcl.exec(TclObj,"list",xb,yb,zb)
        lst4 = Tcl.exec("list",lst2b,lst3b)
        @test typeof(lst4) == Vector{Any}
        @test length(lst4) == 2
        @test all(lst4[1][1] .== xc)
        @test all(lst4[1][2] == "hello")
        @test all(lst4[2][1] .== xc)
        @test all(lst4[2][2] .== yc)
        @test all(lst4[2][3] .== zc)

        lst5 = Tcl.list(π, 1, "hello", 2:6)
        @test length(lst5) == 4
        @test lst5[1] ≈ π
        @test lst5[2] == 1
        @test lst5[3] == "hello"
        @test all(lst5[4] .== 2:6)
        push!(lst5, sqrt(2))
        @test length(lst5) == 5
        @test lst5[end] ≈ sqrt(2)
        @test all(lst5[4:end][1] .== 2:6)
        @test lst5[0] == nothing
        @test lst5[end+1] == nothing
        A = lst5[[1 2; 3 4; 5 6]]
        @test size(A) == (3, 2)
        @test A[1,1] == lst5[1]
        @test A[1,2] == lst5[2]
        @test A[2,1] == lst5[3]
        @test A[2,2] == lst5[4]
        @test A[3,1] == lst5[5]
        @test A[3,2] == lst5[6]

        #yc = Tcl.exec("list",x,z)
        #r = Tcl.exec(TclObj,"list",4.6,pi)
        #q = Tcl.exec(TclObj,"list",4,pi)
        #Tcl.exec("list",4,pi)
        #Tcl.exec("list",x,z,r)
        #Tcl.exec("list",x,z,"hello",r)
        #i = Tcl.getinterp()
        #i[:v] = r
        #Tcl.eval(i, raw"foreach x $v { puts $x }")

    end
end

end


text:
function readSMAT(filename::AbstractString, indextype::Type, valuetype::Type)
    (rows,header) = readdlm(filename;header=true)
    m = parse(indextype,header[1])
    n = parse(indextype,header[2])
    nz = parse(Int64,header[3])
    assert(size(rows,1) == nz)
    A = sparse(
               convert(Array{indextype,1},rows[:,1])+1,
               convert(Array{indextype,1},rows[:,2])+1,
               convert(Array{valuetype,1},rows[:,3]),
               parse(Int,header[1]),
               parse(Int,header[2])
               )
    return A
end

readSMAT(filename) = readSMAT(filename, Int64, Float64)

# readSMATtoCSC
# readSMATtoCSR
# readSMATtoEdges

function _read_bz2(file::AbstractString)
  cf, p = open(`zcat $(file)`, "r", STDOUT)
  return eachline(cf)
end

function _read_xz(file::AbstractString)
  cf, p = open(`zcat $(file)`, "r", STDOUT)
  return eachline(cf)
end

function _read_gz(file::AbstractString)
  cf, p = open(`zcat $(file)`, "r", STDOUT)
  return eachline(cf)
end

function _read_compressed(filename::AbstractString)
  if     endswith(filename,".gz")
    return _read_gz(filename)
  elseif endswith(filename,".bz2")
    return _read_bz2(filename)
  elseif endswith(filename,".xz")
    return _read_xz(filename)
  else
    return eachline(open(filename, "r"))
  end
end

"""
This function reads a list of edges, which is generally some variant of the
form:
  i j val
or
  i j

"""
function load_edges_list(lineiter; indextype::Type = Int, valuetype::Type = Union{}, header::Bool=false, offset::Int = 1)
  if header
    hdr = collect(Base.Iterators.take(lines, 1))[1] # get the first line
  end
  ei = Vector{indextype}()
  ej = similar(ei)
  ev = Vector{valuetype}()
  for line in lineiter
    parts = split(line)
    push!(ei, parse(indextype, parts[1])+offset)
    push!(ej, parse(indextype, parts[2])+offset)
    if valuetype != Union{}
      push!(ev, parse(valuetype, parts[3]))
    end
  end
  return ei, ej, ev
end

function readSMATtoEdges(filename::AbstractString, indextype::Type, valuetype::Type)
  lines = _read_compressed(filename)
  header = collect(take(lines, 1))[1] # get the first line
  m = parse(indextype,header[1])
  n = parse(indextype,header[2])
  nz = parse(Int64,header[3])
  edges = Vector{Tuple{indextype,indextype,valuetype}}()
  for line in lines
    parts = split(line)
    push!(edges, tuple(
                  parse(indextype, parts[1])+1,
                  parse(indextype, parts[2])+1,
                  parse(valuetype, parts[3])))
  end
  return edges, m, n
end

function readSMATtoEdges(filename::AbstractString, indextype::Type)
  lines = _read_compressed(filename)
  header = collect(take(lines, 1))[1] # get the first line
  m = parse(indextype,header[1])
  n = parse(indextype,header[2])
  nz = parse(Int64,header[3])
  edges = Vector{Tuple{indextype,indextype}}()
  for line in lines
    parts = split(line)
    push!(edges, tuple(
                  parse(indextype, parts[1])+1,
                  parse(indextype, parts[2])+1))
  end
  return edges
end

export readSMAT

#=

readSMATtoEdges
readSMATtoCSR(filename, )
function readSMATtoCSR(filename::AbstractString)

=#

function writeSMAT{T}(filename::AbstractString, A::SparseMatrixCSC{T,Int}; values::Bool=true)
    open(filename, "w") do outfile
        write(outfile, join((size(A,1), size(A,2), nnz(A)), " "), "\n")

        rows = rowvals(A)
        vals = nonzeros(A)
        m, n = size(A)
        for j = 1:n
           for nzi in nzrange(A, j)
              row = rows[nzi]
              val = vals[nzi]
              if values
                write(outfile, join((row-1, j-1, val), " "), "\n")
              else
                write(outfile, join((row-1, j-1, 1), " "), "\n")
              end
           end
        end
    end
end

export writeSMAT


text:
using InferableTasks
using Test

@testset "InferableTasks.jl" begin
    # Write your tests here.
end


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

export
    av_des_alloc,
    av_des_init,
    av_des_crypt,
    av_des_mac

function av_des_alloc()
    ccall((:av_des_alloc,libavutil),Ptr{AVDES},())
end

function av_des_init(d,key,key_bits::Integer,decrypt::Integer)
    ccall((:av_des_init,libavutil),Cint,(Ptr{AVDES},Ptr{UInt8},Cint,Cint),d,key,key_bits,decrypt)
end

function av_des_crypt(d,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_des_crypt,libavutil),Void,(Ptr{AVDES},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{UInt8},Cint),d,dst,src,count,iv,decrypt)
end

function av_des_mac(d,dst,src,count::Integer)
    ccall((:av_des_mac,libavutil),Void,(Ptr{AVDES},Ptr{UInt8},Ptr{UInt8},Cint),d,dst,src,count)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module regularization

export regularize, regularize!, AbstractRegularizer, L1, L2, BoxRegularizer

abstract type AbstractRegularizer end

################################################################################
# L1

struct L1 <: AbstractRegularizer
    gamma::Float64

    function L1(gamma)
        gamma <= 0.0 && error("Gamma is $(gamma). It should be strictly positive.")
        new(gamma)
    end
end

function regularize(reg::L1, B, gradB)
    gamma = reg.gamma
    for (i, x) in enumerate(B)
        gradB[i] = x < 0 ? -gamma : gamma
    end
    gamma * sum(abs(x) for x in B)
end

################################################################################
# L2

struct L2 <: AbstractRegularizer
    lambda::Float64

    function L2(lambda)
        lambda <= 0.0 && error("Lambda is $(lambda). It should be strictly positive.")
        new(lambda)
    end
end

function regularize(reg::L2, B, gradB)
    lambda = reg.lambda
    gradB .= lambda .* B
    0.5 * lambda * sum(x^2 for x in B)
end

################################################################################
# BoxRegularizer

"Constrains each parameter in B to be in [lowerbound, upperbound]."
struct BoxRegularizer <: AbstractRegularizer
    lowerbound::Float64
    upperbound::Float64

    function BoxRegularizer(lb, ub)
        lb >= ub && error("The lower bound must be strictly lower than the upper bound.")
        new(lb, ub)
    end
end

# Called by the solver
function regularize!(outB, inB, reg::BoxRegularizer, gradB)
    lb    = reg.lowerbound
    width = reg.upperbound - lb
    for (i, x) in enumerate(inB)
        pw       = width / (1.0 + exp(-x))  # p*w, where p = 1/(1+exp(-x))
        outB[i]  = lb + pw
        gradB[i] = pw - pw*pw/width  # w * p * (1 - p)
    end
end

#=
   Called after the solver, for returning values in the specified box.
   Same as previous version but excludes updating the gradient.
=#
function regularize!(B, reg::BoxRegularizer)
    lb    = reg.lowerbound
    width = reg.upperbound - lb
    for (i, x) in enumerate(B)
        B[i] = lb + width / (1.0 + exp(-x))
    end
end

end


text:
using StrideArrays
using Documenter

makedocs(;
  modules = [StrideArrays],
  authors = "Chris Elrod",
  repo = "https://github.com/JuliaSIMD/StrideArrays.jl/blob/{commit}{path}#L{line}",
  sitename = "StrideArrays.jl",
  format = Documenter.HTML(;
    prettyurls = get(ENV, "CI", "false") == "true",
    canonical = "https://JuliaSIMD.github.io/StrideArrays.jl",
    assets = String[],
  ),
  pages = [
    "Home" => "index.md",
    "Getting Started" => "getting_started.md",
    "Architecture Benchmarks" =>
      ["arches/cascadelake.md", "arches/tigerlake.md", "arches/haswell.md"],
    "Random Number Generation" => "rng.md",
    "Broadcasting" => "broadcasting.md",
    "Stack Allocattion" => "stack_allocation.md",
  ],
  strict = false,
)

deploydocs(; repo = "github.com/JuliaSIMD/StrideArrays.jl")


text:
include("abstract_tree.jl")
include("serial_tree.jl")
include("parallel_tree.jl")
include("parallel.jl")

# Composite type to hold the actual tree in addition to other mesh-related data
# that is not strictly part of the tree.
mutable struct TreeMesh{TreeType<:AbstractTree}
  tree::TreeType
  current_filename::String
  unsaved_changes::Bool
  first_cell_by_rank::OffsetVector{Int, Vector{Int}}
  n_cells_by_rank::OffsetVector{Int, Vector{Int}}

  function TreeMesh{TreeType}(n_cells_max::Integer) where TreeType
    # Create mesh
    m = new()
    m.tree = TreeType(n_cells_max)
    m.current_filename = ""
    m.unsaved_changes = false
    m.first_cell_by_rank = OffsetVector(Int[], 0)
    m.n_cells_by_rank = OffsetVector(Int[], 0)

    return m
  end

  function TreeMesh{TreeType}(n_cells_max::Integer, domain_center::AbstractArray{Float64},
                              domain_length, periodicity=true) where TreeType
    # Create mesh
    m = new()
    m.tree = TreeType(n_cells_max, domain_center, domain_length, periodicity)
    m.current_filename = ""
    m.unsaved_changes = false
    m.first_cell_by_rank = OffsetVector(Int[], 0)
    m.n_cells_by_rank = OffsetVector(Int[], 0)

    return m
  end
end

const TreeMesh1D = TreeMesh{TreeType} where {TreeType <: AbstractTree{1}}
const TreeMesh2D = TreeMesh{TreeType} where {TreeType <: AbstractTree{2}}
const TreeMesh3D = TreeMesh{TreeType} where {TreeType <: AbstractTree{3}}

# Constructor for passing the dimension and mesh type as an argument
TreeMesh(::Type{TreeType}, args...) where TreeType = TreeMesh{TreeType}(args...)

# Constructor accepting a single number as center (as opposed to an array) for 1D
function TreeMesh{TreeType}(n::Int, center::Real, len::Real, periodicity=true) where {TreeType<:AbstractTree{1}}
  return TreeMesh{TreeType}(n, [convert(Float64, center)], len, periodicity)
end

@inline Base.ndims(mesh::TreeMesh) = ndims(mesh.tree)

# Generate initial mesh
function generate_mesh()
  # Get number of spatial dimensions
  ndims_ = parameter("ndims")

  # Get maximum number of cells that should be supported
  n_cells_max = parameter("n_cells_max")

  # Get domain boundaries
  coordinates_min = parameter("coordinates_min")
  coordinates_max = parameter("coordinates_max")

  # Domain length is calculated as the maximum length in any axis direction
  domain_center = @. (coordinates_min + coordinates_max) / 2
  domain_length = maximum(coordinates_max .- coordinates_min)

  # By default, mesh is periodic in all dimensions
  periodicity = parameter("periodicity", true)

  # Create mesh
  if mpi_isparallel()
    tree_type = ParallelTree{ndims_}
  else
    tree_type = SerialTree{ndims_}
  end
  @timeit timer() "creation" mesh = TreeMesh(tree_type, n_cells_max, domain_center,
                                             domain_length, periodicity)

  # Create initial refinement
  initial_refinement_level = parameter("initial_refinement_level")
  @timeit timer() "initial refinement" for l = 1:initial_refinement_level
    refine!(mesh.tree)
  end

  # Apply refinement patches
  @timeit timer() "refinement patches" for patch in parameter("refinement_patches", [])
    mpi_isparallel() && error("non-uniform meshes not supported in parallel")
    if patch["type"] == "box"
      refine_box!(mesh.tree, patch["coordinates_min"], patch["coordinates_max"])
    else
      error("unknown refinement patch type '$(patch["type"])'")
    end
  end

  # Apply coarsening patches
  @timeit timer() "coarsening patches" for patch in parameter("coarsening_patches", [])
    mpi_isparallel() && error("non-uniform meshes not supported in parallel")
    if patch["type"] == "box"
      coarsen_box!(mesh.tree, patch["coordinates_min"], patch["coordinates_max"])
    else
      error("unknown coarsening patch type '$(patch["type"])'")
    end
  end

  # Partition mesh
  if mpi_isparallel()
    partition!(mesh)
  end

  return mesh
end

# Load existing mesh from file
load_mesh(restart_filename) = load_mesh(restart_filename, mpi_parallel())
function load_mesh(restart_filename, mpi_parallel::Val{false})
  # Get number of spatial dimensions
  ndims_ = parameter("ndims")

  # Get maximum number of cells that should be supported
  n_cells_max = parameter("n_cells_max")

  # Create mesh
  @timeit timer() "creation" mesh = TreeMesh(SerialTree{ndims_}, n_cells_max)

  # Determine mesh filename
  filename = get_restart_mesh_filename(restart_filename, Val(false))
  mesh.current_filename = filename
  mesh.unsaved_changes = false

  # Open mesh file
  h5open(filename, "r") do file
    # Set domain information
    mesh.tree.center_level_0 = read(attrs(file)["center_level_0"])
    mesh.tree.length_level_0 = read(attrs(file)["length_level_0"])
    mesh.tree.periodicity    = Tuple(read(attrs(file)["periodicity"]))

    # Set length
    n_cells = read(attrs(file)["n_cells"])
    resize!(mesh.tree, n_cells)

    # Read in data
    mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
    mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
    mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
    mesh.tree.levels[1:n_cells] = read(file["levels"])
    mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
  end

  return mesh
end

# Obtain the mesh filename from a restart file
get_restart_mesh_filename(restart_filename) = get_restart_mesh_filename(restart_filename, mpi_parallel())
function get_restart_mesh_filename(restart_filename, mpi_parallel::Val{false})
  # Get directory name
  dirname, _ = splitdir(restart_filename)

  # Read mesh filename from restart file
  mesh_file = ""
  h5open(restart_filename, "r") do file
    mesh_file = read(attrs(file)["mesh_file"])
  end

  # Construct and return filename
  return joinpath(dirname, mesh_file)
end


text:
using VisualRegressionTests, Base.Test

cd(@__DIR__)
include("visual_regression_funcs.jl")

reference_image_installed = try
    isa(Pkg.installed("ReferenceImages"), VersionNumber)
catch e
    false
end
if !reference_image_installed
    Pkg.clone("https://github.com/SimonDanisch/ReferenceImages.git")
end

refpath(i) = Pkg.dir("ReferenceImages", "reference_images", "img_$i.png")

function record(filepath, func)
    srand(1234)
    scene = func()
    @assert isa(scene, Scene)
    Makie.save(filepath, scene)
end

function record_reference()
    i = 1
    while true
        # functions in visual_regression_funcs.jl all are named after the same sceme
        func = Symbol("testfunc_$i")
        isdefined(Main, func) || break # we only have n functions in visual_regression_funcs
        func_inst = getfield(Main, func)
        record(refpath(i), func_inst)
        i += 1
    end
end
# record_reference()

@testset "Makie visual regression" begin
    i = 1
    while true
        # functions in visual_regression_funcs.jl all are named after the same sceme
        func = Symbol("testfunc_$i")
        isdefined(Main, func) || break # we only have n functions in visual_regression_funcs
        func_inst = getfield(Main, func)
        path = refpath(i)
        f(fn) = record(fn, func_inst)
        @testset "$func" begin
            @test test_images(VisualTest(f, path), popup = false) |> success
        end
        i += 1
    end
end


text:
using Pkg.TOML

#using TOML

using PeaceVote
using PeaceCypher

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:PeaceVote)

function config(demespec)
    config = Dict()

    #metadata = Dict()
    config["name"] = demespec.name
    config["uuid"] = "$(demespec.uuid)"
    config["maintainer"] = demespec.maintainer ### BigInt could be left unchanged

    ### I could have a second layer for this!

    notary = Dict()
    notary["config"] = "$(demespec.crypto)"
    notary["uuid"] = "$(demespec.cryptodep)"

    cypher = Dict()
    cypher["config"] = "$(demespec.cypherconfig)"
    cypher["uuid"] = "$(demespec.cypherdep)"

    peacefounder = Dict()
    peacefounder["uuid"] = "$(demespec.peacefounder)"

    #config["metadata"] = metadata
    config["notary"] = notary
    config["cypher"] = cypher
    config["peacefounder"] = peacefounder

    return config
end

using Base: UUID

import PeaceVote.DemeSpec
function DemeSpec(dict::Dict)
    name = dict["name"]
    uuid = UUID(dict["uuid"])
    maintainer = dict["maintainer"]
    crypto = Symbol(dict["notary"]["config"])
    cryptodep = UUID(dict["notary"]["uuid"])
    cypherconfig = Symbol(dict["cypher"]["config"])
    cypherdep = UUID(dict["cypher"]["uuid"])
    peacefounder = UUID(dict["peacefounder"]["uuid"])

    DemeSpec(uuid,name,maintainer,crypto,cryptodep,cypherconfig,cypherdep,peacefounder)
end

@show DemeSpec(config(demespec)) == demespec

totoml(demespec) = TOML.print(config(demespec))

#dict = TOML.parse(arrbody)


text:
# Various ways to get part of an Associative Array (indexing, diag, or >,<,==)
using LinearAlgebra

# This is the getindex function for Assoc, the Associate Array.
#import Base.getindex
StringOrNumArray  = Union{AbstractString,Array,Number}

#The Base getindex function which most higher level function would call upon.
function getindex(A::Assoc, i::Array{Int64}, j::Array{Int64})
    #Check if A is empty
    if isempty(A.A)
        return Assoc([1],[1],0,(+))
    end

    # Need to check if numeric values- don't move into Vals, keep in A.A
    if A.val == [1.0]
        return condense(Assoc(A.row[i],A.col[j],A.val,A.A[i,j]))
    else
        return deepCondense(Assoc(A.row[i],A.col[j],A.val,A.A[i,j]))
    end
end

#Singular Case
getindex(A::Assoc,i::Any) = getindex(A,i,:)

PreviousTypes = Array{Int64}

#Variations by derivations
#Each of the additional type would call upon the base functions above.
#It does so by building up row and column combinations of the new type with all previous types.

#Get index with basic addressing.
#Variations between Element, Array, Colon, AbstractRange
getindex(A::Assoc,i::Array{Union{AbstractString,Number}},j::PreviousTypes)                          = getindex(A,findall(x-> x in i,A.row),j)
getindex(A::Assoc,i::PreviousTypes,j::Array{Union{AbstractString,Number}})                          = getindex(A,i,findall(x-> x in j,A.col))
getindex(A::Assoc,i::Array{Union{AbstractString,Number}},j::Array{Union{AbstractString,Number}})    = getindex(A,findall(x-> x in i,A.row),findall(x-> x in j,A.col))

PreviousTypes = Union{PreviousTypes,Array{Union{AbstractString,Number}}}

getindex(A::Assoc,i::Int64,j::PreviousTypes)         = getindex(A,[i],j)
getindex(A::Assoc,i::PreviousTypes,j::Int64)         = getindex(A,i,[j])
getindex(A::Assoc,i::Int64,j::Int64)                 = getindex(A,[i],[j])

PreviousTypes = Union{PreviousTypes,Int64}

getindex(A::Assoc,i::Colon,j::PreviousTypes)         = getindex(A,1:size(A.row,1),j)
getindex(A::Assoc,i::PreviousTypes,j::Colon)         = getindex(A,i,1:size(A.col,1))
getindex(A::Assoc,i::Colon,j::Colon)                 = getindex(A,1:size(A.row,1),1:size(A.col,1))

PreviousTypes = Union{PreviousTypes,Colon}

getindex(A::Assoc,i::AbstractRange,j::PreviousTypes)         = getindex(A,collect(i),j)
getindex(A::Assoc,i::PreviousTypes,j::AbstractRange)         = getindex(A,i,collect(j))
getindex(A::Assoc,i::AbstractRange,j::AbstractRange)                 = getindex(A,collect(i),collect(j))

PreviousTypes = Union{PreviousTypes,AbstractRange}

function convertrange(Akeys,r::AbstractString)
    sep = r[end:end]
    idx = 1
    while occursin(sep*":"*sep,r[idx:end])
        idx = findfirst(sep*":"*sep,r)[1]+1
        from = 1:idx-2
        to = idx+2:length(r)

        if occursin(sep,r[1:idx-2])
            from = findprev(sep,r,idx-2)[1]:idx-2
        end
        if occursin(sep,r[idx+1:end-1])
            to = idx+2:findnext(sep,r,idx+2)[1]-1
        end
        newKeys = join(Akeys[searchsortedfirst(Akeys,r[from]):searchsortedlast(Akeys,r[to])],sep)
        if !isempty(newKeys)
            r = r[1:from[1]-1]*newKeys*r[to[end]+1:end]
        else
            r = r[1:from[end]+1]*r[to[1]:end]
        end
            idx = from[end]+length(newKeys)
    end
    return r
end

#Variations of Single Sequence Strings that is separated by a single character separator.
getindex(A::Assoc, i::AbstractString, j::PreviousTypes)   = getindex(A, findall( x -> in(x,StrUnique(convertrange(A.row,i))[1]),A.row), j)
getindex(A::Assoc, i::PreviousTypes ,j::AbstractString)   = getindex(A, i ,findall( x -> in(x,StrUnique(convertrange(A.col,j))[1]),A.col))
getindex(A::Assoc, i::AbstractString ,j::AbstractString)  = getindex(A, findall( x -> in(x,StrUnique(convertrange(A.row,i))[1]),A.row) ,findall( x -> in(x,StrUnique(convertrange(A.col,j))[1]),A.col))

PreviousTypes = Union{PreviousTypes,AbstractString}

#Variations by Regex
getindex(A::Assoc, i::Regex, j::PreviousTypes)  = getindex(A, findall( x -> occursin(i,x),A.row), j)
getindex(A::Assoc, i::PreviousTypes, j::Regex)  = getindex(A, i, findall( x -> occursin(j,x),A.col))
getindex(A::Assoc, i::Regex, j::Regex)          = getindex(A, findall( x -> occursin(i,x),A.row), findall( x -> occursin(j,x),A.col))

PreviousTypes = Union{PreviousTypes,Regex}

#Variation with StartsWith
struct StartsWith
    inputString::AbstractString
end

# returns an array of indices;
# the corresponding subarray of getrow/getcol is what's called for by startswith
function StartsWithHelper(Ar::Array{Union{AbstractString,Number}},S::StartsWith)
    str_list = []
    if S.inputString[end] == ','
        str_list,tmp = StrUnique(S.inputString)
    else
        str_list = [S.inputString]
    end
    result_indice = Array{Int64,1}()
    for str in str_list
        str_result_first = searchsortedfirst(Ar,str)
        str_result_last = searchsortedlast(Ar,str*string(Char(255)))
        if str_result_first <= str_result_last
            [push!(result_indice,x) for x in str_result_first:str_result_last]
        end
    end
    return result_indice
end

getindex(A::Assoc,i::PreviousTypes,j::StartsWith) = getindex(A,i,StartsWithHelper(getcol(A),j))
getindex(A::Assoc,i::StartsWith,j::PreviousTypes) = getindex(A,StartsWithHelper(getrow(A),i),j)
getindex(A::Assoc,i::StartsWith,j::StartsWith) = getindex(A,StartsWithHelper(getrow(A),i),StartsWithHelper(getcol(A),j))

PreviousTypes = Union{PreviousTypes,StartsWith}

#=
> : get a new Assoc where all of the elements of input Assoc mataches the given Element.
=#
function >(A::Assoc, E::Union{AbstractString,Number})
    if (isa(E,Number) & (A.val ==[1.0])  )
        tarIndex = E
    else
        tarIndex = searchsortedlast(getval(A),E)
    end

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        rowkey, colkey, valkey = findnz(SparseMatrixCSC(A.A))
    else
        rowkey, colkey, valkey = findnz(A.A)
    end
    mapping = findall( x-> x > tarIndex, valkey)
    rows, cols, vals = find(A)

    outA = Assoc(rows[mapping],cols[mapping],vals[mapping])

    if A.val==[1.0]
        outA = putVal(outA,A.val)
    end

    return outA
end

>(E::Union{AbstractString,Number},A::Assoc) = (A < E)

#=
< : get a new Assoc where all of the elements of input Assoc mataches the given Element.
=#
function <(A::Assoc, E::Union{AbstractString,Number})
    if (isa(E,Number) & (A.val ==[1.0])  )
        tarIndex = E
    else
        tarIndex = searchsortedfirst(A.val,E)
    end

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        rowkey, colkey, valkey = findnz(SparseMatrixCSC(A.A))
    else
        rowkey, colkey, valkey = findnz(A.A)
    end
    mapping = findall( x-> x < tarIndex, valkey)
    rows, cols, vals = find(A)

    outA = Assoc(rows[mapping],cols[mapping],vals[mapping])

    if A.val==[1.0]
        outA = putVal(outA,A.val)
    end

    return outA
end

<(E::Union{AbstractString,Number},A::Assoc) = (A > E)

#=
== : get a new Assoc where all of the elements of input Assoc matches the given Element.
=#
(==)(A::Assoc,E::Union{AbstractString,Number}) = equal(A::Assoc,E::Union{AbstractString,Number})
function equal(A::Assoc, E::Union{AbstractString,Number})
    if (isa(E,Number) && (A.val == [1.0])  )
        tarIndex = E
    else
        tarIndex = searchsortedfirst(A.val,E)
        if !(E == getval(A)[tarIndex])
            tarIndex = 0
        end
    end

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        rowkey, colkey, valkey = findnz(SparseMatrixCSC(A.A))
    else
        rowkey, colkey, valkey = findnz(A.A)
    end
    mapping = findall( x-> x == tarIndex, valkey)
    rows,cols,vals = find(A)

    Aout = Assoc(rows[mapping],cols[mapping],vals[mapping])

    # May not need this anymore
    if A.val==[1.0]
        Aout = putVal(Aout,A.val)
    end

    return Aout
end

==(E::Union{AbstractString,Number},A::Assoc) = (A == E)

function bounded(A::Assoc, E1::Union{AbstractString,Number}, E2::Union{AbstractString,Number})
    if (isa(E1,Number) & (A.val ==[1.0])  )
        tarIndex1 = E1
    else
        tarIndex1 = searchsortedfirst(A.val,E1)
    end

    if (isa(E2,Number) & (A.val ==[1.0])  )
        tarIndex2 = E2
    else
        tarIndex2 = searchsortedfirst(A.val,E2)
    end

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        rowkey, colkey, valkey = findnz(SparseMatrixCSC(A.A))
    else
        rowkey, colkey, valkey = findnz(A.A)
    end
    mapping = findall( x-> tarIndex1 <= x <= tarIndex2, valkey)
    rows, cols, vals = find(A)

    outA = Assoc(rows[mapping],cols[mapping],vals[mapping])

    if A.val==[1.0]
        outA = putVal(outA,A.val)
    end

    return outA
end

function strictbounded(A::Assoc, E1::Union{AbstractString,Number}, E2::Union{AbstractString,Number})
    if (isa(E1,Number) & (A.val ==[1.0])  )
        tarIndex1 = E1
    else
        tarIndex1 = searchsortedfirst(A.val,E1)
    end

    if (isa(E2,Number) & (A.val ==[1.0])  )
        tarIndex2 = E2
    else
        tarIndex2 = searchsortedfirst(A.val,E2)
    end

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        rowkey, colkey, valkey = findnz(SparseMatrixCSC(A.A))
    else
        rowkey, colkey, valkey = findnz(A.A)
    end
    mapping = findall( x-> tarIndex1 < x < tarIndex2, valkey)
    rows, cols, vals = find(A)

    outA = Assoc(rows[mapping],cols[mapping],vals[mapping])

    if A.val==[1.0]
        outA = putVal(outA,A.val)
    end

    return outA
end

#=
diag : Output the diagonal of input Assoc A.
Outputs the Assoc with only the diagonal elements of A.
=#
function diag(A::Assoc)
    # Check if numeric values first
    if A.val == [1.0]
        return Assoc(A.row,A.col,A.val,dropzeros!(sparse(diagm(0 => diag(A.A)))))
    else
        return deepCondense(Assoc(A.row,A.col,A.val,sparse(diagm(0 => diag(A.A)))))
    end
end

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


text:
using Revise
using NLSolvers
using LinearAlgebra
using SparseDiffTools
using SparseArrays
using IterativeSolvers
using ForwardDiff

# Start with a stupid example
n = 10
A = sprand(10, 10, 0.1)
A = A*A' + I
F(Fx, x) = mul!(Fx, A, x)

x0, = rand(10)
xp = copy(x0)
Fx = copy(xp)

function theta(x)
   if x[1] > 0
       return atan(x[2] / x[1]) / (2.0 * pi)
   else
       return (pi + atan(x[2] / x[1])) / (2.0 * pi)
   end
end

function F_powell!(Fx, x)
        if ( x[1]^2 + x[2]^2 == 0 )
            dtdx1 = 0;
            dtdx2 = 0;
        else
            dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
            dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
        end
        Fx[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
        Fx[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
        Fx[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
    Fx
end

function F_jacobian_powell!(Fx, Jx, x)
    ForwardDiff.jacobian!(Jx, F_powell!, Fx, x)
    Fx, Jx
end
Fc, Jc = zeros(3), zeros(3,3)
F_jacobian_powell!(Fc, Jc, x0)

import NLSolvers: OnceDiffedJv
function OnceDiffedJv(F; seed, autodiff=false)
    JacOp = JacVec(F, seed; autodiff=false)
    OnceDiffedJv(F, JacOp)
end

jv = JacVec(F_powell!, rand(3); autodiff=false)
function jvop(x)
    jv.u .= x
    jv
end
prob_obj = NLSolvers.VectorObjective(F_powell!, nothing, F_jacobian_powell!, jvop)

prob = NEqProblem(prob_obj)
x0 = [-1.0, 0.0, 0.0]
res = solve(prob, copy(x0), LineSearch(Newton(), Backtracking()))
@test norm(res.info.best_residual) < 1e-15

using IterativeSolvers
function inexact_linsolve(x0, JvOp, Fx, ηₖ)
    krylov_iter = IterativeSolvers.gmres_iterable!(x0, JvOp, Fx; maxiter=50)
    res = copy(Fx)
    rhs = ηₖ*norm(Fx, 2)
    for item in krylov_iter
        res = krylov_iter.residual.current
        if res <= rhs
            break
        end
    end
    return x0, res
end
res = solve(prob, copy(x0), InexactNewton(inexact_linsolve, FixedForceTerm(0.001), 1e-4, 300), NEqOptions(maxiter=1000))
@test norm(res.info.best_residual) < 1e-10

function f_2by2!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
    F
end

function g_2by2!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
    J
end

FJacOp2 = OnceDiffedJv(f_2by2!; seed=rand(2))
solve(FJacOp2, rand(2), ResidualKrylov(FixedForceTerm(0.4), 1e-3, 300))

function f_sparse!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function g_sparse!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

function helical_valley()
    tpi = 8*atan(1)
    c7 = 2.5e-1
    c8 = 5e-1

    function f!(fvec, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        fvec[1] = 10(x[3] - 10*temp1)
        fvec[2] = 10(temp2 - 1)
        fvec[3] = x[3]
    end
    function fj!(F, J, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        if !isa(F, Nothing)
            F[1] = 10*(x[3] - 10*temp1)
            F[2] = 10*(temp2 - 1)
            F[3] = x[3]
        end
        if !isa(J, Nothing)
            temp = x[1]^2 + x[2]^2
            temp1 = tpi*temp
            temp2 = sqrt(temp)
            J[1,1] = 100*x[2]/temp1
            J[1,2] = -100*x[1]/temp1
            J[1,3] = 10
            J[2,1] = 10*x[1]/temp2
            J[2,2] = 10*x[2]/temp2
            J[2,3] = 0
            J[3,1] = 0
            J[3,2] = 0
            J[3,3] = 1
            if !isa(F, nothing)
                return F, J
            else
                return J
            end
        end

    end
    f!, fj!, rand(3), [-1.; 0; 0]), [-1.; 0; 0], "Helical Valley")
end


text:
"""
    CopyTo <: CellRule

    CopyTo{W}(from)
    CopyTo{W}(; from)

A simple rule that copies aux array slices to a grid over time.
This can be used for comparing simulation dynamics to aux data dynamics.
"""
struct CopyTo{W,F} <: CellRule{Tuple{},W}
    "An Aux or Grid key for data source or a single value"
    from::F
end
CopyTo(from) = CopyTo{:_default_}(from)
CopyTo(; from) = CopyTo{:_default_}(from)
CopyTo{W}(from) where W = CopyTo{W,typeof(from)}(from)
CopyTo{W}(; from) where W = CopyTo{W,typeof(from)}(from)

ConstructionBase.constructorof(::Type{<:CopyTo{W}}) where W = CopyTo{W}

DynamicGrids.applyrule(data, rule::CopyTo, state, I) = get(data, rule.from, I...)
DynamicGrids.applyrule(data, rule::CopyTo{W}, state, I) where W <: Tuple =
    ntuple(i -> get(data, rule.from, I...), length(_asiterable(W)))


text:
# include -
include("Include.jl")

# How many parameters?
number_of_parameters = 153;
number_of_initial_sets = 5;
number_of_final_sets = 200;

# Simple script to generate ensemble by sampling -
tmp_parameter_array = zeros(number_of_parameters,number_of_initial_sets);
for parameter_set_index = 1:number_of_initial_sets

  # Setup tmp filename -
  tmp_filename = "./parameter_best_v1.dat."*string(parameter_set_index);
  parameter_set = readdlm(tmp_filename);

  for parameter_index = 1:number_of_parameters
    tmp_parameter_array[parameter_index,parameter_set_index] = parameter_set[parameter_index];
  end
end

# Establish bounds -
bounds_array = zeros(number_of_parameters,2);
for parameter_index = 1:number_of_parameters

  # min and max -
  minimum_value = minimum(tmp_parameter_array[parameter_index,:]);
  maximum_value = maximum(tmp_parameter_array[parameter_index,:]);

  # perturbation -
  perturbation_size = 1;
  bounds_array[parameter_index,1] = (1/perturbation_size)*minimum_value;
  bounds_array[parameter_index,2] = perturbation_size*maximum_value;
end

# Generate an ensemble of parameter sets -
ensemble_array = zeros(number_of_parameters,number_of_final_sets);
sobol_sequence = SobolSeq(number_of_parameters,bounds_array[:,1],bounds_array[:,2]);
for parameter_set_index = 1:number_of_final_sets

  # sample the parameter array -
  parameter_array = next(sobol_sequence)

  for parameter_index = 1:number_of_parameters
    ensemble_array[parameter_index,parameter_set_index] = parameter_array[parameter_index];
  end
end

# Solve the model to calculate the error -
error_array = zeros(number_of_final_sets);
for parameter_set_index = 1:number_of_final_sets

  # sample the ensemble_array -
  parameter_array = ensemble_array[:,parameter_set_index]

  # error -
  error_value = sum(objective_function(parameter_array));
  error_array[parameter_set_index] = error_value;

  @show (parameter_set_index,error_value)
end

# Filter -
minimum_error = minimum(error_array);
variance_error = var(error_array);
epsilon = 0.5;
error_cutoff = minimum_error*(1.0+epsilon);
ensemble_index_array = Int[];
for parameter_set_index = 1:number_of_final_sets

  error_value = error_array[parameter_set_index];

  if (error_value<=error_cutoff)
    push!(ensemble_index_array,parameter_set_index);
  else

    # can we accept this set?
    error_difference = (error_value - minimum_error)^2
    prob_keep_set = exp(-1*error_difference/variance_error);
    check_role = rand();
    if (prob_keep_set>check_role)
      push!(ensemble_index_array,parameter_set_index);
    end
  end
end

# Last ... run the model w/these sets, save the data -
for (local_index,parameter_set_index) in enumerate(ensemble_index_array)

  # sample the ensemble_array -
  parameter_array = ensemble_array[:,parameter_set_index]

  # Load the data dictionary -
  data_dictionary = DataDictionary(0.0,0.0,0.0);

  # Update data dictionary to match new parameters before calculating obj
  parameter_mapping_array = data_dictionary["parameter_name_mapping_array"]
  for index = 1:length(parameter_mapping_array)
    if index <= data_dictionary["number_of_binding"]
      data_dictionary["binding_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    else
      data_dictionary["control_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    end
  end

  # Setup the time -
  time_start = 0.0
  time_stop = 120.0
  time_step_size = 0.01

  # Run the simulation -
  (T,X) = Simulation(time_start,time_stop,time_step_size,data_dictionary);

  # dump data to disk -
  local_data = [T X];
  data_filename = "./ensemble/sim_data_EI_"*string(parameter_set_index)*"_LI_"*string(local_index)*".dat"
  writedlm(data_filename,local_data);
end

# Write the simple ensemble -
filename = "./simple_ensemble_parameter_sets.dat"
writedlm(filename,ensemble_array);


text:
function sparse_from_coo(::Type{<:SparseMatrixCSC}, I,J,V,m,n)
  sparse(I,J,V,m,n)
end

@inline function push_coo!(
   ::Type{<:SparseMatrixCSC},I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number)

   push!(I, ik)
   push!(J, jk)
   push!(V, vk)
   nothing
end

@inline function is_entry_stored(::Type{<:SparseMatrixCSC},i::Integer,j::Integer)
  true
end

function push_coo!(I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number)
  push_coo!(SparseMatrixCSC{Float64,Int},I,J,V,ik,jk,vk)
end

function finalize_coo!(
  ::Type{<:SparseMatrixCSC},I::Vector,J::Vector,V::Vector,m::Integer,n::Integer)
  nothing
end

function finalize_coo!(I::Vector,J::Vector,V::Vector,m::Integer,n::Integer)
    finalize_coo!(SparseMatrixCSC,I,J,V,m,n)
end

function add_entry!(A::SparseMatrixCSC{Tv,Ti},v::Number,i::Integer,j::Integer,combine::Function=+) where {Tv,Ti<:Integer}
    cv = convert(Tv, v)
    ci = convert(Ti, i)
    cj = convert(Ti, j)
    m,n = size(A)
    if !((1 <= ci <= m) & (1 <= cj <= n))
        throw(BoundsError(A, (ci,cj)))
    end
    cols = getptr(A)
    indices = getindices(A)
    vals = nonzeros(A)
    coljfirstk = Int(cols[cj])
    coljlastk = Int(cols[cj+1] - 1)
    searchk = searchsortedfirst(indices, ci, coljfirstk, coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && indices[searchk] == ci
        # Column j contains entry A[i,j]. Update
        vold = vals[searchk]
        vals[searchk] = combine(vold, cv)
    else
        throw(ArgumentError("Entry ($i,$j) does not exist in the sparsity pattern"))
    end
    return A
end

function copy_entries!(a::SparseMatrixCSC,b::SparseMatrixCSC)
  _copy_entries_sparse!(a,b)
end

function _copy_entries_sparse!(a,b)
  na = nonzeros(a)
  nb = nonzeros(b)
  if na !== nb
    copyto!(na,nb)
  end
end

hasrowmajororder(::Type{<:SparseMatrixCSC}) = false
hasrowmajororder(a::SparseMatrixCSC) = hasrowmajororder(SparseMatrixCSC)

hascolmajororder(::Type{<:SparseMatrixCSC}) = true
hascolmajororder(a::SparseMatrixCSC) = hascolmajororder(SparseMatrixCSC)

getptr(S::SparseMatrixCSC) = S.colptr

getindices(S::SparseMatrixCSC) = rowvals(S)


text:
using Plots
using LaTeXStrings
using DFWannier
pyplot()
plot_font=font(15,"DejaVu Sans")

pyplot(lab="",yguidefont=plot_font,ytickfont=plot_font,xtickfont=plot_font,legendfont=plot_font)
T=Float32
# x = create_TB_model("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/","/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",[[PhysAtom(0.0,0.0,-0.0239129,-0.155854) for i=1:4]...,[PhysAtom(0.0,0.0,5.5540692,0.318205) for i=1:4]...],Float64)

# x = WannierModel{T}("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/","/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",[[PhysAtom(T[0.0,0.0,-0.0239129,-0.155854]...) for i=1:4]...,[PhysAtom(T[0.0,0.0,5.5540692,0.318205]...) for i=1:4]...]);
x = WannierModel{T}("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/test1","/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",[[PhysAtom(T[0.0,0.0,-0.0239129,-0.155854]...) for i=1:4]...,[PhysAtom(T[0.0,0.0,5.5540692,0.318205]...) for i=1:4]...]);
calculate_angmom(x.wfcs[6],x.wfcs[8])
points = [WfcPoint3D(-p.w,p.p) for p in x.wfcs[8].points]
x.wfcs[8] = Wfc3D(points,x.wfcs[8].cell,x.wfcs[8].atom)
x.wfcs[1].cell
recip = T[1.184212  0.000000  0.403316;-0.592106  1.025558  0.403316;-0.592106 -1.025558  0.403316]
kz = 0.
kz = 0.6049734
ky_max = 1.0
kx_max = 1.0
kx_min,ky_min,kz_min = recip'*x.k_points[90]
kx_max,ky_max,kz_max = recip'*x.k_points[110]
kx_max = 1.184212/1
ky_max = 1.0255575/1
# kx_max = 0.0399672
# ky_max = ky_min
kp_cart = [T[kx,ky,kz] for kx=linspace(-ky_max,ky_max,30),ky=linspace(-ky_max,ky_max,30)]
kp_rec = [inv(recip')*kc for kc in kp_cart]
kp_rec = reshape(kp_rec,(900,1))
k_input = reshape(kp_rec,900,1)
k_plot = reshape(kp_cart,(length(kp_cart)))
tbbands=calculate_eig_cm_angmom(x,k_input);
tbbandssoc=calculate_eig_cm_angmom_soc(x,k_input);

begin
kxs = [k[1] for k in k_plot]
kys = [k[2] for k in k_plot]
# kzs = [k[3] for k in k_input]
u = Float64[]
v = Float64[]
c = Float64[]
# w = Float[]
for (i,(spin,angmom)) in enumerate(zip(tbbandssoc[10].spins,tbbandssoc[10].angmoms))
  sx,sy,sz = Array(spin[2])
  # sx,sy,sz = Array(spin[2])+Array(spin[1])
  lx,ly,lz = Array(angmom[2])
  # lx,ly,lz = Array(angmom[2])+Array(angmom[1])
  push!(u,(lx)/120)
  push!(v,(ly)/120)
  push!(c,sz)
end
println(maximum(c))
quiver(kxs,kys,quiver=(u,v),color=c)
end

6nd
2lot(tbbandssoc[6],:angmom2_x)
test_points = [[ka,kb,kc] for (ka,kb,kc) in zip(linspace(0.5,0.25,201), linspace(1.0,1.0,201),linspace(0.0,0.25,201))]
# dfbands = read_qe_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/bands.out",T)
dfbands = read_qe_bands_file("/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",T)
# dfbandssoc = read_qe_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",T);
dfbandsnsoc = read_qe_bands_file("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/rest/bands.out",T)
# potential=read_potential_file("/home/ponet/Documents/PhD/GeTe/NSOC/p_1.0/density.out")
tbbandssocf=calculate_eig_soc(x)
tbbandsf = calculate_eig_angmom(x)
tbbands=calculate_eig_cm_angmom(x,90:0.2:110);
tbbands=calculate_eig_cm_angmom(x,k_input);
tbbandssoc = calculate_eig_cm_angmom_soc(x,90:0.2:110);
tbbandssoc = calculate_eig_cm_angmom_soc(x);
tbbandssoc = calculate_eig_cm_angmom_soc(x,k_input);
pyplot()
plot(plot(tbbandssoc[9],:angmom2_x),plot(tbbands[5],:angmom2_x),plot(tbbandssoc[9],:angmom2_y),plot(tbbands[5],:angmom2_y))
plot(plot([[real(eig[6]) for eig in tbbands[4].eigvec] [real(eig[8]) for eig in tbbands[4].eigvec] [real(eig[7]) for eig in tbbands[4].eigvec]],title="real no soc"),plot([[imag(eig[6]) for eig in tbbands[4].eigvec] [imag(eig[8]) for eig in tbbands[4].eigvec] [imag(eig[7]) for eig in tbbands[4].eigvec]],title="imag no soc"),plot([[real(eig[6]) for eig in tbbandssoc[7].eigvec] [real(eig[8]) for eig in tbbandssoc[7].eigvec] [real(eig[7]) for eig in tbbandssoc[7].eigvec]],title="real soc"),plot([[imag(eig[6]) for eig in tbbandssoc[7].eigvec] [imag(eig[8]) for eig in tbbandssoc[7].eigvec] [imag(eig[7]) for eig in tbbandssoc[7].eigvec]],title="imag soc"),labels=["pz" "py" "pz"])
#----- paper plots
# two small and 2 big
plot(plot([tbbandssoc[9:10]...,tbbands[5]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot(plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_x,label = ["SOC" "SOC" "No SOC"],leg=false,title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),layout=(2,1)), plot(plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),layout=(2,1)),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"))
plot(plot([tbbandssoc[5:6]...,tbbands[3]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot(plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_x,label = ["SOC" "SOC" "No SOC"],leg=false,title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),layout=(2,1)), plot(plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),layout=(2,1)),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"))
#all same size
plot(plot([tbbandssoc[9:10]...,tbbands[5]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(2.5:0.2:4)),plot([tbbandssoc[9:10]...,tbbands[5]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_x,label = ["SOC" "SOC" "No SOC"],leg=false,title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-1:0.2:0.2)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title="",yticks=collect(-0.2:0.2:0.6)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/oamvseigvalsv1.pdf")
plot(plot([tbbandssoc[5:6]...,tbbands[3]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(3.4:0.2:5)),plot([tbbandssoc[5:6]...,tbbands[3]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)",yticks=collect(-3:0.2:-1)), plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_x,leg=false,label = ["SOC" "SOC" "No SOC"],title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-0.5:0.2:1)),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/oamvseigvalsv3.pdf")
plot(plot([tbbandssoc[7:8]...,tbbands[4]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot([tbbandssoc[7:8]...,tbbands[4]],:angmom2_x,leg=false,label = ["SOC" "SOC" "No SOC"],title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[7:8]...,tbbands[4]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),plot([tbbandssoc[7:8]...,tbbands[4]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/oamvseigvalsv2.pdf")

#-----------plot 3 bandstructures-------------9.2879
plot(plot(dfbands[11:19],tbbandsf,title="No SOC",label="",fermi=9.2879,xticks=([0,101,201],["A","Z","U"]),subplot_index=1,ylim=[-13,9]),plot(dfbands[21:37],tbbandssoc,title="SOC",yl="",fermi=9.2879,xticks=([0,101,201],["A","Z","U"]),subplot_index=2,ylim=[-13,9],leg=false),plot(dfbands[21:37],tbbandssoc,title="SOC",yl="",fermi=9.2879,xticks=([0,101,201],["A","Z","U"]),subplot_index=3,ylim=[-5,5]),titlefont = font(30,"Courier"),size=(1800,600),layout=(1,3))
plot(plot(dfbands[11:19],tbbandsf,title="No SOC",label="",fermi=9.2879,subplot_index=1,ylim=[-13,9]),plot(dfbandssoc[21:37],tbbandssoc,title="SOC",yl="",fermi=9.2879,subplot_index=2,ylim=[-13,9],leg=false),plot(dfbandssoc[21:37],tbbandssoc,title="SOC",yl="",fermi=9.2879,subplot_index=3,ylim=[-5,5]),titlefont = font(30,"Courier"),size=(1800,600),layout=(1,3))
#-----------plot heatmaps-----------------
begin
  T=Float64
  dfbandssoc = read_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",T);
  zmat_te,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/pwo.pdos_atm#1(Te)_wfc#2(p)",fermi=9.2879);
  zmat_ge,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/pwo.pdos_atm#2(Ge)_wfc#2(p)",fermi=9.2879);
  # zmat_ge,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/pwo.pdos_atm#2(Ge)_wfc#2(p)",fermi=9.2879)
  k_lim = [-norm(dfbandssoc[1].k_points[101]-dfbandssoc[1].k_points[1]),0.0,norm(dfbandssoc[1].k_points[101]-dfbandssoc[1].k_points[end])]
end
plot(plot(dfbandssoc,fermi=9.2879,color=:red,xticks=(k_lim,["A","Z","U"]),ylims=(-5-9.2879,18.01-9.2879),width=1,label="",title=""),heatmap(zmat_ge,legend=false,color=:BuPu,yticks=(yticks[1][5:5:end]+10,["","","",""]),xticks=([0,101,201],["A","Z","U"])),heatmap(zmat_te,color=:BuPu,xticks=([0,101,201],["A","Z","U"]),yticks=(yticks[1][5:5:end]+10,["","","",""])),layout=(1,3),size=(1200,500),right_margin=5mm)
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/paperxsf/rest/heatmap_bs.png")
#-----
recip = [1.184212  0.000000  0.403316;
        -0.592106  1.025558  0.403316;
        -0.592106 -1.025558  0.403316]
k_cart = []
for k in x.k_points
  push!(k_cart,recip'*k)
end
x.k_points[end]
k_cart[end]
#------Calculation of H_eff----------------------------------------------------------------#
function Hd(c,Ez,band)
  recip = [1.184212  0.000000  0.403316;
          -0.592106  1.025558  0.403316;
          -0.592106 -1.025558  0.403316]
  out = []
  for i in eachindex(band.k_points)
    push!(out,c*Ez*(recip'*band.k_points[i])[2]*(band.spins[i][2].x+band.angmoms[i][2].x))
  end
  return out
end

function Hso(l_so,band)
  out = []
  for i in eachindex(band.k_points)
    push!(out,l_so*(band.angmoms[i][2].x*band.spins[i][2].x))
  end
  return out
end

H0 = tbbands[5].eigvals

H10=H0+Hd(0.1,-0.5,tbbandssoc[10])+Hso(0.29,tbbandssoc[10])
H9=H0+Hd(0.1,-0.5,tbbandssoc[9])+Hso(0.29,tbbandssoc[9])
H5=tbbands[3].eigvals+Hd(9.0,-0.5,tbbandssoc[5])+Hso(0.29,tbbandssoc[5])
H6=tbbands[3].eigvals+Hd(9.0,-0.5,tbbandssoc[6])+Hso(0.29,tbbandssoc[6])
plot([H10,H9,tbbandssoc[6].eigvals,tbbandssoc[5].eigvals],line=[:blue :red])
plot(H0)
#---------------Proof of J rather than L------------#
\lambda
tbbandssoc[10].angmoms[50][2]
tbbandssoc[9].angmoms[50][2]
tbbandssoc[6].angmoms[50][2]
tbbandssoc[5].angmoms[50][2]
tbbandssoc[10].cms[49]-tbbandssoc[10].cms[50]
tbbandssoc[9].cms[49]-tbbandssoc[9].cms[50]
tbbandssoc[6].cms[49]-tbbandssoc[6].cms[50]
tbbandssoc[5].cms[49]-tbbandssoc[5].cms[50]
tbbandssoc[10].angmoms[50][2]+tbbandssoc[10].spins[50][2]
tbbandssoc[9].angmoms[50][2]+tbbandssoc[9].spins[50][2]
tbbandssoc[6].angmoms[50][2]+tbbandssoc[6].spins[50][2]
tbbandssoc[5].angmoms[50][2]+tbbandssoc[5].spins[50][2]

 #-------------CM stuff-------------------
Ge_s_b = construct_bloch_sum(x.wfcs[1],[0.5,0.5,0.5])
Ge_pz_b = construct_bloch_sum(x.wfcs[2],[0.5,0.5,0.5])
Ge_px_b = construct_bloch_sum(x.wfcs[3],[0.5,0.5,0.5])
Ge_py_b = construct_bloch_sum(x.wfcs[4],[0.5,0.5,0.5])
Te_s_b = construct_bloch_sum(x.wfcs[5],[0.5,0.5,0.5])
Te_pz_b = construct_bloch_sum(x.wfcs[6],[0.5,0.5,0.5])
Te_px_b = construct_bloch_sum(x.wfcs[7],[0.5,0.5,0.5])
Te_py_b = construct_bloch_sum(x.wfcs[8],[0.5,0.5,0.5])
Ge_s_b = construct_bloch_sum(x.wfcs[1],[0.,0.,0.])
Ge_pz_b = construct_bloch_sum(x.wfcs[2],[0.4,0.55,0.5])
Ge_px_b = construct_bloch_sum(x.wfcs[3],[0.4,0.55,0.5])
Ge_py_b = construct_bloch_sum(x.wfcs[4],[0.4,0.55,0.5])
Te_s_b = construct_bloch_sum(x.wfcs[5],[0.4,0.55,0.5])
Te_pz_b = construct_bloch_sum(x.wfcs[6],[0.4,0.55,0.5])
Te_px_b = construct_bloch_sum(x.wfcs[7],[0.4,0.55,0.5])
Te_py_b = construct_bloch_sum(x.wfcs[8],[0.4,0.55,0.5])
write_xsf_file("/Users/ponet/Downloads/test.xsf",Ge_s_b)
Ge_s = x.wfcs[1]
Ge_pz = x.wfcs[2]
Ge_px = x.wfcs[3]
Ge_py = x.wfcs[4]
Te_s = x.wfcs[5]
Te_pz = x.wfcs[6]
Te_px = x.wfcs[7]
Te_py = x.wfcs[8]
calculate_cm(Te_pz,Te_pz)

#---------------Eq 11-------------
term1,term2,wfc = eq11(x,9,98)
write_dipole_mesh("/Users/ponet/Downloads/term2test.xsf",term2,:z)
write_xsf_file("/Users/ponet/Downloads/wfctest.xsf",wfc)
#--------------different polarizations------

x = create_TB_model("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.5/","/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",[[PhysAtom(0.0,0.0,-0.0239129,-0.155854) for i=1:4]...,[PhysAtom(0.0,0.0,5.5540692,0.318205) for i=1:4]...],Float64);
tbbandssocf,tmp=calculate_eig_SOC(x);
tbbands,tmp=calculate_eig_cm(x);
tbbands,dip1=calculate_eig_cm_angmom(x,90:0.2:110);
tbbandssoc,dip= calculate_eig_cm_angmom_SOC(x,90:0.2:110);
dfbands = read_qe_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/bands.out",Float64);
plotlyjs()
dfbandssoc = read_qe_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out");
plot(plot([tbbandssoc[9:10]...,tbbands[5]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(0:0.2:10)),plot([tbbandssoc[9:10]...,tbbands[5]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)"), plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_x,label = ["SOC" "SOC" "No SOC"],leg=false,title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-1:0.2:1)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[9:10]...,tbbands[5]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title="",yticks=collect(-1:0.2:1)),plot([tbbandssoc[9:10]...,tbbands[5]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.5/oamvseigvalsv1.pdf")
plot(plot([tbbandssoc[5:6]...,tbbands[3]],:cm_z,label=["SOC" "SOC" "No SOC"],leg=false,xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(0:0.2:10)),plot([tbbandssoc[5:6]...,tbbands[3]],:eigvals,fermi=9.2879,label=["SOC" "SOC" "No SOC"],xticks=([-0.05,0.0,0.05],["","",""]),yguide=L"E-E_f"*" (eV)",yticks=collect(-3:0.2:-1)), plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_x,leg=false,label = ["SOC" "SOC" "No SOC"],title="OAM around Te",xticks=([-0.05,0.0,0.05],["","",""]),yticks=collect(-1:0.2:1)),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_x,fermi=9.2879,leg=false,title="SAM around Te",xticks=([-0.05,0.0,0.05],["","",""])),plot([tbbandssoc[5:6]...,tbbands[3]],:angmom2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),xguide=L"|$\mathbf{k}_r$|",title=""),plot([tbbandssoc[5:6]...,tbbands[3]],:spin2_y,fermi=9.2879,leg=false,xticks=([-0.05,0.0,0.05],[0.05,"Z",0.05]),title="",xguide=L"|$\mathbf{k}_r$|"),size=(1200,1200),guidefont=font(20,"DejaVu Sans"),titlefont=font(20,"DejaVu Sans"),layout=(3,2))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.5/oamvseigvalsv3.pdf")
plot(dfbandssoc[21:37],tbbandssocf,title="SOC",yl="",fermi=9.2879,subplot_index=2,ylim=[-13,9],leg=false)
dfbands = read_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/bands.out")
plot(tbbandssocf,legend=false,ylim=[2,10])
pyplot()
begin
  T=Float64
  dfbandssoc = read_bands_file("/Users/ponet/Documents/Fysica/PhD/GeTe/fullrel/GeTe_bands.out",T);
  zmat_te,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/pwo.pdos_atm#1(Te)_wfc#2(p)",fermi=9.2879);
  zmat_ge,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/pwo.pdos_atm#2(Ge)_wfc#2(p)",fermi=9.2879);
  # zmat_ge,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/pwo.pdos_atm#2(Ge)_wfc#2(p)",fermi=9.2879)
  k_lim = [-norm(dfbandssoc[1].k_points[101]-dfbandssoc[1].k_points[1]),0.0,norm(dfbandssoc[1].k_points[101]-dfbandssoc[1].k_points[end])]
end
zmat_ges,yticks=read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/pwo.pdos_atm#2(Ge)_wfc#1(s)",fermi=9.2879);
zmat_tes, = read_KPDOS("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/pwo.pdos_atm#1(Te)_wfc#1(s)",fermi=9.2879);
plot(heatmap(zmat_ge),heatmap(zmat_te),heatmap(zmat_ges),heatmap(zmat_tes))
savefig("/Users/ponet/Documents/Fysica/PhD/GeTe/colin/p_0.0/heatmap_p.png")

#----------------total angmom-----------------#
#notes: Weird things happen with the soc calculation. Nonetheless the non-soc calculation seems to work out fine, and results
#       in a much larger Lx than is seen around the atoms. The same dynamics seem to still apply though. The issue with the sign
#       of Ly still remains
x = create_TB_model("/home/ponet/Documents/PhD/GeTe/NSOC/paperxsf/","/home/ponet/Documents/PhD/GeTe/SOC/GeTe_bands.out",[[PhysAtom(0.0,0.0,-0.0239129,-0.155854) for i=1:4]...,[PhysAtom(0.0,0.0,5.5540692,0.318205) for i=1:4]...],Float64);

tbbands=calculate_eig_cm_angmom(x,90:0.2:110);
tbbandssoc= calculate_eig_cm_angmom_SOC(x,90:0.2:110);

total_angmom = []
total_angmom_soc = []
for i = 1:length(tbbands[5].eigvec)
  vec = tbbands[5].eigvec[i]
  vec_soc = tbbandssoc[10].eigvec[i]
  wfc = sum([v*w for (v,w) in zip(vec,x.wfcs)])
  wfc_soc = sum([v*w for (v,w) in zip(vec_soc,[x.wfcs;x.wfcs])])
  center = calculate_cm(wfc,wfc)
  center_soc = calculate_cm(wfc_soc,wfc_soc)
  push!(total_angmom,calculate_angular_momentum(wfc,wfc,center))
  push!(total_angmom_soc,calculate_angular_momentum(wfc_soc,wfc_soc,center_soc))
end
total_angmom[1][4]

plot(map(x->real(x[1]),total_angmom))

gete = load_server_job("GeTe_2/nonrel",phd_dir*"GeTe/NSOC")
remove_flags!(GeTe,:smearing)
remove_flags!(GeTe,Symbol("starting_magnetization(1)"))
change_flags!(GeTe,:occupations=>"'fixed'",:degauss => 0.0f0)
print_flow(GeTe)
change_flow!(GeTe,"scf.in"=>true)
print_flags(gete,"wan")
pull_outputs(GeTe,extras=["*.xsf","*r.dat"])
submit_job(GeTe)
replace_header_word!(GeTe,"frontend","defpart")

#pol calc

te_ez = -3.65563
volume = 360.3090*0.52917724900001^3
ez_rel = te_ez/volume
p_te = ez*((0.52325284+0.002252856)-0.5)*3*6.686168
# = -0.03502832319972468

add_calculation!(GeTe,deepcopy(get_input(GeTe,"GeTe_nscf.in")),filename="berry1.in")
add_flags!(GeTe,"berry1.in",:control, :lberry => true, :gdir => 1 , :nppstr => 10)
change_flags!(GeTe,"berry1.in",:gdir=>1)
add_calculation!(GeTe,deepcopy(get_input(GeTe,"berry1.in")),filename="berry2.in")
change_flags!(GeTe,"berry2.in",:gdir=>2)
add_calculation!(GeTe,deepcopy(get_input(GeTe,"berry2.in")),filename="berry3.in")
change_flags!(GeTe,"berry3.in",:gdir=>3)
change_data_option!(GeTe,"berry3.in",:k_points,:automatic)
change_flow!(GeTe,"scf"=>false)
GeTe.name = "GeTe_NSOC"


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

## basic task functions and TLS

# Container for a captured exception and its backtrace. Can be serialized.
struct CapturedException <: Exception
    ex::Any
    processed_bt::Vector{Any}

    function CapturedException(ex, bt_raw::Vector)
        # bt_raw MUST be a vector that can be processed by StackTraces.stacktrace
        # Typically the result of a catch_backtrace()

        # Process bt_raw so that it can be safely serialized
        bt_lines = process_backtrace(bt_raw, 100) # Limiting this to 100 lines.
        CapturedException(ex, bt_lines)
    end

    CapturedException(ex, processed_bt::Vector{Any}) = new(ex, processed_bt)
end

function showerror(io::IO, ce::CapturedException)
    showerror(io, ce.ex, ce.processed_bt, backtrace=true)
end

"""
    CompositeException

Wrap a `Vector` of exceptions thrown by a [`Task`](@ref) (e.g. generated from a remote worker over a channel
or an asynchronously executing local I/O write or a remote worker under `pmap`) with information about the series of exceptions.
For example, if a group of workers are executing several tasks, and multiple workers fail, the resulting `CompositeException` will
contain a "bundle" of information from each worker indicating where and why the exception(s) occurred.
"""
struct CompositeException <: Exception
    exceptions::Vector{Any}
    CompositeException() = new(Any[])
    CompositeException(exceptions) = new(exceptions)
end
length(c::CompositeException) = length(c.exceptions)
push!(c::CompositeException, ex) = push!(c.exceptions, ex)
isempty(c::CompositeException) = isempty(c.exceptions)
iterate(c::CompositeException, state...) = iterate(c.exceptions, state...)
eltype(::Type{CompositeException}) = Any

function showerror(io::IO, ex::CompositeException)
    if !isempty(ex)
        showerror(io, ex.exceptions[1])
        remaining = length(ex) - 1
        if remaining > 0
            print(io, string("\n\n...and ", remaining, " more exception(s).\n"))
        end
    else
        print(io, "CompositeException()\n")
    end
end

function show(io::IO, t::Task)
    print(io, "Task ($(t.state)) @0x$(string(convert(UInt, pointer_from_objref(t)), base = 16, pad = Sys.WORD_SIZE>>2))")
end

"""
    @task

Wrap an expression in a [`Task`](@ref) without executing it, and return the [`Task`](@ref). This only
creates a task, and does not run it.

# Examples
```jldoctest
julia> a1() = sum(i for i in 1:1000);

julia> b = @task a1();

julia> istaskstarted(b)
false

julia> schedule(b);

julia> yield();

julia> istaskdone(b)
true
```
"""
macro task(ex)
    :(Task(()->$(esc(ex))))
end

"""
    current_task()

Get the currently running [`Task`](@ref).
"""
current_task() = ccall(:jl_get_current_task, Ref{Task}, ())

"""
    istaskdone(t::Task) -> Bool

Determine whether a task has exited.

# Examples
```jldoctest
julia> a2() = sum(i for i in 1:1000);

julia> b = Task(a2);

julia> istaskdone(b)
false

julia> schedule(b);

julia> yield();

julia> istaskdone(b)
true
```
"""
istaskdone(t::Task) = ((t.state == :done) | istaskfailed(t))

"""
    istaskstarted(t::Task) -> Bool

Determine whether a task has started executing.

# Examples
```jldoctest
julia> a3() = sum(i for i in 1:1000);

julia> b = Task(a3);

julia> istaskstarted(b)
false
```
"""
istaskstarted(t::Task) = ccall(:jl_is_task_started, Cint, (Any,), t) != 0

istaskfailed(t::Task) = (t.state == :failed)

task_result(t::Task) = t.result

task_local_storage() = get_task_tls(current_task())
function get_task_tls(t::Task)
    if t.storage === nothing
        t.storage = IdDict()
    end
    (t.storage)::IdDict{Any,Any}
end

"""
    task_local_storage(key)

Look up the value of a key in the current task's task-local storage.
"""
task_local_storage(key) = task_local_storage()[key]

"""
    task_local_storage(key, value)

Assign a value to a key in the current task's task-local storage.
"""
task_local_storage(key, val) = (task_local_storage()[key] = val)

"""
    task_local_storage(body, key, value)

Call the function `body` with a modified task-local storage, in which `value` is assigned to
`key`; the previous value of `key`, or lack thereof, is restored afterwards. Useful
for emulating dynamic scoping.
"""
function task_local_storage(body::Function, key, val)
    tls = task_local_storage()
    hadkey = haskey(tls,key)
    old = get(tls,key,nothing)
    tls[key] = val
    try body()
    finally
        hadkey ? (tls[key] = old) : delete!(tls,key)
    end
end

# NOTE: you can only wait for scheduled tasks
function wait(t::Task)
    if !istaskdone(t)
        if t.donenotify === nothing
            t.donenotify = Condition()
        end
    end
    while !istaskdone(t)
        wait(t.donenotify)
    end
    if istaskfailed(t)
        throw(t.exception)
    end
end

"""
    fetch(t::Task)

Wait for a Task to finish, then return its result value. If the task fails with an
exception, the exception is propagated (re-thrown in the task that called fetch).
"""
function fetch(t::Task)
    wait(t)
    task_result(t)
end

## lexically-scoped waiting for multiple items

function sync_end(refs)
    c_ex = CompositeException()
    for r in refs
        try
            wait(r)
        catch
            if !isa(r, Task) || (isa(r, Task) && !istaskfailed(r))
                rethrow()
            end
        finally
            if isa(r, Task) && istaskfailed(r)
                push!(c_ex, CapturedException(task_result(r), r.backtrace))
            end
        end
    end

    if !isempty(c_ex)
        throw(c_ex)
    end
    nothing
end

const sync_varname = gensym(:sync)

"""
    @sync

Wait until all lexically-enclosed uses of `@async`, `@spawn`, `@spawnat` and `@distributed`
are complete. All exceptions thrown by enclosed async operations are collected and thrown as
a `CompositeException`.
"""
macro sync(block)
    var = esc(sync_varname)
    quote
        let $var = Any[]
            v = $(esc(block))
            sync_end($var)
            v
        end
    end
end

# schedule an expression to run asynchronously

"""
    @async

Wrap an expression in a [`Task`](@ref) and add it to the local machine's scheduler queue.
"""
macro async(expr)
    thunk = esc(:(()->($expr)))
    var = esc(sync_varname)
    quote
        local task = Task($thunk)
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
    end
end

function register_taskdone_hook(t::Task, hook)
    tls = get_task_tls(t)
    push!(get!(tls, :TASKDONE_HOOKS, []), hook)
    t
end

# runtime system hook called when a task finishes
function task_done_hook(t::Task)
    # `finish_task` sets `sigatomic` before entering this function
    err = istaskfailed(t)
    result = task_result(t)
    handled = false
    if err
        t.backtrace = catch_backtrace()
    end

    if isa(t.donenotify, Condition) && !isempty(t.donenotify.waitq)
        handled = true
        notify(t.donenotify, result, true, err)
    end

    # Execute any other hooks registered in the TLS
    if isa(t.storage, IdDict) && haskey(t.storage, :TASKDONE_HOOKS)
        foreach(hook -> hook(t), t.storage[:TASKDONE_HOOKS])
        delete!(t.storage, :TASKDONE_HOOKS)
        handled = true
    end

    if err && !handled
        if isa(result,InterruptException) && isdefined(Base,:active_repl_backend) &&
            active_repl_backend.backend_task.state == :runnable && isempty(Workqueue) &&
            active_repl_backend.in_eval
            throwto(active_repl_backend.backend_task, result) # this terminates the task
        end
    end
    # Clear sigatomic before waiting
    sigatomic_end()
    try
        wait() # this will not return
    catch e
        # If an InterruptException happens while blocked in the event loop, try handing
        # the exception to the REPL task since the current task is done.
        # issue #19467
        if isa(e,InterruptException) && isdefined(Base,:active_repl_backend) &&
            active_repl_backend.backend_task.state == :runnable && isempty(Workqueue) &&
            active_repl_backend.in_eval
            throwto(active_repl_backend.backend_task, e)
        else
            rethrow()
        end
    end
end

"""
    timedwait(testcb::Function, secs::Float64; pollint::Float64=0.1)

Waits until `testcb` returns `true` or for `secs` seconds, whichever is earlier.
`testcb` is polled every `pollint` seconds.
"""
function timedwait(testcb::Function, secs::Float64; pollint::Float64=0.1)
    pollint > 0 || throw(ArgumentError("cannot set pollint to $pollint seconds"))
    start = time()
    done = Channel(1)
    timercb(aw) = begin
        try
            if testcb()
                put!(done, :ok)
            elseif (time() - start) > secs
                put!(done, :timed_out)
            end
        catch e
            put!(done, :error)
        finally
            isready(done) && close(aw)
        end
    end

    if !testcb()
        t = Timer(timercb, pollint, interval = pollint)
        ret = fetch(done)
        close(t)
    else
        ret = :ok
    end
    ret
end

max_stars_repo_name: syntapy/julia

max_issues_repo_name: syntapy/julia

max_forks_repo_name: syntapy/julia


text:
using Pluto
using UUIDs

# TODO create a copy with a uuid of src/normal.jl and load it
newfile = "src/" * string(UUIDs.uuid4()) * ".jl"
tmppath, tmpio = mktemp()
open("src/normal.jl") do io
    for line in eachline(io, keep = true)
        write(tmpio, line)
    end
end
close(tmpio)
mv(tmppath, newfile, force = true)

Pluto.run(notebook = newfile,
          require_secret_for_access = false, require_secret_for_open_links = false,
          port = parse(Int, ARGS[1]), host = "0.0.0.0",
          show_file_system = false, launch_browser = false)

max_stars_repo_name: roualdes/Normal

max_issues_repo_name: roualdes/Normal

max_forks_repo_name: roualdes/Normal


text:
include("default_params.jl")
include("reflecting_model/model_def.jl")

include("../helper_functions.jl")

# initial condition
function map_to_t(d0,B,approx_type,t,h=0.005,limit=false,kwargs...)
    if (approx_type==:dg)&&(limit)
        dt = integrate_time(d0,B,t,StableRK4(h); limiter=GeneralisedMUSCL)
    elseif (approx_type==:dg)&&(!limit)
        dt = integrate_time(d0,B,t,StableRK4(h); limiter=NoLimiter)
    else
        dt = integrate_time(d0,B,t,StableRK4(h))
    end
    return dt
end

point_mass_initial_condition(n,rng) = (X = fill(eps(),n), φ = fill(1,n))
d0_map_point_mass(dq) = interior_point_mass(eps(),1,dq)

function exp_initial_condition(n,rng)
    i = rand(rng,1:2,n)
    X = Vector{Float64}(undef,n)
    for k in 1:n
        x = Inf
        while x>=model.b
            x = -log(rand(rng))
        end
        X[k]=x
    end
    return (X=X,φ=i)
end
d0_map_exp(dq) = SFMDistribution((x,i)->(i∈(1:2))*exp(-x)/(1-exp(-model.b))/2,dq)

function map_to_t(d0,B,approx_type,t,h=0.005,limit=false,kwargs...)
    if (approx_type==:limiter)
        dt = integrate_time(d0,B,t,StableRK4(h); limiter=GeneralisedMUSCL)
    end
    return dt
end

function make_approx_data()
    pth = mkpath((@__DIR__)*"/reflecting_model/data/exp/data")
    mkpath(pth*"/coeffs")
    mkpath(pth*"/cdfs")

    h0 = 1
    for k in 1:2:21
        h = h0/((k+1)/2)
        o = 2
        nodes = 0:h:10

        mesh = DGMesh(nodes,o)

        dq = DiscretisedFluidQueue(model,mesh)

        d0 = d0_map_exp(dq)

        B = build_full_generator(dq)

        dt = map_to_t(d0,B.B,:limiter,transient_time)

        CSV.write(
            pth*"/coeffs/order_"*string(k)*"_model_limiter.csv",
            DataFrame(a=dt.coeffs),
        )
        x = range(0,model.b,length=n_err_evals)
        phase_1_evaluated = cdf(dt).(x,1)
        phase_2_evaluated = cdf(dt).(x,2)
        CSV.write(
            pth*"/cdfs/order_"*string(k)*"_model_limiter.csv",
            DataFrame(x=x,
                phase_1 = phase_1_evaluated,
                phase_2 = phase_2_evaluated,
            ),
        )
    end
    return nothing
end

make_approx_data()


text:
using SeasonalStreamflowForecasts

using Test
using Dates: Date

@testset "SeasonalStreamflowForecasts" begin
    ssf = SSF()
    @testset "SSF()" begin
        nrows, ncols = size(ssf.sites)
        @test nrows > 0 # At least one site.
        @test ncols > 0 # At least one column.
    end

    @testset "get_forecasts()" begin
        site_ids = ["410730"]
        fc_date = Date(2020,8,1)
        for awrc_id in site_ids
            data, header = get_forecasts(ssf, awrc_id, fc_date)

            local nrows, ncols = size(data)
            @test nrows > 0 # At least one data point.
            @test ncols > 0 # At least one column.

            @test length(header) > 0 # At least one header line
        end
        data, header = get_forecasts(ssf, "invalid ID", fc_date)
        @test isempty(data)
        @test isempty(header)
    end
    @testset "show()" begin
        show_str = repr(ssf)
        @test occursin("ID", show_str)
        @test occursin("AWRC", show_str)
    end
    @testset "close()" begin
        close!(ssf)
        @test isempty(ssf.sites)
    end
end


text:
push!(LOAD_PATH, "../src/")

using Documenter, POMDPModelTools

makedocs(
    modules = [POMDPModelTools],
    format = Documenter.HTML(),
    sitename = "POMDPModelTools.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPModelTools.jl.git",
)


text:
@testset "1015.smallest-integer-divisible-by-k.jl" begin
    @test smallest_repunit_div_by_k(787) == 393
    @test smallest_repunit_div_by_k(77) == 6
    @test smallest_repunit_div_by_k(7) == 6
    @test smallest_repunit_div_by_k(541) == 540
    @test smallest_repunit_div_by_k(547) == 91
end


text:
using FunctionalCollections

export Component,
       Action,
       noop,
       setfield,
       childaction,
       indexaction

"""
The abstract type for all Components
"""
abstract Component

"""
The view of a component
returns a UI with annotated intents
"""
function view
end

convert(::Type{Tile}, c::Component) = view(c)

"""
An action by/of a component - can be obtained from
the intent set in the view
"""
abstract Action <: Intent

"""
Given a `Component` and an `Action` transition to
a new `Component` after applying in the given `Action`
"""
function update
end

immutable NoOp <: Action
end

"""
No-op action
"""
const noop = NoOp()
update(c::Component, a::NoOp) = c

"""
Set a field. Update generates code for each component type and field.
Component must not have inner constructors for this to work
"""
immutable SetField{field, T} <: Action
    value::T
end
@compat (::Type{SetField{f}}){f,T}(x::T) = SetField{f, T}(x)

setfield(f::Symbol) = SetField{f}
setfield(f, x) = SetField{f}(x)

@generated function update{field}(c::Component, x::SetField{field})
    Expr(:call, c, [field === f ? :(x.value) : :(c.$f) for f in fieldnames(c)]...)
end

"""
Child action on a field
"""
immutable ChildAction{A<:Action} <: Action
    filed::Symbol
    action::A
end

childaction(f, action) = ChildAction(f, action)
childaction(f, idx, action) = childaction(f, indexaction(idx, action))

function update(c::Component, a::ChildAction)
    update(c, setfield(a.field, update(getfield(c, a.field), a.action)))
end

"""
Action on an index in a collection (vector / dict)
"""
immutable IndexAction{I,A<:Action} <: Action
    index::I
    action::A
end

function update(c::Union{Associative, AbstractArray}, a::IndexAction)
    assoc_(c, a.index, update(c[a.index], a))
end

assoc_(c::Union{PersistentVector, PersistentHashMap}, idx, x) = assoc(c, idx, x)
function assoc_(c, idx, x)
    c′ = copy(c)
    c′[idx] = x
    c′
end


text:
export intersperse

# Utility functions for Elem

boolattr(a, name) = a ? name : nothing

@doc """
Intersperse a value in between elements in a vector

Optionally you can tell it to enclose the result in the seperator element.

e.g.
```
intersperse(0, [1, 2, 3])
# => [1, 0, 2, 0, 3]

intersperse(0, [1, 2, 3], true)
# => [0, 1, 0, 2, 0, 3, 0]
```
""" ->
function intersperse(x, xs, enclose=false)
    if length(xs) > 1
        res = foldl((acc, nxt) -> vcat(acc, x, nxt),
                    Any[xs[1]], xs[2:end])
    else
        res = xs
    end
    enclose ? [x, res, x] : res
end

getproperty(el::Elem, prop, default) =
    hasproperties(el) ? get(properties(el), prop, default) : default

wrapmany(t::TileList, wrap) =
    length(t.tiles) == 1 ?
        render(t.tiles[1]) :
        render(t.tiles, wrap)

@api class => Class <: Tile begin
    arg(class::String)
    curry(content::TileList)
    kwarg(wrap::Symbol=:div)
    kwarg(forcewrap::Bool=false)
end

maybestring(s::String) = s
maybestring(s::TileList) =
    length(s.tiles) == 1 ? maybestring(s.tiles[1]) : render(s.tiles)
maybestring(s) = render(s)

addclasses(t, cs) =
    t & [:className => cs * " " * getproperty(t, :className, "")]

render(c::Class) =
    addclasses(c.forcewrap ? Elem(c.wrap, maybestring(c.content)) :
                           wrapmany(c.content, c.wrap),
               c.class)


text:
#
# HyperCubic lattice
#

# An N-dimensional hypercubic lattice with periodic
# boundary conditions
struct HyperCubic{N}
  dims::NTuple{N,Int}
end
Base.size(l::HyperCubic) = l.dims

const Chain = HyperCubic{1}
const Square = HyperCubic{2}
const Cubic = HyperCubic{3}

struct Edge{N}
  edge::Tuple{NTuple{N,Int},NTuple{N,Int}}
  boundary::Bool
end

Base.getindex(edge::Edge, n::Int) = edge.edge[n]

Base.reverse(edge::Edge) = Edge(reverse(edge.edge), edge.boundary)

sites(lattice::HyperCubic) = (Tuple(s) for s in CartesianIndices(axes(lattice)))

function is_in_edge(site::Tuple, edge::Edge)
  if site == edge[1]
    return false
  elseif site == edge[2]
    return true
  else
    error("Site $site is not incident to edge $edge")
  end
end

function onehot_tuple(n::Integer, length::Val{N}) where {N}
  return ntuple(i -> i == n ? 1 : 0, Val(N))
end

# Obtain the neighbor in dimension `dim` in direction `dir`, for example:
# neighbor((2, 2), 2, 1) == (2, 3)
# neighbor((2, 2), 2, -1) == (2, 1)
# neighbor((2, 2), 1, -1) == (1, 2)
# neighbor((2, 3), 2, 1, lattice_size=(3, 3)) == (2, 1)
function neighbor(
  site::NTuple{N,Int}, dim::Int, dir::Int; lattice_size=typemax(eltype(site))
) where {N}
  return map(
    (a, b, d) -> mod1(a + dir * b, d), site, onehot_tuple(dim, Val(N)), lattice_size
  )
end

# Check if the edge connecting to the specified neighbor of the site
# crosses the boundary of the lattice
function isboundary(
  site::NTuple{N,Int}, dim::Int, dir::Int; lattice_size=typemax(eltype(site))
) where {N}
  return ((dir == 1) && (site[dim] == lattice_size[dim])) || ((dir == -1) && site[dim] == 1)
end

# The neighboring sites of the specified site
function filterneighbors(
  f, lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false
) where {N}
  lattice_size = size(lattice)
  site_neighbors = Vector{NTuple{N,Int}}()
  for dim in 1:N, dir in (-1, 1)
    site_neighbor = neighbor(site, dim, dir; lattice_size=lattice_size)
    bc_condition = periodic || !(isboundary(site, dim, dir; lattice_size=lattice_size))
    if f(site, site_neighbor) && bc_condition
      push!(site_neighbors, site_neighbor)
    end
  end
  return site_neighbors
end

function neighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
  return filterneighbors(≠, lattice, site; periodic=periodic)
end

function inneighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
  return filterneighbors(>, lattice, site; periodic=periodic)
end

function outneighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
  return filterneighbors(<, lattice, site; periodic=periodic)
end

# All of the edges connected to the vertex `site`
function incident_edges(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
  lattice_size = size(lattice)
  site_edges = Vector{Edge{N}}()
  for dim in 1:N, dir in (-1, 1)
    site_neighbor = neighbor(site, dim, dir; lattice_size=lattice_size)
    boundary = false
    if isboundary(site, dim, dir; lattice_size=lattice_size)
      boundary = true
    end
    edge = (site, site_neighbor)
    if dir == -1
      edge = reverse(edge)
    end
    push!(site_edges, Edge(edge, boundary))
  end
  return site_edges
end

function bonds(lattice::HyperCubic; periodic=false)
  return [
    (s, n) for s in sites(lattice) for n in outneighbors(lattice, s; periodic=periodic)
  ]
end

function bonds(lattice::Square, coord::Tuple{Colon,<:Integer})
  rowsize = lattice.dims[1]
  colsites = [(i, coord[2]) for i in 1:(rowsize - 1)]
  return [(s, (s[1] + 1, s[2])) for s in colsites]
end

function bonds(lattice::Square, coord::Tuple{<:Integer,Colon})
  colsize = lattice.dims[2]
  rowsites = [(coord[1], i) for i in 1:(colsize - 1)]
  return [(s, (s[1], s[2] + 1)) for s in rowsites]
end