text:
# SPDX-License-Identifier: MIT

module EnumX

export @enumx

abstract type Enum{T} <: Base.Enum{T} end

@noinline panic(x) = throw(ArgumentError(x))

macro enumx(args...)
    return enumx(__module__, Any[args...])
end

function symbol_map end

function enumx(_module_, args)
    T = :T
    if length(args) > 1 && Meta.isexpr(args[1], :(=), 2) && args[1].args[1] === :T &&
       (args[1].args[2] isa Symbol || args[1].args[2] isa QuoteNode)
        T = args[1].args[2]
        T isa QuoteNode && (T = T.value)
        popfirst!(args) # drop T=...
    end
    name = popfirst!(args)
    if name isa Symbol
        modname = name
        baseT = Int32
    elseif Meta.isexpr(name, :(::), 2) && name.args[1] isa Symbol
        modname = name.args[1]
        baseT = Core.eval(_module_, name.args[2])
    else
        panic("invalid EnumX.@enumx type specification: $(name).")
    end
    name = modname
    if length(args) == 1 && Meta.isexpr(args[1], :block)
        syms = args[1].args
    else
        syms = args
    end
    name_value_map = Vector{Pair{Symbol, baseT}}()
    doc_entries = Vector{Pair{Symbol,Expr}}()
    next = zero(baseT)
    first = true
    for s in syms
        s isa LineNumberNode && continue
        # Handle doc expressions
        doc_expr = nothing
        if Meta.isexpr(s, :macrocall, 4) && s.args[1] isa GlobalRef && s.args[1].mod === Core &&
               s.args[1].name === Symbol("@doc")
            doc_expr = s
            s = s.args[4]
        end
        if s isa Symbol
            if !first && next == typemin(baseT)
                panic("value overflow for Enum $(modname): $(modname).$(s) = $(next).")
            end
            sym = s
        elseif Meta.isexpr(s, :(=), 2) && s.args[1] isa Symbol
            if s.args[2] isa Symbol &&
               (i = findfirst(x -> x.first === s.args[2], name_value_map); i !== nothing)
                @assert name_value_map[i].first === s.args[2]
                nx = name_value_map[i].second
            else
                nx = Core.eval(_module_, s.args[2])
            end
            if !(nx isa Integer && typemin(baseT) <= nx <= typemax(baseT))
                panic(
                    "invalid value for Enum $(modname){$(baseT)}: " *
                    "$(modname).$(s.args[1]) = $(repr(nx))."
                )
            end
            next = convert(baseT, nx)
            sym = s.args[1]
        else
            panic("invalid EnumX.@enumx entry: $(s)")
        end
        if sym === T
            panic("instance name $(modname).$(sym) reserved for the Enum typename.")
        end
        if (idx = findfirst(x -> x.first === sym, name_value_map); idx !== nothing)
            v = name_value_map[idx].second
            panic(
                "duplicate name for Enum $(modname): $(modname).$(sym) = $(next)," *
                " name already used for $(modname).$(sym) = $(v)."
            )
        end
        push!(name_value_map, sym => next)
        if doc_expr !== nothing
            # Replace the documented expression since it might be :(Apple = ...)
            doc_expr.args[4] = sym
            push!(doc_entries, sym => doc_expr)
        end

        next += oneunit(baseT)
        first = false
    end
    value_name_map = Dict{baseT,Symbol}(v => k for (k, v) in reverse(name_value_map))
    module_block = quote
        primitive type $(T) <: Enum{$(baseT)} $(sizeof(baseT) * 8) end
        let value_name_map = $(value_name_map)
            check_valid(x) = x in keys(value_name_map) ||
                throw(ArgumentError("invalid value for Enum $($(QuoteNode(modname))): $(x)."))
            global function $(esc(T))(x::Integer)
                check_valid(x)
                return Base.bitcast($(esc(T)), convert($(baseT), x))
            end
            Base.Enums.namemap(::Base.Type{$(esc(T))}) = value_name_map
            Base.Enums.instances(::Base.Type{$(esc(T))}) =
                ($([esc(k) for (k,v) in name_value_map]...),)
            EnumX.symbol_map(::Base.Type{$(esc(T))}) = $(name_value_map)
        end
    end
    for (k, v) in name_value_map
        push!(module_block.args,
            Expr(:const, Expr(:(=), esc(k), Expr(:call, esc(T), v)))
        )
    end
    for (_, v) in doc_entries
        push!(module_block.args, v)
    end
    # Document the module and the type
    mdoc = Expr(:block, Expr(:meta, :doc), esc(modname))
    # TODO: Attach to the type too?
    # Tdoc = Expr(:block, Expr(:meta, :doc), Expr(:., esc(modname), QuoteNode(T)))
    return Expr(:toplevel, Expr(:module, false, esc(modname), module_block), mdoc, #=Tdoc,=# nothing)
end

function Base.show(io::IO, ::MIME"text/plain", x::E) where E <: Enum
    iob = IOBuffer()
    ix = Integer(x)
    for (k, v) in symbol_map(E)
        if v == ix
            print(iob, "$(nameof(parentmodule(E))).$(k) = ")
        end
    end
    show(iob, ix)
    write(io, seekstart(iob))
    return nothing
end
function Base.show(io::IO, ::MIME"text/plain", ::Type{E}) where E <: Enum
    if !isconcretetype(E) # handle EnumX.Enum and EnumX.Enum{T}
        invoke(show, Tuple{IO, Type}, io, E)
        return
    end
    iob = IOBuffer()
    insts = Base.Enums.instances(E)
    n = length(insts)
    stringmap = Pair{String, Base.Enums.basetype(E)}[
        string("$(nameof(parentmodule(E))).", k) => v for (k, v) in symbol_map(E)
    ]
    mx = maximum(x -> textwidth(x.first), stringmap; init = 0)
    print(iob,
        "Enum type $(nameof(parentmodule(E))).$(nameof(E)) <: ",
        "Enum{$(Base.Enums.basetype(E))} with $(n) instance$(n == 1 ? "" : "s")$(n>0 ? ":" : "")"
    )
    for (k, v) in stringmap
        print(iob, "\n ", rpad(k, mx), " = ")
        show(iob, v)
    end
    write(io, seekstart(iob))
    return nothing
end

end # module EnumX


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 793190e2-fce9-11ea-0f94-8d23919bbda3
using Pkg, DrWatson

# ╔═╡ 7931d048-fce9-11ea-3644-cbdb925b4031
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap, GLM
	using StatisticalRethinking
end

# ╔═╡ 9670408e-fce1-11ea-1d03-51e8376f67e1
md"## Clip-05-15-17s.jl"

# ╔═╡ 32dfbbfe-81e2-11eb-29f8-754d6aebcb49
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
	scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
end;

# ╔═╡ 851c9784-81e2-11eb-01f0-073edd961499
stan5_3 = "
data {
  int N;
  vector[N] divorce_s;
  vector[N] marriage_s;
  vector[N] medianagemarriage_s;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + + bA * medianagemarriage_s + bM * marriage_s;
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  divorce_s ~ normal( mu , sigma );
}
";

# ╔═╡ 851d5eda-81e2-11eb-1174-4fc6ec7295b1
begin
	data = (N = size(df, 1), divorce_s = df.Divorce_s,
		marriage_s = df.Marriage_s, medianagemarriage_s = df.MedianAgeMarriage_s)
	init = (a = 0.0, bM = 0.0, bA = 1.0, sigma = 1.0)
	q5_3s, m5_3s, o5_3s = stan_quap("m5.3s", stan5_3; data, init);
	if !isnothing(q5_3s)
		quap5_3s_df = sample(q5_3s)
	end
	if !isnothing(m5_3s)
		post5_3s_df = read_samples(m5_3s, :dataframe)
		PRECIS(post5_3s_df)
	end
end

# ╔═╡ 952ba332-828b-11eb-2511-bff99eaa1236
md"##### Quadratic approximation:"

# ╔═╡ 851cc740-81e2-11eb-0af3-7f59f78c706e
# Rethinking results
rethinking_results = "
	   mean   sd  5.5% 94.5%
a      0.00 0.10 -0.16  0.16
bM    -0.07 0.15 -0.31  0.18
bA    -0.61 0.15 -0.85 -0.37
sigma  0.79 0.08  0.66  0.91
";

# ╔═╡ 2e11fb74-828b-11eb-012d-67e9445387d5
if !isnothing(q5_3s)
	PRECIS(quap5_3s_df)
end

# ╔═╡ 794c76b4-fce9-11ea-3d4b-cdb6ca10a383
if !isnothing(m5_3s)
	begin
		part5_3s = read_samples(m5_3s, :particles)
		N = size(df, 1)
		plot(xlab="Observed divorce", ylab="Predicted divorce",
			title="Posterior predictive plot")
		v = zeros(size(df, 1), 4);
		for i in 1:N
			mu = mean(part5_3s.bM) * df[i, :Marriage_s] +
				mean(part5_3s.bA) * df[i, :MedianAgeMarriage_s]
			if i == 13
				annotate!([(df[i, :Divorce_s]-0.05, mu,
					Plots.text("ID", 6, :red, :right))])
			end
			if i == 39
				annotate!([(df[i, :Divorce_s]-0.05, mu,
					Plots.text("RI", 6, :red, :right))])
			end
			scatter!([df[i, :Divorce_s]], [mu], color=:red)
			s = rand(Normal(mu, mean(part5_3s.sigma)), 1000)
			v[i, :] = [maximum(s), hpdi(s, alpha=0.11)[2],
				hpdi(s, alpha=0.11)[1], minimum(s)]
		end
		for i in 1:N
			plot!([df[i, :Divorce_s], df[i, :Divorce_s]], [v[i,1], v[i, 4]],
				color=:darkblue, leg=false)
			plot!([df[i, :Divorce_s], df[i, :Divorce_s]], [v[i,2], v[i, 3]],
				line=2, color=:black, leg=false)
		end
		df2 = DataFrame(
			:x => df.Divorce_s,
			:y => [mean(part5_3s.bM) * df[i, :Marriage_s] +
				mean(part5_3s.bA) * df[i, :MedianAgeMarriage_s] for i in 1:N]
		)
		m1 = lm(@formula(y ~ x), df2)
		x = -2.1:0.1:2.2
		y = coef(m1)[2] * x
		plot!(x, y, line=:dash, color=:red)

	end
end

# ╔═╡ 7957c32a-fce9-11ea-09a4-6fdd2e231a7d
md"## End of clip-05-15-17s.jl"

# ╔═╡ Cell order:
# ╟─9670408e-fce1-11ea-1d03-51e8376f67e1
# ╠═793190e2-fce9-11ea-0f94-8d23919bbda3
# ╠═7931d048-fce9-11ea-3644-cbdb925b4031
# ╠═32dfbbfe-81e2-11eb-29f8-754d6aebcb49
# ╠═851c9784-81e2-11eb-01f0-073edd961499
# ╠═851d5eda-81e2-11eb-1174-4fc6ec7295b1
# ╟─952ba332-828b-11eb-2511-bff99eaa1236
# ╠═851cc740-81e2-11eb-0af3-7f59f78c706e
# ╠═2e11fb74-828b-11eb-012d-67e9445387d5
# ╠═794c76b4-fce9-11ea-3d4b-cdb6ca10a383
# ╟─7957c32a-fce9-11ea-09a4-6fdd2e231a7d

ext: jl


text:
function all_spin_one_half(array::Array)
    any_spin_not_one_half = any(x->(x!=1 && x!= -1), array)
    !(any_spin_not_one_half)
end


text:
@testset "Time Marching" begin
    import PotentialFlow.Utils: centraldiff
    T = [1e-1, 1e-2, 1e-3]

    f! = (ẋ, x, t) -> ẋ .= [-x[2], x[1]]
    update = (x₊, x₋, ẋ, Δt) -> @. x₊ = x₋ + Δt*ẋ

    @testset "Forward Euler" begin
        err = map(T) do Δt
            x₋ = [1.0, 2.0]
            x₊ = similar(x₋)
            ẋ  = similar(x₋)

            for t in 0:Δt:5.0-Δt
                forward_euler!(x₊, x₋, t, Δt, f!, update, ẋ)
                x₊, x₋ = x₋, x₊
            end
            norm(x₋ - [cos(5)-2sin(5),2cos(5)+sin(5)])
        end

        @test minimum(centraldiff(log.(err))./centraldiff(log.(T))) > 0.9
    end

    @testset "RK4" begin
        err = map(T) do Δt
            x₋ = [1.0, 2.0]
            x₊ = similar(x₋)
            ẋ  = [similar(x₋) for k in 1:4]

            for t in 0:Δt:5.0-Δt
                rk4!(x₊, x₋, t, Δt, f!, update, ẋ)
                x₊, x₋ = x₋, x₊
            end
            norm(x₋ - [cos(5)-2sin(5),2cos(5)+sin(5)])
        end

        @test minimum(centraldiff(log.(err))./centraldiff(log.(T))) > 3.9
    end
end


text:
using Molecules
import Molecules: Atom

export Molecule, Atom

using Fermi.Options
using Fermi.PhysicalConstants: atomic_number, bohr_to_angstrom

using LinearAlgebra
using Formatting

import Base: show

"""
    Fermi.Molecule

Object storing information about a molecule (group of atoms).

# Fields

    atoms         Array with Fermi.Atom objects
    charge        Charge of the molecule
    multiplicity  Multiplicity ``(2Mₛ + 1)``
    Nα            Number of Alpha electrons
    Nβ            Number of Beta electrons
    Vnuc          Nuclear repulsion energy

# Examples:

A Molecule object can be created by providing the desired keyworded arguments:

    molstring:      A string representing the XYZ of the molecule
    unit:           Unit used for the distance between atoms (Bohr or Angstrom)
    charge:         Molecular charge
    multiplicity    Molecular multiplicity

Any argument not given explicitly will be read from the Options.
```
julia> Molecule()
Molecule:

O    1.209153654800    1.766411818900   -0.017161397200
H    2.198480007500    1.797710062700    0.012116171900
H    0.919788188200    2.458018557000    0.629793883200

Charge: 0   Multiplicity: 1
Nuclear repulsion:    8.8880641737

julia> Molecule(charge=2, multiplicity=3)
Molecule:

O    1.209153654800    1.766411818900   -0.017161397200
H    2.198480007500    1.797710062700    0.012116171900
H    0.919788188200    2.458018557000    0.629793883200

Charge: 2   Multiplicity: 3
Nuclear repulsion:    8.8880641737
```
"""
struct Molecule
    atoms::Vector{Atom}
    charge::Int
    multiplicity::Int
    Nα::Int
    Nβ::Int
    Vnuc::Float64
end

function Molecule(;
    molstring = Options.get("molstring"),
    unit = Symbol(Options.get("unit")),
    charge = Options.get("charge"),
    multiplicity = Options.get("multiplicity")
    )

    atoms = Molecules.parse_string(molstring, unit=unit)
    Molecule(atoms, charge, multiplicity)
end

function Molecule(atoms::Vector{T}, charge::Int, multiplicity::Int) where T <: Atom

    # Compute Nuclear repulsion
    Vnuc = Molecules.nuclear_repulsion(atoms)

    # Compute number of electrons
    nelec = -charge
    for i in eachindex(atoms)
        nelec += atoms[i].Z
    end

    # If the number of electrons turns out negative returns an error
    if nelec ≤ 0
        throw(FermiException("Invalid charge ($charge) for given molecule"))
    end

    # Mult = 2Ms + 1 thus the number of unpaired electrons (taken as α) is Mult-1 = 2Ms
    αexcess = multiplicity-1

    # The number of β electrons must be equal the number of doubly occupied orbitals (Ndo)
    # Ndo = (nelec - n_of_unpaired_elec)/2 this must be an integer
    if isodd(nelec) != isodd(αexcess)
        throw(FermiException("Incompatible charge $(charge) and multiplicity $(multiplicity)"))
    end

    Nβ = (nelec - αexcess)/2
    Nα = nelec - Nβ
    out =  Molecule(atoms, charge, multiplicity, Nα, Nβ, Vnuc)
    return out
end

"""
    Fermi.string_repr(M::Molecule)

Returns a nicely formatted string with all the molecule's information
"""
function string_repr(M::Molecule)
    out = ""
    out = out*format("Molecule:\n\n")
    out = out*format(Molecules.get_xyz(M.atoms))
    out = out*format("\n")
    out = out*format("\nCharge: {}   ", M.charge)
    out = out*format("Multiplicity: {}   \n", M.multiplicity)
    out = out*format("Nuclear repulsion: {:15.10f}", M.Vnuc)
    return out
end

# Pretty printing
function show(io::IO, ::MIME"text/plain", X::Molecule)
    print(io, string_repr(X))
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
#using LinearAlgebra
#using DelimitedFiles
#using Statistics
#using Dates
#using PyCall

struct unit_vector
    x
    y
    xy
end

struct geometry_twisted
    lattice::String
    sites
    inter_vector::unit_vector
    dimension::Int
    twist_angle::Float64
    inter_distance
    inter_alignment
end

function get_theta(m,r,a1,a2)
    theta=acos((3*m^2+3*m*r+r^2/2)/(3*m^2+3*m*r+r^2))
    A1=m*a2+(m+r)*a1
    A2=-(m+r)*a2+(2*m+r)*a1                                     # the a1,a2 in literature is the opposite definition from here
    n=3*m^2+3*m*r+r^2
    return theta, A1, A2, n
end

function get_parameter_for_theta(theta)
    m=1
    r=1
    theta2=acos((3*m^2+3*m*r+r^2/2)/(3*m^2+3*m*r+r^2))

    while theta2>theta
        m+=1
        theta2=acos((3*m^2+3*m*r+r^2/2)/(3*m^2+3*m*r+r^2))
    end
    m-=1
    theta2=acos((3*m^2+3*m*r+r^2/2)/(3*m^2+3*m*r+r^2))
    factor=theta2/theta

    return m, factor
end

function rotation_along_z_3D(theta)
    U=[[cos(theta),sin(theta),0] [-sin(theta),cos(theta),0] [0,0,1]]     #in this way U=anti-clock rotation of theta
    return U
end

function geometry_simple_3D_rotated(n,m,a1,a2,z,theta,r0)
    U=rotation_along_z_3D(theta)
    a1=U*a1
    a2=U*a2
    R0=[0,0,z]-(n-1)/2*a1-(m-1)/2*a2+r0                                  # s.t. the second layer is rotated w.r.t. [0,0]
    R=Array{Any}(undef,n*m)
    ii=1
    for i=0:n-1
        for j=0:m-1
            R[ii]=i*a1+j*a2+R0
            ii=ii+1
        end
    end
    return R,a1,a2
end

function geometry_bipartite_3D_rotated(n,m,a1,a2,a3,z,theta,r0)
    U=rotation_along_z_3D(theta)
    a1=U*a1
    a2=U*a2
    a3=U*a3
    R0=[0,0,z]-(n-1)/2*a1-(m-1)/2*a2+r0
    R=Array{Any}(undef,2*n*m)
    ii=1
    for i=0:n-1
        for j=0:m-1
            for k=0:1
                R[ii]=i*a1+j*a2+k*a3+R0
                ii=ii+1
            end
        end
    end
    return R,a1,a2
end

function nearest(site_a,site_b)
    return 0.1<norm(site_a-site_b)<1.1
end

function next_nearest(site_a,site_b)
    return 1.1<norm(site_a-site_b)<1.9
end

function next_next_nearest(site_a,site_b)
    return 1.9<norm(site_a-site_b)<2.1
end

function plot_R(R)
    plt.figure(figsize=(6,6),dpi=80)
    N=size(R,1)
    x=Array{Float64}(undef,Int(N))
    y=Array{Float64}(undef,Int(N))
    x2=Array{Float64}(undef,Int(N))
    y2=Array{Float64}(undef,Int(N))
    j1=1
    j2=1
    for i=1:N
        if R[i][3]==0
            x[j1]=R[i][1]
            y[j1]=R[i][2]
            j1+=1
        else
            x2[j2]=R[i][1]
            y2[j2]=R[i][2]
            j2+=1
        end
    end
    if j1>1
        x=x[1:j1-1]
        y=y[1:j1-1]
        plt.scatter(x,y,color="red",s=10)
    end
    if j2>1
        x2=x2[1:j2-1]
        y2=y2[1:j2-1]
        plt.scatter(x2,y2,color="blue",s=10)
    end
    plt.axis([-40,40,-40,40])
    plt.xlabel("")
    plt.ylabel("")
    plt.show()
end

function add_vacancy(R,vacant_site)
    deleteat!(R,vacant_site)
    return R
end

function merge_R(R,R2)
    M=size(R,1)
    N=size(R2,1)
    RR=Array{Any}(undef,M+N)
    for i=1:M
        RR[i]=R[i]
    end
    for j=1:N
        RR[N+j]=R2[j]
    end
    return R
end

function check_site_in_unit_cell(site,A1,A2)      #check whether site is in unitcell [-1/2,1/2]*A1+[-1/2,1/2]*A2, site can be 3D
    num=A1[1]*A2[2]-A1[2]*A2[1]
    num2=site[1]*A1[2]-site[2]*A1[1]
    num3=site[1]*A2[2]-site[2]*A2[1]
    alpha=num3/num
    beta=-num2/num                                # site=alpha*a1+beta*a2
    if -0.50<round(alpha,digits=10)<=0.50 && -0.50<round(beta,digits=10)<=0.50
        return true
    else
        return false
    end
end

function check_R(R_unitcell,inter_vector,R)
    n_inter=size(inter_vector,1)
    R=[round.(x,digits=3) for x in R]
    #println(R)
    N=size(R_unitcell,1)
    for ii=2:n_inter
        r=inter_vector[ii]
        for i=1:N
            R0=R_unitcell[i]+r
            R0=round.(R0,digits=3)
            #println(R0)
            if !(R0 in R)
                println(R0)
                #println(R)
                return false
            end
        end
    end

    return true
end

function plot_R_unitcell(R_unitcell,inter_vector)
    n_inter=size(inter_vector,1)
    N=size(R_unitcell,1)
    for ii=2:n_inter
        r=inter_vector[ii]
        for i=1:N
            push!(R_unitcell,R_unitcell[i]+r)
        end
    end
    plot_R(R_unitcell)
end

function get_reciprocal_vector(A1,A2)
    A3=[0,0,1]
    B1=2*pi*cross(A2,A3)/dot(A1,cross(A2,A3))
    B2=2*pi*cross(A1,A3)/dot(A2,cross(A1,A3))
    #b1=[B1[1],B1[2]]
    #b2=[B2[1],B2[2]]
    return B1,B2
end

function get_inter_vector_3D(n,m,a1,a2)
    inter_vector=[[0,0,0],n*a1,m*a2,m*a2-n*a1]
    return inter_vector
end

function get_lattice_3D_rotated(lattice,n_a1,n_a2,z,theta,r0,a)   #theta in units of 1, z: layer-index  #add confinement in unitcell, r0: parallel alignment of layer, a=lattice constant
    if lattice=="triangular"
        a1=[1,0,0]*a
        a2=[1/2,-sqrt(3)/2,0]*a
        #println(a1,a2)
        R,a1,a2=geometry_simple_3D_rotated(n_a1,n_a2,a1,a2,z,theta,r0)
    elseif lattice=="square"
        a1=[1,0,0]*a
        a2=[0,1,0]*a
        R,a1,a2=geometry_simple_3D_rotated(n_a1,n_a2,a1,a2,z,theta,r0)
    elseif lattice=="honeycomb"
        a1=[sqrt(3),0,0]*a
        a2=[sqrt(3)/2,-3/2,0]*a
        a3=[sqrt(3)/2,-1/2,0]*a
        R,a1,a2=geometry_bipartite_3D_rotated(n_a1,n_a2,a1,a2,a3,z,theta,r0)
    else
        println("invalid lattice")
    end

    return R,a1,a2
end

function get_twisted_lattice_2(lattice,n_a1,n_a2,m,r=1,d=1,r0=[0,0,0],mode="pi_flux")
    mode=="pi_flux" ? a=2 : a=1
    R,a1,a2=get_lattice_3D_rotated(lattice,n_a1,n_a2,0,0,[0,0,0],a)
    theta, A1, A2, n_sites=get_theta(m,r,a1,a2)
    R2,a1_2,a2_2=get_lattice_3D_rotated(lattice,n_a1,n_a2,d,theta,r0,a)
    R_unitcell,n1,n2=(mode=="pi_flux" ? get_unit_cell_2(R,R2,A1,A2,a1,a2,a1_2,a2_2,n_sites) : get_unit_cell(R,R2,A1,A2,n_sites))
    N=size(R_unitcell,1)

    @assert n1==n2 "number of sites on different layers mismatch",n1,n2
    @assert n1+n2==N "primitive cell too small"
    inter_vector=unit_vector(A1,A2,A1-A2)
    g=geometry_twisted(lattice, R_unitcell, inter_vector, 2, theta, d, r0)

    return g
end

function get_twisted_lattice_ribbon(lattice,n_a1,n_a2,m,r=1,d=1,r0=[0,0,0],mode="pi_flux",n_A1=1,n_A2=1)
    mode=="pi_flux" ? a=2 : a=1
    R,a1,a2=get_lattice_3D_rotated(lattice,n_a1,n_a2,0,0,[0,0,0],a)
    theta, A1, A2, n_sites=get_theta(m,r,a1,a2)
    R2,a1_2,a2_2=get_lattice_3D_rotated(lattice,n_a1,n_a2,d,theta,r0,a)
    R_unitcell,n1,n2=(mode=="pi_flux" ? get_unit_cell_2(R,R2,A1,A2,a1,a2,a1_2,a2_2,n_sites,n_A1,n_A2) : get_unit_cell(R,R2,A1,A2,n_sites,n_A1,n_A2))
    N=size(R_unitcell,1)

    @assert n1==n2 "number of sites on different layers mismatch",n1,n2
    @assert n1+n2==N "primitive cell too small"
    inter_vector=unit_vector(n_A1*A1,n_A2*A2,n_A1*A1-n_A2*A2)
    g=geometry_twisted(lattice, R_unitcell, inter_vector, 2, theta, d, r0)

    return g
end

function get_unit_cell_2(R,R2,A1,A2,a1,a2,a1_2,a2_2,n_sites,n_A1=1,n_A2=1)
    A1=n_A1*A1
    A2=n_A2*A2
    n_sites=n_sites*8*n_A1*n_A2   # for triangular, the factor is 2*4=8
    R_unitcell=Array{Any}(undef,n_sites)
    N=size(R,1)
    n1=0
    n2=0
    n=1
    for i=1:N
        if check_site_in_unit_cell(R[i],A1,A2)
            n1=n1+4
            R_unitcell[n]=R[i]
            R_unitcell[n+1]=R[i]+a1/2
            R_unitcell[n+2]=R[i]+a2/2
            R_unitcell[n+3]=R[i]+a1/2+a2/2
            n=n+4
        end
        if check_site_in_unit_cell(R2[i],A1,A2)
            n2=n2+4
            R_unitcell[n]=R2[i]
            R_unitcell[n+1]=R2[i]+a1_2/2
            R_unitcell[n+2]=R2[i]+a2_2/2
            R_unitcell[n+3]=R2[i]+a1_2/2+a2_2/2
            n=n+4
        end
    end

    return R_unitcell,n1,n2
end

function get_unit_cell(R,R2,A1,A2,n_sites,n_A1=1,n_A2=1)
    A1=n_A1*A1
    A2=n_A2*A2
    n_sites=n_sites*4*n_A1*n_A2 # 4 for honeycomb (bipartite)
    R_unitcell=Array{Any}(undef,n_sites)
    N=size(R,1)
    n1=0
    n2=0
    n=1
    for i=1:N
        if check_site_in_unit_cell(R[i],A1,A2)
            n1=n1+1
            R_unitcell[n]=R[i]
            n=n+1
        end
        if check_site_in_unit_cell(R2[i],A1,A2)
            n2=n2+1
            R_unitcell[n]=R2[i]
            n=n+1
        end
    end

    return R_unitcell,n1,n2
end

function remove_layer_R(R,layer_index)
    N=size(R,1)
    for i=1:N
        j=N+1-i
        if R[j][3]==layer_index
            deleteat!(R,j)
        end
    end
    return R
end

function remove_layer_g(g,layer_index)
    R=remove_layer_R(g.sites,layer_index)
    return geometry_twisted(g.lattice,R,g.inter_vector,g.dimension,g.twist_angle,g.inter_distance,g.inter_alignment)
end

function test_geometry_twisted()
    n_a1=50
    n_a2=50
    m=1
    r=1
    r0=[0,0,0]
    d=1
    lattice,mode="triangular","pi_flux"
    #lattice,mode="honeycomb","uniform_hop"

    g=get_twisted_lattice_2(lattice,n_a1,n_a2,m,r,d,r0,mode)
    println(g.twist_angle*180/pi)
    println(size(g.sites,1))
    #println(typeof(g))
    remove_layer_g(g,1)
    println(size(g.sites,1))
    plot_R(g.sites)
end

#end  #end module


text: @inline A_mul_Bc(tvec::TransposedVector, B::Union{LinAlg.QRCompactWYQ,LinAlg.QRPackedQ}) = ctranspose(B*ctranspose(tvec))


text:
struct HessianNull <: Hessian end
struct HessianD21{H1,H2} <: Hessian
    h1::H1
    h2::H2
    ref::MyRef{Float64}
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
end
struct Hessian11{F,H1,H11} <: Hessian
    h1::H1
    h11::H11
    fref::MyRef{Float64}
    ref::MyRef{Float64}

    function Hessian11(e::Expression1{F,E},d, indexer = nothing) where {F,E}
        h1 = Hessian(e.e1,d.d1,indexer)
        h11= Hessian(d.d1,d.d1,indexer)
        return new{F,typeof(h1),typeof(h11)}(h1,h11,ref(e.e1),ref(d))
    end
end
struct Hessian11F1{F,H1,H11,R} <: Hessian
    a::R
    h1::H1
    h11::H11
    fref::MyRef{Float64}
    ref::MyRef{Float64}
    function Hessian11F1(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F,E1 <: Real,E2}
        h1 = Hessian(e.e2,d.d1,indexer)
        h11= Hessian(d.d1,d.d1,indexer)
        return new{F,typeof(h1),typeof(h11),typeof(e.e1)}(e.e1,h1,h11,ref(e.e2),ref(d))
    end
end
struct Hessian11F2{F,H1,H11,R} <: Hessian
    a::R
    h1::H1
    h11::H11
    fref::MyRef{Float64}
    ref::MyRef{Float64}
    function Hessian11F2(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F,E1,E2 <: Real}
        h1 = Hessian(e.e1,d.d1,indexer)
        h11= Hessian(d.d1,d.d1,indexer)
        return new{F,typeof(h1),typeof(h11),typeof(e.e2)}(e.e2,h1,h11,ref(e.e1),ref(d))
    end
end
struct Hessian02{H11,H12,H21,H22} <: Hessian
    h11::H11
    h12::H12
    h21::H21
    h22::H22
    ref11::MyRef{Float64}
    ref12::MyRef{Float64}
    ref21::MyRef{Float64}
    ref22::MyRef{Float64}
end
struct Hessian22{F,H1,H2,H11,H12,H21,H22} <: Hessian
    h1::H1
    h2::H2
    h11::H11
    h12::H12
    h21::H21
    h22::H22
    fref1::MyRef{Float64}
    fref2::MyRef{Float64}
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}

    function Hessian22(e::Expression2{F,E1,E2},d,indexer = nothing) where {F,E1,E2}
        h1 = Hessian(e.e1,d.d1,indexer)
        h2 = Hessian(e.e2,d.d2,indexer)
        h11= Hessian(d.d1,d.d1,indexer)
        h12= Hessian(d.d1,d.d2,indexer)
        h21= Hessian(d.d2,d.d1,indexer)
        h22= Hessian(d.d2,d.d2,indexer)
        new{F,typeof(h1),typeof(h2),typeof(h11),typeof(h12),typeof(h21),typeof(h22)}(h1,h2,h11,h12,h21,h22,ref(e.e1),ref(e.e2),ref1(d),ref2(d))
    end
end
# struct Hessian20{H1,H2} <: Hessian
#     h1::H1
#     h2::H2
#     ref1::MyRef{Float64}
#     ref2::MyRef{Float64}
# end
struct HessianD00 <: Hessian
    index1::Int
    index2::Int
    islower::Bool
end
struct HessianD00S <: Hessian
    index::Int
    islower::Bool
end
struct HessianD20{H1,H2} <: Hessian
    h1::H1
    h2::H2
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
end
struct HessianD10{H} <: Hessian
    h::H
    ref::MyRef{Float64}
end
struct HessianD11{H} <: Hessian
    h::H
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
end
struct Hessian22a{H1,H2} <: Hessian
    h1::H1
    h2::H2
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
end
struct Hessian22m{H1,H2,H12,H21} <: Hessian
    h1::H1
    h2::H2
    h12::H12
    h21::H21
    ref1::MyRef{Float64}
    ref2::MyRef{Float64}
end
struct Hessian11a{H} <: Hessian
    h1::H
    ref::MyRef{Float64}
end
struct HessianSum{I,H} <: Hessian
    inner::I
    hs::Vector{H}
end

const HESSIAN_NULL = HessianNull()

@inline function (h::HessianSum{I,H})(z,x,p=nothing,h0=1) where {I,H}
    inner(h)(z,x,p,h0)
    @simd for i in eachindex(h.hs)
        @inbounds h.hs[i](z,x,p,h0)
    end
end
@inline function (h::HessianSum{Nothing,H})(z,x,p=nothing,h0=1) where H
    @simd for i in eachindex(h.hs)
        @inbounds h.hs[i](z,x,p,h0)
    end
end

@inline (::HessianNull)(z,x,p=nothing,h0=1) = nothing
@inline function (h::HessianD00)(z,x,p=nothing,h0 = 1)
    islower(h) && @inbounds z[index1(h)::Int,index2(h)::Int] += h0
    return
end
@inline function (h::HessianD00S)(z,x,p=nothing,h0 = 1)
    islower(h) && @inbounds z[index(h)::Int] += h0
    return
end
@inline function (h::HessianD10{H})(z,x,p=nothing,h0 = 1) where H
    h.h(z,x,p,h0*refval(h))
    return
end
@inline function (h::HessianD11{H})(z,x,p=nothing,h0 = 1) where H
    h.h(z,x,p,h0*refval1(h)*refval2(h))
    return
end
@inline function (h::HessianD21{H1,H2})(z,x,p=nothing,h0=1) where {H1,H2}
    h.h1(z,x,p,h0*refval(h)*refval1(h))
    h.h2(z,x,p,h0*refval(h)*refval2(h))
    return
end
@inline function (h::HessianD20{H1,H2})(z,x,p=nothing,h0=1) where {H1,H2}
    h.h1(z,x,p,h0*refval1(h))
    h.h2(z,x,p,h0*refval2(h))
end
@inline function (h::Hessian02{H11,H12,H21,H22})(z,x,p=nothing,h0=1) where {H11,H12,H21,H22}
    h.h11(z,x,p,h0*refval11(h)*refval21(h))
    h.h12(z,x,p,h0*refval11(h)*refval22(h))
    h.h21(z,x,p,h0*refval12(h)*refval21(h))
    h.h22(z,x,p,h0*refval12(h)*refval22(h))
    return
end
@inline function (h::Hessian11a{H})(z,x,p=nothing,h0=1) where {H}
    h.h1(z,x,p,h0*refval(h))
    return
end
@inline function (h::Hessian22a{H1,H2})(z,x,p=nothing,h0=1) where {H1,H2}
    h.h1(z,x,p,h0*refval1(h))
    h.h2(z,x,p,h0*refval2(h))
    return
end
@inline function (h::Hessian22m{H1,H2,H12})(z,x,p=nothing,h0=1) where {H1,H2,H12}
    h.h1(z,x,p,h0*refval1(h))
    h.h2(z,x,p,h0*refval2(h))
    h.h12(z,x,p,h0)
    h.h21(z,x,p,h0)
    return
end

Hessian(e::ExpressionSum{E,I1},d::GradientSum{D,I2},indexer = nothing) where {E,D,I1,I2} = HessianSum(Hessian(inner(e),inner(d),indexer),[Hessian(e,d,indexer) for (e,d) in zip(e.es,d.ds)])
Hessian(e::ExpressionSum{E,Nothing},d::GradientSum{D,Nothing},indexer = nothing) where {E,D} = HessianSum(nothing,[Hessian(e,d,indexer) for (e,d) in zip(e.es,d.ds)])
Hessian(e::Variable,::G,indexer) where G <: Gradient = HESSIAN_NULL
Hessian(e::Parameter,::G,indexer) where G <: Gradient = HESSIAN_NULL
Hessian(e::Constant,::G,indexer) where G <: Gradient = HESSIAN_NULL
Hessian(d1::G,d2::GradientNull,indexer) where G <: Gradient = HESSIAN_NULL
Hessian(d1::GradientNull,d2::G,indexer) where G <: Gradient = HESSIAN_NULL
Hessian(d1::GradientNull,d2::GradientNull,indexer) = HESSIAN_NULL
Hessian(d1::Gradient0,d2::Gradient0,indexer)= HessianD00S(index(d1)>=index(d2) ? set_indexer!(indexer,index(d1),index(d2)) : 0,index(d1) >= index(d2))
Hessian(d1::Gradient0,d2::Gradient0,::Nothing) = HessianD00(index(d1),index(d2),index(d1) >= index(d2))
Hessian(d1::Gradient0,d2::G1, indexer = nothing) where G1 <: Union{Gradient1,Gradient2F1,Gradient2F2} = HessianD10(Hessian(d1,d2.d1,indexer),ref(d2))
Hessian(d1::G1,d2::Gradient0, indexer = nothing) where G1 <: Union{Gradient1,Gradient2F1,Gradient2F2} = HessianD10(Hessian(d1.d1,d2,indexer),ref(d1))
Hessian(d1::G1,d2::G2, indexer = nothing) where {G1 <: Union{Gradient1,Gradient2F1,Gradient2F2}, G2 <: Union{Gradient1,Gradient2F1,Gradient2F2}}= HessianD11(Hessian(d1.d1,d2.d1,indexer),ref(d1),ref(d2))
Hessian(d1::G1,d2::G2, indexer = nothing) where {G1 <: Union{Gradient1,Gradient2F1,Gradient2F2}, G2 <: Gradient2} = HessianD21(Hessian(d1.d1,d2.d1,indexer),Hessian(d1.d1,d2.d2,indexer),ref(d1),ref1(d2),ref2(d2))
Hessian(d1::G1,d2::G2, indexer = nothing) where {G1 <: Gradient2, G2 <: Union{Gradient1,Gradient2F1,Gradient2F2}} = HessianD21(Hessian(d1.d1,d2.d1,indexer),Hessian(d1.d2,d2.d1,indexer),ref(d2),ref1(d1),ref2(d1))
Hessian(d1::Gradient0,d2::Gradient2{F,F1,F2}, indexer = nothing) where {F,F1,F2} = HessianD20(Hessian(d1,d2.d1,indexer),Hessian(d1,d2.d2,indexer),ref1(d2),ref2(d2))
Hessian(d1::Gradient2{F,F1,F2},d2::Gradient0, indexer = nothing) where {F,F1,F2} = HessianD20(Hessian(d1.d1,d2,indexer),Hessian(d1.d2,d2,indexer),ref1(d1),ref2(d1))
Hessian(d1::Gradient2{F,F1,F2} where {F,F1,F2},d2::Gradient2{F,F1,F2} where {F,F1,F2}, indexer = nothing) = Hessian02(Hessian(d1.d1,d2.d1,indexer),Hessian(d1.d1,d2.d2,indexer),Hessian(d1.d2,d2.d1,indexer),Hessian(d1.d2,d2.d2,indexer),ref1(d1),ref2(d1),ref1(d2),ref2(d2))
Hessian(e::Expression1{F,E},d, indexer = nothing) where {F,E} = Hessian11(e,d,indexer)
Hessian(e::Expression2{F,E1,E2},d,indexer = nothing) where {F,E1,E2} = Hessian22(e,d,indexer)
Hessian(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F,E1 <: Real,E2} = Hessian11F1(e,d,indexer)
Hessian(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F,E1,E2 <: Real} = Hessian11F2(e,d,indexer)
Hessian(e::Expression2{typeof(*),E1,E2},d,indexer = nothing) where {E1,E2} = Hessian22m(Hessian(e.e1,d.d1,indexer),Hessian(e.e2,d.d2,indexer),Hessian(d.d1,d.d2,indexer),Hessian(d.d2,d.d1,indexer),ref1(d),ref2(d))
Hessian(e::Expression2{F,E1,E2},d,indexer = nothing) where {F<:Union{typeof(+),typeof(-)},E1,E2} = Hessian22a(Hessian(e.e1,d.d1,indexer),Hessian(e.e2,d.d2,indexer),ref1(d),ref2(d))
Hessian(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F<:Union{typeof(+),typeof(-),typeof(*)},E1 <: Real,E2 <: Expression} = Hessian11a(Hessian(e.e2,d.d1,indexer),ref(d))
Hessian(e::Expression2{F,E1,E2}, d,indexer = nothing) where {F<:Union{typeof(+),typeof(-),typeof(*),typeof(/)},E1 <: Expression,E2 <: Real} = Hessian11a(Hessian(e.e1,d.d1,indexer),ref(d))

# performance killer ---------------
function Hessian(d1::GradientSum{D1,I1},d2::GradientSum{D2,I2},indexer = nothing) where {D1,D2,I1,I2}
    @warn "This operation is expensive"
    hinner = Hessian(inner(d1),d2,indexer)
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        hinner(z,x,p,h0)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::GradientSum{D1,Nothing},d2::GradientSum{D2,I2},indexer = nothing) where {D1,D2,I2}
    @warn "This operation is expensive"
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::GradientSum{D1,I1},d2::GradientSum{D2,Nothing},indexer = nothing) where {D1,D2,I1}
    @warn "This operation is expensive"
    hinner = Hessian(inner(d1),d2,indexer)
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        hinner(z,x,p,h0)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::GradientSum{D1,Nothing},d2::GradientSum{D2,Nothing},indexer = nothing) where {D1,D2}
    @warn "This operation is expensive"
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::GradientSum{D1,I1},d2::G,indexer = nothing) where {D1,D2,I1,I2,G <: Gradient}
    @warn "This operation is expensive"
    hinner = Hessian(inner(d1),d2,indexer)
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        hinner(z,x,p,h0)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::GradientSum{D1,Nothing},d2::G,indexer = nothing) where {D1,D2,I2,G <: Gradient}
    @warn "This operation is expensive"
    hs = [Hessian(d,d2,indexer) for d in d1.ds]
    @inline function (z,x,p=nothing,h0=1)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::G,d2::GradientSum{D2,I2},indexer = nothing) where {G <: Gradient,D2,I2}
    @warn "This operation is expensive"
    hinner = Hessian(d1,inner(d2),indexer)
    hs = [Hessian(d1,d,indexer) for d in d2.ds]
    @inline function (z,x,p=nothing,h0=1)
        hinner(z,x,p,h0)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end
function Hessian(d1::G,d2::GradientSum{D2,Nothing},indexer = nothing) where {G <: Gradient, D2}
    @warn "This operation is expensive"
    hs = [Hessian(d1,d,indexer) for d in d2.ds]
    @inline function (z,x,p=nothing,h0=1)
        @simd for i in eachindex(hs)
            @inbounds hs[i](z,x,p,h0)
        end
    end
end

islower(h) = h.islower


text:
function intersect_box_box(	lbx1::Float64,
							ubx1::Float64,
							lby1::Float64,
							uby1::Float64,
							lbz1::Float64,
							ubz1::Float64,
							lbx2::Float64,
							ubx2::Float64,
							lby2::Float64,
							uby2::Float64,
							lbz2::Float64,
							ubz2::Float64,
							Lx::Float64,
							Ly::Float64,
							Lz::Float64)

	if 0.5 * (lbx1 + ubx1) - 0.5 * (lbx2 + ubx2) < -0.5*Lx
		lbx1 = lbx1 + Lx
		ubx1 = ubx1 + Lx
	elseif 0.5 * (lbx1 + ubx1) - 0.5 * (lbx2 + ubx2) > 0.5*Lx
		lbx1 = lbx1 - Lx
		ubx1 = ubx1 - Lx
	end

	if 0.5 * (lby1 + uby1) - 0.5 * (lby2 + uby2) < -0.5*Ly
		lby1 = lby1 + Ly
		uby1 = uby1 + Ly
	elseif 0.5 * (lby1 + uby1) - 0.5 * (lby2 + uby2) > 0.5*Ly
		lby1 = lby1 - Ly
		uby1 = uby1 - Ly
	end

	if 0.5 * (lbz1 + ubz1) - 0.5 * (lbz2 + ubz2) < -0.5*Lz
		lbz1 = lbz1 + Lz
		ubz1 = ubz1 + Lz
	elseif 0.5 * (lbz1 + ubz1) - 0.5 * (lbz2 + ubz2) > 0.5*Lz
		lbz1 = lbz1 - Lz
		ubz1 = ubz1 - Lz
	end

	return 	(abs( (lbx1 + ubx1) - (lbx2 + ubx2)) <= ubx1 - lbx1 + ubx2 - lbx2 ) &&
			(abs( (lby1 + uby1) - (lby2 + uby2)) <= uby1 - lby1 + uby2 - lby2 ) &&
			(abs( (lbz1 + ubz1) - (lbz2 + ubz2)) <= ubz1 - lbz1 + ubz2 - lbz2 )
end


text:
using Makie
using AbstractPlotting: textslider

function main()
    sa, _a = textslider(range(-1.,stop=1.,length=100),"a", start=1.)
    sb, _b = textslider(range(-1.,stop=1.,length=100),"b", start=0.)
    xs = -π:0.01:π
    N = length(xs)
    ab = lift(tuple,_a,_b)
    pos = lift(ab) do (a,b)
        map(LinRange(-2π, 2π, 100)) do x
            Point2f0(x,sin(a*x+b))
        end
    end

    p = lines(pos)

    RecordEvents(
        hbox(p, vbox(sa,sb), parent=Scene()),
        "output"
    )
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
struct CostFunctionUE <: CostFunction
    tfc::TimeFunctionContainer

    function CostFunctionUE(network::AbstractNetwork, fn::Function)
        new(TimeFunctionContainer(network, fn))
    end
end

function (f::CostFunctionUE)(x::Array{<:Real,1}, ids=nothing; returnitems=[:costs], tolls=nothing)
    t, dt = nothing, nothing

    if (:costs ∈ returnitems) || (:tsc ∈ returnitems)
        t = times(f.tfc, x, ids)
    end
    if (:derivs ∈ returnitems)
        dt = dtimes(f.tfc, x, ids)
    end
    tolls = (tolls == nothing) ? zero(x) : (ids == nothing) ? tolls : tolls[ids]

    rv = Dict()
    if :costs ∈ returnitems
        rv[:costs] = t + tolls
    end
    if :derivs ∈ returnitems
        rv[:derivs] = dt
    end
    if :tsc ∈ returnitems
        rv[:tsc] = sum(x .* t)
    end

    return NamedTuple(rv)
end

struct CostFunctionSO <: CostFunction
    tfc::TimeFunctionContainer

    function CostFunctionSO(network::AbstractNetwork, fn::Function)
        new(TimeFunctionContainer(network, fn))
    end
end

function (f::CostFunctionSO)(x::Array{<:Real,1}, ids=nothing; returnitems=[:costs], tolls=nothing)
    t, dt, ddt = nothing, nothing, nothing

    if (:costs ∈ returnitems) || (:tsc ∈ returnitems)
        t = times(f.tfc, x, ids)
    end
    if (:costs ∈ returnitems) || (:derivs ∈ returnitems)
        dt = dtimes(f.tfc, x, ids)
    end
    if (:derivs ∈ returnitems)
        ddt = ddtimes(f.tfc, x, ids)
    end
    tolls = (tolls == nothing) ? zero(x) : (ids == nothing) ? tolls : tolls[ids]

    rv = Dict()
    if :costs ∈ returnitems
        rv[:costs] = t + x .* dt
    end
    if :derivs ∈ returnitems
        rv[:derivs] = 2 * dt + x .* ddt
    end
    if :tsc ∈ returnitems
        rv[:tsc] = sum(x .* t)
    end
    return NamedTuple(rv)
end

(f::CostFunction)(x::AbstractMatrix, args...; kwargs...) = f(sum(x, dims=2)[:,1], args...; kwargs...)


text:
using DataStreams
using DataFrames

filename = Pkg.dir( "ExcelDataStreams", "test", "test.xlsx" )

columnindices = [0,2]
names = [:date, :Fibonacci]
eltypes = [Date, Int]
skipstart = 2

ed = ExcelDataStreams.ExcelDataStream(
    filename,
    columnindices = columnindices,
    names = names,
    eltypes = eltypes,
    skipstart = skipstart,
)

df = Data.close!( Data.stream!( ed, DataFrame ) )
assert( size(df)==(9,2) )
assert( unique(diff(df[:date])) == [Dates.Day(1)] )
assert( all(df[3:end,:Fibonacci] .== df[1:end-2,:Fibonacci] + df[2:end-1,:Fibonacci]) )


text:
using NamedTuples
using Base.Test

@test @NT( a => 1 ).a == 1
@test @NT( a => 1 )[1] == 1
@test @NT( a => 1 )[:a] == 1

@test @NT( :a => 1 ).a == 1
@test @NT( :a => 1 )[1] == 1
@test @NT( :a => 1 )[:a] == 1

@test length( @NT( a => 1)) == 1
@test length( @NT( a => 1, b => 2.0)) == 2

@test first( @NT( a => 1, b => 2.0 )) == 1
@test last( @NT( a => 1, b => "hello", c => 2.0 )) == 2.0
@test [ v for v in @NT( a => 1.0, b => 2.0 ) ] == [ 1.0, 2.0 ]

@test ( x = @NT( a::Int64, b::Float64 )( 1, 2.0 ) ; typeof(x.a) == Int64 && typeof(x.b) == Float64 )
@test @NT( a => 1, b => "hello")  ==  @NT( a, b )( 1, "hello")
@test @NT( a => 1) != @NT( b => 1 )

@test hash( @NT( a => 1, b => "hello"))  ==  hash( @NT( a, b )( 1, "hello") )
@test hash( @NT( a => 1, b => "hello")) != hash( @NT( a => 1, b => 2.0 ))

@test @NT( a ) ==  @NT( a )
@test @NT( a ) !=  @NT( b )

@test typeof( @NT( a::Int64, b::Float64 )(1, 3.0) ) == typeof( @NT( a => 1, b => 2.0 ))

# Syntax tests, including anon named tuples
@test @NT( a, b ) <: NamedTuple
@test @NT( ::Int64, ::Float64 ) <: NamedTuple
@test typeof( @NT( 1, 2, "3" )) <: NamedTuple
@test typeof( @NT( 1 + 2, "hello")) <: NamedTuple
@test length( @NT( 1, 2, 3, "hello")[ 1:2 ]) == 2

@test isbits( @NT( ::Int64, ::Float64)) == true

nt = @NT( a=>1, b=>2, c=>3 )
@test nt.a == 1
@test nt.b == 2
@test nt.c == 3
@test haskey( nt, :x ) == false
x = setindex( nt, :x, 123 )
@test x.x == 123
@test x.a == 1
@test x.b == 2
@test x.c == 3
@test nt[[:c,:b]] == @NT( c=>3, b=>2 )

y = delete( x, :a)
@test x != y
@test haskey( nt, :a ) == true
@test haskey( y, :a ) == false
@test x.x == 123
@test x.a == 1
@test x.b == 2
@test x.c == 3

@test merge( nt, @NT( d => "hello", e => "world"))  == @NT( a=>1,b=>2,c=>3,d=>"hello",e=>"world")


text:
@json struct PipelineSchedule
    id::Int
    description::String
    ref::String
    cron::String
    cron_timezone::String
    next_run_at::DateTime
    active::Bool
    created_at::DateTime
    updated_at::DateTime
    owner::User
end

## Pipeline Schedules

# Get All Pipeline Schedules
# https://docs.gitlab.com/ee/api/pipeline_schedules.html#get-all-pipeline-schedules
function endpoint(
    ::GitLabAPI,
    ::typeof(list_pipeline_schedules),
    project::Integer,
)
    return Endpoint(:GET, "/projects/$project/pipeline_schedules")
end

function endpoint(
    ::GitLabAPI,
    ::typeof(list_pipeline_schedules),
    owner::AStr,
    repo::AStr,
)
    return Endpoint(:GET, "/projects/$(encode(owner, repo))/pipeline_schedules")
end

into(::GitLabAPI, ::typeof(list_pipeline_schedules)) = Vector{PipelineSchedule}


text:
facts("Search space") do
  context("isinspace") do
    for(i in 1:NumTestRepetitions)
      reps = rand(1:10)
      ss1 = symmetric_search_space(reps, (0.0, 1.0))
      ind = rand_individual(ss1)
      for(j in 1:reps)
        @fact (mins(ss1)[j] <= ind[j] <= maxs(ss1)[j])  => true
      end
    end
  end

  context("Symmetric search space with default range") do
    ss1 = symmetric_search_space(1)
    @fact numdims(ss1) => 1
    @fact ranges(ss1) => [(0.0, 1.0)]
    @fact range_for_dim(ss1,1) => (0.0, 1.0)

    for(i in 1:NumTestRepetitions)
      ind = rand_individual(ss1)
      @fact size(ind) => (1, 1)
      @fact isinspace(ind, ss1) => true
    end

    ss3 = symmetric_search_space(3)
    @fact numdims(ss3) => 3
    @fact ranges(ss3) => [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]
    @fact range_for_dim(ss3,1) => (0.0, 1.0)
    @fact range_for_dim(ss3,2) => (0.0, 1.0)
    @fact range_for_dim(ss3,3) => (0.0, 1.0)

    for(i in 1:NumTestRepetitions)
      ind = rand_individual(ss3)
      @fact size(ind) => (3, 1)
      @fact isinspace(ind, ss3) => true
    end
  end

  context("SymmetricSearchSpace with given range") do
    ss1 = symmetric_search_space(1, (-1.0, 1.0))
    @fact numdims(ss1) => 1
    @fact ranges(ss1) => [(-1.0, 1.0)]
    @fact range_for_dim(ss1,1) => (-1.0, 1.0)

    for(i in 1:NumTestRepetitions)
      reps = rand(1:100)
      range = (rand(), rand())
      ss = symmetric_search_space(reps, range)
      @fact numdims(ss) => reps
      @fact all([(dr == range) for dr in ranges(ss)]) => true
    end
  end

  context("rand_individual is within the search space") do
    for(i in 1:NumTestRepetitions)
      reps = rand(1:100)
      mm = sort(rand(2,1), 1)
      range = (mm[1], mm[2])
      ss = symmetric_search_space(reps, range)
      ind = rand_individual(ss)
      @fact length(ind) => numdims(ss)
      @fact isinspace(ind, ss) => true
    end
  end

  context("rand_individuals creates many individuals and all are within the search space") do
    for(i in 1:NumTestRepetitions)
      reps = rand(1:10)
      mm = sort(rand(2,1), 1)
      range = (mm[1], mm[2])
      ss = symmetric_search_space(reps, range)
      numinds = rand(1:10)
      inds = rand_individuals(ss, numinds)
      @fact size(inds,1) => numdims(ss)
      @fact size(inds,2) => numinds
      for(j in 1:numinds)
        @fact isinspace(inds[j], ss) => true
      end
    end
  end

  context("RangePerDimSearchSpace") do
    ss = RangePerDimSearchSpace([(0.0, 1.0)])
    @fact mins(ss) => hcat([0.0])
    @fact maxs(ss) => hcat([1.0])
    @fact deltas(ss) => hcat([1.0])

    ss = RangePerDimSearchSpace([(0.0, 1.0), (0.5, 10.0)])
    @fact mins(ss) => [0.0, 0.5]''
    @fact maxs(ss) => [1.0, 10.0]''
    @fact deltas(ss) => [1.0, 9.5]''
  end

  context("rand_individuals_lhs samples in LHS intervals") do
    ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

    inds = rand_individuals_lhs(ss, 2)
    @fact size(inds, 1) => 2
    @fact size(inds, 2) => 3

    sorted = sort(inds, 1) # Sort per column => in their ordered intervals
    @fact 0.0 <= sorted[1,1] <= 0.5 => true
    @fact 0.5 <= sorted[2,1] <= 1.0 => true

    @fact 2.0 <= sorted[1,2] <= 2.5 => true
    @fact 2.5 <= sorted[2,2] <= 3.0 => true

    @fact 4.0 <= sorted[1,3] <= 4.5 => true
    @fact 4.5 <= sorted[2,3] <= 5.0 => true
  end

  context("feasible finds feasible points in the search space") do
    ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

    # We use the double transpose below to ensure the actual and expected
    # values have the same type (matrices, not vectors).
    @fact BlackBoxOptim.feasible([1.1, 2.0, 4.0], ss) == [1.0, 2.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([1.1, 3.0, 4.0], ss) == [1.0, 3.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([1.1, 2.0, 5.0], ss) == [1.0, 2.0, 5.0]'' => true
    @fact BlackBoxOptim.feasible([1.1, 3.0, 5.0], ss) == [1.0, 3.0, 5.0]'' => true

    @fact BlackBoxOptim.feasible([-0.1, 2.0, 4.0], ss) == [0.0, 2.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([-0.1, 3.0, 4.0], ss) == [0.0, 3.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([-0.1, 2.0, 5.0], ss) == [0.0, 2.0, 5.0]'' => true
    @fact BlackBoxOptim.feasible([-0.1, 3.0, 5.0], ss) == [0.0, 3.0, 5.0]'' => true

    @fact BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]'' => true
    @fact BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]'' => true

    @fact BlackBoxOptim.feasible([0.0, 3.3, 4.0], ss) == [0.0, 3.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([0.0, 3.2, 4.0], ss) == [0.0, 3.0, 4.0]'' => true
    @fact BlackBoxOptim.feasible([1.0, 3.1, 5.0], ss) == [1.0, 3.0, 5.0]'' => true
    @fact BlackBoxOptim.feasible([1.0, 3.9, 5.0], ss) == [1.0, 3.0, 5.0]'' => true

    @fact BlackBoxOptim.feasible([-0.4, 3.3, 14.5], ss) == [0.0, 3.0, 5.0]'' => true
  end

  context("diameters") do
    ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])
    diams = diameters(ss)

    @fact length(diams) => 3
    @fact diams[:] => [1.0, 1.0, 1.0]
  end

end


text:
# # Basic Regression
# This tutorial is a port of [hello_world.py](https://github.com/aamini/evidential-deep-learning/blob/d1d8e395fb083308d14fa92c5ce766e97b2a066a/hello_world.py) to julia.
# We create a synthetic dataset, train a simple evidential regression model and plot the results
using Flux, Plots
import EvidentialDeepLearning as EDL
using Distributions
# Lets create some data
function my_data(x_min, x_max, n; train=true)
    x = range(Float32(x_min), Float32(x_max), length=n)
    x = reshape(x, (1,n))
    σ = train ? 3f0 : 0f0
    y = x .^3 .+ rand(Normal(0f0, σ), size(x))
    return x, y
end
x_train, y_train = my_data(-4, 4, 1000)
x_test, y_test = my_data(-7, 7, 1000, train=false)
plot(title="Dataset", ylims=(-150,150), xlims=(-7,7))
scatter!(vec(x_train), vec(y_train), label="train", markersize=1)
plot!(vec(x_test), vec(y_test), label="truth")

# Our model is a basic single hidden layer dense network, where the output is lifted to the evidential world:
model = Flux.Chain(
    Dense(1, 64, Flux.relu),
    Dense(64, 4, identity),
    EDL.NIGRegressor(),
)
loss(x, y) = EDL.regression_loss(model(x), y, λ=1f-2)
# Lets train the model
data = Flux.Data.DataLoader((x_train, y_train), batchsize=100, shuffle=true)
opt = ADAM(5e-4) # training is quite sensitive to learning rate
for epoch in 1:500 # takes ~ 30s
    Flux.train!(loss, Flux.params(model), data, opt)
end
# Lets plot the results.
plot(legend=:bottom)
h = model(x_test)
preds = vec(EDL.predict.(h))
plot!(vec(x_test), preds, label="Prediction", color=:blue, xlims=(-7,7), ylims=(-150,150))
σ_pp = vec(EDL.std_predict.(h)) # std of the predictive posterior
lo = preds - σ_pp
hi = preds + σ_pp
plot!(vec(x_test), hi, fillrange=lo, label="Predicted error", color=:green, alpha=0.5)
scatter!(vec(x_train), vec(y_train), label="training data", markersize=1, color=:red)
plot!(vec(x_test), vec(y_test), label="truth", color=:black)


text:
module RTM_64x200

    using JUDI.TimeModeling, SegyIO, LinearAlgebra, PyPlot, JUDI.SLIM_optim, Statistics, Random
    using ImageFiltering, Random, DelimitedFiles

    export rtm_isic, plot_rtm

    function rtm_isic(m, m0)

        m = adjoint(m)
        m0 = adjoint(m0)

        n = size(m)
        dm = vec(m - m0)

        # Setup info and model structure
        d = (25., 25.)
        o = (0., 0.)
        nsrc = 21	# number of sources

        ## Set up receiver geometry in meters
        nxrec = 201
        xrec = range(0f0, stop=5000f0, length=nxrec)
        yrec = 0f0
        zrec = range(0f0, stop=0f0, length=nxrec)

        # receiver sampling and recording time
        timeR = 2000f0   # receiver recording time [ms]
        dtR = 1f0    # receiver sampling interval [ms]

        # Set up receiver structure
        recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

        ## Set up source geometry (cell array with source locations for each shot)
        xsrc = convertToCell(range(0f0, stop=5000f0, length=nsrc))
        ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
        zsrc = convertToCell(range(0f0, stop=0f0, length=nsrc))

        # source sampling and number of time steps
        timeS = timeR  # ms
        dtS = dtR   # ms

        # Set up source structure
        srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

        # setup wavelet
        f0 = 0.006f0     # KHz
        wavelet = ricker_wavelet(timeS, dtS, f0)
        q = judiVector(srcGeometry, wavelet);

        # Set up model and info structure for linear operators
        model = Model(n, d, o, m; nb=200)
        model0 = Model(n, d, o, m0; nb=200)
        ntComp = get_computational_nt(srcGeometry, recGeometry, model)
        info = Info(prod(n), nsrc, ntComp)

        opt = Options(space_order=16, isic=true) # Use linearized inverse scattering imaging condition for the Jacobian

        # Setup operators
        Pr = judiProjection(info, recGeometry)
        F = judiModeling(info, model; options=opt)
        F0 = judiModeling(info, model0; options=opt)
        Ps = judiProjection(info, srcGeometry)
        J = judiJacobian(Pr*F0*adjoint(Ps), q) # Born forward modeling operator

        ####################################################################################################
        ## Forward modeling and RTM
        d_obs = Pr*F*adjoint(Ps)*q  # Practical observation data generated on m
        # d_syn = Pr*F0*adjoint(Ps)*q # Synthetic observation data generated on m0

        # When the migarion velocity is close to the true velocity
        idx_wb = 9 # index of the water bottom
        m1 = ones(Float32, n) .* m[1,1] # constant velocity
        m1[:, idx_wb+1:end] .= m[:, idx_wb+1]
        model1 = Model(n, d, o, m1; nb=200)
        F1 = judiModeling(info, model1; options=opt)
        d_syn = Pr*F1*adjoint(Ps)*q # forward modeling on the constant velocity to mute the direct wave

        rtm = adjoint(J) * (d_obs - d_syn)

        # Alternatively using fwi_objective() or lsrtm_objective()
        # fval, rtm = fwi_objective(model0, q, d_obs; options=opt)
        # x = zeros(Float32, model0.n)
        # fval, rtm = lsrtm_objective(model0, q, d_obs, Mr*x; nlind=true, options=opt) # Mr is a preconditioner like model topmute

        return rtm
    end

    """
    plot_rtm(image, spacing)

    Plots an rtm image. Image needs to be (x, z).

    # Arguments

    - `image`: RTM image shape (nx, nz)
    - `spacing`: Grid spacing (dx, dz)
    - `perc`: perc for vmin/vmax, default at 98
    - `cmap`: Colormap, default is "PuOr"
    - `o`: Origin of image, default is (0, 0)
    - `interp`: matplotlib interplolation, default is "hanning"
    - `aspect`: Aspect ratio of plot, defaults to "auto"
    - `d_scale`: Depth scaling power, depth^d_scale, default is 1.5
    - `name`: Plot title, default is "RTM"
    - `units`: Unit of spacing, default is "m"
    - `new_fig`: Whether to create a new figure, default is true
    - `save`: Whether to save the figure, default name is "name" if just true, or save if save is a string
    """

    function plot_rtm(image, spacing; perc=98, cmap="PuOr", o=(0, 0),
                    interp="hanning", aspect=nothing, d_scale=1.5,
                    name="RTM", units="m", new_fig=true, save=nothing)
        nx, nz = size(image)
        dx, dz = spacing
        ox, oz = o
        depth = range(oz, oz + (nz - 1)*spacing[2], length=nz).^d_scale

        scaled = image .* depth'

        a = quantile(abs.(vec(scaled)), perc/100)
        extent = [ox, ox+ (nx-1)*dx, oz+(nz-1)*dz, oz]
        isnothing(aspect) && (aspect = .5 * nx/nz)

        if new_fig
            figure()
        end
        imshow(scaled', vmin=-a, vmax=a, cmap=cmap, aspect=aspect,
            interpolation=interp, extent=extent)
        xlabel("X ($units)")
        ylabel("Depth ($units)")
        title("$name")

        if ~isnothing(save)
            save == True ? filename=name : filename=save
            savefig(filename, bbox_inches="tight", dpi=150)
        end
    end

end


text:
function event_occurred(obj::ODEEventHandler, arg0::ODEStateAndDerivative, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (ODEStateAndDerivative, jboolean), arg0, arg1)
end

function g(obj::ODEEventHandler, arg0::ODEStateAndDerivative)
    return jcall(obj, "g", jdouble, (ODEStateAndDerivative,), arg0)
end

function init(obj::ODEEventHandler, arg0::ODEStateAndDerivative, arg1::jdouble)
    return jcall(obj, "init", void, (ODEStateAndDerivative, jdouble), arg0, arg1)
end

function reset_state(obj::ODEEventHandler, arg0::ODEStateAndDerivative)
    return jcall(obj, "resetState", ODEState, (ODEStateAndDerivative,), arg0)
end


text:
using Contour, GLVisualize, GeometryTypes, GLAbstraction, Colors, FileIO

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0.0, 1.0, 360))
end

description = """
This example uses the Contour.jl library to generate
contour lines from an heightfield.
"""

# create a rotation from the time signal
rotation = map(timesignal) do t
    rotationmatrix_z(Float32(t*2pi)) # -> 4x4 Float32 rotation matrix
end

xrange = -5f0:0.02f0:5f0
yrange = -5f0:0.02f0:5f0

z = Float32[sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x) for x in xrange, y in yrange]
mini = minimum(z)
maxi = maximum(z)
color_ramp = map(x->RGBA{Float32}(x, 1.0), colormap("Blues"))
points = Point3f0[]
colors = RGBA{Float32}[]

for h in mini:0.2f0:maxi
    c = Contour.contour(xrange, yrange, z, h)
    for elem in c.lines
        c = color_lookup(color_ramp, h, mini, maxi)
        for p in elem.vertices
            push!(points, Point3f0(p[1], p[2], h))
            push!(colors, c)
        end
        # we can seperate colors with NaN's (might change soon!)
        push!(points, Point3f0(NaN32))
        push!(colors, c)
    end
end

line_renderable = visualize(
    points, :lines,
    color = colors,
    model = rotation
)
_view(line_renderable, window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using NQCModels.DiabaticModels: LargeDiabaticModel

struct LargeDiabaticCalculator{T,M} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::Hermitian{T,Matrix{T}}
    derivative::Matrix{Hermitian{T,Matrix{T}}}
    eigen::LinearAlgebra.Eigen{T,T,Matrix{T},Vector{T}}
    adiabatic_derivative::Matrix{Matrix{T}}
    nonadiabatic_coupling::Matrix{Matrix{T}}
    tmp_mat::Matrix{T}
    tmp_mat_complex1::Matrix{Complex{T}}
    tmp_mat_complex2::Matrix{Complex{T}}
    function LargeDiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        n = nstates(model)
        potential = Hermitian(zeros(n, n))
        derivative = [Hermitian(zeros(n, n)) for i=1:ndofs(model), j=1:atoms]
        eigen = Eigen(zeros(n), zeros(n,n)+I)
        adiabatic_derivative = [zeros(n, n) for i=1:ndofs(model), j=1:atoms]
        nonadiabatic_coupling = [zeros(n, n) for i=1:ndofs(model), j=1:atoms]
        tmp_mat = zeros(T, n, n)
        tmp_mat_complex1 = zeros(Complex{T}, n, n)
        tmp_mat_complex2 = zeros(Complex{T}, n, n)
        new{T,M}(model, potential, derivative, eigen, adiabatic_derivative, nonadiabatic_coupling,
            tmp_mat, tmp_mat_complex1, tmp_mat_complex2)
    end
end

function Calculator(model::LargeDiabaticModel, atoms::Integer, T::Type=Float64)
    LargeDiabaticCalculator{T}(model, atoms)
end

function evaluate_potential!(calc::LargeDiabaticCalculator, R)
    NQCModels.potential!(calc.model, calc.potential, R)
end

function evaluate_potential!(calc::LargeDiabaticCalculator, R::AbstractArray{T,3}) where {T}
    @views for i in axes(R, 3)
        NQCModels.potential!(calc.model, calc.potential[i], R[:,:,i])
    end
end

function transform_derivative!(calc::LargeDiabaticCalculator)
    for I in eachindex(calc.derivative)
        LinearAlgebra.mul!(calc.tmp_mat, calc.derivative[I], calc.eigen.vectors)
        LinearAlgebra.mul!(calc.adiabatic_derivative[I], calc.eigen.vectors', calc.tmp_mat)
    end
end

function eigen!(calc::LargeDiabaticCalculator)
    eig = LinearAlgebra.eigen(calc.potential)
    correct_phase!(eig, calc.eigen.vectors)
    copyto!(calc.eigen.vectors, eig.vectors)
    copyto!(calc.eigen.values, eig.values)
end

function correct_phase!(eig::LinearAlgebra.Eigen, old_eigenvectors::AbstractMatrix)
    @views for i=1:length(eig.values)
        if LinearAlgebra.dot(eig.vectors[:,i], old_eigenvectors[:,i]) < 0
            eig.vectors[:,i] .*= -1
        end
    end
end

function evaluate_nonadiabatic_coupling!(calc::LargeDiabaticCalculator)
    for I in eachindex(calc.adiabatic_derivative)
        evaluate_nonadiabatic_coupling!(calc.nonadiabatic_coupling[I], calc.adiabatic_derivative[I], calc.eigen.values)
    end
end

function evaluate_nonadiabatic_coupling!(coupling::Matrix, adiabatic_derivative::Matrix, eigenvalues::Vector)
    for i=1:length(eigenvalues)
        for j=i+1:length(eigenvalues)
            coupling[j,i] = -adiabatic_derivative[j,i] / (eigenvalues[j]-eigenvalues[i])
            coupling[i,j] = -coupling[j,i]
        end
    end
end


text:
mutable struct KalmanFilter <: AbstractEstimator
    belief_::MvNormal{Float64}
    motion_noise_stds::Dict{String,Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    pose_::Vector{Float64}
    map_::Map
    function KalmanFilter(
        envmap::Map,
        initial_pose::Vector{Float64},
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05;
        glob = false, # GlobalKf
        xlim = [-5.0, 5.0],
        ylim = [-5.0, 5.0],
    )
        @assert length(initial_pose) == 3
        if glob
            pose_distrib = PoseUniform(xlim, ylim)
            rand_pose = uniform(pose_distrib)
            new(
                MvNormal(rand_pose, Diagonal([1e+4, 1e+4, 1e+4])),
                motion_noise_stds,
                distance_dev_rate,
                direction_dev,
                copy(initial_pose),
                envmap,
            )
        else
            new(
                MvNormal(initial_pose, Diagonal([1e-10, 1e-10, 1e-10])),
                motion_noise_stds,
                distance_dev_rate,
                direction_dev,
                copy(initial_pose),
                envmap,
            )
        end
    end
end

function matM(v::Float64, ω::Float64, dt::Float64, stds::Dict{String,Float64})
    return Diagonal([
        stds["vv"]^2 * abs(v) / dt + stds["vω"]^2 * abs(ω) / dt,
        stds["ωv"]^2 * abs(v) / dt + stds["ωω"]^2 * abs(ω) / dt,
    ])
end

function matA(v::Float64, ω::Float64, dt::Float64, θ::Float64)
    st, ct = sin(θ), cos(θ)
    stw, ctw = sin(θ + ω * dt), cos(θ + ω * dt)
    a11 = (stw - st) / ω
    a12 = -v / (ω^2) * (stw - st) + v / ω * dt * ctw
    a21 = (-ctw + ct) / ω
    a22 = -v / (ω^2) * (-ctw + ct) + v / ω * dt * stw
    a31 = 0.0
    a32 = dt
    return [
        a11 a12
        a21 a22
        a31 a32
    ]
end

function matF(v::Float64, ω::Float64, dt::Float64, θ::Float64)
    F = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    F[1, 3] = v / ω * (cos(θ + ω * dt) - cos(θ))
    F[2, 3] = v / ω * (sin(θ + ω * dt) - sin(θ))
    return F
end

function matH(pose::Vector{Float64}, landmark_pos::Vector{Float64})
    mx, my = landmark_pos[1], landmark_pos[2]
    mux, muy, muθ = pose[1], pose[2], pose[3]
    q = (mux - mx)^2 + (muy - my)^2
    q11 = (mux - mx) / sqrt(q)
    q12 = (muy - my) / sqrt(q)
    q21 = (my - muy) / q
    q22 = (mux - mx) / q
    return [
        q11 q12 0.0
        q21 q22 -1.0
    ]
end

function matQ(distance_dev::Float64, direction_dev::Float64)
    return Diagonal([distance_dev^2, direction_dev^2])
end

function motion_update(kf::KalmanFilter, v::Float64, ω::Float64, dt::Float64; kwargs...)
    if abs(ω) < 1e-5
        ω = 1e-5
    end
    μ = mean(kf.belief_)
    Σ = cov(kf.belief_)
    M = matM(v, ω, dt, kf.motion_noise_stds)
    A = matA(v, ω, dt, μ[3])
    F = matF(v, ω, dt, μ[3])
    cov_mat = F * Σ * transpose(F) + A * M * transpose(A)
    cov_mat = (cov_mat + transpose(cov_mat)) / 2.0
    kf.belief_ = MvNormal(state_transition(μ, v, ω, dt), cov_mat)
    kf.pose_ = mean(kf.belief_)
end

function observation_update(
    kf::KalmanFilter,
    observation::Vector{Vector{Float64}},
    envmap = Map();
    kwargs...,
)
    μ = mean(kf.belief_)
    Σ = cov(kf.belief_)
    envmap = kf.map_
    for obsv in observation
        z = obsv[1:2]
        obs_index = convert(Int64, obsv[3]) + 1 # landmark id starts from 0!

        H = matH(μ, envmap[obs_index].pos_)
        estimated_z = observation_function(μ, envmap[obs_index].pos_)
        Q = matQ(estimated_z[1] * kf.distance_dev_rate, kf.direction_dev)
        K = Σ * transpose(H) * inv(Q + H * Σ * transpose(H))
        μ += K * (z - estimated_z)
        Σ = (Matrix(1.0I, 3, 3) - K * H) * Σ
    end
    Σ = (Σ + transpose(Σ)) / 2.0
    kf.belief_ = MvNormal(μ, Σ)
    kf.pose_ = μ
end

function draw(kf::KalmanFilter, p::Plot{T}) where {T}
    pose = kf.pose_
    cov_mat = cov(kf.belief_)
    p = covellipse!(pose[1:2], cov_mat[1:2, 1:2], n_std = 3, aspect_ratio = 1)

    x, y, θ = pose[1], pose[2], pose[3]
    sigma3 = sqrt(cov_mat[3, 3]) * 3.0
    xs = [x + cos(θ - sigma3), x, x + cos(θ + sigma3)]
    ys = [y + sin(θ - sigma3), y, y + sin(θ + sigma3)]
    p = plot!(xs, ys, color = "blue", alpha = 0.5)
    annota = "($(round(pose[1], sigdigits=3)), $(round(pose[2], sigdigits=3)), $(round(pose[3], sigdigits=3)))"
    p = annotate!(pose[1] + 1.0, pose[2] + 1.0, text(annota, 10))
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
"""
A context structure for holding state related to serializing Julia
objects. A key component is an `IOBuffer` used to hold the serialized
result.
"""
struct SerializeContext
    io::IOBuffer
    store::Dict{Any,Any}   # Meant to map Julia object to variable name
    init::Vector{Any}      # Expressions to run initially
end
SerializeContext(io::IOBuffer = IOBuffer()) = SerializeContext(io, Dict(), Vector{Expr}())

const _td = IdDict(
    Any => :jl_any_type,
    Float64 => :jl_float64_type,
    Float32 => :jl_float32_type,
    Int64 => :jl_int64_type,
    Int32 => :jl_int32_type,
    Int16 => :jl_int16_type,
    Int8 => :jl_int8_type,
    UInt64 => :jl_uint64_type,
    UInt32 => :jl_uint32_type,
    UInt16 => :jl_uint16_type,
    UInt8 => :jl_uint8_type,
    Cint => :jl_int32_type,
    Cvoid => :jl_any_type,
    Array => :jl_array_type,
    Array{Any,1} => :jl_array_any_type,
    Array{Int32,1} => :jl_array_int32_type,
    Array{UInt8,1} => :jl_array_uint8_type,
    ErrorException => :jl_errorexception_type,
    DataType => :jl_datatype_type,
    UnionAll => :jl_unionall_type,
    Union => :jl_union_type,
    Core.TypeofBottom => :jl_typeofbottom_type,
    TypeVar => :jl_tvar_type,
)

const _t = IdDict()

for (t,s) in _td
    _t[t] = :(unsafe_load(cglobal($(QuoteNode(s)), Type)))
end

const _gd = IdDict(
    Core => :jl_core_module,
    Main => :jl_main_module,
    nothing => :jl_nothing,
    () => :jl_emptytuple,
    Core.svec() => :jl_emptysvec,
    UndefRefError() => :jl_undefref_exception,
)

const _g = IdDict()

for (x,s) in _gd
    _g[x] = :(unsafe_load(cglobal($(QuoteNode(s)), Any)))
end

"""
    serialize(ctx::SerializeContext, x)

Serialize `x` into the context object `ctx`. `ctx.io` is the `IOBuffer` where the
serialized results are stored. Get the result with `take!(ctx.io)`.

This function returns an expression that will deserialize the object. Several `serialize`
methods can be called recursively to build up deserialization code for nested objects.
The expression returned is meant to be `eval`ed into a function that can be called
to do the serialization.

The deserialization code should be pretty low-level code that can be compiled
relatively easily. It especially shouldn't use global variables.

Serialization / deserialization code can use `ctx` to hold state information.

Some simple types like boxed variables do not need to write anything to `ctx.io`.
They can return an expression that directly creates the object.
"""
function serialize(ctx::SerializeContext, @nospecialize(x))
    haskey(_g, x) && return _g[x]
    # TODO: fix this major kludge.
    if nfields(x) > 0
        return Expr(:tuple, (serialize(ctx, getfield(x,i)) for i in 1:nfields(x))...)
    end
    return :(unsafe_load(cglobal(:jl_emptytuple, Any)))
end

function serialize(ctx::SerializeContext, @nospecialize(t::DataType))
    if haskey(_t, t)
        return _t[t]
    elseif haskey(ctx.store, t)
        return ctx.store[t]
    else
        # primary = unwrap_unionall(t.wrapper)
        name = gensym(Symbol(:type, "-", t.name.name))
        ctx.store[t] = name
        e = quote
            $name = let
                local tn    = $(serialize(ctx, t.name))
                # names = $(serialize(ctx, t.names))
                local super = $(serialize(ctx, t.super))
                local parameters = $(serialize(ctx, t.parameters))
                local types = $(serialize(ctx, t.types))
                local ndt = ccall(:jl_new_datatype, Any,
                            (Any, Any, Any, Any, Any, Any, Cint, Cint, Cint),
                            tn, tn.module, super, parameters, #=names=# unsafe_load(cglobal(:jl_any_type, Any)), types,
                            $(t.abstract), $(t.mutable), $(t.ninitialized))
                # tn.wrapper = ndt.name.wrapper
                # ccall(:jl_set_const, Cvoid, (Any, Any, Any), tn.module, tn.name, tn.wrapper)
                ndt
                # ty = tn.wrapper
                # $(ctx.types[string(t)]) = ndt
                # hasinstance = serialize(ctx, )
                # $(if isdefined(primary, :instance) && !isdefined(t, :instance)
                #     # use setfield! directly to avoid `fieldtype` lowering expecting to see a Singleton object already on ty
                #     :(Core.setfield!(ty, :instance, ccall(:jl_new_struct, Any, (Any, Any...), ty)))
                # end)
            end
        end
        push!(ctx.init, e)
        return name
    end
end

function serialize(ctx::SerializeContext, tn::Core.TypeName)
    haskey(ctx.store, tn) && return ctx.store[tn]
    name = gensym(Symbol(:typename, "-", tn.name))
    ctx.store[tn] = name
    e = quote
        $name = ccall(:jl_new_typename_in, Ref{Core.TypeName}, (Any, Any),
            #   $(serialize(ctx, tn.name)), Main  #=__deserialized_types__ =# )
              $(serialize(ctx, tn.name)), unsafe_load(cglobal(:jl_main_module, Any))  #=__deserialized_types__ =# )
    end
    push!(ctx.init, e)
    return name
end

function serialize(ctx::SerializeContext, mi::Core.MethodInstance)
    return :(unsafe_load(cglobal(:jl_emptytuple, Any)))
end

function serialize(ctx::SerializeContext, x::String)
    advance!(ctx.io)
    v = Vector{UInt8}(x)
    ioptr = ctx.io.ptr
    write(ctx.io, v)
    quote
        unsafe_string(Vptr + $(ioptr - 1), $(length(v)))
    end
end

function serialize(ctx::SerializeContext, x::Symbol)
    haskey(ctx.store, x) && return ctx.store[x]
    name = gensym(Symbol(:symbol, "-", x))
    ctx.store[x] = name
    e = quote
        $name = ccall(:jl_symbol_n, Any, (Ptr{UInt8}, Csize_t), $(serialize(ctx, string(x))), $(length(string(x))))
        # ccall(:jl_set_global, Cvoid, (Any, Any, Any), unsafe_load(cglobal(:jl_main_module, Any)), $(QuoteNode(name)), x)
    end
    push!(ctx.init, e)
    return name
end

# Define functions that return an expression. Example:
#    serialize(ctx::SerializeContext, x::Int) = :(ccall(:jl_box_int64, Any, (Int,), $x))
for (fun, type) in (:jl_box_int64 => Int64,     :jl_box_int32 => Int32,    :jl_box_int8 => Int16,    :jl_box_int8 => Int8,
                    :jl_box_uint64 => UInt64,   :jl_box_uint32 => UInt32,  :jl_box_uint8 => UInt16,  :jl_box_uint8 => UInt8,
                    :jl_box_voidpointer => Ptr{Cvoid},
                    :jl_box_float64 => Float64, :jl_box_float32 => Float32)
    @eval serialize(ctx::SerializeContext, x::$type) = Expr(:call, :ccall, QuoteNode($(QuoteNode(fun))), Any, Expr(:tuple, $type), x)
end
serialize(ctx::SerializeContext, x::Char) = :(ccall(:jl_box_char, Any, (UInt32,), $x))
serialize(ctx::SerializeContext, x::Bool) = :(ccall(:jl_box_bool, Any, (UInt8,), $x))

function serialize(ctx::SerializeContext, a::Tuple)
    length(a) == 0 && return :(unsafe_load(cglobal(:jl_emptytuple, Any)))
    Expr(:tuple, (serialize(ctx, x) for x in a)...)
end

function serialize(ctx::SerializeContext, a::Core.SimpleVector)
    length(a) == 0 && return :(unsafe_load(cglobal(:jl_emptysvec, Any)))
    Expr(:call, Expr(:., :Core, QuoteNode(:svec)), (serialize(ctx, x) for x in a)...)
end

advance!(io) = write(io, repeat('\0', -rem(io.ptr - 1, 8, RoundUp)))  # Align data to 8 bytes

function serialize(ctx::SerializeContext, a::Array{T,N}) where {T,N}
    elty = eltype(a)
    aty = typeof(a)
    dims = size(a)
    atys = serialize(ctx, aty)
    if isbitstype(elty)
        advance!(ctx.io)
        ioptr = ctx.io.ptr
        write(ctx.io, a)
        if N == 1
            advance!(ctx.io)
            ioptr = ctx.io.ptr
            write(ctx.io, a)
            quote
                p = Vptr + $ioptr - 1
                ccall(:jl_ptr_to_array_1d, $aty, (Any, Ptr{Cvoid}, Csize_t, Cint), $atys, p, $(length(a)), false)
            end
        else
            dms = serialize(ctx, dims)
            advance!(ctx.io)
            ioptr = ctx.io.ptr
            write(ctx.io, a)
            quote
                p = Vptr + $ioptr - 1
                ccall(:jl_ptr_to_array, $aty, (Any, Ptr{Cvoid}, Any, Int32), $atys, p, $dms, false)
            end
        end
    else
        idx = Int[]
        e = Array{Any}(undef, length(a))
        @inbounds for i in eachindex(a)
            if isassigned(a, i)
                e[i] = serialize(ctx, a[i])
                push!(idx, i)
            end
        end
        aname = gensym()
        resulte = [quote
            # $aname = Array{$elty, $(length(dims))}(undef, $dims)
            $aname = ccall(:jl_new_array, $aty, (Any, Any), $atys, $(serialize(ctx, dims)))
        end]
        for i in idx
            push!(resulte, quote
                # unsafe_store!(pointer($aname), $(e[i]), $i)
                unsafe_store!(convert(Ptr{Any}, pointer($aname)), $(e[i]), $i)
                # unsafe_store!(convert(Ptr{Csize_t}, pointer($aname)), pointer_from_objref($(e[i])), $i)
                # @inbounds $aname[$i] = $(e[i])
            end)
        end
        push!(resulte, :($aname = $aname))
        Expr(:block, resulte...)
    end
end


text:
using REPL
import REPL: Terminals
using RCall

mutable struct FakeTerminal <: Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    FakeTerminal(stdin,stdout,stderr,hascolor=true) =
        new(stdin,stdout,stderr,hascolor,false)
end

Terminals.hascolor(t::FakeTerminal) = t.hascolor
Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
Terminals.size(t::FakeTerminal) = (24, 80)

# fake repl

input = Pipe()
output = Pipe()
err = Pipe()
Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
Base.link_pipe!(output, reader_supports_async=true, writer_supports_async=true)
Base.link_pipe!(err, reader_supports_async=true, writer_supports_async=true)

repl = REPL.LineEditREPL(FakeTerminal(input.out, output.in, err.in), true)

repltask = @async begin
    REPL.run_repl(repl)
end

send_repl(x, enter=true) = write(input, enter ? "$x\n" : x)

function read_repl(io::IO, x)
    cache = Ref{Any}("")
    read_task = @task cache[] = readuntil(io, x)
    t = Base.Timer((_) -> Base.throwto(read_task,
                ErrorException("Expect \"$x\", but wait too long.")), 5)
    schedule(read_task)
    fetch(read_task)
    close(t)
    cache[]
end

check_repl_stdout(x) = length(read_repl(output, x)) > 0
check_repl_stderr(x) = length(read_repl(err, x)) > 0

# waiting for the repl
send_repl("using RCall")

RCall.RPrompt.repl_init(repl)

send_repl("\$", false)
@test check_repl_stdout("R> ")

send_repl("bar = 'apple'")
@test check_repl_stdout("bar = 'apple'")

send_repl("paste0('pine', bar)")
@test check_repl_stdout("pineapple")

send_repl("mtca\t", false)
@test check_repl_stdout("mtcars")

send_repl("")
@test check_repl_stdout("\n")

send_repl("\\alp\t", false)
@test check_repl_stdout("\\alpha")

send_repl("\t", false)
@test check_repl_stdout("α")

send_repl("")
@test check_repl_stdout("\n")

send_repl("foo]")
@test check_repl_stderr("unexpected")

send_repl("stop('something is wrong')")
@test check_repl_stderr("something is wrong")

send_repl("not_found")
@test check_repl_stderr("not found")

send_repl("\b", false)
@test check_repl_stdout("julia> ")

send_repl("x = \"orange\"")
@test check_repl_stdout("x = \"orange\"")

send_repl("\$", false)
@test check_repl_stdout("R> ")

send_repl("\$x")
@test check_repl_stdout("orange")

send_repl("\$not_found")
@test check_repl_stderr("UndefVarError")

# check visibility

send_repl("'apple'")
send_repl("paste0('check', 'point')")
@test occursin("\"apple\"", read_repl(output, "checkpoint"))

send_repl("invisible('apple')")
send_repl("paste0('check', 'point')")
@test !occursin("\"apple\"", read_repl(output, "checkpoint"))


text:
using Images, FileIO

h, w = 50, 70
img = zeros(RGB{N0f8}, h, w)
img[10:40, 5:35] = colorant"skyblue"
for i in 26:50, j in (i-25):40
    img[i, j] = colorant"sienna1"
end

save("data/bitmapWrite.ppm", img)
save("data/bitmapWrite.png", img)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
macro goldensectiontrace()
    quote
        if tracing
            dt = Dict()
            if extended_trace
                dt["x_minimum"] = x_minimum
                dt["x_lower"] = x_lower
                dt["x_upper"] = x_upper
            end
            update!(tr,
                    iteration,
                    f_minimum,
                    NaN,
                    dt,
                    store_trace,
                    show_trace,
                    show_every,
                    callback)
        end
    end
end

immutable GoldenSection <: Optimizer end

function optimize{T <: AbstractFloat}(f::Function, x_lower::T, x_upper::T,
                                      mo::GoldenSection;
                                      rel_tol::T = sqrt(eps(T)),
                                      abs_tol::T = eps(T),
                                      iterations::Integer = 1_000,
                                      store_trace::Bool = false,
                                      show_trace::Bool = false,
                                      callback = nothing,
                                      show_every = 1,
                                      extended_trace::Bool = false,
                                      nargs...)
    if !(x_lower < x_upper)
        error("x_lower must be less than x_upper")
    end

    # Save for later
    initial_lower = x_lower
    initial_upper = x_upper

    const golden_ratio::T = 0.5 * (3.0 - sqrt(5.0))

    x_minimum = x_lower + golden_ratio*(x_upper-x_lower)
    f_minimum = f(x_minimum)
    f_calls = 1 # Number of calls to f

    iteration = 0
    converged = false

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = store_trace || show_trace || extended_trace || callback != nothing
    @goldensectiontrace

    while iteration < iterations

        tolx = rel_tol * abs(x_minimum) + abs_tol

        x_midpoint = (x_upper+x_lower)/2

        if abs(x_minimum - x_midpoint) <= 2*tolx - (x_upper-x_lower)/2
            converged = true
            break
        end

        iteration += 1

        if x_upper - x_minimum > x_minimum - x_lower
            x_new = x_minimum + golden_ratio*(x_upper - x_minimum)
            f_new = f(x_new)
            f_calls += 1
            if f_new < f_minimum
                x_lower = x_minimum
                x_minimum = x_new
                f_minimum = f_new
            else
                x_upper = x_new
            end
        else
            x_new = x_minimum - golden_ratio*(x_minimum - x_lower)
            f_new = f(x_new)
            f_calls += 1
            if f_new < f_minimum
                x_upper = x_minimum
                x_minimum = x_new
                f_minimum = f_new
            else
                x_lower = x_new
            end
        end

        @goldensectiontrace
    end

    return UnivariateOptimizationResults("Golden Section Search",
                                         initial_lower,
                                         initial_upper,
                                         x_minimum,
                                         Float64(f_minimum),
                                         iteration,
                                         iteration == iterations,
                                         converged,
                                         rel_tol,
                                         abs_tol,
                                         tr,
                                         f_calls)
end


text:
using Test
using CliOptions

include("testutils.jl")

@testset "CounterOption()" begin
    @testset "ctor" begin
        @testset "names" begin
            @test_throws MethodError CounterOption()
            @test_throws MethodError CounterOption("-a", "--foo", "--bar")
            @test_throws ArgumentError CounterOption("")
            @test_throws ArgumentError CounterOption("a")
            @test_throws ArgumentError CounterOption("-")
            @test_throws ArgumentError CounterOption("--")
            @test_throws ArgumentError CounterOption(""; decrementers = "-a")
            @test_throws ArgumentError CounterOption("-a"; decrementers = "")
            @test_throws ArgumentError CounterOption("-a"; decrementers = [""])
            @test_throws ArgumentError CounterOption("-a"; decrementers = ["a"])
            @test_throws ArgumentError CounterOption("-a"; decrementers = ["-"])
            @test_throws ArgumentError CounterOption("-a"; decrementers = ["--"])
            @test_throws ArgumentError CounterOption(UInt8, "-a")

            option = CounterOption("-a")
            @test option.names == ("-a",)
            @test option.decrementers == String[]
        end

        @testset "decrementers" begin
            option = CounterOption("-a", "-b", decrementers = "-c")
            @test option.names == ("-a", "-b")
            @test option.decrementers == ["-c"]

            option = CounterOption("-a", "-b", decrementers = ["-c", "-d"])
            @test option.names == ("-a", "-b")
            @test option.decrementers == ["-c", "-d"]
        end

        @testset "default" begin
            let result = CliOptions.ParseResult()
                @test_throws InexactError CounterOption(Int8, "-v", default = -129)
                CounterOption(Int8, "-v", default = -128)
                CounterOption(Int8, "-v", default = 127)
                @test_throws InexactError CounterOption(Int8, "-v", default = 128)
            end
        end
    end

    @testset "ctor; duplicates, $(v[1]) and $(v[2])" for v in [
        (["-f", "--foo"], ["-b", "--bar"], (nothing, nothing)),
        (["-f", "-f"], ["-b", "--bar"], (ArgumentError, "-f")),
        (["-f", "--foo"], ["-f", "--bar"], (ArgumentError, "-f")),
        (["-f", "--foo"], ["-b", "--foo"], (ArgumentError, "--foo")),
        (["-f", "--foo"], ["-b", "-b"], (ArgumentError, "-b")),
    ]
        names, decrementers, expected = v
        if expected[1] isa Type
            tr = @test_throws expected[1] CounterOption(names...;
                                                        decrementers = decrementers)
            if tr isa Test.Pass
                buf = IOBuffer()
                showerror(buf, tr.value)
                msg = String(take!(buf))
                @test msg == ("ArgumentError: Duplicate names for a CounterOption found: " *
                              expected[2])
            end
        else
            @test CounterOption(names...; decrementers = decrementers) !== nothing
        end
    end

    @testset "show(x); $(join(v[1],','))" for v in [
        (["-a"], "CounterOption(:a)"),
        (["-a", "--foo-bar"], "CounterOption(:a,:foo_bar)"),
    ]
        names, expected_repr = v
        option = CounterOption(names...)
        @test repr(option) == expected_repr
    end

    @testset "show(io, x)" begin
        let option = CounterOption("-a")
            buf = IOBuffer()
            redirect_stdout(buf) do
                show(option)
            end
            @test String(take!(buf)) == "CounterOption(:a)"
        end
    end

    @testset "consume!(); $(v[1])" for v in [
        (["-v"], (1, 1)),
        (["--verbose"], (1, 1)),
        (["-q"], (1, -1)),
        (["--quiet"], (1, -1)),
    ]
        args, expected = v
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        option = CounterOption("-v", "--verbose"; decrementers = ["-q", "--quiet"])
        num_consumed = CliOptions.consume!(d, option, args, ctx)
        @test num_consumed == expected[1]
        @test d["verbose"] == expected[2]
    end

    @testset "consume!(); type, $v" for v in [
        Int, Int8, Int128,
    ]
        T = v
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        option = CounterOption(T, "-v")
        CliOptions.consume!(d, option, ["-v"], ctx)
        @test typeof(d["v"]) == T
    end

    @testset "consume!(); overflow, $(v[1][1]) × $(v[2])" for v in [
        (["-v"], 127, 127),
        (["-v"], 128, CliOptionError),
        (["-q"], 128, -128),
        (["-q"], 129, CliOptionError),
    ]
        args, count, expected = v
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        option = CounterOption(Int8, "-v"; decrementers = "-q")
        for _ in 1:count-1
            CliOptions.consume!(d, option, args, ctx)
        end
        if expected isa Type
            @test_throws expected CliOptions.consume!(d, option, args, ctx)
        else
            CliOptions.consume!(d, option, args, ctx)
            @test d["v"] == expected
        end
    end
end


text:
using SARSOP
using POMDPs
using POMDPLinter
using POMDPModels
using BeliefUpdaters
using POMDPSimulators
using POMDPModelTools
using QuickPOMDPs
using Test

using SARSOP: evaluate # to avoid ambiguity

@testset "POMDPFile" begin
    pomdp_file = POMDPFile(joinpath(@__DIR__, "tiger.pomdpx"))
    @test isfile(pomdp_file.filename)
    tiger_file = POMDPFile(TigerPOMDP())
    @test isfile(tiger_file.filename)
end

@testset "Basic interface" begin
    solver = SARSOPSolver(precision=0.01, timeout=10.0, memory=100.0, fast=true,
                          randomization=true,
                          trial_improvement_factor=1e-3)
    pomdp = TigerPOMDP()
    policy = solve(solver, pomdp)
    mdp = SimpleGridWorld()
    @test_throws ErrorException solve(solver, mdp)
    POMDPLinter.@show_requirements solve(solver, pomdp)
    sim = RolloutSimulator(max_steps=100)
    r = simulate(sim, pomdp, policy, DiscreteUpdater(pomdp))
end

@testset "Simulator" begin
    simulator = SARSOPSimulator(sim_len=5, sim_num=5, fast=true, srand = 1)
    simulate(simulator)
    evaluator = SARSOPEvaluator(sim_len=5, sim_num=5, fast=true, srand = 1, memory=200.0)
    evaluate(evaluator)
end

@testset "Policy Graph" begin
    graph_generator = PolicyGraphGenerator(graph_filename="tiger.dot",
                                           pomdp_filename="model.pomdpx",
                                           fast = true)
    generate_graph(graph_generator)
    @test isfile(graph_generator.graph_filename)
end

@testset "Load policy" begin
    pomdp = TigerPOMDP()
    policy = solve(SARSOPSolver(verbose=false), pomdp)
    policy2 = load_policy(pomdp, "policy.out")
    @test policy.alphas == policy2.alphas
    @test_throws ErrorException load_policy(pomdp, "dum.out")
end

@testset "Issue #39" begin
    solve(SARSOPSolver(), BabyPOMDP())
end

@testset "Issue #40" begin
    include("cancer.jl")
end


text:
import Pkg; Pkg.activate("../../")

include("GPs.jl")

using Distributions
using PyPlot
using BenchmarkTools
import Random

Random.seed!(4)
n = 10
x = randn(n)
f(x) =  sin(3 * x) * sin(x) * (-1)^(x > 0)
sigma = 0
y = f.(x) + randn(n) * sigma;

xgrid = range(minimum(x) - 0.3, maximum(x) + 0.3, length=100)
rhos = rand(Uniform(0.4, 0.6), 50)

p = [let
    gp = GPs.GP(y, x[:,:], GPs.SqExpKernel(1, rho), sigma=0);
    post = GPs.posterior(gp, xgrid[:,:])
    rand(post)
end for rho in rhos]
p = hcat(p...)

plt.scatter(x, y)
p_lower, p_upper = GPs.ci(p)
plt.plot(xgrid, mean(p, dims=2))
plt.fill_between(xgrid, p_lower, p_upper, alpha=0.2, color="blue")

for i in sample(1:size(p, 2), 10)
    plt.plot(xgrid, p[:, i], alpha=0.2, color="black")
end

plt.savefig("out/post.pdf")
plt.close()

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
export ExternalSamplingMCCFRPolicy

using Random
using StatsBase: sample, Weights

"""
    ExternalSamplingMCCFRPolicy

This implementation uses stochasticaly-weighted averaging.

Ref:

- [MONTE CARLO SAMPLING AND REGRET MINIMIZATION FOR EQUILIBRIUM COMPUTATION AND DECISION-MAKING IN LARGE EXTENSIVE FORM GAMES](http://mlanctot.info/files/papers/PhD_Thesis_MarcLanctot.pdf)
- [Monte Carlo Sampling for Regret Minimization in Extensive Games](https://papers.nips.cc/paper/3713-monte-carlo-sampling-for-regret-minimization-in-extensive-games.pdf)
"""
struct ExternalSamplingMCCFRPolicy{S,T,R<:AbstractRNG} <: AbstractPolicy
    nodes::Dict{S,InfoStateNode}
    behavior_policy::QBasedPolicy{TabularLearner{S,T},WeightedExplorer{true,R}}
end

(p::ExternalSamplingMCCFRPolicy)(env::AbstractEnv) = p.behavior_policy(env)

RLBase.get_prob(p::ExternalSamplingMCCFRPolicy, env::AbstractEnv) =
    get_prob(p.behavior_policy, env)

function ExternalSamplingMCCFRPolicy(;
    env::AbstractEnv,
    n_iter::Int,
    rng = Random.GLOBAL_RNG,
)
    @assert NumAgentStyle(env) isa MultiAgent
    @assert DynamicStyle(env) === SEQUENTIAL
    @assert RewardStyle(env) === TERMINAL_REWARD
    @assert ChanceStyle(env) === EXPLICIT_STOCHASTIC
    @assert DefaultStateStyle(env) === Information{String}()

    nodes = init_info_state_nodes(env)

    for i in 1:n_iter
        for p in get_players(env)
            if p != get_chance_player(env)
                external_sampling(copy(env), p, nodes, rng)
            end
        end
    end

    behavior_policy = QBasedPolicy(;
        learner = TabularLearner{String}(),
        explorer = WeightedExplorer(; is_normalized = true, rng = rng),
    )

    for (k, v) in nodes
        s = sum(v.cumulative_strategy)
        if s != 0
            update!(behavior_policy, k => v.cumulative_strategy ./ s)
        end
    end

    ExternalSamplingMCCFRPolicy(nodes, behavior_policy)
end

function external_sampling(env, i, nodes, rng)
    current_player = get_current_player(env)

    if get_terminal(env)
        get_reward(env, i)
    elseif current_player == get_chance_player(env)
        env(rand(rng, get_actions(env)))
        external_sampling(env, i, nodes, rng)
    else
        I = get_state(env)
        node = nodes[I]
        regret_matching!(node)
        σ, rI, sI = node.strategy, node.cumulative_regret, node.cumulative_strategy
        n = length(node.strategy)

        if i == current_player
            u = zeros(n)
            uσ = 0
            for (aᵢ, a) in enumerate(get_legal_actions(env))
                u[aᵢ] = external_sampling(child(env, a), i, nodes, rng)
                uσ += σ[aᵢ] * u[aᵢ]
            end
            rI .+= u .- uσ
            uσ
        else
            a′ = sample(rng, Weights(σ, 1.0))
            env(get_legal_actions(env)[a′])
            u = external_sampling(env, i, nodes, rng)
            sI .+= σ
            u
        end
    end
end


text:
"""
Signature(𝑟, 𝑠) represents a Signature for 𝑧 in which
```𝑠 = (𝑧 + 𝑟𝑒) / 𝑘```
𝑘 being a random integer.
"""
struct Signature
    𝑟::BigInt
    𝑠::BigInt
    Signature(𝑟, 𝑠) = new(𝑟, 𝑠)
end

"Formats Signature as (r, s) in hexadecimal format"
function show(io::IO, z::Signature)
    print(io, "scep256k1 signature(𝑟, 𝑠):\n", string(z.𝑟, base = 16), ",\n", string(z.𝑠, base = 16))
end

==(x::Signature, y::Signature) = x.𝑟 == y.𝑟 && x.𝑠 == y.𝑠

"""
Serialize a Signature to DER format

sig2der(x::Signature) -> Array{UInt8,1}
"""
function sig2der(x::Signature)
    rbin = int2bytes(x.𝑟)
    # if rbin has a high bit, add a 00
    if rbin[1] >= 128
        rbin = pushfirst!(rbin, 0x00)
    end
    result = cat([0x02], int2bytes(length(rbin)), rbin; dims=1)
    sbin = int2bytes(x.𝑠)
    # if sbin has a high bit, add a 00
    if sbin[1] >= 128
        sbin = pushfirst!(sbin, 0x00)
    end
    result = cat(result, [0x02], int2bytes(length(sbin)), sbin; dims=1)
    return cat([0x30], int2bytes(length(result)), result; dims=1)
end

"""
Parse a DER binary to a Signature

der2sig(signature_bin::AbstractArray{UInt8}) -> Signature
"""
function der2sig(signature_bin::AbstractArray{UInt8})
    s = IOBuffer(signature_bin)
    bytes = UInt8[]
    readbytes!(s, bytes, 1)
    if bytes[1] != 0x30
        throw(DomainError("Bad Signature"))
    end
    readbytes!(s, bytes, 1)
    if bytes[1] + 2 != length(signature_bin)
        throw(DomainError("Bad Signature Length"))
    end
    readbytes!(s, bytes, 1)
    if bytes[1] != 0x02
        throw(DomainError("Bad Signature"))
    end
    readbytes!(s, bytes, 1)
    rlength = Int(bytes[1])
    readbytes!(s, bytes, rlength)
    r = bytes2hex(bytes)
    readbytes!(s, bytes, 1)
    if bytes[1] != 0x02
        throw(DomainError("Bad Signature"))
    end
    readbytes!(s, bytes, 1)
    slength = Int(bytes[1])
    bytes = UInt8[]
    readbytes!(s, bytes, slength)
    s = bytes2hex(bytes)
    if length(signature_bin) != 6 + rlength + slength
        throw(DomainError("Signature too long"))
    end
    return Signature(parse(BigInt, r, base=16),
                     parse(BigInt, s, base=16))
end


text:
using Base.Iterators
using Distributions
using LinearAlgebra
using Random
using RCall

# Define the DGP
function simulate_sample(N0, N1)
    # Set the conditional density
    function dy(xc)
        atoms = [(-xc, 0.2), (xc, 0.2)]
        weights = [0.5, 0.5]
        # weights = [exp(-2xc), 1 - exp(-2xc)]
        MixtureModel(Normal, atoms, weights)
    end

    # Simulate a sample
    x0 = rand(N0)
    y0 = [rand(dy(xc)) for xc in x0]
    X0 = rcopy(R"cbind(1, splines::bs($x0, df = 6, Boundary.knots = c(0, 1)))")

    # Generate the grid points
    y_grid = LinRange(-3.0, 3.0, N1)
    x_grid = [0.25, 0.5, 0.75]
    y1, x1 = Iterators.product(y_grid, x_grid) |> x -> zip(x...) .|> collect
    X1 = rcopy(R"cbind(1, predict(splines::bs($x0, df = 6, Boundary.knots = c(0, 1)), $x1))")
    return dy, y0, x0, X0, y1, x1, X1
end

# # Not run
# Random.seed!(1);
# N0, N1 = 500, 50;
# dy, y0, x0, X0, y1, x1, X1 = simulate_sample(N0, N1);


text:
struct AE
	encoder
	decoder
end

Flux.@functor AE

function (ae::AE)(x)
	z = ae.encoder(x)
	x̂ = ae.decoder(z)
	return x̂
end

# ae_loss already in svdd

function conv_ae_constructor(
	;idim=(2,2,1),
	zdim::Int=1,
	activation="relu",
	kernelsizes=(1,1),
	channels=(1,1),
	scalings=(1,1),
	init_seed=nothing,
	batchnorm=false,
	kwargs...
)
	# if seed is given, set it
	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing
	# encoder
	encoder = Chain(
		conv_encoder(idim, zdim, kernelsizes, channels, scalings; activation=activation, batchnorm=batchnorm)...
	)
	# decoder
	decoder = Chain(
		conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels), reverse(scalings),
			activation=activation, vec_output=false, vec_output_dim=nothing, batchnorm=batchnorm)...,
			x->sigmoid.(x) # our dataset images are scaled from 0-1
	)
	# reset seed
	(init_seed !== nothing) ? Random.seed!() : nothing

	model = AE(encoder, decoder)
end

function ae_constructor(
	;idim::Int=1,
	zdim::Int=1,
	hdim::Int=64,
	activation = "relu",
	nlayers::Int=3,
	init_seed=nothing,
	kwargs...
	)
	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing

	encoder = build_mlp(idim, hdim, zdim, nlayers, activation=activation, lastlayer="linear")

	decoder = build_mlp(zdim, hdim, idim, nlayers, activation=activation, lastlayer="linear")

	(init_seed !== nothing) ? Random.seed!() : nothing

	return AE(encoder, decoder)
end

function anomaly_score(ae::AE, real_input; dims=(1,2,3), to_testmode::Bool=true)
	(to_testmode == true) ? Flux.testmode!(ae) : nothing
	score = mse_score(ae(real_input), real_input, dims=dims)
	(to_testmode == true) ? Flux.testmode!(ae, false) : nothing
	return score
end

function anomaly_score_gpu(ae::AE, real_input; batch_size=64, dims=(1,2,3), to_testmode::Bool=true)
	real_input = Flux.Data.DataLoader(real_input, batchsize=batch_size)
	(to_testmode == true) ? Flux.testmode!(ae) : nothing
	ae = ae |> gpu
	scores = Array{Float32}([])
	for x in real_input
		x = x |> gpu
		score = mse_score(ae(x), x, dims=dims) #vec(Flux.sum((ae(x) .- x).^2, dims=dims))
		scores = cat(scores, score |> cpu, dims=1)
	end
	(to_testmode == true) ? Flux.testmode!(ae, false) : nothing
	return scores
end

function StatsBase.fit!(ae::AE, data, params)
	# prepare batches & loaders
	train_loader, val_loader = prepare_dataloaders(data, params)
	history = MVHistory()
	# prepare for early stopping
	best_val_loss = Inf

	patience = params.patience
	val_batches = length(val_loader)

	if ndims(first(train_loader)) == 2
		dims = 1
	elseif ndims(first(train_loader)) == 4
		dims = (1,2,3)
	else
		error("unknown data type (no 2D or 4D tensor)")
	end

	#model to gpu
	ae = ae |> gpu
	# ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) = Optimiser(ADAM(η, β), WeightDecay(decay))
	opt_ae = haskey(params, :decay) ? ADAMW(params.lr, (0.9, 0.999), params.decay) : ADAM(params.lr)

	best_ae = deepcopy(ae)

	# train "AE"
	progress = Progress(params.iters)
	ps_ae = Flux.params(AE)

	for (iter, X) in enumerate(train_loader)
		X = getobs(X)|>gpu
		loss_ae, back = Flux.pullback(ps_ae) do
			ae_loss(ae(X), X, dims=dims)
		end
		grad_ae = back(1f0)
		Flux.Optimise.update!(opt_ae, ps_ae, grad_ae)

		push!(history, :loss_ae, iter, loss_ae)

		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:ae_loss, loss_ae)
			])

		if mod(iter, params.check_every) == 0
			total_val_loss = 0
			Flux.testmode!(ae)
			for X_val in val_loader
				X_val = X_val |> gpu
				total_val_loss += ae_loss(ae(X), X)
			end
			Flux.testmode!(ae, false)
			push!(history, :validation_ae_loss, iter, total_val_loss/val_batches)
			if total_val_loss < best_val_loss
				best_val_loss = total_val_loss
				patience = params.patience
				best_ae = deepcopy(ae)
			else
				patience -= 1
				if patience == 0
					@info "Stopped training of Autoencoder after $(iter) iterations"
					global iters = iter - params.check_every*params.patience
					break
				end
			end
		end
		if iter == params.iters
			global iters = params.iters
		end
	end
	return history, best_ae, sum(map(p->length(p), Flux.params(ae))), iters
end


text:
# regression test for 5442a4e

@testset "if statement" begin
    function bar1()
        if false
            return 1
        else
            return 0
        end
    end
    @test TinyCassette.execute(nothing, bar1) == 0

    @noinline foo() = false
    function bar2()
        if foo()
            return 1
        else
            return 0
        end
    end
    @test TinyCassette.execute(nothing, bar2) == 0
end


text:
@testset "variableoutputcpnn.jl" begin

    @testset "VariableOutputCPNN" begin
        modelNN = SeaPearl.VariableOutputCPNN(
            graphChain = Flux.Chain(
                SeaPearl.GraphConv(3 => 3),
                SeaPearl.GraphConv(3 => 3),
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(3, 3),
                Flux.Dense(3, 3),
            ),
            outputChain = Flux.Chain(
                Flux.Dense(6, 1)
            )
        )

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(3, 10) for i = 1:4]
        featuredGraphs = [SeaPearl.FeaturedGraph(g; nf=nf) for (g, nf) in zip(graphs, nodeFeatures)]

        trajectoryVector = SeaPearl.TsptwTrajectoryState.(featuredGraphs, rand(1:10, 4), [rand(1:10, rand(1:5)) for i = 1:4])
        actionSize = map(ts -> length(ts.possibleValuesIdx), trajectoryVector)
        nnInput = trajectoryVector |> cpu

        output = modelNN.(nnInput)

        @test all(length.(output) .== actionSize)
    end
end