text:
abstract type AbstractOrbital end

"""
    const MQ = Union{Int,Symbol}

Defines the possible types that may represent the main quantum number. It can either be an
non-negative integer or a `Symbol` value (generally used to label continuum electrons).
"""
const MQ = Union{Int,Symbol}

nisless(an::T, bn::T) where T = an < bn
# Our convention is that symbolic main quantum numbers are always
# greater than numeric ones, such that ks appears after 2p, etc.
nisless(an::I, bn::Symbol) where {I<:Integer} = true
nisless(an::Symbol, bn::I) where {I<:Integer} = false

# * Non-relativistic orbital

"""
    struct Orbital{N <: AtomicLevels.MQ} <: AbstractOrbital

Represents an orbital with just the angular quantum number `ℓ`.

The type parameter has to be such that it can represent a proper main quantum number (i.e. a
subtype of [`AtomicLevels.MQ`](@ref)).
"""
struct Orbital{N<:MQ} <: AbstractOrbital
    n::N
    ℓ::Int
    function Orbital(n::Int, ℓ::Int)
        n ≥ 1 || throw(ArgumentError("Invalid main quantum number $(n)"))
        0 ≤ ℓ && ℓ < n || throw(ArgumentError("Angular quantum number has to be ∈ [0,$(n-1)] when n = $(n)"))
        new{Int}(n, ℓ)
    end
    function Orbital(n::Symbol, ℓ::Int)
        new{Symbol}(n, ℓ)
    end
end

Base.show(io::IO, orb::Orbital{N}) where N =
    write(io, "$(orb.n)$(spectroscopic_label(orb.ℓ))")

"""
    degeneracy(orbital::Orbital)

Returns the degeneracy of `orbital` which is `2(2ℓ+1)`

# Examples

```jldoctest
julia> degeneracy(o"1s")
2

julia> degeneracy(o"2p")
6
```
"""
degeneracy(orb::Orbital) = 2*(2orb.ℓ + 1)

"""
    isless(a::Orbital, b::Orbital)

Compares the orbitals `a` and `b` to decide which one comes before the
other in a configuration.

# Examples

```jldoctest
julia> o"1s" < o"2s"
true

julia> o"1s" < o"2p"
true

julia> o"ks" < o"2p"
false
```
"""
function Base.isless(a::Orbital, b::Orbital)
    nisless(a.n, b.n) && return true
    a.n == b.n && a.ℓ < b.ℓ && return true
    false
end

"""
    parity(orbital::Orbital)

Returns the parity of `orbital`, defined as `(-1)^ℓ`.

# Examples

```jldoctest
julia> parity(o"1s")
even

julia> parity(o"2p")
odd
```

"""
parity(orb::Orbital) = p"odd"^orb.ℓ

"""
    symmetry(orbital::Orbital)

Returns the symmetry for `orbital` which is simply `ℓ`.
"""
symmetry(orb::Orbital) = orb.ℓ

"""
    isbound(::Orbital)

Returns `true` is the main quantum number is an integer, `false`
otherwise.

```jldoctest
julia> isbound(o"1s")
true

julia> isbound(o"ks")
false
```
"""
isbound(::Orbital{Int}) = true
isbound(::Orbital{Symbol}) = false

"""
    mℓrange(orbital::Orbital)

Returns the range of valid values of `mℓ` for `orbital`.

# Examples

```jldoctest
julia> mℓrange(o"2p")
-1:1
```
"""
mℓrange(orb::Orbital) = (-orb.ℓ:orb.ℓ)

"""
    struct SpinOrbital{O<:Orbital} <: AbstractOrbital

Spin orbitals are fully characterized orbitals, i.e. the quantum numbers `n`, `ℓ`, `mℓ` and
`ms` are all specified.
"""
struct SpinOrbital{O<:Orbital} <: AbstractOrbital
    orb::O
    mℓ::Int
    spin::Bool
    function SpinOrbital(orb::O, mℓ::Int, spin::Bool) where {O<:Orbital}
        abs(mℓ) ≤ orb.ℓ ||
            throw(ArgumentError("Magnetic quantum number not in valid range -$(orb.ℓ)..$(orb.ℓ)"))
        new{O}(orb, mℓ, spin)
    end
end
function Base.show(io::IO, so::SpinOrbital)
    show(io, so.orb)
    write(io, to_subscript(so.mℓ))
    write(io, so.spin ? "α" : "β")
end

degeneracy(::SpinOrbital) = 1

Base.isless(a::SpinOrbital, b::SpinOrbital) =
    a.orb < b.orb ||
    a.orb == b.orb && a.mℓ < b.mℓ ||
    a.orb == b.orb && a.mℓ == b.mℓ && a.spin > b.spin # We prefer α < β

parity(so::SpinOrbital) = parity(so.orb)
symmetry(so::SpinOrbital) = (symmetry(so.orb), so.mℓ, so.spin)

isbound(so::SpinOrbital) = isbound(so.orb)

Base.promote_type(::Type{SpinOrbital{O}}, ::Type{SpinOrbital}) where O = SpinOrbital
Base.promote_type(::Type{SpinOrbital}, ::Type{SpinOrbital{O}}) where O = SpinOrbital
Base.promote_type(::Type{SpinOrbital{Orbital{I}}}, ::Type{SpinOrbital{Orbital{Symbol}}}) where {I<:Integer} = SpinOrbital
Base.promote_type(::Type{SpinOrbital{Orbital{Symbol}}}, ::Type{SpinOrbital{Orbital{I}}}) where {I<:Integer} = SpinOrbital

"""
    spin_orbitals(orbital)

Generate all permissible spin-orbitals for a given `orbital`, e.g. 2p
-> 2p ⊗ mℓ = {-1,0,1} ⊗ s = {α,β}

# Examples

```jldoctest
julia> spin_orbitals(o"2p")
6-element Array{SpinOrbital{Orbital{Int64}},1}:
 2p₋₁α
 2p₋₁β
 2p₀α
 2p₀β
 2p₁α
 2p₁β
```

"""
function spin_orbitals(orb::O) where {O<:Orbital}
    map([true,false]) do spin
        map(mℓrange(orb)) do mℓ
            SpinOrbital(orb, mℓ, spin)
        end
    end |> so -> vcat(so...) |> sort
end

# * Relativistic orbital

"""
    kappa_to_ℓ(κ::Integer) :: Integer

Calculate the `ℓ` quantum number corresponding to the `κ` quantum number.

Note: `κ` and `ℓ` values are always integers.
"""
function kappa_to_ℓ(kappa::Integer)
    kappa == zero(kappa) && throw(ArgumentError("κ can not be zero"))
    (kappa < 0) ? -(kappa+1) : kappa
end

"""
    kappa_to_j(κ::Integer) :: HalfInteger

Calculate the `j` quantum number corresponding to the `κ` quantum number.

Note: `κ` is always an integer.
"""
function kappa_to_j(kappa::Integer)
    kappa == zero(kappa) && throw(ArgumentError("κ can not be zero"))
    HalfInteger(2*abs(kappa) - 1, 2)
end

"""
    ℓj_to_kappa(ℓ::Integer, j::Real) :: Integer

Converts a valid `(ℓ, j)` pair to the corresponding `κ` value.

**Note:** there is a one-to-one correspondence between valid `(ℓ,j)` pairs and `κ` values
such that for `j = ℓ ± 1/2`, `κ = ∓(j + 1/2)`.
"""
function ℓj_to_kappa(ℓ::Integer, j::Real)
    assert_orbital_ℓj(ℓ, j)
    (j < ℓ) ? ℓ : -(ℓ + 1)
end

function assert_orbital_ℓj(ℓ::Integer, j::Real)
    j = HalfInteger(j)
    s = hi"1/2"
    (ℓ == j + s) || (ℓ == j - s) ||
        throw(ArgumentError("Invalid (ℓ, j) = $(ℓ), $(j) pair, expected j = ℓ ± 1/2."))
    return
end

struct RelativisticOrbital{N<:MQ} <: AbstractOrbital
    n::N
    κ::Int
    function RelativisticOrbital(n::Integer, κ::Integer)
        n ≥ 1 || throw(ArgumentError("Invalid main quantum number $(n)"))
        κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
        ℓ = kappa_to_ℓ(κ)
        0 ≤ ℓ && ℓ < n || throw(ArgumentError("Angular quantum number has to be ∈ [0,$(n-1)] when n = $(n)"))
        new{Int}(n, κ)
    end
    function RelativisticOrbital(n::Symbol, κ::Integer)
        κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
        new{Symbol}(n, κ)
    end
end
RelativisticOrbital(n::MQ, ℓ::Integer, j::Real) = RelativisticOrbital(n, ℓj_to_kappa(ℓ, j))

function Base.show(io::IO, orb::RelativisticOrbital)
    write(io, "$(orb.n)$(spectroscopic_label(kappa_to_ℓ(orb.κ)))")
    orb.κ > 0 && write(io, "⁻")
end

function flip_j(orb::RelativisticOrbital)
    orb.κ == -1 && return RelativisticOrbital(orb.n, -1) # nothing to flip for s-orbitals
    RelativisticOrbital(orb.n, orb.κ < 0 ? abs(orb.κ) - 1 : -(orb.κ + 1))
end

degeneracy(orb::RelativisticOrbital{N}) where N = 2*abs(orb.κ) # 2j + 1 = 2|κ|

function Base.isless(a::RelativisticOrbital, b::RelativisticOrbital)
    nisless(a.n, b.n) && return true
    aℓ, bℓ = kappa_to_ℓ(a.κ), kappa_to_ℓ(b.κ)
    a.n == b.n && aℓ < bℓ && return true
    a.n == b.n && aℓ == bℓ && abs(a.κ) < abs(b.κ) && return true
    false
end

parity(orb::RelativisticOrbital) = p"odd"^kappa_to_ℓ(orb.κ)
symmetry(orb::RelativisticOrbital) = orb.κ

isbound(::RelativisticOrbital{Int}) = true
isbound(::RelativisticOrbital{Symbol}) = false

# * Orbital construction from strings

parse_orbital_n(m::RegexMatch,i=1) =
    isnumeric(m[i][1]) ? parse(Int, m[i]) : Symbol(m[i])

function parse_orbital_ℓ(m::RegexMatch,i=2)
    ℓs = strip(m[i], ['[',']'])
    if isnumeric(ℓs[1])
        parse(Int, ℓs)
    else
        ℓi = findfirst(ℓs, spectroscopic)
        isnothing(ℓi) && throw(ArgumentError("Invalid spectroscopic label: $(m[i])"))
        first(ℓi) - 1
    end
end

function orbital_from_string(::Type{O}, orb_str::AbstractString) where {O<:AbstractOrbital}
    m = match(r"^([0-9]+|.)([a-z]|\[[0-9]+\])([-]{0,1})$", orb_str)
    m === nothing && throw(ArgumentError("Invalid orbital string: $(orb_str)"))
    n = parse_orbital_n(m)
    ℓ = parse_orbital_ℓ(m)
    if O == RelativisticOrbital
        j = ℓ + (m[3] == "-" ? -1 : 1)*1//2
        O(n, ℓ, j)
    else
        m[3] == "" || throw(ArgumentError("Non-relativistic orbitals cannot have their spins explicitly specified"))
        O(n, ℓ)
    end
end

macro o_str(orb_str)
    orbital_from_string(Orbital, orb_str)
end

macro ro_str(orb_str)
    orbital_from_string(RelativisticOrbital, orb_str)
end

function orbitals_from_string(::Type{O}, orbs_str::AbstractString) where {O<:AbstractOrbital}
    map(split(orbs_str)) do orb_str
        m = match(r"^([0-9]+|.)\[([a-z]|[0-9]+)(-([a-z]|[0-9]+)){0,1}\]$", strip(orb_str))
        m === nothing && throw(ArgumentError("Invalid orbitals string: $(orb_str)"))
        n = parse_orbital_n(m)
        ℓs = map(filter(i -> !isnothing(m[i]), [2,4])) do i
            parse_orbital_ℓ(m, i)
        end
        orbs = if O == RelativisticOrbital
            orbs = map(ℓ -> O(n, ℓ, ℓ-1//2), max(first(ℓs),1):last(ℓs))
            append!(orbs, map(ℓ -> O(n, ℓ, ℓ+1//2), first(ℓs):last(ℓs)))
        else
            map(ℓ -> O(n, ℓ), first(ℓs):last(ℓs))
        end
        sort(orbs)
    end |> o -> vcat(o...) |> sort
end

macro os_str(orbs_str)
    orbitals_from_string(Orbital, orbs_str)
end

macro ros_str(orbs_str)
    orbitals_from_string(RelativisticOrbital, orbs_str)
end

function kappa_from_string(κ_str)
    m = match(r"^([a-z]|\[[0-9]+\])([-]{0,1})$", κ_str)
    m === nothing && throw(ArgumentError("Invalid κ string: $(κ_str)"))
    ℓ = parse_orbital_ℓ(m, 1)
    j = ℓ + (m[2] == "-" ? -1 : 1)//2
    ℓj_to_kappa(ℓ, j)
end

macro κ_str(κ_str)
    kappa_from_string(κ_str)
end

export Orbital, SpinOrbital, RelativisticOrbital, @o_str, @ro_str, @os_str, @ros_str, degeneracy, symmetry, isbound, mℓrange, spin_orbitals, @κ_str


text:
# import Rasters._sourcetype, Rasters._filterkw
    # Rasters.defaultmappedcrs,
    # Rasters.defaultcrs;
import Rasters: FileArray, getraster, layerkeys, RasterDataSource,
    defaultcrs, defaultmappedcrs, _sourcetype,
    haslayers, _open

# `lazy`: If true Lazy DiskArray will be used; if false, all data will be read
# into memory.
function Rasters.Raster(ds::NCD.CFVariable, filename::AbstractString, key=nothing;
    crs=nothing, mappedcrs=nothing, dims=nothing, refdims=(),
    name=Symbol(key isa Nothing ? "" : string(key)),
    metadata=metadata(ds), missingval=missingval(ds), write=false,
    source=_sourcetype(filename), lazy = true)

    crs = defaultcrs(source, crs)
    mappedcrs = defaultmappedcrs(source, mappedcrs)
    dims = dims isa Nothing ? DD.dims(ds, crs, mappedcrs) : dims

    data = lazy ? FileArray(ds, filename; key, write) : ds[:]
    return Raster(data, dims, refdims, name, metadata, missingval)
end

raster = Raster;
brick = RasterStack;

export raster, brick, dims

ext: jl


text:
#"""
#Functions to processes arrays for AbcPmcSteps
#"""

#{A <: UnivariateAbcPmc}
#{A <: MultivariateAbcPmc}
#{A <: SingleMeasureAbcPmc}
#{A <: MultiMeasureAbcPmc

"""
`totalsamples_acceptbws(Array{AbcPmcStep, 1})`

Takes array of AbcPmcStep and returns a tuple `(labels, [samples, acceptbws])` where row `i` corresponds to the `ith` entry of the input Array.  The first column will list the total number of samples drawn in the corresponding step while subsequent columns will list acceptbw information.  In the case of a `MultiMeasureAbcPmc` array labels are given as 'acceptbw1', 'acceptbw2', ...
"""
function totalsamples_acceptbws(x::Array{A, 1}) where A <: SingleMeasureAbcPmc
    T = eltype(x[1].acceptbw)
    totalsamples = zeros(T, length(x))
    acceptbws = zeros(T, length(x))
    for (ii, a) in enumerate(x)
        totalsamples[ii] = sum(T, a.nsampled)
        acceptbws[ii] = a.acceptbw
    end
    return (["samples", "acceptbw"], [cumsum(totalsamples) acceptbws])
end

function totalsamples_acceptbws(x::Array{A, 1}) where A <: MultiMeasureAbcPmc
    T = eltype(x[1].acceptbw)
    totalsamples = zeros(T, length(x))
    acceptbws = zeros(T, length(x), length(x[1].acceptbw))
    for (ii, a) in enumerate(x)
        totalsamples[ii] = sum(T, a.nsampled)
        acceptbws[ii, :] = a.acceptbw
    end
    return (["samples"; map(x -> string("acceptbw", x), 1:size(acceptbws, 2))], [cumsum(totalsamples) acceptbws])
end

"""
`particles_weights(Array{AbcPmcStep, 1})`

Takes array of AbcPmcStep and returns a tuple `(labels, [step, particles, weights])` where row `step` corresponds to the entry of the input Array.  In the case of a `MultivariateAbcPmc` array labels are given as 'particle1', 'particle2', ....  The last column contains the weights associated with the particles.
"""
function particles_weights(x::Array{A, 1}) where A <: UnivariateAbcPmc
    n = length(x)
    m = length(x[1].particles)
    out = Array{promote_type(eltype(x[1].particles), eltype(x[1].weights), Int64)}(n * m, 3)
    out[:, 1] = repeat(1:n, inner = m)
    for ii in 1:n
        rng = range((ii - 1) * m + 1, 1, m)
        out[rng, 2] = x[ii].particles
        out[rng, 3] = x[ii].weights.values
    end
    return (["step", "particle", "weight"], out)
end

function particles_weights(x::Array{A, 1}) where A <: MultivariateAbcPmc
    n = length(x)
    m, k = size(x[1].particles)
    out = Array{promote_type(eltype(x[1].particles), eltype(x[1].weights), Int64)}(n * m, k + 2)
    out[:, 1] = repeat(1:n, inner = m)
    for ii in 1:n
        rng = range((ii - 1) * m + 1, 1, m)
        out[rng, 2:(k + 1)] = x[ii].particles
        out[rng, k + 2] = x[ii].weights.values
    end
    return (["step"; map(x -> string("particle", x), 1:k); "weight"], out)
end

#workspace()
#whos()


text:
function sumdigits(number)
    rest = abs(number)
    sum = 0
    while rest > 0
      sum += rest % 10
      rest = rest ÷ 10
    end
    return sum
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module SLICOTMath
using SLICOT_jll
using LinearAlgebra
using LinearAlgebra: BlasInt
using DocStringExtensions

function chkargsok(ret::BlasInt)
    if ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to SLICOT call"))
    end
end
include("slicot_m.jl")
include("m_docs.jl")

end


text:
@testset "Scitypes" begin
  @test scitype(Composition(a=1, b=2)) === Compositional{2}
  @test scitype(Composition(a=1, b=2, c=3)) === Compositional{3}

  table = (a=rand(100), b=rand(100), c=rand(100), d=1:100)
  coda  = compose(table, (:a,:b,:c))
  @test schema(coda).scitypes == (Count, Compositional{3})
end


text:
using Circuitscape

const FILE_CELLMAP_ASC = joinpath("sample1", "cellmap.asc")
const FILE_POINTS_ASC = joinpath("sample1", "points.asc")

h = Circuitscape.HabitatRasterMap(FILE_CELLMAP_ASC, true)
hab = Circuitscape.Habitat(h)
fp = Circuitscape.FocalPoints(FILE_POINTS_ASC)
Circuitscape.solve_pairwise(hab, fp)


text:
const NOTE_NAMES =
    ("C", "C♯/D♭", "D", "D♯/E♭", "E", "F", "F♯/G♭", "G", "G♯/A♭", "A", "A♯/B♭", "B")

const SHARPS = (
    ("C", "0"),
    ("C", "1"),
    ("D", "0"),
    ("D", "1"),
    ("E", "0"),
    ("F", "0"),
    ("F", "1"),
    ("G", "0"),
    ("G", "1"),
    ("A", "0"),
    ("A", "1"),
    ("B", "0"),
)

const FLATS = (
    ("C", "0"),
    ("D", "-1"),
    ("D", "0"),
    ("E", "-1"),
    ("E", "0"),
    ("F", "0"),
    ("G", "-1"),
    ("G", "0"),
    ("A", "-1"),
    ("A", "0"),
    ("B", "-1"),
    ("B", "0"),
)

const A440_MIDI_CODE = 69
const C0_MIDI_CODE = 12

function get_nearest_midi_code(frequency)
    A440_MIDI_CODE + round(Int, 12 * log(frequency / 440Hz) / log(2))
end

function get_note_parts(midi_code)
    fldmod(round(Int, midi_code) - C0_MIDI_CODE, 12)
end

function round_nearest_midi(frequency)
    get_frequency(get_nearest_midi_code(frequency))
end

function get_frequency(midi_code)
    2.0^((midi_code - A440_MIDI_CODE) / 12) * 440Hz
end

const NOTE_TYPES = (
    "maxima",
    "long",
    "breve",
    "whole",
    "half",
    "quarter",
    "eighth",
    "16th",
    "32nd",
    "64th",
    "128th",
    "256th",
    "512th",
    "1024th",
)

function skip!(measure, direction, beats)
    add_text(new_child(new_child(measure, direction), "duration"), string(beats))
end

function skip!(measure, beats)
    if beats > 0
        skip!(measure, "forward", beats)
    elseif beats < 0
        skip!(measure, "backup", -beats)
    end
end

function new_measure!(part, index)
    measure = new_child(part, "measure")
    set_attribute(measure, "number", string(index))
    measure
end

function maybe_new_measure!(part, measures, index)
    if index == length(measures) + 1
        measure = new_measure!(part, index)
        push!(measures, measure)
        measure
    else
        measures[index]
    end
end

function make_notations!(xml_note)
    new_child(xml_note, "notations")
end

function tie!(xml_note, start_or_stop)
    set_attribute(new_child(xml_note, "tie"), "type", start_or_stop)
    maybe_notations = find_element(xml_note, "notations")
    set_attribute(
        new_child(if maybe_notations === nothing
            new_child(xml_note, "notations")
        else
            maybe_notations
        end, "tied"),
        "type",
        start_or_stop,
    )
end

struct ParsedPitch
    staff::String
    step::String
    alter::String
    octave::String
end

function ParsedPitch(scale, frequency)
    midi_code = get_nearest_midi_code(frequency)
    octave, degree = get_note_parts(midi_code)
    step, alter = scale[degree + 1]
    # middle C
    staff = if midi_code >= 60
        "1"
    else
        "2"
    end
    ParsedPitch(staff, step, alter, string(octave))
end

function add_sub_note!(
    measure,
    beat_type_offset,
    parsed_pitch,
    beat_type_index;
    dotted = false,
)
    xml_note = new_child(measure, "note")
    xml_pitch = new_child(xml_note, "pitch")
    add_text(new_child(xml_pitch, "step"), parsed_pitch.step)
    add_text(new_child(xml_pitch, "alter"), parsed_pitch.alter)
    add_text(new_child(xml_pitch, "octave"), parsed_pitch.octave)
    add_text(
        new_child(xml_note, "duration"),
        string(Int(if dotted
            1.5
        else
            1.0
        end * 2.0^(beat_type_offset - 1))),
    )
    add_text(
        new_child(xml_note, "type"),
        NOTE_TYPES[beat_type_index - (beat_type_offset - 1)],
    )
    if dotted
        new_child(xml_note, "dot")
    end
    add_text(new_child(xml_note, "staff"), parsed_pitch.staff)
    xml_note
end

function add_measure_note!(
    measure,
    beat_type_index,
    parsed_pitch,
    beats;
    end_previous_tie = false,
    start_new_tie = false,
)
    # TODO: avoid allocations
    duration_offsets = map(
        first,
        Iterators.filter(((_, digit),) -> Bool(digit), enumerate(digits(beats; base = 2))),
    )
    if length(duration_offsets) == 2 && duration_offsets[1] + 1 == duration_offsets[2]
        xml_note = add_sub_note!(
            measure,
            duration_offsets[2],
            parsed_pitch,
            beat_type_index;
            dotted = true,
        )
        if end_previous_tie
            tie!(xml_note, "stop")
        end
        if start_new_tie
            tie!(xml_note, "start")
        end
    else
        number_of_sub_notes = length(duration_offsets)
        for (sub_index, beat_type_offset) in enumerate(duration_offsets)
            xml_note =
                add_sub_note!(measure, beat_type_offset, parsed_pitch, beat_type_index)
            if end_previous_tie || sub_index != 1
                tie!(xml_note, "stop")
            end
            if start_new_tie || sub_index != number_of_sub_notes
                tie!(xml_note, "start")
            end
        end
    end
    skip!(measure, -beats)
end

function split_measures(beats_per_measure, chord_from_beat, beats)
    beginning = beats_per_measure - chord_from_beat
    middle_measures, ending = fldmod(beats - beginning, beats_per_measure)
    beginning, middle_measures, ending
end

function add_whole_note!(
    part,
    measures,
    measure_index,
    key,
    scale,
    note,
    beats_per_measure,
    chord_from_beat,
    beat_type_index,
)
    parsed_pitch = ParsedPitch(scale, key * Rational(note.interval))
    note_beats = note.beats
    destination = chord_from_beat + note_beats
    if destination > beats_per_measure
        # TODO: avoid allocations
        beginning, middle_measures, ending =
            split_measures(beats_per_measure, chord_from_beat, note_beats)
        measure_index_beats =
            map((index -> (index, beats_per_measure)), (measure_index + 1):(measure_index + middle_measures))
        if beginning > 0
            pushfirst!(measure_index_beats, (measure_index, beginning))
        end
        if ending > 0
            push!(measure_index_beats, (measure_index + middle_measures + 1, ending))
        end
        number_of_note_measures = length(measure_index_beats)
        for (sub_index, (measure_index, beats)) in enumerate(measure_index_beats)
            add_measure_note!(
                maybe_new_measure!(part, measures, measure_index),
                beat_type_index,
                parsed_pitch,
                beats;
                end_previous_tie = sub_index != 1,
                start_new_tie = sub_index != number_of_note_measures,
            )
        end
    else
        add_measure_note!(
            measures[measure_index],
            beat_type_index,
            parsed_pitch,
            note_beats,
        )
    end
end

function add_chord!(
    part,
    measures,
    key,
    scale,
    chord,
    measure_index,
    beats_per_measure,
    beat_type_index,
    chord_from_beat,
)
    key = round_nearest_midi(Rational(chord.modulation.interval) * key)
    for note in chord.notes
        add_whole_note!(
            part,
            measures,
            measure_index,
            key,
            scale,
            note,
            beats_per_measure,
            chord_from_beat,
            beat_type_index,
        )
    end
    chord_beats = chord.beats
    if chord_beats > 0
        destination = chord_from_beat + chord_beats
        if destination > beats_per_measure
            beginning, middle_measures, ending =
                split_measures(beats_per_measure, chord_from_beat, chord_beats)
            if beginning > 0
                skip!(measures[measure_index], beginning)
            end
            for middle_index in 1:middle_measures
                skip!(maybe_new_measure!(part, measures, measure_index + middle_index), beats_per_measure)
            end
            if ending > 0
                measure_index = measure_index + middle_measures + 1
                skip!(maybe_new_measure!(part, measures, measure_index), ending)
                chord_from_beat = ending
            else
                chord_from_beat = beats_per_measure
                measure_index = measure_index + middle_measures
            end
        else
            skip!(measures[measure_index], chord_beats)
            chord_from_beat = chord_from_beat + chord_beats
        end
    elseif chord_beats < 0
        if chord_beats < -chord_from_beat
            throw(
                ErrorException(
                    "Skipping chord backwards into a previous measure is not unsupported",
                ),
            )
        end
        skip!(measures[measure_index], chord_beats)
        chord_from_beat = chord_from_beat + chord_beats
    end
    key, measure_index, chord_from_beat
end

"""
    make_music_xml(
        song;
        key_fifths = 0,
        beats_per_quarter_note = 1,
        time_signature_numerator = 4,
        time_signature_denominator = 4,
    )

Make a music xml file from a song. You can save it to a new file with
`LightXML.save_file`.

Key fifths is the number of sharps and flats in the key signature. So if it's
`+1`, the key is G major/E minor, and if it's `-1`, the key is F major/D minor.

```jldoctest
julia> using Justly

julia> using LightXML: save_file

julia> cd(joinpath(pkgdir(Justly), "test"))

julia> song = read_justly("song.justly");

julia> xml_song = make_music_xml(song, beats_per_quarter_note = 2);

julia> save_file(xml_song, "test_song.musicxml");

julia> read("test_song.musicxml", String) == read("song.musicxml", String)
true

julia> rm("test_song.musicxml")
```
"""
function make_music_xml(
    song;
    key_fifths = 0,
    beats_per_quarter_note = 1,
    time_signature_numerator = 4,
    time_signature_denominator = 4,
)
    beats_per_measure =
        Int(4 * time_signature_numerator / time_signature_denominator * beats_per_quarter_note)
    scale = if key_fifths >= 0
        SHARPS
    else
        FLATS
    end

    beat_type_index = 6 + Int(log2(beats_per_quarter_note))

    part_id = "P1"
    score = XMLDocument()
    score_partwise = create_root(score, "score-partwise")
    set_attribute(score_partwise, "version", "4.0")
    part_list = new_child(score_partwise, "part-list")
    first_part = new_child(part_list, "score-part")
    set_attribute(first_part, "id", part_id)
    add_text(new_child(first_part, "part-name"), "Part 1")
    part = new_child(score_partwise, "part")
    set_attribute(part, "id", part_id)
    measure_index = 1
    first_measure = new_measure!(part, measure_index)
    measures = [first_measure]
    attributes = new_child(first_measure, "attributes")
    add_text(new_child(attributes, "divisions"), string(beats_per_quarter_note))
    add_text(new_child(new_child(attributes, "key"), "fifths"), string(key_fifths))
    time_signature = new_child(attributes, "time")
    add_text(new_child(time_signature, "beats"), string(time_signature_numerator))
    add_text(new_child(time_signature, "beat-type"), string(time_signature_denominator))
    add_text(new_child(attributes, "staves"), "2")
    treble_clef = new_child(attributes, "clef")
    set_attribute(treble_clef, "number", "1")
    add_text(new_child(treble_clef, "sign"), "G")
    add_text(new_child(treble_clef, "line"), "2")
    bass_clef = new_child(attributes, "clef")
    set_attribute(bass_clef, "number", "2")
    add_text(new_child(bass_clef, "sign"), "F")
    add_text(new_child(bass_clef, "line"), "4")

    chord_from_beat = 0
    key = round_nearest_midi(song.initial_key)
    chord = song.chords[1]
    for chord in song.chords
        key, measure_index, chord_from_beat = add_chord!(
            part,
            measures,
            key,
            scale,
            chord,
            measure_index,
            beats_per_measure,
            beat_type_index,
            chord_from_beat,
        )
    end

    score
end
export make_music_xml

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
# ---
# title: 522. Longest Uncommon Subsequence II
# id: problem522
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/longest-uncommon-subsequence-ii/description/>
# hidden: true
# ---
#
# Given a list of strings, you need to find the longest uncommon subsequence
# among them. The longest uncommon subsequence is defined as the longest
# subsequence of one of these strings and this subsequence should not be **any**
# subsequence of the other strings.
#
# A **subsequence** is a sequence that can be derived from one sequence by
# deleting some characters without changing the order of the remaining elements.
# Trivially, any string is a subsequence of itself and an empty string is a
# subsequence of any string.
#
# The input will be a list of strings, and the output needs to be the length of
# the longest uncommon subsequence. If the longest uncommon subsequence doesn't
# exist, return -1.
#
# **Example 1:**
#
#
#
#     Input: "aba", "cdc", "eae"
#     Output: 3
#
#
# **Note:**
#
#   1. All the given strings' lengths will not exceed 10.
#   2. The length of the given list will be in the range of [2, 50].
#
#
## @lc code=start
using LeetCode

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


text:
#
#   PropagateThreeBody.jl
#
#   Includes functions and structures for propagating orbits
#   within the circular restricted three-body problem.
#

"""
Uses `DifferentialEquations` solvers to propagate normalized,
Synodic CR3BP states Δt into the future. All keyword arguments are passed
directly to `DifferentialEquations` solvers.
"""
propagate(orbit::NormalizedSynodicCR3BPOrbit, Δt::Unitful.Time; kwargs...) = propagate(
    orbit, nondimensionalize(Δt, normalized_time_unit(orbit.system)); kwargs...
)


text:
abstract type AbstractTransformer <: Layer end

# struct PwFFN
#     din::Dense
#     dout::Dense
# end
#
# @treelike PwFFN
#
#
# "just a wrapper for two dense layer."
# PwFFN(size::Int, h::Int, act = relu) = PwFFN(
#     Dense(input=size, output=h, activation=act),
#     Dense(h, size)
# )
#
# function (pw::PwFFN)(x)::AbstractMatrix
#     # size(x) == (dims, seq_len)
#     pw.dout(pw.din(x))
# end

struct Transformer <: AbstractTransformer
    mh::MultiheadAttention
    mhn::LayerNorm
    pw::MLP
    pwn::LayerNorm
    drop::Dropout
end

"""
    Transformer(size::Int, head::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1)
    Transformer(size::Int, head::Int, hs::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1)
Transformer layer.
`size` is the input size. if `hs` is not specify, use `div(size, head)` as the hidden size of multi-head attention.
`ps` is the hidden size & `act` is the activation function of the positionwise feedforward layer.
When `future` is `false`, the k-th token can't see the j-th tokens where j > k. `pdrop` is the dropout rate.
"""
function Transformer(size::Int, head::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    Transformer(size, head, div(size, head), ps;future=future, act=act, pdrop=pdrop)
end

Transformer(size::Int, head::Int, hs::Int, ps::Int; future::Bool = true, act = ReLU(), pdrop = 0.1) = Transformer(
    MultiheadAttention(head, size, hs, size; future=future, pdrop=pdrop),
    LayerNorm(size),
    MLP(size, ps, size; activation=act),
    LayerNorm(size),
    Dropout(pdrop),
)

function (t::Transformer)(x, mask=nothing)
    a = t.mh(x, x, x; mask=mask)
    a = t.drop(a)
    res_a = x .+ a
    if ndims(x) == 3
        insize = size(res_a)
        res_a = reshape(res_a, insize[1], :)
    end
    res_a = t.mhn(res_a)
    pwffn = t.pw(res_a)
    pwffn = t.drop(pwffn)
    res_pwffn = res_a .+ pwffn
    res_pwffn = t.pwn(res_pwffn)
    if ndims(x) == 3
        res_pwffn = reshape(res_pwffn, :, Base.tail(insize)...)
    end
    res_pwffn
end

function Base.show(io::IO, t::Transformer)
    hs = div(size(t.mh.iqproj)[1], t.mh.head)
    h, ps = size(t.pw.layers[end])

    print(io, "Transformer(")
    print(io, "head=$(t.mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h)")
    print(io, ", dropout=$(t.drop.p))")
    print(io, ")")
end

struct TransformerDecoder <: AbstractTransformer
    mh::MultiheadAttention
    mhn::LayerNorm
    imh::MultiheadAttention
    imhn::LayerNorm
    pw::MLP
    pwn::LayerNorm
    drop::Dropout
end

"""
    TransformerDecoder(size::Int, head::Int, ps::Int; act = ReLU(), pdrop = 0.1)
    TransformerDecoder(size::Int, head::Int, hs::Int, ps::Int; act = ReLU(), pdrop = 0.1)
TransformerDecoder layer. Decode the value from a Encoder.
`size` is the input size. if `hs` is not specify, use `div(size, head)` as the hidden size of multi-head attention.
`ps` is the hidden size & `act` is the activation function of the positionwise feedforward layer.
`pdrop` is the dropout rate.
"""
function TransformerDecoder(size::Int, head::Int, ps::Int; act = ReLU(), pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    TransformerDecoder(size, head, div(size, head), ps; act=act, pdrop=pdrop)
end

TransformerDecoder(size::Int, head::Int, hs::Int, ps::Int; act = ReLU(), pdrop = 0.1) = TransformerDecoder(
    MultiheadAttention(head, size, hs, size; future=false, pdrop=pdrop),
    LayerNorm(size),
    MultiheadAttention(head, size, hs, size; future=true, pdrop=pdrop),
    LayerNorm(size),
    MLP(size, ps, size; activation=act),
    LayerNorm(size),
    Dropout(pdrop),
)

function (td::TransformerDecoder)(x, m, mask=nothing)
    a = td.mh(x,x,x)
    a = td.drop(a)
    res_a = x .+ a
    res_a = td.mhn(res_a)

    ia = td.imh(res_a, m, m, mask=mask)
    ia = td.drop(ia)
    res_ia = res_a .+ ia
    if ndims(x) == 3
        insize = size(res_ia)
        res_ia = reshape(res_ia, insize[1], :)
    end
    res_ia = td.imhn(res_ia)
    pwffn = td.pw(res_ia)
    pwffn = td.drop(pwffn)
    res_pwffn = res_ia .+ pwffn
    res_pwffn = td.pwn(res_pwffn)
    if N == 3
        res_pwffn = reshape(res_pwffn, :, Base.tail(insize)...)
    end
    res_pwffn
end

function Base.show(io::IO, td::TransformerDecoder)
    hs = div(size(td.imh.iqproj)[1], td.imh.head)
    h, ps = size(td.pw.layers[end])

    print(io, "TransformerDecoder(")
    print(io, "head=$(td.mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "pwffn_size=$(ps), ")
    print(io, "size=$(h)")
    print(io, ", dropout=$(td.drop.p))")
end

mutable struct PositionEmbedding <: Layer
    trainable::Bool
    embedding::Embed
end
get_value(e::PositionEmbedding, name::Symbol, xs::NamedTuple) = pe(first(xs))

function PE(size, pos, i::Int)
    if rem(i, 2) == 0
        sin(pos/1e4^(i/size))
    else
        cos(pos/1e4^((i-1)/size))
    end
end

function PositionEmbedding(size::Int, max_len::Int = 1024; trainable::Bool = false, atype=arrtype)
    if trainable
        embedding = param(size, max_len; atype=atype, init=randn)
    else
        embedding = Matrix{Float64}(undef, size, max_len)
        for l = 1:max_len
            map!(i->PE(size, l, i), selectdim(embedding, 2, l), 1:size)
        end
        embedding = convert(atype,embedding)
    end
    PositionEmbedding(trainable, Embed(embedding))
end

function (pe::PositionEmbedding)(len::Int)
    max_len = size(pe.embedding)[2]
    if len > max_len
       if pe.trainable
           error("position embedding length exceeded")
       else
           over = similar(weight, D, len)
           copyto!(over, 1, weight, 1, length(weight))
           for l = cur_len+1:len
               over[:,l] .=  convert(atype, map(i->PE(D, l, i), 1:D))
           end
           pe.embedding = Embed(over)
       end
   end
   return pe.embedding.weight[:,1:len]
end

(pe::PositionEmbedding)(x::AbstractArray{Int}) = pe(size(x, 1))
(pe::PositionEmbedding)(x) = x .+ pe(size(x, 2))

function Base.show(io::IO, pe::PositionEmbedding)
    s, max_len = size(pe.embedding)
    if pe.trainable
        print(io, "PositionEmbedding($(s), max_len=$(max_len))")
    else
        print(io, "PositionEmbedding($(s))")
    end
end

struct TransformerModel{E, T <: AbstractTransformer, C} <: Layer
    embed::E
    transformers::T
    classifer::C
end

TransformerModel(embed, transformers) = TransformerModel(embed, transformers, identity)

function Base.show(io::IO, model::TransformerModel)
    print(io, "TransformerModel{")
    print(io, typeof(model.transformers))
    print(io, "}(")
    print(io, model.embed)
    print(io, ", ")
    print(io, model.transformers)
    if model.classifer !== identity
        print(io, ", ")
        print(io, model.classifer)
    end
    print(io, ")")
end


text:
# Julia wrapper for header: /usr/include/scip/scip_expr.h
# Automatically generated using Clang.jl wrap_c

function SCIPgetExprtreeTransformedVars(scip, tree)
    ccall((:SCIPgetExprtreeTransformedVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_EXPRTREE}), scip, tree)
end

function SCIPevalExprtreeSol(scip, tree, sol, val)
    ccall((:SCIPevalExprtreeSol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_EXPRTREE}, Ptr{SCIP_SOL}, Ptr{Cdouble}), scip, tree, sol, val)
end

function SCIPevalExprtreeGlobalBounds(scip, tree, infinity, val)
    ccall((:SCIPevalExprtreeGlobalBounds, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_EXPRTREE}, Cdouble, Ptr{Cint}), scip, tree, infinity, val)
end

function SCIPevalExprtreeLocalBounds(scip, tree, infinity, val)
    ccall((:SCIPevalExprtreeLocalBounds, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_EXPRTREE}, Cdouble, Ptr{Cint}), scip, tree, infinity, val)
end


text:
typealias SCIP_Stage Int8

SCIP_STAGE_INIT         = int8(  0)
SCIP_STAGE_PROBLEM      = int8(  1)
SCIP_STAGE_TRANSFORMING = int8(  2)
SCIP_STAGE_TRANSFORMED  = int8(  3)
SCIP_STAGE_INITPRESOLVE = int8(  4)
SCIP_STAGE_PRESOLVING   = int8(  5)
SCIP_STAGE_EXITPRESOLVE = int8(  6)
SCIP_STAGE_PRESOLVED    = int8(  7)
SCIP_STAGE_INITSOLVE    = int8(  8)
SCIP_STAGE_SOLVING      = int8(  9)
SCIP_STAGE_SOLVED       = int8( 10)
SCIP_STAGE_EXITSOLVE    = int8( 11)
SCIP_STAGE_FREETRANS    = int8( 12)
SCIP_STAGE_FREE         = int8( 13)

SCIP_STAGE = Dict{SCIP_Stage, String}()
SCIP_STAGE[SCIP_STAGE_INIT]         = "SCIP data structures are initialized, no problem exists"
SCIP_STAGE[SCIP_STAGE_PROBLEM]      = "the problem is being created and modified"
SCIP_STAGE[SCIP_STAGE_TRANSFORMING] = "the problem is being transformed into solving data space"
SCIP_STAGE[SCIP_STAGE_TRANSFORMED]  = "the problem was transformed into solving data space"
SCIP_STAGE[SCIP_STAGE_INITPRESOLVE] = "presolving is initialized"
SCIP_STAGE[SCIP_STAGE_PRESOLVING]   = "the problem is being presolved"
SCIP_STAGE[SCIP_STAGE_EXITPRESOLVE] = "presolving is exited"
SCIP_STAGE[SCIP_STAGE_PRESOLVED]    = "the problem was presolved"
SCIP_STAGE[SCIP_STAGE_INITSOLVE]    = "the solving process data is being initialized"
SCIP_STAGE[SCIP_STAGE_SOLVING]      = "the problem is being solved"
SCIP_STAGE[SCIP_STAGE_SOLVED]       = "the problem was solved"
SCIP_STAGE[SCIP_STAGE_EXITSOLVE]    = "the solving process data is being freed"
SCIP_STAGE[SCIP_STAGE_FREETRANS]    = "the transformed problem is being freed"
SCIP_STAGE[SCIP_STAGE_FREE]         = "SCIP data structures are being freed"


text:
using Turing, ContinuousBenchmarks.TuringTools, ContinuousBenchmarks.Reporter

data = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

log_report("Dummy Benchmark started!")

@model constrained_test(obs) = begin
    p ~ Beta(2,2)
    for i = 1:length(obs)
        obs[i] ~ Bernoulli(p)
    end
    p
end

log_report("Dummy model constructed!")

LOG_DATA = @tbenchmark(HMC(1000, 1.5, 3), constrained_test, data)

log_report("Dummy benchmark finished!")

print_log(LOG_DATA)

log_report("Dummy benchmark printed!")


text:
using RipQP, QPSReader, QuadraticModels, Quadmath
using Test

function createQuadraticModel128(qpdata; name="qp_pb")
    return QuadraticModel(convert(Array{Float128}, qpdata.c), qpdata.qrows, qpdata.qcols,
            convert(Array{Float128}, qpdata.qvals),
            Arows=qpdata.arows, Acols=qpdata.acols,
            Avals=convert(Array{Float128}, qpdata.avals),
            lcon=convert(Array{Float128}, qpdata.lcon),
            ucon=convert(Array{Float128}, qpdata.ucon),
            lvar=convert(Array{Float128}, qpdata.lvar),
            uvar=convert(Array{Float128}, qpdata.uvar),
            c0=Float128(qpdata.c0), x0 = zeros(Float128, length(qpdata.c)), name=name)
end

@testset "mono_mode" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1))
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats1.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3))
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "multi_mode" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(mode = :multi), display=false)
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(mode = :multi), display=false)
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(mode = :multi), display=false)
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "dynamic_regularization" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(sp = K2LDLParams(regul=:dynamic)), display=false)
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(sp = K2LDLParams(regul=:dynamic)), display=false)
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(sp = K2LDLParams(regul=:dynamic)), display=false)
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "centrality_corrections" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(kc=-1), display=false) # automatic centrality corrections computation
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(kc=2), display=false)
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(kc=2), display=false)
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "Float128" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    qm128_1 = createQuadraticModel128(qps1)
    stats1 = ripqp(qm128_1, itol = InputTol(ϵ_rb32=0.1, ϵ_rb64=0.01), iconf = InputConfig(mode=:multi, normalize_rtol=false), display=false)
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    qm128_2 = createQuadraticModel128(qps2)
    stats2 = ripqp(qm128_2, display=false)
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable
end

@testset "refinement" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement=:zoom), display=false) # automatic centrality corrections computation
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement=:ref), display=false) # automatic centrality corrections computation
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(mode=:multi, refinement=:multizoom),  display=false)
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(mode=:multi, refinement=:multiref),  display=false)
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "K2_5" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), display=false, iconf = InputConfig(sp = K2_5LDLParams()))
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), display=false, iconf = InputConfig(sp = K2_5LDLParams(), mode = :multi))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps3 = readqps("HS52.SIF") # free bounds
    stats3 = ripqp(QuadraticModel(qps3), display=false, iconf = InputConfig(sp = K2_5LDLParams(regul=:dynamic)))
    @test isapprox(stats3.objective, 5.32664756, atol=1e-2)
    @test stats3.status == :acceptable
end

@testset "IPF" begin
    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), display=false, iconf = InputConfig(solve_method = :IPF, mode=:multi))
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), display=false, iconf = InputConfig(solve_method = :IPF))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), display=false, iconf = InputConfig(solve_method = :IPF, sp = K2_5LDLParams(), refinement=:zoom))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps1 = readqps("QAFIRO.SIF") #lower bounds
    stats1 = ripqp(QuadraticModel(qps1), display=false, iconf = InputConfig(solve_method = :IPF, mode=:multi, refinement=:multiref))
    @test isapprox(stats1.objective, -1.59078179, atol=1e-2)
    @test stats1.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), display=false, iconf = InputConfig(solve_method = :IPF, mode=:multi,
                                                                            refinement=:multizoom, sp = K2_5LDLParams()))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable

    qps2 = readqps("HS21.SIF") # low/upp bounds
    stats2 = ripqp(QuadraticModel(qps2), display=false, iconf = InputConfig(solve_method = :IPF, refinement=:zoom,
                                                                            sp = K2LDLParams(regul=:dynamic)))
    @test isapprox(stats2.objective, -9.99599999e1, atol=1e-2)
    @test stats2.status == :acceptable
end


text:
using POMDPGym, POMDPs, Test
pomdp = PendulumPOMDP()
actions(pomdp)
@test length(actions(pomdp)) == 2

s = rand(initialstate(pomdp))
@test size(s) == (2,)

c = render(pomdp, s, -2)

mdp = InvertedPendulumMDP()

max_forks_repo_name: ancorso/POMDPGym


text:
using Distributed
addprocs(4)
@everywhere include("reproducibility.jl")
##

@everywhere using QSWalk
@everywhere using LightGraphs
@everywhere using NPZ
@everywhere using LinearAlgebra, SparseArrays
data_folder_name = "propagation_qsw"

## ######################################################################
# ######################## prob_dist ####################################
# #######################################################################

n_ex = 121
@assert isodd(n_ex)
t_ex = 15

g_ex = PathGraph(n_ex)
mid_ex = div(n_ex+1, 2)

ϱ = proj(mid_ex, n_ex)
adj = adjacency_matrix(g_ex)

s_global = evolve_generator(adj, [adj])
ϱ_global = evolve(s_global, ϱ, t_ex)
npzwrite("$data_folder_name/global.npz", real.(diag(ϱ_global)))

s_ctqw = evolve_generator(adj, [spzeros(size(adj)...)])
ϱ_ctqw = evolve(s_ctqw, ϱ, t_ex)
npzwrite("$data_folder_name/ctqw.npz", real.(diag(ϱ_ctqw)))

s_local = evolve_generator(adj, local_lind(adj))
ϱ_local = evolve(s_local, ϱ, t_ex)
npzwrite("$data_folder_name/local.npz", real.(diag(ϱ_local)))

## ######################################################################
# ####################### α calculator new ##############################
# #######################################################################
ωs = 0.:0.2:1.
max_n = 300
t_scales = 0.02:0.02:1
t_scales |> collect

##
println("Global")
αs_global = []
for t_scale = t_scales
    t = t_scale*max_n
    println("t: $t")
    n = ceil(Int, t/0.2)
    n += iseven(n)
    mid = ceil(Int, div(n+1, 2))

    g = PathGraph(n)
    h = adjacency_matrix(g)
    s = pmap(ω->evolve_generator(h, [h], ω), ωs)
    initial_state = proj(mid, n)
    println("n: $n")
    @time prob_dists = pmap(x-> diag(evolve(x, initial_state, t)), s)

    positions = -(mid-1):(mid-1)
    sec_moments = pmap(pdist -> sum(pdist .* positions.^2), prob_dists)
    push!(αs_global, log.(sec_moments))
end
npzwrite("$data_folder_name/sec_moments_global.npz", hcat(αs_global...))
##
println("Local")

αs_global = []
for t_scale = t_scales
    t = t_scale*max_n
    println("t: $t")
    n = ceil(Int, t/0.2)
    n += iseven(n)
    mid = ceil(Int, div(n+1, 2))

    g = PathGraph(n)
    h = adjacency_matrix(g)
    s = pmap(ω->evolve_generator(h, local_lind(h), ω), ωs)
    initial_state = proj(mid, n)
    println("n: $n")
    @time prob_dists = pmap(x-> diag(evolve(x, initial_state, t)), s)

    positions = -(mid-1):(mid-1)
    sec_moments = pmap(pdist -> sum(pdist .* positions.^2), prob_dists)
    push!(αs_global, log.(sec_moments))
end
npzwrite("$data_folder_name/sec_moments_local.npz", hcat(αs_global...))


text:
"""
Provides a basic interface for defining and solving MDPs/POMDPs
"""
module POMDPs

using Random
import Base: rand
import Statistics
import Distributions: rand, pdf, mode, mean, support
import Pkg
import LibGit2

export
    # Abstract type
    POMDP,
    MDP,

    # Discrete Functions
    n_states,
    n_actions,
    n_observations,

    # Model functions
    discount,
    states,
    actions,
    observations,
    transition,
    observation,
    reward,
    isterminal,
    isterminal_act,
    isterminal_obs,

    # Generative model functions
    generate_s,
    generate_o,
    generate_sr,
    generate_so,
    generate_or,
    generate_sor,
    initialstate,

    # Discrete Functions
    length,
    stateindex,
    actionindex,
    obsindex,
    weight,

    # Common Functions
    rand,
    pdf,
    mode,
    mean,
    dimensions,
    lowerbound,
    upperbound,
    getindex,
    support,
    sampletype,
    initialstate_distribution,

    max_possible_weight,
    new_particle,
    reset_distribution,

    # Solver types
    Solver,
    solve,

    # Beliefs
    Updater,
    update,
    initialize_belief,

    # Policy
    Policy,
    action,
    updater,
    value,

    # Simulation
    Simulator,
    simulate,

    # Utilities
    implemented,
    @implemented,
    convert_s,
    convert_a,
    convert_o,
    statetype,
    actiontype,
    obstype,

    # Requirements checking
    RequirementSet,
    check_requirements,
    show_requirements,
    get_requirements,
    requirements_info,
    @POMDP_require,
    @POMDP_requirements,
    @requirements_info,
    @get_requirements,
    @show_requirements,
    @warn_requirements,
    @req,
    @subreq,

    # Deprecated
    state_index,
    action_index,
    obs_index,
    state_type,
    action_type,
    obs_type,
    initial_state,
    initial_state_distribution,
    iterator

include("requirements_internals.jl")
include("requirements_printing.jl")
include("pomdp.jl")
include("solver.jl")
include("simulator.jl")
include("requirements_interface.jl")
include("distribution.jl")
include("belief.jl")
include("space.jl")
include("policy.jl")
include("generative.jl")
include("generative_impl.jl")
include("type_inferrence.jl")
include("constants.jl")
include("utils.jl")

end


text:
"""
Implementation of `Paxos.Transports.Common.Transport` for in-memory communication
within the current process.
"""
module Memory

export memory

using UUIDs
using Base

using ..Common

struct MemoryTransport <: Transport
    listeners::Dict{Any,Channel}
end

struct MemoryConnection <: Connection
    sent::Channel
    received::Channel
end

"""
Create a new transport for passing messages within
"""
function memory()
    MemoryTransport(Dict())
end

# -----------------

function Common.connectTo(transport::MemoryTransport, recipient)
    backlog = get!(transport.listeners,recipient, Channel())
    sent = Channel()
    received = Channel()
    # put the channels in the reverse order of thow the client sees them,
    # because when the serer sends it just writes to the received of the client,
    # and vice versa for client sending to server
    put!(backlog,MemoryConnection(received, sent))
    MemoryConnection(sent, received)
end

function Common.listenOn(handler::Function, transport::MemoryTransport, address)
    @sync try
        @debug "Beginning to listen for memory connections on $address"
        connections = get!(transport.listeners,address, Channel())
        while true
            client = take!(connections)
            @async handler(client)
        end
    finally
        delete!(transport.listeners, address)
        @debug "Finished listening for memory connections on $address"
    end
end

function Common.sendTo(connection::MemoryConnection, message)
    put!(connection.sent, message)
end

function Common.receiveFrom(connection::MemoryConnection)
    take!(connection.received)
end

function Base.close(connection::MemoryConnection)
    close(connection.sent)
    close(connection.received)
end

end


text:
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random

function R_sir_rej!(rate,xc,xd,parms,t,issum::Bool)
	(S,I,R,~) = xd
	(beta,mu) = parms
	infection = beta*S*I
	recovery = mu*I
	rate_display = parms[1]
	if issum == false
		rate[1] = infection
		rate[2] = recovery
		rate[3] = rate_display
		return 0., rate_display + 3.5
	else
		return infection+recovery + rate_display, rate_display + 3.5
	end
end

xc0 = [0.0]
xd0 = [99,10,0,0]
nu = [[-1 1 0 0];[0 -1 1 0];[0 0 0 1]]
parms = [0.1/100.0,0.01]
tf = 150.0

Random.seed!(1234)
	println("--> rejection algorithm for SSA")
	problem = PDMP.PDMPProblem(PDMP.F_dummy,R_sir_rej!,nu, xc0, xd0, parms, (0.0, tf))
	result = PDMP.solve(problem, Rejection(Tsit5()); n_jumps = 1000)

# using Plots
# gr()
# plot(result.time,result.xd[1,:])
#	 plot!(result.time,result.xd[2,:])
#	 plot!(result.time,result.xd[3,:])
#	 plot!(result_chv.time,result_chv.xd[1,:],marker=:d,color=:blue)
#	 plot!(result_chv.time,result_chv.xd[2,:],marker=:d,color=:red)
#	 plot!(result_chv.time,result_chv.xd[3,:],marker=:d,color=:green,line=:step)


text:
include("system.jl")

@doc raw"""
    KineticPropagator(Δt)

Creates a propagator that updates positions according to the canonical equation

``\boldsymbol q\leftarrow\boldsymbol q+\frac{\partial H}{\partial\boldsymbol p}\Delta t``
"""
struct KineticPropagator <: Propagator
    Δt::Real
end
(propagator!::KineticPropagator)(phase::Phase, system::ClassicalSystem) = begin
    phase.q .+= gradient(p -> kinetic(system, p), phase.p) * propagator!.Δt
end

@doc raw"""
    KineticPropagator(Δt)

Creates a propagator that updates momenta according to the canonical equation

``\boldsymbol p\leftarrow\boldsymbol p-\frac{\partial H}{\partial\boldsymbol q}\Delta t``
"""
struct PotentialPropagator <: Propagator
    Δt::Real
end
(propagator!::PotentialPropagator)(phase::Phase, system::ClassicalSystem) = begin
    phase.p .-= gradient(q -> potential(system, q), phase.q) * propagator!.Δt
end

@doc raw"""
    ExternalKineticPropagator(Δt)

Creates a propagator that updates momenta according to the canonical equation, with PIMD external potential only

``\boldsymbol p\leftarrow\boldsymbol p-\frac{\partial V_{\rm ext}}{\partial\boldsymbol q}\Delta t``
"""
struct ExternalPotentialPropagator <: Propagator
    Δt::Real
end
(propagator!::ExternalPotentialPropagator)(phase::Phase, system::ClassicalSystem) = begin
    phase.p .-= gradient(q -> potentialExternal(system, q), phase.q) * propagator!.Δt
end

@doc raw"""
    LangevinPropagator(system, Δt)

Creates a propagator that thermostats momenta using Langevin equation

``\boldsymbol p\leftarrow e^{-\Gamma\Delta t}\boldsymbol p+\sqrt{1-e^{-2\Gamma\Delta t}}\sqrt{\frac{M}{\beta}}\boldsymbol{\xi}``

where ``\boldsymbol{\xi}`` is sampled from standard normal distribution and ``\Gamma`` is determined by the characteristic frequency.
"""
struct LangevinPropagator <: Propagator
    dissipation::Operator
    fluctuation::Operator
    LangevinPropagator(system::ClassicalSystem, Δt::Real) = begin
        Γ = friction(system)
        dissipation = exp(-Δt * Γ)
        fluctuation = sqrt(I - dissipation^2) * sqrt(system.M / system.β)
        new(dissipation, fluctuation)
    end
end
(propagator!::LangevinPropagator)(phase::Phase, system::ClassicalSystem) = begin
    phase.p .= propagator!.dissipation * phase.p + propagator!.fluctuation * randn(system.n)
end

friction(system::Primitive) = √(inv(system.M) * system.A)
friction(system::Union{Normal, Staging}) = Diagonal(√system.n / system.β * I, system.n)

@doc raw"""
    GeneralizedCayleyPropagator(system, generalizedCayleyFunction, Δt)

Creates a generalized Cayley propagator that  PIMD internal potential using a generalized cayley function ``\varphi``

``\boldsymbol z\leftarrow\varphi(\mathsf J\mathsf H)\boldsymbol z``

Examples of generalized Cayley functions are `exp()`, `cayley()` and `sqrtcayley()` (recommended).
"""
struct GeneralizedCayleyPropagator <: Propagator
    matrixRepresentation::Matrix
    GeneralizedCayleyPropagator(system::IsomorphicClassicalSystem, generalizedCayleyFunction::Function, Δt::Real) = begin
        n = system.n
        B = Matrix([0I inv(system.M); -system.A 0I]) * Δt
        new(generalizedCayleyFunction(B))
    end
end
(propagator!::GeneralizedCayleyPropagator)(phase::Phase, system::ClassicalSystem) = begin
    z = vcat(phase.q, phase.p)
    z .= propagator!.matrixRepresentation * z
    phase.q .= z[1:system.n]
    phase.p .= z[system.n+1:2(system.n)]
end

const cayley = A::Matrix -> (I + A / 2) / (I - A / 2)
const sqrtcayley = A::Matrix -> (real ∘ √)((I + A)/(I - A))


text:
module LDAPClient
using OpenLDAPClient_jll

include("enums.jl")
include("types.jl")
include("libldap.jl")
include("liblber.jl")
include("messages.jl")
include("attributes.jl")

function is_ldap_url(url::AbstractString) :: Bool
    ldap_is_ldap_url(url) != 0
end

function get_error_message(error_code::Integer)
    return unsafe_string(ldap_err2string(error_code))
end

function get_error_message(err::AuthErr)
    return get_error_message(err.err_code)
end

function error_check(err::Integer)
    if err != 0
        error("LDAP Error Code $err: $(get_error_message(err))")
    end
end

function LDAPConnection(uri::AbstractString; protocol::LDAPVersion=LDAP_VERSION3)
    @assert is_ldap_url(uri) "$uri is not a valid LDAP URL."
    ldp_handle_ref = Ref{Ptr{Cvoid}}()
    err = ldap_initialize(ldp_handle_ref, uri)
    error_check(err)
    result = LDAPConnection(ldp_handle_ref[], uri)
    set_protocol_version(result, protocol)
    return result
end

function set_protocol_version(ldap::LDAPConnection, protocol::LDAPVersion)
    err = ldap_set_option(ldap.handle, LDAP_OPT_PROTOCOL_VERSION, protocol)
    error_check(err)
end

function get_protocol_version(ldap::LDAPConnection)
    out_int = Ref{Cint}(0)
    err = ldap_get_option(ldap.handle, LDAP_OPT_PROTOCOL_VERSION, out_int)
    return LDAPVersion(out_int[])
end

@inline parse_opt_string(str::Cstring) = str == C_NULL ? nothing : unsafe_string(str)

@inline function parse_null_terminated_string_vector(ptr::Ptr{Ptr{UInt8}}) :: Vector{String}
    result = Vector{String}()

    if ptr == C_NULL
        return result
    end

    local i = 1
    while true
        cstr_ptr = unsafe_load(ptr, i)
        if cstr_ptr == C_NULL
            break
        else
            push!(result, unsafe_string(cstr_ptr))
        end

        i += 1
    end

    return result
end

function URL(url::AbstractString)

    ldap_url_desc_handle_ref = Ref{Ptr{LDAPURLDesc}}()
    err = ldap_url_parse(url, ldap_url_desc_handle_ref)
    error_check(err)
    handle = ldap_url_desc_handle_ref[]

    local result::URL

    try
        ldap_url_desc = unsafe_load(ldap_url_desc_handle_ref[])

        lud_scheme = unsafe_string(ldap_url_desc.lud_scheme)
        lud_host = unsafe_string(ldap_url_desc.lud_host)
        lud_port = ldap_url_desc.lud_port
        lud_dn = parse_opt_string(ldap_url_desc.lud_dn)
        lud_attrs = parse_null_terminated_string_vector(ldap_url_desc.lud_attrs)
        lud_scope = ldap_url_desc.lud_scope
        lud_filter = parse_opt_string(ldap_url_desc.lud_filter)
        lud_exts = parse_null_terminated_string_vector(ldap_url_desc.lud_exts)
        lud_crit_exts = ldap_url_desc.lud_crit_exts

        result = URL(lud_scheme, lud_host, lud_port, lud_dn, lud_attrs, lud_scope, lud_filter, lud_exts, lud_crit_exts)

    finally
        ldap_free_urldesc(handle)
    end

    return result
end

function simple_bind(ldap::LDAPConnection, who::AbstractString, password::AbstractString)
    err = ldap_simple_bind_s(ldap.handle, who, password)
    error_check(err)
    nothing
end

function unbind(ldap::LDAPConnection)
    err = ldap_unbind_s(ldap.handle)
    error_check(err)
end

function authenticate(uri::AbstractString, who::AbstractString, password::AbstractString;
        protocol::LDAPVersion=LDAP_VERSION3) :: AuthenticationResult

    ldap = LDAPConnection(uri, protocol=protocol)
    err = ldap_simple_bind_s(ldap.handle, who, password)
    result = new_authentication_result(uri, who, err)
    unbind(ldap)
    return result
end

function new_authentication_result(uri::AbstractString, who::AbstractString, err::Integer)
    if err == 0
        return AuthOk(uri, who)
    else
        return AuthErr(uri, who, err)
    end
end

end # module


text:
homodyne(mag, phase; kwargs...) = homodyne(mag .* exp.(1im .* phase); kwargs...)
homodyne(I; kwargs...) = homodyne!(copy(I); kwargs...)
homodyne!(I; dims=1:2, σ=8 .* ones(length(dims))) = I ./= gaussiansmooth3d(I, σ; dims)

function NumART2star(image, TE)
    neco = length(TE)
    t2star(m) = (TE[end]-TE[1]) / 2(neco-1) * (m[1]+m[end]+2sum(m[2:end-1])) / (m[1]-m[end])
    return [t2star(image[I,:]) for I in CartesianIndices(size(image)[1:3])]
end

r2s_from_t2s(t2s) = 1000 ./ t2s


text:
using Statistics:quantile, mean, median, std
using SpecialFunctions:erf

const _shared_params = """    scores_train::AbstractVector{<:Real}
A vector of training scores, typically the result of [`fit`](@ref) with a detector.

    scores_test::AbstractVector{<:Real}
A vector of test scores, typically the result of [`score`](@ref) using a previously fitted detector."""

"""
    combine(scores_mat,
            strategy = :mean)

Combination method to merge outlier scores from multiple detectors by using a combination strateg. This function is
typically used within a [`Evaluator`](@ref). TODO: Add AOM/MOA/LSCP combination strategies

Parameters
----------
    scores_mat::AbstractMatrix{T}
A matrix, with each row representing the scores for a specific instance and each column representing a detector.

    strategy::Symbol=:mean
Determines how to combine the scores of multiple detectors, e.g. maximum or mean of all scores.

Returns
----------
    combined_scores::AbstractVector{T}
The combined scores.e.g. the maximum of scores from different detectors.

Examples
----------
    scores = [1 2; 3 4; 5 6]
    combine(scores) # [1.5, 3.5, 5.5]
"""
function combine(scores_mat::AbstractMatrix{<:Real}, strategy::Symbol=:mean)::Scores
    @assert strategy in (:mean, :maximum, :median)

    if strategy == :mean
        return  dropdims(mean(scores_mat, dims=2), dims=2)
    elseif strategy == :maximum
        return dropdims(maximum(scores_mat, dims=2), dims=2)
    elseif strategy == :median
        return dropdims(median(scores_mat, dims=2), dims=2)
    end
end

"""
    normalize(scoresTrain,
              scoresTest)

Transform an array of scores into a range between [0,1] using min-max scaling.

Parameters
----------
$_shared_params

Returns
----------
    normalized_scores::Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}}
The normalized train and test scores.

Examples
----------
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    normalize(scores_train, scores_test) # ([0.0, 0.5, 1.0], [1.0, 1.0, 0.5, 0.0, 0.0])
    normalize(scores_train) # [0.0, 0.5, 1.0]
"""
function normalize(scores_train::Scores, scores_test::Scores)::Result
    minTrain = minimum(scores_train)
    maxTrain = maximum(scores_train)
    @assert minTrain < maxTrain "Cannot normalize scores if they are all equal"
    f = scores -> clamp.((scores .- minTrain) ./ (maxTrain - minTrain), 0, 1)
    f(scores_train), f(scores_test)
end
normalize(scores) = normalize(scores, scores)[1]  # only extract scores train because they are the same

"""
    unify(scores_train,
          scores_test)

Transform an array of scores into a range between [0,1] using unifying scores as described in [1].

Parameters
----------
$_shared_params

Returns
----------
    unified_scores::Tuple{AbstractVector{<:Real}, AbstractVector{<:Real}}
The unified train and test scores.

Examples
----------
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    unify(scores_train, scores_test) # ([0.0, 0.0, 0.68..], [0.95.., 0.68.., 0.0, 0.0, 0.0])
    unify(scores_train) # [0.0, 0.0, 0.68..]

References
----------
Kriegel, Hans-Peter; Kroger, Peer; Schubert, Erich; Zimek, Arthur (2011): Interpreting and Unifying Outlier Scores.
"""
function unify(scores_train::Scores, scores_test::Scores)::Result
    μ = mean(scores_train)
    σ = std(scores_train)
    @assert σ > 0 "Cannot normalize scores if they are all equal"
    f = scores -> clamp.(erf.((scores .- μ) ./ (σ * √2)), 0, 1)
    f(scores_train), f(scores_test)
end
unify(scores) = unify(scores, scores)[1] # only extract scores train because they are the same

"""
    classify(outlier_fraction,
             scores_train,
             scores_test)

Convert an array of scores to an array of classes with `1` indicating normal data and `-1` indicating outliers. The
conversion is based on percentiles of the training data, i.e. all datapoints above the `1 - outlier_fraction` percentile
are considered outliers.

Parameters
----------
    outlier_fraction::Real
The fraction of outliers (number between 0 and 1) in the data used to determine the score threshold to classify the
samples into inliers and outliers.

$_shared_params

Returns
----------
    classes::AbstractVector{<:Integer}
The vector of classes consisting of `-1` (outlier) and `1` (inlier) elements.

Examples
----------
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    classify(0.3, scores_train, scores_test) # [-1, -1, 1, 1, 1]
    classify(0.3, scores_train) # [1, 1, -1]
"""
function classify(outlier_fraction::Real, scores_train::Scores, scores_test::Scores)::Labels
    @assert 0 < outlier_fraction < 1
    ifelse.(scores_test .> quantile(scores_train, 1 - outlier_fraction), CLASS_OUTLIER, CLASS_NORMAL)
end
classify(outlier_fraction, scores) = classify(outlier_fraction, scores, scores)


text:
using Base.Graphics
using Cairo
using Tk

using Lasem
using LaTeXStrings

function view_svg(filename)
    # use lasem to open svg file
    doc = Lasem.lsm_dom_document_new_from_path(filename)
    doc != C_NULL || error("Opening and/or processing of file $filename failed.")

    view = Lasem.lsm_dom_document_create_view(doc)
    view != C_NULL || error("Creating lasem view failed.")

    Lasem.lsm_dom_view_set_resolution(view, 96.0)
    w,h = Lasem.lsm_dom_view_get_size_pixels(view)

    # create Cairo surface and render image into it
    surface = CairoARGBSurface(w, h)
    Lasem.lsm_dom_view_render(view, CairoContext(surface))

    # create Window to display image
    win = Toplevel("SVG", w, h)
    c = Canvas(win)
    pack(c, expand=true, fill="both")

    ctx = getgc(c)

    image(ctx, surface, 0, 0, int(w), int(h))

    reveal(c)
    Tk.update()

    Lasem.destroy(view)
    Lasem.destroy(doc)
end

function view_tex(expr::ASCIIString)
    const margin = 10

    # use lasem to render equation
    doc = Lasem.lsm_mathml_document_new_from_itex(expr)
    doc != C_NULL || error("Processing of $expr failed.")

    view = Lasem.lsm_dom_document_create_view(doc)
    view != C_NULL || error("Creating lasem view failed.")

    Lasem.lsm_dom_view_set_resolution(view, 96.0)
    w,h = Lasem.lsm_dom_view_get_size_pixels(view)

    # create Cairo surface and render image into it
    surface = CairoARGBSurface(w, h)
    Lasem.lsm_dom_view_render(view, CairoContext(surface))

    # create Window to display image
    win = Toplevel("iTeX", w+2*margin, h+2*margin)
    c = Canvas(win)
    pack(c, expand=true, fill="both")

    ctx = getgc(c)

    rectangle(ctx, margin, margin, int(w), int(h))
    set_source_rgb(ctx, 1, 0, 1)
    stroke(ctx)
    image(ctx, surface, margin, margin, int(w), int(h))

    reveal(c)
    Tk.update()

    Lasem.destroy(view)
    Lasem.destroy(doc)
end

view_tex(ls::LaTeXString) = view_tex(bytestring(ls))


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear

Random.seed!(100)

dts = 1 .//2 .^(8:-1:4)
testTol = 0.25

f = (u,p,t)->cos(t)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = sin(t)
prob_ode_sin = ODEProblem(f, 0.,(0.0,1.0))

f = (du,u,p,t)->du[1]=cos(t)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = [sin(t)]
prob_ode_sin_inplace = ODEProblem(f, [0.], (0.0,1.0))

f = (u,p,t)->sin(u)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = 2*acot(exp(-t)*cot(0.5))
prob_ode_nonlinear = ODEProblem(f, 1.,(0.,0.5))

f = (du,u,p,t)->du[1]=sin(u[1])
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = [2*acot(exp(-t)*cot(0.5))]
prob_ode_nonlinear_inplace = ODEProblem(f,[1.],(0.,0.5))

test_problems_only_time = [prob_ode_sin, prob_ode_sin_inplace]
test_problems_linear = [prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear]
test_problems_nonlinear = [prob_ode_nonlinear, prob_ode_nonlinear_inplace]

f_ssp = (u,p,t) -> begin
  sin(10t) * u * (1-u)
end
test_problem_ssp = ODEProblem(f_ssp, 0.1, (0., 8.))
test_problem_ssp_long = ODEProblem(f_ssp, 0.1, (0., 1.e3))

f_ssp_inplace = (du,u,p,t) -> begin
  @. du = sin(10t) * u * (1-u)
end
test_problem_ssp_inplace = ODEProblem(f_ssp_inplace, rand(3,3), (0., 8.))

# test SSP coefficient for explicit Euler
alg = Euler()
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg)+1.e-3, dense=false)
@test any(sol.u .< 0)

alg = SSPRK22()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)
# test SSP property of dense output
sol = solve(test_problem_ssp, alg, dt=1.)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)
sol = solve(test_problem_ssp_inplace, alg, dt=1.)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)

alg = SSPRK33()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # This corresponds to Simpson's rule; due to symmetric quadrature nodes,
  # it is of degree 4 instead of 3, as would be expected.
  @test abs(sim.𝒪est[:final]-1-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)
# test SSP property of dense output
sol = solve(test_problem_ssp, alg, dt=1.)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)
sol = solve(test_problem_ssp_inplace, alg, dt=1.)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)

alg = SSPRK53()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)

alg = SSPRK63()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)

alg = SSPRK73()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)

alg = SSPRK83()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)

alg = SSPRK432()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # higher order as pure quadrature
  @test abs(sim.𝒪est[:final]-1-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)
# test SSP property of dense output
sol = solve(test_problem_ssp, alg, dt=8/5, adaptive=false)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)
sol = solve(test_problem_ssp_inplace, alg, dt=8/5, adaptive=false)
@test mapreduce(t->all(0 .<= sol(t) .<= 1), (u,v)->u&&v, range(0, stop=8, length=50), init=true)

alg = SSPRK932()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false,maxiters=1e7)
@test all(sol.u .>= 0)

alg = SSPRK54()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # convergence order seems to be worse for this problem
  @test abs(sim.𝒪est[:final]+0.25-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  # convergence order seems to be better for this problem
  @test abs(sim.𝒪est[:final]-0.5-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)

alg = SSPRK104()
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test abs(sim.𝒪est[:final]-OrdinaryDiffEq.alg_order(alg)) < testTol
end
# test SSP coefficient
sol = solve(test_problem_ssp_long, alg, dt=OrdinaryDiffEq.ssp_coefficient(alg), dense=false)
@test all(sol.u .>= 0)


text:
################################################################################
#
#   The unified Module for all different sorts of lattices
#   and lattice based calculations.
#
#   Module structure (mostly included in subfiles):
#
#   1) TYPE DEFINITION
#       - Definition of UNITCELL TYPE
#       - Definition of LATTICE TYPE
#       - saving / loading types
#       - information on types
#       - connectivity of tpyes
#
#   2) PRE-IMPLEMENTED UNITCELLS (various 2D and 3D stuff)
#
#   3) CONSTRUCTION OF UNITCELLS FROM SITES / CONSTRUCTION HELPER
#
#   4) CONSTRUCTION OF LATTICES
#       - Building Lattices as (periodic/open/semiperiodic) repeated patterns of Unitcells
#       - Building Lattices by bond distance to an origin site
#       - Building Lattices in a shape around an origin site
#
#   5) MODIFICATION OF LATTICES
#       - Adding / Removing Connections
#       - Adding / Removing sites
#       - Interaction Strengths
#       - Scaling / Rotating / Shifting in real space
#       - Lattice vector optimization
#
#   6) PLAQUETTE RELATED FUNCTIONS
#       - obtaining plaquettes of lattice
#       - printing plaquette statistics
#
#   7) INTERACTION MATRICES (REAL AND MOMENTUM SPACE)
#
#   8) PATHS (IN MOMENTUM SPACE)
#       - type definition
#       - modification of paths
#       - pre-implemented paths
#       - TODO create paths based on unitcell
#
#   9) BRILLOUIN ZONES
#       - type definition
#       - pre-implemented BZs
#       - construction of default BZs (2D)
#       - construction of default BZs (3D)
#       - plotting of BZ
#
#  10) SVG PLOTTING
#       - 2D lattices
#       - TODO 3D lattices
#       - 2D plaquettes
#
#  11) BLENDER PLOTTING
#
#  12) BANDSTRUCTURES (CALCULATION AND PLOTTING)
#
#  13) FERMI SURFACES (CALCULATION AND PLOTTING)
#       - 2D unitcells
#       - 3D unitcells
#
#  14) LUTTINGER TISZA (CALCULATION AND PLOTTING)
#       - LTBandstructure type
#       - Spin interaction matrices
#       - calculation and plotting of LT Bandstructures (along path in k space)
#       - calculation and plotting of LT Groundstates (in k space)
#
################################################################################

# Start of module
module LatticePhysics

################################################################################
#
#   IMPORTED LIBRARIES THAT ARE USED IN LATTICEPHYSICS.JL
#
################################################################################

# JLD is used for all save / load applications
using JLD

# PyPlot is used for all plotting of band structures (plotting of lattices is done via SVG file creation)
using PyCall
using PyPlot

# Optim is used for minimizing the band structure to find the ground state energy of the system
using Optim

################################################################################
#
#   DEFAULT FILENAMES FOR FOLDERS AND OBJECTS
#
################################################################################

# THE FOLDER FOR UNITCELLS
FOLDER_UNITCELLS = "unitcells/"
# THE FOLDER FOR LATICES
FOLDER_LATTICES = "lattices/"

# export the folders
export FOLDER_UNITCELLS
export FOLDER_LATTICES

# Dummy filenames
UNITCELL_DUMMY_FILENAME = "THIS IS NOT A UNIT CELL"
LATTICE_DUMMY_FILENAME = "THIS IS NOT A LATTICE"
export UNITCELL_DUMMY_FILENAME
export LATTICE_DUMMY_FILENAME

# THE FOLDER FOR BAND STRUCTURES
FOLDER_SPECTRA = "band_structures/"

# export the folder
export FOLDER_SPECTRA

################################################################################
#
#   FUNCTIONS TO BUILD DEFAULT FOLDERS
#
################################################################################
function buildFolderUnitcells()
    mkpath(FOLDER_UNITCELLS)
end
export buildFolderUnitcells
function buildFolderLattices()
    mkpath(FOLDER_LATTICES)
end
export buildFolderLattices
function buildFolderSpectra()
    mkpath(FOLDER_SPECTRA)
end
export buildFolderSpectra

function buildFolders()
    buildFolderUnitcells()
    buildFolderLattices()
    buildFolderSpectra()
end
export buildFolders

################################################################################
#
#   TYPE DEFINITIONS OF OBJECT CLASSES IN JULIA
#
#   STRUCTURE OF THE FILE:
#
#   1) Definition of UNITCELL OBJECT
#       - Saving
#       - Loading
#
#   2) Definition of LATTICE OBJECT
#       - Saving
#       - Loading
#
#   3) CONVERSION LATTICE --> UNITCELL
#
#   4) OBTAINING GENERAL INFORMATION
#       - Testing a Unitcell object
#       - Print Information (printInfo(..))
#       - obtain the connection strengths of an object
#
#   5) OBTAINING CONNECTIVITY INFORMATION
#       - getConnectivityList(..)
#       - getConnectionList(..)
#
################################################################################

# included in subfile (1153 lines)
include("LatticePhysics_type_definitions.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT UNITCELLS in 2D and 3D
#
#   STRUCTURE OF THE FILE
#
#   2D UNITCELLS
#   - SQUARE / RECTANGLE
#       - getUnitcellSquare
#       - getUnitcellExtendedSquare
#       - getUnitcellCheckerboard
#       - getUnitcellShastrySutherland
#       - getUnitcellAdvancedSquare
#       - getUnitcellSquareOctagon
#       - getUnitcellBCC2D
#       - getUnitcellFullyConnectedSquare
#   - TRIANGULAR
#       - getUnitcellTriangular
#       - getUnitcellHoneycomb
#       - getUnitcellKagome
#       - getUnitcellKagomeMinus
#       - getUnitcellHoneycombXXX
#
#   3D UNITCELLS
#   - CUBIC / FCC
#       - getUnitcellDiamond
#       - getUnitcellCubic
#       - getUnitcellSC
#       - getUnitcellFCC
#       - getUnitcellBCC
#       - getUnitcellPyrochlore
#   - (X,3)y FAMILY
#       - getUnitcell_8_3_a
#       - getUnitcell_8_3_b
#       - getUnitcell_8_3_c
#       - getUnitcell_8_3_n
#       - getUnitcell_9_3_a
#       - getUnitcell_10_3_a / getUnitcellHyperoctagon
#       - getUnitcell_10_3_b / getUnitcellHyperhoneycomb
#       - getUnitcell_10_3_c
#       - getUnitcell_10_3_d
#
################################################################################

# included in subfile (4315 lines)
include("LatticePhysics_unitcell_implementations.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT UNITCELL CONSTRUCTION RELATED STUFF
#
#   STRUCTURE OF THE FILE
#
#   1) UNITCELL FROM SITE COLLECTION (automatically determine the connections)
#       - 2D (not exported)
#       - 3D (not exported)
#       - independent of the dimension
#
#   2) UNITCELL GENERATING CODE
#
################################################################################

# included in subfile (554 lines)
include("LatticePhysics_unitcell_construction.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE CONSTRUCTION RELATED STUFF
#
#   STRUCTURE OF THE FILE
#
#   1) CONSTRUCTION BASED ON UNITCELLS (2D & 3D)
#   - Periodic Boundaries
#   - Open Boundaries
#   - Semi-Periodic Boundaries
#   - General Boundaries
#
#   2) CONSTRUCTION BASED ON BOND DISTANCE (2D & 3D)
#
#   3) CONSTRUCTION BASED ON SHAPE (2D & 3D)
#   - General Shape
#   - Sphere
#   - Box
#
################################################################################

# included in subfile (2079 lines)
include("LatticePhysics_lattice_construction.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE MODIFICATION FUNCTIONS
#   (MOST OF THEM ALSO WORK FOR UNITCELLS)
#
#   STRUCTURE OF THE FILE
#
#   1) CONNECTIONS ADDING REMOVING
#       - Add a new connection
#       - remove connections (based on indices)
#       - remove connections (based on strength)
#
#   2) SITES ADDING REMOVING
#       - add a new site
#       - remove a site (based on index)
#
#   3) CONNECTION STRENGTH MODIFICATION
#       - all connections
#       - mapping of strengths
#       - evaluate strengths
#       - optimize connections
#
#   4) REAL SPACE ROTATION / SCALING / SHIFTING
#       - Rotation around some axis
#       - Scaling along some axis
#       - Global (isotropic) scaling
#       - Global (isotropic) scaling to mean / minimum bond length
#       - Shifting along some axis
#       - Shifting along some lattice vector
#       - Shifting center of lattice to origin
#
#   5) LATTICE VECTOR OPTIMIZATION
#
################################################################################

# included in subfile (1451 lines)
include("LatticePhysics_lattice_modification.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT PLAQUETTE RELATED FUNCTIONS
#
#   STRUCTURE OF THE FILE
#
#   1) FINDING PLAQUETTES IN LATTICE OBJECTS
#   - of a site
#   - of complete lattice
#
#   2) PRINTING PLAQUETTE STATISTICS
#   - of a site
#   - of complete lattice
#
#
#   BUG: Plaquettes dont implement periodic boundaries so far
#
################################################################################

# included in subfile (333 lines)
include("LatticePhysics_plaquette_operations.jl")

################################################################################
#
#   METHODS FOR CONSTRUCTION INTERACTION MATRICES FOR LATTICES
#
#   STRUCTURE OF THE FILE
#
#   1) INTERACTION MATRICES IN REAL SPACE
#
#   2) INTERACTION MATRICES IN MOMENTUM SPACE
#
#   NOTE: No Majorana fermions so far on the level of matrices, because of some
#   problems with the various gauges that can be tuned. Want: sublattice
#   idenfitication in bipartite systems.
#
################################################################################

# included in subfile (246 lines)
include("LatticePhysics_interaction_matrices.jl")

################################################################################
#
#   METHODS FOR CONSTRUCTION OF PATHS INSIDE THE BZ OF A UNITCELL
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE PATH
#       - type definition
#       - printInfo function
#       - path String function (e.g. X-G-M-X)
#
#   2) CONSTRUCTION FUNCTIONS PATH
#       - function to add points
#       - TODO function to remove points
#       - function to scale the total resolution
#       - function to set the total resolution
#
#   3) DEFAULT PATHS
#       - getDefaultPathTriangular
#       - getDefaultPathSquare
#       - getDefaultPathSquareOctagon
#       - getDefaultPathFCC
#
#   4) TODO FUNCTION TO CREATE A DEFAULT PATH BASED ON A UNITCELL OBJECT
#
#   5) FUNCTION TO PLOT A PATH (in both 2D and 3D)
#
################################################################################

# included in subfile (799 lines)
include("LatticePhysics_paths.jl")

################################################################################
#
#   METHODS FOR CONSTRUCTION AND PLOTTING OF BRILLOUIN ZONES (BZ) OF A UNITCELL
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE BRILLOUINZONE
#       - type definition
#       - TODO printInfo function
#
#   2) DEFAULT BZ
#       - TODO getDefaultBZSquare
#       - getDefaultBZFCC
#
#   3) FUNCTION TO CREATE A DEFAULT BZ BASED ON A UNITCELL OBJECT
#       - for 2D unitcells
#       - for 3D unitcells
#
#   4) PLOTTING OF A BZ
#
################################################################################

# included in subfile (1265 lines)
include("LatticePhysics_brillouin_zones.jl")

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE PLOTTING FUNCTIONS FOR
#   PLOTTING TO SVG IMAGES
#
#   STRUCTURE OF THE FILE
#
#   1) HELPER SVG FUNCTIONS (NOTE: not exported)
#       - Header / Footer String
#       - Ellipse (stroked) String
#       - Line String
#       - Plaquette String
#       - Color conversion
#       - Color collections / sequences
#
#   2) PLOTTING LATTICES
#       - plot in 2D
#       - TODO plot in 3D
#       - TODO plot independent of dimension
#
#   3) PLOTTING PLAQUETTES (only 2D)
#
################################################################################

# included in subfile (2094 lines)
include("LatticePhysics_SVG_plotting.jl")

################################################################################
#
#   METHODS FOR DUMPING A LATTICE TO BLENDER INPUT FILE
#   (input file can be used in Blender with the provided AddOn)
#
#   STRUCTURE OF THE FILE
#
#   1) SAVING LATTICE DATA TO TEXT FILE FOR BLENDER USE
#
################################################################################

# included in subfile (450 lines)
include("LatticePhysics_blender.jl")

################################################################################
#
#   METHODS FOR CONSTRUCTION OF BANDSTRUCTURES (ALONG PATHS)
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE BANDSTRUCTURE
#       - type definition
#       - printInfo function
#
#   2) CALCULATION OF BAND STRUCTURES OF UNTICELL OBJECTS
#
#   3) PLOTTING OF BAND STRUCTURES
#       - plotting of Bandstructure objects
#       - plotting of bandstructures of unitcells along paths
#
################################################################################

# included in subfile (529 lines)
include("LatticePhysics_bandstructures.jl")

################################################################################
#
#   METHODS FOR CONSTRUCTION OF FERMI SURFACES (2D & 3D)
#
#   STRUCTURE OF THE FILE
#
#   1) CALCULATION OF FERMI SURFACE
#
#   2) PLOTTING OF FERMI SURFACE
#       - plotting from points
#       - plotting from unitcell
#
################################################################################

# included in subfile (788 lines)
include("LatticePhysics_fermi_surfaces.jl")

################################################################################
#
#   METHODS FOR LUTTINGER TISZA CALCULATION
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE LTBANDSTRUCTURE
#       - type definition
#       - printInfo function
#
#   2) CALCULATION OF LT BAND STRUCTURES OF UNTICELL OBJECTS
#      (- LT constraint and deviation functions (NOT EXPORTED) )
#       - spin interaction matrices
#       - calculation of band structures
#
#   3) PLOTTING OF LT BAND STRUCTURES
#       - plotting of Bandstructure objects
#       - plotting of bandstructures of unitcells along paths
#
#   4) CALCULATION OF LT GROUND STATES (k space manifold)
#
#   5) PLOTTING OF LT GROUND STATES (k space manifold)
#       - plotting from points
#       - plotting from unitcell
#
#   BUG: So far, the constraint in all calculation ONLY works for k=0
#        For all other k values, the constraint is checked also for k=0
#
#   NOTE: Constraint has to be check with exp(ikr) factors on other sites as well
#
#   NOTE: Factor of 1/2 is missing for correct minimal energy
#
################################################################################

# included in subfile (1871 lines)
include("LatticePhysics_luttinger_tisza.jl")

# MODULE END
# total lines: 608 + 1153 + 4315 + 554 + 2079 + 1451 + 333 + 246 + 799 + 1265 + 2094 + 450 + 529 + 788 + 1871
# = 18535 lines
end


text:
@testset "NumField/NfAbs/NfAbs" begin
  cyclo_expl = function(n, m)
    Fn, zn = CyclotomicField(n)
    Fnm, znm = CyclotomicField(n*m)
    x = zn
    x_up = Hecke.force_coerce_cyclo(Fnm, x)
    x_down = Hecke.force_coerce_cyclo(Fn, x_up)
    return (x, x_up, x_down)
  end

  res = cyclo_expl(3, 4)
  @test (res[1]^3, res[2]^3) == (1, 1)

  res = cyclo_expl(3, 2)
  z6 = gen(parent(res[2]))
  # Test that z3 is mapped to z6^2
  @test z6^2 == res[2]
end

@testset "Splitting Field" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^3-2
  K = splitting_field(f)
  @test typeof(K) == AnticNumberField
  K1 = number_field([x^3-2, x^2+x+1])[1]
  K1abs = simple_extension(K1)[1]
  @test isisomorphic(K, K1abs)[1]
  K, R = splitting_field(f, do_roots = true)
  for r in R
    @test iszero(f(r))
  end

  K, a = number_field(x^2+1)
  Kt, t = PolynomialRing(K, "t")
  g = t^4-2
  L = splitting_field(g)
  @test typeof(L) == Hecke.NfRel{nf_elem}
  @test absolute_degree(L) == 8
end


text:
# arnoldi.jl

mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S}
    k::Int # current Krylov dimension
    V::OrthonormalBasis{T} # basis of length k
    H::Vector{S} # stores the Hessenberg matrix in packed form
    r::T # residual
end

Base.length(F::ArnoldiFactorization) = F.k
Base.sizehint!(F::ArnoldiFactorization, n) = begin
    sizehint!(F.V, n)
    sizehint!(F.H, (n*n + 3*n) >> 1)
    return F
end
Base.eltype(F::ArnoldiFactorization) = eltype(typeof(F))
Base.eltype(::Type{<:ArnoldiFactorization{<:Any,S}}) where {S} = S

basis(F::ArnoldiFactorization) = F.V
rayleighquotient(F::ArnoldiFactorization) = PackedHessenberg(F.H, F.k)
residual(F::ArnoldiFactorization) = F.r
@inbounds normres(F::ArnoldiFactorization) = abs(F.H[end])
rayleighextension(F::ArnoldiFactorization) = SimpleBasisVector(F.k, F.k)

# Arnoldi iteration for constructing the orthonormal basis of a Krylov subspace.
struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
    operator::F
    x₀::T
    orth::O
end
ArnoldiIterator(A, x₀) = ArnoldiIterator(A, x₀, KrylovDefaults.orth)
ArnoldiIterator(A::AbstractMatrix, x₀::AbstractVector, orth::Orthogonalizer =
                KrylovDefaults.orth) = ArnoldiIterator(x->A*x, x₀, orth)

Base.IteratorSize(::Type{<:ArnoldiIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:ArnoldiIterator}) = Base.EltypeUnknown()

function Base.iterate(iter::ArnoldiIterator)
    state = initialize(iter)
    return state, state
end
function Base.iterate(iter::ArnoldiIterator, state)
    nr = normres(state)
    if nr < eps(typeof(nr))
        return nothing
    else
        state = expand!(iter, deepcopy(state))
        return state, state
    end
end

function initialize(iter::ArnoldiIterator; verbosity::Int = 0)
    β₀ = norm(iter.x₀)
    β₀ == 0 && throw(ArgumentError("initial vector should not have norm zero"))
    invβ₀ = one(eltype(iter.x₀))/β₀
    T = typeof(invβ₀) # division might change eltype
    x₀ = mul!(similar(iter.x₀, T), iter.x₀, invβ₀)
    w = iter.operator(x₀) # applying the operator might change eltype
    v = eltype(x₀) == eltype(w) ? x₀ : copyto!(similar(w), x₀)
    r, α = orthogonalize!(w, v, iter.orth)
    β = norm(r)
    V = OrthonormalBasis([v])
    H = [α, β]
    if verbosity > 0
        @info "Arnoldi iteration step 1: normres = $β"
    end

    state = ArnoldiFactorization(1, V, H, r)
end
function initialize!(iter::ArnoldiIterator, state::ArnoldiFactorization; verbosity::Int = 0)
    x₀ = iter.x₀
    V = state.V
    while length(V) > 1
        pop!(V)
    end
    H = empty!(state.H)

    v = mul!(V[1], x₀, 1/norm(x₀))
    w = iter.operator(v)
    r, α = orthogonalize!(w, v, iter.orth)
    β = norm(r)
    state.k = 1
    push!(H, α, β)
    state.r = r
    if verbosity > 0
        @info "Arnoldi iteration step 1: normres = $β"
    end
    return state
end
function expand!(iter::ArnoldiIterator, state::ArnoldiFactorization; verbosity::Int = 0)
    state.k += 1
    k = state.k
    V = state.V
    H = state.H
    r = state.r
    β = normres(state)
    push!(V, rmul!(r, 1/β))
    m = length(H)
    resize!(H, m+k+1)
    r, β = arnoldirecurrence!(iter.operator, V, view(H, (m+1):(m+k)), iter.orth)
    H[m+k+1] = β
    state.r = r
    if verbosity > 0
        @info "Arnoldi iteration step $k: normres = $β"
    end
    return state
end
function shrink!(state::ArnoldiFactorization, k)
    length(state) <= k && return state
    V = state.V
    H = state.H
    while length(V) > k+1
        pop!(V)
    end
    r = pop!(V)
    resize!(H, (k*k + 3*k) >> 1)
    state.k = k
    state.r = rmul!(r, normres(state))
    return state
end

# Arnoldi recurrence: simply use provided orthonormalization routines
function arnoldirecurrence!(operator, V::OrthonormalBasis, h::AbstractVector,
                            orth::Orthogonalizer)
    w = operator(last(V))
    r, h = orthogonalize!(w, V, h, orth)
    return r, norm(r)
end


text:
function solve()
    n = parse(Int, readline())
    zoro = 0
    for i=1:n
        d1, d2 = split(readline())
        if d1 == d2
            zoro += 1
        else
            zoro = 0
        end
        if zoro ≥ 3
            return "Yes"
        end
    end
    "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia