text:
using Plots
plotlyjs()

function total_time(dic; names = dic["names"], labels = dic["names"])
    plt = Plots.plot(title = "Seconds Distribution")
    for i in 1:length(names)
        name = names[i]
        Plots.plot!(sort(dic["$(name)_tot"]), label=labels[i], linewidth=3)
    end
    display(plt)
    return plt
end

function total_time_relative(dic; names = dic["names"], labels=names, cols=distinguishable_colors(length(names)))

    if isa(cols,Dict)
        cols = [cols[x] for x in names];
    end

  plt = Plots.plot(title = "(Time / Time for $(labels[1])), Distribution")
  rel = dic["$(names[1])_tot"]

  mn = 1
  mx = 1

    for i in 2:length(names)
        name = names[i]
      y = sort(min.(dic["$(name)_tot"] ./ rel, 10))
      mn = min(mn,minimum(y))
      mx = max(mx,maximum(y))
      Plots.plot!(y, label=labels[i], linewidth=3, color=cols[i])
  end
  Plots.plot!(ylimits=(0:10),yticks=[0;round(mn,2);1:10])

  plot!(ytickfont = font(16))
  plot!(legendfont = font(16))
  plot!(guidefont = font(16))
    plot!(xticks=[])

  display(plt)
  return plt
end

function total_time_nnz(dic; names = dic["names"],
                        labels = names,
                        cols=distinguishable_colors(length(names)),
                        ref=[])
    if isa(cols,Dict)
        cols = [cols[x] for x in names];
    end

    if ~isempty(ref)
        yref = 10*dic["$(ref)_tot"] ./ dic["ne"]
    end

    plt = Plots.plot()
    mn = Inf
    mx = 0
    for name in names
        y = dic["$(name)_tot"] ./ dic["ne"]
        if ~isempty(ref)
            y = min.(y, yref)
        end

        mn = min(mn,minimum(y))
        mx = max(mx,maximum(y[y .< Inf]))
    end
    b = floor(log(mx)/log(10))
    for i in 1:length(names)
        name = names[i]
               y = dic["$(name)_tot"] ./ dic["ne"]
        if ~isempty(ref)
            y = min.(y, yref)
        end
      y = sort(y)
      y[y .== Inf] = mx
      y = y * 10^(-b)
      Plots.plot!(y, label=labels[i], linewidth=3, color=cols[i])
    end
  plot!(ytickfont = font(12))
  plot!(legendfont = font(16))
  plot!(guidefont = font(16))
    plot!(xticks=[])
  title!("$(10.0^(round.(Int,b))) Seconds / Edge, Distribution")
  display(plt)
  return mn, mx
end

function time_vs_n(dic; names = dic["names"], labels=names, cols=distinguishable_colors(length(names)))
    if isa(cols,Dict)
        cols = [cols[x] for x in names];
    end

    plt = Plots.plot()
    for i in 1:length(names)
        name = names[i]
      Plots.plot!(dic["nv"],dic["$(name)_tot"], label=labels[i], linewidth=3, color=cols[i])
  end
    Plots.plot!(xaxis=(:log), yaxis=(:log))
    Plots.plot!(xlabel="Number Vertices",ylabel="Seconds")
      plot!(ytickfont = font(12))
      plot!(xtickfont = font(12))
  plot!(legendfont = font(16))
  plot!(guidefont = font(16))
  display(plt)
  return plt
end

"""
    Given the list of solvers that were used, and a dic containing the results.
    Plot some graphs explaining how it went.
"""
function plots_and_stats(solvers,dic)

  plots_and_stats_raw(solvers,dic,"tot","Total Time")
  plots_and_stats_raw(solvers,dic,"build","Build Time")
  plots_and_stats_raw(solvers,dic,"solve","Solve Time")
  plots_and_stats_raw(solvers,dic,"its","Iterations")

end

function plots_and_stats(dic)

  plots_and_stats_raw(dic,"tot","Total Time")
  plots_and_stats_raw(dic,"build","Build Time")
  plots_and_stats_raw(dic,"solve","Solve Time")
  plots_and_stats_raw(dic,"its","Iterations")

end

function plot_relative(dic)

  plot_relative(dic,"tot","Total Time")
  plot_relative(dic,"build","Build Time")
  plot_relative(dic,"solve","Solve Time")
  plot_relative(dic,"its","Iterations")

end

function vec_stats(v)
    v = sort(v)
    n = length(v)
    print("mean: ",round(mean(v),3))
    pts = round.(Int, [1, n/4, n/2, 3n/4, n])
    print(" quarts: ", v[pts])
    println()
end

function vec_stats(v,cnt)
    v = sort(v)
    n = length(v)
    print("mean: ",round.(mean(v),3))
    pts = round.(Int, [1, n/4, n/2, 3n/4, n])
    print(" quarts: ", v[pts])

    print("  (Inf: ",sum(isinf.(cnt)), ")")

    println()
end

function plots_and_stats_raw(solvers,dic,str,title)

  plt = Plots.plot(title = title)
  sol1 = solvers[1]
  p = sortperm(dic["$(sol1.name)_$(str)"])
  for solver in solvers
      Plots.plot!(dic["$(solver.name)_$(str)"][p], label=solver.name)
  end
  display(plt)
end

function plots_and_stats_raw(dic,str,title)

  plt = Plots.plot(title = title)
    sol1 = dic["names"][1]
    key = "$(sol1)_$(str)"
    p = sortperm(dic[key])

    baserow = dic["$(sol1)_$(str)"][p]

    for name in dic["names"]
        key = "$(name)_$(str)"
        thisrow = dic[key][p]
        threshrow = min.(thisrow, 10*baserow)
        Plots.plot!(threshrow, label=name)

        print(key, " ")
        vec_stats(threshrow, thisrow)
  end
  display(plt)
end

function plot_relative(dic,str,title)

  plt = Plots.plot(title = title)
    sol1 = dic["names"][1]
    key = "$(sol1)_$(str)"

    baserow = dic["$(sol1)_$(str)"]

    for i in 2:length(dic["names"])
        name = dic["names"][i]
        rowname = "$name / $(sol1)"
        key = "$(name)_$(str)"
        thisrow = dic[key]
        threshrow = min.(thisrow, 10*baserow) ./ baserow
        Plots.plot!(sort(threshrow), label=rowname)

        print(str, " " ,rowname, " ")
        vec_stats(threshrow, thisrow)
  end
  display(plt)
end

function compare_two(dic,sol1,sol2)

  plt = Plots.plot(title = "Total Time, Relative")
    key1 = "$(sol1)_tot"
    key2 = "$(sol2)_tot"

    tot1 = dic[key1]
    tot2 = dic[key2]

    Plots.plot!(sort(tot1./tot2), label="$(sol1) / $(sol2)")
    Plots.plot!(sort(tot2./tot1), label="$(sol2) / $(sol1)")

  display(plt)

  plt = Plots.plot(title = "Total Time / nnz")
    key1 = "$(sol1)_tot"
    key2 = "$(sol2)_tot"

    tot1 = dic[key1]
    tot2 = dic[key2]
    ne = dic["ne"]

    Plots.plot!(sort(tot1./ne), label="$(sol1) / nnz")
    Plots.plot!(sort(tot2./ne), label="$(sol2) / nnz")

  display(plt)
end

"""
    sd = mask_dict(mcf,mask)

For example:
```
mask = mcf["nv"] .> 20000
@show sum(mask)
sd = mask_dict(mcf,mask)
```
"""
function mask_dict(dic,mask)
    sol = dic["names"][1]
    ntests = length(dic["$(sol)_tot"])

    subdict = Dict()

    for (key, value) in dic
        if length(dic[key]) == ntests
            subdict[key] = dic[key][mask]
        else
            subdict[key] = dic[key]
        end
    end
    return subdict
end


text:
module SVDD

include("svdd_base.jl")
include("svdd_util.jl")

include("classifiers/classifier.jl")
include("classifiers/classifier_svdd.jl")
include("classifiers/classifier_svdd_vanilla.jl")
include("classifiers/classifier_svdd_neg.jl")
include("classifiers/classifier_ssad.jl")
include("classifiers/classifier_randomOC.jl")

include("init_strategies/strategies_C.jl")
include("init_strategies/strategies_gamma.jl")
include("init_strategies/strategies_combined.jl")

include("solvers/solver_base.jl")
include("solvers/smo_svdd.jl")

using Memento
using LinearAlgebra, Random

const LOGGER = getlogger(@__MODULE__)

function __init__()
    Memento.register(LOGGER)
end

export
    OCClassifier,
    VanillaSVDD,
    SVDDneg,
    SSAD,

    SMOSolver,

    ModelSolverException, ModelStateException,

    FixedParameterInitialization,
    FixedGammaStrategy,
    FixedCStrategy,
    SimpleCombinedStrategy,
    GammaFirstCombinedStrategy,
    BinarySearchCStrategy,
    TaxErrorEstimate,
    BoundedTaxErrorEstimate,
    RuleOfThumbSilverman,
    RuleOfThumbScott,

    initialize!,
    fit!,
    predict,
    classify,

    instantiate,
    set_adjust_K!,
    set_data!,
    set_pools!,
    set_kappa!,
    merge_pools,
    get_kernel,
    get_model_params,

    adjust_kernel_matrix!
end

ext: jl


text:
@testset "Parallel.Betweenness" begin

    s2 = SimpleDiGraph(3)
    add_edge!(s2, 1, 2); add_edge!(s2, 2, 3); add_edge!(s2, 3, 3)
    s1 = SimpleGraph(s2)
    g3 = PathGraph(5)
    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")
    for g in testdigraphs(gint)
        z  = @inferred(LightGraphs.betweenness_centrality(g))
        zp = @inferred(Parallel.betweenness_centrality(g))
        @test all(isapprox.(z, zp))

        y  = LightGraphs.betweenness_centrality(g, endpoints=true, normalize=false)
        yp = Parallel.betweenness_centrality(g, endpoints=true, normalize=false)
        @test all(isapprox(y, yp))

        xp = Parallel.betweenness_centrality(g, 3)
        xp2 = Parallel.betweenness_centrality(g, collect(1:20))

        @test length(xp) == 50
        @test length(xp2) == 50
    end
    @test @inferred(Parallel.betweenness_centrality(s1)) == [0, 1, 0]
    @test Parallel.betweenness_centrality(s2)   == [0, 0.5, 0]

    g = SimpleGraph(2)
    add_edge!(g, 1, 2)
    z  = LightGraphs.betweenness_centrality(g; normalize=true)
    zp = Parallel.betweenness_centrality(g; normalize=true)
    @test all(isapprox(z, zp))
    @test z[1] == z[2] == 0.0
    z2  = LightGraphs.betweenness_centrality(g, vertices(g))
    zp2 = Parallel.betweenness_centrality(g, vertices(g))
    @test all(isapprox(z2, zp2))
    z3  = LightGraphs.betweenness_centrality(g, [vertices(g);])
    zp3 = Parallel.betweenness_centrality(g, [vertices(g);])
    @test all(isapprox(z3, zp3))

    @test z == z2 == z3

    z  = LightGraphs.betweenness_centrality(g3; normalize=false)
    zp = Parallel.betweenness_centrality(g3; normalize=false)
    @test all(isapprox(z, zp))
    @test z[1] == z[5] == 0.0

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Access <: SwaggerModel

    function Access(;)
        o = new()
        o
    end
end # type Access

const _property_map_Access = Dict{Symbol,Symbol}()
const _property_types_Access = Dict{Symbol,String}()
Base.propertynames(::Type{ Access }) = collect(keys(_property_map_Access))
Swagger.property_type(::Type{ Access }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Access[name]))}
Swagger.field_name(::Type{ Access }, property_name::Symbol) =  _property_map_Access[property_name]

function check_required(o::Access)
    true
end

function validate_property(::Type{ Access }, name::Symbol, val)
end


text:
using Test
using StaticArrays: SVector
using LinearAlgebra: dot
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Essentials.DegreesOfFreedom
using QuantumLattices.Essentials.Spatials: PID,Point,pidtype,rcoord,icoord
using QuantumLattices.Mathematics.AlgebraOverFields: ID,idpropertynames
import QuantumLattices.Interfaces: dimension,decompose,update!,sequence,reset!
import QuantumLattices.Essentials.DegreesOfFreedom: script,latexsuperscript,latexsubscript

struct DID <: IID nambu::Int end
Base.adjoint(sl::DID)=DID(3-sl.nambu)

struct DIndex{S} <: Index{PID{S},DID}
    scope::S
    site::Int
    nambu::Int
end
Base.fieldnames(::Type{<:DIndex})=(:scope,:site,:nambu)
Base.union(::Type{P},::Type{I}) where {P<:PID,I<:DID}=DIndex{fieldtype(P,:scope)}
function Base.angle(id::OID{<:DIndex},vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float})
    phase=  length(vectors)==1 ? 2pi*dot(decompose(id.icoord,vectors[1]),values) :
            length(vectors)==2 ? 2pi*dot(decompose(id.icoord,vectors[1],vectors[2]),values) :
            length(vectors)==3 ? 2pi*dot(decompose(id.icoord,vectors[1],vectors[2],vectors[3]),values) :
            error("angle error: not supported number of input basis vectors.")
    id.index.nambu==1 ? phase : -phase
end

script(oid::OID{<:DIndex},::Val{:site})=oid.index.site
script(oid::OID{<:DIndex},::Val{:nambu})=oid.index.nambu==2 ? "\\dagger" : ""

struct DFock <: Internal{DID}
    atom::Int
    nnambu::Int
end
dimension(f::DFock)=f.nnambu
Base.getindex(f::DFock,i::Int)=DID(i)

struct DOperator{V<:Number,I<:ID{OID}} <: Operator{V,I}
    value::V
    id::I
end
latexformat(DOperator,LaTeX{(:nambu,),(:site,)}('d'))

@testset "Index" begin
    index=DIndex(PID('S',4),DID(1))
    @test index|>pidtype==PID{Char}
    @test index|>typeof|>pidtype==PID{Char}
    @test index|>iidtype==DID
    @test index|>typeof|>iidtype==DID
    @test index|>pid==PID('S',4)
    @test index|>iid==DID(1)
    @test index|>adjoint==DIndex('S',4,2)
    @test union(PID,IID)==Index{PID,IID}
    @test union(index|>pidtype,index|>iidtype)==index|>typeof
end

@testset "IndexToTuple" begin
    index=DIndex(PID('S',4),DID(1))
    @test directindextotuple(index)==('S',4,1)
    filteredindextotuple=FilteredAttributes(DIndex)
    @test filteredindextotuple==FilteredAttributes(:scope,:site,:nambu)
    @test isequal(filteredindextotuple,FilteredAttributes(:scope,:site,:nambu))
    @test filteredindextotuple|>length==3
    @test filteredindextotuple|>typeof|>length==3
    @test filteredindextotuple(index)==('S',4,1)
    @test filter(attr->attr≠:scope,filteredindextotuple)(index)==(4,1)
    @test filter(attr->attr≠:nambu,filteredindextotuple)(index)==('S',4)
    @test filter(attr->attr∉(:site,:nambu),filteredindextotuple)(index)==('S',)
end

@testset "Internal" begin
    it=DFock(1,2)
    @test it|>eltype==DID
    @test it|>typeof|>eltype==DID
    @test it==deepcopy(it)
    @test isequal(it,deepcopy(it))
    @test it|>string=="DFock(atom=1,nnambu=2)"
    @test it|>collect==[DID(1),DID(2)]
end

@testset "IDFConfig" begin
    config=IDFConfig{DFock}(pid->DFock((pid.site-1)%2+1,2),[PID(1,1),PID(1,2)])
    @test convert(Dict,config)==Dict(PID(1,1)=>DFock(1,2),PID(1,2)=>DFock(2,2))
    reset!(config,(PID(2,1),PID(2,2)))
    @test convert(Dict,config)==Dict(PID(2,1)=>DFock(1,2),PID(2,2)=>DFock(2,2))
end

@testset "Table" begin
    by=filter(attr->attr≠:nambu,FilteredAttributes(DIndex))

    table=Table([DIndex(1,1,1),DIndex(1,1,2)])
    @test table==Dict(DIndex(1,1,1)=>1,DIndex(1,1,2)=>2)
    table=Table([DIndex(1,1,1),DIndex(1,1,2)],by)
    @test table==Dict(DIndex(1,1,1)=>1,DIndex(1,1,2)=>1)

    config=IDFConfig{DFock}(pid->DFock((pid.site-1)%2+1,2),[PID(1,1),PID(1,2)])
    inds1=(DIndex(PID(1,1),iid) for iid in DFock(1,2))|>collect
    inds2=(DIndex(PID(1,2),iid) for iid in DFock(2,2))|>collect
    @test Table(config)==Table([inds1;inds2])
    @test Table(config,by)==Table([inds1;inds2],by)
    @test Table(config)==union(Table(inds1),Table(inds2))
    @test Table(config,by)|>reverse==Dict(1=>Set([DIndex(1,1,1),DIndex(1,1,2)]),2=>Set([DIndex(1,2,1),DIndex(1,2,2)]))

    table=Table(config)
    @test reset!(empty(table),[inds1;inds2])==table
    @test reset!(empty(table),config)==table
end

@testset "OID" begin
    oid=OID(DIndex(1,1,1),rcoord=SVector(0.0,-0.0),icoord=SVector(0.0,0.0),seq=1)
    @test oid'==OID(DIndex(1,1,2),rcoord=SVector(0.0,0.0),icoord=SVector(0.0,0.0),seq=1)
    @test hash(oid,UInt(1))==hash(OID(DIndex(1,1,1),rcoord=SVector(0.0,0.0),icoord=SVector(0.0,0.0),seq=1),UInt(1))
    @test idpropertynames(ID{OID})==(:indexes,:rcoords,:icoords,:seqs)
    @test fieldnames(OID)==(:index,:rcoord,:icoord,:seq)
    @test string(oid)=="OID(DIndex(1,1,1),[0.0,0.0],[0.0,0.0],1)"
    @test ID(oid',oid)'==ID(oid',oid)
    @test isHermitian(ID(oid',oid))==true
    @test isHermitian(ID(oid,oid))==false
    @test oidtype(DID,Point{2,PID{Char}},Nothing,Val(true))==OID{DIndex{Char},SVector{2,Float},SVector{2,Float},Nothing}
    @test oidtype(DID,Point{2,PID{Char}},Table{DIndex{Char}},Val(true))==OID{DIndex{Char},SVector{2,Float},SVector{2,Float},Int}
    @test oidtype(DID,Point{2,PID{Char}},Nothing,Val(false))==OID{DIndex{Char},Nothing,Nothing,Nothing}
    @test oidtype(DID,Point{2,PID{Char}},Table{DIndex{Char}},Val(false))==OID{DIndex{Char},Nothing,Nothing,Int}
end

@testset "Operator" begin
    @test latexformat(DOperator)==LaTeX{(:nambu,),(:site,)}('d')

    opt=DOperator(1.0im,(DIndex(1,2,2),DIndex(1,1,1)),rcoords=(SVector(1.0,0.0),SVector(0.0,0.0)),icoords=(SVector(2.0,0.0),SVector(0.0,0.0)),seqs=(2,1))
    @test opt'==DOperator(-1.0im,(DIndex(1,1,2),DIndex(1,2,1)),rcoords=(SVector(0.0,0.0),SVector(1.0,0.0)),icoords=(SVector(0.0,0.0),SVector(2.0,0.0)),seqs=(1,2))
    @test isHermitian(opt)==false
    @test string(opt)=="DOperator(1.0im,ID(OID(DIndex(1,2,2),[1.0,0.0],[2.0,0.0],2),OID(DIndex(1,1,1),[0.0,0.0],[0.0,0.0],1)))"
    @test twist(opt,[[1.0,0.0],[0.0,1.0]],[0.1,0.0])≈replace(opt,value=1.0im*conj(exp(2im*pi*0.2)))
    @test sequence(opt)==(2,1)
    @test sequence(opt,Dict(DIndex(1,2,2)=>3,DIndex(1,1,1)=>4))==(3,4)

    opt=DOperator(1.0,(DIndex(1,1,2),DIndex(1,1,1)),rcoords=(SVector(0.5,0.5),SVector(0.5,0.5)),icoords=(SVector(1.0,1.0),SVector(1.0,1.0)),seqs=(1,1))
    @test opt'==opt
    @test isHermitian(opt)==true

    opt=DOperator(1.0,(DIndex(1,1,2),DIndex(1,1,1)),rcoords=(SVector(0.5,0.5),SVector(0.0,0.5)),icoords=(SVector(1.0,1.0),SVector(0.0,1.0)),seqs=(1,1))
    @test rcoord(opt)==SVector(0.5,0.0)
    @test icoord(opt)==SVector(1.0,0.0)

    opt=DOperator(1.0,ID(OID(DIndex(1,1,2),SVector(0.5,0.0),SVector(1.0,0.0),1)))
    @test rcoord(opt)==SVector(0.5,0.0)
    @test icoord(opt)==SVector(1.0,0.0)
end

@testset "Operators" begin
    opt1=DOperator(1.0im,(DIndex(1,2,2),DIndex(1,1,1)),rcoords=(SVector(1.0,0.0),SVector(0.0,0.0)),icoords=(SVector(2.0,0.0),SVector(0.0,0.0)),seqs=(2,1))
    opt2=DOperator(1.0,(DIndex(1,1,2),DIndex(1,1,1)),rcoords=(SVector(0.0,0.0),SVector(0.0,0.0)),icoords=(SVector(0.0,0.0),SVector(0.0,0.0)),seqs=(1,1))
    opts=Operators(opt1,opt2)
    @test summary(opts)=="Operators{$(opts|>valtype)}"
    @test opts'==Operators(opt1',opt2')
    @test opts'+opts==Operators(opt1,opt1',opt2*2)
    @test isHermitian(opts)==false
    @test isHermitian(opts'+opts)==true
end

@testset "LaTeX" begin
    latex=LaTeX{(:nambu,),(:site,)}('c')
    @test latexsuperscript(latex|>typeof)==(:nambu,)
    @test latexsubscript(latex|>typeof)==(:site,)

    oid=OID(DIndex('d',1,2),rcoord=SVector(0.0,0.0),icoord=SVector(1.0,0.0),seq=1)
    @test script(oid,Val(:rcoord))=="[0.0,0.0]"
    @test script(oid,Val(:icoord))=="[1.0,0.0]"
    @test script(oid,latex,Val(:B))=='c'
    @test script(oid,latex,Val(:SP))==("\\dagger",)
    @test script(oid,latex,Val(:SB))==(1,)
    @test repr(oid,latex)=="c^{\\dagger}_{1}"

    opt=DOperator(1.0,(DIndex('d',2,2),DIndex('d',1,1)),rcoords=(SVector(1.0,0.0),SVector(0.0,0.0)),icoords=(SVector(2.0,0.0),SVector(0.0,0.0)),seqs=(2,1))
    @test repr(opt,LaTeX{(:nambu,),(:site,)}('c'))=="c^{\\dagger}_{2}c^{}_{1}"
    @test repr(opt,LaTeX{(:nambu,),(:site,)}())=="d^{\\dagger}_{2}d^{}_{1}"
    @test repr(opt,LaTeX{(:nambu,),(:rcoord,)}())=="d^{\\dagger}_{[1.0,0.0]}d^{}_{[0.0,0.0]}"
    @test repr(opt)=="d^{\\dagger}_{2}d^{}_{1}"

    opts=Operators(DOperator(1.0-1.0im,(DIndex('d',2,2),DIndex('d',1,1))),DOperator(-1.0,(DIndex('d',1,2),DIndex('d',1,1))))
    @test repr(opts,LaTeX{(:nambu,),(:site,)}('c'))=="(1.0-1.0im)c^{\\dagger}_{2}c^{}_{1}-c^{\\dagger}_{1}c^{}_{1}"

    io=IOBuffer()
    show(io,MIME"text/latex"(),opt)
    @test String(take!(io))=="\$d^{\\dagger}_{2}d^{}_{1}\$"
    show(io,MIME"text/latex"(),opts)
    @test String(take!(io))=="\$(1.0-1.0im)d^{\\dagger}_{2}d^{}_{1}-d^{\\dagger}_{1}d^{}_{1}\$"
end

@testset "Boundary" begin
    opt=DOperator(4.5,(DIndex('a',1,2),DIndex('b',2,1)),rcoords=(SVector(0.5,0.5),SVector(1.5,1.5)),icoords=(SVector(0.0,0.0),SVector(1.0,1.0)),seqs=(1,2))
    bound=Boundary{(:θ₁,:θ₂)}([0.1,0.2],[[1.0,0.0],[0.0,1.0]])
    @test bound==deepcopy(bound) && isequal(bound,deepcopy(bound))
    @test angle(bound,opt)≈0.6pi
    @test bound(opt)≈replace(opt,value=4.5*exp(2im*pi*0.3))
    update!(bound,θ₁=0.3)
    @test bound(opt)≈replace(opt,value=4.5*exp(2im*pi*0.5))
    bound=Boundary()
    @test angle(bound,opt)==0
    @test bound(opt)==opt
    @test update!(bound)==bound
end


text:
export PointMassFormConstraint, call_optimizer, call_starting_point

using Distributions
using Optim

"""
    PointMassFormConstraint

One of the form constraint objects. Constraint a message to be in a form of dirac's delta point mass.
By default uses `Optim.jl` package to find argmin of -logpdf(x).
Accepts custom `optimizer` callback which might be used to customise optimisation procedure with different packages
or different arguments for `Optim.jl` package.

See also: [`constrain_form`](@ref), [`DistProduct`](@ref)
"""
struct PointMassFormConstraint{F, P}
    optimizer      :: F
    starting_point :: P
end

PointMassFormConstraint(; optimizer = default_point_mass_form_constraint_optimizer, starting_point = default_point_mass_form_constraint_starting_point) = PointMassFormConstraint(optimizer, starting_point)

default_form_check_strategy(::PointMassFormConstraint) = FormConstraintCheckLast()

is_point_mass_form_constraint(::PointMassFormConstraint) = true

call_optimizer(pmconstraint::PointMassFormConstraint, data)      = pmconstraint.optimizer(variate_form(data), value_support(data), pmconstraint, data)
call_starting_point(pmconstraint::PointMassFormConstraint, data) = pmconstraint.starting_point(variate_form(data), value_support(data), pmconstraint, data)

function constrain_form(pmconstraint::PointMassFormConstraint, message::Message)
    data       = ReactiveMP.getdata(message)
    is_clamped = ReactiveMP.is_clamped(message)
    is_initial = ReactiveMP.is_initial(message)
    return Message(call_optimizer(pmconstraint, data), is_clamped, is_initial)
end

function default_point_mass_form_constraint_optimizer(::Type{ Univariate }, ::Type{ Continuous }, constraint::PointMassFormConstraint, distribution)

    target = let distribution = distribution
        (x) -> -logpdf(distribution, x[1])
    end

    support = Distributions.support(distribution)

    result = if isinf(Distributions.minimum(support)) && isinf(Distributions.maximum(support))
        optimize(target, call_starting_point(constraint, distribution), LBFGS())
    else
        lb = [ Distributions.minimum(support) ]
        rb = [ Distributions.maximum(support) ]
        optimize(target, lb, rb, call_starting_point(constraint, distribution), Fminbox(GradientDescent()))
    end

    if Optim.converged(result)
        return PointMass(Optim.minimizer(result)[1])
    else
        error("Optimisation procedure for point mass estimation did not converge", result)
    end
end

function default_point_mass_form_constraint_starting_point(::Type{ Univariate }, ::Type{ Continuous }, constraint::PointMassFormConstraint, distribution)
    support = Distributions.support(distribution)
    lb      = Distributions.minimum(support)
    rb      = Distributions.maximum(support)
    return if isinf(lb) && isinf(rb)
        return zeros(1)
    else
        error("No default starting point specified for a range [ $(lb), $(rb) ]")
    end
end


text:
push!(LOAD_PATH,".")
using CodeGenerator

buffer=CodeGenerator.read_grammar()
g=CodeGenerator.Grammar(buffer)
println(g)

CodeGenerator.prettyshow(stdout,g)

G=CodeGenerator.Generator(g)

CodeGenerator.bind(G,Symbol("<Number>"),()->rand(1:9))

res=CodeGenerator.generate(G,Symbol("Expr");max_depth=6)

for each in res
      println(each)
end

res=CodeGenerator.generate(G,Symbol("Expr");max_depth=2,generator=CodeGenerator.generate_bfs)

for each in res
      println(each)
end

max_stars_repo_name: yeruoforever/SimpleGrammar

max_issues_repo_name: yeruoforever/SimpleGrammar

max_forks_repo_name: yeruoforever/SimpleGrammar


text:
using BQCESubroutine: BQCESubroutine, Utilities
using BenchmarkTools
using LinearAlgebra
using YaoLocations
using Test

suite = BenchmarkGroup()
suite["generic broutine!"] = BenchmarkGroup()

macro suite(title, ex)
    return esc(suite_m(__module__, title, ex))
end

function suite_m(m::Module, title, ex)
    quote
        if !isdefined($m, :__suite_stack__)
            if !isdefined($m, :SUITE)
                global const SUITE = $BenchmarkGroup()
            end
            global const __suite_stack__ = [SUITE]
        end

        __suite__ = last(__suite_stack__)
        __suite__ = __suite__[$title] = $BenchmarkGroup()
        push!(__suite_stack__, __suite__)

        let
           $(record_benchmarkable(ex))
        end

        pop!(__suite_stack__)

        SUITE
    end
end

function record_benchmarkable(ex)
    ex isa Expr || return ex

    if ex.head === :block
        body = Expr(:block)
        count = 0
        for each in ex.args
            line = if each isa Expr && each.head === :macrocall && each.args[1] === Symbol("@benchmarkable")
                count += 1
                :(__suite__[$count] = $each)
            else
                each
            end
            push!(body.args, line)
        end
        return quote
            __suite__ = last(__suite_stack__)
            __suite__ = __suite__[$title] = $BenchmarkGroup()
            push!(__suite_stack__, __suite__)

            let
            $(record_benchmarkable(ex))
            end

            pop!(__suite_stack__)
        end
    elseif ex.head === :for
        error("not implemented")
    end

    return ret
end

@suite "generic broutine!" begin
    @suite "rand(ComplexF64, 1<<3)" begin
        st = rand(ComplexF64, 1 << 8)
        U = rand(ComplexF64, 1 << 3, 1 << 3)
        locs = Locations((1, 3, 4, 5, 6))[1:3]
        @benchmarkable broutine!(r, $U, $locs) setup=(r=copy($st))
    end
end

SUITE["generic broutine!"]["rand(ComplexF64, 1<<3)"][1]

tune!(SUITE)
results = run(SUITE, verbose = true, seconds = 1)
results["generic broutine!"]["rand(ComplexF64, 1<<3)"][1]


text:
function get_occ(ev::Vector{Float64}, nocc::Int64, Tbeta::Float64)
    #
    # usage: [occ,efermi] = get_occ(ev, nocc, Tbeta);
    #
    # Obtained from KSSOLV
    #
    # ZX: Transferred to julia language.
    nev = length(ev)
    tol = 1e-15
    maxiter = 100

    if (nev > nocc)
       #
       # use bisection to find efermi such that
       #       sum_i fermidirac(ev(i)) = nocc
       #
       ilb = nocc-1
       iub = nocc+1
       lb = ev[ilb]
       ub = ev[iub]
       #
       # make sure flb < nocc and fub > nocc
       #
       flb = sum(fermidirac(ev,lb,Tbeta))
       fub = sum(fermidirac(ev,ub,Tbeta))
       while ( (nocc-flb)*(fub-nocc) < 0 )
          #fprintf('getocc: initial bounds are off:\n');
          #fprintf('flb = #11.3e, fub = #11.3e, nocc = #d\n', flb,fub,nocc);
          if (flb > nocc)
             if (ilb > 1)
                ilb = ilb - 1
                lb = ev[ilb]
                flb = sum(fermidirac(ev,lb,Tbeta))
             else
                error("getocc: cannot find a lower bound for efermi, something is wrong")
                break
             end
          end;
          if (fub < nocc)
             if (iub < nev)
                iub = iub + 1
                ub  = ev[iub]
                fub = sum(fermidirac(ev,ub,Tbeta))
             else
                error("getocc: cannot find an upper bound for efermi, something is wrong, try increasing the number of wavefunctions in X0")
                break
             end
          end
       end
       #fprintf('flb = #11.3e, fub = #11.3e\n', flb, fub);
       #pause;
       efermi = (lb+ub)/2
       occ = fermidirac(ev,efermi,Tbeta)
       occsum = sum(occ)
       iter = 1
       while ( abs(occsum-nocc) > tol ) && ( iter < maxiter )
          # fprintf('iter = #d, efermi = #11.3e, sum = #11.3e\n', ...
                   # iter, efermi, occsum);
          # fprintf('lb = #11.3e, ub = #11.3e\n', lb, ub);
          if (occsum < nocc)
             lb = efermi
          else
             ub = efermi
          end
          efermi = (lb+ub)/2
          occ = fermidirac(ev,efermi,Tbeta)
          occsum = sum(occ)
          iter = iter + 1
       end;
    elseif (nev == nocc)
       occ    = ones(nocc,1)
       efermi = ev[nocc]
    else
       occ = []
       efermi = 0
       error("The number of eigenvalues in ev should be larger than nocc")
    end
    return (occ, efermi)
end

function fermidirac(ev,efermi,Tbeta::Float64)
#
# usage: [f,g] = fermidirac(ev, efermi, Tbeta);

f = 1 ./(1 .+exp.(Tbeta*(ev .- efermi)))

return f
    # if( nargout > 1 )
    # # Compute the derivative
    # g = Tbeta./(1+exp(Tbeta*(-ev+efermi))).*f;
    # end
    #
    # return [f,g]
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
export PolynomialZonotope, dim, σ,
       polynomial_order,
       order,
       linear_map,
       scale,
       minkowsi_sum

"""
    PolynomialZonotope{N}

Type that represents a polynomial zonotope.

### Fields

- `c`  -- starting point
- `E`  -- matrix of multi-indexed generators such that *all* indices have the same value
- `F`  -- matrix of multi-indexed generators such that *not all* indices have the same value
- `G`  -- matrix of single-indexed generators

### Notes

Polynomial zonotopes were introduced by M. Althoff in [1] and have been applied
as a non-convex set representation in the reachability problem of nonlinear ODEs.

Mathematically, a polynomial zonotope is defined as the tuple ``(c, E, F, G)``,
where:

- ``c ∈ \\mathbb{R}^n`` is the starting point (in some particular cases it corresponds
  to the center of a usual zonotope),

- ``E = [E^{[1]} ⋯ E^{[η]}]`` is an ``n × p × η(η+1)/2`` matrix with column-blocks

```math
E^{[i]} = [f^{([i], 1, 1, …, 1)} ⋯ f^{([i], p, p, …, p)}], \\qquad i = 1,…, η
```
called the matrix of *multi-indexed generators with equal indices*, where each
``f^{([i], k_1, k_2, …, k_i)}`` is an ``n``-vector,

- ``F = [F^{[2]} ⋯ F^{[η]}]`` is a matrix with column-blocks

```math
F^{[i]} = [f^{([i], 1, 1, …, 1, 2)} f^{([i], 1, 1, …, 1, 3)} ⋯ f^{([i], 1, 1, …, 1, p)} \\\\
f^{([i], 1, 1, …, 2, 2)} f^{([i], 1, 1, …, 2, 3)} ⋯ f^{([i], 1, 1, …, 2, p)} \\\\
f^{([i], 1, 1, …, 3, 3)} ⋯], \\qquad i = 1,…, η
```
called the matrix of *multi-indexed generators with unequal indices* (or, more accurately,
not-all-equal indices), where each ``f^{([i], k_1, k_2, …, k_i)}`` is an ``n``-vector,

- ``G = [G^{[1]} ⋯ G^{[q]}]`` is an ``n × q`` matrix with columns

```math
G^{[i]} = g^{(i)}, \\qquad i = 1,…, q
```
called the matrix of *single-indexed generators*, where each ``g^{(i)}`` is an
``n``-vector.

The polynomial zonotope ``(c, E, F, G)`` defines the set:

```math
\\mathcal{PZ} = \\left\\{ c + ∑_{j=1}^p β_j f^{([1], j)} + ∑_{j=1}^p ∑_{k=j}^p β_j β_k f^{([2], j, k)} + \\\\
+ … + ∑_{j=1}^p ∑_{k=j}^p ⋯ ∑_{m=ℓ}^p β_j β_k ⋯ β_m f^{([η], j, k, …, m)} + \\\\
+ ∑_{i=1}^q γ_i g^{(i)}, \\qquad β_i, γ_i ∈ [-1, 1] \\right\\},
```
where the number of factors in the final product, ``β_j β_k ⋯ β_m``, corresponds to
the polynomial order ``η``.

[1] M. Althoff in *Reachability analysis of nonlinear systems using conservative
    polynomialization and non-convex sets*, Hybrid Systems: Computation and
    Control, 2013, pp. 173–182.
"""
struct PolynomialZonotope{N}
    c::Vector{N}
    E::Vector{Matrix{N}}
    F::Vector{Matrix{N}}
    G::Matrix{N}

    # default constructor with dimension check
    function PolynomialZonotope{N}(c::Vector{N},
                                   E::Vector{Matrix{N}},
                                   F::Vector{Matrix{N}},
                                   G::Matrix{N}) where {N}

        # check polynomial order
        @assert length(E) == 1 + length(F)

        return new(c, E, F, G)
    end
end

# type-less convenience constructor
PolynomialZonotope(c::Vector{N},
                   E::Vector{Matrix{N}},
                   F::Vector{Matrix{N}},
                   G::Matrix{N}) where {N} = PolynomialZonotope{N}(c, E, F, G)

"""
    dim(pz::PolynomialZonotope)::Int

Return the ambient dimension of a polynomial zonotope.

### Input

- `pz` -- polynomial zonotope

### Output

An integer representing the ambient dimension of the polynomial zonotope.
"""
dim(pz::PolynomialZonotope)::Int = length(pz.c)

"""
    σ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::AbstractVector{N} where {N}

Return the support vector of a polynomial zonotope along direction `d`.

### Input

- `d`  -- direction
- `pz` -- polynomial zonotope

### Output

Vector representing the support vector.
"""
function σ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::AbstractVector{N} where {N}
    error("this function is not yet implemented")
end

"""
    ρ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::AbstractVector{N} where {N}

Return the support function of a polynomial zonotope along direction `d`.

### Input

- `d`  -- direction
- `pz` -- polynomial zonotope

### Output

Value of the support function.
"""
function ρ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::N where {N}
    error("this function is not yet implemented")
end

"""
    polynomial_order(pz::PolynomialZonotope)::Int

Polynomial order of a polynomial zonotope.

### Input

- `pz` -- polynomial zonotope

## Output

The polynomial order, defined as the maximal power of the scale factors ``β_i``.
Usually denoted ``η``.
"""
polynomial_order(pz::PolynomialZonotope)::Int = length(pz.E)

"""
    order(pz::PolynomialZonotope)::Rational{Int}

Order of a polynomial zonotope.

### Input

- `pz` -- polynomial zonotope

## Output

The order, a rational number defined as the total number of generators divided
by the ambient dimension.
"""
function order(pz::PolynomialZonotope)::Rational{Int}
    η = polynomial_order(pz)  # polynomial order
    p = size(pz.E[1], 2)  # number of dependent factors
    q = size(pz.G, 2)  # number of independent factors
    ξ = sum(i -> binomial(p+i-1, i), 1:η) + q
    n = dim(pz)
    return ξ//n
end

"""
    linear_map(M::Matrix, pz::PolynomialZonotope)

Return the linear map of a polynomial zonotope.

### Input

- `M`  -- matrix
- `pz` -- polynomial zonotope

## Output

Polynomial zonotope such that its starting point and generators are those of `pz`
multiplied by the matrix `M`.
"""
function linear_map(M::Matrix, pz::PolynomialZonotope)
    c = M * pz.c
    E = [M*Ei for Ei in pz.E]
    F = [M*Fi for Fi in pz.F]
    G = M * pz.G
    return PolynomialZonotope(c, E, F, G)
end

"""
    scale(α::Number, pz::PolynomialZonotope)

Return a polynomial zonotope modified by a scale factor.

### Input

- `α`  -- polynomial zonotope
- `pz` -- polynomial zonotope

## Output

Polynomial zonotope such that its center and generators are multiples of those
of `pz` by a factor ``α``.
"""
function scale(α::Number, pz::PolynomialZonotope)
    c = α * pz.c
    E = [α*Ei for Ei in pz.E]
    F = [α*Fi for Fi in pz.F]
    G = α * pz.G
    return PolynomialZonotope(c, E, F, G)
end

"""
    minkowski_sum(pz::PolynomialZonotope, z::Zonotope)

Return the Minkowski sum of a polynomial zonotope and a usual zonotope.

### Input

- `pz` -- polynomial zonotope
- `z`  -- usual zonotope

## Output

Polynomial zonotope such that its center is the sum of the centers of `pz` and `z`
and its generators is the matrix that results by concatenating them.
"""
function minkowski_sum(pz::PolynomialZonotope, z::Zonotope)
    c = pz.c + z.center
    G = [pz.G z.generators]
    return PolynomialZonotope(c, pz.E, pz.F, G)
end

# symmetric method
minkowski_sum(z::Zonotope, pz::PolynomialZonotope) = minkowski_sum(pz, z)


text:
"""
Neural network constructor
"""
mutable struct network
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
    transfers::Array{Function,1}
    function network(sizes, activation_functions)
        num_layers = length(sizes)
        sizearr = sizes
        biases = [randn(y) for y in sizes[2:end]]
        weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
        transfers = [activation_functions[i] for i in 1:length(activation_functions)]
        return new(num_layers, sizearr, biases, weights, transfers)
    end
end

"""
Evaluate network on input x
"""
function nn_forward(net::network, x)
    for (w, b, h) in zip(net.weights, net.biases, net.transfers[2:end])
        x = h.(w*x .+ b)
    end
    return x
end

"""
Log posterior for binary classification (loglikelihood + logprior)
Loglikelihood can be obtained from the loss function implemented in Flux.
"""
function log_posterior(net::network, x, y, prior_sigma)
    return sum(Flux.binarycrossentropy.(nn_forward(net,x), y)) +
                0.5*norm(Flux.params([net.weights, net.biases]))^2 / prior_sigma^2
end

"""
Function that implements one iteration of the
Hamiltonian Monte Carlo algorithm to update the weights.
Leapfrog acts on matrices and arrays
"""
function hmc_bnn!(net::network, X, Y, prior_sigma, step_size, num_steps)

    p_for_biases = [randn(y) for y in net.sizearr[2:end]]
    p_for_weights = [randn(y, x) for (x, y) in zip(net.sizearr[1:end-1], net.sizearr[2:end])]
    current_K = 0.5 * norm([p_for_biases, p_for_weights])^2

    current_weights = copy(net.weights)
    current_biases = copy(net.biases)
    current_U = log_posterior(net, X, Y, prior_sigma)

    grads = Zygote.gradient(m -> log_posterior(m, X, Y, prior_sigma), net)[1][]
    # new state according to leapfrog for hamilton dynamics
    p_for_weights =  p_for_weights .- step_size .* grads.weights / 2.0
    p_for_biases = p_for_biases .- step_size .* grads.biases / 2.0
    for jump in 1:1:num_steps-1
        net.weights = net.weights .+ step_size .* p_for_weights# gradK.(p)
        net.biases = net.biases .+ step_size .* p_for_biases

        # update gradients to use with new  p
        grads = Zygote.gradient(m -> log_posterior(m, X, Y, prior_sigma), net)[1][]
        p_for_weights = p_for_weights .- step_size .* grads.weights
        p_for_biases = p_for_biases .- step_size .* grads.biases
    end
    net.weights = net.weights .+ step_size .* p_for_weights # gradK.(p)
    net.biases = net.biases .+ step_size .* p_for_biases

    # update gradients to use with new  p
    grads = Zygote.gradient(m -> log_posterior(m, X, Y, prior_sigma), net)[1][]

    # the second half for the kinetic variables
    p_for_weights = p_for_weights .- step_size .* grads.weights / 2.0
    p_for_biases = p_for_biases .- step_size .* grads.biases / 2.0   # the second half for the kinetic variables

    p_for_weights = -p_for_weights
    p_for_biases = -p_for_biases

    # test wethear we accept the new state or not
    # Evaluate potential and kinetic energies at start and end of trajectory
    proposed_U = log_posterior(net, X, Y, prior_sigma)
    proposed_K = 0.5 * norm([p_for_biases, p_for_weights])^2
    if  rand(Float64) .< exp(current_U - proposed_U + current_K - proposed_K)
        nothing
    else
        net.weights = current_weights
        net.biases = current_biases
    end
end

"""
Some helper functions.
"""
normalize_data(x) = (x .- mean(x)) ./ std(x);

function bnn_posterior!(net::network, weights, biases)
    net.weights = weights
    net.biases = biases
end

function cell2array(Z)
    m, n = size(Z)
    ZZ = zeros(m, n)
    for i in 1:1:m
        for j in 1:1:n
            ZZ[i,j] = Z[i,j][1]
        end
    end
    return ZZ
end


text:
module FindDifferences

using Images

"""
2つの画像を見つけ出す
2つの画像は重ならない
widthを半分にして，探す
差分率が最小となる二つの矩形を探索する
"""
function find_two_images(fn_img; max_count=3)
    dn_img = dirname(fn_img)
    fn_img1 = joinpath(dn_img, "img1.png")
    fn_img2 = joinpath(dn_img, "img2.png")
    fn_img_diff = joinpath(dn_img, "img_diff.png")
    # size(img) => (height, width)
    img = load(fn_img)
    img1t = nothing
    img2t = nothing

    height_0, width_0 = size(img)
    height = height_0
    width = floor(Int, width_0 / 2)

    diff_rate_min = Inf
    count_no_update = 0

    # widthを小さくしていく
    # widthを小さくしていき，小さくしても最小値が更新されない場合はstop
    while width > 0 && count_no_update <= max_count
        @show width
        begin_idx_width_1 = 0

        while begin_idx_width_1 + width < width_0
            # img1をずらしていく begin_idx_width_1を増やしていく
            begin_idx_width_1 += 1
            end_idx_width_1 = begin_idx_width_1 + width - 1
            # end_idx_width_1 + widthが元の画像サイズを超えていたらcontinue
            if end_idx_width_1 + width > width_0; continue; end
            img1 = img[1:height, begin_idx_width_1:end_idx_width_1]

            begin_idx_width_2 = 0 + end_idx_width_1
            while begin_idx_width_2 + width < width_0
                # img2をずらしていく
                begin_idx_width_2 += 1
                end_idx_width_2 = begin_idx_width_2 + width - 1
                img2 = img[1:height, begin_idx_width_2:end_idx_width_2]
                img_diff = img1 .- img2
                rgb_zero = zero(img_diff[1])
                # 差分が存在する面積を計算する
                is_diffs::BitArray{2} = img_diff .!= rgb_zero
                # 差分率 0 ~ 1
                diff_rate = sum(is_diffs) / length(is_diffs)
                if diff_rate < diff_rate_min
                    count_no_update = 0
                    diff_rate_min = diff_rate
                    @show diff_rate_min
                    # 画像を保存する
                    Images.save(fn_img1, img1)
                    Images.save(fn_img2, img2)
                    save_diff_image(fn_img_diff, img1, img2)
                    img1t = img1
                    img2t = img2
                end
            end
        end
        width -= 1
        count_no_update += 1
    end
    return img1t, img2t
end

"""
差分があるところを分かりやすくして保存する
"""
function save_diff_image(fn_img, img1, img2)
    img_diff = img1 .- img2
    rgb_zero = zero(img_diff[1])
    rgb_red = isa(rgb_zero, RGB) ? RGB(1, 0, 0) : RGBA(1, 0, 0, 0)
    # 差分があるところを塗り替える
    idx_diff = img_diff .!= rgb_zero
    img = copy(img1)
    img[idx_diff] .= rgb_red
    Images.save(fn_img, img)
end

end # module

max_stars_repo_name: Soyukke/FindDifferences

max_issues_repo_name: Soyukke/FindDifferences

max_forks_repo_name: Soyukke/FindDifferences


text:
using DFTShims: ColinearSpin, SpinDegenerate, is_spin_polarized
using AxisArrays
using Unitful

const DH = DFTShims.Dispatch.Hartree
const Dρ = DH.Scalars.ρ
const Dϵ = DH.Scalars.ϵ
const SIZES = 10, 2
const AXES = Axis{:radius}(1:SIZES[1]), Axis{:bb}((:α, :β))
const polarized = ColinearSpin()
const unpolarized = SpinDegenerate()

@testset "SpinDegenerate" begin
    ρ = zeros(Dρ{Int64}, false, SIZES...)
    @test typeof(ρ) <: AxisArray{Dρ{Int64}}
    @test size(ρ) == SIZES

    ρ = zeros(Dρ{Int64}, unpolarized, SIZES[1], AXES[1], SIZES[2])
    @test typeof(ρ) <: AxisArray{Dρ{Int64}}
    @test size(ρ) == SIZES
    @test axes(ρ, 1) == AXES[1]

    ρ = zeros(Dρ{Int64}, unpolarized, SIZES..., AXES[1], AXES[2])
    @test typeof(ρ) <: AxisArray{Dρ{Int64}}
    @test size(ρ) == SIZES
    @test axes(ρ, 1) == AXES[1]
    @test axes(ρ, 2) == AXES[2]

    ρ = @inferred reinterpret(Dρ{Int64}, SpinDegenerate(), [1 1 1; 2 2 2])
    @test typeof(ρ) <: AxisArray
    @test eltype(ρ) == Dρ{Int64}
    @test !is_spin_polarized(ρ)

    @test_throws ArgumentError zeros(Dρ{Int64}, unpolarized, SIZES..., AXES..., AXES[1])

    @inferred zeros(Dρ{Int64}, unpolarized, SIZES)
    @inferred zeros(Dρ{Int64}, unpolarized, SIZES, AXES)
    @inferred zeros(Dρ{Int64}, unpolarized, SIZES, AXES[1:1])

    ϵ = @inferred zeros(Dϵ{Int64}, ColinearSpin(), SIZES)
    @test is_spin_polarized(ϵ)
    @test :spin ∈ axisnames(ϵ)
    @test size(ϵ) == (SIZES..., 2)

    @test @inferred(SpinCategory(ρ)) == SpinDegenerate()
end

@testset "Axis Manipulations" begin
    a = AxisArray(zeros(2, 3))
    const add_spin_axis = DFTShims.ArrayInitialization.add_spin_axis
    const replace_spin_axis = DFTShims.ArrayInitialization.replace_spin_axis
    @test @inferred(add_spin_axis(polarized, (2, 3), 5)) == (2, 3, 5)
    @test @inferred(add_spin_axis(polarized, AXES, AXES[1])) == (AXES..., AXES[1])
    @test @inferred(add_spin_axis(ColinearSpinFirst(), AXES, AXES[1])) ==
        (AXES[1], AXES...)

    saxes = AXES..., Axis{:spin}((:u, :d))
    actual = @inferred replace_spin_axis(Dρ, polarized, (4, 3, 2), saxes)
    @test actual == ((4, 3, 2), (AXES..., Axis{:spin}((:α, :β))))
    actual = @inferred replace_spin_axis(Dρ, ColinearSpinFirst(), (4, 3, 2), saxes)
    @test actual == ((2, 4, 3), (Axis{:spin}((:α, :β)), AXES...))
    actual = @inferred replace_spin_axis(Dρ, ColinearSpinLast(), (4, 3, 2), saxes)
    @test actual == ((4, 3, 2), (AXES..., Axis{:spin}((:α, :β))))

    saxes = Axis{:spin}((:u, :d)), AXES...
    actual = @inferred replace_spin_axis(Dρ, ColinearSpinFirst(), (2, 3, 4), saxes)
    @test actual == ((2, 3, 4), (Axis{:spin}((:α, :β)), AXES...))
    actual = @inferred replace_spin_axis(Dρ, ColinearSpinLast(), (2, 3, 4), saxes)
    @test actual == ((3, 4, 2), (AXES..., Axis{:spin}((:α, :β))))
end

@testset "ColinearSpin" begin
    ρ = zeros(Dρ{Int64}, true, SIZES...)
    @test typeof(ρ) <: AxisArray{Dρ{Int64}}
    @test size(ρ) == (SIZES..., 2)
    @test axes(ρ, 3) == Axis{:spin}((:α, :β))
    @test @inferred(SpinCategory(ρ)) == ColinearSpin()

    ρ = @inferred zeros(Dρ{Int64}, polarized, AXES)
    @test typeof(ρ) <: AxisArray{Dρ{Int64}}
    @test size(ρ) == (SIZES..., 2)
    @test axes(ρ, 1) == AXES[1]
    @test axes(ρ, 2) == AXES[2]
    @test axes(ρ, 3) == Axis{:spin}((:α, :β))

    @test_throws(ArgumentError,
                 zeros(Dρ{Int64}, polarized, SIZES..., AXES..., Axis{:aa}((1, 2))))

    const D∂³ϵ_∂ρ∂σ² = DH.Scalars.∂³ϵ_∂ρ∂σ²
    ∂³ϵ_∂ρ∂σ² = @inferred zeros(D∂³ϵ_∂ρ∂σ²{Int16}, polarized, AXES)
    @test size(∂³ϵ_∂ρ∂σ²) == (SIZES..., 12)
    @test axes(∂³ϵ_∂ρ∂σ², 1) == AXES[1]
    @test axes(∂³ϵ_∂ρ∂σ², 2) == AXES[2]

    ρ = @inferred reinterpret(Dρ{Int64}, ColinearSpinFirst(), [1 1 1; 2 2 2])
    @test typeof(ρ) <: AxisArray
    @test eltype(ρ) == Dρ{Int64}
    @test is_spin_polarized(ρ)
    @test @inferred(SpinCategory(ρ)) == ColinearSpinFirst()

    ∂ϵ_∂ρ = @inferred rand(typeof(1u"eV*nm^3"), ColinearSpin{length(AXES)}(), AXES)
    @test axes(∂ϵ_∂ρ, length(AXES)) isa Axis{:spin}
    @test ndims(∂ϵ_∂ρ) == length(AXES) + 1

    @test_throws(ArgumentError,
                 reinterpret(Dρ{Int64}, ColinearSpinLast(), [1 1 1; 2 2 2]))

    @inferred zeros(Dρ{Int64}, polarized, SIZES)
    @inferred zeros(Dρ{Int64}, polarized, AXES)
    @inferred rand(typeof(1u"eV*nm^3"), ColinearSpin{length(AXES)}(), AXES)
    @inferred zeros(Dρ{Int64}, polarized, SIZES, AXES)
    @inferred zeros(Dρ{Int64}, polarized, SIZES, AXES[1:1])
    @inferred zeros(Dρ{Int64}, polarized, SIZES, AXES[1:1])
end

@testset "From template array" begin
    ρ = zeros(Dρ{Int32}, true, SIZES..., AXES...)
    @test all(zeros(ρ, DH.Scalars.∂²ϵ_∂σ², SpinAware()) .== 0u"∂²ϵ_∂σ²")
    @test !(typeof(zeros(ρ, DH.Scalars.∂²ϵ_∂σ², SpinAware()).data) <: AxisArray)
    @test eltype(zeros(ρ, DD.Scalars.∂²ϵ_∂σ², SpinAware())) == typeof(Int32(0)u"∂²ϵ_∂σ²")
    @test eltype(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinAware())) == typeof(0u"∂²ϵ_∂σ²")
    @test is_spin_polarized(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinAware()))
    @test length(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinAware())[Axis{:spin}]) == 6
    @test find(x -> typeof(x) <: Axis{:spin},
               axes(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinAware()))) == [length(AXES) + 1]
    @test axes(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinAware()))[1:end - 1] == AXES
    @test axes(zeros(ρ, DH.Scalars.∂²ϵ_∂σ², ColinearSpinFirst()))[2:end] == AXES

    @test !is_spin_polarized(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinDegenerate()))
    @test !(typeof(zeros(DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinDegenerate()).data) <: AxisArray)
    @test ndims(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinDegenerate())) == length(SIZES)
    @test axes(zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int64}, SpinDegenerate())) == AXES

    ϵ = @inferred zeros(ρ, Dϵ, SpinDegenerate())
    @test !is_spin_polarized(ϵ)
    @test :spin ∉ axisnames(ϵ)
    @test size(ϵ) == SIZES
    @test !(typeof(ϵ.data) <: AxisArray)

    ρ = zeros(ρ, Dρ{Int64}, SpinDegenerate())
    @test axes(zeros(ρ, DD.Scalars.∂²ϵ_∂σ², ColinearSpinFirst()))[2:end] == AXES
    @test is_spin_polarized(zeros(ρ, DH.Scalars.∂²ϵ_∂σ², ColinearSpinFirst()))

    ϵ = @inferred zeros(ρ, Dϵ, ColinearSpin())
    @test is_spin_polarized(ϵ)
    @test :spin ∈ axisnames(ϵ)
    @test size(ϵ) == (SIZES..., 2)

    @inferred zeros(ρ, DH.Scalars.∂²ϵ_∂σ², SpinAware())
    @inferred zeros(ρ, DH.Scalars.∂²ϵ_∂σ²{Int32}, SpinAware())
    @inferred zeros(ρ, DD.Scalars.∂²ϵ_∂σ², ColinearSpin())
end

@testset "Convert between arrays" begin
    ρₙ = zeros(Dρ{Int32}, true, SIZES..., AXES...)
    ρₙ[:] = (1:length(ρₙ)) * oneunit(eltype(ρₙ))
    @test @inferred(convert(ColinearSpinLast(), ρₙ)) === ρₙ
    @test !(typeof(convert(ColinearSpinLast(), ρₙ).data) <: AxisArray)
    ρ₀  = @inferred convert(ColinearSpinFirst(), ρₙ)
    @test axes(ρ₀) == (axes(ρₙ, Axis{:spin}), AXES...)
    @test !(typeof(ρ₀.data) <: AxisArray)

    ρ = @inferred convert(ColinearSpinLast(), ρ₀)
    @test axes(ρ) == (AXES..., axes(ρₙ, Axis{:spin}))
    @test all(ρ .== ρₙ)
    @test @inferred(convert(ColinearSpinLast(), ρ)) === ρ

    @test_throws InexactError convert(typeof(1u"m^-3"), ColinearSpinFirst(), ρₙ)
    ρₘ = @inferred convert(typeof(1.0u"m^-3"), ColinearSpinFirst(), ρₙ)
    @test axes(ρₘ) == axes(ρ₀)
    @test ρₘ ≈ ρ₀
    @test !(ustrip.(ρₘ) ≈ ustrip.(ρ₀))

    ϵ = similar(ρ₀, DH.Scalars.ϵ, SpinAware())
    @test is_spin_polarized(ϵ)
    @test !(typeof(ϵ.data) <: AxisArray)
end

@testset "Creation via keyword arguments" begin
    axis = 0u"a₀":0.1u"a₀":4u"a₀"
    ρ = zeros(DH.Scalars.ρ{Float64}, false, radius=axis)
    @test eltype(ρ) == DH.Scalars.ρ{Float64}
    @test size(ρ) == (length(axis), )
    @test axes(ρ, 1) == Axis{:radius}(axis)

    ρ = rand(DH.Scalars.ρ{Float64}, false, radius=axis)
    @test eltype(ρ) == DH.Scalars.ρ{Float64}
    @test size(ρ) == (length(axis), )
    @test axes(ρ, 1) == Axis{:radius}(axis)
end

@testset "Wrap a dimensionless array" begin
    @test is_spin_polarized(wrap([1, 2, 3]u"m^-3")) == false
    @test eltype(wrap([1, 2, 3]u"m^-3")) <: DD.Scalars.ρ
    @test_throws ArgumentError wrap(ColinearSpin(), [1, 2, 3]u"m^-3")
    @test is_spin_polarized(wrap(ColinearSpinFirst(), [1 2 3; 4 5 6]u"m^-3"))
    @test is_spin_polarized(wrap(DD.Scalars.ρ, ColinearSpinFirst(), [1 2 3; 4 5 6]))
    @test_throws ArgumentError wrap(DD.Scalars.ρ, ColinearSpin(), [1 2 3; 4 5 6])
end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
module ParameTeX

using CSV
using DataFrames

function readParamFile( fileName )
  return DataFrame( CSV.File( fileName ) );
end

function generatePdfs( df , folder, baseTeX )
    dfheaders = names( df );
    rows = eachrow( df )
    numCols = length( dfheaders )
    paramsFile = joinpath( folder, "parameters.tex" )

    currFolder = pwd()

    print(" Starts row by row generation...\n")
    for row in rows

        # open file of parameters and write values
        paramsFileFid = open( paramsFile , "w" )
        for col in 1:numCols
            if typeof(row[col]) <:AbstractString
                valueAsString = replace( row[col], "&"=>"\\&")

            elseif typeof(row[col]) == Int64
                valueAsString = string( row[col] )
            else
                error("type not implemented yet, create and issue!")
            end
            write( paramsFileFid, "\\newcommand{\\" * dfheaders[col] * "}{" * valueAsString * "}\n")
        end
        close( paramsFileFid );

        # compile pdf
        fileExtension = baseTeX[end-3:end]
        fileName      = baseTeX[1:end-4]
        outputPdfFilename = fileName * "_" * replace( row[1]," "=>"_" ) * ".pdf"
        cd(folder)
        print(outputPdfFilename,"\n")
        run(`pdflatex $baseTeX`)
        run(`mv $fileName.pdf $outputPdfFilename`)
        cd(currFolder)
    end

end

export readParamFile
export generatePdfs

end # module


text:
# This shows how to plot a 3D seismic tomography model
# The paper that describes it is:
#
# Zhao, L., Paul, A., Malusà, M.G., Xu, X., Zheng, T., Solarino, S., Guillot, S., Schwartz, S., Dumont, T., Salimbeni, S., Aubert, C., Pondrelli, S., Wang, Q., Zhu, R., 2016. Continuity of the Alpine slab unraveled by high-resolution P wave tomography. Journal of Geophysical Research: Solid Earth 121, 8720–8737. doi:10.1002/2016JB013310
#
# You will need to download the file `Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt` from:
#   https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/
#
# And make sure that you are in the same directory as the data file

using DelimitedFiles, GeophysicalModelGenerator

# Load data:
data=readdlm("Zhao_etal_JGR_2016_Pwave_Alps_3D_k60.txt",' ',Float64,'\n', skipstart=0,header=false)
lon        = data[:,1];
lat        = data[:,2];
depth      = data[:,3];
dVp_perc   = data[:,4];

# Create 3D regular grid:
resolution =  (length(unique(lon)), length(unique(lat)), length(unique(depth)))
Lon          = reshape(lon,      resolution);
Lat          = reshape(lat,      resolution);
Depth        = reshape(depth,    resolution);
dVp_perc_3D  = reshape(dVp_perc, resolution);

# save paraview file
Data_set    =   GeoData(Lon,Lat,Depth,(dVp_Percentage=dVp_perc_3D,))
Write_Paraview(Data_set, "Zhao_etal_2016_dVp_percentage")

# extract cross-sections
Data_cross  =   CrossSection(Data_set, Depth_level=-100km)
Write_Paraview(Data_cross, "Zhao_CrossSection_100km")

Data_cross  =   CrossSection(Data_set, Lon_level=10)
Write_Paraview(Data_cross, "Zhao_CrossSection_Lon10")

Data_cross  =   CrossSection(Data_set, Lon_level=10, Interpolate=true)
Write_Paraview(Data_cross, "Zhao_CrossSection_Lon10_interpolated");

Data_cross  =   CrossSection(Data_set, Start=(1.0,39), End=(18,50))
Write_Paraview(Data_cross, "Zhao_CrossSection_diagonal")

# Extract a 3D subset of the data
Data_subset     =   ExtractSubvolume(Data_set,Lon_level=(5,12), Lat_level=(40,45))
Write_Paraview(Data_subset, "Zhao_Subset")

Data_subset_interp     =   ExtractSubvolume(Data_set,Lon_level=(5,12), Lat_level=(40,45), Interpolate=true)
Write_Paraview(Data_subset, "Zhao_Subset_interp")


text:
#reads in the ibis data, th
function prepibes(;
  ibesname::String = PARAM[:ibesname],
  ibespath::String = PARAM[:ibespath],
  cilname::String = PARAM[:cilname],
  cilpath::String = PARAM[:cilpath],
  incsvstream::Function = IN_CSV_STREAM,
  csvextension::String = CSV_EXTENSION,
  ibescolumns::Union{Vector{Symbol}, Nothing} = PARAM[:ibescolumns])::DataFrame

  ibes::DataFrame = incsvstream("$ibespath\\$ibesname.$csvextension") |> CSV.File |> DataFrame
  rename!(ibes, (s::String->lowercase(s)).(names(ibes)))
  ibes.actdats = parseibes.(Date, ibes.actdats)
  ibes.anndats = parseibes.(Date, ibes.anndats)
  ibes.perioddate = parseibes.(Date, ibes.pends)

  #according to WRDS, announce date is used after 1/1/1999
  cutoffdate::Date = Date(1999,1,1)
  function pickdate(act,ann)
    #check if one is missing
    (ann === missing || act === missing) && return coalesce(ann, act, missing)

    ann ≥ cutoffdate && return ann
    return act
  end

  ibes.ibesadate = pickdate.(ibes.actdats, ibes.anndats)

  #read the cil linking file
  cil::DataFrame = CSV.read("$cilpath\\$cilname.csv", DataFrame)
  rename!(cil, (s::String->lowercase(s)).(names(cil)))
  cil = cil[cil.score .≤ 4, :] #only want likely matches- see WRDS documentation
  cil.edate = parsecil.(Date, cil.edate)
  cil.sdate = parsecil.(Date, cil.sdate)

  #NOTE: finish the linkage
  ibes = innerjoin(ibes, cil, on=:ticker)
  validlink(d, sdate, ::Missing) = true
  validlink(d, sdate, edate) = (d≥sdate) & (d≤edate)

  ibes = ibes[validlink.(ibes.ibesadate, ibes.sdate, ibes.edate), :]

  #we want the US currency row if there are multiple rows due to fx restatement
  ibes.uscurr = ibes.curr_act .=== "USD"
  sort!(ibes, [:permno, :perioddate, :uscurr])
  ibes = combine(last, groupby(ibes, [:permno, :perioddate]))

  #drop the suffix
  rename!(ibes, (n->n => (n[(end-4):end] === "_last" ? n[1:(end-5)] : n)).(names(ibes)))
  rename!(ibes, [:sdate=>:cilsdate, :edate=>:ciledate, :ticker=>:ibesticker])

  #now patch in a gvkey
  #onyl care that the link is valid on the proposed ibes adate
  ccm = prepccm(ccmcolumns=[:gvkey, :permno, :ccmeffdate, :ccmenddate])
  ibes = innerjoin(ibes, ccm, on=:permno)
  ibes = ibes[(ibes.ccmeffdate .≤ ibes.ibesadate) .& (ibes.ibesadate .≤ ibes.ccmenddate), :]

  #knock out dups by taking the first adate for any gvkey-period
  sort!(ibes, [:gvkey, :perioddate, :ibesadate])
  ibes = combine(first, groupby(ibes, [:gvkey, :perioddate]))

  #validate integrity of link
  @assert all((ibes.perioddate .=== lastdayofmonth.(ibes.perioddate)) .|
    (ibes.perioddate .=== missing))
  @assert ibes[nonunique(ibes,[:permno, :perioddate]),:] |> isempty
  @assert ibes[nonunique(ibes,[:gvkey, :perioddate]),:] |> isempty
  @assert ibes[nonunique(ibes,[:ibesticker, :perioddate]),:] |> isempty

  (!(ibescolumns === nothing)) && select!(ibes, ibescolumns)

  return ibes
end

parseibes(::Type{Date}, s::String, ibesdateformat::DateFormat=dateformat"yyyymmdd"
  ) = Dates.Date(s, ibesdateformat)
parseibes(::Type{<:Any}, ::Missing) = missing
parseibes(::Type{Date}, i::Int, ibesdateformat::DateFormat=dateformat"yyyymmdd") = parseibes(
  Date, "$i", ibesdateformat)

parsecil(::Type{<:Any}, ::Missing) = missing
parsecil(::Type{Date}, s::String, cildateformat::DateFormat=dateformat"dduuuyyyy"
  ) = Dates.Date(s, cildateformat)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "Test likelihood function" begin
    distribution = truncated(Normal(5.0, 1.0), 0, Inf)
    image = generate_image((1000, 1000), 0.15, distribution)
    # 1.18 is a factor of safety since all this thing is purely random
    # and can never achieve its "true" maximum value.
    max_likelihood = 1.18*Sandcorn.likelihood(image, distribution; cutoff = 20, ncloud = 150)

    for _ in 1:5
        new_distribution = truncated(Normal(5 + 1.0 - 2rand(Float64),
                                            1 + 0.5 -  rand(Float64)),
                                     0, Inf)
        shifted_likelihood = Sandcorn.likelihood(image, new_distribution;
                                                 cutoff = 20,
                                                 ncloud = 150)
        @test shifted_likelihood < max_likelihood
    end
end


text:
import ArgParse

function logo(version::Int32)

    border_print_()
    print(
        "                   __   __ ___________                 _                        \n" *
        "                   \\ \\ / ##  ___| ___ \\               | |                       \n" *
        "                    \\ V / \\ `--.| |_/ / ___ _ __   ___| |__                     \n" *
        "                    /   \\  `--. \\ ___ \\/ _ \\ '_ \\ / __| '_ \\                    \n" *
        "                   / /^\\ \\/\\__/ / |_/ /  __/ | | | (__| | | |                   \n" *
        "                   \\/   \\/\\____/\\____/ \\___|_| |_|\\___|_| |_|                   \n\n",
    )
    border_print_()
    center_print_("Developed at Argonne National Laboratory", 79)
    center_print_("XSBench.jl is developed at Oak Ridge National Laboratory", 79)

    v::String = string("Version: %d", version)
    center_print_(v, 79)
    border_print_()

end

function border_print_()
    printf(
        "===================================================================" *
        "=============\n",
    )
end

# Prints Section titles in center of 80 char terminal
function center_print_(s::String, width::Int32)

    length::Int32 = length(s)
    for i = 0:(width-length)/2
        print(" ")
    end
    println(s)

end

function read_CLI(args::Array{String,1})::Inputs

    ## Check if user sets these
    user_g = 0
    default_lookups = 1
    default_particles = 1

    parse_settings = ArgParse.ArgParseSettings(
        "XSBench.jl argument settings.",
        exc_handler = ArgParse.debug_handler,
    )

    ArgParse.@add_arg_table! parse_settings begin
        "-g"
        arg_type = Int64
        default = Int64(11303)
        range_tester = x -> x > 0
        help = "Number of Grid Points, must be larger than 0"
        "-m"
        help = "Simulation Method"
        arg_type = String
        default = "history"
        "-l"
        help = "Lookups"
        arg_type = Int32
        default = Int32(34)
        range_tester = x -> x > 0
        # using upper case H since h is reserved when add_help is true
        # https://argparsejl.readthedocs.io/en/latest/argparse.html#general-settings
        "-H"
        help = "Hash Bins"
        arg_type = Int32
        default = Int32(10000)
        "-p"
        help = "Particles"
        arg_type = Int32
        default = Int32(500000)
        "-s"
        help = "HM"
        arg_type = String
        default = "large"
        "-G"
        help = "Grid Type"
        arg_type = String
        default = "unionized"
        "-b"
        help = "Binary Mode"
        arg_type = String
        default = "none"
        "-k"
        help = "Kernel Optimization Selection"
        arg_type = Int32
        default = Int32(0)
    end

    parsed_args = ArgParse.parse_args(args, parse_settings)

    default_lookups::Int32 = 1
    default_particles::Int32 = 1

    # Set up returning inputs struct
    inputs = Inputs()

    for (key, value) in parsed_args

        if key == "g"
            inputs.n_gridpoints = value

        elseif key == "m"
            inputs.simulation_method = value
            if value == "event"
                if default_lookups != 0 && default_particles != 0
                    inputs.lookups = inputs.lookups * inputs.particles
                    inputs.particles = 0
                end

            end

        elseif key == "l"
            inputs.lookups = value
            default_lookups = 0

        elseif key == "H"
            inputs.hash_bins = value

        elseif key == "p"
            inputs.particles = value
            default_particles = 0

        elseif key == "s"
            inputs.HM = value

        elseif key == "G"
            inputs.grid_type = value

        elseif key == "b"
            inputs.binary_mode = value

        elseif key == "k"
            inputs.kernel_id = value

        end

    end

    # Validate inputs, throw exceptions
    return inputs
end


text:
module SBMLCases

using JSON, SimSolver, DataFrames, CSV, DataStructures, StatsPlots

# files containing settings and tags of sbml models
const cases_db = "./cases.json"
const results_db = "./results.json"
# models' dirs and output dir paths
const cases_path = "./cases/semantic"
const output_path = "./cases/output"

# default backend
const default_backend = Val{:SimSolver}

include("../$cases_path/julia/model.jl")

### Metelkin

function case_build_errors(
    case::AbstractDict;
    build_dict::Vector{Any}
)
    id = case["name"]
    filename_regex = Regex("$id-sbml-l2v4\\.xml\$")

    # to get only errors which refer to cases id
    f = x -> x["level"] == "error" && ( # select only errors
        ( haskey(x["opt"], "space") && x["opt"]["space"] == "x$id" ) ||  # search id in space
        ( haskey(x["opt"], "filename") && occursin(filename_regex, x["opt"]["filename"]) ) # search id in filename
    )

    return filter(f, build_dict)
end

function case_sim_result(
    case::AbstractDict;
    cases_path::AbstractString,
    output_path::AbstractString
)
    result = OrderedDict()
    try
        status = solve_case(case; cases_path = cases_path, output_path = output_path)
        if status
            result["status"] = "SUCCESS"
            result["message"] = "Simulations meet the criteria"
            println("$(case["name"])...................success")
        else
            result["status"] = "TOLERANCE_FAIL"
            result["message"] = "Simulation tolerance test not passed"
            println("$(case["name"])...................failure")
        end
    catch e
        result["status"] = "ERROR"
        result["message"] = "Error while running model. $e"
        println("$(case["name"])...................error")
    end

    return result
end

########################## Upload cases from cases_db ########################

"""
    upload_cases(;db_path::AbstractString=cases_db)

Upload cases from `cases_db`.
"""
function upload_cases(;db_path::AbstractString=cases_db)
    f = open(db_path, "r")
    dict = JSON.parse(f, dicttype=OrderedDict)
    close(f)
    return dict
end

########################## Add new cases #####################################

"""
    add_cases(;
        cases_path::AbstractString=cases_path,
        cases_db::AbstractString=cases_db
    )

Add new cases from `cases_path` to `cases_db`.
"""
function add_cases(;
    cases_path::AbstractString=cases_path,
    cases_db::AbstractString=cases_db
)
    cases_dict = upload_cases(db_path=cases_db)
    new_cases = 0
    foreach(readdir(cases_path)) do case
        if !haskey(cases_dict, case)
            @show case
            cases_dict[case] = add_single_case(cases_path, case)
            new_cases += 1
        end
    end
    new_cases > 0 && save_as_json(cases_dict, cases_db)
    return nothing
end

function add_single_case(path::AbstractString, case::AbstractString)
    OrderedDict(
        "name" => case,
        "settings" => add_settings(path, case),
        "tags" => add_tags(path, case)
    )
end

function add_settings(path::AbstractString, case::AbstractString)
    settings = OrderedDict{String,Any}()
    open("$path/$case/$case-settings.txt", "r") do f
        for line in eachline(f)
            split_line = split(line, ":")

            k = first(split_line)
            v = last(split_line)

            if k in ["start", "duration", "absolute", "relative"]
                settings[k] = !isempty(v) ? parse(Float64, v) : nothing
            elseif k == "steps"
                settings[k] = !isempty(v) ? parse(Int64, v) : nothing
            elseif k in ["variables", "amount", "concentration"]
                settings[k] = String[]
                for v_i in split(v, ",")
                    push!(settings[k], strip(v_i))
                end
            elseif !isempty(k)
                @warn "Case $case: setting $k currently not supported"
            end
        end
    end
    return settings
end

function add_tags(path::AbstractString, case::AbstractString)
    tags = OrderedDict{String,Vector{String}}()
    tags_parsed = 0
    open("$path/$case/$case-model.m", "r") do f
        for line in eachline(f)
            split_line = split(line, ": ")
            k = strip(first(split_line))
            v = strip(last(split_line))

            if k in ["componentTags", "testTags"]
                tags[k] = String[]
                for v_i in split(v, ",")
                    push!(tags[k], strip(v_i))
                end
                tags_parsed += 1
            end
            tags_parsed == 2 && break
        end
    end
    return tags
end

########################## Save cases or results ############################

function save_as_json(dict::AbstractDict, json::AbstractString)
    #stringdata = JSON.json(dict)
    open(json, "w") do f
        JSON.print(f, dict, 4)
    end
    return nothing
 end

############# Run simulations and update results in results_db ##############

"""
    update_results(
        case::AbstractString,
        cases_dict::AbstractDict=upload_cases(db_path=results_db);
        cases_path::AbstractString=cases_path,
        cases_db::AbstractString=cases_db,
        results_db::AbstractString=results_db,
        backend::DataType=default_backend,
        kwargs...
    )

Reads `cases_db` to `cases_dict`, accesses the `case`,
solves it with the chosen `backend` solver and writes result to `results_db`.
"""
function update_results(
    case::AbstractString,
    cases_dict::AbstractDict=upload_cases(db_path=results_db);
    cases_path::AbstractString=cases_path,
    cases_db::AbstractString=cases_db,
    results_db::AbstractString=results_db,
    backend::DataType=default_backend,
    kwargs...
)
    cases_dict["cases"][case]["result"] = OrderedDict()
    try
        status = solve_case(cases_dict["cases"][case], backend; kwargs...)
        if status
            cases_dict["cases"][case]["result"]["status"] = "success"
            cases_dict["cases"][case]["result"]["message"] = ""
            println("$case...................success")
        else
            cases_dict["cases"][case]["result"]["status"] = "failure"
            cases_dict["cases"][case]["result"]["message"] = "tolerance test not passed"
            println("$case...................failure")
        end
    catch e
        cases_dict["cases"][case]["result"]["status"] = "error"
        cases_dict["cases"][case]["result"]["message"] = "Check the model: $e"
        println("$case...................error")
    finally
        save_as_json(cases_dict, results_db)
    end
    return nothing
end

"""
    update_results(
        cases_vec::Vector{String},
        cases_dict::AbstractDict=upload_cases(db_path=results_db);
        cases_path::AbstractString=cases_path,
        cases_db::AbstractString=cases_db,
        results_db::AbstractString=results_db,
        backend::DataType=default_backend,
        kwargs...
    )

Reads `cases_db` to `cases_dict`, accesses the selected cases from `cases_vec`,
solves it with the chosen `backend` solver and writes results to `results_db`.
"""
function update_results(
    cases_vec::Vector{String},
    cases_dict::AbstractDict=upload_cases(db_path=results_db);
    cases_path::AbstractString=cases_path,
    cases_db::AbstractString=cases_db,
    results_db::AbstractString=results_db,
    backend::DataType=default_backend,
    kwargs...
)

    for case in cases_vec
        update_results(
            case,
            cases_dict;
            cases_path=cases_path,
            cases_db=cases_db,
            results_db=results_db,
            backend=default_backend,
            kwargs...
        )
    end
    return nothing
end

"""
    update_results(
        cases_range::UnitRange,
        cases_dict::AbstractDict=upload_cases(db_path=results_db);
        cases_path::AbstractString=cases_path,
        cases_db::AbstractString=cases_db,
        results_db::AbstractString=results_db,
        backend::DataType=default_backend,
        kwargs...
    )

Reads `cases_db` to `cases_dict`, accesses the selected cases from `cases_range`,
solves it with the chosen `backend` solver and writes results to `results_db`.
"""
function update_results(
    cases_range::UnitRange,
    cases_dict::AbstractDict=upload_cases(db_path=results_db);
    cases_path::AbstractString=cases_path,
    cases_db::AbstractString=cases_db,
    results_db::AbstractString=results_db,
    backend::DataType=default_backend,
    kwargs...
)
    for case in cases_range
        update_results(
            lpad(case, 5, "0"),
            cases_dict;
            cases_path=cases_path,
            cases_db=cases_db,
            results_db=results_db,
            backend=default_backend,
            kwargs...
        )
    end
    return nothing
end

"""
    update_results(
        cases_dict::AbstractDict=upload_cases(db_path=results_db);
        include_test_tags::Vector{String}=String[],
        include_component_tags::Vector{String}=String[],
        exclude_test_tags::Vector{String}=String[],
        exclude_component_tags::Vector{String}=String[],
        cases_path::AbstractString=cases_path,
        cases_db::AbstractString=cases_db,
        results_db::AbstractString=results_db,
        backend::DataType=default_backend,
        kwargs...
    )

Reads `cases_db` to `cases_dict`, filters the cases according to `include` and `exclude` tags,
solves it with the chosen `backend` solver and writes results to `results_db`.
"""
function update_results(
    cases_dict::AbstractDict=upload_cases(db_path=results_db);
    include_test_tags::Vector{String}=String[],
    include_component_tags::Vector{String}=String[],
    exclude_test_tags::Vector{String}=String[],
    exclude_component_tags::Vector{String}=String[],
    cases_path::AbstractString=cases_path,
    cases_db::AbstractString=cases_db,
    results_db::AbstractString=results_db,
    backend::DataType=default_backend,
    kwargs...
)
    cases_vec = isempty(include_test_tags) &&
                isempty(include_component_tags) &&
                isempty(exclude_test_tags) &&
                isempty(exclude_component_tags) ? cases_dict.keys :
                    filter_cases(cases_dict;
                        include_test_tags=include_test_tags,
                        include_component_tags=include_component_tags,
                        exclude_test_tags=exclude_test_tags,
                        exclude_component_tags=exclude_component_tags,
                    )
    update_results(
        cases_vec,
        cases_dict;
        cases_path=cases_path,
        cases_db=cases_db,
        results_db=results_db,
        backend=default_backend,
        kwargs...
    )
end

# solve single model with SimSolver backend
function solve_case(
    case::AbstractDict,
    backend::Type{Val{:SimSolver}}=default_backend;
    cases_path::AbstractString=cases_path,
    output_path::AbstractString=output_path,
    alg::Symbol=:Vern9
)

    # include heta models code
    #=
    file_path = "$cases_path/julia/model.jl"
    isfile(file_path) ? include(file_path) : "Model file doesn't exist"
    eval(quote using SBMLCases.SimSolverPlatform end)
    =#
    #model = eval_model(SBMLCases.julia.models) invokelatest not needed ?

    case_name = case["name"]
    #!haskey(SimSolverPlatform.models, case_name) && throw("Model $case_name is not compiled")

    tspan = (case["settings"]["start"],case["settings"]["duration"])
    step = (case["settings"]["duration"]-case["settings"]["start"])/case["settings"]["steps"]

    #outputs
    outputs = Symbol[]
    for v in case["settings"]["variables"]
        if v in case["settings"]["amount"]
            push!(outputs, Symbol(v*"_"))
        else
            push!(outputs, Symbol(v))
        end
    end

    # add step option to SimSolver
    saveat = collect(range(case["settings"]["start"],case["settings"]["duration"]; step=step))
    subtask = SubTask(saveat, outputs)

    solver = Dict(
        :alg=>alg,
        :reltol=>1e-7, #case["settings"]["relative"],
        :abstol=>1e-14 #case["settings"]["absolute"],
        #:maxiters => 10^5,
        #:dtmax => step/2
    )

    stask = SimpleSTask(
        SBMLCases.SimSolverPlatform.models[Symbol(join(["x",case_name]))],
        NamedTuple(),
        subtask,
        tspan,
        solver=solver,
        evt_save=(false,false)
    )

    res = solve_task(stask)

    df_sim = sol_as_df(res)
    CSV.write("$output_path/$case_name.csv",  df_sim)
    df_ans = DataFrame!(CSV.File("$cases_path/$case_name/$case_name-results.csv"))

    # create plot
    p_ref = plot_results(df_sim, df_ans)
    savefig(p_ref, "$output_path/$case_name")

    return compare_results(case, df_sim, df_ans)

end

# tolerance test according to:
# https://github.com/sbmlteam/sbml-test-suite/blob/master/cases/semantic/README.md#tolerances-and-errors-for-timecourse-tests
function compare_results(case::AbstractDict, df_sim::DataFrame, df_ans::DataFrame)
    (i_length,j_length) = size(df_sim)
    for i in 1:i_length
        for j in 2:j_length
            if abs(df_sim[i,j]-df_ans[i,j]) > (case["settings"]["absolute"] + case["settings"]["relative"] * abs(df_ans[i,j]))
                return false
            end
        end
    end
    return true
end

"""
    filter_cases(
        cases_dict::AbstractDict;
        include_test_tags::Vector{String}=String[],
        include_component_tags::Vector{String}=String[],
        exclude_test_tags::Vector{String}=String[],
        exclude_component_tags::Vector{String}=String[],
    )

Functions replicates the behavior of SBML Runner filter.
It filters out names of the cases which satisfy `include` and `exclude` tags.
"""
function filter_cases(
    cases_dict::AbstractDict;
    include_test_tags::Vector{String}=String[],
    include_component_tags::Vector{String}=String[],
    exclude_test_tags::Vector{String}=String[],
    exclude_component_tags::Vector{String}=String[],
)

    function include_filter_func(x)
        (!isempty(include_component_tags) ? any(in(cases_dict[x]["tags"]["componentTags"]).(include_component_tags)) : true) ||
        (!isempty(include_test_tags) ? any(in(cases_dict[x]["tags"]["testTags"]).(include_test_tags)) : true)
    end

    function exclude_filter_func(x)
        (!isempty(exclude_component_tags) ? !any(in(cases_dict[x]["tags"]["componentTags"]).(exclude_component_tags)) : true) &&
        (!isempty(exclude_test_tags) ? !any(in(cases_dict[x]["tags"]["testTags"]).(exclude_test_tags)) : true)
    end

    return filter(x->include_filter_func(x) && exclude_filter_func(x), keys(cases_dict))
end

# update functions after uploading a new model
function eval_model(model_code)
    Model(
        (cons)->Base.invokelatest(model_code.start, cons),
        (du, u, p, t)->Base.invokelatest(model_code.ode, du, u, p, t),
        [eval_event(evt) for evt in model_code.events],
        (outputIds)->Base.invokelatest(model_code.saving, outputIds),
        model_code.default_constants
    )
end

eval_event(evt::TimeEvent) = TimeEvent((cons)->Base.invokelatest(evt.condition_func, cons), (integrator)->Base.invokelatest(evt.affect_func, integrator))
eval_event(evt::DEvent) = DEvent((u, t, integrator)->Base.invokelatest(evt.condition_func, u, t, integrator), (integrator)->Base.invokelatest(evt.affect_func, integrator))
eval_event(evt::CEvent) = CEvent((u, t, integrator)->Base.invokelatest(evt.condition_func, u, t, integrator), (integrator)->Base.invokelatest(evt.affect_func, integrator))

function results_to_df(res::AbstractDict, ran::UnitRange{Int64}=1:955)
    num = length(ran)
    id_str = Vector{String}(undef, num)
    comp_tags = Vector(undef, num)
    test_tags = Vector(undef, num)
    status = Vector{String}(undef, num)
    message = Vector{String}(undef, num)

    for j in ran
        i = j-first(ran)+1
        id_str[i] = lpad(j, 5, "0")
        comp_tags[i] = res[id_str[i]]["tags"]["componentTags"]
        test_tags[i] = res[id_str[i]]["tags"]["testTags"]
        status[i] = String(res[id_str[i]]["result"]["status"])
        message[i] = String(res[id_str[i]]["result"]["message"])
    end

    DataFrame(
        id = id_str,
        comp_tags = comp_tags,
        test_tags = test_tags,
        status = status,
        message = message
    )
end

function plot_results(df_sim, df_ans)
    names_sim = names(df_sim)[2:end]
    names_ans = names(df_ans)[2:end]
    cl = size(df_sim)[2]

    time_sim = :time
    time_ans = Symbol(names(df_ans)[1])

    p_sim = StatsPlots.@df df_sim plot(
        cols(time_sim),
        cols(2:cl),
        title = "Simulations",
        legend = false)
    p_ans = StatsPlots.@df df_ans plot(
        cols(time_ans),
        cols(2:cl),
        title = "Answers",
        legend = false)

    df_diff = copy(df_ans)
    for (col_sim,col_ans) in zip(names_sim,names_ans)
        df_diff[!,col_ans] .= abs.(df_sim[!,col_sim] - df_ans[!,col_ans])
    end
    p_diff = StatsPlots.@df df_diff plot(
        cols(time_ans),
        cols(2:cl),
        title = "Difference",
        legend = false)
    legend = plot(
        permutedims(zeros(Int, cl-1)),
        showaxis = false,
        grid = false,
        label = permutedims(names_ans))
    plot(p_sim, p_ans, p_diff, legend, dpi=300)
end

###  ###
function run_and_update_status!(
    cases_dict::OrderedDict;
    build_dict::Vector{Any} = OrderedDict(),
    range::UnitRange{Int64} = 1:1780
)
    for (id, value) in collect(cases_dict)[range]
        build_errors = case_build_errors(
            value;
            build_dict = build_dict
        )
        value["build_errors"] = build_errors

        if (length(build_errors) == 0)
            sim_report = case_sim_result(
                value;
                cases_path = "./cases/semantic",
                output_path = "./cases/output"
            )
        else
            sim_report = Dict(
                "status" => "SKIPPED",
                "message" => "Model was not simulated because of build errors"
            )
        end

        value["result"] = sim_report
    end
end
### ###

export upload_cases, filter_cases, add_cases, update_results,
    case_build_errors, case_sim_result, run_and_update_status! # Metelkin

end #module

max_stars_repo_name: insysbio/sbml-heta-cases

max_issues_repo_name: insysbio/sbml-heta-cases

max_forks_repo_name: insysbio/sbml-heta-cases


text:
VERIFICATION_DIR = "../verification/"
EXPERIMENTS = ["stratified_couette_flow"]

for exp in EXPERIMENTS
    script_filepath = joinpath(VERIFICATION_DIR, exp, exp * ".jl")
    try
        include(script_filepath)
    catch err
        @error sprint(showerror, err)
    end
end

function run_stratified_couette_flow_verification(arch)
    simulate_stratified_couette_flow(Nxy=16, Nz=8, arch=arch, Ri=0.01, Ni=1, end_time=1e-5)
    return true  # We're just checking to make sure the script runs with no errors.
end

@testset "Verification" begin
    @info "Testing verification scripts..."

    for arch in archs
        @testset "Stratified Couette flow verification [$(typeof(arch))]" begin
            @info "  Testing stratified Couette flow verification [$(typeof(arch))]"
            @test run_stratified_couette_flow_verification(arch)
        end
    end
end


text:
# ECC protocol implementation in Julia
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# PyECC (https://github.com/amintos/PyECC) as Pythonic reference

importall Base
using Krypto

abstract ECCKey

type ECCPubKey <: ECCKey
end

type ECCPrivKey <: ECCKey
end


text:
typealias SCIP_Objsense Int8

SCIP_OBJSENSE_MAXIMIZE = int8(-1)
SCIP_OBJSENSE_MINIMIZE = int8(+1)

SCIP_OBJSENSE = Dict{SCIP_Objsense, String}()
SCIP_OBJSENSE[SCIP_OBJSENSE_MAXIMIZE] = "maximize"
SCIP_OBJSENSE[SCIP_OBJSENSE_MINIMIZE] = "minimize"