text:
mutable struct RBMBatchedCache{VT,VS,VST} <: NNBatchedCache{RBMSplit}
    θ::VT
    θ_tmp::VT
    logℒθ::VT
    ∂logℒθ::VT

    # complex sigmas
    res::VS #batch
    res_tmp::VST #batch

    # states
    σ::VT

    valid::Bool # = false
end

cache(net::RBM, batch_sz) = begin
    n_v = length(net.a)
    n_h = length(net.b)
    return RBMBatchedCache(
                  similar(net.b, n_h, batch_sz),
                  similar(net.b, n_h, batch_sz),
                  similar(net.b, n_h, batch_sz),
                  similar(net.b, n_h, batch_sz),
                  similar(net.b, 1, batch_sz),
                  similar(net.b, 1, batch_sz),
                  similar(net.b, n_v, batch_sz),
                  false)
end
batch_size(c::RBMBatchedCache) = size(c.θ, 2)

function Base.show(io::IO, m::RBMBatchedCache)
    print(io, "RBMBatchedCache with batch-size = $(batch_size(m))")
end

function logψ!(out::AbstractArray, net::RBM, c::RBMBatchedCache, σ_r::AStateBatch)
    θ = c.θ
    θ_tmp = c.θ_tmp
    logℒθ = c.logℒθ
    T = eltype(θ)

    # copy the states to complex valued states for the computations.
    σ = c.σ;  σ.=σ_r #σ = copy!(c.σ, σ_r)

    #θ .= net.b .+ net.W * σ
    mul!(θ, net.W, σ)
    θ .+= net.b
    logℒθ .= net.f.(θ)

    #res = σᵗ*net.a + sum(logℒθ, dims=1)
    mul!(out, transpose(net.a), σ)
    Base.mapreducedim!(identity, +, out, logℒθ)

    return out
end

function logψ_and_∇logψ!(∇logψ, out, net::RBM, c::RBMBatchedCache, σ_r::AStateBatch)
    θ = c.θ
    θ_tmp = c.θ_tmp
    logℒθ = c.logℒθ
    ∂logℒθ = c.∂logℒθ
    res = out # c.res
    res_tmp = c.res_tmp
    T = eltype(θ)

    # copy the states to complex valued states for the computations.
    σ = c.σ;  σ.=σ_r #σ = copy!(c.σ, σ_r)

    #θ .= net.b .+ net.W * σ
    mul!(θ, net.W, σ)
    θ .+= net.b
    logℒθ  .= net.f.(θ)
    ∂logℒθ .= fwd_der.(net.f, θ)

    #res = σ'*net.a + sum(logℒθ, dims=1)
    mul!(out, transpose(net.a), σ)
    #sum!(logℒθ, res, init=false)
    Base.mapreducedim!(identity, +, res, logℒθ)

    ∇logψ.a   .= σ
    ∇logψ.b   .= ∂logℒθ
    #∇logψ.W  .= ∂logℒθ .* transpose(σ)

    _batched_outer_prod!(∇logψ.W, ∂logℒθ, σ)

    # TODO make this better
    #copyto!(out, 1, res, 1, length(out))

    return out
end


text:
function get_maxima_minima(amplitudes)
    # Find all local maxima
    maxima = zeros(size(amplitudes))
    begin
        maxima_idx = findall(diff(sign.(diff(amplitudes))) .== -2) .+ 1
        maxima[maxima_idx] .= 1.0
        if amplitudes[1] > amplitudes[2]
            maxima[1] = 1.0
        end
        if amplitudes[end] > amplitudes[end-1]
            maxima[end] = 1.0
        end
    end
    maxima .*= amplitudes

    minima = ones(size(amplitudes)) .* Inf
    begin
        minima_idx = findall(diff(sign.(diff(amplitudes))) .== 2) .+ 1
        if amplitudes[1] < amplitudes[2]
            minima[1] = amplitudes[1]
        end
        if amplitudes[end] < amplitudes[end-1]
            maxima[end] = amplitudes[end]
        end
        minima[minima_idx] .= amplitudes[minima_idx]
    end

    return maxima, minima
end

function get_regions(amplitudes, minima, pitch, bin_size)

    region_size = 0.5 * pitch / bin_size + 1

    regions = []
    for k in 1:200
        idx = pitch * k / bin_size + 1
        idx_int = Int(round(idx))
        if idx_int >= length(amplitudes)
            break
        end

        region_start = max(1, Int(round(idx - region_size)))
        region_min = findmin(minima[region_start:idx_int-1])
        if region_min[1] == Inf
            region_min = findmin(amplitudes[region_start:idx_int-1])
        end
        region_start = region_min[2] + region_start - 1

        region_end = min(Int(round(idx + region_size)), length(amplitudes))
        region_min = findmin(minima[idx_int+1:region_end])
        if region_min[1] == Inf
            region_min = findmin(amplitudes[idx_int+1:region_end])
        end
        region_end = region_min[2] + idx_int - 1

        push!(regions, (idx, region_start:region_end))
    end

    idx = 1
    noise = []
    for i in 1:length(regions)
        push!(noise, idx:(regions[i][2][1]-1))
        idx = regions[i][2][end] + 1
    end
    push!(noise, (regions[end][2][end]+1):length(amplitudes))

    return regions, noise
end

function interpolate_FFT(k, s)
    N = length(s)
    sum([s[n] * exp(-im * 2 * π * k * n / N) for n in 1:N])
end

function linear_interpolation(s, k)
    x_1 = Int.(ceil.(k))
    x_0 = Int.(floor.(k))
    return (s[x_0] .* (x_1 .- k) .+ s[x_1] .* (k .- x_0)) / (x_1[1] .- x_0[1])
end

##################################################################################
using GLMakie, Makie.Colors
set_theme!(theme_black())
using WAV
using DSP
include("Utils.jl")

function get_song(file)
    y, Fs = wavread(string("Inputs/", file, ".wav"))
    low_note = Utils.getNote(1) #hz

    # Signal pre-processing
    y = [(y[i, 1] + y[i, 2]) / 2 for i = 1:length(y[:, 1])]
    function highpassfilter(signals, Fs, cutoff, order = 4)
        wdo = 2.0 * cutoff / Fs
        filth = digitalfilter(Highpass(wdo), Butterworth(order))
        filtfilt(filth, signals)
    end
    y = highpassfilter(y, Fs, low_note)
    return y, Fs
end

file = "Arpeggio"
y, Fs = get_song(file)

##################################################################################
include("Utils.jl")
using .Utils
include("SPA.jl")
using FFTW
using WAV

Wl = 2048
bin_size = Fs / Wl
hop_size = Int(floor(0.25 * Wl))
windows = Utils.hannWindow.(Utils.getWindows(y, Wl, hop_size = hop_size))

nyquist = Int(Wl / 2) + 1
fft_in = [fft(window)[1:nyquist] for window in windows]
fft_out = [zeros(ComplexF64, Wl) for n in 1:length(fft_in)]

######################################################################
# Do stuff with fft_in
for n in 1:length(fft_in)
    # Convert to polar
    r = abs.(fft_in[n])

    # Get frequencies in hz
    ω_0 = SPA.detect_pitch(fft_in[n], Fs)
    ω_end = Utils.getNote(Int(round(Utils.getNoteFromFreq(ω_0))))

    # Get ratio of difference
    Δω = (ω_end - ω_0) / ω_0

    # Decompose into regions
    maxima, minima = get_maxima_minima(r)
    regions, noise = get_regions(r, minima, ω_0, bin_size)

    # Shift frequencies
    for region in noise
        fft_out[n][region] .+= fft_in[n][region]
    end

    for region in regions
        Δbins = region[1] * Δω
        Δbins_int = Int(floor(Δbins))
        Δbins_frac = Δbins - Δbins_int

        region_intshift = region[2] .+ Δbins_int

        # Remove invalid indexes
        valid_idx = findall(x -> x >= 1 && x <= length(fft_in[n]), region_intshift)
        if length(valid_idx) > 0
            region_intshift = region_intshift[valid_idx]

            fft_region = linear_interpolation(r, region[2][valid_idx] .- Δbins_frac)

            # Phase correction
            Δϕ = ℯ^(im * Δω * hop_size)
            ϕ = ℯ .^ (im * angle.(fft_in[n][region[2][valid_idx]])) .* Δϕ

            fft_out[n][region_intshift] += fft_region .* ϕ
        end
    end
end
#####################################################################################

# Convert back to signal
for n in 1:length(fft_out)
    copyto!(fft_out[n], nyquist + 1, conj.(fft_out[n][nyquist-1:-1:2]))
end
signal = zeros(length(y) + (Wl - length(y) % Wl))
for n in 1:length(fft_out)
    range_start = (n - 1) * hop_step + 1
    range = range_start:range_start+Wl-1
    signal[range] .+= real.(ifft(fft_out[n]))
end
signal = signal[1:length(y)]
wavwrite(signal, string("Outputs/", file, "_modified.wav"), Fs = Fs)
lines(signal)

ext: jl


text:
# unit test the set up

include("../src/DMD_Util.jl");
include("../src/DMD_Type.jl");

using LowRankApprox

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
@time X, t, at, Bt = simDMD(m, n, k, T; seed=123456);

# obtain exact dmd and trap dmd estimates

dmdexactestimate(m,n,k,X,t)

@time a_exact, B_exact = dmdexactestimate(m,n,k,X,t,dmdtype="exact")
@time a_trap, B_trap = dmdexactestimate(m,n,k,X,t,dmdtype="trap")

@time a_exact_approx, B_exact_approx = dmdexactestimate(m,n,k,X,
                                                        t,dmdtype="exact",
                                                        uselowrank=true)
@time a_trap_approx, B_trap_approx = dmdexactestimate(m,n,k,X,
                                                        t,dmdtype="trap",
                                                        uselowrank=true)

@printf "%7.4e abs err evals w/ exactDMD\n" besterrperm(a_exact,at)
@printf "%7.4e abs err evals w/ trapDMD\n" besterrperm(a_trap,at)
@printf "%7.4e abs err evals w/ exactDMD (approximate)\n" besterrperm(
    a_exact_approx,at)
@printf "%7.4e abs err evals w/ trapDMD (approximate)\n" besterrperm(
    a_trap_approx,at)


text:
# Autogenerated wrapper script for SSGraphBLAS_jll for x86_64-apple-darwin
export libgraphblas

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SSGraphBLAS")
JLLWrappers.@declare_library_product(libgraphblas, "@rpath/libgraphblas.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libgraphblas,
        "lib/libgraphblas.6.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export barycenter,
       centroid,
       bounding_box

"""
    barycenter(wps::AbstractVector{WeightedPoint23})
    barycenter(wp₁::WeightedPoint23, wp₂::WeightedPoint23, wps::WeightedPoint23...)
    barycenter(ps::AbstractVector{Point23}, ws::AbstractVector{Real})

Computes the barycenter of a non-empty set of 2D or 3D weighted points.

Returns [`Point2`](@ref) or [`Point3`](@ref) depending on the dimension of the
input values.

!!! info "Precondition"

    `length(ps) > 1`, and the sum of the weights is non-zero.
"""
barycenter(wps::AbstractVector{WeightedPoint23})

barycenter(ps::AbstractVector) =
    length(ps) > 1 ?
        barycenter(collect(CxxRef{_pointfor(ps[1])}, CxxRef.(ps))) :
        error("length(ps) ≯ 1")

for D ∈ 2:3
    P = Symbol(:Point, D)
    WP = Symbol(:Weighted, P)
    @eval begin
        @cxxdereference barycenter(wp₁::$WP, wp₂::$WP,
                                   wps::reference_type_union($WP)...) =
            barycenter(CxxRef.([wp₁, wp₂, wps...]))
    end
end

barycenter(ps::AbstractVector, ws::AbstractVector) =
    length(ws) > 1 && length(ps) > 1 ?
        barycenter(collect(CxxRef{_pointfor(ps[1])}, CxxRef.(ps)),
                   iscxxtype(FT) ?
                       CxxRef.(collect(FT, convert.(FT, ws))) :
                       collect(FT, convert.(FT, ws))) :
        barycenter(ps)

"""
    bounding_box(ps::AbstractVector{Point23})
    bounding_box(p::Point23, q::Point23, ps::Point23...)

Computes the bounding box of a non-empty set of 2D or 3D points.

Returns [`IsoRectangle2`](@ref) or [`IsoCuboid3`](@ref) depending on the
dimension of the input values.

!!! info "Precondition"

    `length(ps) > 1`
"""
bounding_box(ps::AbstractVector{Point23})

bounding_box(ps::AbstractVector) =
    length(ps) > 1 ?
        bounding_box(collect(CxxRef{_pointfor(ps[1])}, CxxRef.(ps))) :
        error("length(ps) ≯ 1")

for P ∈ (Point2, Point3)
    @eval @cxxdereference bounding_box(p::$P, q::$P,
                                       ps::reference_type_union($P)...) =
        bounding_box(CxxRef.([p, q, ps...]))
end

const _centroid_T = Union{Point23
                        , Segment23
                        , Triangle23
                        , IsoBox23
                        , Tetrahedron3
                        , Circle2
                        , Sphere3}
"""
    centroid(xs::AbstractVector{T})
    centroid(x::T, xs::T...)

Computes the centroid of a non-empty set of 2D or 3D objects.

Returns [`Point2`](@ref) or [`Point3`](@ref) dependencing on the dimension of
the input objects.

For two dimensional inputs, `T` must be either [`Point2`](@ref),
[`Segment2`](@ref), [`Triangle2`](@ref), [`IsoRectangle2`](@ref), or
[`Circle2`](@ref).

For three dimensional inputs, `T` must be either [`Point3`](@ref),
[`Segment3`](@ref), [`Triangle3`](@ref), [`IsoCuboid3`](@ref),
[`Sphere3`](@ref), or [`Tetrahedron3`](@ref).

!!! info "Precondition"

    `!isempty(xs)`
"""
centroid(xs::AbstractVector{_centroid_T})

centroid(xs::AbstractVector) =
    !isempty(xs) ?
        centroid(CxxRef.(xs)) :
        error("isempty(xs)")

for T ∈ (:Segment2, :Segment3
       , :IsoRectangle2, :IsoCuboid3
       , :Circle2, :Sphere3)
    @eval @cxxdereference centroid(x::$T,
                                   xs::reference_type_union($T)...) =
        centroid(CxxRef.([x, xs...]))
end
for T ∈ (:Triangle2, :Triangle3, :Tetrahedron3)
    @eval @cxxdereference centroid(x::$T, y::$T,
                                   xs::reference_type_union($T)...) =
        centroid(CxxRef.([x, y, xs...]))
end
for P ∈ (:Point2, :Point3)
    @eval @cxxdereference centroid(a::$P, b::$P, c::$P, d::$P, e::$P,
                                   xs::reference_type_union($P)...) =
        centroid(CxxRef.([a, b, c, d, e, xs...]))
end


text:
include("../ReconstructionEvaluations.jl");

# file write convention
dt = 170227
author = "tm"
fn = "u19_graph_clustering"
dir = joinpath(homedir(), "seungmount/research/tommy/s1/data");

# load edges and seg list
edges_fn = joinpath(dir, "170215_nt_combined_edges_cons.csv")
seg_list_fn = joinpath(dir, "170202_kl_seg_list.csv");
edges = load_edges(edges_fn);
# seg_list: id, size(vx), no synapses, pre/post category (post then pre)
seg_list = load_seg_list(seg_list_fn);
seg_list = quantify_pre_post(seg_list, edges);

# filter for pre-post counts
pre_post = count_segs_both_pre_and_post(edges[:,2])
post_only_count = 10
edge_prepost_filter = create_prepost_edge_mask(edges, pre_post, post_only_count);
filtered_edges = edges[edge_prepost_filter,:];

# filter for segment size
pre_min = 1e5
pre_max = 5e7
post_min = 7e5
post_max = 5e8
pre_filter = (seg_list[:,4] .== 1) & (pre_min .<= seg_list[:,2] .<= pre_max);
post_filter = (seg_list[:,4] .== 2) & (post_min .<= seg_list[:,2] .<= post_max);
seg_filter = pre_filter | post_filter
syn_to_segs, segs_to_syn, syn_coords, syn_size = create_graph_dicts(filtered_edges, seg_list[seg_filter,1])
segs = sort(unique(vcat(values(syn_to_segs)...)))
seg_to_index, index_to_seg = create_index_dict(segs);
pre_to_post, post_to_pre, pre_to_syn, post_to_syn = create_post_pre_syn_dicts(syn_to_segs);

# cluster
adj = create_adjacency_matrix(seg_to_index, syn_to_segs, syn_size);
reweight_adj!(adj) # apply log to the synapse weights
rand_order = randperm(size(adj,1))
radj = adj[rand_order,rand_order]
# Randomize adjacency matrix for better clustering results
radj_fn = joinpath(dir, "$(dt)_$(author)_rand_adj_sparse_$(fn).csv")
write_sparse(radj_fn, radj);

# Cluster the randomized matrix
cluster, louvain_order = louvain_clustering(radj);

# Read/write the permutations of the matrix sorting so the original segment IDs can be recovered
perm_fn = joinpath(dir, "$(dt)_$(author)_cluster_perm_$(fn).csv")
perm_tbl = hcat(collect(1:size(radj,1)), rand_order, louvain_order, round(Int64, cluster))
writedlm(perm_fn, perm_tbl)

# Sort the randomized matrix by its clustering order
cadj = radj[louvain_order, louvain_order]
cadj_fn = joinpath(dir, "$(dt)_$(author)_clustered_adj_sparse_$(fn).csv")
write_sparse(cadj_fn, cadj)


text:
using Documenter, TestThings

makedocs(;
    modules=[TestThings],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Foo" => "foo.md"
    ],
    repo="https://github.com/under-Peter/TestThings.jl/blob/{commit}{path}#L{line}",
    sitename="TestThings.jl",
    authors="Andreas Peter",
)

deploydocs(;
    repo="github.com/under-Peter/TestThings.jl",
)


text:
# simple test that things get created
@testset "RSS gen" begin
    f = joinpath(p, "basic", "feed.xml")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""<rss version="2.0" xmlns:atom="http://www.w3.org/2005/Atom">""", fc)
    @test occursin(raw"""<title>JuDoc Template</title>""", fc)
    @test occursin(raw"""<description><![CDATA[Example website using JuDoc
  ]]></description>""", fc)
    @test !occursin(raw"""<author>""", fc)
    @test occursin(raw"""<link>https://tlienart.github.io/JuDocTemplates.jl/pub/menu1.html</link>""", fc)
    @test occursin(raw"""<description><![CDATA[A short description of the page which would serve as blurb in a RSS feed;""", fc)
end


text:
function solveMUMPS(A::SparseMatrixCSC{T1}, rhs::AbstractArray{T2}, sym::Int = 0, ooc::Int = 0, tr::Int = 0) where {T1, T2}
    x = zeros(promote_type(T1, T2), size(rhs))
    return solveMUMPS!(A, rhs, x, sym, ooc, tr)
end

function solveMUMPS!(A::SparseMatrixCSC, rhs::AbstractArray, x::Array, sym::Int=0, ooc::Int=0, tr::Int=0)
    # factorization
    factor = factorMUMPS(A, sym, ooc)
    # solve system
    x = applyMUMPS!(factor, rhs, x, tr)
    # free memory
    destroyMUMPS(factor)
    return x
end

function factorMUMPS(A::SparseMatrixCSC{ComplexF64}, sym::Int=0, ooc::Int=0)
    # Generate LU-factorization of complex matrix A, 'sym' can be: 0=unsymmetric, 1=symm. pos def, 2=general symmetric
    if size(A, 1) != size(A, 2)
    error("factorMUMPS: Matrix must be square!")
    end
    n  = size(A, 1)
    mumpsstat = [0]
    start_time = time_ns()
    p  = ccall((:factor_mumps_cmplx_, MUMPSlibPath),
               Int64,
               (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}),
               Ref(n), Ref(sym), Ref(ooc),  convert(Ptr{ComplexF64}, pointer(A.nzval)), A.rowval, A.colptr, mumpsstat)
    checkMUMPSerror(mumpsstat)
    facTime = float(time_ns() - start_time)
    factor = MUMPSfactorization(p, myid(), n, A.nzval[1], facTime)
    return factor
end

function factorMUMPS(A::SparseMatrixCSC{Float64}, sym::Int=0, ooc::Int=0)
    # Generate LU-factorization of a real matrix A, 'sym' can be: 0=unsymmetric, 1=symm. pos def, 2=general symmetric
    if size(A, 1) != size(A, 2)
        error("factorMUMPS: Matrix must be square!")
    end
    n  = size(A, 1)
    mumpsstat = [0]
    start_time = time_ns()
    p  = ccall((:factor_mumps_, MUMPSlibPath),
               Int64,
               (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}),
               Ref(n), Ref(sym), Ref(ooc), A.nzval, A.rowval, A.colptr, mumpsstat)
    checkMUMPSerror(mumpsstat)
    facTime = float(time_ns() - start_time)
    factor = MUMPSfactorization(p, myid(), n, A.nzval[1], facTime)
    return factor
end

function checkMUMPSerror(mumpsstat)
    # check if MUMPS reported error
    if mumpsstat[1] == -10
        error("MUMPS: Numerically singular matrix.")
    elseif mumpsstat[1] == -13
        error("MUMPS: memory allocation error")
    elseif mumpsstat[1] == -40
        error("MUMPS: matrix is not positive definite")
    elseif mumpsstat[1] == -90
        error("MUMPS: Error in out-of-core management.Probably there is not enough disk space.")
    elseif mumpsstat[1] < 0
        error(@sprintf("MUMPS: error --> %d <--. Please refer to Ch. 7 of MUMPS User's guide!", mumpsstat[1]))
    end
end

function applyMUMPS(factor::MUMPSfactorization{T1}, rhs::AbstractArray{T2, N}, x::Array=zeros(promote_type(T1, T2), size(rhs)), tr::Int=0) where {T1, T2, N}
    id1 = myid()
    id2 = factor.worker
    if id1 != id2
        @warn "Worker $id1 has no access to MUMPS factorization stored on $id2. Trying to remotecall!"
        return remotecall_fetch(applyMUMPS, factor.worker, factor, rhs, x, tr)::Array{promote_type(T1, T2), N}
    end

    if size(rhs, 1) != factor.n
        error("applyMUMPS: wrong size of rhs, size(A)=$(factor.n), size(rhs)=$(size(rhs, 1)) x $(size(rhs, 2)).")
    end

    if size(x) != size(rhs)
        if isempty(x)
            x = zeros(promote_type(T1, T2), size(rhs))
        else
            error("applyMUMPS: wrong size of x, size(A)=$(factor.n), size(rhs)=$(size(rhs)), size(x)=$(size(x)) provided")
        end
    end

    return applyMUMPS!(factor, rhs, x, tr)::Array{promote_type(T1, T2), N}
end

function applyMUMPS!(factor::MUMPSfactorization{Float64}, rhs::Array{Float64}, x::Array{Float64}, tr::Int=0)
    nrhs = size(rhs, 2)
    ptr = factor.ptr
    ccall((:solve_mumps_, MUMPSlibPath),
          Int64,
          (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}),
          Ref(ptr), Ref(nrhs), rhs, x, Ref(tr))
    return x
end

function applyMUMPS!(factor::MUMPSfactorization{Float64}, rhs::SparseMatrixCSC{Float64}, x::Array{Float64, 2}, tr::Int=0)
    nrhs = size(rhs, 2)
    nzrhs = nnz(rhs)
    ptr = factor.ptr
    ccall((:solve_mumps_sparse_rhs_, MUMPSlibPath),
          Nothing,
          (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}),
          Ref(ptr), Ref(nzrhs), Ref(nrhs), rhs.nzval, rhs.rowval, rhs.colptr, x, Ref(tr))
    return x
end

function applyMUMPS!(factor::MUMPSfactorization{Float64}, rhs::SparseVector{Float64}, x::Array{Float64, 1}, tr::Int=0)
    nrhs = 1
    nzrhs = nnz(rhs)
    colptr = [1; nzrhs+1]
    ptr = factor.ptr
    ccall((:solve_mumps_sparse_rhs_, MUMPSlibPath),
          Nothing,
          (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}),
          Ref(ptr), Ref(nzrhs), Ref(nrhs), rhs.nzval, rhs.nzind, colptr, x, Ref(tr))
    return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64}, rhs::Array{ComplexF64}, x::Array{ComplexF64}, tr::Int=0)
    n = size(rhs, 1)
    nrhs = size(rhs, 2)
    ptr = factor.ptr
    ccall((:solve_mumps_cmplx_, MUMPSlibPath),
          Int64,
          (Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF32}, Ptr{Int64}),
          Ref(ptr), Ref(nrhs), convert(Ptr{ComplexF64}, pointer(rhs)), convert(Ptr{ComplexF64}, pointer(x)), Ref(tr))
    return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64}, rhs::SparseMatrixCSC{ComplexF64}, x::Array{ComplexF64, 2}, tr::Int=0)
    nrhs = size(rhs, 2)
    nzrhs = nnz(rhs)
    ptr = factor.ptr
    ccall((:solve_mumps_cmplx_sparse_rhs_, MUMPSlibPath),
          Nothing,
          (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{Int64}),
          Ref(ptr), Ref(nzrhs),Ref(nrhs), convert(Ptr{ComplexF64}, pointer(rhs.nzval)), rhs.rowval, rhs.colptr, convert(Ptr{ComplexF64}, pointer(x)), Ref(tr))
    return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64}, rhs::SparseVector{ComplexF64}, x::Array{ComplexF64, 1}, tr::Int=0)
    nrhs = 1
    nzrhs = nnz(rhs)
    colptr = [1; nzrhs + 1]
    ptr = factor.ptr
    ccall((:solve_mumps_cmplx_sparse_rhs_, MUMPSlibPath),
          Nothing,
          (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{Int64}),
          Ref(ptr), Ref(nzrhs),Ref(nrhs), convert(Ptr{ComplexF64}, pointer(rhs.nzval)), rhs.nzind, colptr, convert(Ptr{ComplexF64}, pointer(x)), Ref(tr))
    return x
end

function destroyMUMPS(factor::MUMPSfactorization{Float64})
     #  free memory
    id1 = myid()
    id2 = factor.worker
    if myid() != factor.worker
        remotecall_fetch(destroyMUMPS, id2, factor)
        return
    end
    ccall((:destroy_mumps_, MUMPSlibPath), Int64, (Ptr{Int64}, ), Ref(factor.ptr))
    factor.ptr = -1
    factor.n = -1
    factor.time = -1.0
end

function destroyMUMPS(factor::MUMPSfactorization{ComplexF64})
     #  free memory
    id1 = myid()
    id2 = factor.worker
    if myid() != factor.worker
        remotecall_fetch(destroyMUMPS, id2, factor)
        return
    end
    ccall((:destroy_mumps_cmplx_, MUMPSlibPath), Int64, (Ptr{Int64}, ), Ref(factor.ptr))
    factor.ptr = -1
    factor.n = -1
    factor.time = -1.0
end


text:
# default polynomial type
struct MPoly{T} <: AbstractPolynomial
    terms::T
end
#MPoly() = MPoly(EMPTY_TERMS)
#MPoly(x::AbstractTerm) = MPoly([x])
#MPoly(x::M) where {M<:AbstractMonomial} = MPoly(Term(x))
MPoly{T}(x::AbstractTerm) where T = MPoly([x])
MPoly{T}(x::M) where {T,M<:AbstractMonomial} = MPoly(Term(x))
MPoly{T}(x::CoeffType) where {T} = MPoly(eltype(T)(x))
terms(x::MPoly) = x.terms

Base.copy(x::MPoly) = MPoly(map(copy, terms(x)))
Base.one(::Type{<:MPoly{T}}) where T = MPoly(eltype(T)(1))

Base.:*(x::AbstractTerm, p::MPoly) = p * x
function Base.:*(p::MPoly, x::T) where {T<:AbstractTerm}
    if iszero(p)
        return p
    elseif iszero(x)
        return MPoly(emptyterm(T))
    elseif isone(x)
        return p
    else
        MPoly(T[t * x for t in terms(p) if !iszero(t)])
    end
end
function Base.:*(p::MPoly, x::MPoly)
    ts = terms(x)
    out = similar(terms(p), length(ts) * length(terms(p)))
    current_len = 0
    s = MPoly(out)
    for tp in terms(p)
        start = 1
        for tx in terms(x)
            _, (current_len, start) = _add!(s, tp * tx, current_len, start)
        end
    end
    resize!(out, current_len)
    (debugmode() && !issorted(terms(s), rev=true)) && throw("Polynomial not sorted!")

    return s
end
Base.:+(p::MPoly, x::AbstractTerm) = add!(copy(p), x)
Base.:-(p::MPoly, x::AbstractTerm) = sub!(copy(p), x)

Base.:-(p::MPoly) = -1 * p
sub!(p::MPoly, x::AbstractTerm) = add!(p, -x)

addcoef(x::T, c) where {T<:AbstractTerm} = (c += coeff(x); return iszero(c), T(c, monomial(x)))
addcoef(x::T, c::T) where {T<:AbstractTerm} = addcoef(x, c.coeff)
subcoef(x::T, c) where {T<:AbstractTerm} = (c = coeff(x) - c; return iszero(c), T(c, monomial(x)))
subcoef(x::T, c::T) where {T<:AbstractTerm} = subcoef(x, c.coeff)

function add!(p::MPoly, x::AbstractTerm)
    q, (current_len, _) = _add!(p, x, length(terms(p)), 1);
    resize!(terms(q), current_len)
    (debugmode() && !issorted(terms(q), rev=true)) && throw("Polynomial not sorted!")
    q
end

@inline function _add!(p::MPoly, x::AbstractTerm, current_len, start)
    iszero(x) && return p, (current_len, 1)
    ts = terms(p)
    i = start
    @inbounds for i in start:current_len
        t = ts[i]
        if t < x
            current_len += 1
            if current_len > length(ts)
                insert!(ts, i, x)
            else
                @simd for j in current_len:-1:i+1
                    ts[j] = ts[j-1]
                end
                ts[i] = x
            end
            return p, (current_len, i)
        elseif ismatch(t, x)
            return p, _add_term_matched(ts, i, t, x, current_len)
        end
    end
    current_len += 1
    if current_len > length(ts)
        push!(ts, x)
    else
        ts[current_len] = x
    end
    return p, (current_len, i)
end

@inline function _add_term_matched(ts, i, t, x, current_len)
    iz, t = addcoef(t, x)
    if iz
        tsp = pointer(ts) + i*Base.aligned_sizeof(eltype(ts))
        GC.@preserve ts Base.unsafe_copyto!(tsp - Base.aligned_sizeof(eltype(ts)), tsp, current_len - i)
        current_len -= 1
        return current_len, max(1, i-1)
    else
        @inbounds ts[i] = t
        return current_len, i
    end
end

@inline function _add_rev!(p::MPoly, x::AbstractTerm, _end=1)
    iszero(x) && return p, 1
    ts = terms(p)
    i = _end
    @inbounds for i in length(ts):-1:_end
        t = ts[i]
        if t > x
            insert!(ts, i+1, x)
            return p, i
        elseif ismatch(t, x)
            return p, _add_term_matched_rev(ts, i, t, x)
        end
    end
    push!(ts, x)
    return p, i
end

function _add_term_matched_rev(ts, i, t, x)
    iz, t = addcoef(t, x)
    if iz
        deleteat!(ts, i)
        return max(1, i-1)
    else
        ts[i] = t
        return i
    end
end

function add!(p::AbstractPolynomial, x::AbstractPolynomial)
    tp = terms(p)
    tx = terms(x)
    if tp isa Vector
        sizehint!(tp, length(tp) + length(tx))
    end
    for t in tx
        add!(p, t)
    end
    return p
end
function sub!(p::AbstractPolynomial, x::AbstractPolynomial)
    tp = terms(p)
    tx = terms(x)
    if tp isa Vector
        sizehint!(tp, length(tp) + length(tx))
    end
    for t in tx
        sub!(p, t)
    end
    return p
end
Base.:+(p::AbstractPolynomial, x::AbstractPolynomial) = add!(copy(p), x)
Base.:-(p::AbstractPolynomial, x::AbstractPolynomial) = sub!(copy(p), x)

lt(p::AbstractPolynomial) = first(terms(p))
lc(p::AbstractPolynomial) = coeff(lt(p))
function rmlt!(p::MPoly)
    ts = terms(p)
    popfirst!(ts)
    return p
end
function takelt!(p::MPoly, x::MPoly)
    add!(p, lt(x))
    rmlt!(x)
    return p
end

function fnmadd!(p::MPoly, l, s)
    for li in terms(l)
        sub!(p, li * s)
    end
end

function Base.divrem(p::MPoly, d::MPoly)
    p = copy(p)
    q = MPoly(similar(terms(p), 0))
    r = MPoly(similar(terms(p), 0))
    while !isempty(terms(p))
        nx, fail = lt(p) / lt(d)
        if fail
            takelt!(r, p)
        else
            #p -= d * nx
            #q += nx
            fnmadd!(p, d, nx)
            add!(q, nx)
        end
    end
    return q, r
end

function divexact(p::MPoly, d::MPoly)
    p = copy(p)
    ts = similar(terms(p), 0)
    sizehint!(ts, length(terms(d))-1)
    q = MPoly(ts)
    _end = 1
    while !isempty(terms(p))
        nx, fail = lt(p) / lt(d)
        @assert !fail
        fnmadd!(p, d, nx)
        _, _end = _add_rev!(q, nx, _end)
    end
    (debugmode() && !issorted(terms(q), rev=true)) && throw("Polynomial not sorted!")
    return q
end

function Base.rem(p::AbstractPolynomial, d::AbstractPolynomial)
    p = copy(p)
    r = MPoly(similar(terms(p), 0))
    while !isempty(terms(p))
        nx, fail = lt(p) / lt(d)
        if fail
            takelt!(r, p)
        else
            sub!(p, d * nx)
        end
    end
    return r
end

# TODO
Base.:(/)(x::MPoly, y::MPoly) = divexact(x, y)

Base.gcd(x::Union{AbstractTerm,AbstractMonomial}, y::MPoly) = gcd(MPoly(x), y)
Base.gcd(x::MPoly, y::Union{AbstractTerm,AbstractMonomial}) = gcd(y, x)
function Base.gcd(x::MPoly, y::MPoly)
    # trival case
    if iszero(x) || isone(y)
        return y
    elseif iszero(y) || isone(x) || x == y
        return x
    end

    v1, p1 = to_univariate(x)
    v2, p2 = to_univariate(y)
    if v1 < v2
        x, y = y, x
        v1, v2 = v2, v1
        p1, p2 = p2, p1
    end
    # v2 < v1
    # both are constants
    v2 == NOT_A_VAR && return MPoly(gcd(lt(x), lt(y)))
    if v2 < v1
        # `v2` in p2 doesn't exist in `x`, so the gcd at this level is 1 and we
        # just move on to the next level
        return gcd(x, content(p2))
    end
    v1 == v2 || error("unreachable")

    g = gcd(p1, p2)
    return univariate_to_multivariate(g)
end

function pick_var(x::MPoly)
    ts = terms(x)
    v = NOT_A_VAR
    for (i, t) in enumerate(ts)
        if degree(t) > 0
            m = monomial(t)
            vv = firstid(m) # get the minvar
            if vv < v
                v = vv
            end
        end
    end
    return IDType(v)
end

function to_univariate(x::MPoly)
    v = pick_var(x)
    v, (v == NOT_A_VAR ? nothing : SparsePoly(x, v))
end


text:
using BinaryBuilder

name = "Rmath"
version = v"0.3.0"

sources = [
    GitSource("https://github.com/JuliaStats/Rmath-julia.git",
              "c71948a8ce7b8b1aca7b94a3e97840bc06c51749"),
]

script = raw"""
cd $WORKSPACE/srcdir/Rmath-julia*
make -j${nproc}
mkdir -p "${libdir}"
mv src/libRmath-julia.* "${libdir}"
"""

platforms = supported_platforms(;experimental=true)

products = [
    LibraryProduct("libRmath-julia", :libRmath),
]

dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
__precompile__()

module Binning
    using Missings

    # import Base: ==, !=, >, <, >=, <=, +, -, *, !, &, |, ⊻, ^, /

    export equaldistancelimits, find_bin, binindex, binsupport

    include("BinningUtils.jl")
end


text:
using NeighbourLists
using NeighbourLists: SMat, SVec
using Test
using LinearAlgebra

include("nn_list.jl")

Ns = [3, 3, 4, 4, 5, 5]
@info("Testing PairList Correctness: ")
for N in Ns
   C = SMat( diagm(0 => (2.0 .+ 0.2 * rand(3))) * N )
   X = [ C' * rand(SVec)   for i = 1:ceil(Int, abs(det(C))) ÷ 4 + 2 ]
   pbc = SVec(rand(Bool, 3))
   cutoff = 2.0

   # compute a cell list
   nlist = PairList(X, cutoff, C, pbc)
   @show typeof(nlist)

   # compute a NearestNeighbors list
   i, j, r, R = nn_list(X, cutoff, C, pbc)
   R = X[j] - X[i]
   first = NeighbourLists.get_first(i, length(X))
   NeighbourLists.sort_neigs!(j, (r, R), first)

   println(@test (nlist.i == i) && (nlist.j == j))

   # check that they are sorted
   print("sorted: ");
   println(@test all( issorted(j)  for (_i, j, _R) in NeighbourLists.sites(nlist) ))

   # check that the neighbourhoods produced by `neigs` are correct
   pass_neigs = true
   for (i, j, R) in NeighbourLists.sites(nlist)
      j_, R_ = NeighbourLists.neigs(nlist, i)
      if j != j_ || !(R ≈ R_)
         pass_neigs = false
         break
      end
      js, Rs, Ss = NeighbourLists.neigss(nlist, i)
      if j != js || !(R ≈ Rs)
         pass_neigs = false
         break
      end
   end
   print("neigs: "); println(@test pass_neigs)
end
println()


text:
"""
    SlackCallRecord

This object has fields

* `call::Symbol`: either `:DummyThread` or `:send_exception_message`
* `args::Tuple`: the arguments to the call
* `kwargs::NamedTuple`: any keyword arguments to the call

corresponding to function calls made with [`DummyThread`](@ref)s.
"""
struct SlackCallRecord
    call::Symbol
    args::Tuple
    kwargs::NamedTuple
end

StructTypes.StructType(::Type{SlackCallRecord}) = StructTypes.Struct()
# This makes it so JSON3 doesn't omit tuples/namedtuples if they are empty
@inline StructTypes.isempty(::Type{SlackCallRecord}, ::NamedTuple) = false
@inline StructTypes.isempty(::Type{SlackCallRecord}, ::Tuple) = false

"""
    DummyThread <: AbstractSlackThread

Provides a "dummy" SlackThread which does not make any API calls
to Slack, and instead logs messages to it's `logged` field.

This can be used for testing or to pass to code expecting an `AbstractSlackThread`
when you don't want to log anything to Slack.

The `logged` field is a `Vector{SlackCallRecord}` corresponding to logging calls.

See also: [`SlackCallRecord`](@ref).
"""
struct DummyThread <: AbstractSlackThread
    logged::Vector{SlackCallRecord}
end

# Allow "setting" `ts` or `channel` for interchangability with `SlackThread`
function Base.setproperty!(d::DummyThread, name::Symbol, x::Any)
    if name === :ts || name === :channel
        if x !== nothing
            # throw error if we can't convert, just like a SlackThread would
            convert(String, x)
        end
        return x
    end
    return setfield!(d, name, x) # get the usual error
end

# Return `nothing` if asked for `ts` or `channel`
function Base.getproperty(d::DummyThread, name::Symbol)
    if name === :ts || name === :channel
        return nothing
    end
    return getfield(d, name)
end

Base.propertynames(::DummyThread) = (:channel, :ts, :logged)

DummyThread() = DummyThread([])

StructTypes.StructType(::Type{DummyThread}) = StructTypes.Struct()

function (d::DummyThread)(args...; kwargs...)
    push!(d.logged, SlackCallRecord(:DummyThread, args, NamedTuple(kwargs)))
    return nothing
end

function send_exception_message(d::DummyThread, msg)
    push!(d.logged, SlackCallRecord(:send_exception_message, tuple(msg), NamedTuple()))
    return nothing
end


text:
# NOTE: we now have an AbstractODEProblem
abstract type AbstractODEProblem end

struct ODEProblem{F,T<:Tuple{Number,Number},U<:AbstractVector,P<:AbstractVector} <: AbstractODEProblem
    f::F
    tspan::T
    u0::U
    θ::P
end

abstract type ODESolver end

struct Euler{T} <: ODESolver
    dt::T
end

function (solver::Euler)(prob::AbstractODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    (u + dt*f(u,θ), t+dt)
end
struct RK2{T} <: ODESolver
    dt::T
end
function (solver::RK2)(prob::AbstractODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    du = f(u,θ)
    uh = u + du*dt
    u + dt/2*(du + f(uh,θ)), t+dt
end

function solve(prob::AbstractODEProblem, solver::ODESolver)
    t = prob.tspan[1]; u = prob.u0
    us = [u]; ts = [t]
    while t < prob.tspan[2]
        (u,t) = solver(prob, u, t)
        push!(us,u)
        push!(ts,t)
    end
    # NOTE: we are not `reduce(hcat,us)`ing any more
    ts, us
end


text:
using NSGAIII
using vOptGeneric, GLPKMathProgInterface

# using PyPlot

function plot_pop(P)
    clf()
    pop = filter(x -> x.CV ≈ 0, P)
    non_dom = NSGAIII.fast_non_dominated_sort!(P)[1]
    pop = setdiff(pop, non_dom)
    p = plot3D(map(x -> x.y[1], pop), map(x -> x.y[2], pop), map(x -> x.y[3], pop), "bo", markersize=1)
    p = plot3D(map(x -> x.y[1], non_dom), map(x -> x.y[2], non_dom), map(x -> x.y[3], non_dom), "go", markersize=1)
    ax = gca()
    ax[:set_xlim]([-95, 110])
    ax[:set_ylim]([-60, 45])
    ax[:set_zlim]([-50, 110])
    !isinteractive() && show()
    sleep(0.1)
end

m = vModel(solver = GLPKSolverMIP())

@variable(m, 0 <=x[1:5] <= 10)
@variable(m, δ[1:3], Bin)

@addobjective(m, Max, dot([17,-12,-12,-19,-6], x) + dot([-73, -99, -81], δ))
@addobjective(m, Max, dot([2,-6,0,-12,13], x) + dot([-61,-79,-53], δ))
@addobjective(m, Max, dot([-20,7,-16,0,-1], x) + dot([-72,-54,-79], δ))

@constraint(m, sum(δ) <= 1)
@constraint(m, -x[2] + 6x[5] + 25δ[1] <= 52)
@constraint(m, -x[1] + 18x[4] + 18x[5] + 8δ[2] <= 77)
@constraint(m, 7x[4] + 9x[5] + 19δ[3] <= 66)
@constraint(m, 16x[1] + 20x[5] <= 86)
@constraint(m, 13x[2] + 7x[4] <= 86)

print(m)

res = nsga(0, 100, 10, m,  pmut=0.3);

solve(m, method=:lex)
println("\nRésolution lexico-graphique : ")
display(getY_N(m));

println()
println("Meilleur individu sur le premier objectif")
x1 = sort(res, by = x -> x[2][1])[end];
println("x = $(x1[1][1:5]) , δ = $(x1[1][6:8])")
println("z = $(x1[2])")
println("CV : $(x1[3])")

println("Meilleur individu sur le deuxième objectif")
x2 = sort(res, by = x -> x[2][2])[end];
println("x = $(x2[1][1:5]) , δ = $(x2[1][6:8])")
println("z = $(x2[2])")
println("CV : $(x2[3])")

println("Meilleur individu sur le troisième objectif")
x3 = sort(res, by = x -> x[2][3])[end];
println("x = $(x3[1][1:5]) , δ = $(x3[1][6:8])")
println("z = $(x3[2])")
println("CV : $(x3[3])")


text:
#!/usr/bin/env python
# coding: utf-8

# # Spontaneous moralization on path graph

# ## Loading modules

using QSWalk
using LightGraphs
using PyPlot

# ## Non-symmetric case

# Here we provide more advanced functions corresponding to sponatneous moralization. Below we start with simple parametrizations. Note *dim* should be odd for unique middle-point.

dim = 101
midpoint = ceil(Int, dim/2)
w = 0.5
timepoint = 40.
adjacency = adjacency_matrix(PathGraph(dim));

# We generate all of the operators needed for the evolution, including initial state.

lind, vset = nm_lind(adjacency)
hglobal = nm_glob_ham(adjacency)
hlocal = nm_loc_ham(vset)
opnonsymmetric = evolve_generator(hglobal, [lind], hlocal, w)

rhoinit = nm_init(vset[[midpoint]], vset);

# Finally we make an evolution.

rho = evolve(opnonsymmetric, rhoinit, timepoint);

# Note that first moment of natural measurement deviates from zero.

positions = (collect(1:dim).-midpoint)
measurement_nonsymmetric = nm_measurement(rho, vset)
println("First moment centralized in midpoint: $(sum(positions .* nm_measurement(rho, vset)))")

# This is because of the non-symmetrices *lind* choice (analysis shows, that even removing *hlocal* and *hglobal* operators does not result in symmetric evolution). To confirm this, we plot natural measurement and its reflection around *n=midpoint*.

plot(positions, measurement_nonsymmetric, "k")
plot(positions, reverse(measurement_nonsymmetric), "b")
xlabel("position")
ylabel("probability")
axis([positions[1], positions[end], 0, maximum(measurement_nonsymmetric)])
vlines(0, 0., maximum(measurement_nonsymmetric), linestyles="--")

# ## Symmetric case

# The way to correct this is to choose another, symmetric Lindblad operator. While standard suage of *nm_lind* will always output the same result, we can choose different basic orthogonal matrices to form different operators. In following example. We choose dictionary, which for different vertex degree chooses different matrix.

linddescription1 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; 1 -1])
linddescription2 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; -1 1])
lind1, vset = nm_lind(adjacency, linddescription1)
lind2, vset = nm_lind(adjacency, linddescription2);

# We can make similar creation for each vertex. For example one can choose. We restrict ourselves to *lind1* and *lind2*, as those guarantees symmetrization.

vset = make_vertex_set(adjacency)
linddescription3 = Dict(v=>rand(length(v), length(v)) for v = vlist(vset))
lind3, _ = nm_lind(adjacency, linddescription3);

# Other functions should be adjusted to use both *lind1* and *lind2*.

hglobal = nm_glob_ham(adjacency)
hlocal = nm_loc_ham(vset)
opsymmetric = evolve_generator(hglobal, [lind1, lind2], hlocal, w)

rhoinit = nm_init(vset[[midpoint]], vset)
rho = evolve(opsymmetric, rhoinit, timepoint);

# Now both first momement and the distribution confirms symmetric evolution.

positions = (collect(1:dim).-midpoint)
measurement_symmetric = nm_measurement(rho, vset)
println("First moment centralized in midpoint: $(sum(positions .* measurement_symmetric))")

plot(positions, measurement_symmetric, "k")
plot(positions, reverse(measurement_symmetric), "b")
xlabel("position")
ylabel("probability")
axis([positions[1], positions[end], 0, maximum(measurement_symmetric)])
vlines(0, 0., maximum(measurement_symmetric), linestyles="--")


text:
module task

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
v = [1, 2, 3, 4]
dv = DataArray(v, falses(size(v)))

m = [1 2; 3 4]
dm = DataArray(m, falses(size(m)))

t = Array(Int, 2, 2, 2)
t[1:2, 1:2, 1:2] = 1
dt = DataArray(t, falses(size(t)))

dv = DataArray(v)
dv = DataArray(v, [false, false, false, false])

dv = DataArray(dv)
dm = DataArray(dm)
dt = DataArray(dt)

dv = DataArray(Int, 2)
dm = DataArray(Int, 2, 2)
dt = DataArray(Int, 2, 2, 2)

similar(dv)
similar(dm)
similar(dt)

similar(dv, 2)
similar(dm, 2, 2)
similar(dt, 2, 2, 2)

# Test vecbind
a = [1:4]
d = DataArray(a)
@assert isequal(vecbind(a,a), [a,a])
@assert isequal(vecbind(a,1.0 * a), 1.0 * [a,a])
@assert isequal(vecbind(d,a), [d,d])
@assert isequal(vecbind(a,d), [d,d])
@assert isequal(vecbind(1.0 * a,d), 1.0 * [d,d])
@assert isequal(vecbind(d,d), [d,d])
@assert isequal(vecbind(a,IndexedVector(a)), [a,a])
@assert isequal(vecbind(a,IndexedVector(d)), [d,d])
@assert isequal(vecbind(PooledDataArray(a),IndexedVector(d)), [d,d])
@assert isequal(vecbind(PooledDataArray(a),IndexedVector(a)), [a,a])
@assert isequal(vecbind(a,RepeatedVector(a,2)), [a,a,a])
@assert isequal(vecbind(a,StackedVector({a,1.0*a})), [a,a,a])


text:
using MultidimensionalSpectroscopy, PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles, Random

# make sure to set script directory as pwd()
cd(@__DIR__)

pygui(true)

calc_2d = true

cmp = create_colormap("bright");

Δ = -0.2
wc_k(θ) = (2.0 + Δ + sin(θ))  * 2 * pi # cavity frequency
wc = wc_k(0)
wa = 2.0  * 2 * pi  # atom frequency
g  = 0.15 * 2 * pi   # coupling strength
κ = 0.5         # cavity dissipation rate
γ = 0.00005         # atom dissipation rate
N = 1               # number of cavity fock states
n_th_a = 0.0        # avg number of thermal bath excitation ???
#use_rwa = True

tlist = [0:0.05:50;]
b_fock = FockBasis(N)
b_atom = NLevelBasis(2)
Psi0 = fockstate(b_fock,1) ⊗ nlevelstate(b_atom,1)    # start with an excited cavity

a = destroy(b_fock) ⊗ one(b_atom)
at = dagger(a)
sm = one(b_fock) ⊗ transition(b_atom,1,2)
sp = dagger(sm)

H = wc * at * a + wa * sp * sm + g * (at * sm + a * sp)

H_k(θ) = wc_k(θ) * at * a + wa * sp * sm + g * (at * sm + a * sp)

println(H)

#J = [sqrt(γ)*sm,sqrt(κ)*a]
# the same ?
J     = [a,sm]
rates = [sqrt(κ), sqrt(γ)]
#rates  = [1, 1]
tout, rhot = timeevolution.master(tlist,Psi0,H,rates .* J)
nc = real(expect(at * a, rhot))
na = real(expect(sp * sm, rhot))

figure(figsize=(6,6))
subplot(211)
plot(tout,nc)
plot(tout,na)

#ω, spec = timecorrelations.correlation2spectrum(tlist, nc; normalize_spec=true)
#figure(figsize=(6,3))
#plot(ω, spec)

rho0 = dm(Psi0)

corr = timecorrelations.correlation(tlist, rho0, H, J, at, a; rates=rates)
ω, spec = timecorrelations.correlation2spectrum(tlist, corr; normalize_spec=true)
plot(tlist,corr)

subplot(212)
plot(ω, spec)

width = .5
ωq_TDBC     = [H.data[3,3] - H.data[2,2]]
ampl_TDBC   = [1] / (N)
gam = .1

try
    local m = @which spectral_density(1)
    Base.delete_method(m)
catch
end
function spectral_density(ω) # power spectral density, thermal bath noise spectrum
    w = [(ω .- ωq_TDBC[i]) ./ (.5 * width) for i in 1:length(ωq_TDBC)]
    if ω == 0
        #Jω = 0
        Jω = gam
    elseif !isless(real(ω),0)
        Jω = sum([ampl_TDBC[i] ./ (1 .+ w[i].^2) for i in 1:length(ωq_TDBC)])  # probably most realistic ... ???
        #Jω = gam / 2 * (ω / (2 * pi)) * (ω > 0.0)
    elseif  isless(real(ω),0)
        Jω = 0
        #Jω = gam / 2 * (-ω / (2 * pi)) * (ω < 0.0)
    end
    return Jω
end

Γ_R = .01
a_ops = [Γ_R * sp * sm, spectral_density]# plot correlation function and spectrum
a_ops = [Γ_R * at * a,  spectral_density]# plot correlation function and spectrum

R, ekets = timeevolution.bloch_redfield_tensor(H, [a_ops]; J=rates .* J)

μ12 = at * sm*sp + a * sm*sp
μ23 = at * a * sp + at * a * sm

tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0,R,H)
corr = expect(μ12,rhot)

corr    = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)
ωnew, spec = timecorrelations.correlation2spectrum(tnew, corr; normalize_spec=true);
plot(ωnew,spec)

rhot = nothing
tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0*μ12,R,H)

vecs   = eigvecs(dense(H).data)
eivecs = [Ket(H.basis_l,vecs[:,i]) for i in 1:length(vecs[1,:])]

#Operator(H.basis_l, H.basis_l, eivecs)

transf = copy(H)
transf.data = vecs

rhot_coupl = [transf' * i * transf for i in rhot]

H_coupl = transf' * H * transf

tout, rhot = timeevolution.master_bloch_redfield(tlist,μ12*rho0*μ12,R,H)

vecs   = eigvecs(dense(H_coupl).data)
eivecs = [Ket(H_coupl.basis_l,vecs[:,i]) for i in 1:length(vecs[1,:])]

figure()
subplot(132);
title("eigen/coupled basis") # levels: 1 -> GS, 2 -> LP, 3 - N-1 -> ..., N -> UP
for i in 1:length(H.basis_l)
    es_n = real(expect(dm(eivecs[i]),rhot))
    plot(tout,es_n,linestyle="dashed", label=string(i)); #append!(leg,[string(i)])
    ylim((-0.1, 1))
end
legend()

vecs   = eigvecs(dense(H).data)
eivecs = [Ket(H.basis_l,vecs[:,i]) for i in 1:length(vecs[1,:])]

es_mat = zeros(length(tout))
subplot(133); title("uncoupled basis") # levels: 1 -> GS, 2 -> cav., 3 - N -> TLSs
for i in 1:length(H.basis_l)
    es_n = real(expect(dm(eivecs[i]),rhot))
    plot(tout,es_n,linestyle="dashed", label=string(i)); #append!(leg,[string(i)])
    if i > 2
        global es_mat = es_mat .+ es_n
    end
end
#plot(tout, es_mat, label="Σ(TLSs)")
legend()

corr = []
spec = []
angles = collect(-pi/4:0.001*pi:pi/2)
for i in angles
    #temp1 = timecorrelations.correlation(tlist, rho0, H_k(i), J, at, a; rates=rates)
    local R, ekets = timeevolution.bloch_redfield_tensor(H_k(i), [a_ops]; J=rates .* J)
    local tout, temp = timeevolution.master_bloch_redfield(tlist,μ12*rho0,R,H_k(i))
    temp1 = expect(μ12,temp)
    append!(corr, [temp1])
    local ω, temp2 = timecorrelations.correlation2spectrum(tlist, temp1; normalize_spec=false)
    append!(spec,[temp2])
end

spec = hcat(spec...)
figure()
contourf(angles ./ (2*pi) .* 360,ω, spec)

Psi0 = fockstate(b_fock,0) ⊗ nlevelstate(b_atom,1)    # start with system in gs
rho0 = dm(Psi0)

use_sub = false

if calc_2d
        ## calculate (complex) 3rd order corr function (with T=0)
        zp = 11 # zeropad up to 2^zp

        F = R

        ## calculate 2d spectra at
        T = [0., 10., 100.]

        spectra2d = Array{out2d}(undef, length(T))
        for i = 1:length(T)
            spectra2d[i] = make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],
                                                "redfield";debug=false,use_sub=use_sub,zp=zp);
        end

        ## crop 2D data and increase dw
            spectra2d = [crop2d(spectra2d[i],3.5;w_max=20,step=2) for i = 1:length(T)]

        ## assign ω-axis from output
        ω = spectra2d[1].ω

        ## plot 2D spectra for each(?) T
        # what to plot
        rep="absorptive"
        scal="lin"

        ## make  subplot layout
        nplots = length(T); ncols = ceil(sqrt(nplots)); nrows = ceil(nplots / ncols)

        ## create figure
        figure(figsize=(ncols*3+0.25,nrows*3+0.5));
        suptitle(rep * " 2D spectrum (" * scal * ". scaling)")
        for i = 1:nplots
                subplot(nrows,ncols,i,aspect="equal")
                plot2d(ω,spectra2d[i].full2d;repr=rep,scaling=scal)
                title("2D spectrum at $(T[i]) fs")
                #xlim([0, E₁+E₂]); ylim([0, E₁+E₂]);
        end
        tight_layout()

        #plot([E₁, E₁], [ω[1], ω[end]],"k--",linewidth=1,alpha=0.25)
        #plot([E₂, E₂], [ω[1], ω[end]],"k--",linewidth=1,alpha=0.25)
        #plot([ω[1], ω[end]],[E₁, E₁],"k--",linewidth=1,alpha=0.25)
        #plot([ω[1], ω[end]],[E₂, E₂],"k--",linewidth=1,alpha=0.25)
        #plot([energies[2:3], energies[2:3]], [ω[1], ω[end]],"g--",linewidth=1,alpha=0.25)
        #plot([ω[1], ω[end]],[energies[2:3], energies[2:3]],"g--",linewidth=1,alpha=0.25)

        # choose one of the following plotting ranges
        #xlim([ω[1], ω[end]]); ylim([ω[1], ω[end]]);
        #xlim([0, ω[end]]); ylim([0, ω[end]]);
        #xlim([0, wa+wc]); ylim([0, wa+wc]);
end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
immutable DynamicSystem
    f::Function # Transition function (function of time, state, control, disturbance)
    U::Function # "Gain" function (function of time, state, control, disturbance)
    Ubar::Function # Terminal value (function of state)
    T::Int # Terminal time
    statedim::Int # Dimension of state
    amin::Function  # Control minimum value (function of time and state)
    amax::Function  # Control maximum value (function of time and state)
    W::Tuple     # Disturbance (draws independent random variables)
end

immutable DynamicSystem1D{T1<:Real}
    f::Function    # Transition function
    U::Function    # "Gain" function
    Ubar::Function # Terminal value
    T::Int         # Terminal time
    amin::T1       # Control minimum value
    amax::T1       # Control maximum value
end

type DynamicSystemTrajectory1D{T1<:Real}
    system::DynamicSystem1D
    state::Vector{T1} # Observed states (t=0,\dots,T)
    control::Vector{T1} # Applied controls (t=0,\dots,T-1)
    value::Vector{T1} # Achieved value (t=0,\dots,T), value[0] = 0
end

function DynamicSystemTrajectory1D{T<:Real}(system::DynamicSystem1D{T})
    DynamicSystemTrajectory1D(system, Vector{T}(), Vector{T}(), Vector{T}())
end

type DynamicSystemTrajectory{T1<:Real}
    system::DynamicSystem
    state::Array{T1,2} # Observed states (t=0,\dots,T)
    control::Array{T1,2} # Applied controls (t=0,\dots,T-1)
    value::Vector{T1} # Achieved value (t=0,\dots,T), value[t=0] = 0
end

function DynamicSystemTrajectory(system::DynamicSystem)
    DynamicSystemTrajectory(system, Array{Float64,2}(system.statedim, system.T+1),
                            Array{Float64,2}(system.statedim, system.T),
                            Vector{Float64}(system.T+1))
end

abstract type AbstractSystemController1D end
abstract type OnlineSystemController1D <: AbstractSystemController1D end

immutable OfflineSystemControl1D{T1<:Real} <: AbstractSystemController1D
    system::DynamicSystem1D{T1}
    initialstate::T1
    policy::Function
end

function OfflineSystemControl1D{T1<:Real}(system::DynamicSystem1D{T1},
                                          initialstate::T1,
                                          xarr::Vector{T1},
                                          αvec::Array{T1,2})
    αinterp = interpolate((xarr,0:system.T-1), αvec, Gridded(Linear()))
    α(t, x) = αinterp[x, t]
    OfflineSystemControl1D(system, initialstate, α)
end

immutable MPCSystem1D{T1<:Real} <: OnlineSystemController1D
    system::DynamicSystem1D{T1}
    initialstate::T1
    numscenarios::Vector{Int} # Number of scenarios for given time step
end

function MPCSystem1D{T1<:Real}(system::DynamicSystem1D{T1}, initialstate::T1, numscenarios::Int)
    MPCSystem1D(system, initialstate, fill(numscenarios,system.T))
end

immutable OLFCSystem1D{T1<:Real} <: OnlineSystemController1D
    system::DynamicSystem1D{T1}
    initialstate::T1
    numscenarios::Vector{Int} # Number of scenarios for given time step
end

function OLFCSystem1D{T1<:Real}(system::DynamicSystem1D{T1}, initialstate::T1, numscenarios::Int)
    OLFCSystem1D(system, initialstate, fill(numscenarios,system.T))
end

immutable OfflineSystemControl
    system::DynamicSystem
    policy::Function
end

function OfflineSystemControl{T1<:Real,N<:Int}(system::DynamicSystem,
                                               xtup,
                                               αvec::Array{T1,N})
    αinterp = interpolate(reverse((0:system.T-1, xtup...)),
                          αvec, Gridded(Linear()))
    α(t, x) = αinterp[reverse(x)..., t]
    OfflineSystemControl1D(system, α)
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Function: calculate_transcription_control_array
# Description: Calculate the transcriptional control array at time t
# Generated on: 2020-02-05T11:38:24.149
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{String,Any} => Dictionary holding model parameters
#
# Output arguments:
# control_array::Array{Float64,1} => Transcriptional control array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_transcription_control_array(t::Float64,x::Array{Float64,1},data_dictionary::Dict{String,Any})

	# initialize the control -
	control_array = zeros(2)

	# Alias the species -
	deGFP = x[1]
	sigma_70 = x[2]
	mRNA_deGFP = x[3]
	mRNA_sigma_70 = x[4]
	protein_deGFP = x[5]
	protein_sigma_70 = x[6]

	# Alias the binding parameters -
	binding_parameter_dictionary = data_dictionary["binding_parameter_dictionary"]
	n_deGFP_sigma_70 = binding_parameter_dictionary["n_deGFP_sigma_70"]
	K_deGFP_sigma_70 = binding_parameter_dictionary["K_deGFP_sigma_70"]

	# Alias the control function parameters -
	control_parameter_dictionary = data_dictionary["control_parameter_dictionary"]
	W_deGFP_RNAP = control_parameter_dictionary["W_deGFP_RNAP"]
	W_deGFP_sigma_70 = control_parameter_dictionary["W_deGFP_sigma_70"]
	W_sigma_70_RNAP = control_parameter_dictionary["W_sigma_70_RNAP"]

	# basline -
	BLW = data_dictionary["base_line_weight"]

	# Transfer function target:deGFP actor:sigma_70
	actor_set_deGFP_sigma_70 = [
		protein_sigma_70
	]
	actor = prod(actor_set_deGFP_sigma_70)
	b_deGFP_sigma_70 = (actor^(n_deGFP_sigma_70))/(K_deGFP_sigma_70^(n_deGFP_sigma_70)+actor^(n_deGFP_sigma_70))

	# Control function for deGFP -
	control_array[1] = (W_deGFP_RNAP+W_deGFP_sigma_70*b_deGFP_sigma_70)/(1.0+W_deGFP_RNAP+W_deGFP_sigma_70*b_deGFP_sigma_70)

	# Control function for sigma_70 -
	control_array[2] = (W_sigma_70_RNAP)/(1.0+W_sigma_70_RNAP)

	# return -
	return control_array
end

#
# ----------------------------------------------------------------------------------- #
# Function: calculate_translation_control_array
# Description: Calculate the translation control array at time t
# Generated on: 2020-02-05T11:38:24.152
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{String,Any} => Dictionary holding model parameters
#
# Output arguments:
# control_array::Array{Float64,1} => Translation control array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_translation_control_array(t::Float64,x::Array{Float64,1},data_dictionary::Dict{String,Any})

	# initialize the control -
	control_array = ones(2)

	# correct for "resource"?
    correction_term = (x[7]/100.0)
    control_array = control_array*correction_term

	# return -
	return control_array
end

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
import SpecialFunctions
import SpecialFunctions: besselj

@testset "Bessel Roots" begin
    tol = 1e-11

    # Check if besselj(nu, besselroots(nu, n) ) is small:

    for nu = 0.:0.1:5.
        let n = 10
            @test norm( besselj.(nu, besselroots(nu, n) ), Inf ) < tol
        end
    end
end


text:
__precompile__()
module Spelling
using DataStructures
using StringDistances

export score, spellcheck, spellcheck_sentence, correct, correct_sentence

const dir = Pkg.dir("Spelling")
const linuxwords = Set(readlines("$dir/data/words"))
include("$dir/data/guten.jl")
include("$dir/data/tv.jl")

#Uppercase/lowercase isn't very well-supported, though there are some cased words in the
# lists, such as country names and pronouns, hence why no automatic lowercase is being performed

const formality = 0.6
const cutoff    = 0.5
"""
score(s)

Scores a given word in terms of how common it is. A value of 1 means it is a member of the
linux.words file provided with fedora 26, a higher value is its relative frequency per
billion words in text, based on the ~40,000 most common words from project Gutenberg
and the ~40,000 most common words from TV show scripts (both from Wikipedia).
"""
function score(s)
    g = s∈ keys(guten) ? guten[s] : 0.0
    t = s∈ keys(tv) ? tv[s] : 0.0
    l = s∈ linuxwords ? 1 : 0
    sqrt(formality*g^2 + (1-formality)*t^2 + l)
end

"""
spellcheck(s, c=cutoff)

Returns true if a word scores above the cutoff value (i.e. is a word) and false otherwise.
"""
spellcheck(s, c=cutoff) = score(s) >= c

"""
spellcheck_sentence(s::String, c=cutoff)

Spellchecks a sentence and returns regex matches for the failing words.
"""
function spellcheck_sentence(s::String, c=cutoff)
    words = collect(eachmatch(r"\w+", s))
    words[(w->!spellcheck(w.match, c)).(words)]
end

function correct_tv(s::String)
    ks = [keys(tv)...]
    ns = (x->StringDistances.compare(DamerauLevenshtein(), x, s)).(ks)
    ks[indmax(ns)]
end

function correct_guten(s::String)
    ks = [keys(guten)...]
    ns = (x->StringDistances.compare(DamerauLevenshtein(), x, s)).(ks)
    ks[indmax(ns)]
end

function correct(s::String, f=formality)
    tks = [keys(tv)...]
    gks = [keys(guten)...]
    tns = (x->StringDistances.compare(DamerauLevenshtein(), x, s)).(tks)
    gns = (x->StringDistances.compare(DamerauLevenshtein(), x, s)).(gks)
    tm, gm = indmax(tns), indmax(gns)
    if (1-f)*tns[tm] == f*gns[gm]
        tv[tks[tm]] > guten[gks[gm]] ? tks[tm] : gks[gm]
    else
        (1-f)*tns[tm] > f*gns[gm] ? tks[tm] : gks[gm]
    end
end

function correct_sentence(s::String, f=formality)
    words = String.((x->x.match).(collect(eachmatch(r"\w+", s))))
    sn = split(s, r"\w+")
    nw = [correct.(words, f); ""]
    join(sn.*nw)
end

end


text:
function categorical_ticks!(ax, x, y)
	cat_xticks = categorical_ticks(x)
	cat_yticks = categorical_ticks(y)

	if !isnothing(cat_xticks)
		ax.xticks[] = cat_xticks
	end
	if !isnothing(cat_yticks)
		ax.yticks[] = cat_yticks
	end
	ax
end

function categorical_ticks(var)
	ticks  = categorical_range(var)
	labels = categorical_labels(var)

	if ticks isa Automatic
		nothing
	else
		ticks, labels
	end
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
using Revise
using Random
using SpecialFunctions
using PyCall
using PyPlot
using StatsPlots

function mySliceSampler(pdf_log, x_0, w, m, N, burnIn = 1000)
    # REQUIRES MODULES:
    # Random
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_init --> inital vector of x (dimension 1xD), where D dimension of pdf
    # w --> typical slice size, for a window (dimension 1xD)
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    # burnIn --> number of burn-in samples (optional, results will be discarded)
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations
    pdflog_x_1 = 0.0;

    # to be able to input scalars:
    if typeof(x_0) == Float64
        x_0 = [x_0]
    end

    for ii in 1:(N+burnIn)
        # Update the new x progressively for each dimension:
        L   = 1*x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   = 1*x_0
        x_1 = 1*x_0  # "1*" needed, otherwise x_0 suddendly takes mysterious values
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            if (ii == 1)
                #vertical = rand() * pdf_log(x_0) # in normal space
                z = pdf_log(x_0) + log(rand()) # in log-space
            else
                # saving evaluations, pdflog_x_1 was updated for current permutation in the end
                z = pdflog_x_1 + log(rand())
            end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            L[d] = x_0[d] - rand() * w[d]
            R[d] = L[d] + w[d]

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L[d] -= w[d]
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R[d] += w[d]
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                x_1[d] = L[d] + rand() * (R[d] - L[d])

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0[d] = x_1[d]     # update value for this dimension
                    break
                # Value was not within slice, shrink the interval:
                elseif (x_1[d] < x_0[d])
                    L[d] = x_1[d]
                elseif (x_1[d] > x_0[d])
                    R[d] = x_1[d]
                else
                    throw(ErrorException("Error during shrinking the interval"))
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        # Just overwrite the burnIn by using modular arithmetic:
        i = 1 + (ii-1)%N    # gives 1,2,3,...N,1,2,3,...,(N+burnIn)
        x_s[i,:] = x_1    # sample point, after one round of permutations
        pdflog_x_s[i] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations

    end

    return x_s, pdflog_x_s
end

#%%
################################################################

function log_beta(θ,a::Any,b::Any)
    #using SpecialFunctions
    θ = θ[1]
    if 0 < θ < 1
        # Numerator
        num = log(θ)*(a-1) + log(1-θ)*(b-1)
        # The denominator (Beta-function) can be integrated or made up of Gamma-functions
        # Since we only treat integer values here, Γ(a) = (a-1)!
        #den = log(factorial(a-1) * factorial(b-1) / factorial(a+b-1))
        den = log(SpecialFunctions.beta(a,b))
        return num - den
    else
        return -Inf
    end
end

################################################################

# Log-Bernoulli PDF:
function log_bernoulli_pdf(θ,y)
    if 0 < θ < 1 && (y==0 || y==1)    # Only admissible values
        return y * log(θ) + (1-y) * log(1-θ)
    else
        return -Inf
    end
end

# Log-Likelihood fct of a (series of) coin tosses with bias:
function log_likelihood_fct(θ, y)
    # Make array s.t. multiple theta values can be input:
    #likelihood = Array{Float64}(undef,length(θ))
    # Do that for every theta value requested:
    #for ii in 1:length(θ)
    #    # Elementwise application of the log_beronoulli_pdf over outcome vector y, then SUM together:
    #    likelihood[ii] = sum(log_bernoulli_pdf.(θ[ii],y))
    #end

    # Un array-ise the output:
    #if length(θ) == 1
    #    likelihood = likelihood[1]
    #end

    θ = θ[1]
    likelihood = sum(log_bernoulli_pdf.(θ,y))

    return likelihood
end

#%%
################################################################

# Chain inputs:
N = 10^6
x_start = 0.5;
w = .5;  # typical window width
m = 100; # maximum number of window widths

# Final posteriors in one array:
posteriors = Array{Float64}(undef,N,3);

#%%
# First Column of Figure 6.4
a = 100
b = 100
y = [ones(1,17) zeros(1,3)]

log_prior(θ) = log_beta(θ,a,b)
log_lklhd(θ) = log_likelihood_fct(θ,y)
log_posterior(θ) = log_lklhd(θ) + log_prior(θ)

# Need to use "StatsPlots" or "Plots" histogram to add the density curve, PyPlot won't work
chain, pdf = mySliceSampler(log_prior,x_start,w,m,N);
histogram(chain, bins=100, normalize=:pdf, label="prior")   # Comes from StatsPlots now
density!(chain,linewidth=3,label="prior")

chain, pdf = mySliceSampler(log_lklhd,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="likelihood")
density!(chain,linewidth=3,label="likelihood")

chain, pdf = mySliceSampler(log_posterior,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="posterior")
density!(chain,linewidth=3,label="posterior")

posteriors[:,1] = chain # save the posterior

Plots.savefig("A4-A-C1_Slice.pdf")

#%%
# Second Column
a = 18.25
b = 6.75
y = [ones(1,17) zeros(1,3)]

log_prior(θ) = log_beta(θ,a,b)
log_lklhd(θ) = log_likelihood_fct(θ,y)
log_posterior(θ) = log_lklhd(θ) + log_prior(θ)

chain, pdf = mySliceSampler(log_prior,x_start,w,m,N);
histogram(chain, bins=100, normalize=:pdf, label="prior")   # Comes from StatsPlots now
density!(chain,linewidth=3,label="prior")

chain, pdf = mySliceSampler(log_lklhd,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="likelihood")
density!(chain,linewidth=3,label="likelihood")

chain, pdf = mySliceSampler(log_posterior,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="posterior")
density!(chain,linewidth=3,label="posterior")

posteriors[:,2] = chain # save the posterior

Plots.savefig("A4-A-C2_Slice.pdf")

#%%
# Third Column
a = 1
b = 1
y = [ones(1,17) zeros(1,3)]

log_prior(θ) = log_beta(θ,a,b)
log_lklhd(θ) = log_likelihood_fct(θ,y)
log_posterior(θ) = log_lklhd(θ) + log_prior(θ)

chain, pdf = mySliceSampler(log_prior,x_start,w,m,N);
histogram(chain, bins=100, normalize=:pdf, label="prior")   # Comes from StatsPlots now
density!(chain,linewidth=3,label="prior")

chain, pdf = mySliceSampler(log_lklhd,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="likelihood")
density!(chain,linewidth=3,label="likelihood")

chain, pdf = mySliceSampler(log_posterior,x_start,w,m,N)
histogram!(chain, bins=100, normalize=:pdf,label="posterior")
density!(chain,linewidth=3,label="posterior")

posteriors[:,3] = chain # save the posterior

Plots.savefig("A4-A-C3_Slice.pdf")

#%%

density(posteriors, label=["a=100, b=100","a=18.25, b=6.75","a=1, b=1"])
Plots.savefig("A4-A-all_Slice.pdf")


text:
module MusicManipulations

using Reexport
@reexport using MIDI
import MIDI: Notes, Note

include("general.jl")
include("midifiles.jl")
include("quantize.jl")
include("scale.jl")
include("humanize.jl")

include("data_handling/data_extraction.jl")
include("data_handling/timeseries.jl")

include("specific_modules/jazz.jl")
include("specific_modules/drums.jl")

include("motifs/notes_with_motifs.jl")

export Jazz, Drums

end


text:
module TestChebyshev

using IterativeSolvers
using Test
using Random
using LinearAlgebra

function randSPD(T, n)
    A = rand(T, n, n) + n * I
    A' * A
end

function approx_eigenvalue_bounds(A)
    λs = eigvals(A)
    mnv, mxv = minimum(real(λs)), maximum(real(λs))
    Δ = (mxv - mnv) / 100
    mnv - Δ, mxv + Δ
end

#Chebyshev
@testset "Chebyshev" begin

n = 10
Random.seed!(1234321)

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = randSPD(T, n)
    b = rand(T, n)
    reltol = √(eps(real(T)))
    abstol = reltol

    @testset "Without preconditioner" begin
        λ_min, λ_max = approx_eigenvalue_bounds(A)
        x0 = rand(n)
        x, history = chebyshev(A, b, λ_min, λ_max, reltol=reltol, maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ reltol
    end

    @testset "With an initial guess" begin
        λ_min, λ_max = approx_eigenvalue_bounds(A)
        x0 = rand(T, n)
        initial_residual = norm(A * x0 - b)
        x, history = chebyshev!(x0, A, b, λ_min, λ_max, reltol=reltol, maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test x == x0
        @test norm(A * x - b) ≤ reltol * initial_residual

        x0 = rand(T, n)
        x, history = chebyshev!(x0, A, b, λ_min, λ_max, abstol=abstol, reltol=zero(real(T)), maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test x == x0
        @test norm(A * x - b) ≤ 2*abstol
    end

    @testset "With a preconditioner" begin
        B = randSPD(T, n)
        B_fact = cholesky!(B, Val(false))
        λ_min, λ_max = approx_eigenvalue_bounds(B_fact \ A)
        x, history = chebyshev(A, b, λ_min, λ_max, Pl = B_fact, reltol=reltol, maxiter=10n, log=true)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ reltol
    end
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]
        λ_min, λ_max = approx_eigenvalue_bounds(A)

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = chebyshev!(x, A, b, λ_min, λ_max, log=true)
        @test 2 ≤ niters(ch) ≤ n

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = chebyshev!(x, A, b, λ_min, λ_max, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end
end

end # module


text:
using Stats
using Base.Test

@test_approx_eq mean([1, 2, 3]) 2.0
@test_approx_eq mean(1:3) 2.0

@test_approx_eq geomean([1, 2, 3]) (6.0)^(1/3)
@test_approx_eq geomean(1:3) (6.0)^(1/3)
@test_approx_eq geomean([2, 8]) 4.0
@test_approx_eq geomean([4, 1, 1/32]) 0.5
@test geomean([1, 0, 2]) == 0.0

@test_approx_eq harmmean([1, 2, 3]) 3 / (1 + 1/2 + 1/3)
@test_approx_eq harmmean(1:3) 3 / (1 + 1/2 + 1/3)
@test_approx_eq harmmean([1, 2, 4]) 12 / 7

@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.0) 22.4
@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.4) 4.0
@test_approx_eq trimmean([-100, 2, 3, 7, 200], 0.8) 3.0

@test_approx_eq wmean([1.0, 2.0, 3.0], [1/3, 1/3, 1/3]) 2.0
@test_approx_eq wmean([1.0, 2.0, 3.0], [1.0, 0.0, 0.0]) 1.0
@test_approx_eq wmean([1.0, 2.0, 3.0], [0.0, 1.0, 0.0]) 2.0
@test_approx_eq wmean([1.0, 2.0, 3.0], [0.0, 0.0, 1.0]) 3.0
@test_approx_eq wmean([1.0, 2.0, 3.0], [0.5, 0.0, 0.5]) 2.0
@test_approx_eq wmean([1.0, 2.0, 3.0], [0.5, 0.5, 0.0]) 1.5
@test_approx_eq wmean([1.0, 2.0, 3.0], [0.0, 0.5, 0.5]) 2.5

@test_approx_eq wmean(1:3, [1/3, 1/3, 1/3]) 2.0
@test_approx_eq wmean(1:3, [1.0, 0.0, 0.0]) 1.0
@test_approx_eq wmean(1:3, [0.0, 1.0, 0.0]) 2.0
@test_approx_eq wmean(1:3, [0.0, 0.0, 1.0]) 3.0
@test_approx_eq wmean(1:3, [0.5, 0.0, 0.5]) 2.0
@test_approx_eq wmean(1:3, [0.5, 0.5, 0.0]) 1.5
@test_approx_eq wmean(1:3, [0.0, 0.5, 0.5]) 2.5

@test_approx_eq mean(1:3, weights([1.0, 1.0, 0.5])) 1.8


text:
using HorizonSideRobots

function mark_cells(r)
    for side in 0:3
        digits = move_through_barrier!(r, side)
        putmarker!(r)
        back!(r, digits, side)
    end
end

function move_through_barrier!(r, side)
    digits = []
    while (length(digits)==0 || digits[length(digits)]!=0)
        count1 = 0
        count2 = 0
        count3 = 0
        while !isborder(r, HorizonSide(side))
            move!(r, HorizonSide(side))
            count1+=1
        end
        while (isborder(r, HorizonSide(side)) && !isborder(r, HorizonSide((side+1)%4)))
            move!(r, HorizonSide((side+1)%4))
            count2 += 1
        end
        if !isborder(r, HorizonSide(side))
            move!(r, HorizonSide(side))
            count3 += 1
            while isborder(r, HorizonSide((side+3)%4))
                move!(r, HorizonSide(side))
                count3 += 1
            end
            for i in 1:count2
                move!(r, HorizonSide((side+3)%4))
            end
        else
            for i in 1:count2
                move!(r, HorizonSide((side+3)%4))
            end
            count2 = 0
        end
        push!(digits,count1)
        push!(digits,count2)
        push!(digits,count3)
        push!(digits,count2)
    end
    return digits
end

function movements!(r, side, count)
    for i in 0:count-1
        move!(r, HorizonSide(side))
    end
end

function back!(r, digits, side)
    digits = reverse(digits)
    for i in 1:length(digits)
        rside = side
        if i%2 == 1
            rside = (rside+i)%4
        else
            rside = (rside+2)%4
        end
        movements!(r, rside, digits[i])
    end
end

mark_cells(r)

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
include("pristine_graphene.jl")

"""
    Defect

Abstract type for all defect types.
"""
abstract type Defect end

"""
    LocalSpin(x::Float64, y::Float64, z::Float64, coord::GrapheneCoord) <: Defect

A local spin with components `x`, `y`, and `z` located at `coord`. See
[`GrapheneCoord`](@ref).
"""
struct LocalSpin <: Defect
    x::Float64
    y::Float64
    z::Float64
    coord::GrapheneCoord
end

"""
    ImpurityState(ϵ::Float64, coupling::Vector{Tuple{Float64,GrapheneCoord}}) <: Defect

An impurity state of energy `ϵ` (in eV) coupled to the graphene system. The
tuples in the `coupling` field contain all the coupling energies (in eV) and the
 corresponding [`GrapheneCoord`](@ref)'s.
"""
struct ImpurityState <: Defect
    ϵ::Float64
    coupling::Vector{Tuple{Float64,GrapheneCoord}}
end

"""
    Hopping(c1::GrapheneCoord, c2::GrapheneCoord, Δ::ComplexF64) <: Defect

Hopping modification (in eV) between two [`GrapheneCoord`](@ref)'s. If `c1==c2`,
 this quantity corresponds to a local energy modification.
"""
struct Hopping <: Defect
    c1::GrapheneCoord
    c2::GrapheneCoord
    Δ::ComplexF64
end

"""
    GrapheneSystem(
        μ::Float64,
        T::Float64,
        Δ::Array{ComplexF64,2},
        V::Array{Float64,2},
        scattering_states::Vector{GrapheneState},
        imps::Vector{Float64},
    )

A structure describing the perturbed graphene system.

See also [`mkGrapheneSystem`](@ref).
"""
struct GrapheneSystem
    μ::Float64                                  # Chemical potential
    T::Float64                                  # Temperature
    Δ::Array{ComplexF64,2}                      # Δ matrix
    V::Array{Float64,2}                         # V Matrix
    scattering_states::Vector{GrapheneState}    # List of all perturbed states
    imps::Vector{Float64}                       # Impurity energies
end

"""
    mkGrapheneSystem(
        μ::Float64,
        T::Float64,
        defects::Vector{Defect},
    )

Construct [`GrapheneSystem`](@ref).

# Arguments
* `μ`: chemical potential
* `T`: temperature
* `defects`: a list of [`Defect`]'s

When supplying []`Hopping`](@ref) in `defects`, ensure that each coordinate pair
appears only once as repeated pairs with different couplings will cause eariler
values to be overwritten. The order of the coordinates does not matter.

The function constructs a [`GrapheneSystem`](@ref) with the `Δ` and `V` matrices.
 In addition, a list of all
[`GrapheneState`](@ref) that are perturbed (`scattering_states` field in
[`GrapheneCoord`](@ref)) and a list of
impurity energies (`imps` field in [`GrapheneCoord`](@ref)) are included.
"""
function mkGrapheneSystem(μ::Float64, T::Float64, defects::Vector{Defect})

    hops = filter(x -> typeof(x) == Hopping, defects)
    spins = filter(x -> typeof(x) == LocalSpin, defects)
    imps = filter(x -> typeof(x) == ImpurityState, defects)

    # Create a coupling dictionary
    if isempty(hops)
        hops_dict = Dict{Tuple{GrapheneState,GrapheneState},ComplexF64}()
    else
        hops_dict =
            reduce(
                vcat,
                [
                    [
                        (
                            (
                                GrapheneState(x.c1, SpinUp),
                                GrapheneState(x.c2, SpinUp),
                            ),
                            x.Δ,
                        ),
                        (
                            (
                                GrapheneState(x.c2, SpinUp),
                                GrapheneState(x.c1, SpinUp),
                            ),
                            conj(x.Δ),
                        ),
                        (
                            (
                                GrapheneState(x.c1, SpinDown),
                                GrapheneState(x.c2, SpinDown),
                            ),
                            x.Δ,
                        ),
                        (
                            (
                                GrapheneState(x.c2, SpinDown),
                                GrapheneState(x.c1, SpinDown),
                            ),
                            conj(x.Δ),
                        ),
                    ] for x in hops
                ],
            ) |> Dict
    end

    if isempty(spins)
        spins_dict = Dict{Tuple{GrapheneState,GrapheneState},ComplexF64}()
    else
        spins_dict =
            reduce(
                vcat,
                [
                    [
                        (
                            (
                                GrapheneState(s.coord, SpinUp),
                                GrapheneState(s.coord, SpinUp),
                            ),
                            s.z,
                        ),
                        (
                            (
                                GrapheneState(s.coord, SpinDown),
                                GrapheneState(s.coord, SpinDown),
                            ),
                            -s.z,
                        ),
                        (
                            (
                                GrapheneState(s.coord, SpinUp),
                                GrapheneState(s.coord, SpinDown),
                            ),
                            s.x - 1im * s.y,
                        ),
                        (
                            (
                                GrapheneState(s.coord, SpinDown),
                                GrapheneState(s.coord, SpinUp),
                            ),
                            s.x + 1im * s.y,
                        ),
                    ] for s in spins
                ],
            ) |> Dict
    end

    pert_dict = mergewith(+, hops_dict, spins_dict)

    # Get the coordinates of all the directly-perturbed atoms
    perturbed_atoms =
        [x[1].coord for x in pert_dict |> keys |> collect] |> unique
    # Get the coordinates of atoms coupled to impurities
    coupled_atoms = [
        y[2] for y in reduce(
            vcat,
            [x.coupling for x in imps],
            init = Vector{Tuple{Float64,GrapheneCoord}}[],
        )
    ]
    # Combine the two sets of coordinates and keep only unique entries
    all_atoms = vcat(perturbed_atoms, coupled_atoms) |> unique |> sort
    all_states = vcat(
        [GrapheneState(x, SpinUp) for x in all_atoms],
        [GrapheneState(x, SpinDown) for x in all_atoms],
    )

    # Assemble Δ
    Δ = [
        get!(pert_dict, (x, y), 0.0 + 0.0im) for x in all_states,
        y in all_states
    ]

    # Assemble V
    V_per_spin = [
        [
            (
                sum([
                    ((atom == c[2]) * c[1]) for c in imp.coupling
                ])::Float64
            ) for atom in all_atoms
        ] for imp in imps
    ]

    if !isempty(V_per_spin)
        V_per_spin = reduce(hcat, V_per_spin)
        V = [
            1 0
            0 1
        ] ⊗ V_per_spin |> collect
    else
        V = Array{Float64}(undef, 0, 0)
    end

    imp_energies = [x.ϵ for x in imps]
    return GrapheneSystem(μ, T, Δ, V, all_states, imp_energies)
end

"""
    peierls_phase(vec_pot, a1::GrapheneCoord, a2::GrapheneCoord)

Calculate the phase used in the Peierls substitution to include the effects of
the magnetic field.

To make the units work out better, the magnetic
field ``\\mathbf{B}(x, y, z) = \\Phi_0/\\mathcal{V} \\mathbf{f}(x, y, z)``,
where ``\\mathcal{V}`` is the area of the graphene unit cell in Å², ``\\Phi_0 = h / 2e``
is the magnetic flux quantum, and ``\\mathbf{f}(x, y, z) =
\\nabla \\times \\mathbf{g}(x,y,z)`` is a dimensionless vector function. Note that
``\\mathbf{f} = 1`` produces a field of about 40000 T.

For the vector potential, we have
``\\mathbf{A}(x, y, z) =\\Phi_0 / \\mathcal{V} \\mathbf{g}(x, y, z)``, where
``\\mathbf{g}(x, y, z)`` has the units of Å. Using the definition of the Peierls
phase, one gets

```math
\\phi = -\\frac{\\pi}{ \\Phi_0} \\int \\mathbf{A}\\cdot d\\mathbf{l} =
-\\frac{\\pi}{ \\mathcal{V}} \\int  \\mathbf{g}(x, y, z) \\cdot d\\mathbf{l}
\\rightarrow
-\\frac{\\pi}{ \\mathcal{V}} \\int  \\mathbf{g}_{xy}(x, y) \\cdot d\\mathbf{l}\\,.
```

The last step follows from the fact that the graphene system resides in the ``xy``
plane, so one needs to retain only the ``x`` and ``y`` componends of ``\\mathbf{g}``,
as denoted by the subscript ``xy``.

# Arguments
* `vec_pot(x,y)`: ``\\mathbf{g}_{xy}(x, y)`` in Å with `x` and `y` in Å.
* `a1`: [`GrapheneCoord`](@ref) of the "from" atom.
* `a2`: [`GrapheneCoord`](@ref) of the "to" atom.

`vec_pot(x,y)` needs to return a tuple corresponding to the vector potential in
    ``x`` and ``y`` directions.
"""
function peierls_phase(vec_pot, a1::GrapheneCoord, a2::GrapheneCoord)
    a1_loc = crystal_to_cartesian(a1)
    a2_loc = crystal_to_cartesian(a2)
    disp = a2_loc - a1_loc
    res = quadgk(
        u ->
            vec_pot(a1_loc[1] + u * disp[1], a1_loc[2] + u * disp[2]) .* disp |> sum,
        0,
        1,
    )[1]
    return (-π * res / UC_area)
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# AffineDecisionFunction #
# ========================== #
MA.mutability(::Type{<:TypedDecisionLike}) = MA.IsMutable()

function MA.mutable_copy(f::AffineDecisionFunction)
    return AffineDecisionFunction(MA.mutable_copy(f.variable_part),
                                  MA.mutable_copy(f.decision_part))
end
function MA.mutable_copy(f::QuadraticDecisionFunction)
    return QuadraticDecisionFunction(
        MA.mutable_copy(f.variable_part),
        MA.mutable_copy(f.decision_part),
        MA.mutable_copy(f.cross_terms))
end

function MA.isequal_canonical(f::F, g::F) where F<:Union{AffineDecisionFunction, VectorAffineDecisionFunction}
    return MA.isequal_canonical(f.variable_part, g.variable_part) &&
        MA.isequal_canonical(f.decision_part, g.decision_part)
end

function MA.isequal_canonical(f::F, g::F) where F <: QuadraticDecisionFunction
    return MA.isequal_canonical(f.variable_part, g.variable_part) &&
        MA.isequal_canonical(f.decision_part, g.decision_part) &&
        MA.isequal_canonical(f.cross_terms, g.cross_terms)
end

function MA.iszero!(f::TypedScalarDecisionLike)
    return iszero!(MOI.constant(f)) && MOIU._is_constant(MOIU.canonicalize!(f))
end

function MA.scaling(f::TypedScalarDecisionLike)
    return MA.scaling(f.variable_part)
end

MA.promote_operation(::Union{typeof(zero), typeof(one)}, F::Type{<:TypedScalarDecisionLike}) = F

const PROMOTE_IMPLEMENTED_OP = Union{typeof(+), typeof(-), typeof(*), typeof(/)}
function MA.promote_operation(op::PROMOTE_IMPLEMENTED_OP,
                              F::Type{<:ScalarDecisionLike{T}},
                              G::Type{<:ScalarDecisionLike{T}}) where T
    MOIU.promote_operation(op, T, F, G)
end
function MA.promote_operation(op::PROMOTE_IMPLEMENTED_OP, F::Type{T},
                              G::Type{<:TypedDecisionLike{T}}) where T
    MOIU.promote_operation(op, T, F, G)
end
function MA.promote_operation(op::PROMOTE_IMPLEMENTED_OP, F::Type{<:TypedDecisionLike{T}},
                              G::Type{T}) where T
    MOIU.promote_operation(op, T, F, G)
end
function MA.promote_operation(op::PROMOTE_IMPLEMENTED_OP, F::Type{<:Number},
                              G::Type{<:Union{SingleDecision, VectorOfDecisions}})
    MOIU.promote_operation(op, F, F, G)
end
function MA.promote_operation(op::PROMOTE_IMPLEMENTED_OP,
                              F::Type{<:Union{SingleDecision, VectorOfDecisions}},
                              G::Type{<:Number})
    MOIU.promote_operation(op, G, F, G)
end

function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, f::AffineDecisionFunction)
    MA.mutable_operate!(op, f.variable_part)
    MA.mutable_operate!(zero, f.decision_part)
    return f
end
function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, f::QuadraticDecisionFunction{T}) where T
    MA.mutable_operate!(op, f.variable_part)
    MA.mutable_operate!(zero, f.decision_part)
    MA.mutable_operate!(zero, f.cross_terms)
    return f
end

function MA.mutable_operate!(::typeof(-), f::AffineDecisionFunction)
    MA.mutable_operate!(op, f.variable_part)
    MA.mutable_operate!(op, f.decision_part)
    return f
end
function MA.mutable_operate!(::typeof(-), f::QuadraticDecisionFunction{T}) where T
    MA.mutable_operate!(op, f.variable_part)
    MA.mutable_operate!(op, f.decision_part)
    MA.mutable_operate!(op, f.cross_terms)
    return f
end

function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::AffineDecisionFunction{T},
                             g::T) where T
    MA.mutable_operate!(op, f.variable_part, g)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::AffineDecisionFunction{T},
                             g::MOI.SingleVariable) where T
    _operate_terms!(op, f.variable_part, AffineDecisionFunction{T}(g).variable_part)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::AffineDecisionFunction{T},
                             g::SingleDecision) where T
    _operate_terms!(op, f.decision_part, AffineDecisionFunction{T}(g).decision_part)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::AffineDecisionFunction{T},
                             g::MOI.ScalarAffineFunction{T}) where T
    _operate_terms!(op, f.variable_part, g)
    f.variable_part.constant = op(f.variable_part.constant, g.constant)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::AffineDecisionFunction{T},
                             g::AffineDecisionFunction{T}) where T
    # Variable part
    _operate_terms!(op, f.variable_part, g.variable_part)
    f.variable_part.constant = op(f.variable_part.constant, g.variable_part.constant)
    # Decision part
    _operate_terms!(op, f.decision_part, g.decision_part)
    f.decision_part.constant = op(f.decision_part.constant, g.decision_part.constant)
    return f
end

function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::T) where T
    MA.mutable_operate!(op, f.variable_part, g)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::MOI.SingleVariable) where T
    _operate_terms!(op, f.variable_part, convert(MOI.ScalarAffineFunction{T}, g))
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::SingleDecision) where T
    _operate_terms!(op, f.decision_part, convert(AffineDecisionFunction{T}, g).decision_part)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::MOI.ScalarAffineFunction{T}) where T
    _operate_terms!(op, f.variable_part, g)
    f.variable_part.constant = op(f.variable_part.constant, g.constant)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::MOI.ScalarQuadraticFunction{T}) where T
    _operate_terms!(op, f.variable_part, g)
    f.variable_part.constant = op(f.variable_part.constant, g.constant)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::AffineDecisionFunction{T}) where T
    # Variable_Part
    _operate_terms!(op, f.variable_part, g.variable_part)
    f.variable_part.constant = op(f.variable_part.constant, g.variable_part.constant)
    # Decision part
    _operate_terms!(op, f.decision_part, g.decision_part)
    f.decision_part.constant = op(f.decision_part.constant, f.decision_part.constant)
    return f
end
function MA.mutable_operate!(op::Union{typeof(+), typeof(-)},
                             f::QuadraticDecisionFunction{T},
                             g::QuadraticDecisionFunction{T}) where T
    # Variable part
    _operate_terms!(op, f.variable_part, g.variable_part)
    f.variable_part.constant = op(f.variable_part.constant, g.variable_part.constant)
    # Decision part
    _operate_terms!(op, f.decision_part, g.decision_part)
    f.decision_part.constant = op(f.decision_part.constant, g.decision_part.constant)
    # Cross terms
    _operate_terms!(op, f.cross_terms, g.cross_terms)
    return f
end

MOIU.similar_type(::Type{<:AffineDecisionFunction}, ::Type{T}) where T = AffineDecisionFunction{T}
MOIU.similar_type(::Type{<:QuadraticDecisionFunction}, ::Type{T}) where T = QuadraticDecisionFunction{T}
MOIU.similar_type(::Type{<:VectorAffineDecisionFunction}, ::Type{T}) where T = VectorAffineDecisionFunction{T}

MA.promote_operation(::typeof(real), T::Type{<:Union{SingleDecision, VectorOfDecisions}}) = T
function MA.promote_operation(::typeof(imag), ::Type{T}, ::Type{F}) where {T, F <: SingleDecision}
    return AffineDecisionFunction{T}
end
function MA.promote_operation(::typeof(imag), ::Type{T}, ::Type{F}) where {T, F <: VectorOfDecisions}
    return VectorAffineDecisionFunction{T}
end
function MA.promote_operation(::typeof(conj), ::Type{T}, ::Type{F}) where {T, F <: Union{SingleDecision, VectorOfDecisions}}
    return T
end
function MA.promote_operation(op::Union{typeof(real), typeof(imag), typeof(conj)}, F::Type{<:TypedDecisionLike{T}}) where T
    return MOIU.similar_type(F, MA.promote_operation(op, T))
end

function MOIU.operate(::typeof(imag), ::Type{T}, ::SingleDecision) where T
    return zero(AffineDecisionFunction{T})
end
function MOIU.operate(::typeof(imag), ::Type{T}, f::VectorOfDecisions) where T
    zero_with_output_dimension(VectorAffineDecisionFunction{T}, MOIU.output_dimension(f))
end
function MOIU.operate(::typeof(imag), ::Type, f::TypedDecisionLike)
    imag(f)
end

Base.real(f::Union{SingleDecision, VectorOfDecisions}) = f
Base.real(f::TypedDecisionLike) = operate_coefficients(real, f)
Base.imag(f::TypedDecisionLike) = operate_coefficients(imag, f)
Base.conj(f::Union{SingleDecision, VectorOfDecisions}) = f
Base.conj(f::TypedDecisionLike) = operate_coefficients(conj, f)

function _operate_terms!(::typeof(+),
                         f::MOI.ScalarAffineFunction,
                         g::MOI.ScalarAffineFunction)
    for term in g.terms
        add_term!(f.terms, term)
    end
    return nothing
end

function _operate_terms!(::typeof(-),
                         f::MOI.ScalarAffineFunction,
                         g::MOI.ScalarAffineFunction)
    for term in g.terms
        remove_term!(f.terms, term)
    end
    return nothing
end

function _operate_terms!(::typeof(+),
                         f::MOI.ScalarQuadraticFunction,
                         g::MOI.ScalarAffineFunction)
    for term in g.terms
        add_term!(f.affine_terms, term)
    end
    return nothing
end

function _operate_terms!(::typeof(-),
                         f::MOI.ScalarQuadraticFunction,
                         g::MOI.ScalarAffineFunction)
    for term in g.terms
        remove_term!(f.affine_terms, term)
    end
    return nothing
end

function _operate_terms!(::typeof(+),
                         f::MOI.ScalarQuadraticFunction,
                         g::MOI.ScalarQuadraticFunction)
    for term in g.affine_terms
        add_term!(f.affine_terms, term)
    end
    for term in g.quadratic_terms
        add_term!(f.quadratic_terms, term)
    end
    return nothing
end

function _operate_terms!(::typeof(-),
                         f::MOI.ScalarQuadraticFunction,
                         g::MOI.ScalarQuadraticFunction)
    for term in g.affine_terms
        remove_term!(f.affine_terms, term)
    end
    for term in g.quadratic_terms
        remove_term!(f.quadratic_terms, term)
    end
    return nothing
end