text:
using ReinforcementLearningBase

mutable struct TicTacToe <: AbstractEnv
    board::BitArray{3}
    player::Int
end

function TicTacToe()
    board = BitArray{3}(undef, 3, 3, 3)
    fill!(board, false)
    board[:, :, 3] .= true
    TicTacToe(board, 1)
end

function RLBase.reset!(env::TicTacToe)
    fill!(env.board, false)
    env.board[:, :, 3] .= true
    env.player = 1
    return env.board
end

const TIC_TAC_TOE_STATE_INFO = Dict{
    TicTacToe,
    NamedTuple{
        (:index, :is_terminated, :winner),
        Tuple{Int,Bool,Union{Nothing,Int}},
    },
}()

Base.hash(env::TicTacToe, h::UInt) = hash(env.board, h)
Base.isequal(a::TicTacToe, b::TicTacToe) = isequal(a.board, b.board)

RLBase.action_space(::TicTacToe) = Base.OneTo(9)

RLBase.legal_action_space(env::TicTacToe, p) = findall(legal_action_space_mask(env))

function RLBase.legal_action_space_mask(env::TicTacToe, p)
    if is_win(env, 1) || is_win(env, 2)
        zeros(false, 9)
    else
        vec(view(env.board, :, :, 3))
    end
end

(env::TicTacToe)(action::Int) = env(CartesianIndices((3, 3))[action])

function (env::TicTacToe)(action::CartesianIndex{2})
    env.board[action, 3] = false
    env.board[action, env.player] = true
    env.player = mod1(env.player+1,2)
	return env.board
end

RLBase.current_player(env::TicTacToe) = env.player
RLBase.players(env::TicTacToe) = (1, 2)

RLBase.state(env::TicTacToe, ::Observation{BitArray{3}}, p) = env.board
RLBase.state_space(env::TicTacToe, ::Observation{BitArray{3}}, p) =
    Space(fill(p, 3, 3, 3))
RLBase.state(env::TicTacToe, ::Observation{Int}, p) =
    get_tic_tac_toe_state_info()[env].index
RLBase.state_space(env::TicTacToe, ::Observation{Int}, p) =
    Base.OneTo(length(get_tic_tac_toe_state_info()))

RLBase.state_space(env::TicTacToe, ::Observation{String}, p) = WorldSpace{String}()

function RLBase.state(env::TicTacToe, ::Observation{String}, p)
    buff = IOBuffer()
    for i in 1:3
        for j in 1:3
            if env.board[i, j, 1]
                x = 'o'
            elseif env.board[i, j, 2]
                x = 'x'
            else
                x = '.'
            end
            print(buff, x)
        end
        print(buff, '\n')
    end
    String(take!(buff))
end

RLBase.is_terminated(env::TicTacToe) = get_tic_tac_toe_state_info()[env].is_terminated

function RLBase.reward(env::TicTacToe, player)
    if is_terminated(env)
        winner = get_tic_tac_toe_state_info()[env].winner
        if isnothing(winner)
            0
        elseif winner === player
            1
        else
            -1
        end
    else
        0
    end
end

function is_win(env::TicTacToe, player)
    b = env.board
    p = env.player
    @inbounds begin
        b[1, 1, p] & b[1, 2, p] & b[1, 3, p] ||
            b[2, 1, p] & b[2, 2, p] & b[2, 3, p] ||
            b[3, 1, p] & b[3, 2, p] & b[3, 3, p] ||
            b[1, 1, p] & b[2, 1, p] & b[3, 1, p] ||
            b[1, 2, p] & b[2, 2, p] & b[3, 2, p] ||
            b[1, 3, p] & b[2, 3, p] & b[3, 3, p] ||
            b[1, 1, p] & b[2, 2, p] & b[3, 3, p] ||
            b[1, 3, p] & b[2, 2, p] & b[3, 1, p]
    end
end

function get_tic_tac_toe_state_info()
    if isempty(TIC_TAC_TOE_STATE_INFO)
        # @info "initializing state info..."
        t = @elapsed begin
            n = 1
            root = TicTacToe()
            TIC_TAC_TOE_STATE_INFO[root] =
                (index = n, is_terminated = false, winner = nothing)
            walk(root) do env
                if !haskey(TIC_TAC_TOE_STATE_INFO, env)
                    n += 1
                    has_empty_pos = any(view(env.board, :, :, 3))
                    w = if is_win(env, 1)
                        1
                    elseif is_win(env, 2)
                        2
                    else
                        nothing
                    end
                    TIC_TAC_TOE_STATE_INFO[env] = (
                        index = n,
                        is_terminated = !(has_empty_pos && isnothing(w)),
                        winner = w,
                    )
                end
            end
        end
        # @info "finished initializing state info in $t seconds"
    end
    TIC_TAC_TOE_STATE_INFO
end

RLBase.NumAgentStyle(::TicTacToe) = MultiAgent(2)
RLBase.DynamicStyle(::TicTacToe) = SEQUENTIAL
RLBase.ActionStyle(::TicTacToe) = FULL_ACTION_SET
RLBase.InformationStyle(::TicTacToe) = PERFECT_INFORMATION
RLBase.StateStyle(::TicTacToe) =
    (Observation{String}(), Observation{Int}(), Observation{BitArray{3}}())
RLBase.RewardStyle(::TicTacToe) = TERMINAL_REWARD
RLBase.UtilityStyle(::TicTacToe) = ZERO_SUM
RLBase.ChanceStyle(::TicTacToe) = DETERMINISTIC

using Crayons

player_color(p) = p == 1 ? crayon"light_red" : crayon"light_blue"
player_name(p)  = p == 1 ? "1" : "2"
player_mark(p)  = p == 1 ? "o" : "x"

function render_game(env::AbstractEnv; with_position_names=true, botmargin=true)
	p=current_player(env)
	pname = player_name(p)
	pcol = player_color(p)
	print(pcol, pname, " plays:", crayon"reset", "\n\n")
	board_red=env.board[:,:,1]
	board_blue=env.board[:,:,2]
	indices=LinearIndices(board_blue)
	for x in 1:3
		for y in 1:3
			if board_red[x,y]==0 && board_blue[x,y]==0
				print(".")
			elseif board_red[x,y]==1
				print(player_color(1), player_mark(1), crayon"reset")
			elseif board_blue[x,y]==1
				print(player_color(2), player_mark(2), crayon"reset")
			end

			print(" ")
		end
		if with_position_names
			print(" | ")
			for y in 1:3
				print("$(indices[CartesianIndex(x,y)]) ")
			end
		end
		print("\n")
	end
  	botmargin && print("\n")
end

function human_input()::Int
	print("Please take an action:\n")
	action_taken = parse(Int, readline())
	return action_taken
end


text:
abstract type ClassificationModel{Class} end
labels(model::Type{<:ClassificationModel{Class}}) where {Class} = labels(Class)
labels(model::ClassificationModel{Class}) where {Class} = labels(Class)
Base.broadcastable(m::ClassificationModel) = Ref(m)

ext: jl


text:
module SchwarzChristoffel

using Reexport

include("MapTypes.jl")
@reexport using .MapTypes

include("Polygons.jl")
@reexport using .Polygons

include("Exterior.jl")
@reexport using .Exterior

export Polygons, Exterior

# plotting stuff
include("plot_recipes.jl")

end


text:
# Possible actions:
#   :preprocess
#   :bibentry
#   :textentry
abstract type CitationStyle <: Selectors.AbstractSelector end

abstract type DefaultStyle <: CitationStyle end
Selectors.order(::Type{DefaultStyle}) = 1.0
Selectors.matcher(::Type{DefaultStyle}, fmt, doc, node) = fmt === :default
Selectors.runner(::Type{DefaultStyle}, fmt, doc, node) = format_bibentries(node)

function format_bibentries(node::BibliographyNode)
    unique!(node.keys)

    for key in node.keys
        if !haskey(node.bibdatabase, key)
            @warn "Could not find information for key $key."
            node.database[key] = (key, key)
            continue
        end
        node.database[key] = format_entry(node.bibdatabase[key])
    end
end

format_entry(entry) = (format_bibentry(entry), format_textentry(entry))

function format_bibentry(entry)
    authors = strip.(split(entry["author"], "and"))
    authorstr = ""
    for author in authors
        if authorstr != ""
            authorstr *= ", "
        end

        name = split(author)
        for piece in name[1:end-1]
            authorstr *= string(uppercase(piece[1]), ". ")
        end
        authorstr *= name[end]
    end

    journal = get(entry, "journal", "")

    citation = string(authorstr, ", \"", entry["title"], "\", ",
        journal, " (", entry["year"], ").")

    return citation
end

function format_textentry(entry)
    authors = strip.(split(entry["author"], "and"))
    lastname = split(authors[1])[end]
    year = entry["year"]

    return string("[", lastname, year, "]")
end


text:
using Kwant, PyPlot, PyCall
kwant = Kwant; pyplot = PyPlot

include("graphene.jl")

pot = 0.1
syst, leads = make_system(pot=pot)

# To highlight the two sublattices of graphene, we plot one with
# a filled, and the other one with an open circle:
py"""
    def family_colors(site):
        return 0 if site.family == $(a.o) else 1
    """

# Plot the closed system without leads.
kwant.plot(syst, site_color=py"family_colors", site_lw=0.1, colorbar=false)

# Compute some eigenvalues.
compute_evs(syst.finalized())

    # Attach the leads to the system.
for lead in leads
    syst.attach_lead(lead)
end

    # Then, plot the system with leads.
kwant.plot(syst, site_color=py"family_colors", site_lw=0.1,
           lead_site_lw=0, colorbar=false)

# Finalize the system.
system = syst.finalized()

    # Compute the band structure of lead 0.
momenta = [-pi + 0.02 * pi * i for i in range(0,length=101)]
plot_bandstructure(system.leads[1], momenta)

# Plot conductance.
energies = [-2 * pot + 4. / 50. * pot * i for i in range(0,length=51)]
plot_conductance(system, energies)


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    struct ModifyConstraintNotAllowed{F<:AbstractFunction, S<:AbstractSet,
                                      C<:AbstractFunctionModification} <: NotAllowedError
        constraint_index::ConstraintIndex{F, S}
        change::C
        message::String
    end

An error indicating that the constraint modification `change` cannot be applied
to the constraint of index `ci`.
"""
struct ModifyConstraintNotAllowed{
    F<:AbstractFunction,
    S<:AbstractSet,
    C<:AbstractFunctionModification,
} <: NotAllowedError
    constraint_index::ConstraintIndex{F,S}
    change::C
    message::String
end
function ModifyConstraintNotAllowed(
    ci::ConstraintIndex{F,S},
    change::AbstractFunctionModification,
    message = "",
) where {F<:AbstractFunction,S<:AbstractSet}
    return ModifyConstraintNotAllowed{F,S,typeof(change)}(ci, change, message)
end
function throw_modify_not_allowed(ci::ConstraintIndex, args...)
    return throw(ModifyConstraintNotAllowed(ci, args...))
end

function operation_name(err::ModifyConstraintNotAllowed{F,S}) where {F,S}
    return "Modifying the constraints $(err.constraint_index) with $(err.change)"
end

"""
    struct ModifyObjectiveNotAllowed{C<:AbstractFunctionModification} <: NotAllowedError
        change::C
        message::String
    end

An error indicating that the objective modification `change` cannot be applied
to the objective.
"""
struct ModifyObjectiveNotAllowed{C<:AbstractFunctionModification} <:
       NotAllowedError
    change::C
    message::String
end
function ModifyObjectiveNotAllowed(change::AbstractFunctionModification)
    return ModifyObjectiveNotAllowed(change, "")
end
function throw_modify_not_allowed(::ObjectiveFunction, args...)
    return throw(ModifyObjectiveNotAllowed(args...))
end

function operation_name(err::ModifyObjectiveNotAllowed)
    return "Modifying the objective function with $(err.change)"
end

"""
## Constraint Function

    modify(model::ModelLike, ci::ConstraintIndex, change::AbstractFunctionModification)

Apply the modification specified by `change` to the function of constraint `ci`.

An [`ModifyConstraintNotAllowed`](@ref) error is thrown if modifying
constraints is not supported by the model `model`.

### Examples

```julia
modify(model, ci, ScalarConstantChange(10.0))
```

## Objective Function

    modify(model::ModelLike, ::ObjectiveFunction, change::AbstractFunctionModification)

Apply the modification specified by `change` to the objective function of
`model`. To change the function completely, call `set` instead.

An [`ModifyObjectiveNotAllowed`](@ref) error is thrown if modifying
objectives is not supported by the model `model`.

### Examples

```julia
modify(model, ObjectiveFunction{ScalarAffineFunction{Float64}}(), ScalarConstantChange(10.0))
```

## Multiple modifications in Constraint Functions

    modify(
        model::ModelLike,
        cis::AbstractVector{<:ConstraintIndex},
        changes::AbstractVector{<:AbstractFunctionModification},
    )

Apply multiple modifications specified by `changes` to the functions of constraints `cis`.

A [`ModifyConstraintNotAllowed`](@ref) error is thrown if modifying
constraints is not supported by `model`.

### Examples

```julia
modify(
    model,
    [ci, ci],
    [
        ScalarCoefficientChange{Float64}(VariableIndex(1), 1.0),
        ScalarCoefficientChange{Float64}(VariableIndex(2), 0.5),
    ],
)
```

## Multiple modifications in the Objective Function

    modify(
        model::ModelLike,
        attr::ObjectiveFunction,
        changes::AbstractVector{<:AbstractFunctionModification},
    )

Apply multiple modifications specified by `changes` to the functions of constraints `cis`.

A [`ModifyObjectiveNotAllowed`](@ref) error is thrown if modifying
objective coefficients is not supported by `model`.

### Examples

```julia
modify(
    model,
    ObjectiveFunction{ScalarAffineFunction{Float64}}(),
    [
        ScalarCoefficientChange{Float64}(VariableIndex(1), 1.0),
        ScalarCoefficientChange{Float64}(VariableIndex(2), 0.5),
    ],
)
```
"""
function modify end

function modify(
    model::ModelLike,
    ci::ConstraintIndex,
    change::AbstractFunctionModification,
)
    return throw_modify_not_allowed(ci, change)
end

function modify(
    model::ModelLike,
    cis::AbstractVector{<:ConstraintIndex},
    changes::AbstractVector{<:AbstractFunctionModification},
)
    @assert length(cis) == length(changes)
    for (ci, change) in zip(cis, changes)
        modify(model, ci, change)
    end
    return
end

function modify(
    model::ModelLike,
    attr::ObjectiveFunction,
    change::AbstractFunctionModification,
)
    return throw_modify_not_allowed(attr, change)
end

function modify(
    model::ModelLike,
    attr::ObjectiveFunction,
    changes::AbstractVector{<:AbstractFunctionModification},
)
    for change in changes
        modify(model, attr, change)
    end
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
#Hermite-Gaussian Modes
HG(n,m,x,y) = hermiteh(n,√2*x)*hermiteh(m,√2*y)*exp(-x^2-y^2);

function θ(x,y)
	#θ ∈ [-π,π]
	if x ≥ 0
		atan(y/x)
	else
		if y≥ 0
			atan(y/x) + π
		else
			atan(y/x) - π
		end
	end
end;

function LG(p,l,x,y)
	#Laguerre-Gaussian Modes
    if x==0 && y==0
        0.0
    else
        (2*(x^2+y^2))^(0.5*abs(l))*laguerrel(p,abs(l),2*(x^2+y^2))*exp(-(x^2+y^2))*exp(im*l*θ(x,y))
    end
end

function airy_beam(x,α)
	#Airy beam (in a single direction)
	airyai(x)*exp(α*x)
end

#Series of obstacles to simulate difraction
function vertical_obstacle(x,y,l)
	if abs(x)<l
		0.0
	else
		1.0
	end
end

function horizontal_obstacle(x,y,l)
	if abs(y)<l
		0.0
	else
		1.0
	end
end

function circular_obstacle(x,y,r)
	if x^2+y^2<r
		0.0
	else
		1.0
	end
end

function rectangular_obstacle(x,y,a,b,c,d)
	if (a<x<b) & (c<y<d)
		0.0
	else
		1.0
	end
end


text:
using CSV, DataFrames, Statistics

const minimum_similarity = 0.8
const movies = CSV.read("top_10_movies_user_rankings.tsv", delim = '\t')

function user_similarity(target_user)
  similarity = Dict{Symbol,Float64}()

  for user in names(movies[:, 2:end])
    user == target_user && continue

    ratings = movies[:, [user, target_user]]
    common_movies = ratings[(ratings[user] .> 0) .& (ratings[target_user] .> 0), :]
    # common_movies = ratings[(ratings[user] .> 7) .& (ratings[target_user] .> 0), :]

    correlation = try
      cor(common_movies[user], common_movies[target_user])
    catch
      0.0
    end

    similarity[user] = correlation
  end

  similarity
end

function recommendations(target_user)
  recommended = Dict{String,Float64}()

  for (user,similarity) in user_similarity(target_user)
    similarity > minimum_similarity || continue

    ratings = movies[:, [Symbol("Movie title"), user, target_user]]
    recommended_movies = ratings[(ratings[user] .>= 7) .& (ratings[target_user] .== 0), :]

    for movie in eachrow(recommended_movies)
      recommended[movie[Symbol("Movie title")]] = movie[user] * similarity
    end
  end

  recommended
end

for user in names(movies)[2:end]
  println("Recommendations for $user: $(recommendations(user))")
end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
"""
Vertex (`Vector{Int64}`) -> `DeterministicStrategy`

    convert(
        ::Type{DeterministicStrategy},
        vertex  :: Vector{Int64},
        scenario :: Scenario;
        rep = "normalized" :: String
    )
"""
function convert(::Type{DeterministicStrategy}, vertex::Vector{Int64}, scenario::Union{LocalSignaling,BlackBox}; rep="normalized" :: String)
    if !(rep in ("normalized", "generalized"))
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'"))
    end

    s_dims = strategy_dims(scenario)

    strategy_matrix = (rep == "normalized") ? cat(
        reshape(vertex, (s_dims[1]-1, s_dims[2])),
        zeros(Int64, (1,s_dims[2])),
        dims = 1
    ) : reshape(vertex, s_dims)

    # for vertices in the normalized representation, if a column does not contain
    # a one, add it to the last row.
    if rep == "normalized"
        for col_id in 1:s_dims[2]
            if sum(strategy_matrix[:,col_id]) == 0
                strategy_matrix[s_dims[1], col_id] = 1
            end
        end
    end

    DeterministicStrategy(strategy_matrix, scenario)
end

"""
`DeterministicStrategy` -> Vertex (`Vector{Int64}`)

    convert(
        ::Type{Vector{Int64}},
        strategy :: DeterministicStrategy;
        rep = "normalized" :: String
    )
"""
function convert(::Type{Vector{Int64}}, strategy::DeterministicStrategy; rep = "normalized"::String)
    if !(rep in ("normalized", "generalized"))
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'"))
    end

    max_row = (rep == "normalized") ? size(strategy,1) - 1 : size(strategy,1)

    strategy[1:max_row,:][:]
end

"""
    convert(
        S :: Type{<:AbstractStrategy}, vertex::Vector{<:Real}, scenario::BipartiteNonSignaling;
        rep="non-signaling" :: String
    )

Transforms a behavior vector or vertex in to either a [`Strategy`](@ref) or [`DeterministicStrategy`](@ref).
If converting into a `DeterministicStrategy`, the vertex must contain `Int64` values.
Valid representations are "non-signaling", "normalized", and "generalized".
"""
function convert(S::Type{<:AbstractStrategy}, vertex::Vector{<:Real}, scenario::BipartiteNonSignaling;
    rep="non-signaling" :: String
)
    if !(rep in ["non-signaling","normalized","generalized"])
        throw(DomainError(rep, "input `rep` must be in [\"non-signaling\",\"normalized\",\"generalized\"]"))
    end

    num_type = (S == Strategy) ? Float64 : Int64

    if S == DeterministicStrategy && eltype(vertex) !== num_type
        throw(DomainError((num_type), "for type $S, vertex must be `Vector{$num_type}`."))
    end

    gen_strat_dims = strategy_dims(scenario)
    gen_strat = (rep == "generalized") ? reshape(vertex, gen_strat_dims) : zeros(num_type, gen_strat_dims)

    if rep == "non-signaling"
        # code to get non-signaling to normalized
        α_dim = (scenario.A-1)*scenario.X
        β_dim = (scenario.B-1)*scenario.Y

        α_strat = reshape(vertex[1:α_dim], (scenario.A-1, scenario.X))
        β_strat = reshape(vertex[α_dim+1:α_dim+β_dim], (scenario.B-1, scenario.Y))
        αβ_strat = reshape(vertex[α_dim+β_dim+1:end], ((scenario.A-1)*(scenario.B-1), scenario.X*scenario.Y))

        # reversing non-signaling constraints for Alice
        for i in 0:scenario.A-2
            gen_strat[i*scenario.B+1:i*scenario.B + scenario.B-1,:] = αβ_strat[i*(scenario.B-1)+1:(i+1)*(scenario.B-1),:]

            p_ax_row = kron(α_strat[i+1,:],ones(Int64, scenario.Y))'
            gen_strat[i*scenario.B+scenario.B,:] = p_ax_row - ones(Int64, (1,scenario.B-1)) * αβ_strat[i*(scenario.B-1)+1:(i+1)*(scenario.B-1),:]
        end

        for i in 1:scenario.B-1
            αβ_row_ids = i:scenario.B-1:(scenario.A-1)*(scenario.B-1)

            p_by_row = kron(ones(Int64,scenario.X),β_strat[i,:])'
            gen_strat[(scenario.A-1)*scenario.B + i,:] = p_by_row - ones(Int64, (1,scenario.A-1)) * αβ_strat[αβ_row_ids,:]
        end

    elseif rep == "normalized"
        gen_strat[1:gen_strat_dims[1]-1,:] = reshape(vertex, (gen_strat_dims[1]-1, gen_strat_dims[2]))
    end

    # strategy is in "normalized" rep transform it to "generalized" rep
    if rep in ("normalized", "non-signaling")
        gen_strat[end,:] = 1 .- ones(Int64, (1,gen_strat_dims[1]-1)) * gen_strat[1:end-1,:  ]
    end

    (S == Strategy) ? Strategy(gen_strat) : DeterministicStrategy(gen_strat)
end


text:
using OrdinaryDiffEq, ParameterizedFunctions, ODE, ODEInterfaceDiffEq,
      LSODA, Sundials, DiffEqDevTools

f = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

p = [0.7,0.8,1/12.5,0.5]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);

sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()

plot(sol)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) #fails
          Dict(:alg=>dopri5())
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) # fails
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=10000,error_estimate=:L2,dense_errors=true)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>CVODE_Adams())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
          #Dict(:alg=>lsoda())
          Dict(:alg=>odex())
          Dict(:alg=>ddeabm())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100,maxiters=1000,error_estimate=:L2,dense_errors=true)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using ForwardDiff

println("Testing Partials...")
t = @elapsed include("PartialsTest.jl")
println("done (took $t seconds).")

println("Testing Dual...")
t = @elapsed include("DualTest.jl")
println("done (took $t seconds).")

println("Testing derivative functionality...")
t = @elapsed include("DerivativeTest.jl")
println("done (took $t seconds).")

println("Testing gradient functionality...")
t = @elapsed include("GradientTest.jl")
println("done (took $t seconds).")

println("Testing jacobian functionality...")
t = @elapsed include("JacobianTest.jl")
println("done (took $t seconds).")

println("Testing hessian functionality...")
t = @elapsed include("HessianTest.jl")
println("done (took $t seconds).")

println("Testing perturbation confusion functionality...")
t = @elapsed include("ConfusionTest.jl")
println("done (took $t seconds).")

println("Testing miscellaneous functionality...")
t = @elapsed include("MiscTest.jl")
println("done (took $t seconds).")

# These tests need to be run in a process where bounds checking is not explicitly enabled
# (like they are with Pkg.test)
println("Testing SIMD vectorization...")
project = Base.active_project()
simdfile = joinpath(@__DIR__, "SIMDTest.jl")
t = @elapsed run(`$(Base.julia_cmd()) --check-bounds=no --code-coverage=none -O2 --project=$project $simdfile`)
println("done (took $t seconds).")


text:
##### Beginning of file

module Types # Begin submodule MirrorUpdater.Types

__precompile__(true)

abstract type AbstractInterval end

function _name_with_jl(x::AbstractString)::String
    name_without_jl::String = _name_without_jl(x)
    name_with_jl::String = string(name_without_jl, ".jl")
    return name_with_jl
end

function _name_without_jl(x::AbstractString)::String
    temp::String = strip(convert(String, x))
    if endswith(lowercase(temp), ".jl")
        result = strip(temp[1:end-3])
    else
        result = temp
    end
    return result
end

function _name_with_git(x::AbstractString)::String
    name_without_git::String = _name_without_git(x)
    name_with_git::String = string(name_without_git, ".git")
    return name_with_git
end

function _name_without_git(x::AbstractString)::String
    temp::String = strip(convert(String, x))
    if endswith(lowercase(temp), ".git")
        result = strip(temp[1:end-4])
    else
        result = temp
    end
    return result
end

function _is_interval(x::String)::Bool
    if _is_no_bounds_interval(x)
        return true
    elseif _is_lower_bound_only_interval(x)
        return true
    elseif _is_upper_bound_only_interval(x)
        return true
    elseif _is_lower_and_upper_bound_interval(x)
        return true
    else
        return false
    end
end

function _get_lower_and_upper_bound_interval_regex()::Regex
    lower_and_upper_bound_interval_regex::Regex =
        r"\[(\w\w*?)\,(\w\w*?)\)"
    return lower_and_upper_bound_interval_regex
end

function _get_lower_bound_only_interval_regex()::Regex
    lower_bound_only_interval_regex::Regex =
        r"\[(\w\w*?)\,\)"
    return lower_bound_only_interval_regex
end

function _get_upper_bound_only_interval_regex()::Regex
    upper_bound_only_interval_regex::Regex =
        r"\[\,(\w\w*?)\)"
    return upper_bound_only_interval_regex
end

function _get_no_bounds_interval_regex()::Regex
    no_bounds_interval_regex::Regex =
        r"\[\,\)"
    return no_bounds_interval_regex
end

function _is_no_bounds_interval(x::String)::Bool
    result::Bool = occursin(
        _get_no_bounds_interval_regex(),
        x,
        )
    return result
end

function _is_lower_and_upper_bound_interval(x::String)::Bool
    result::Bool = occursin(
        _get_lower_and_upper_bound_interval_regex(),
        x,
        )
    return result
end

function _is_lower_bound_only_interval(x::String)::Bool
    result::Bool = occursin(
        _get_lower_bound_only_interval_regex(),
        x,
        )
    return result
end

function _is_upper_bound_only_interval(x::String)::Bool
    result::Bool = occursin(
        _get_upper_bound_only_interval_regex(),
        x,
        )
    return result
end

struct Package
    name::String
    uuid::String
    source_url::String
    function Package(
            name::String,
            uuid::String,
            source_url::String,
            )::Package
        correct_name = _name_with_jl(name)
        correct_uuid = strip(uuid)
        correct_source_url = strip(source_url)
        result::Package = new(
            correct_name,
            correct_uuid,
            correct_source_url,
            )
        return result
    end
end

function Package(
        ;
        name::AbstractString,
        uuid::AbstractString,
        source_url::AbstractString,
        )::Package
    result::Package = Package(
        convert(String,name),
        convert(String,uuid),
        convert(String,source_url),
        )
    return result
end

struct Registry
    owner::String
    name::String
    uuid::String
    url::String
    function Registry(
            owner::String,
            name::String,
            uuid::String,
            url::String,
            )::Registry
        correct_owner = strip(owner)
        correct_name = _name_without_jl(name)
        correct_uuid = strip(uuid)
        correct_url = strip(url)
        result::Registry = new(
            correct_owner,
            correct_name,
            correct_uuid,
            correct_url,
            )
        return result
    end
end

function Registry(
        ;
        owner::AbstractString,
        name::AbstractString,
        uuid::AbstractString,
        url::AbstractString,
        )::Registry
    result::Registry = Registry(
        convert(String,owner),
        convert(String,name),
        convert(String,uuid),
        convert(String,url),
        )
    return result
end

struct SrcDestPair
    source_url::String
    destination_repo_name::String
    function SrcDestPair(
            source_url::String,
            destination_repo_name::String,
            )::SrcDestPair
        correct_source_url = strip(source_url)
        correct_destination_repo_name = strip(destination_repo_name)
        result::SrcDestPair = new(
            correct_source_url,
            correct_destination_repo_name,
            )
        return result
    end
end

function SrcDestPair(
    ;
    source_url::AbstractString,
    destination_repo_name::String,
    )::SrcDestPair
    result::SrcDestPair = SrcDestPair(
        convert(String, source_url),
        convert(String, destination_repo_name),
        )
    return result
end

struct NoBoundsInterval <: AbstractInterval
end

struct LowerAndUpperBoundInterval <: AbstractInterval
    left::String
    right::String
    function LowerAndUpperBoundInterval(
            left::String,
            right::String,
            )::LowerAndUpperBoundInterval
        correct_left = strip(left)
        correct_right = strip(right)
        result::LowerAndUpperBoundInterval = new(
            correct_left,
            correct_right,
            )
        return result
    end
end

struct LowerBoundOnlyInterval <: AbstractInterval
    left::String
    function LowerBoundOnlyInterval(
            left::String,
            )::LowerBoundOnlyInterval
        correct_left = strip(left)
        result::LowerBoundOnlyInterval = new(
            correct_left,
            )
        return result
    end
end

struct UpperBoundOnlyInterval <: AbstractInterval
    right::String
    function UpperBoundOnlyInterval(
            right::String,
            )::UpperBoundOnlyInterval
        correct_right = strip(right)
        result::UpperBoundOnlyInterval = new(
            correct_right,
            )
        return result
    end
end

function _construct_interval(x::String)::AbstractInterval
    if _is_no_bounds_interval(x)
        result = NoBoundsInterval()
    elseif _is_lower_bound_only_interval(x)
        loweronly_regexmatch::RegexMatch = match(
            _get_lower_bound_only_interval_regex(),
            x,
            )
        loweronly_left::String = strip(
            convert(String, loweronly_regexmatch[1])
            )
        result = LowerBoundOnlyInterval(loweronly_left)
    elseif _is_upper_bound_only_interval(x)
        upperonly_regexmatch::RegexMatch = match(
            _get_upper_bound_only_interval_regex(),
            x,
            )
        upperonly_right::String = strip(
            convert(String, upperonly_regexmatch[1])
            )
        result = UpperBoundOnlyInterval(upperonly_right)
    elseif _is_lower_and_upper_bound_interval(x)
        lowerandupper_regexmatch::RegexMatch = match(
            _get_lower_and_upper_bound_interval_regex(),
            x,
            )
        lowerandupper_left::String = strip(
            convert(String, lowerandupper_regexmatch[1])
            )
        lowerandupper_right::String = strip(
            convert(String, lowerandupper_regexmatch[2])
            )
        result = LowerAndUpperBoundInterval(
            lowerandupper_left,
            lowerandupper_right,
            )
    else
        error("argument is not a valid interval")
    end
    return result
end

function Base.isless(
        x::SrcDestPair,
        y::SrcDestPair,
        )::Bool
    x_destination_repo_name::String = strip(x.destination_repo_name)
    y_destination_repo_name::String = strip(y.destination_repo_name)
    result::Bool = Base.isless(
        x_destination_repo_name,
        y_destination_repo_name,
        )
    return result
end

end # End submodule MirrorUpdater.Types

##### End of file


text:
using .Meta

export Func, NullFunc, nullfunc, @func

abstract type Func end

struct NullFunc <: Func end
const nullfunc = NullFunc()

struct ParamCreator <: Func end
const paramcreator = ParamCreator()

const SymbolFunc = Union{NullFunc, ParamCreator}

"""
    (f::Func)(args...)
"""
function (f::Func)(args...)
    args = Tuple(asvar(x) for x in args)
    xs = [x.data for x in args]
    ys = forward(f, xs...)
    ys isa Tuple || (ys = (ys,))
    outputs = enable_backprop[] ? [asvar(y, f) for y in ys] : [asvar(y) for y in ys]
    f.args = args
    f.outputs = Tuple(WeakRef(output) for output in outputs)
    return length(outputs) > 1 ? outputs : outputs[1]
end

"""
    @func
"""
macro func(obj)
    @assert obj.head == :struct "@func is for struct type."
    if obj.args[2] isa Symbol
        name = obj.args[2]
        obj.args[2] = :($name <: Func)
    elseif isexpr(obj.args[2], :curly)
        name = obj.args[2].args[1]
        header = obj.args[2]
        if obj.args[2].args[end] != :Func obj.args[2] = :($header <: Func) end
    elseif isexpr(obj.args[2], :<:)
        name = obj.args[2].args[1]
    else
        @error "unexpected args."
    end
    definedfield = obj.args[3]
    obj.args[3] = quote
        $definedfield
        args::Tuple
        outputs::Tuple
    end
    if length(definedfield.args) == 1
        push!(obj.args[3].args, :($name() = new()))
    end
    return esc(obj)
end


text:
function readcondnuc!(read::AbstractVector{Float64},nuc::Int,q::UInt8)
	p = phred2prob(q)
	read[:] = p/3*ones(4)
	read[nuc] = 1-p
	read
end

# Using the implementations in codonfrequency.jl. They work for nucleotides too.
# TODO: Refactor.

function setupproblem(observations::AbstractVector{NucQualityDict})
	@assert size(observations)==(4,) # allow trailing ones?

	# for each of the 4 nucleotides and each unique quality value in them
	# create one 4-vector with p(read|nucleotide)
	n = sum(length,observations)
	rcn = zeros(4,n)
	counts = zeros(1,n)

	nucCounts = zeros(4) # number of observations supporting each nucleotide, used for initial guess only

	i = 1
	for nuc=1:4
		d = observations[nuc]
		s = 0
		for (k,v) in d
			readcondnuc!(view(rcn,:,i),nuc,k)
			counts[1,i] = v
			nucCounts[nuc] += v
			s += v
			i += 1
		end
	end

	s = sum(nucCounts)
	theta0 = nucCounts/(s>0 ? s : 1)

	theta0,rcn,counts
end

# TODO: add logging!

function mlnucfreqs(observations::AbstractVector{NucQualityDict};
					method=:Newton, newtonRegularization=1e-6)
	theta0,rcn,counts = setupproblem(observations)
	coverage = sum(counts)

	all(theta0.==0) && return zeros(4), true, coverage # no (useful) reads - set all freqs to zero

	# early out if we can detect that there is support for a single nucleotide only?

	if method==:Newton
		theta,converged,nbrIter = newtonsolve(theta0,rcn,counts,regularization=newtonRegularization)
	elseif method==:GradientDescent
		theta,converged,nbrIter = gradientdescentsolve(theta0,rcn,counts)
	end

	#theta,obj,nbrIter
	theta, converged, coverage
end

# for a single sequence
function mlnucfreqs(m::NucQualityMap; log=stdout, method=:Newton, newtonRegularization=1e-6)
	len = size(m,2)
	freqs = zeros(4,len)
	coverage = zeros(Int,len)
	for i=1:len
		freqs[:,i], converged, coverage[i] = mlnucfreqs(m[:,i],method=method,newtonRegularization=newtonRegularization)
		converged || println(log, "WARNING: Solver did not converge (position=", i, ')')
	end
	freqs, 1:len, coverage
end

# solve for the whole genome
function mlnucfreqs(nqa::NucQualityAccumulator; log=stdout, method=:Newton, newtonRegularization=1e-6)
	N = length(nqa)
	freqs = Vector{Matrix{Float64}}(undef,N)
	positions = Vector{Vector{Int}}(undef,N)
	coverage = Vector{Vector{Int}}(undef,N)
	for (i,m) in enumerate(nqa)
		freqs[i],positions[i],coverage[i] = mlnucfreqs(m,log=log,method=method,
                                                       newtonRegularization=newtonRegularization)
	end

	freqs,positions,coverage
end

# solve for bam file
function mlnucfreqs(bamFile::BamFile; log=stdout, strands=:both,
                    mappingQualityThreshold=30, baseQualityThreshold=30,
                    ignoreChimeric::Bool=true,
                    method=:Newton, newtonRegularization=1e-6)
	nqa = nucqualitycount(bamFile, strands=strands, mappingQualityThreshold=mappingQualityThreshold,
	                      ignoreChimeric=ignoreChimeric)
	qualityfilter!(nqa, baseQualityThreshold)

	mlnucfreqs(nqa,log=log,method=method,newtonRegularization=newtonRegularization)
end


text:
# ---
# title: 1595. Minimum Cost to Connect Two Groups of Points
# id: problem1595
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Graph
# link: <https://leetcode.com/problems/minimum-cost-to-connect-two-groups-of-points/description/>
# hidden: true
# ---
#
# You are given two groups of points where the first group has `size1` points,
# the second group has `size2` points, and `size1 >= size2`.
#
# The `cost` of the connection between any two points are given in an `size1 x
# size2` matrix where `cost[i][j]` is the cost of connecting point `i` of the
# first group and point `j` of the second group. The groups are connected if
# **each point in both  groups is connected to one or more points in the
# opposite group**. In other words, each point in the first group must be
# connected to at least one point in the second group, and each point in the
# second group must be connected to at least one point in the first group.
#
# Return  _the minimum cost it takes to connect the two groups_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/03/ex1.jpg)
#
#
#
#     Input: cost = [[15, 96], [36, 2]]
#     Output: 17
#     **Explanation** : The optimal way of connecting the groups is:
#     1--A
#     2--B
#     This results in a total cost of 17.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/03/ex2.jpg)
#
#
#
#     Input: cost = [[1, 3, 5], [4, 1, 1], [1, 5, 3]]
#     Output: 4
#     **Explanation** : The optimal way of connecting the groups is:
#     1--A
#     2--B
#     2--C
#     3--A
#     This results in a total cost of 4.
#     Note that there are multiple points connected to point 2 in the first group and point A in the second group. This does not matter as there is no limit to the number of points that can be connected. We only care about the minimum total cost.
#
#
# **Example 3:**
#
#
#
#     Input: cost = [[2, 5, 1], [3, 4, 7], [8, 1, 2], [6, 2, 4], [3, 8, 8]]
#     Output: 10
#
#
#
#
# **Constraints:**
#
#   * `size1 == cost.length`
#   * `size2 == cost[i].length`
#   * `1 <= size1, size2 <= 12`
#   * `size1 >= size2`
#   * `0 <= cost[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Base.Test
using QuantumOptics

ωc = 1.2
ωa = 0.9
g = 1.0
Ncutoff = 10

T = Float64[0.,1.]

fockbasis = FockBasis(Ncutoff)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = g/2*(sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis))

H = Ha + Hc + Hint

Omega_n(n::Int) = sqrt((ωa-ωc)^2 + g^2*(n+1))
Alpha_n(n::Int) = atan2(g*sqrt(n+1),(ωa-ωc))

Em(n::Int) = ωc*(n+0.5) - 0.5*Omega_n(n)
Ep(n::Int) = ωc*(n+0.5) + 0.5*Omega_n(n)

psi_e(n::Int) = spinup(spinbasis) ⊗ fockstate(fockbasis, n)
psi_g(n::Int) = spindown(spinbasis) ⊗ fockstate(fockbasis, n+1)

psi_p(n::Int) = cos(Alpha_n(n)/2)*psi_e(n) + sin(Alpha_n(n)/2)*psi_g(n)
psi_m(n::Int) = -sin(Alpha_n(n)/2)*psi_e(n) + cos(Alpha_n(n)/2)*psi_g(n)
basisstates_p = Ket[psi_p(n) for n=0:Ncutoff-1]
basisstates_m = Ket[psi_m(n) for n=0:Ncutoff-1]

Emvec = Float64[Em(n) for n=0:Ncutoff-1]
Epvec = Float64[Ep(n) for n=-1:Ncutoff-1]

P = sortperm([Emvec; Epvec])
basisstates = Ket[basisstates_m; psi_g(-1); basisstates_p][P][1:(Ncutoff+1)]
E = [Emvec; Epvec][P][1:(Ncutoff+1)]

@test norm(operatorspectrum(full(H))[1:(Ncutoff+1)] - E) < 1e-12
@test norm(operatorspectrum_hermitian(full(H))[1:(Ncutoff+1)] - E) < 1e-12
@test norm(operatorspectrum(H)[1:(Ncutoff+1)] - E) < 1e-12
@test norm(operatorspectrum_hermitian(H)[1:(Ncutoff+1)] - E) < 1e-12

b = eigenstates(full(H))[1:(Ncutoff+1)]
for i=1:length(b)
    @test 1-abs(dagger(b[i])*basisstates[i])<1e-12
end

b = eigenstates_hermitian(full(H))[1:(Ncutoff+1)]
for i=1:length(b)
    @test 1-abs(dagger(b[i])*basisstates[i])<1e-12
end

b = eigenstates(H)[1:(Ncutoff+1)]
for i=1:length(b)
    @test 1-abs(dagger(b[i])*basisstates[i])<1e-12
end

b = eigenstates_hermitian(H)[1:(Ncutoff+1)]
for i=1:length(b)
    @test 1-abs(dagger(b[i])*basisstates[i])<1e-12
end


text:
#=
# Hyperdimensional computing interface

date: 21 October 2021

author: Michiel Stock

In this document, I present a Julia interface for computing with hyperdimensional
vectors (HDV). These are vectors of large dimensionality (typically 10,000) that
allow for distributed and holographic representing, storing and computing with patterns.

## Hypervectors

We first define a root abstract `AbstractHDV` type as a subtype of `AbstractVector`,
meaning we inherit all behaviour expected from a vector.
=#

using LinearAlgebra

abstract type AbstractHDV{T} <: AbstractVector{T} end

Base.getindex(hdv::AbstractHDV, i) = hdv.v[(i+hdv.offset)%length(hdv)+1]
Base.size(hdv::AbstractHDV) = size(hdv.v)
Base.setindex!(hdv::AbstractHDV, val, i) = (hdv.v[(i+hdv.offset)%length(hdv)+1] = val)
normalize!(::AbstractHDV) = nothing  ## vectors have no normalization by default

# The first concrete type is `BipolarHDV`, which stores patterns with values $\{-1, 0, 1\}$.
# All concrete type store a vector with the values and an offset, allowing for inplace
# permutations.

mutable struct BipolarHDV <: AbstractHDV{Int}
    v::Vector{Int}
    offset::Int
    BipolarHDV(v::Vector, offset=0) = new(v, offset)
end

# We always provide a constructor with optinal dimensionality (n=10,000 by default) and
# a method `similar`.

BipolarHDV(n::Int=10_000) = BipolarHDV(rand((-1, 1), n))
Base.similar(hdv::BipolarHDV) = BipolarHDV(similar(hdv.v))

normalize!(hdv::BipolarHDV) = (hdv.v .= sign.(hdv.v))

# `BinaryHDV` contain binary vectors.

mutable struct BinaryHDV <: AbstractHDV{Bool}
    v::Vector{Bool}
    offset::Int
    BinaryHDV(v::AbstractVector, offset=0) = new(v, offset)
end

BinaryHDV(n::Int=10_000) = BinaryHDV(rand(Bool, n))
Base.similar(hdv::BinaryHDV) = BinaryHDV(similar(hdv.v))

# GradedBipolarHDV are vectors in $[-1, 1]^n$, allowing for graded relations.

mutable struct GradedBipolarHDV{T<:Real} <: AbstractHDV{T}
    v::Vector{T}
    offset::Int
    GradedBipolarHDV(v::AbstractVector, offset=0) = new{eltype(v)}(v, offset)
end

GradedBipolarHDV(T::Type, n::Int=10_000) = GradedBipolarHDV(2rand(T, n).-1)
GradedBipolarHDV(n::Int=10_000) = GradedBipolarHDV(Float32, n)

Base.similar(hdv::GradedBipolarHDV) = GradedBipolarHDV(similar(hdv.v))

normalize!(hdv::GradedBipolarHDV)= (hdv.v .= camp.(hdv.v, -1, 1))

# Finally, `RealHDV` contain real values, drawn from a standard normal distribution
# by default.

mutable struct RealHDV{T<:Real} <: AbstractHDV{T}
    v::Vector{T}
    offset::Int
    RealHDV(v::Vector{T}, offset=0) where {T} = new{T}(v, offset)
end

RealHDV(n::Int=10_000) = RealHDV(randn(n), 1.0)

normalize!(hdv::RealHDV) = (hdv.v .*=  √(length(hdv) / sum(abs2, hdv.v)))

Base.similar(hdv::RealHDV) = RealHDV(similar(hdv.v))

#=
## Computing with HDV

The following table contains the basic operations (two binary and one unary) that
are used to calculate with HDV. Important is that all these representations need to
be *reversible*.

| Operation            | symbol | remark                                                                                |
| -------------------- | ------ | --------------------------------------------------------------------------------------------------------------- |
| Bundling/aggregating | `+`    | combines the information of two vectors into a new vector similar to both                                       |
| Binding              | `*`    | mapping, combines the two vectors in something different from both, preserves distance, distributes of bundling |
| Shifting             | `Π`    | Permutation (in practice cyclic shifting), distributes over addition, conserves distance                        |

There are many options possible, so we use type-based dispatch to give the user control.

### Aggregation

Aggregation is at the highest level represented by `+`. (Question: should we introduce a specific symbol?)
=#

aggr(hdv1::AbstractHDV, hdv2::AbstractHDV) = aggr!(similar(hdv1), hdv1, hdv2)

Base.:+(hdv1::AbstractHDV, hdv2::AbstractHDV) = aggr(hdv1, hdv2)

aggr!(r::BipolarHDV, hdv1::BipolarHDV, hdv2::BipolarHDV) = (r .= hdv1 .+ hdv2)
aggr!(r::BinaryHDV, hdv1::BinaryHDV, hdv2::BinaryHDV) = (r .= hdv1 .& hdv2)  ## Conservative

function aggr!(r::RealHDV, hdv1::RealHDV, hdv2::RealHDV)
    r .=  hdv1 .+ hdv2
    r.norm = hdv1.norm + hdv2.norm
    return r
end

aggr!(r::GradedBipolarHDV, x::GradedBipolarHDV, y::GradedBipolarHDV) = (@. r = x * y / (x * y) + (1-x) * (1-y))

#=

### Binding

=#

bind(hdv1::AbstractHDV, hdv2::AbstractHDV) = bind!(similar(hdv1), hdv1, hdv2)
Base.:*(hdv1::AbstractHDV, hdv2::AbstractHDV) = bind(hdv1, hdv2)

bind!(r::BipolarHDV, hdv1::BipolarHDV, hdv2::BipolarHDV) = (r .= hdv1 .* hdv2)
bind!(r::RealHDV, hdv1::RealHDV, hdv2::RealHDV) = (r .= hdv1 .* hdv2)
bind!(r::BinaryHDV, hdv1::BinaryHDV, hdv2::BinaryHDV) = (r .= hdv1 .⊻ hdv2)
bind!(r::GradedBipolarHDV, x::GradedBipolarHDV, y::GradedBipolarHDV) = (@. r = x + y - x * y)

#=
### Permutation

Vectors can be permuted by performing a cyclic shift.
=#

Base.circshift!(hdv::AbstractHDV, k) = (hdv.offset += k)

function Base.circshift(hdv::AbstractHDV, k)
    hdv = copy(hdv)
    hdv.offset += k
    return hdv
end

Π(hdv::AbstractHDV, k) = circshift(hdv, k)

#=

## Inference and comparision

To perform inference, one has to compute (dis)similarity between vectors.
=#

cos_sim(x::HDV, y::HDV) where {HDV<:AbstractHDV} = (x ⋅ y) / (norm(x.v) * norm(x.v))

jaccard(x::HDV, y::HDV) where {HDV<:BinaryHDV} = sum(xᵢ & yᵢ for (xᵢ, yᵢ) in zip(x, y)) / sum(xᵢ | yᵢ for (xᵢ, yᵢ) in zip(x, y))


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iotdeviceadvisor
using AWS.Compat
using AWS.UUIDs

"""
    create_suite_definition()
    create_suite_definition(params::Dict{String,<:Any})

Creates a Device Advisor test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionConfiguration"`: Creates a Device Advisor test suite with suite
  definition configuration.
- `"tags"`: The tags to be attached to the suite definition.
"""
create_suite_definition(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions"; aws_config=aws_config)
create_suite_definition(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions", params; aws_config=aws_config)

"""
    delete_suite_definition(suite_definition_id)
    delete_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Deletes a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite to be deleted.

"""
delete_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
delete_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)

"""
    get_suite_definition(suite_definition_id)
    get_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Gets information about a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite to get.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionVersion"`: Suite definition version of the test suite to get.
"""
get_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
get_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)

"""
    get_suite_run(suite_definition_id, suite_run_id)
    get_suite_run(suite_definition_id, suite_run_id, params::Dict{String,<:Any})

Gets information about a Device Advisor test suite run.

# Arguments
- `suite_definition_id`: Suite definition Id for the test suite run.
- `suite_run_id`: Suite run Id for the test suite run.

"""
get_suite_run(suiteDefinitionId, suiteRunId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)"; aws_config=aws_config)
get_suite_run(suiteDefinitionId, suiteRunId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)", params; aws_config=aws_config)

"""
    get_suite_run_report(suite_definition_id, suite_run_id)
    get_suite_run_report(suite_definition_id, suite_run_id, params::Dict{String,<:Any})

Gets a report download link for a successful Device Advisor qualifying test suite run.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite.
- `suite_run_id`: Suite run Id of the test suite run.

"""
get_suite_run_report(suiteDefinitionId, suiteRunId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/report"; aws_config=aws_config)
get_suite_run_report(suiteDefinitionId, suiteRunId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/report", params; aws_config=aws_config)

"""
    list_suite_definitions()
    list_suite_definitions(params::Dict{String,<:Any})

Lists the Device Advisor test suites you have created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return at once.
- `"nextToken"`: A token used to get the next set of results.
"""
list_suite_definitions(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions"; aws_config=aws_config)
list_suite_definitions(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions", params; aws_config=aws_config)

"""
    list_suite_runs()
    list_suite_runs(params::Dict{String,<:Any})

Lists the runs of the specified Device Advisor test suite. You can list all runs of the
test suite, or the runs of a specific version of the test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return at once.
- `"nextToken"`: A token to retrieve the next set of results.
- `"suiteDefinitionId"`: Lists the test suite runs of the specified test suite based on
  suite definition Id.
- `"suiteDefinitionVersion"`: Must be passed along with suiteDefinitionId. Lists the test
  suite runs of the specified test suite based on suite definition version.
"""
list_suite_runs(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteRuns"; aws_config=aws_config)
list_suite_runs(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteRuns", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags attached to an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The ARN of the IoT Device Advisor resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    start_suite_run(suite_definition_id)
    start_suite_run(suite_definition_id, params::Dict{String,<:Any})

Starts a Device Advisor test suite run.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionVersion"`: Suite definition version of the test suite.
- `"suiteRunConfiguration"`: Suite run configuration.
- `"tags"`: The tags to be attached to the suite run.
"""
start_suite_run(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns"; aws_config=aws_config)
start_suite_run(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns", params; aws_config=aws_config)

"""
    stop_suite_run(suite_definition_id, suite_run_id)
    stop_suite_run(suite_definition_id, suite_run_id, params::Dict{String,<:Any})

Stops a Device Advisor test suite run that is currently running.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite run to be stopped.
- `suite_run_id`: Suite run Id of the test suite run to be stopped.

"""
stop_suite_run(suiteDefinitionId, suiteRunId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/stop"; aws_config=aws_config)
stop_suite_run(suiteDefinitionId, suiteRunId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/stop", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds to and modifies existing tags of an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The resource ARN of an IoT Device Advisor resource.
- `tags`: The tags to be attached to the IoT Device Advisor resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The resource ARN of an IoT Device Advisor resource.
- `tag_keys`: List of tag keys to remove from the IoT Device Advisor resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_suite_definition(suite_definition_id)
    update_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Updates a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Suite definition Id of the test suite to be updated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionConfiguration"`: Updates a Device Advisor test suite with suite
  definition configuration.
"""
update_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("PATCH", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
update_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("PATCH", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)


text:
### CIndex

"""
    CIndex

Data holding column names and map to linear indices for RowTable.
"""
struct CIndex
    names::Vector{Symbol}
    smap::Dict{Symbol, Int}
end

"""
    Base.names(c::CIndex) = c.names

Return the names of the index `c`.
"""
Base.names(c::CIndex) = c.names

@inline _names(c::CIndex) = c.names
@inline _smap(c::CIndex) = c.smap

### CIndex constructors

CIndex() = CIndex(Symbol[])
CIndex(names::Tuple, smap) = CIndex([names...], smap)
CIndex(names::Tuple) = CIndex([names...])

"""
    CIndex(names::Vector{Symbol})::CIndex

Create a `Cindex` with names and their ordinal position as in `names`.
"""
CIndex(names::Vector{Symbol}) = CIndex(names, Dict(s => i for (i, s) in enumerate(names)))

# Hack around unknown bug that shows up only in v0.7
CIndex(names::Set{Symbol}) = CIndex([x for x in names])

Base.:(==)(c1::CIndex, c2::CIndex) = (_names(c1) == _names(c2))

@inline Base.getindex(c::CIndex, inds) = inds
@inline Base.getindex(c::CIndex, s::Symbol) = _smap(c)[s]  # returns Int
# can we avoid allocation by returning an iterator ?
@inline Base.getindex(c::CIndex, syms::AbstractVector{T}) where T <: Symbol = [_smap(c)[s] for s in syms]
@inline Base.getindex(c::CIndex, inds::AbstractVector{T}) where T <: Integer = inds
@inline Base.getindex(c::CIndex, inds::AbstractVector) = [isa(s, Symbol) ?  _smap(c)[s] : s  for s in inds]

### Copy

@inline Base.copy(ci::CIndex) = CIndex(copy(ci.names), copy(ci.smap))
Base.deepcopy(ci::CIndex) = CIndex(deepcopy(_names(ci)), deepcopy(ci.smap))

### Rename

DataFrames.rename!(c::CIndex, a::Array{T}) where T <: Pair = rename!(c::CIndex, Dict(a))

function DataFrames.rename!(c::CIndex, d::AbstractDict)
    newnames = copy(_names(c))
    for (from, to) in d
        haskey(_smap(c), from) || throw(ErrorException("There is no existing name $from"))
        newnames[_smap(c)[from]] =  to
    end
    length(newnames) == length(unique(newnames)) || throw(ArgumentError("names must be unique"))
 @inbounds for (i, k) in enumerate(newnames)
        _names(c)[i] = k
        _smap(c)[k] = i
    end
end


text:
using Test
using Libxc

# Wrap in an outer testset to get a full report if one test fails
@testset "Libxc.jl" begin

@testset "Metadata" begin
    @test Libxc.libxc_version ≥ v"5.1.0"
    @test Libxc.libxc_doi isa AbstractString
end

@testset "Functional list" begin
    available = available_functionals()
    @test :lda_x     in available
    @test :lda_c_vwn in available
    @test :lda_xc_teter93 in available
    @test :gga_x_pbe in available
    @test :gga_c_pbe in available
end

@testset "Functional construction" begin
    lda = Functional(:lda_x, n_spin=2)
    @test is_lda(lda)
    @test !is_gga(lda)
    @test !is_mgga(lda)
    @test lda.family == :lda
    @test sort(lda.flags) == sort([:vxc, :dim3, :fxc, :exc])
    @test lda.kind == :exchange
    @test lda.name == "Slater exchange"
    @test lda.n_spin == 2
    @test all(i in supported_derivatives(lda) for i in 0:2)
end

@testset "LDA and GGA evaluate" begin
    rho = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma = [0.2, 0.3, 0.4, 0.5, 0.6]
    result = zeros(Float64, 5)

    # LDA
    lda_x = Functional(:lda_x)
    result = evaluate(lda_x, rho=rho, derivatives=0).zk
    @test result ≈ [-0.342809, -0.431912, -0.494416, -0.544175, -0.586194] atol=1e-5

    # GGA
    gga_x = Functional(:gga_x_pbe)
    result = evaluate(gga_x, rho=rho, sigma=sigma, derivatives=0).zk
    @test result ≈ [-0.452598, -0.478878, -0.520674, -0.561428, -0.598661] atol=1e-5
end

@testset "LDA and GGA evaluate!" begin
    rho = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma = [0.2, 0.3, 0.4, 0.5, 0.6]
    result = zeros(Float64, 5)

    # LDA
    func = Functional(:lda_x)
    evaluate!(func, rho=rho, zk=result)
    @test result ≈ [-0.342809, -0.431912, -0.494416, -0.544175, -0.586194] atol=1e-5

    # GGA
    func = Functional(:gga_x_pbe)
    evaluate!(func, rho=rho, sigma=sigma, zk=result)
    @test result ≈ [-0.452598, -0.478878, -0.520674, -0.561428, -0.598661] atol=1e-5
end

@testset "LDA / GGA evaluate! with spin" begin
    rho = abs.(randn(5))
    sigma = abs.(randn(5))

    # Duplicate rho and sigma for spin = 2 tests
    rho2 = zeros(Float64, 2length(rho))
    rho2[1:2:end] = rho2[2:2:end] = rho ./ 2

    sigma2 = zeros(Float64, 3length(rho))
    sigma2[1:3:end] = sigma2[2:3:end] = sigma2[3:3:end] = sigma ./ 4

    # LSDA
    for sym in (:lda_x, :lda_c_vwn)
        E = zeros(Float64, length(rho))
        func = Functional(sym, n_spin=1)
        evaluate!(func, rho=rho, zk=E)

        E2 = zeros(Float64, length(rho))
        func = Functional(sym, n_spin=2)
        evaluate!(func, rho=rho2, zk=E2)

        @test E ≈ E2
    end

    # GGA
    for sym in (:gga_x_pbe, :gga_c_pbe)
        E = zeros(Float64, length(rho))
        func = Functional(sym, n_spin=1)
        evaluate!(func, rho=rho, sigma=sigma, zk=E)

        E2 = zeros(Float64, length(rho))
        func = Functional(sym, n_spin=2)
        evaluate!(func, rho=rho2, sigma=sigma2, zk=E2)

        @test E ≈ E2
    end
end

@testset "LDA / GGA evaluate with spin" begin
    shape = (2, 3, 4)
    rho   = reshape(abs.(randn(shape)), 1, shape...)
    sigma = reshape(abs.(randn(shape)), 1, shape...)

    # Duplicate rho and sigma for spin = 2 tests
    rho2 = 0.5vcat(rho, rho)
    sigma2 = 0.25vcat(sigma, sigma, sigma)

    # LSDA
    for sym in (:lda_x, :lda_c_vwn)
        res = evaluate(Functional(sym, n_spin=1), rho=rho, zk=zeros(shape))
        @test size(res.zk) == shape
        @test size(res.vrho) == (1, shape...)

        res2 = evaluate(Functional(sym, n_spin=2), rho=rho2)
        @test size(res2.zk) == shape
        @test size(res2.vrho) == (2, shape...)

        @test res.zk ≈ res2.zk
        @test res.vrho[1, :, :, :] ≈ res2.vrho[1, :, :, :]
        @test res.vrho[1, :, :, :] ≈ res2.vrho[2, :, :, :]
    end

    # GGA
    for sym in (:gga_x_pbe, :gga_c_pbe)
        res = evaluate(Functional(sym, n_spin=1), rho=rho, sigma=sigma)
        @test size(res.zk) == shape
        @test size(res.vrho) == (1, shape...)
        @test size(res.vsigma) == (1, shape...)

        res2 = evaluate(Functional(sym, n_spin=2), rho=rho2, sigma=sigma2)
        @test size(res2.zk) == shape
        @test size(res2.vrho) == (2, shape...)
        @test size(res2.vsigma) == (3, shape...)

        @test res.zk ≈ res2.zk
        @test res.vrho[1, :, :, :] ≈ res2.vrho[1, :, :, :]
        @test res.vrho[1, :, :, :] ≈ res2.vrho[2, :, :, :]
        @test 4res.vsigma ≈ sum(res2.vsigma, dims=1)
    end
end

@testset "mGGA evaluate! without Laplacian" begin
    rho    = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma  = [0.2, 0.3, 0.4, 0.5, 0.6]
    tau    = [0.5, 0.4, 0.3, 0.2, 0.1]
    result = zeros(Float64, 5)

    tpss_x = Functional(:mgga_x_tpss, n_spin=1)
    evaluate!(tpss_x, rho=rho, sigma=sigma, tau=tau, zk=result)
    @test result ≈ [-0.4254894, -0.4685985, -0.5341998, -0.6017734, -0.6646443] atol=1e-6
end

@testset "mGGA evaluate! with Laplacian" begin
    rho    = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma  = [0.2, 0.3, 0.4, 0.5, 0.6]
    lapl   = [-0.2, -0.3, -0.4, -0.5, -0.6]
    tau    = [0.5, 0.4, 0.3, 0.2, 0.1]
    vrho   = randn(Float64, 5)
    vsigma = randn(Float64, 5)
    vtau   = randn(Float64, 5)
    vlapl  = randn(Float64, 5)

    # TB09 is a potential functional, construct and check this
    tb09_x = Functional(:mgga_x_tb09, n_spin=1)
    @test all(i in supported_derivatives(tb09_x) for i in 1:2)
    @test !(0 in supported_derivatives(tb09_x))

    evaluate!(tb09_x; rho, sigma, lapl, tau, vrho, vsigma, vtau, vlapl)
    @test vrho ≈ [ 0.07015148, -0.4085103, -0.73564854, -0.97043302, -1.15268102] atol=1e-6
    @test iszero(vsigma)
    @test iszero(vtau)
    @test iszero(vlapl)
end

@testset "evaluate! with invalid arguments" begin
    rho = [0.1, 0.2, 0.3, 0.4, 0.5]
    func = Functional(:lda_x)

    # Derivatives ≥ 3 not compiled into libxc at the moment
    @test_throws ArgumentError evaluate(func, rho=rho, derivatives=0:3)
    @test_throws ArgumentError evaluate!(func, rho=rho, v3rho3=randn(5))

    @test_throws DimensionMismatch evaluate!(func, rho=rho, zk=randn(2))
end

@testset "Custom evaluate! dispatch" begin
    # A dummy version of Slater exchange for Float32
    function Libxc.evaluate!(func::Functional, ::Val{:lda}, rho::Array{Float32};
                             zk::Array{Float32})
        @assert func.identifier == :lda_x
        @. zk = -Float32(3/4) * cbrt(Float32(3/π) * rho)
    end
    rho = Float32[0.1, 0.2, 0.3, 0.4, 0.5]
    result = zeros(Float32, 5)
    func = Functional(:lda_x)

    res = evaluate(func, rho=rho, derivatives=0)
    @test eltype(res.zk) == Float32
    @test eltype(res.zk) == Float32
    @test res.zk ≈ [-0.342809, -0.431912, -0.494416, -0.544175, -0.586194] atol=1e-5
end

@testset "Reference extraction" begin
    lda = Functional(:lda_x)
    @test length(lda.references) == 2
    @test lda.references[1].doi  == "10.1017/S0305004100016108"
    @test lda.references[2].doi  == "10.1007/BF01340281"
end

@testset "Setting thresholds" begin
    lda = Functional(:lda_x)

    lda.density_threshold = 1e-4
    @test lda.density_threshold == 1e-4
    lda.density_threshold = 1e-6
    @test lda.density_threshold == 1e-6

    for field in (:zeta_threshold, :sigma_threshold, :tau_threshold)
        setproperty!(lda, field, 1e-4)
        @test getproperty(lda, field) == 1e-4
        setproperty!(lda, field, 1e-6)
        @test getproperty(lda, field) == 1e-6
    end
end

@testset "Read-only properties" begin
    lda = Functional(:lda_x)

    # LDA has none of these things
    allprops = [:exx_coefficient, :cam_alpha, :cam_beta, :cam_omega, :nlc_b, :nlc_C]
    for s in allprops
        @test isnothing(getproperty(lda, s))
    end

    # Global hybrid functional
    b3lyp = Functional(:hyb_gga_xc_b3lyp)
    @test b3lyp.exx_coefficient == 0.2
    @test isnothing(b3lyp.cam_alpha)
    @test isnothing(b3lyp.cam_beta)
    @test isnothing(b3lyp.nlc_b)

    # SCAN + VV10 correlation functional
    scan_vv10 = Functional(:xc_mgga_c_scan_vv10)
    @test isnothing(scan_vv10.exx_coefficient)
    @test isnothing(scan_vv10.cam_alpha)
    @test scan_vv10.nlc_b == 14.0
    @test scan_vv10.nlc_C == 0.0093
    @test !needs_laplacian(scan_vv10)

    # ωB97 functional
    ωB97 = Functional(:xc_hyb_gga_xc_wb97)
    @test isnothing(ωB97.exx_coefficient)
    @test ωB97.cam_alpha == 1.0
    @test ωB97.cam_beta  == -1.0
    @test ωB97.cam_omega == 0.4
    @test isnothing(ωB97.nlc_b)

    # ωB97x_v functional
    ωB97x_v = Functional(:xc_hyb_gga_xc_wb97x_v)
    @test isnothing(ωB97.exx_coefficient)
    @test ωB97x_v.cam_alpha == 1.0
    @test ωB97x_v.cam_beta  == -(1.0 - 0.167)
    @test ωB97x_v.cam_omega == 0.3
    @test ωB97x_v.nlc_b     == 6.0
    @test ωB97x_v.nlc_C     == 0.01
end
end  # outer wrapper


text:
module atn_coupled_module
export time_series_par
export time_series_eta
export time_series
using Statistics

function infected_c(
    AgentState_1,
    AgentState_2,
    AgentAct_1,
    AgentAct_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
)

    N = length(AgentAct_1)
    AgentStateNext_1 = zeros(Bool, N)
    AgentStateNext_2 = zeros(Bool, N)
    for k = 1:N
        AgentStateNext_1[k] = AgentState_1[k]
        AgentStateNext_2[k] = AgentState_2[k]
    end
    for i = 1:N
        #Change someone sick in group 1 recovers, mu_1
        if AgentState_1[i] == 1
            if rand() < mu_1
                AgentStateNext_1[i] = 0
            end
        end
        #Change someone sick in group 2 recovers, mu_2
        if AgentState_2[i] == 1
            if rand() < mu_2
                AgentStateNext_2[i] = 0
            end
        end

        #If node is active
        if rand() < AgentAct_1[i]
            #it makes m connections
            for j = 1:m
                #choose a random index
                rInt = rand(1:N)
                #probability of staying is beta_1

                if rand() < beta_1
                    if AgentState_2[i] == 1
                        #the active node gets the new neigbor sick in POP_1
                        if ((AgentState_1[rInt] == 0) && (rand() < lambda_1))
                            AgentStateNext_1[rInt] = 1
                        end
                    else
                        #the active node gets sick from the new neighbor node
                        if ((AgentState_1[rInt] == 1) && (rand() < lambda_1))
                            AgentStateNext_1[i] = 1
                        end
                    end
                else
                    #the active node gets the new neigbor sick in POP_2
                    if AgentState_1[i] == 1
                        if ((AgentState_1[rInt] == 0) && (rand() < lambda_2))
                            AgentStateNext_2[rInt] = 1
                        end
                    else
                        if ((AgentState_2[rInt] == 1) && (rand() < lambda_1))
                            AgentStateNext_1[i] = 1
                        end
                    end

                end
            end
        end

        if rand() < AgentAct_2[i]
            #it makes m connections
            for j = 1:m
                #choose a random index
                rInt = rand(1:N)
                #probability of staying is beta_2

                if rand() < beta_2
                    if AgentState_2[i] == 1
                        #the active node gets the new neigbor sick in POP_1
                        if ((AgentState_2[rInt] == 0) && (rand() < lambda_2))
                            AgentStateNext_2[rInt] = 1
                        end
                    else
                        #the active node gets sick from the new neighbor node
                        if ((AgentState_2[rInt] == 1) && (rand() < lambda_2))
                            AgentStateNext_2[i] = 1
                        end
                    end
                else
                    #the active node gets the new neigbor sick in POP_2
                    if AgentState_2[i] == 1
                        if ((AgentState_2[rInt] == 0) && (rand() < lambda_1))
                            AgentStateNext_1[rInt] = 1
                        end
                    else
                        if ((AgentState_1[rInt] == 1) && (rand() < lambda_2))
                            AgentStateNext_2[i] = 1
                        end
                    end

                end
            end

        end
    end
    return AgentStateNext_1, AgentStateNext_2

end

function avg_infected(
    N,
    NTrials,
    m,
    mu,
    eta,
    lambda,
    etaS,
    AgentAct,
    tmax,
    TWindow,
)
    InfectedOverTrials = 0.0
    for nt = 1:NTrials
        AgentState = zeros(Bool, N)
        AgentState[1:Int(0.01 * N)] .= 1
        InfectedOverATrial = 0.0
        for i_time = 1:tmax
            NumInf, AgentState = infected_c(
                AgentState,
                AgentAct,
                m,
                mu,
                lambda,
                etaS,
                eta,
                InfectedOverATrial,
                tmax,
                TWindow,
            )
            if i_time > (tmax - TWindow)
                NumInf = sum(AgentState)
                InfectedOverATrial = InfectedOverATrial + NumInf
            end
        end
        InfectedOverATrial = InfectedOverATrial / TWindow
        InfectedOverTrials = InfectedOverTrials + InfectedOverATrial
    end
    return (InfectedOverTrials / NTrials) / N
end

function avg_infected_par(
    N,
    NTrials,
    m,
    mu,
    grid,
    etaS,
    AgentAct,
    tmax,
    TWindow,
)
    eta = grid[1]
    lambda = grid[2]
    InfectedOverTrials = 0.0
    for nt = 1:NTrials
        AgentState = zeros(Bool, N)
        for k = 1:100
            AgentState[k] = 1
        end
        InfectedOverATrial = 0.0
        for i_time = 1:tmax
            AgentState = infected_c(
                AgentState,
                AgentAct,
                m,
                mu,
                lambda,
                etaS,
                eta,
                InfectedOverATrial,
                tmax,
                TWindow,
            )
            if i_time > (tmax - TWindow)
                NumInf = sum(AgentState)
                InfectedOverATrial = InfectedOverATrial + NumInf
            end
        end
        InfectedOverATrial = InfectedOverATrial / TWindow
        InfectedOverTrials = InfectedOverTrials + InfectedOverATrial
    end
    return (InfectedOverTrials / NTrials) / N
end
function infected_totals(
    tmax,
    agent_act_1,
    agent_act_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
)
    N = length(agent_act_1)
    agent_state_1 = zeros(Bool, N)
    agent_state_2 = zeros(Bool, N)

    Iend_1 = zeros(tmax)
    Iend_2 = zeros(tmax)
    for k = 1:100
        agent_state_1[k] = 1
        # agent_state_2[k] = 1
    end

    for i_time = 1:tmax
        agent_state_1, agent_state_2 = infected_c(
            agent_state_1,
            agent_state_2,
            agent_act_1,
            agent_act_2,
            m,
            mu_1,
            mu_2,
            lambda_1,
            lambda_2,
            beta_1,
            beta_2,
        )
        Iend_1[i_time] = sum(agent_state_1)
        Iend_2[i_time] = sum(agent_state_2)
    end
    return Iend_1, Iend_2
end

function infected_totals(
    tmax,
    agent_act_1,
    agent_act_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
)
    N = length(agent_act_1)
    agent_state_1 = zeros(Bool, N)
    agent_state_2 = zeros(Bool, N)

    # # Iend_1 = zeros(tmax)
    # Iend_2 = zeros(tmax)
    for k = 1:100
        # agent_state_1[k] = 1
        agent_state_2[k] = 1
    end
    t_window = 500
    infected = 0.0
    for i_time = 1:tmax
        agent_state_1, agent_state_2 = infected_c(
            agent_state_1,
            agent_state_2,
            agent_act_1,
            agent_act_2,
            m,
            mu_1,
            mu_2,
            lambda_1,
            lambda_2,
            beta_1,
            beta_2,
        )
        # Iend_1[i_time] = sum(agent_state_1)
        if i_time > (tmax - t_window)
            infected = infected + sum(agent_state_1)
        end
    end
    #
    return infected / t_window
end

function parameter_sweep(
    n_trials,
    agent_act_1,
    agent_act_2,
    mu_1,
    mu_2,
    lambda_1,
    # beta_1,
    grid,
)
    N = length(agent_act_1)
    beta_2 = 1 - grid[1]
    beta_1 = beta_2
    lambda_2 = grid[2]
    m = 5
    tmax = 2000
    total = 0.0
    for k = 1:n_trials
        infected = infected_totals(
            tmax,
            agent_act_1,
            agent_act_2,
            m,
            mu_1,
            mu_2,
            lambda_1,
            lambda_2,
            beta_1,
            beta_2,
        )
        total = total + infected / N
    end
    return total / n_trials
end
end


text:
abstract type Triangle{T <: Real} end

struct Incremental{T} <: Triangle{T}
    claims::Matrix{T}
    function Incremental(claims::Matrix{T}) where T <: Real
        if size(claims)[1] != size(claims)[2]
            error("The number of rows and the number of columns must be the same.")
        elseif siniestros_bajo_diagonal(claims::Matrix{T})
            error("Having claims under the diagonal is not supported.")
        end
        new{T}(claims)
    end
end

struct Accumulated{T} <: Triangle{T}
    claims::Matrix{T}
    function Accumulated(claims::Matrix{T}) where T <: Real
        if size(claims)[1] != size(claims)[2]
            error("The number of rows and the number of columns must be the same.")
        elseif siniestros_bajo_diagonal(claims::Matrix{T})
            error("NHaving claims under the diagonal is not supported.")
        end
        new{T}(claims)
    end
end

struct IncrementalCompleted{T} <: Triangle{T}
    claims::Matrix{T}
    function IncrementalCompleted(claims::Matrix{T}) where T <: Real
        if size(claims)[1] != size(claims)[2]
            error("The number of rows and the number of columns must be the same.")
        end
        new{T}(claims)
    end
end

struct AccumulatedCompleted{T} <: Triangle{T}
    claims::Matrix{T}
    function AccumulatedCompleted(claims::Matrix{T}) where T <: Real
        if size(claims)[1] != size(claims)[2]
            error("The number of rows and the number of columns must be the same.")
        end
        new{T}(claims)
    end
end


text:
using Images, ImageFeatures, CoordinateTransformations, ImageDraw

img1 = Gray.(load("cat-3417184_640.jpg"))
img2 = Gray.(load("cat-3417184_640_watermarked.jpg"))

# rotate img2 around the center and resize it
rot = recenter(RotMatrix(5pi/6), [size(img2)...] .÷ 2)  # a rotation around the center
tform = rot ∘ Translation(-50, -40)
img2 = warp(img2, tform, size(img2))
img2 = imresize(img2, Int.(trunc.(size(img2) .* 0.7)))

features_1 = Features(fastcorners(img1, 12, 0.35));
features_2 = Features(fastcorners(img2, 12, 0.35));

brisk_params = BRISK()
desc_1, ret_features_1 = create_descriptor(img1, features_1, brisk_params);
desc_2, ret_features_2 = create_descriptor(img2, features_2, brisk_params);
matches = match_keypoints(Keypoints(ret_features_1), Keypoints(ret_features_2), desc_1, desc_2, 0.2)

# create
img3 = zeros(Gray, size(img1))
img3[1:size(img2, 1), 1:size(img2, 2)] = img2

grid = hcat(img1, img3)
offset = CartesianIndex(0, size(img1, 2))
map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)

imshow(grid)


text:
using PlotUtils: zscale
using PyCall
using SiriusB

fits = pyimport("astropy.io.fits")
pro = pyimport("proplot")

# rcParams
pro.rc["image.origin"] = "lower"
pro.rc["image.cmap"] = "inferno"
pro.rc["grid"] = false

data_cube = fits.getdata(datadir("epoch_2020nov21", "processed", "2020nov21_sirius-b_cube_calib.fits"))
frame = Float32.(data_cube[169, :, :])

fig, axs = pro.subplots(figwidth="4in")
ax = first(axs)
vmin, vmax = zscale(frame)
ax.imshow(frame; vmin, vmax, colorbar="r")
ax.format(xlabel="x [px]", ylabel="y [px]")
fig.save(figuredir("spike.pdf"))

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
module SimpleRandom
using LinearAlgebra
import Random.randperm

include("RV.jl")

export random_unit_vector, random_subset

"""
`random_unit_vector(d)` returns a random `d`-dimensional unit vector.
"""
function random_unit_vector(d::Int)
  v = randn(d)
  return v / norm(v)
end

"""
`random_subset` is used to create random subsets as follows:

+ `random_subset(A)`: random subset of `A` with each element
chosen with probability 1/2.
+ `random_subset(A,k)`: random `k`-element subset of `A`.
+ `random_subset(n)`: random subset of `1:n`.
+ `random_subset(n,k)`: random `k`-element subset of `1:n`.
"""
function random_subset(A::Union{Set,BitSet})
  T = typeof(A)
  B = T()
  for a in A
    if rand() < 0.5
      push!(B,a)
    end
  end
  return B
end

random_subset(n::Int) = random_subset(Set(1:n))

function random_subset(A::Union{Set,BitSet}, k::Int)
  n = length(A)
  if k<0 || k>n
    error("k = $k is out of range")
  end
  T = typeof(A)
  B = T()
  elements = collect(A)
  p = randperm(n)
  for j=1:k
    push!(B,elements[p[j]])
  end
  return B
end

function random_subset(n::Int, k::Int)
  if n<0 || k<0 || k>n
    error("n = $n and/or k = $k invalid")
  end
  x = randperm(n)
  y = x[1:k]
  return Set(y)
end

export random_choice

"""
`random_choice(weights)` randomly chooses a value from `1` to `n`,
where `n` is the number of elements in `weights`. The probability
that `k` is chosen is proportional to `weights[k]`. The `weights`
must be nonnegative and not all zero.

`random_choice(dict)` choose a random key `k` from `dict` with weight
proportional to `dict[k]`. Thus, `dict` must be of type
`Dict{S, T<:Real}`.
"""
function random_choice(weights::Vector{T}) where {T<:Real}
  vals = cumsum(weights)
  vals /= vals[end]
  idx = rand()
  for k=1:length(vals)
    @inbounds if idx <= vals[k]
      return k
    end
  end
  error("Impropper input")
end

function random_choice(d::Dict{S,T}) where {S,T<:Real}
  ks = collect(keys(d))
  n = length(ks)
  wts = [ d[ks[j]] for j=1:n ]
  idx = random_choice(wts)
  return ks[idx]
end

import Distributions

export binom_rv, poisson_rv, exp_rv

"""
`binom_rv(n,p)` generates a random binomial random value.
`p` defaults to `0.5`.
"""
binom_rv(n::Int,p::Real=0.5) = rand(Distributions.Binomial(n,p))

"""
`poisson_rv(lambda)` generates a Poisson random value with
mean `lamba` (which defaults to `1.0`).
"""
poisson_rv(lambda::Real=1.0) = rand(Distributions.Poisson(lambda))

"""
`exp_rv(theta)` returns an exponential random value with
mean `theta` (which defaults to `1.0`).
"""
exp_rv(theta::Real=1.0) = rand(Distributions.Exponential(theta))

end  # end of module


text:
module Renderer

export respond, html, json, redirect

using Nullables, JSON, HTTP, Reexport, Markdown
using Genie, Genie.Util, Genie.Configuration, Genie.Loggers

@reexport using Genie.Flax

const CONTENT_TYPES = Dict{Symbol,String}(
  :html       => "text/html",
  :plain      => "text/plain",
  :text       => "text/plain",
  :json       => "application/json",
  :js         => "application/javascript",
  :javascript => "application/javascript",
  :xml        => "text/xml",
  :markdown   => "text/markdown"
)
const DEFAULT_CONTENT_TYPE = :html
const ResourcePath = Union{String,Symbol}
const HTTPHeaders = Dict{String,String}

"""
    mutable struct WebResource

Represents a resource that can be resolved by the view layer
"""
mutable struct WebResource
  resource::ResourcePath
  action::ResourcePath
  layout::ResourcePath
end

"""
    mutable struct WebRenderable

Represents an object that can be rendered on the web
"""
mutable struct WebRenderable
  body::String
  content_type::Symbol
  status::Int
  headers::HTTPHeaders
end

WebRenderable(body::String) = WebRenderable(body, DEFAULT_CONTENT_TYPE, 200, HTTPHeaders())
WebRenderable(body::String, content_type::Symbol) = WebRenderable(body, content_type, 200, HTTPHeaders())
WebRenderable(; body::String = "", content_type::Symbol = DEFAULT_CONTENT_TYPE,
                status::Int = 200, headers::HTTPHeaders = HTTPHeaders()) = WebRenderable(body, content_type, status, headers)
function WebRenderable(wr::WebRenderable, status::Int, headers::HTTPHeaders)
  wr.status = status
  wr.headers = headers

  wr
end

"""
"""
function tohtml(resource::ResourcePath, action::ResourcePath;
                  layout::ResourcePath = Genie.config.renderer_default_layout_file, context::Module = @__MODULE__, vars...) :: WebRenderable
  WebRenderable(Flax.html_renderer(resource, action; layout = layout, mod = context, vars...) |> Base.invokelatest)
end
function tohtml(data::String; context::Module = @__MODULE__, layout::Union{ResourcePath,Nothing} = nothing, vars...) :: WebRenderable
  WebRenderable(Flax.html_renderer(data; mod = context, layout = layout, vars...) |> Base.invokelatest)
end
function tohtml(restful_resource::WebResource; context::Module = @__MODULE__, vars...) :: WebRenderable
  WebRenderable(restful_resource.resource, restful_resource.action, layout = restful_resource.layout, context = context, vars...)
end

"""
"""
function html(resource::ResourcePath, action::ResourcePath; layout::ResourcePath = Genie.config.renderer_default_layout_file,
                context::Module = @__MODULE__, status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), vars...) :: HTTP.Response
  WebRenderable(tohtml(resource, action; layout = layout, context = context, vars...), status, headers) |> respond
end
function html(data::String; context::Module = @__MODULE__, status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), layout::Union{ResourcePath,Nothing} = nothing, vars...) :: HTTP.Response
  WebRenderable(tohtml(data; context = context, layout = layout, vars...), status, headers) |> respond
end
function html(data::HTML; context::Module = @__MODULE__, status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), layout::Union{ResourcePath,Nothing} = nothing, vars...) :: HTTP.Response
  html(data.content, context = context, status = status, headers = headers, layout = layout, vars...)
end

### JSON RENDERING ###

"""
Invokes the JSON renderer of the underlying configured templating library.
"""
function tojson(resource::ResourcePath, action::ResourcePath; context::Module = @__MODULE__, vars...) :: WebRenderable
  WebRenderable(Flax.json_renderer(resource, action; mod = context, vars...) |> Base.invokelatest, :json)
end
function tojson(data) :: WebRenderable
  WebRenderable(JSON.json(data), :json)
end

"""
"""
function json(resource::ResourcePath, action::ResourcePath; context::Module = @__MODULE__,
              status::Int = 200, headers::HTTPHeaders = HTTPHeaders(), vars...) :: HTTP.Response
  WebRenderable(tojson(resource, action; mod = context, vars...), status, headers) |> respond
end
function json(data; status::Int = 200, headers::HTTPHeaders = HTTPHeaders()) :: HTTP.Response
  WebRenderable(tojson(data), status, headers) |> respond
end

### REDIRECT RESPONSES ###

"""
Sets redirect headers and prepares the `Response`.
"""
function redirect(location::String, code::Int = 302, headers::HTTPHeaders = HTTPHeaders()) :: HTTP.Response
  headers["Location"] = location
  WebRenderable("Redirecting you to $location", :text, code, headers) |> respond
end
function redirect(named_route::Symbol, code::Int = 302, headers::HTTPHeaders = HTTPHeaders()) :: HTTP.Response
  redirect(Genie.Router.linkto(named_route), code, headers)
end

"""
    hasrequested(content_type::Symbol) :: Bool

Checks wheter or not the requested content type matches `content_type`.
"""
function hasrequested(content_type::Symbol) :: Bool
  task_local_storage(:__params)[:response_type] == content_type
end

### RESPONSES ###

"""
Constructs a `Response` corresponding to the content-type sof the request.
"""
function respond(r::WebRenderable) :: HTTP.Response
  haskey(r.headers, "Content-Type") || (r.headers["Content-Type"] = CONTENT_TYPES[r.content_type])

  HTTP.Response(r.status, [h for h in r.headers], body = r.body)
end
function respond(response::HTTP.Response) :: HTTP.Response
  response
end
function respond(body::String, params::Dict{Symbol,T})::HTTP.Response where {T}
  r = params[:RESPONSE]
  r.data = body

  r |> respond
end
function respond(err::T, content_type::Union{Symbol,String} = Genie.Router.responsetype(), code::Int = 500)::HTTP.Response where {T<:Exception}
  HTTP.Response(code, (isa(content_type, Symbol) ? ["Content-Type" => CONTENT_TYPES[content_type]] : ["Content-Type" => content_type]), body = string(err))
end
function respond(body::String, content_type::Union{Symbol,String} = Genie.Router.responsetype(), code::Int = 200) :: HTTP.Response
  HTTP.Response(code, (isa(content_type, Symbol) ? ["Content-Type" => CONTENT_TYPES[content_type]] : ["Content-Type" => content_type]), body = body)
end
function respond(body, code::Int = 200, headers::HTTPHeaders = HTTPHeaders())
  HTTP.Response(code, [h for h in headers], body = string(body))
end
function respond(f::Function, code::Int = 200, headers::HTTPHeaders = HTTPHeaders())
  respond(f(), code, headers)
end

end


text:
using LatPhysBase
using LatPhysLatticeConstruction
using Test

@testset "LatPhysLatticeConstruction.jl" begin

    #------------------
    # create unitcells
    #------------------

    # in 2d
    unitcell_2d = newUnitcell(
        # Type of the unitcell
        Unitcell{Site{Int64,2}, Bond{Symbol,2}},
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[sqrt(3.0)/2, -0.5],
            Float64[sqrt(3.0)/2, +0.5]
        ],
        # Sites
        Site{Int64,2}[
            newSite(Site{Int64,2}, Float64[0.0, 0.0],         getDefaultLabel(Int64)),
            newSite(Site{Int64,2}, Float64[1/sqrt(3.0), 0.0], getDefaultLabel(Int64))
        ],
        # Bonds
        Bond{Symbol,2}[
            newBond(Bond{Symbol,2}, 1, 2, getDefaultLabel(Symbol), (0, 0)),
            newBond(Bond{Symbol,2}, 1, 2, getDefaultLabel(Symbol), (-1, 0)),
            newBond(Bond{Symbol,2}, 1, 2, getDefaultLabel(Symbol), (0, -1)),
            newBond(Bond{Symbol,2}, 2, 1, getDefaultLabel(Symbol), (0, 0)),
            newBond(Bond{Symbol,2}, 2, 1, getDefaultLabel(Symbol), (1, 0)),
            newBond(Bond{Symbol,2}, 2, 1, getDefaultLabel(Symbol), (0, 1))
        ]
    )

    # in 3d
    unitcell_3d = newUnitcell(
        # Type of the unitcell
        Unitcell{Site{Int64,3}, Bond{Symbol,3}},
        # lattice vectors
        Vector{Float64}[
            Float64[1, 0, 0],
            Float64[0, 1, 0],
            Float64[0, 0, 1]
        ],
        # sites
        Site{Int64,3}[
            newSite(Site{Int64,3}, Float64[0,0,0], getDefaultLabel(Int64))
        ],
        # bonds
        Bond{Symbol,3}[
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (+1,0,0)),
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (-1,0,0)),
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (0,+1,0)),
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (0,-1,0)),
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (0,0,+1)),
            newBond(Bond{Symbol,3}, 1,1, getDefaultLabel(Symbol), (0,0,-1))
        ]
    )

    #------------------
    # periodic lattice
    #------------------
    @testset "Periodic lattices" begin

        # 2d given all extents
        @test_nowarn getLatticePeriodic(unitcell_2d, (2,2))
        # 2d given only 1 extent
        @test_nowarn getLatticePeriodic(unitcell_2d, 3)

        # 3d given all extents
        @test_nowarn getLatticePeriodic(unitcell_3d, (2,2,2))
        # 3d given only 1 extent
        @test_nowarn getLatticePeriodic(unitcell_3d, 3)

    end

    #------------------
    # by bond distance
    #------------------
    @testset "Lattices by bond distance" begin

        # 2d given only the bond distance
        @test_nowarn getLatticeByBondDistance(unitcell_2d, 3)
        # 2d given bond distance and origin
        @test_nowarn getLatticeByBondDistance(unitcell_2d, 3, 2)

        # 3d given only the bond distance
        @test_nowarn getLatticeByBondDistance(unitcell_3d, 3)
        # 3d given bond distance and origin
        @test_nowarn getLatticeByBondDistance(unitcell_3d, 3, 1)

    end

end


text:
# enumerate primes
#https://www.kaggle.com/dserdyuk/julia-primes

function main()
    primes = Int64[]
    i = 2
    while i < 1e7
        found = false
        for p in primes
            if i % p == 0
                found = true
                break
            end
            if p > √i # p > sqrt(i)
                break
            end
        end
        if !found
            push!(primes, i)
        end
        i += 1
    end

    f = open("jloutput.csv", "w")
    write(f, "primes\n")
    for p in primes
        write(f, "$(p)\n")
    end
end

tic()
main()
toc()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# ---
# title: 941. Valid Mountain Array
# id: problem941
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/valid-mountain-array/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, return _`true` if and only if it is a valid
# mountain array_.
#
# Recall that arr is a mountain array if and only if:
#
#   * `arr.length >= 3`
#   * There exists some `i` with `0 < i < arr.length - 1` such that:
#     * `arr[0] < arr[1] < ... < arr[i - 1] < A[i] `
#     * `arr[i] > arr[i + 1] > ... > arr[arr.length - 1]`
#
# ![](https://assets.leetcode.com/uploads/2019/10/20/hint_valid_mountain_array.png)
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [2,1]
#     Output: false
#
#
# **Example 2:**
#
#
#
#     Input: arr = [3,5,5]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: arr = [0,3,2,1]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 104`
#   * `0 <= arr[i] <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end