text:
mutable struct PolarGrid <: AbstractGrid
    rows::Vector{Vector{Cell}}
    wrap::Bool
end

function PolarGrid(ncircles; wrap = true)
    rows = [Cell[] for r in 1:ncircles]
    rowheight = 1 / ncircles
    rows[1] = [Cell(1, 1)]
    for row in 2:ncircles
        radius = (row - 1) / ncircles
        circ = 2π * radius
        prevcount = length(rows[row - 1])
        cellwidth = circ / prevcount
        ratio = round(Int, cellwidth / rowheight)
        ncells = prevcount * ratio
        rows[row] = [Cell(row, col) for col in 1:ncells]
    end
    return PolarGrid(rows, wrap)
end

rows(g::PolarGrid) = g.rows

iswrapped(g::PolarGrid) = g.wrap

Base.size(g::PolarGrid) = length(rows(g)), length(last(rows(g)))

function Base.getindex(g::PolarGrid, r, c)
    roob = r < 1 || r > length(rows(g))
    coob = !iswrapped(g) && (c < 1 || c > length(rows(g)[r]))
    (roob || coob) && return nothing
    return rows(g)[r][mod1(c, length(rows(g)[r]))]
end

function Base.iterate(g::PolarGrid)
    isempty(rows(g)) && return nothing
    isempty(first(rows(g))) && return nothing
    return first(first(rows(g))), (2, 1)
end

Base.iterate(g::PolarGrid, s) = iterate(vcat(rows(g)...), s)

max_stars_repo_name: basp/ariadne

max_issues_repo_name: basp/ariadne

max_forks_repo_name: basp/ariadne


text:
# RECOMMENDED: Use macOS or linux, Gtk on Windows is too slow. I made this on Windows :(
using  Gtk, Graphics, Cairo

# Create canvas on which man is drawn
c = @GtkCanvas(600, 600)

# Render a line on the canvas: No documentation found, but thanks to Casey Kneale!
# Original snippet by Casey Kneale here: https://gist.github.com/caseykneale/49e447f41427cfdcc1efbd681c8f6833
function render_line(ctx::CairoContext, start, finish ; dash = nothing, RGB = (0,0,0))
    @assert(length(start) == 2); @assert(length(finish) == 2); @assert(length(RGB) == 3)
    set_source_rgb(ctx, RGB...);
    if !isa(dash, Nothing)
        set_dash(ctx, dash)
    end
    move_to(ctx, start...);
    line_to(ctx, finish...);
    Gtk.stroke(ctx)
    return nothing
end

# Get input recursively: If error, call same function again
# needs_word specifies if the input required is a word or letter: true if word, or false if letter
function get_input(text_to_print, needs_word, in_same_line)
    in_same_line ? print(text_to_print) : println(text_to_print)
    txtin = lowercase(readline());
    if needs_word && length(split(txtin, " ")) > 1
        println("Please enter valid input (1 word only)")
        get_input(text_to_print, needs_word, in_same_line)
    elseif !needs_word && length(split(txtin, "")) > 1
        println("Please enter valid input (1 letter only)")
        get_input(text_to_print, needs_word, in_same_line)
    else
        return txtin
    end
end

# Draw the initial setup for hangman
@guarded draw(c) do widget
    ctx = getgc(c)
    h = height(c)
    w = width(c)

    render_line(ctx, [200, 80], [200, 30])
    render_line(ctx, [200, 30], [300, 30])
    render_line(ctx, [300, 500], [300, 30])
    render_line(ctx, [250, 500], [350, 500])
end

# Draw the man according to the number of guesses left
function render_man(lefta)
    abc = [([200, 180], [200, 350]), ([200, 180], [150, 280]), ([200, 180], [250, 280]), ([200, 350], [150, 450]), ([200, 350], [250, 450])]
    @guarded draw(c) do widget
        ctx = getgc(c)
        h = height(c)
        w = width(c)

        ind = 6 - lefta

        if ind == 1
            arc(ctx, 200, 130, 50, 0, 2pi)
            render_line(ctx, [200, 80], [200, 80])
        elseif ind > 1
            render_line(ctx, abc[ind-1][1], abc[ind-1][2])
        end
    end
end

# Print start screen
println("""

██   ██   █████   ███    ██   ██████   ███    ███   █████   ███    ██
██   ██  ██   ██  ████   ██  ██        ████  ████  ██   ██  ████   ██
███████  ███████  ██ ██  ██  ██   ███  ██ ████ ██  ███████  ██ ██  ██
██   ██  ██   ██  ██  ██ ██  ██    ██  ██  ██  ██  ██   ██  ██  ██ ██
██   ██  ██   ██  ██   ████   ██████   ██      ██  ██   ██  ██   ████
""")

# Hide previous input line and overwrite with current partially guessed word
function show_word()
    print("\u1b[1F")
    print(join(current_array, " "))
    print("\u1b[0K")
end

# Start up game using get_input
word = get_input("Player 1, please enter your word (word will disappear once entered): ", true, true)

# Open and setup window for man to be drawn
game_window = GtkWindow("Hangman", 600, 600)

g = GtkGrid()
g[1, 1] = c

set_gtk_property!(g, :column_homogeneous, true)
set_gtk_property!(g, :column_spacing, 15)  # introduce a 15-pixel gap between columns
push!(game_window, g)
showall(game_window)

# Create a empty word (not revealed at all), and display it
current_array = ["__" for x in 1:length(word)]
show_word()

# Define number of wrong guesses left
left = 6

# Run game until end game state is reached
while true
    # Make wrong guesses left global so we can use it inside the while loop. If this is removed, it stops working for some reason.
    global left

    # Print the prompt for guessing the letter, and get its input
    println()
    letter = get_input("Player 2, guess a letter in the word: You have $left guesses left!", false, false)
    print("\u1b[1F")

    # Check if letter guessed is in the word
    for (index, check_letter) in enumerate(split(word, ""))
        if letter == check_letter
            current_array[index] = letter
        end
    end

    # If letter is not in word and is not guessed already, reduce number of wrong left chances by one, and render the man again
    if !(letter in current_array)
        left -= 1
        render_man(left)
    end

    # Show updated word with partially guessed letters
    show_word()

    # If game is in ending state, break out of the loop and declare winner
    if !("__" in current_array) || (left == 0)
        break
    end
end

# Declare winner
println()
if left == 0
    println("Player 1 wins! Sorry Player 2, you took the L here. The word was \"$word\"")
else
    println("Player 2 wins! Sorry, Player 1, let's have a harder one next time.")
end

# Create recursive function to ask if player wnats to play again.
function replay_prompt()
    print("Play again? (y/n): ")
    reply = lowercase(readline());
    if reply == "y"
        visible(game_window, false)
        run(`julia $PROGRAM_FILE`)
        exit()
    elseif reply == "n"
        exit()
    else
        println("Please enter a valid response: ")
        replay_prompt()
    end
end

# Run replay_prompt to ask if player wants to replay.
replay_prompt()

ext: jl

max_stars_repo_name: soumitradev/Hangman

max_issues_repo_name: soumitradev/Hangman

max_forks_repo_name: soumitradev/Hangman


text:
include(normpath(@__DIR__, "setup.jl"))

@testset "EscapeAnalysis" begin
    include(normpath(@__DIR__, "EscapeAnalysis.jl"))
end


text:
### This file contains convenience functions for constructing LoopSets.

# function strip_unneeded_const_deps!(ls::LoopSet)
#     for op ∈ operations(ls)
#         if isconstant(op) && iszero(length(reducedchildren(op)))
#             op.dependencies = NODEPENDENCY
#         end
#     end
#     ls
# end

function Base.copyto!(ls::LoopSet, q::Expr)
    q.head === :for || throw("Expression must be a for loop.")
    add_loop!(ls, q, 8)
    # strip_unneeded_const_deps!(ls)
end

function add_ci_call!(q::Expr, @nospecialize(f), args, syms, i, mod = nothing)
    call = if f isa Core.SSAValue
        Expr(:call, syms[f.id])
    else
        Expr(:call, f)
    end
    for arg ∈ @view(args[2:end])
        if arg isa Core.SSAValue
            push!(call.args, syms[arg.id])
        else
            push!(call.args, arg)
        end
    end
    if mod !== nothing # indicates it's `vmaterialize(!)`
        reg_size = Expr(:call, lv(:register_size))
        reg_count = Expr(:call, lv(:available_registers))
        cache_lnsze = Expr(:call, lv(:cache_linesize))
        push!(call.args, Expr(:call, Expr(:curly, :Val, QuoteNode(mod))), reg_size, reg_count, cache_lnsze)
    end
    push!(q.args, Expr(:(=), syms[i], call))
end

function substitute_broadcast(q::Expr, mod::Symbol)
    ci = first(Meta.lower(LoopVectorization, q).args).code
    nargs = length(ci)-1
    ex = Expr(:block,)
    syms = [gensym() for _ ∈ 1:nargs]
    for n ∈ 1:nargs
        ciₙ = ci[n]
        ciₙargs = ciₙ.args
        f = first(ciₙargs)
        if ciₙ.head === :(=)
            push!(ex.args, Expr(:(=), f, syms[((ciₙargs[2])::Core.SSAValue).id]))
        elseif isglobalref(f, Base, :materialize!)
            add_ci_call!(ex, lv(:vmaterialize!), ciₙargs, syms, n, mod)
        elseif isglobalref(f, Base, :materialize)
            add_ci_call!(ex, lv(:vmaterialize), ciₙargs, syms, n, mod)
        else
            add_ci_call!(ex, f, ciₙargs, syms, n)
        end
    end
    ex
end

function LoopSet(q::Expr, mod::Symbol = :Main)
    contract_pass!(q)
    ls = LoopSet(mod)
    copyto!(ls, q)
    resize!(ls.loop_order, num_loops(ls))
    ls
end
LoopSet(q::Expr, m::Module) = LoopSet(macroexpand(m, q)::Expr, Symbol(m))

function loopset(q::Expr) # for interactive use only
    ls = LoopSet(q)
    set_hw!(ls)
    ls
end

"""
    @avx

Annotate a `for` loop, or a set of nested `for` loops whose bounds are constant across iterations, to optimize the computation. For example:

    function AmulBavx!(C, A, B)
        @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ 1:size(A,2)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end

The macro models the set of nested loops, and chooses an ordering of the three loops to
minimize predicted computation time.

It may also apply to broadcasts:

```jldoctest
julia> using LoopVectorization

julia> a = rand(100);

julia> b = @avx exp.(2 .* a);

julia> c = similar(b);

julia> @avx @. c = exp(2a);

julia> b ≈ c
true
```

# Extended help

Advanced users can customize the implementation of the `@avx`-annotated block
using keyword arguments:

```
@avx inline=false unroll=2 body
```

where `body` is the code of the block (e.g., `for ... end`).

`inline` is a Boolean. When `true`, `body` will be directly inlined
into the function (via a forced-inlining call to `_avx_!`).
When `false`, it wont force inlining of the call to `_avx_!` instead, letting Julia's own inlining engine
determine whether the call to `_avx_!` should be inlined. (Typically, it won't.)
Sometimes not inlining can lead to substantially worse code generation, and >40% regressions, even in very
large problems (2-d convolutions are a case where this has been observed).
One can find some circumstances where `inline=true` is faster, and other circumstances
where `inline=false` is faster, so the best setting may require experimentation. By default, the macro
tries to guess. Currently the algorithm is simple: roughly, if there are more than two dynamically sized loops
or and no convolutions, it will probably not force inlining. Otherwise, it probably will.

`check_empty` (default is `false`) determines whether or not it will check if any of the iterators are empty.
If false, you must ensure yourself that they are not empty, else the behavior of the loop is undefined and
(like with `@inbounds`) segmentation faults are likely.

`unroll` is an integer that specifies the loop unrolling factor, or a
tuple `(u₁, u₂) = (4, 2)` signaling that the generated code should unroll more than
one loop. `u₁` is the unrolling factor for the first unrolled loop and `u₂` for the next (if present),
but it applies to the loop ordering and unrolling that will be chosen by LoopVectorization,
*not* the order in `body`.
`uᵢ=0` (the default) indicates that LoopVectorization should pick its own value,
and `uᵢ=-1` disables unrolling for the correspond loop.

The `@avx` macro also checks the array arguments using `LoopVectorization.check_args` to try and determine
if they are compatible with the macro. If `check_args` returns false, a fall back loop annotated with `@inbounds`
and `@fastmath` is generated. Note that `VectorizationBase` provides functions such as `vadd` and `vmul` that will
ignore `@fastmath`, preserving IEEE semantics both within `@avx` and `@fastmath`.
`check_args` currently returns false for some wrapper types like `LinearAlgebra.UpperTriangular`, requiring you to
use their `parent`. Triangular loops aren't yet supported.
"""
macro avx(q)
    q = macroexpand(__module__, q)
    isa(q, Expr) || return q
    q2 = if q.head === :for
        setup_call(LoopSet(q, __module__), q, __source__)
    else# assume broadcast
        substitute_broadcast(q, Symbol(__module__))
    end
    esc(q2)
end

function check_inline(arg)
    a1 = (arg.args[1])::Symbol
    a1 === :inline || return zero(Int8)
    i = (arg.args[2])::Bool % Int8
    i + i - one(Int8)
end
function check_unroll(arg)
    a1 = (arg.args[1])::Symbol
    default = (zero(Int8),zero(Int8))
    a1 === :unroll || return default
    tup = arg.args[2]
    u₂ = -one(Int8)
    if tup isa Integer
        u₁ = convert(Int8, tup)
    elseif isa(tup, Expr)
        if length(tup.args) == 1
            u₁ = convert(Int8, tup.args[1])
        elseif length(tup.args) == 2
            u₁ = convert(Int8, tup.args[1])
            u₂ = convert(Int8, tup.args[2])
        else
            return default
        end
    else
        return default
    end
    u₁, u₂
end
function check_checkempty(arg)
    arg.args[1] === :check_empty ? (arg.args[2])::Bool : nothing
end
function check_macro_kwarg(arg, inline::Int8 = zero(Int8), check_empty::Bool = false, u₁::Int8 = zero(Int8), u₂::Int8 = zero(Int8))
    @assert arg.head === :(=)
    i = check_inline(arg)
    if iszero(i)
        ce = check_checkempty(arg)
        if isnothing(ce)
            u₁, u₂ = check_unroll(arg)
        else
            check_empty = ce
        end
    else
        inline = i
    end
    inline, check_empty, u₁, u₂
end
macro avx(arg, q)
    @assert q.head === :for
    @assert arg.head === :(=)
    q = macroexpand(__module__, q)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg)
    ls = LoopSet(q, __module__)
    esc(setup_call(ls, q, __source__, inline, check_empty, u₁, u₂))
end
macro avx(arg1, arg2, q)
    @assert q.head === :for
    q = macroexpand(__module__, q)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg1)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg2, inline, check_empty, u₁, u₂)
    esc(setup_call(LoopSet(q, __module__), q, __source__, inline, check_empty, u₁, u₂))
end
macro avx(arg1, arg2, arg3, q)
    @assert q.head === :for
    q = macroexpand(__module__, q)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg1)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg2, inline, check_empty, u₁, u₂)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg3, inline, check_empty, u₁, u₂)
    esc(setup_call(LoopSet(q, __module__), q, __source__, inline, check_empty, u₁, u₂))
end

"""
    @_avx

This macro transforms loops similarly to [`@avx`](@ref).
While `@avx` punts to a generated function to enable type-based analysis, `_@avx`
works on just the expressions. This requires that it makes a number of default assumptions. Use of `@avx` is preferred.

This macro accepts the `inline` and `unroll` keyword arguments like `@avx`, but ignores the `check_empty` argument.
"""
macro _avx(q)
    q = macroexpand(__module__, q)
    ls = LoopSet(q, __module__)
    set_hw!(ls)
    esc(Expr(:block, ls.prepreamble, lower_and_split_loops(ls, -1)))
end
macro _avx(arg, q)
    @assert q.head === :for
    q = macroexpand(__module__, q)
    inline, check_empty, u₁, u₂ = check_macro_kwarg(arg)
    ls = LoopSet(q, __module__)
    set_hw!(ls)
    esc(Expr(:block, ls.prepreamble, lower(ls, u₁ % Int, u₂ % Int, -1)))
end

macro avx_debug(q)
    q = macroexpand(__module__, q)
    ls = LoopSet(q, __module__)
    esc(LoopVectorization.setup_call_debug(ls))
end


text:
# OLS estimation using (\)
# https://github.com/giob1994/Alistair.jl/blob/3a11c19150169695581b46e4d1895f0641a4c29d/src/linregress.jl#L38-L41
function ols(X, Y)
    # Ax = B
    # x = B \ A
    # https://github.com/JuliaLang/julia/blob/d789231e9985537686052db9b2314c0d51656308/stdlib/LinearAlgebra/src/generic.jl#L827-L855
    return (X' * X) \ (X' * Y)
end

# OLS estimation using inv()
# inv() is a more computationally-intensive procedure than \ operator
# https://github.com/mcreel/Econometrics/blob/b29cffa8ce7f845fa73788522f9c7753d40684f4/Examples/OLS/OlsFit.jl#L13-L14
function ols_(X, Y)
    return inv(X'*X)*X'*Y
end


text:
"""
    ILP(optimizer, max_iter)

ILP iteratively solves a linearized primal optimization to compute maximum allowable disturbance.
It iteratively adds the linear constraint to the problem.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: halfspace

# Return
`AdversarialResult`

# Method
Iteratively solve a linear encoding of the problem.
It only considers the linear piece of the network that has the same activation pattern as the reference input.
Default `optimizer` is `GLPKSolverMIP()`.
We provide both iterative method and non-iterative method to solve the LP problem.
Default `iterative` is `true`.

# Property
Sound but not complete.

# Reference
[O. Bastani, Y. Ioannou, L. Lampropoulos, D. Vytiniotis, A. Nori, and A. Criminisi,
"Measuring Neural Net Robustness with Constraints,"
in *Advances in Neural Information Processing Systems*, 2016.](https://arxiv.org/abs/1605.07262)
"""
@with_kw struct ILP{O<:AbstractMathProgSolver}
    optimizer::O    = GLPKSolverMIP()
    iterative::Bool = true
end

function solve(solver::ILP, problem::Problem)
    nnet = problem.network
    x = problem.input.center
    model = Model(solver = solver.optimizer)
    δ = get_activation(nnet, x)
    neurons = init_neurons(model, nnet)
    add_complementary_set_constraint!(model, problem.output, last(neurons))
    o = max_disturbance!(model, first(neurons) - problem.input.center)

    if !solver.iterative
        encode_lp!(model, nnet, neurons, δ)
        status = solve(model, suppress_warnings = true)
        status != :Optimal && return AdversarialResult(:Unknown)
        return interpret_result(solver, getvalue(o), problem.input)
    end

    encode_relaxed_lp!(model, nnet, neurons, δ)
    while true
        status = solve(model, suppress_warnings = true)
        status != :Optimal && return AdversarialResult(:Unknown)
        x = getvalue(first(neurons))
        matched, index = match_activation(nnet, x, δ)
        if matched
            return interpret_result(solver, getvalue(o), problem.input)
        end
        add_constraint!(model, nnet, neurons, δ, index)
    end
end

function interpret_result(solver::ILP, o, input)
    if o >= maximum(input.radius)
        return AdversarialResult(:SAT, o)
    else
        return AdversarialResult(:UNSAT, o)
    end
end

function add_constraint!(model::Model,
                         nnet::Network,
                         z::Vector{Vector{Variable}},
                         δ::Vector{Vector{Bool}},
                         index::Tuple{Int64, Int64})
    i, j = index
    layer = nnet.layers[i]
    val = layer.weights[j, :]' * z[i] + layer.bias[j]
    if δ[i][j]
        @constraint(model, val >= 0.0)
    else
        @constraint(model, val <= 0.0)
    end
end

function match_activation(nnet::Network, x::Vector{Float64}, δ::Vector{Vector{Bool}})
    curr_value = x
    for (i, layer) in enumerate(nnet.layers)
        curr_value = layer.weights * curr_value + layer.bias
        for (j, val) in enumerate(curr_value)
            act = δ[i][j]
            if act && val < -0.0001 # Floating point operation
                return (false, (i, j))
            end
            if !act && val > 0.0001 # Floating point operation
                return (false, (i, j))
            end
        end
        curr_value = layer.activation(curr_value)
    end
    return (true, (0, 0))
end


text:
function makeH_Ising(Nx::Int, Ny::Int, J::Real, hx, hz; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    X = ITensor(s, s')
    X[s(1), s'(2)] = 0.5
    X[s(2), s'(1)] = 0.5

    hxs = hx isa Matrix ? hx : fill(hx, Ny, Nx)
    hzs = hz isa Matrix ? hz : fill(hz, Ny, Nx)

    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            if J != 0.0
                op_a  = J * copy(Z)
                op_b  = copy(Z)
                sites = [row=>col, row+1=>col]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))
            end
        end
        if col < Nx
            if J != 0.0
                op_a  = J * Z
                op_b  = copy(Z)
                sites = [row=>col, row=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))
            end
        end
        if hzs[ row, col ] != 0.0
            push!(H[row, col], Operator([row=>row], [hzs[row,col]*copy(Z)], s, Field))
        end
        if hxs[ row, col ] != 0.0
            push!(H[row, col], Operator([row=>row], [hxs[row,col]*copy(X)], s, Field))
        end
    end
    return H
end

function makeCuH_Ising(Nx::Int, Ny::Int, J::Real, hx, hz; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    X = ITensor(s, s')
    X[s(1), s'(2)] = 0.5
    X[s(2), s'(1)] = 0.5

    hxs = hx isa Matrix ? hx : fill(hx, Ny, Nx)
    hzs = hz isa Matrix ? hz : fill(hz, Ny, Nx)

    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            if J != 0.0
                op_a  = -J * Z
                op_b  = copy(Z)
                sites = [row=>col, row+1=>col]
                push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Vertical))
            end
        end
        if col < Nx
            if J != 0.0
                op_a  = -J * Z
                op_b  = copy(Z)
                sites = [row=>col, row=>col+1]
                push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Horizontal))
            end
        end
        if hzs[ row, col ] != 0.0
            push!(H[row, col], Operator([row=>col], [cuITensor(hzs[row,col]*copy(Z))], s, Field))
        end
        if hxs[ row, col ] != 0.0
            push!(H[row, col], Operator([row=>col], [cuITensor(hxs[row,col]*copy(X))], s, Field))
        end
    end
    return H
end

function makeH_XXZ(Nx::Int, Ny::Int, J::Real; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            op_a  = 0.5 * P
            op_b  = copy(M)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            op_a  = 0.5 * M
            op_b  = copy(P)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            if J != 0.0
                op_a  = J * Z
                op_b  = copy(Z)
                sites = [row=>col, row+1=>col]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))
            end
        end
        if col < Nx
            op_a  = 0.5 * P
            op_b  = copy(M)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            op_a  = 0.5 * M
            op_b  = copy(P)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            if J != 0.0
                op_a  = J * Z
                op_b  = copy(Z)
                sites = [row=>col, row=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))
            end
        end
    end
    # pinning fields
    J_ = 1.;
    @inbounds for row in 1:Ny
        op = isodd(row-1) ? (J_/2.0) * Z : (-J_/2.0) * Z
        push!(H[row, 1], Operator([row=>1], [op], s, Field))
        op = isodd(row-1) ? (-J_/2.0) * Z : (J_/2.0) * Z
        push!(H[row, Nx], Operator([row=>Nx], [op], s, Field))
    end
    @inbounds for col in 2:Nx-1
        op = isodd(col-1) ? (-J_/2.0) * Z : (J_/2.0) * Z
        push!(H[Ny, col], Operator([Ny=>col], [op], s, Field))
        op = isodd(col-1) ? (J_/2.0) * Z : (-J_/2.0) * Z
        push!(H[1, col], Operator([1=>col], [op], s, Field))
    end
    return H
end

function makeH_J1J2(Nx::Int, Ny::Int, J1::Real, J2::Real; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            op_a  = 0.5 * J1 * P
            op_b  = copy(M)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            op_a  = 0.5 * J1 * M
            op_b  = copy(P)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            op_a  = J1 * Z
            op_b  = copy(Z)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))
        end
        if col < Nx
            op_a  = 0.5 * J1 * P
            op_b  = copy(M)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            op_a  = 0.5 * J1 * M
            op_b  = copy(P)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            op_a  = J1 * Z
            op_b  = copy(Z)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))
        end
        # handle J2
        if col < Nx
            if row < Ny
                op_a  = 0.5 * J2 * P
                op_b  = copy(M)
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = 0.5 * J2 * M
                op_b  = copy(P)
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = J2 * Z
                op_b  = copy(Z)
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))
            end
            if row > 1
                op_a  = 0.5 * J2 * P
                op_b  = copy(M)
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = 0.5 * J2 * M
                op_b  = copy(P)
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = J2 * Z
                op_b  = copy(Z)
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))
            end
        end
    end
    return H
end

function makeCuH_J1J2(Nx::Int, Ny::Int, J1::Real, J2::Real; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            op_a  = cuITensor(0.5 * J1 * P)
            op_b  = cuITensor(copy(M))
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            op_a  = cuITensor(0.5 * J1 * M)
            op_b  = cuITensor(copy(P))
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))

            op_a  = cuITensor(J1 * Z)
            op_b  = cuITensor(copy(Z))
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Vertical))
        end
        if col < Nx
            op_a  = cuITensor(0.5 * J1 * P)
            op_b  = cuITensor(copy(M))
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            op_a  = cuITensor(0.5 * J1 * M)
            op_b  = cuITensor(copy(P))
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))

            op_a  = cuITensor(J1 * Z)
            op_b  = cuITensor(copy(Z))
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [op_a; op_b], s, Horizontal))
        end
        # handle J2
        if col < Nx
            if row < Ny
                op_a  = cuITensor(0.5 * J2 * P)
                op_b  = cuITensor(copy(M))
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = cuITensor(0.5 * J2 * M)
                op_b  = cuITensor(copy(P))
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = cuITensor(J2 * Z)
                op_b  = cuITensor(copy(Z))
                sites = [row=>col, row+1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))
            end
            if row > 1
                op_a  = cuITensor(0.5 * J2 * P)
                op_b  = cuITensor(copy(M))
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = cuITensor(0.5 * J2 * M)
                op_b  = cuITensor(copy(P))
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))

                op_a  = cuITensor(J2 * Z)
                op_b  = cuITensor(copy(Z))
                sites = [row=>col, row-1=>col+1]
                push!(H[row, col], Operator(sites, [op_a; op_b], s, Diag))
            end
        end
    end
    return H
end

function makeCuH_XXZ(Nx::Int, Ny::Int, J::Real; pinning::Bool=false)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    #Ident = spinI(s)
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    H = Matrix{Vector{Operator}}(undef, Ny, Nx)
    @inbounds for col in 1:Nx, row in 1:Ny
        H[row, col] = Vector{Operator}()
        if row < Ny
            op_a  = 0.5 * P
            op_b  = copy(M)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Vertical))

            op_a  = 0.5 * M
            op_b  = copy(P)
            sites = [row=>col, row+1=>col]
            push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Vertical))

            if J != 0.0
                op_a  = J * Z
                op_b  = copy(Z)
                sites = [row=>col, row+1=>col]
                push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Vertical))
            end
        end
        if col < Nx
            op_a  = 0.5 * P
            op_b  = copy(M)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Horizontal))

            op_a  = 0.5 * M
            op_b  = copy(P)
            sites = [row=>col, row=>col+1]
            push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Horizontal))

            if J != 0.0
                op_a  = J * Z
                op_b  = copy(Z)
                sites = [row=>col, row=>col+1]
                push!(H[row, col], Operator(sites, [cuITensor(op_a); cuITensor(op_b)], s, Horizontal))
            end
        end
    end
    # pinning fields
    J_ = 1.
    @inbounds for row in 1:Ny
        op = isodd(row-1) ? (J_/2.0) * Z : (-J_/2.0) * Z
        push!(H[row, 1], Operator([row=>1], [cuITensor(op)], s, Field))
        op = isodd(row-1) ? (-J_/2.0) * Z : (J_/2.0) * Z
        push!(H[row, Nx], Operator([row=>Nx], [cuITensor(op)], s, Field))
    end
    @inbounds for col in 2:Nx-1
        op = isodd(col-1) ? (-J_/2.0) * Z : (J_/2.0) * Z
        push!(H[Ny, col], Operator([Ny=>col], [cuITensor(op)], s, Field))
        op = isodd(col-1) ? (J_/2.0) * Z : (-J_/2.0) * Z
        push!(H[1, col], Operator([1=>col], [cuITensor(op)], s, Field))
    end
    return H
end


text:
# This contains several objectives

##
# TODO make objectives have ! names
##
# Objective: L∞ norm of the disturbance
function max_disturbance!(model::Model, var)
    o = symbolic_infty_norm(var)
    @objective(model, Min, o)
    return o
end

function min_sum!(model::Model, var)
    o = sum(sum.(var))
    @objective(model, Min, o)
    return o
end

function max_sum!(model::Model, var)
    o = sum(sum.(var))
    @objective(model, Max, o)
    return o
end

function zero_objective!(model::Model)
    @objective(model, Max, 0.0)
end

function linear_objective!(mode::Model, map::HPolytope, var)
    c, d = tosimplehrep(map)
    o = c * var - d
    @objective(model, Min, o)
    return o
end


text:
const push_format! = PetscViewerPushFormat
const set_mode! = PetscViewerFileSetMode
const set_name! = PetscViewerFileSetName
const set_type! = PetscViewerSetType

"""
    PetscViewer(comm::MPI.Comm, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii", mode::PetscFileMode = FILE_MODE_WRITE)

Constructor for a PetscViewer intended to read/write a matrix or a vector with the supplied type and format.
"""
function PetscViewer(comm::MPI.Comm, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii", mode::PetscFileMode = FILE_MODE_WRITE)
    viewer = PetscViewerCreate(comm)
    set_type!(viewer, type)
    set_mode!(viewer, FILE_MODE_WRITE)
    push_format!(viewer, format)
    set_name!(viewer, filename)
    return viewer
end

destroy!(viewer::PetscViewer) = PetscViewerDestroy(viewer)

max_stars_repo_name: bmxam/JuliaPetsc

max_issues_repo_name: bmxam/JuliaPetsc

max_forks_repo_name: bmxam/JuliaPetsc


text:
@time begin
    @time using GaussianBeamLens
    @time using GR

    @time X = GaussianBeamLens.GaussianBeamLensPropagate(0, 25, 1000, 1000);
    # See https://gr-framework.org/workstations.html
    ENV["GKSwstype"] = 41 # Use quick native Win32 GUI
    @time GR.init()
    @time w, h = size(X)
    @time GR.setwindow(0, w, 0, h)
    @time GR.setcolormap(GR.COLORMAP_VIRIDIS)
    @time cmap = colormap();
    @time cmin, cmax = extrema(X)
    #@time data = map(x -> GR.jlgr.normalize_color(x, cmin, cmax), X)
    #@time rgba = [GR.jlgr.to_rgba(value, cmap) for value = data]
    data = GR.jlgr.normalize_color.(X, cmin, cmax)
    rgba = GR.jlgr.to_rgba.(data, Ref(cmap))
    @time GR.drawimage(0, w, 0, h, w, h, rgba)
    @time GR.updatews()
    println("Total time:") #  2.238569 seconds
end

max_stars_repo_name: mkitti/GaussianBeamLens

max_issues_repo_name: mkitti/GaussianBeamLens

max_forks_repo_name: mkitti/GaussianBeamLens


text:
using AxisKeys
using LinearAlgebra
using MAT
using NamedDims
using Plots
using Statistics

"""
    load_besa_av(file_name)

Produces a dictionary of arrays (using AxisKeys and NamedDims), based on the input .mat file
of a single average. The input .mat file is produced from the BESA®-MATLAB® interface and
then saved as a separate .mat file containing the average of a single condition.

Each array has the following format: [time(ms), channels]

"""
function load_besa_av(file_name)
    ## Loading the data of all epochs
    data_file = matread(file_name)
    dict_content = collect(keys(data_file))
    if length(dict_content) > 1
        error("Contents of file contain more than a single struct")
        return
    end

    data = data_file[dict_content[1]]
    averaged_data  = data["data"]["amplitudes"]

    averaged_data_array= wrapdims(
        averaged_data,
        time = dropdims(data["data"]["latencies"],dims=1),
        #WARNING: last channel called REF is being indexed out
        channels = Symbol.(data["channellabels"][1:351])
    )

    return averaged_data_array

end

"""
    load_cont_epochs(file_name)

Produces a dictionary of arrays (using AxisKeys and NamedDims), based on the different
conditions classified from the input .mat file. The input .mat file is produced from the
BESA®-MATLAB® interface when the "Epochs around triggers" option is selected.

Each array has the following format: [time(ms), channels, trials]

"""
function load_cont_epochs(file_name)
## Loading the data of all epochs
cont_epochs = matread(file_name)
cont_epochs = cont_epochs["besa_channels"]
all_epochs  = cont_epochs["data"]["amplitudes"]

## Making sure all epochs are the same length
one_latency = cont_epochs["data"]["latencies"][1]
all_latencies = dropdims(cont_epochs["data"]["latencies"],dims=1)
test_all = map(x -> x==one_latency, all_latencies)
if false ∈ test_all
throw("Input Error: all epochs are not of the same size")
return
end

## Extracting all the conditons via their trigger codes

# Proallocating total number of conditions
stacked_conditions = Array{String,1}(undef, length(all_epochs))

# Going through all epochs to determine the type of event at t=0

for epoch = 1:length(all_epochs)
    # determine index at t=0 in each event
    epoch_latencies  = cont_epochs["data"]["event"][epoch]["latency"]

    if length(epoch_latencies) > 1
        t_zero   = findfirst(isequal(0),epoch_latencies[:])
    else
        t_zero   = findfirst(isequal(0),epoch_latencies)
    end

    # Get the condition label/trigger at t=0
    stacked_conditions[epoch]=string(cont_epochs["data"]["event"][epoch]["label"][t_zero][1])
end

# Getting all the different conditions and their counts so they can be nested (under a subject)
# by name. The count is so that we can preallocate the number of trials
unique_conditions = unique(stacked_conditions)
unique_counts = Dict([(condition,count(x->x==condition, stacked_conditions))
for condition in unique_conditions])

# Determine size of a specific epoch
epoch_dims = size(cont_epochs["data"]["amplitudes"][1])

# Making a dictionary to contain all epochs in each condition
# Preallocating based on known details
# DOC: epoch size of all conditions has to be the same
stacked_epochs = Dict()
for (condition,unique_count) in unique_counts
    stacked_epochs[condition] = Array{Float64,3}(undef,epoch_dims[1],epoch_dims[2],unique_count )
end
for condition in unique_conditions

    for epoch = 1:length(unique_counts[condition])
        # Placing the epoch data based on the condition
        condition_index = condition.==stacked_conditions
        relevant_epochs = all_epochs[condition_index]
        for trial = 1:length(relevant_epochs)
            stacked_epochs[condition][:,:,trial] = relevant_epochs[trial]
        end
        # Adding named dimentions. Getting the timing from from the first epoch, i.e
        # assuming that all epoch sizes are the same
         stacked_epochs[condition] = wrapdims(
            stacked_epochs[condition],
            time = dropdims(cont_epochs["data"]["latencies"][1],dims=1),
            channels = dropdims(Symbol.(cont_epochs["channellabels"]),dims=1),
            trials = 1:(size(stacked_epochs["1"],3)),
            )

    end
end

return stacked_epochs

end

"""
    load_BSepochs(subject_path::String)

Reads epoch data (trials) from the outputs generated by Brainstorm and bundles them together.
The general output of Brainstorm consists of .mat data files for each different epoch
(trial`number`), in each different condition (data`condition_label`). An example file name
is `data_15_trial045.mat`. The output is a Dict type with labels for each condition, with
data being present as KeyedArrays.

For more information on Brainstorm:
Tadel F, Baillet S, Mosher JC, Pantazis D, Leahy RM (2011). Brainstorm: A User-Friendly
Application for MEG/EEG Analysis Computational Intelligence and Neuroscience, vol. 2011,
ID 879716. Website: https://neuroimage.usc.edu/brainstorm/

# Examples
```julia-repl
julia> load_BSepochs("reg_soi_ext1_ya10_sss")
Dict{Any,Any} with 19 entries:
  "4"  => [showing 3 of 100 slices]…
  "1"  => [showing 3 of 100 slices]…
  "12" => [showing 3 of 100 slices]…
  "2"  => [showing 3 of 100 slices]…
  ⋮    => ⋮
```
"""
function load_BSepochs(subject_path::String)

    # Switch working directory to the location of the subjects data
    # TODO: Correct to actual file structure as designated by Brainstorm
    path = subject_path

    # Collecting all unique conditions present
    dir_contents = readdir(path)

    # Getting channel Data
    channel_labels = matread(joinpath(path, "channel_vectorview306_acc1.mat"))["Channel"]["Name"]

    # Getting names of bad trials
    bad_trials = matread(joinpath(path, "brainstormstudy.mat"))["BadTrials"]

    # Getting epoch files
    epoch_files = filter(x -> x[1:4] == "data", dir_contents)
    epoch_conditions = Vector{String}(undef, length(epoch_files))
    for (idx, file) in enumerate(epoch_files)

        condition = split(file, "_")[2] # Get condition label
        epoch_conditions[idx] = condition
    end

    unique_conditions = unique(epoch_conditions)

    # Going through the conditions (one by one) and collecting the data
    data = Dict()

    for condition in unique_conditions

        # WARNING: This assumes that there are not more than 99 epoch_conditions
        to_read = filter( x-> split(x, "_")[2] == condition, epoch_files)
        trial_full = matread(joinpath(path, to_read[1])) # Contains all sorts of metadata + signal data
        one_trial = trial_full["F"]' # Reading signal data to get expected dimentions
        total_timepoints = size(one_trial, 1)
        n_channels = size(one_trial, 2)
        # Removing bad trials
        n_trials = length(to_read) - length(intersect(bad_trials, to_read))
        stacked_epochs = Array{Float64}(undef, total_timepoints, n_channels, n_trials)
        clean_trials = setdiff(to_read, bad_trials)

        # Bundling together all the trials in a single array
        for (idx, trial) in enumerate(clean_trials)

            trial_data = matread(joinpath(path, trial))["F"]' # F is the label given to the data
            stacked_epochs[:,:,idx] = trial_data .* 10^15 # Converting from tesla to femto tesla

        end

        # Converting the array into an AxisKeys Array for better indexing
        stacked_epochs = wrapdims(
           stacked_epochs,
           # Time is converted to ms (from seconds) to ms
           # This also means that epoch selection should always be integers in proprocessing
           time = round.(Int, trial_full["Time"][:] * 1000),
           channels = dropdims(Symbol.(channel_labels), dims=1),
           trials = 1:n_trials,
           )

        # Assigning the condition data to the central data dictionary in order to keep all
        # conditions accessible from a single subject
        data[condition] = stacked_epochs

    end

    return data

end

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text:
using Test, Psycho

for file in readdir("tests/")
    include("tests/" * file)
end


text:
@testset "vertices" begin
    boundary = HyperRectangle(SVector(-1., -2, -3), SVector(2., 4, 6))
    v = vertices(boundary)
    @test v[1,1,1] == [-1, -2, -3]
    @test v[2,1,1] == [1, -2, -3]
    @test v[1,2,1] == [-1, 2, -3]
    @test v[2,2,1] == [1, 2, -3]
    @test v[1,1,2] == [-1, -2, 3]
    @test v[2,1,2] == [1, -2, 3]
    @test v[1,2,2] == [-1, 2, 3]
    @test v[2,2,2] == [1, 2, 3]

    @test center(boundary) == [0, 0, 0]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule AzStorage_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("AzStorage")
JLLWrappers.@generate_main_file("AzStorage", UUID("00c928b4-b5f3-54d8-b38d-afd4635c4ad2"))
end  # module AzStorage_jll


text:
using Test
using Quon
using Quon: check

@testset "string-genus Rule" begin
    c = contract!(tait_copy(), tait_copy(), [2, 3], [1, 3])
    m = Quon.match(Quon.Rule{:string_genus}(), c)
    @test length(m) == 1
    @test check(c, m[1])
    g = m[1].vertices[1]
    rewrite!(c, m[1])
    @test !(g in c.genuses)
end

@testset "perm_rz" begin
    rz2 = contract!(tait_rz(0*im), tait_rz(0*im))
    m = Quon.match(Quon.Rule{:perm_rz}(), rz2)
    @test length(m) == 2
    @test check(rz2, m[1]) && check(rz2, m[2])

    rz = tait_rz(0*im)
    m = Quon.match(Quon.Rule{:perm_rz}(), rz)
    @test check(rz, m[1])
    Quon.rewrite!(rz, m[1])
    @test rz.g.f_max == 4
end

@testset "Yang-Baxter Rule" begin
    ry = contract!(contract!(tait_rz(pi/2*im), tait_rx(-π/2*im)), tait_rz(π/2*im))
    m1 = match(Rule{:yang_baxter_triangle}(), ry)
    @test length(m1) == 1
    @test check(ry, m1[1])
    rewrite!(ry, m1[1])

    m2 = match(Rule{:yang_baxter_star}(), ry)
    @test length(m2) == 1
    @test check(ry, m2[1])
    Quon.rewrite!(ry, m2[1])
    m3 = match(Rule{:yang_baxter_triangle}(), ry)
    @test length(m3) == 1
end

@testset "Fusion Rule" begin
    rx2 = contract!(tait_rx(π/2*im), tait_rx(π/2*im))
    mx = match(Rule{:x_fusion}(), rx2)
    @test check(rx2, mx[1])
    rewrite!(rx2, mx[1])
    @test length(rx2.quon_params) == 2

    rz2 = contract!(tait_rz(π/2*im), tait_rz(π/2*im))
    mz = match(Rule{:z_fusion}(), rz2)
    @test check(rz2, mz[1])
    rewrite!(rz2, mz[1])
    @test length(rz2.quon_params) == 2
end

@testset "charge remove" begin
    rx2 = contract!(tait_rx(pi*im), tait_rx(pi*im))
    m_charge_rm_v = match(Rule(:charge_rm_v), rx2)
    @test check(rx2, m_charge_rm_v[1])
    rewrite!(rx2, m_charge_rm_v[1])
    @test nv(rx2) == 5

    rz2 = contract!(tait_rz(pi*im), tait_rz(pi*im))
    m_charge_rm_f = match(Rule(:charge_rm_f), rz2)
    @test check(rz2, m_charge_rm_f[1])
    rewrite!(rz2, m_charge_rm_f[1])
    @test nv(rz2) == 5
end

@testset "SWAP-genus Rule" begin
    s2 = contract!(tait_swap(), tait_swap())
    m_swap_genus = match(Rule(:swap_genus), s2)
    @test check(s2, m_swap_genus[1])
    rewrite!(s2, m_swap_genus[1])
    @test !check(s2, m_swap_genus[1])
end

# TODO: fix locations for cz
# @testset "CZ^2 = I" begin
#     cz = tait_cz()
#     contract!(cz, tait_cz())
#     m = match(Rule{:perm_rz}(), cz)
#     Quon.rewrite!(cz, m[3])
#     Quon.rewrite!(cz, m[9])
#     m_string_genus = match(Quon.Rule{:string_genus}(), cz)
#     @test length(m_string_genus) == 1
#     Quon.rewrite!(cz, m_string_genus[1])
#     m_fusion = match(Quon.Rule{:z_fusion}(), cz)
#     @test length(m_fusion) == 3
#     Quon.rewrite!(cz, m_fusion[1])
#     Quon.rewrite!(cz, m_fusion[2])
#     Quon.rewrite!(cz, m_fusion[3])
#     m_yb = match(Quon.Rule{:yang_baxter_triangle}(), cz)
#     @test length(m_yb) == 1
#     rewrite!(cz, m_yb[1])
#     m_id = match(Quon.Rule{:identity}(), cz)
#     Quon.rewrite!(cz, m_id[1])
#     Quon.rewrite!(cz, m_id[3])
#     Quon.rewrite!(cz, m_id[5])
#     @test length(m_id) == 6
#     m_gf = match(Rule(:genus_fusion), cz)
#     @test length(m_gf) == 2
#     rewrite!(cz, m_gf[2])
# end


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

using Markdown
using InteractiveUtils

# ╔═╡ 31e93532-d915-4638-ac42-6d5a8b51472e
using PlutoUI

# ╔═╡ 9ba84189-80cd-4760-9e48-feca71f313c5
using DualNumbers

# ╔═╡ a37e3c58-bd22-4bc6-b003-bb8d350d49ed
using Unitful

# ╔═╡ b4b177dd-97db-4a6d-987c-64426b90bd51
bigbreak = HTML("<br>"^4);

# ╔═╡ 47e2cc32-8290-11ec-2253-e528a7b8198e
md"""
# Multiple dispatch ⭐
"""

# ╔═╡ ab103ac0-dba1-4f4f-91d8-3872f234436f
md"""### Starting with the familiar: Generic functions

If we write generic functions, Julia will figure out the types and try to make sense of it (if possible):
"""

# ╔═╡ 808b1e07-2912-4417-aef7-1b5de503e244
f(x) = x.^2

# ╔═╡ 849bd1c6-969b-479e-b519-f5d444bcb840
md"Hence, all of these inputs work:"

# ╔═╡ 67437820-a50b-49f2-af6b-7ba80277103e
(
	f_int    = f(1),
	f_array  = f([1,2]),
	f_string = f("hi")
)

# ╔═╡ e0fd7d97-187e-4937-9ca1-6473e167de4a
bigbreak

# ╔═╡ 19e8f83e-ca03-4677-9c6a-44f392d92a12
md"Indeed, under the hood, Julia does compile for all of these inputs specialised functions:"

# ╔═╡ 3e6b096d-58d3-46a6-8936-93186042c9b7
with_terminal() do
	@code_llvm f( (1.0, 1.0,1.0, 1.0) )
end

# ╔═╡ a213f381-52bf-4008-815e-b3c007b08b50
with_terminal() do
	@code_llvm f( [1.0, 1.0, 1.0] )
end

# ╔═╡ e6fc063f-19ef-4af4-b476-2c28ba161a93
md"""
### Short interlude

If you want to deep dive into how Julia works, it is often very interesting to check how Julia actually converts simple commands.
It needs some time to get used to the output, but it is worth it!

The `@code_*` commands are
```julia
@code_lowered f(x)   # Julia's internal representation after removing syntax sugar
@code_warntype f(x)  # highlighted lines mark type instabilities
@code_llvm f(x)      # the final product of the Julia JIT compiler
@code_native f(x)    # after optimisation with the LLVM for your machine
```

"""

# ╔═╡ 08714061-f90d-45f7-83ad-191b8d309e74

# ╔═╡ 2a255075-0bf6-4347-bf2a-1ded087d2514

# ╔═╡ 08d4c62e-1265-4925-8465-9ff1904057a3

# ╔═╡ 93051f8a-7802-4a2b-891b-46037a7f97b3
md"## Specifying the types

Sometimes the same function needs different implementation depending on the type:
"

# ╔═╡ 1c2b3808-ce65-4313-9380-58a42844eba9
begin
	# this function picks randomly between a and b:
	random(a, b) = rand([a,b])

	# integers between a and b
	random(a::Int64, b::Int64) = rand(a:b)

	# floats between a and b
	random(a::Float64, b::Float64) = a + (b-a)*rand()

	random() = 1
end

# ╔═╡ 1ff517d6-2145-4ab5-93e2-faeb0882dcdc
methods(random)

# ╔═╡ 4292041a-db6d-4235-809e-3bac423a7342
random(1.0, 2.0)

# ╔═╡ c0834dbb-0128-4219-a119-6df81bf876e4
random(1, 3)

# ╔═╡ a92b084e-ede0-408b-8bee-f8508fc2a470
random(1, 6.32)

# ╔═╡ 6682b551-f8da-4bc9-ad77-d41444ba6022

# ╔═╡ a9c019a5-a08e-467a-b635-fab0debd158d
md"#### Definition:

- **A function** (say `random`) is identified by the function name.
- The **different implementations** (say `random(::Int64,::Int64), ...`) of the function **are called methods**.
Hence, each function can have several methods.
"

# ╔═╡ 97036832-f4ca-4b34-b2da-6857439d529b

# ╔═╡ 88d8106c-2008-455b-a42a-fa3dd582c831
md"### → Julia always applies the most specific method for the given inputs!

These methods could overlap in crazy ways! _(It (usually) just works!™)_"

# ╔═╡ e53b709a-dd90-4a8c-bc59-781c685f4954

# ╔═╡ 2f381972-5972-4369-bb2d-67d01199eaf2

# ╔═╡ 2161f348-7efe-4cd0-b998-2d301fac22e4
md"## The key point of multiple dispatch

We want maximal control about how functions look for different types.

Multiple dispatch is
- **dynamic**
  - deciding at runtime what to do
- **modular**
  - reducing duplicate code
- **extensible**
  - allowing to add something to existing code
- **~ fast**
  - trying to keep amount of recomplilation minimal

(see: [Expression problem](https://en.wikipedia.org/wiki/Expression_problem))
"

# ╔═╡ 54eac40f-31aa-4df5-a1e6-3b863c12198f

# ╔═╡ 433fca14-845a-4267-b6be-b4ea1d9100c4
bigbreak

# ╔═╡ dd260866-8e10-4b36-a62f-2f930617e579
md"""
## The type hiearachy

How does Julia decide which method to call?

#### With a partial ordering `<:` on all types! 🧮 😎
"""

# ╔═╡ ba329815-e6bd-410f-82b4-bc7f89655661
Float64 <: Real <: Number

# ╔═╡ b7751fd9-6d65-4d5f-a570-4d71af818132
Int64 <: Integer <: Real <: Number

# ╔═╡ 72476e34-5ff3-48bf-848b-3e9be57faca3
Float64 <: Integer

# ╔═╡ 2431f98b-4e98-4528-acf6-65e27abb08c7
String <: AbstractString <: Any

# ╔═╡ f38fccdb-167c-44bd-8d9d-c3ae2fb79db7
Vector{Float64} <: Vector <: AbstractArray

# ╔═╡ ba7efcc8-479f-409b-846f-2696001d7a6b

# ╔═╡ 1ca91764-5f29-403e-8417-7d7cbb55042f
md"For each input with types `T1, T2, T3, ...`, Julia calls basically the method with input types

$\inf \{ \mathbf{S} \mid \text{method with argument types } \mathbf{S} \text{ exists and } \mathbf{T} <: \mathbf{S} \}$

with respect to the type ordering `<:`."

# ╔═╡ de4767b1-e882-4c55-9d22-0b91cc0d49d9
begin
	g(x::Number) = x
	g(x::Real) = x^2
	g(x::Float64) = x^3
end

# ╔═╡ 4d2f8eaa-fd3f-4b65-a804-7c1658a319c8
g( 1 + 2im )

# ╔═╡ 8d6623d2-60cf-43e7-b040-6787e8b7edad
g(3.0)

# ╔═╡ 352f54bf-3520-4f1c-9124-8d2871466e47
bigbreak

# ╔═╡ ce1a54c5-1ff0-4458-82f5-a1984b5be4fc
md"**Quick exercise:**

Find a case in which Julia cannot determine which method to call!"

# ╔═╡ 684ca6d9-ba49-446c-b394-3d4ca985e1d2
begin
	fnc(x::Int64,  y::Int64) = 1
	fnc(x::Number, y::Number) = 2
end

# ╔═╡ c76d5344-ebae-4bef-9d7c-09940be5c492
fnc(1,1.0)

# ╔═╡ b31fdae8-76f1-4998-945b-ca37f6173c0f

# ╔═╡ 82a2733a-da60-485f-8482-6189101933aa

# ╔═╡ 32a1a877-2146-4e8e-90c1-3b021b6d9348
md"**Note:** A case which is important to keep in mind is this one:"

# ╔═╡ 0a9e5987-869d-4b10-987d-4408769894c1
Vector{Int64} <: Vector

# ╔═╡ 9eb5aa13-7df1-4f2d-8f9b-d2b8b7e13a84
md"However, with so called _parametric types_ we can handle that:"

# ╔═╡ 57e75a99-1af9-452e-af11-81157551c12e
Vector{Int64} <: Vector{T} where T <: Integer

# ╔═╡ 61b1d73d-f7c3-4136-b92b-3a95b5143d2c
bigbreak

# ╔═╡ e97156d1-bea8-491f-b190-5c95d9a5b536
md"### How is that related to Object Oriented Programming?

> Yes and No: Multiple dispatch is more flexible 😃, but less ordered 😢.

In other words: OOP ensures that the programmer does not too stupid things. In Julia you can literally be a **type pirate** ☠.
"

# ╔═╡ c000aa68-b1ab-40ed-9376-32be9c51dc1a
md"
#### Short interlude: Let's be a type pirate ☠
_Run for example: `Base.:+(x::Float64, y::Float64) = 0.0` in a Julia REPL and compute `sin(1)`. (Afterwards, do something like that never again!)_
"

# ╔═╡ 4737b51a-539e-45ec-bf32-e3798d5ae5cd

# ╔═╡ c7a2a4bd-3e1d-438f-9270-725c3c9616f8
md" → You should usually only define methods for your custom types. Don't mess too much with other peoples types or it will backfire!"

# ╔═╡ 19412b07-a5ce-4c7a-b193-2b5e60c521f7
bigbreak

# ╔═╡ e4b689a9-f03a-4e28-8739-e018041d3c79
md"""
## Putting it all together: Being fast and generic is cool!

#### Example: (Derivative of solutions with respect to inital data)
"""

# ╔═╡ 3e537363-610e-40e0-8b8e-44e4713a17cd
md"Let's consider this simple implementation of a (scalar) explicit Euler method."

# ╔═╡ 4e965956-3633-4f4c-80ad-ff3f29958b54
function euler(f, x0::T, dt, n) where T
	x = Array{T}(undef, n)
	x[1] = x0

	for i in 2:n
		x[i] = x[i-1] + dt * f(x[i-1])
	end
	return x
end

# ╔═╡ ae2eedbd-6cb8-4356-924f-eed0926f388b
sol = euler( x -> -sin(x), 1.0, 0.1, 100)

# ╔═╡ 9228a5bc-84d5-4e5e-a971-211ac9b47fe3

# ╔═╡ 65c828b1-ff10-4428-afd9-ec1748054a5b
md"Now, let's see what happens if we plug in another type of `Numbers`!"

# ╔═╡ 71449f89-4c20-4eb6-b41b-5ac39f2ab829

# ╔═╡ c8639b54-fef0-485a-8e8d-10d90ab6ef95
md"### Automatic differentiation"

# ╔═╡ aba97bf0-2519-4ebb-8ead-d7abb4042adb
dual_sol = euler( x -> -sin(x), Dual(1.0, 1.0), 0.1, 100)

# ╔═╡ b3a96ab1-2285-4db3-b67e-040278da2312

# ╔═╡ 6890676a-37eb-4808-8967-5ac17ec286f4
md"##### It just works!™

We have implemented a method to compute the gradient of the numerical solution with respect to initial data. (aka Automatic Differentiation)"

# ╔═╡ 39072267-77e1-4b7e-9efc-f0a904ebf587
md"""
The solution $x(t, x_0)$ of

$\dot x = -\sin(x), \quad x(0) = x_0$

has approximately the derivative
$\frac{\partial x}{\partial x_0}(10, 1) ≈$ $(dualpart(dual_sol[end])).

_Notice that in order to do automatic differentiation, not a single line of the `euler` function had to change!_

"""

# ╔═╡ 6bf56af9-ad4f-417b-9356-9b582f9655e9

# ╔═╡ 91c8ea46-9fc4-4d94-94ec-0824555d4a2d
md"#### Units"

# ╔═╡ f5113d03-754f-4ace-91f3-27bf015cce9c
unit_sol = euler( x -> -x/u"hr", 1.0u"μm", 0.1u"s", 100)

# ╔═╡ 384f9598-abca-4727-8f26-037bc5f377bb

# ╔═╡ 23c8d96c-314d-40d8-a14e-a317378ad705

# ╔═╡ 59ca17c8-98bc-467b-9abc-495c978298d2

# ╔═╡ 5683bafe-96a6-404d-9ef3-4b493ac7e5e4

# ╔═╡ 7b92479f-2638-4bcc-822f-ce638a073d3c

# ╔═╡ 0a97a970-199e-4853-895e-fefc6d7753e6
md"### Unlimited possibilities to hijack this Euler method!"

# ╔═╡ fe5ad85b-f1d6-4d3e-a382-3d050c60c4cf
md"_Reminder:_ All of that compiled to almost C-like machine code!"

# ╔═╡ 8beec153-5a22-44b2-ae98-4cff82ec263c
bigbreak

# ╔═╡ e101a76d-d2ca-43c7-ab01-953c8397b7c9
md"""
# Exercises

7.1

A famous example for multiple dispatch is the
`Base.show(io, x)` function.

It is responsible for displaying Julia objects.
_(It works for text, numbers, images, html, ... you name it, then you can implement it)._

Overwrite the show function for `Tuple{Float64,Float64}` such that it displays
> 1 ❤ 2
when you call `print( (1,2) )`.

Hint: Use `print(io, ...)` to write something directly to input-output `io`.
"""

# ╔═╡ c1ed7bc2-2936-452c-b5b7-06de12adbe40
Base.show(io::IO, x::Tuple{Int64, Int64}) = nothing

# ╔═╡ 35f7e19f-bca1-4d5a-8f8d-dc4c4c77592e
with_terminal() do
	print( (1, 2+5) )
end

# ╔═╡ c79e90d9-d5e0-4f37-8abe-f49731a5d830
bigbreak

# ╔═╡ 716a277e-525b-4fed-81be-c8ceb50f9aec
md"""
7.2

Implement the product for `Tuple{Number,Number}` as
$(x,y) * (a,b) = xb - ya$

_Hint: The multiplication is called `Base.:*`_
"""

# ╔═╡ 8d48eb42-bbca-45c1-8c32-be9da3712835
Tup2f = Tuple{Number,Number} # if you want, you can use this type

# ╔═╡ 7e9c3310-e098-4e98-a02e-fdd905292f64
# enter your code here

# ╔═╡ 1fe523b2-f9b2-4a90-b45f-54dd9cb9dcd4

# ╔═╡ 214287e1-ac41-458a-aedf-16150f0b1da3
# (1.0,2) * (4,5) == -3.0  # test

# ╔═╡ 9b3fc16a-bee4-4917-92cf-41e64ab6ce39
bigbreak

# ╔═╡ e1befdce-9555-4a14-ad8e-48c64b4a74a6
md"""
7.3

Define a function `h` which returns the sum of integer arrays, but the product for all non-integer arrays.
"""

# ╔═╡ 838b3265-4724-4e56-95f9-56566862ffff
# enter your code here

# ╔═╡ 35e901dc-93f1-4927-a6c7-c01a6934ca83
# h([1, 2]) == 3
# h([1.0, 2.0]) == 2.0

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[compat]
DualNumbers = "~0.6.6"
PlutoUI = "~0.7.32"
Unitful = "~1.10.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f9982ef575e19b0e5c7a98c6e75ee496c0f73a93"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.12.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "84f04fe68a3176a583b864e492578b9466d87f1e"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.7"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "e6bf188613555c78062842777b116905a9f9dd49"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.0"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "b95e0b8a8d1b6a6c3e0b3ca393a7a285af47c264"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.10.1"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─31e93532-d915-4638-ac42-6d5a8b51472e
# ╟─b4b177dd-97db-4a6d-987c-64426b90bd51
# ╟─47e2cc32-8290-11ec-2253-e528a7b8198e
# ╟─ab103ac0-dba1-4f4f-91d8-3872f234436f
# ╠═808b1e07-2912-4417-aef7-1b5de503e244
# ╟─849bd1c6-969b-479e-b519-f5d444bcb840
# ╠═67437820-a50b-49f2-af6b-7ba80277103e
# ╟─e0fd7d97-187e-4937-9ca1-6473e167de4a
# ╟─19e8f83e-ca03-4677-9c6a-44f392d92a12
# ╠═3e6b096d-58d3-46a6-8936-93186042c9b7
# ╠═a213f381-52bf-4008-815e-b3c007b08b50
# ╟─e6fc063f-19ef-4af4-b476-2c28ba161a93
# ╟─08714061-f90d-45f7-83ad-191b8d309e74
# ╟─2a255075-0bf6-4347-bf2a-1ded087d2514
# ╟─08d4c62e-1265-4925-8465-9ff1904057a3
# ╟─93051f8a-7802-4a2b-891b-46037a7f97b3
# ╠═1c2b3808-ce65-4313-9380-58a42844eba9
# ╟─1ff517d6-2145-4ab5-93e2-faeb0882dcdc
# ╠═4292041a-db6d-4235-809e-3bac423a7342
# ╠═c0834dbb-0128-4219-a119-6df81bf876e4
# ╠═a92b084e-ede0-408b-8bee-f8508fc2a470
# ╟─6682b551-f8da-4bc9-ad77-d41444ba6022
# ╟─a9c019a5-a08e-467a-b635-fab0debd158d
# ╟─97036832-f4ca-4b34-b2da-6857439d529b
# ╟─88d8106c-2008-455b-a42a-fa3dd582c831
# ╟─e53b709a-dd90-4a8c-bc59-781c685f4954
# ╟─2f381972-5972-4369-bb2d-67d01199eaf2
# ╟─2161f348-7efe-4cd0-b998-2d301fac22e4
# ╟─54eac40f-31aa-4df5-a1e6-3b863c12198f
# ╟─433fca14-845a-4267-b6be-b4ea1d9100c4
# ╟─dd260866-8e10-4b36-a62f-2f930617e579
# ╠═ba329815-e6bd-410f-82b4-bc7f89655661
# ╠═b7751fd9-6d65-4d5f-a570-4d71af818132
# ╠═72476e34-5ff3-48bf-848b-3e9be57faca3
# ╠═2431f98b-4e98-4528-acf6-65e27abb08c7
# ╠═f38fccdb-167c-44bd-8d9d-c3ae2fb79db7
# ╟─ba7efcc8-479f-409b-846f-2696001d7a6b
# ╟─1ca91764-5f29-403e-8417-7d7cbb55042f
# ╠═de4767b1-e882-4c55-9d22-0b91cc0d49d9
# ╠═4d2f8eaa-fd3f-4b65-a804-7c1658a319c8
# ╠═8d6623d2-60cf-43e7-b040-6787e8b7edad
# ╟─352f54bf-3520-4f1c-9124-8d2871466e47
# ╟─ce1a54c5-1ff0-4458-82f5-a1984b5be4fc
# ╠═684ca6d9-ba49-446c-b394-3d4ca985e1d2
# ╠═c76d5344-ebae-4bef-9d7c-09940be5c492
# ╟─b31fdae8-76f1-4998-945b-ca37f6173c0f
# ╟─82a2733a-da60-485f-8482-6189101933aa
# ╟─32a1a877-2146-4e8e-90c1-3b021b6d9348
# ╠═0a9e5987-869d-4b10-987d-4408769894c1
# ╟─9eb5aa13-7df1-4f2d-8f9b-d2b8b7e13a84
# ╠═57e75a99-1af9-452e-af11-81157551c12e
# ╟─61b1d73d-f7c3-4136-b92b-3a95b5143d2c
# ╟─e97156d1-bea8-491f-b190-5c95d9a5b536
# ╟─c000aa68-b1ab-40ed-9376-32be9c51dc1a
# ╟─4737b51a-539e-45ec-bf32-e3798d5ae5cd
# ╟─c7a2a4bd-3e1d-438f-9270-725c3c9616f8
# ╟─19412b07-a5ce-4c7a-b193-2b5e60c521f7
# ╟─e4b689a9-f03a-4e28-8739-e018041d3c79
# ╟─3e537363-610e-40e0-8b8e-44e4713a17cd
# ╠═4e965956-3633-4f4c-80ad-ff3f29958b54
# ╠═ae2eedbd-6cb8-4356-924f-eed0926f388b
# ╟─9228a5bc-84d5-4e5e-a971-211ac9b47fe3
# ╟─65c828b1-ff10-4428-afd9-ec1748054a5b
# ╟─71449f89-4c20-4eb6-b41b-5ac39f2ab829
# ╟─c8639b54-fef0-485a-8e8d-10d90ab6ef95
# ╠═9ba84189-80cd-4760-9e48-feca71f313c5
# ╠═aba97bf0-2519-4ebb-8ead-d7abb4042adb
# ╟─b3a96ab1-2285-4db3-b67e-040278da2312
# ╟─6890676a-37eb-4808-8967-5ac17ec286f4
# ╟─39072267-77e1-4b7e-9efc-f0a904ebf587
# ╟─6bf56af9-ad4f-417b-9356-9b582f9655e9
# ╟─91c8ea46-9fc4-4d94-94ec-0824555d4a2d
# ╠═a37e3c58-bd22-4bc6-b003-bb8d350d49ed
# ╠═f5113d03-754f-4ace-91f3-27bf015cce9c
# ╟─384f9598-abca-4727-8f26-037bc5f377bb
# ╠═23c8d96c-314d-40d8-a14e-a317378ad705
# ╠═59ca17c8-98bc-467b-9abc-495c978298d2
# ╟─5683bafe-96a6-404d-9ef3-4b493ac7e5e4
# ╟─7b92479f-2638-4bcc-822f-ce638a073d3c
# ╟─0a97a970-199e-4853-895e-fefc6d7753e6
# ╟─fe5ad85b-f1d6-4d3e-a382-3d050c60c4cf
# ╟─8beec153-5a22-44b2-ae98-4cff82ec263c
# ╟─e101a76d-d2ca-43c7-ab01-953c8397b7c9
# ╠═c1ed7bc2-2936-452c-b5b7-06de12adbe40
# ╠═35f7e19f-bca1-4d5a-8f8d-dc4c4c77592e
# ╟─c79e90d9-d5e0-4f37-8abe-f49731a5d830
# ╟─716a277e-525b-4fed-81be-c8ceb50f9aec
# ╠═8d48eb42-bbca-45c1-8c32-be9da3712835
# ╠═7e9c3310-e098-4e98-a02e-fdd905292f64
# ╟─1fe523b2-f9b2-4a90-b45f-54dd9cb9dcd4
# ╠═214287e1-ac41-458a-aedf-16150f0b1da3
# ╟─9b3fc16a-bee4-4917-92cf-41e64ab6ce39
# ╟─e1befdce-9555-4a14-ad8e-48c64b4a74a6
# ╠═838b3265-4724-4e56-95f9-56566862ffff
# ╠═35e901dc-93f1-4927-a6c7-c01a6934ca83
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/7 - Multiple dispatch.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/7 - Multiple dispatch.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/7 - Multiple dispatch.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
###### beam splitter and special counter example matrices ######

function beam_splitter(transmission_amplitude = sqrt(0.5))

	"""2d beam_splitter matrix, follows the conventions of Leonardo"""

	# |t|^2 is the "transmission probability"

	t = transmission_amplitude
	r = sqrt(1-t^2)

	bs = [[t -r]; [r t]]

end

function beam_splitter_modes(;in_up, in_down, out_up, out_down, transmission_amplitude, n)

	"""beam splitter incorrporated for connecting specific input modes, output modes in a dimension n interferometer"""

	if in_up == in_down || out_up == out_down
		throw(ArgumentError())
	end

	bs = Matrix{ComplexF64}(I, n, n)

	sub_bs = beam_splitter(transmission_amplitude)

	bs[in_up, out_up] = sub_bs[1,1]
	bs[in_down, out_down] = sub_bs[2,2]
	bs[in_down, out_up] = sub_bs[2,1]
	bs[in_up, out_down] = sub_bs[1,2]
	#
	# # before convention checks:
	# bs[in_down, out_up] = sub_bs[1,2] ####### this change is a bit ad hoc and needs to be checked carefully with the conventions
	# bs[in_up, out_down] = sub_bs[2,1]
	#

	bs

end

function rotation_matrix(angle::Float64)
	[cos(angle) -sin(angle);
	 sin(angle) coss(angle)]
 end

 function rotation_matrix_modes(;in_up, in_dow, out_up, out_down, angle, n)

	 """beam splitter incorrporated for connecting specific input modes, output modes in a dimension n interferometer"""

 	if in_up == in_down || out_up == out_down
 		throw(ArgumentError())
 	end

 	r = Matrix{ComplexF64}(I, n, n)

 	sub_r = rotation_matrix(angle)

 	r[in_up, out_up] = sub_r[1,1]
 	r[in_down, out_down] = sub_r[2,2]
 	# bs[in_down, out_up] = sub_bs[2,1]
 	# bs[in_up, out_down] = sub_bs[1,2]
 	r[in_down, out_up] = sub_r[1,2] ####### this change is a bit ad hoc and needs to be checked carefully with the conventions
 	r[in_up, out_down] = sub_r[2,1]

 	r

 end

 function phase_shift(shifted_modes::Array, param_::Array)

	U = zeros(length(shifted_modes), length(shifted_modes))

 	for i in length(shifted_modes)
		for j in length(shifted_modes)
			i == j && i in shifted_modes ? U[i,j] = exp(1im * param_[i]) : continue
		end
	end

	  U

end

# function build_circuit(mode_occ::ModeOccupation, circuit_elements::Vector{Vector{T}}, layers::Vector) where T <: Interferometer
#
# 	m = mode_occ.m
# 	id = Matrix{ComplexF64}(LinearAlgebra.I, m, m)
# 	res = id
#
# 	for i in 1:length(circuit_elements)
# 		u = id
# 		for j in 1:length(circuit_elements[i])
# 			U = circuit_elements[i][j].U
# 			u[layers[i][j], layers[i][j]] = U
# 		end
#
# 		res *= u
# 	end
#
# 	res
# end
#
# mode_occ = first_modes(3,5)
#
# m = mode_occ.m
#
# b1 = BeamSplitter(1/sqrt(6))
#
# is_unitary(b1.U)
#
# U1 = RandHaar(1)
#
# is_unitary(U1.U)
#
# b2 = BeamSplitter(1/sqrt(2))
#
# is_unitary(b2.U)
#
# U2 = RandHaar(1)
#
# is_unitary(U2.U)
#
# circuit_elem = [[b1,U1],[b2,U2]]
#
# v_ = [[[1,2],[3]], [[1,3],[2]]]
#
# res = build_circuit(mode_occ, circuit_elem, v_)
#
# is_unitary(res)
#
# res * conj(transpose(res))
#


text:
# ---
# title: 464. Can I Win
# id: problem464
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Minimax
# link: <https://leetcode.com/problems/can-i-win/description/>
# hidden: true
# ---
#
# In the "100 game" two players take turns adding, to a running total, any
# integer from `1` to `10`. The player who first causes the running total to
# **reach or exceed** 100 wins.
#
# What if we change the game so that players **cannot** re-use integers?
#
# For example, two players might take turns drawing from a common pool of
# numbers from 1 to 15 without replacement until they reach a total >= 100.
#
# Given two integers maxChoosableInteger and desiredTotal, return `true` if the
# first player to move can force a win, otherwise return `false`. Assume both
# players play **optimally**.
#
#
#
# **Example 1:**
#
#
#
#     Input: maxChoosableInteger = 10, desiredTotal = 11
#     Output: false
#     Explanation:
#     No matter which integer the first player choose, the first player will lose.
#     The first player can choose an integer from 1 up to 10.
#     If the first player choose 1, the second player can only choose integers from 2 up to 10.
#     The second player will win by choosing 10 and get a total = 11, which is >= desiredTotal.
#     Same with other integers chosen by the first player, the second player will always win.
#
#
# **Example 2:**
#
#
#
#     Input: maxChoosableInteger = 10, desiredTotal = 0
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: maxChoosableInteger = 10, desiredTotal = 1
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= maxChoosableInteger <= 20`
#   * `0 <= desiredTotal <= 300`
#
#
## @lc code=start
using LeetCode

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


text:
module TestAqua

import Aqua
import ThreadedScans

test() = Aqua.test_all(ThreadedScans)

end  # module


text:
using PlotlyLight
using JSON3
using Test

#-----------------------------------------------------------------------------# simple checks
@testset "simple checks" begin
    @test Plot(Config(x = 1:10)) isa Plot
    @test Plot(Config(x = 1:10), Config(title="Title")) isa Plot
    @test Plot(Config(x = 1:10), Config(title="Title"), Config(displaylogo=true)) isa Plot
end
#-----------------------------------------------------------------------------# src
@testset "src" begin
    p = Plot(Config(y=1:10))

    PlotlyLight.src!(:cdn)
    @test occursin("cdn", repr("text/html", p))

    PlotlyLight.src!(:none)
    @test !occursin("cdn", repr("text/html", p))

    PlotlyLight.src!(:standalone)
    @test length(repr("text/html", p)) > 1000

    PlotlyLight.src!(:local)
    @test occursin("artifacts", repr("text/html", p))
end
#-----------------------------------------------------------------------------# templates
@testset "templates" begin
    PlotlyLight.template!(:plotly_dark)
    p = Plot(Config(y=1:10))
    s = repr("text/html", p)

    PlotlyLight.template!(:plotly_white)
    p2 = Plot(Config(y=1:10))
    s2 = repr("text/html", p2)

    @test s != s2
end
#-----------------------------------------------------------------------------# vecvec
@testset "collectrows" begin
    x = rand(3, 2)
    v = collectrows(x)
    @test v[1] == x[1, :]
    @test v[2] == x[2, :]
    @test v[3] == x[3, :]
end


text:
println("\n Testing solver build \n")

x = Variable(10)
A = randn(5, 10)
y = Variable(7)
B = randn(5, 7)
b = randn(5)

prob = problem(ls(A*x + b), norm(x, 2) <= 1.0)
built_slv = build(prob, StructuredOptimization.PG())
solve!(built_slv)

~x .= 0.
prob = problem(ls(A*x - B*y + b) + norm(y, 1), norm(x, 2) <= 1.0)
built_slv = build(prob, FPG())
solve!(built_slv)

println("\n Testing @minimize \n")
~x .= 0.
~y .= 0.
slv, = @minimize ls(A*x - B*y + b) st norm(x, 2) <= 1e4, norm(y, 1) <= 1.0 with PG()
~x .= 0.
slv, = @minimize ls(A*x - b) st norm(x, 1) <= 1.0 with PG()
~x .= 0.
slv, = @minimize ls(A*x - b) st norm(x, 1) <= 1.0
~x .= 0.
slv, = @minimize ls(A*x - b) + norm(x, 1) with PG()
~x .= 0.
slv, = @minimize ls(A*x - b) + norm(x, 1)
~x .= 0.
slv, = @minimize ls(A*x - b)

#TODO many many more tests
x = Variable(5)
A = randn(10, 5)
b = randn(10)

println("\nTesting @minimize nonlinear \n")
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with PG()
xpg = copy(~x)
~x .= 0.
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with ZeroFPR()
xz = copy(~x)
~x .= 0.
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with PANOC()
xp = copy(~x)
~x .= 0.

@test norm(xz-xpg) <1e-4
@test norm(xp-xpg) <1e-4

# test nonconvex Rosenbrock function with known minimum
solvers = ["ZeroFPR(tol = 1e-6)","PANOC(tol = 1e-6)"]
for slv in solvers
    solver = eval(Meta.parse(slv))
    x = Variable(1)
    y = Variable(1)
    a,b = 2.0, 100.0

    cf = norm(x-a)^2+b*norm(pow(x,2)-y)^2
    @minimize cf+1e-10*norm(x,1)+1e-10*norm(y,1) with solver

    @test norm(~x-[a]) < 1e-4
    @test norm(~y-[a^2]) < 1e-4
end


text:
#########################################################
# These are helper functions for the DisjointSets type. #
#########################################################

# This takes a DisjointSets object and returns its ground set.
function _ground_set{T}(DS::DisjointSets{T})
    G = Set{T}()
    for item in keys(DS.intmap)
        push!(G,item)
    end
    return G
end

# Create a set of sets containing the parts of this partition.
function set_of_sets{T}(DS::DisjointSets{T})
    n = num_groups(DS)
    GS = _ground_set(DS)

    # Get root elements
    roots = Set{Int}()
    for item in GS
        r = find_root(DS,item)
        push!(roots,r)
    end
    roots = collect(roots)

    # Map root numbers to [1:n]
    rootmap = Dict{Int,Int}()
    for k=1:n
        rootmap[roots[k]] = k
    end

    # Create an array of sets to hold the parts
    parts = Array(Set{T},n)
    for k=1:n
        parts[k] = Set{T}()
    end

    # place each ground set item in its appropriate part
    for item in GS
        index = rootmap[find_root(DS,item)]
        push!(parts[index], item)
    end

    # now take those parts and pack them into a set
    P = Set{Set{T}}() # This is a set of sets
    for item in parts
        push!(P,item)
    end

    return P
end


text:
module baffled_piston_examples
using FinEtools
using FinEtoolsAcoustics
using LinearAlgebra

function baffled_piston_H8_ABC_example_algo()
    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  7500*phun("rev/s");      # frequency of the piston
    a_piston =  -1.0*phun("mm/s")     # amplitude of the piston acceleration
    R = 50.0*phun("mm");# radius of the piston
    Ro = 150.0*phun("mm"); # radius of the enclosure
    nref = 4;#number of refinements of the sphere around the piston
    nlayers = 35;                     # number of layers of elements surrounding the piston
    tolerance = R/(2^nref)/100

    println("""

    Baffled piston in a half-sphere domain with ABC.

    Hexahedral mesh. Algorithm version.
    """)

    t0  =  time()

    # Hexahedral mesh
    fens,fes  =  H8sphere(R,nref);
    bfes  =  meshboundary(fes)
    File  =   "baffledabc_boundary.vtk"
    vtkexportmesh(File, connasarray(bfes), fens.xyz, FinEtools.MeshExportModule.Q4)
    @async run(`"paraview.exe" $File`)

    l = selectelem(fens,bfes,facing = true,direction = [1.0 1.0  1.0], dotmin= 0.001)
    ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
    fens1,fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)

    # Piston surface mesh
    bfes  =  meshboundary(fes)
    l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
    l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
    piston_fes = subset(bfes,intersect(l1,l2));

    # Outer spherical boundary
    louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
    outer_fes = subset(bfes,louter);

    println("Pre-processing time elapsed  =  ",time() - t0,"s")

    t1  =  time()

    material = MatAcoustFluid(bulk, rho)
    # Region of the fluid
    region1 =  FDataDict("femm"=>FEMMAcoust(IntegDomain(fes, GaussRule(3, 2)), material))

    # Surface for the ABC
    abc1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(outer_fes, GaussRule(2, 2)),
    material))

    # Surface of the piston
    flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(piston_fes, GaussRule(2, 2)),
    material),  "normal_flux"=> -rho*a_piston+0.0im);

    # Make model data
    modeldata =  FDataDict("fens"=>  fens,
    "omega"=>omega,
    "regions"=>[region1],
    "flux_bcs"=>[flux1], "ABCs"=>[abc1])

    # Call the solver
    modeldata = FinEtoolsAcoustics.AlgoAcoustModule.steadystate(modeldata)

    println("Computing time elapsed  =  ",time() - t1,"s")
    println("Total time elapsed  =  ",time() - t0,"s")

    geom = modeldata["geom"]
    P = modeldata["P"]

    File  =   "baffledabc.vtk"
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H8;
    scalars = [("absP", abs.(P.values))])
    @async run(`"paraview.exe" $File`)
end # baffled_piston_H8_ABC_example_algo

function allrun()
    println("#####################################################")
    println("# baffled_piston_H8_ABC_example_algo ")
    baffled_piston_H8_ABC_example_algo()
    return true
end # function allrun

end # module baffled_piston_examples


text:
module Soft
using Spec
using ..Omega
using ..Omega: TaggedΩ, tag
import ..Space: combinetags
using ..NonDet: apl
import Omega
import ForwardDiff
import Cassette
using LinearAlgebra: norm

using DocStringExtensions
export  d,
        SoftBool,
        softeq,
        softlt,
        softgt,
        >ₛ,
        >=ₛ,
        <=ₛ,
        <ₛ,
        ==ₛ,
        err,
        logerr,
        anyₛ,
        allₛ,

        # Kernels
        kse,
        kseα,
        kf1,
        kf1β,
        withkernel,
        atα,
        @atα,

        indomain,
        indomainₛ,
        applynotrackerr,
        applytrackerr

abstract type AbstractSoftBool <: Real end
const AbstractBool = Union{AbstractSoftBool, Bool}

include("inf.jl")            # Infinity
include("kernels.jl")        # Kernels
include("kernelctx.jl")      # Kernels context
include("softbool.jl")       # Soft Boolean
include("dualsoftbool.jl")   # Dual Soft Boolean
include("distances.jl")      # Standard Distance Functions
include("trackerror.jl")     # Tracking error
include("any.jl")            # any, all

# Using Dual Soft Bools
softeq(x, y) = dsofteq(x, y)
softgt(x, y) = dsoftgt(x, y)
softlt(x, y) = dsoftlt(x, y)

const softtrue = dsofttrue
const softfalse = dsoftfalse

# Using Not Dual Soft Bools
# const softeq = ssofteq
# const softgt = ssoftgt
# const softlt = ssoftlt

# const softtrue = ssofttrue
# const softfalse = ssoftfalse

const >ₛ = softgt
const >=ₛ = softgt
const <=ₛ = softlt
const <ₛ = softlt
const ==ₛ = softeq

# Lifts #
  # Omega.lift(:softeq, 2)
  # Omega.lift(:softeq, 3)
  # Omega.lift(:softgt, 2)
  # Omega.lift(:softlt, 2)
  # Omega.lift(:logerr, 1)
  # Omega.lift(:err, 1)
  # Omega.lift(:kf1β, 1)
  # Omega.lift(:kseα, 1)
  # Omega.lift(:logkseα, 1)

end


text:
module TestList

using DocumentationOverview
using Test

function test_default()
    list = DocumentationOverview.list(DocumentationOverview)
    text = sprint(show, "text/plain", list)
    @test occursin("•  ", text)
    @test occursin("DocumentationOverview.table", text)
    @test occursin("DocumentationOverview.find", text)
    @test occursin("DocumentationOverview.API", text)
end

end  # module


text:
struct LwplsQda
    X::Array{Float64}
    y::AbstractMatrix
    fm
    metric::String
    h::Real
    k::Int
    nlv::Int
    prior::String
    tol::Real
    verbose::Bool
    lev::AbstractVector
    ni::AbstractVector
end

"""
    lwplsqda(X, y; nlvdis, metric, h, k, nlv, prior = "unif",
        tol = 1e-4, verbose = false)
kNN-LWPLS-QDA models.
* `X` : X-data.
* `y` : y-data (class membership).
* `nlvdis` : Number of latent variables (LVs) to consider in the
    global PLS used for the dimension reduction before
    calculating the dissimilarities. If `nlvdis = 0`, there is no dimension reduction.
* `metric` : Type of dissimilarity used to select the neighbors.
    Possible values are "eucl" (default; Euclidean distance)
    and "mahal" (Mahalanobis distance).
* `h` : A scalar defining the shape of the weight function. Lower is h,
    sharper is the function. See function `wdist`.
* `k` : The number of nearest neighbors to select for each observation to predict.
* `nlv` : Nb. latent variables (LVs).
* `prior` : Type of prior probabilities for class membership
    (`unif`: uniform; `prop`: proportional).
* `tol` : For stabilization when very close neighbors.
* `verbose` : If true, fitting information are printed.

This is the same methodology as for `lwplsr` except that
PLSR is replaced by PLS-QDA.
"""
function lwplsqda(X, y; nlvdis, metric, h, k, nlv,
    prior = "unif", tol = 1e-4, verbose = false)
    X = ensure_mat(X)
    y = ensure_mat(y)
    ztab = tab(y)
    if nlvdis == 0
        fm = nothing
    else
        fm = plskern(X, dummy(y).Y; nlv = nlvdis)
    end
    return LwplsQda(X, y, fm, metric, h, k, nlv, prior, tol, verbose,
        ztab.keys, ztab.vals)
end

"""
    predict(object::LwplsQda, X)
Compute the y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::LwplsQda, X; nlv = nothing)
    X = ensure_mat(X)
    m = size(X, 1)
    a = object.nlv
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    # Getknn
    if isnothing(object.fm)
        res = getknn(object.X, X; k = object.k, metric = object.metric)
    else
        res = getknn(object.fm.T, transform(object.fm, X);
            k = object.k, metric = object.metric)
    end
    listw = copy(res.d)
    for i = 1:m
        w = wdist(res.d[i]; h = object.h)
        w[w .< object.tol] .= object.tol
        listw[i] = w
    end
    # End
    pred = locwlv(object.X, object.y, X;
        listnn = res.ind, listw = listw, fun = plsqda, nlv = nlv,
        prior = object.prior, verbose = object.verbose).pred
    (pred = pred, listnn = res.ind, listd = res.d, listw = listw)
end


text:
using WordCloud
using Random
words = [randstring(rand(1:8)) for i in 1:300]
weights = randexp(length(words)) .* 1000 .+ rand(1:100, length(words))
wordcloud(words, weights) |> generate!


text:
using PiecewiseLinear
using Base.Test

using JuMP

let d = linspace(1,2π,8), f = sin
    mCC = Model()
    @variable(mCC, xCC)
    zCC = piecewiselinear(mCC, xCC, d, sin, method=:CC)
    @objective(mCC, Max, zCC)
    @test solve(mCC) == :Optimal
    @test isapprox(getvalue(xCC), 1.75474, rtol=1e-4)
    @test isapprox(getvalue(zCC), 0.98313, rtol=1e-4)

    @constraint(mCC, xCC ≤ 1.5zCC)
    @test solve(mCC) == :Optimal
    @test isapprox(getvalue(xCC), 1.36495, rtol=1e-4)
    @test isapprox(getvalue(zCC), 0.90997, rtol=1e-4)
end

let d = linspace(1,2π,8), f = sin
    mMC = Model()
    @variable(mMC, xMC)
    zMC = piecewiselinear(mMC, xMC, d, sin, method=:MC)
    @objective(mMC, Max, zMC)
    @test solve(mMC) == :Optimal
    @test isapprox(getvalue(xMC), 1.75474, rtol=1e-4)
    @test isapprox(getvalue(zMC), 0.98313, rtol=1e-4)

    @constraint(mMC, xMC ≤ 1.5zMC)
    @test solve(mMC) == :Optimal
    @test isapprox(getvalue(xMC), 1.36495, rtol=1e-4)
    @test isapprox(getvalue(zMC), 0.90997, rtol=1e-4)
end

let d = linspace(1,2π,8), f = sin
    mInc = Model()
    @variable(mInc, xInc)
    zInc = piecewiselinear(mInc, xInc, d, sin, method=:Incremental)
    @objective(mInc, Max, zInc)
    @test solve(mInc) == :Optimal
    @test isapprox(getvalue(xInc), 1.75474, rtol=1e-4)
    @test isapprox(getvalue(zInc), 0.98313, rtol=1e-4)

    @constraint(mInc, xInc ≤ 1.5zInc)
    @test solve(mInc) == :Optimal
    @test isapprox(getvalue(xInc), 1.36495, rtol=1e-4)
    @test isapprox(getvalue(zInc), 0.90997, rtol=1e-4)
end

let d = linspace(1,2π,8), f = sin
    mLog = Model()
    @variable(mLog, xLog)
    zLog = piecewiselinear(mLog, xLog, d, sin, method=:Logarithmic)
    @objective(mLog, Max, zLog)
    @test solve(mLog) == :Optimal
    @test isapprox(getvalue(xLog), 1.75474, rtol=1e-4)
    @test isapprox(getvalue(zLog), 0.98313, rtol=1e-4)

    @constraint(mLog, xLog ≤ 1.5zLog)
    @test solve(mLog) == :Optimal
    @test isapprox(getvalue(xLog), 1.36495, rtol=1e-4)
    @test isapprox(getvalue(zLog), 0.90997, rtol=1e-4)
end


text:
function compress_then_write(b::Vector{SubString{String}}, io)
    @warn "JDF is write a SubString vector. When loaded back it will be a String vector not a SubString vector"
    compress_then_write(String.(b), io)
end

function column_loader(::Type{SubString{String}}, io, metadata)
    @warn "JDF is loading SubString vector. It will be loaded as a String vector not a SubString vector"
    column_loader(String, io, metadata)
end