text:
#
# This difference is important because for Orthorhombic cells it is
# possible to run over only half of the cells, and wrapping coordinates
# in Orthorhombic cells is slightly cheaper.
#
struct TriclinicCell end
struct OrthorhombicCell end

struct UnitCell{UnitCellType,N,T,M}
    matrix::SMatrix{N,N,T,M}
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Structure that contains the maximum lengths on each direction,
to dispatch on the construction of boxes without periodic boundary
conditions.

"""
struct Limits{T<:AbstractVector}
    limits::T
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Structure that contains some data required to compute the linked cells. To
be initialized with the box size and cutoff.

## Examples

```julia-repl
julia> using CellListMap

julia> sides = [250,250,250];

julia> cutoff = 10;

julia> box = Box(sides,cutoff)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [250.0 0.0 0.0; 0.0 250.0 0.0; 0.0 0.0 250.0]
  cutoff: 10.0
  number of computing cells on each dimension: [27, 27, 27]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 19683

```

```julia-repl
julia> box = Box([ 10  0  0
                    0 10  5
                    0  0 10 ], 1)
Box{TriclinicCell, 3, Float64, 9}
unit cell matrix: [10.0 0.0 0.0; 0.0 10.0 5.0; 0.0 0.0 10.0]
cutoff: 1.0
number of computing cells on each dimension: [12, 17, 12]
computing cell sizes: [1.0, 1.0, 1.0] (lcell: 1)
Total number of cells: 2448

```

"""
Base.@kwdef struct Box{UnitCellType,N,T,TSQ,M}
    unit_cell::UnitCell{UnitCellType,N,T,M}
    lcell::Int
    nc::SVector{N,Int}
    cutoff::T
    cutoff_sq::TSQ
    ranges::SVector{N,UnitRange{Int}}
    cell_size::SVector{N,T}
    unit_cell_max::SVector{N,T}
end

"""

```
_promote_types(cell,cutoff)
```

Internal function or structure - interface may change.

# Extended help

Promotes the types of the unit cell matrix (or sides) and cutoff to floats if one or both were input as integers.

"""
function _promote_types(cell,cutoff)
    if eltype(cell) <: Integer || typeof(cutoff) <: Integer
        if cutoff isa Integer && eltype(cell) <: Integer
            cell = convert.(Float64,cell)
            cutoff = convert(eltype(cell),cutoff)
        elseif eltype(cell) <: Integer
            cell = convert.(typeof(cutoff),cell)
        else
            cutoff = convert(eltype(cell),cutoff)
        end
    end
    return cell, cutoff
end

"""

```
Box(
  unit_cell_matrix::AbstractMatrix,
  cutoff,
  lcell::Int=1,
  UnitCellType=TriclinicCell
)
```

Construct box structure given the cell matrix of lattice vectors. This
constructor will always return a `TriclinicCell` box type, unless the
`UnitCellType` parameter is set manually to `OrthorhombicCell`

## Example
```julia-repl
julia> unit_cell = [ 100   50    0
                       0  120    0
                       0    0  130 ];

julia> box = Box(unit_cell,10)
Box{TriclinicCell, 3, Float64, 9}
  unit cell matrix: [100.0 50.0 0.0; 0.0 120.0 0.0; 0.0 0.0 130.0]
  cutoff: 10.0
  number of computing cells on each dimension: [17, 14, 15]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 3570

```

"""
function Box(
    unit_cell_matrix::AbstractMatrix,
    cutoff,
    lcell::Int,
    ::Type{UnitCellType}
) where {UnitCellType}
    unit_cell_matrix, cutoff = _promote_types(unit_cell_matrix, cutoff)
    T = eltype(unit_cell_matrix)

    s = size(unit_cell_matrix)
    unit_cell_matrix = SMatrix{s[1],s[2],T,s[1]*s[2]}(unit_cell_matrix)

    @assert lcell >= 1 "lcell must be greater or equal to 1"

    N = size(unit_cell_matrix)[1]
    @assert N == size(unit_cell_matrix)[2] "Unit cell matrix must be square."
    @assert check_unit_cell(unit_cell_matrix,cutoff) " Unit cell matrix does not satisfy required conditions."

    unit_cell = UnitCell{UnitCellType,N,T,N*N}(unit_cell_matrix)
    unit_cell_max = sum(unit_cell_matrix[:,i] for i in 1:N)

    # To use the advantages of Orthorhombic cells, box size must
    # be  multiple of the cell size. In some pathological cases this
    # may be bad for performance, because we are increasing the
    # effective cell cutoff
    if UnitCellType == OrthorhombicCell
        nc = floor.(Int,lcell*unit_cell_max/cutoff)
        cell_size = unit_cell_max ./ nc
        nc = nc .+ 2*lcell
    # For triclinic cells we have to use the ceil here, because we need
    # to major the number of cells, when the box size is not a multiple of
    # the cell size
    else
        cell_size = SVector{N,T}(ntuple(i->cutoff/lcell,N)...)
        nc = ceil.(Int,(unit_cell_max .+ 2*cutoff) ./ cell_size)
    end

    #ranges = ranges_of_replicas(cell_size, lcell, nc, unit_cell_matrix)
    ranges = SVector{N,UnitRange{Int}}(ntuple(i->-1:1,N))
    return Box{UnitCellType,N,T,typeof(cutoff^2),N*N}(
        unit_cell,
        lcell,
        nc,
        cutoff,
        cutoff^2,
        ranges,
        cell_size,
        unit_cell_max
    )
end
Box(
    unit_cell_matrix::AbstractMatrix,
    cutoff;
    lcell::Int=1,
    UnitCellType=TriclinicCell
) = Box(unit_cell_matrix,cutoff,lcell,UnitCellType)

function Base.show(io::IO,::MIME"text/plain",box::Box{UnitCellType,N}) where {UnitCellType,N}
    println(io,typeof(box))
    print(io,"  unit cell matrix = [ ")
    print(io,join(box.unit_cell.matrix[1:N,1],", "))
    for i in 2:N
        print(io,"; ",join(box.unit_cell.matrix[1:N,i],", "))
    end
    println(io," ]")
    println(io,"  cutoff = ", box.cutoff)
    println(io,"  number of computing cells on each dimension = ",box.nc)
    println(io,"  computing cell sizes = [", join(box.cell_size,", "), "] (lcell: ",box.lcell,")")
    print(io,"  Total number of cells = ", prod(box.nc))
end

"""

```
cell_matrix_from_sides(sides::AbstractVector)
```

Internal function or structure - interface may change.

# Extended help

Function that returns the Orthorhombic unit cell matrix given a sides vector.

## Example

```
julia> CellListMap.cell_matrix_from_sides([1,1,1])
3×3 SMatrix{3, 3, Int64, 9} with indices SOneTo(3)×SOneTo(3):
 1  0  0
 0  1  0
 0  0  1
```

"""
function cell_matrix_from_sides(sides::AbstractVector)
    T = eltype(sides)
    N = length(sides)
    cart_idxs = CartesianIndices((1:N,1:N))
    unit_cell_matrix = SMatrix{N,N,T,N*N}(
        ntuple(N*N) do i
            c = cart_idxs[i]
            if c[1] == c[2]
                return sides[c[1]]
            else
                return zero(T)
            end
        end
    )
    return unit_cell_matrix
end

"""

```
Box(
  sides::AbstractVector,
  cutoff,
  lcell::Int=1,
  UnitCellType=OrthorhombicCell
)
```

For orthorhombic unit cells, `Box` can be initialized with a vector of the
length of each side.

## Example
```julia-repl
julia> box = Box([120,150,100],10)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [120.0 0.0 0.0; 0.0 150.0 0.0; 0.0 0.0 100.0]
  cutoff: 10.0
  number of computing cells on each dimension: [14, 17, 12]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 2856

```

"""
function Box(
    sides::AbstractVector,
    cutoff,
    lcell::Int,
    ::Type{UnitCellType}
) where {UnitCellType}
    sides, cutoff = _promote_types(sides, cutoff)
    unit_cell_matrix = cell_matrix_from_sides(sides)
    return Box(unit_cell_matrix,cutoff,lcell,UnitCellType)
end
Box(
    sides::AbstractVector,
    cutoff;
    lcell::Int=1,
    UnitCellType=OrthorhombicCell
) = Box(sides,cutoff,lcell,UnitCellType)

"""

```
Box(
    limits::Limits,
    cutoff;
    lcell::Int=1
)
```

This constructor receives the output of `limits(x)` or `limits(x,y)` where `x` and `y` are
the coordinates of the particles involved, and constructs a `Box` with size larger than
the maximum coordinates ranges of all particles plus the cutoff. This is used to
emulate pairwise interactions in non-periodic boxes. The output box is always an `Orthorhombic`
cell.

## Examples

```julia-repl
julia> x = [ [100,100,100] .* rand(3) for i in 1:100_000 ];

julia> box = Box(limits(x),10)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [109.99633932875878 0.0 0.0; 0.0 109.99780283179763 0.0; 0.0 0.0 109.99587254766517]
  cutoff: 10.0
  number of computing cells on each dimension: [12, 12, 12]
  computing cell sizes: [10.999633932875877, 10.999780283179764, 10.999587254766517] (lcell: 1)
  Total number of cells: 1728

julia> y = [ [150,150,50] .* rand(3) for i in 1:100_000 ];

julia> box = Box(limits(x,y),10)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [159.99787690924168 0.0 0.0; 0.0 159.98878289444897 0.0; 0.0 0.0 109.99587254766517]
  cutoff: 10.0
  number of computing cells on each dimension: [18, 17, 12]
  computing cell sizes: [10.666525127282778, 10.665918859629931, 10.999587254766517] (lcell: 1)
  Total number of cells: 3672

```

"""
function Box(limits::Limits, cutoff::T; lcell::Int=1) where T
    sides = limits.limits .+ cutoff
    return Box(sides, cutoff, lcell, OrthorhombicCell)
end


text:
# CUDA module management
# This uses the driver API

function checkdrv(code::Integer)
    if code != 0
        error(driver_error_descriptions[int(code)])
    end
    nothing
end

function checkdrv(code::Integer, msg::String)
    if code != 0
        warn(msg)
        error(driver_error_descriptions[int(code)])
    end
    nothing
end

immutable CuModule
    handle::Ptr{Void}

    function CuModule(filename::String, finalize::Bool = true)
        a = Array(Ptr{Void}, 1)
        checkdrv(ccall((:cuModuleLoad, libcuda), Cint, (Ptr{Ptr{Void}}, Ptr{Cchar}), a, filename), filename)
        md = new(a[1])
        if finalize
            finalizer(md, unload)
        end
        md
    end
    CuModule() = new(C_NULL)
end

# do syntax, f(module)
function CuModule(f::Function, filename::String)
    md = CuModule(filename, false)
    local ret
    try
        ret = f(md)
    finally
        unload(md)
    end
    ret
end

function unload(md::CuModule)
    checkdrv(ccall((:cuModuleUnload, libcuda), Cint, (Ptr{Void},), md.handle))
end

immutable CuFunction
    handle::Ptr{Void}

    function CuFunction(md::CuModule, name::String)
        a = Array(Ptr{Void}, 1)
        checkdrv(ccall((:cuModuleGetFunction, libcuda), Cint,
                       (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Cchar}),
                       a, md.handle, name))
        new(a[1])
    end
end

const driver_error_descriptions = @compat Dict{Int,String}(
        0 => "Success",
        1 => "Invalid value",
        2 => "Out of memory",
        3 => "Driver not initialized",
        4 => "Driver being shutdown",
        5 => "Profiler disabled",
        6 => "Profiler not initialized",
        7 => "Profiler already started",
        8 => "Profiler already stopped",
        100 => "No CUDA-capable device",
        101 => "Invalid device ordinal",
        200 => "Invalid kernel image",
        201 => "Invalid context",
        202 => "Context already current",
        205 => "Map operation failed",
        206 => "Unmap operation failed",
        207 => "Array mapped",
        208 => "Resource already mapped",
        209 => "No kernel image available/suitable for GPU",
        210 => "Resource already acquired",
        211 => "Resource not mapped",
        212 => "Resource not mapped as array",
        213 => "Resource not mapped as pointer",
        214 => "Uncorrectable ECC error detected",
        215 => "Unsupported limit",
        216 => "Context already in use",
        217 => "Peer access not supported",
        300 => "Invalid kernel source",
        301 => "File not found",
        302 => "Shared object symbol not found",
        303 => "Shared object initialization failed",
        304 => "OS call failed",
        400 => "Invalid handle",
        500 => "Named symbol not found",
        600 => "Not ready",
        700 => "Kernel launch failed",
        701 => "Launch out of resources",
        702 => "Launch timeout",
        703 => "Incompatible texturing mode",
        704 => "Peer access already enabled",
        705 => "Peer access not enabled",
        708 => "Primary context already active",
        709 => "Context destroyed",
        710 => "Assertion triggered failure",
        711 => "Too many peers",
        712 => "Host memory already registered",
        713 => "Host memory not registered",
        719 => "Launch failed, perhaps due to an invalid pointer",
        800 => "Operation not permitted",
        801 => "Operation not supported",
        999 => "Unknown error"
)

ext: jl


text:
export SpelledInterval, spelled, spelledp
export SpelledIC, sic, spc
export octaves, internalocts, fifths
export degree, generic, diasteps, alteration, letter
export parsespelled, parsespelledpitch, @i_str, @p_str

import Base: +, -, *, ==

# helpers

const perfectints = Set{Int}([0,3,4])
accstr(n, u, d) =
    if n > 0; repeat(u,n) elseif n < 0; repeat(d,-n) else "" end
qualpf(n, a, p, d) =
    if n > 0; repeat(a,n) elseif n < 0; repeat(d,-n) else p end
qualimpf(n, a, mj, mn, d) =
    if n > 0; repeat(a,n) elseif n < -1; repeat(d,-n-1) elseif n == -1; mn else mj end

# common functions (special interface)

"""
    degree(i)
    degree(p)

Return the "relative scale degree" (0-6) to which the interval points
(unison=`0`, 2nd=`1`, octave=`0`, 2nd down=`6`, etc.).
For pitches, return the integer that corresponds to the letter (C=`0`, D=`1`, ...).

See also: [`generic`](@ref), [`diasteps`](@ref), [`letter`](@ref)
"""
function degree end
degree(fifths::Int) = mod(fifths*4,7)
degree(p::Pitch) = degree(p.pitch)

"""
    generic(i)

Return the generic interval, i.e. the number of diatonic steps modulo octave.
Unlike [`degree`](@ref), the result respects the sign of the interval
(unison=`0`, 2nd up=`1`, 2nd down=`-1`).
For pitches, use [`degree`](@ref).

See also: [`degree`](@ref), [`diasteps`](@ref)
"""
function generic end

"""
    letter(p)

Returns the letter of a pitch as a character.

See also: [`degree`](@ref).
"""
letter(p::Pitch) = 'A' + mod(degree(p)+2, 7)

"""
    diasteps(i)

Return the diatonic steps of the interval (unison=`0`, 2nd=`1`, ..., octave=`7`).
Respects both direction and octaves.

See also [`degree`](@ref) and [`generic`](@ref).
"""
function diasteps end

"""
    alteration(i)
    alteration(p)

Return the number of semitones by which the interval is altered from its the perfect or major variant.
Positive alteration always indicates augmentation,
negative alteration indicates diminution (minor or smaller) of the interval.
For pitches, return the accidentals (positive=sharps, negative=flats, `0`=natural).
"""
function alteration end
alteration(p::Pitch) = alteration(ic(p.pitch))

"""
    octaves(i)
    octaves(p)

Return the number of octaves the interval spans.
Positive intervals start at 0 octaves, increasing.
Negative intervals start at -1 octaves, decreasing.
(You might want to use `octaves(abs(i))` instead).

For a pitch, return its octave.
"""
function octaves end
octaves(p::Pitch) = octaves(p.pitch)

"""
    internalocts(i)

Return the internal octaves (i.e. dependent on the 5ths dimension) of an interval.
"""
function internalocts end

"""
    fifths(i)
    fifths(p)

Return the octave-invariant part of the interval in fifths
(unison=`0`, 5th up=`1`, 4th up/5th down=`-1`).
For a pitch, return the pitch class on the line of fifths.
"""
function fifths end
fifths(p::Pitch) = fifths(p.pitch)

# spelled interval
# ----------------

"""
    SpelledInterval <: Interval

Spelled intervals represented as pairs of fifths and octaves.
E.g., `SpelledInterval(-3, 2)` represents a minor 3rd upwards
(3 fifths down, 2 octaves up).
"""
struct SpelledInterval <: Interval
    fifths :: Int
    octaves :: Int
end

"""
    spelled(fifths, octaves)

Creates a spelled interval from `fifths` and `octaves`.
"""
spelled(f, o) = SpelledInterval(f, o)

"""
    spelledp(fifths, octaves)

Creates a spelled pitch from `fifths` and `octaves`.
"""
spelledp(f, o) = Pitch(spelled(f, o))

# accessors

degree(i::SpelledInterval) = degree(i.fifths)
generic(i::SpelledInterval) = if sign(i) < 0; -degree(-i.fifths) else degree(i.fifths) end
diasteps(i::SpelledInterval) = i.fifths*4 + i.octaves*7
alteration(i::SpelledInterval) = let absi = abs(i); fld(absi.fifths+1, 7) end
octaves(i::SpelledInterval) = i.octaves + fld(i.fifths*4, 7)
internalocts(i::SpelledInterval) = i.octaves
fifths(i::SpelledInterval) = i.fifths

# general interface

function Base.show(io::IO, i::SpelledInterval)
   # negative? print as -abs(i)
    if sign(i) == -1
        print(io, "-")
        print(io, abs(i))
        return
    end

    dia = degree(i)
    diff = alteration(i) # interval is always positive, so direction is correct
    qual = if dia ∈ perfectints
        qualpf(diff, 'a', 'P', 'd')
    else
        qualimpf(diff, 'a', 'M', 'm', 'd')
    end

    octs = octaves(i)

    print(io, qual * string(dia+1) * ':' * string(octs))
end

function Base.show(io::IO, p::Pitch{SpelledInterval})
    i = p.pitch
    dia = degree(i)
    alter = alteration(i)
    print(io, string('A' + mod(dia + 2, 7), accstr(alter, '♯', '♭'), string(octaves(i))))
end

Base.isless(i1::SpelledInterval, i2::SpelledInterval) =
    (diasteps(i1), alteration(i1)) <
    (diasteps(i2), alteration(i2))

Base.isequal(i1::SpelledInterval, i2::SpelledInterval) =
    i1.octaves == i2.octaves && i1.fifths == i2.fifths

Base.hash(i::SpelledInterval, x::UInt) = hash(i.octaves, hash(i.fifths, x))

+(i1::SpelledInterval, i2::SpelledInterval) =
    spelled(i1.fifths + i2.fifths, i1.octaves + i2.octaves)
-(i1::SpelledInterval, i2::SpelledInterval) =
    spelled(i1.fifths - i2.fifths, i1.octaves - i2.octaves)
-(i::SpelledInterval) = spelled(-i.fifths, -i.octaves)
Base.zero(::Type{SpelledInterval}) = spelled(0,0)
Base.zero(::SpelledInterval) = spelled(0,0)
*(i::SpelledInterval, n::Integer) = spelled(i.fifths*n, i.octaves*n)
*(n::Integer,i::SpelledInterval) = spelled(i.fifths*n, i.octaves*n)

octave(::Type{SpelledInterval}) = spelled(0,1)

Base.sign(i::SpelledInterval) = sign(diasteps(i))
Base.abs(i::SpelledInterval) = if sign(i) < 0; -i else i end

ic(i::SpelledInterval) = sic(i.fifths)
embed(i::SpelledInterval) = i
intervaltype(::Type{SpelledInterval}) = SpelledInterval
intervalclasstype(::Type{SpelledInterval}) = SpelledIC

isstep(i::SpelledInterval) = abs(diasteps(i)) <= 1
chromsemi(::Type{SpelledInterval}) = spelled(7,-4)

# conversion

tomidi(i::SpelledInterval) = midi(7*i.fifths + 12*i.octaves)
tomidi(p::Pitch{SpelledInterval}) = Pitch(tomidi(p.pitch) + midi(12))  # C4 = 48 semitones above C0 = midi(60)

# spelled interval class
# ----------------------

"""
    SpelledIC <: IntervalClass

Spelled interval class represented on the line of 5ths with `0 = C`.
E.g., `SpelledIC(3)` represents a major 6th upwards or minor 3rd downwards
(i.e., three 5ths up modulo octave).
"""
struct SpelledIC <: IntervalClass
    fifths :: Int
end

"""
    sic(fifths)

Creates a spelled interval class going `fifths` 5ths upwards.
"""
sic(fs) = SpelledIC(fs)

"""
    spc(fifths)

Creates a spelled pitch class.
In analogy to `sic`, this function takes a number of 5ths.
"""
spc(fs) = Pitch(sic(fs))

# accessors

degree(i::SpelledIC) = degree(i.fifths)
generic(i::SpelledIC) = degree(i.fifths)
diasteps(i::SpelledIC) = degree(i.fifths)
alteration(i::SpelledIC) = fld(i.fifths + 1, 7)
octaves(i::SpelledIC) = 0
internalocts(i::SpelledIC) = 0
fifths(i::SpelledIC) = i.fifths

# interface functions

function Base.show(io::IO, ic::SpelledIC)
    dia = degree(ic)
    diff = alteration(ic)
    qual = if dia ∈ perfectints
        qualpf(diff, 'a', 'P', 'd')
    else
        qualimpf(diff, 'a', 'M', 'm', 'd')
    end

    print(io, qual * string(dia+1))
end

function Base.show(io::IO, p::Pitch{SpelledIC})
    i = p.pitch
    dia = degree(i)
    alter = alteration(i)
    print(io, ('A' + mod(dia+2, 7)) * accstr(alter, '♯', '♭'))
end

Base.isless(i1::SpelledIC, i2::SpelledIC) = isless(i1.fifths,i2.fifths)
Base.isequal(i1::SpelledIC, i2::SpelledIC) = isequal(i1.fifths,i2.fifths)
Base.hash(i::SpelledIC, x::UInt) = hash(i.fifths, x)

+(i1::SpelledIC, i2::SpelledIC) = sic(i1.fifths + i2.fifths)
-(i1::SpelledIC, i2::SpelledIC) = sic(i1.fifths - i2.fifths)
-(i::SpelledIC) = sic(-i.fifths)
Base.zero(::Type{SpelledIC}) = sic(0)
Base.zero(::SpelledIC) = sic(0)
*(i::SpelledIC, n::Integer) = sic(i.fifths * n)
*(n::Integer,i::SpelledIC) = sic(i.fifths * n)

octave(::Type{SpelledIC}) = sic(0)
function Base.sign(i::SpelledIC)
    dia = degree(i)
    if dia == 0; 0 elseif dia > 3; -1 else 1 end
end
Base.abs(i::SpelledIC) = if sign(i) < 0; -i else i end

ic(i::SpelledIC) = i
embed(i::SpelledIC) = spelled(i.fifths, -fld(i.fifths*4, 7))
intervaltype(::Type{SpelledIC}) = SpelledInterval
intervalclasstype(::Type{SpelledIC}) = SpelledIC

isstep(i::SpelledIC) = degree(i) ∈ [0,1,6]
chromsemi(::Type{SpelledIC}) = sic(7)

# conversion

tomidi(i::SpelledIC) = midic(i.fifths * 7)
tomidi(p::Pitch{SpelledIC}) = midipc(p.pitch.fifths * 7)

# parsing
# -------

const rgsic = r"^(-?)(a+|d+|[MPm])([1-7])$"
const rgspelled = r"^(-?)(a+|d+|[MPm])([1-7]):(-?)(\d+)$"

function matchinterval(modifier, num)
    dia = parse(Int, num) - 1
    perfect = dia ∈ perfectints

    alt = if modifier == "M" && !perfect
        0
    elseif modifier == "m" && !perfect
        -1
    elseif lowercase(modifier) == "p" && perfect
        0
    elseif occursin(r"^a+$", modifier)
        length(modifier)
    elseif occursin(r"^d+$", modifier)
        -length(modifier) - (perfect ? 0 : 1)
    else
        error("cannot parse interval \"$modifier$num\"")
    end

    mod(dia*2+1, 7) - 1 + 7*alt
end

"""
    parsespelled(str)

Parse a spelled interval or interval class string.
The type is determined from the string,
so `i"M3:0"` returns an interval while `i"M3"` returns an interval class.

See also: [`@i_str`](@ref), [`parsespelledpitch`](@ref).
"""
function parsespelled(str)
    m = match(rgsic, str)
    if m != nothing
        int = sic(matchinterval(m[2], m[3]))
    else
        m = match(rgspelled, str)
        if m != nothing
            fifths = matchinterval(m[2], m[3])
            octs = parse(Int, m[4]*m[5])
            int = spelled(fifths, octs - fld(fifths*4, 7))
        else
            error("cannot parse interval \"$str\"")
        end
    end

    # invert if necessary
    if m[1] == "-"
        -int
    else
        int
    end
end

"""
    i"str"

Parse a spelled interval or interval class string.
The type is determined from the string,
so `i"M3:0"` returns an interval while `i"M3"` returns an interval class.

See also: [`@p_str`](@ref), [`parsespelled`](@ref), [`parsespelledpitch`](@ref).
"""
macro i_str(str)
    parsespelled(str)
end

const rgspelledpc = r"^([a-g])(♭+|♯+|b+|#+)?$"i
const rgspelledp = r"^([a-g])(♭+|♯+|b+|#+)?(-?\d+)$"i

function matchpitch(letter, accs)
    letter = uppercase(letter)[1]
    if letter >= 'A' && letter <= 'G'
        dia = mod(letter - 'A' - 2, 7)
    else
        error("cannot parse pitch letter \"$letter\"")
    end

    alt = if accs == nothing || accs == ""
        0
    elseif occursin(r"^♭+|b+$"i, accs)
        -length(accs)
    elseif occursin(r"^♯+|#+$"i, accs)
        length(accs)
    else
        error("cannot parse accidentals \"$accs\"")
    end

    mod(dia*2 + 1, 7) - 1 + 7*alt
end

"""
    parsespelledpitch(str)

Parse a spelled pitch or pitch class string.
The type is determined from the string,
so `p"G4"` returns a pitch while `p"G"` returns a pitch class.

See also: [`@p_str`](@ref), [`parsespelled`](@ref).
"""
function parsespelledpitch(str)
    m = match(rgspelledpc, str)
    if m != nothing
        spc(matchpitch(m[1], m[2]))
    else
        m = match(rgspelledp, str)
        if m != nothing
            octs = parse(Int, m[3])
            fifths = matchpitch(m[1], m[2])
            spelledp(fifths, octs - fld(fifths*4, 7))
        else
            error("cannot parse pitch \"$str\"")
        end
    end
end

"""
    p"str"

Parse a spelled pitch or pitch class string.
The type is determined from the string,
so `p"G4"` returns a pitch while `p"G"` returns a pitch class.

See also: [`@i_str`](@ref), [`parsespelledpitch`](@ref), [`parsespelled`](@ref).
"""
macro p_str(str)
    parsespelledpitch(str)
end


text:
"""
    version_major()

Return the *major* version of the used SEAL library as an integer.

See also: [`version_minor`](@ref), [`version_patch`](@ref), [`version`](@ref)
"""
function version_major()
  val = Ref{UInt8}(0)
  retval = ccall((:Version_Major, libsealc), Clong, (Ref{UInt8},), val)
  @check_return_value retval
  return convert(Int, val[])
end

"""
    version_minor()

Return the *minor* version of the used SEAL library as an integer.

See also: [`version_major`](@ref), [`version_patch`](@ref), [`version`](@ref)
"""
function version_minor()
  val = Ref{UInt8}(0)
  retval = ccall((:Version_Minor, libsealc), Clong, (Ref{UInt8},), val)
  @check_return_value retval
  return convert(Int, val[])
end

"""
    version_patch()

Return the *patch* version of the used SEAL library as an integer.

See also: [`version_major`](@ref), [`version_minor`](@ref), [`version`](@ref)
"""
function version_patch()
  val = Ref{UInt8}(0)
  retval = ccall((:Version_Patch, libsealc), Clong, (Ref{UInt8},), val)
  @check_return_value retval
  return convert(Int, val[])
end

"""
    version()

Return the version of the used SEAL library as a `VersionNumber` in the format
`v"major.minor.patch"`..

See also: [`version_major`](@ref), [`version_minor`](@ref), [`version_patch`](@ref)
"""
function version()
  major = version_major()
  minor = version_minor()
  patch = version_patch()
  return VersionNumber("$major.$minor.$patch")
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

n = 10
plot(x=rand(n), y=rand(n),
     Scale.x_continuous(minvalue=-10, maxvalue=10),
     Scale.y_continuous(minvalue=-10, maxvalue=10))


text:
module ProgenyTestingTools

using Random
using Printf
using Statistics
using SparseArrays
using LinearAlgebra

using CSV
using QMSimFiles
using DataFrames
using OffsetArrays
using SparseMatrixDicts
using IterativeSolvers
using StatsBase: sample, Weights

# genetic parameters
"""
    par = PTParameters(var_p, h2_poly, h2_qtl, rep)

It defines the genetic parameters used in data generation supporting a single trait.
The parameters include `var_p` for the phenotypic variance, `h2_poly` for
the polygenic heritability, `h2_qtl` for the qtl heritability, and `rep` for
the repeatability (for repeated records).
"""
struct PTParameters
   mu::Float64
   var_p::Float64
   h2_poly::Float64
   h2_qtl::Float64     # >0 if using genotypes
   rep::Float64
   var_g_sim::Float64  # true genetic variance in QMSim
end

# animal information
# multiple records but a single BV
mutable struct PTAnimal
   cg::Vector{Int}
   y::Vector{Float64}
   pe::Vector{Float64}
   e::Vector{Float64}
end

# central database for all animals
# will be updated when a new animal is born.
#
# maxCG: group x generation: increments when phenotype is assigned
"""
It defines a set of animals in a population.
The individual information will be stored in a dataframe `df` with the following columns.

- `id::Int` = sequential code (from 1)
- `male::Bool` = `true` for male, `false` for female
- `year::Int` = birth year
- `alive::Bool` = `true` for living individual
- `pregnant::Bool` = `true` for pregnant individual
- `genotyped::Bool` = `true` for genotyped individual
- `candidate::Bool` = `true` for candidate for selection
- `status::Int` = user-defined status code (default = 0)
- `sire::Int` = sire code
- `dam::Int` = dam code
- `siregroup::Int` = group ID of sire
- `damgroup::Int` = group ID of dam
- `nprog::Int` = number of progeny
- `nrecprog::Int` = number of progeny with record(s)
- `inb::Float64` = inbreeding coefficient (from 0.0 to 1.0)
- `pbv::Float64` = polygenic breeding value (PBV)
- `qbv::Float64` = QTL breeding value (QBV)
- `tbv::Float64` = true (total) breeding value: TBV = PBV + QBV
- `ebv::Float64` = estimated breeding value by pedigree BLUP
- `gebv::Float64` = genomic EBV by "pseudo" genomic prediction
- `rel::Float64` = user-supplied reliability of GEBV
- `nrec::Int` = number of record(s)
- `firsty::Float64` = the first observation
- `lasty::Float64` = the last observation
- `avgy::Float64` = average of observations
"""
mutable struct PTPopulation
   par::PTParameters
   hp::Bool
   maxAnimal::Int
   maxGroup::Int
   maxCG::Int
   df::DataFrame
   animal::Vector{PTAnimal}
   map::Union{Nothing,QMSimMap}
   gfile::String
end

# group: a subset of the population
# sires, dams, id: linked to IDs in a population
mutable struct PTGroup
   pop::PTPopulation
   groupid::Int
   n::Int
   maxSire::Int
   maxDam::Int
   sires::Vector{Int}
   dams::Vector{Int}
   id::Vector{Int}
   generation::Vector{Int}
end

import Base: show, write
export PTParameters, PTAnimal, PTPopulation, PTGroup
export check_parameters, get_var_poly, get_var_qtl, get_var_pe, get_var_error
export assign_phenotype!
export generate_population, migrate_from_hp!, generate_group, copy_group, add_new_animal!,
       selectid, random_sampling, add_sires!, add_dams!, vacancy_for_sires, vacancy_for_dams,
       mark_candidate!, change_status!, cull!
export update_inbreeding!
export mating!, calving!
export genetic_evaluation!, update_approximated_reliability!, write_files_for_blup

export test_mating!, regular_mating!, cull_old_bulls!, cull_old_cows!,
       mean_inbreeding, number_of_active_females, phenotyped_cows, expected_frequency_of_mating

include("parameter.jl")
include("population.jl")
include("phenotype.jl")
include("pedigree.jl")
include("mating.jl")
include("io.jl")
include("evaluation.jl")
include("batch.jl")

end


text:
@testset "WasmTable" begin
    wasm_code = wat"""
    (module
        (import "js" "tbl" (table 2 anyfunc))
            (func $f42 (result i32) i32.const 42)
            (func $f83 (result i32) i32.const 83)
            (elem (i32.const 0) $f42 $f83)
    )"""

    engine = WasmEngine()
    store = WasmStore(engine)

    wasm_module = WasmModule(store, wasm_code)

    table = Wasmtime.WasmTable(store, 2 => 5)

    @test Base.unsafe_convert(Ptr{Wasmtime.wasm_func_t}, table[1]) == C_NULL
    @test Base.unsafe_convert(Ptr{Wasmtime.wasm_func_t}, table[2]) == C_NULL

    instance = WasmInstance(store, wasm_module, [table])

    @test Base.unsafe_convert(Ptr{Wasmtime.wasm_func_t}, table[1]) != C_NULL
    @test Base.unsafe_convert(Ptr{Wasmtime.wasm_func_t}, table[2]) != C_NULL

    @test Base.convert(Int32, table[1]()[1]) == 42
    @test Base.convert(Int32, table[2]()[1]) == 83
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ProximityPlacementGroupListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{ProximityPlacementGroup} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function ProximityPlacementGroupListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(ProximityPlacementGroupListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(ProximityPlacementGroupListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type ProximityPlacementGroupListResult

const _property_map_ProximityPlacementGroupListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_ProximityPlacementGroupListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{ProximityPlacementGroup}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ ProximityPlacementGroupListResult }) = collect(keys(_property_map_ProximityPlacementGroupListResult))
Swagger.property_type(::Type{ ProximityPlacementGroupListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ProximityPlacementGroupListResult[name]))}
Swagger.field_name(::Type{ ProximityPlacementGroupListResult }, property_name::Symbol) =  _property_map_ProximityPlacementGroupListResult[property_name]

function check_required(o::ProximityPlacementGroupListResult)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ProximityPlacementGroupListResult }, name::Symbol, val)
end


text:
export RainbowLearner

"""
    RainbowLearner(;kwargs...)

See paper: [Rainbow: Combining Improvements in Deep Reinforcement Learning](https://arxiv.org/abs/1710.02298)

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get Q-values of a state.
- `target_approximator`::[`AbstractApproximator`](@ref): similar to `approximator`, but used to estimate the target (the next state).
- `loss_func`: the loss function.
- `Vₘₐₓ::Float32`: the maximum value of distribution.
- `Vₘᵢₙ::Float32`: the minimum value of distribution.
- `n_actions::Int`: number of possible actions.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=4`: the frequency of updating the `approximator`.
- `target_update_freq::Int=500`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `default_priority::Float32=1.0f2.`: the default priority for newly added transitions. It must be `>= 1`.
- `n_atoms::Int=51`: the number of buckets of the value function distribution.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `rng = Random.GLOBAL_RNG`
"""
mutable struct RainbowLearner{
    Tq<:AbstractApproximator,
    Tt<:AbstractApproximator,
    Tf,
    Ts,
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    sampler::NStepBatchSampler
    Vₘₐₓ::Float32
    Vₘᵢₙ::Float32
    n_actions::Int
    n_atoms::Int
    support::Ts
    delta_z::Float32
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float32
    β_priority::Float32
    rng::R
    loss::Float32
end

Flux.functor(x::RainbowLearner) =
    (Q = x.approximator, Qₜ = x.target_approximator, S = x.support),
    y -> begin
        x = @set x.approximator = y.Q
        x = @set x.target_approximator = y.Qₜ
        x = @set x.support = y.S
        x
    end

function RainbowLearner(;
    approximator,
    target_approximator,
    loss_func,
    Vₘₐₓ,
    Vₘᵢₙ,
    n_actions,
    n_atoms = 51,
    support = collect(range(Float32(-Vₘₐₓ), Float32(Vₘₐₓ), length = n_atoms)),
    stack_size = 4,
    delta_z = Float32(support[2] - support[1]),
    γ = 0.99,
    batch_size = 32,
    update_horizon = 1,
    min_replay_history = 32,
    update_freq = 1,
    target_update_freq = 500,
    update_step = 0,
    default_priority = 1.0f2,
    β_priority = 0.5f0,
    traces = SARTS,
    rng = Random.GLOBAL_RNG,
)
    default_priority >= 1.0f0 || error("default value must be >= 1.0f0")
    copyto!(approximator, target_approximator)  # force sync
    support = send_to_device(device(approximator), support)
    sampler = NStepBatchSampler{traces}(;
        γ = γ,
        n = update_horizon,
        stack_size = stack_size,
        batch_size = batch_size,
    )
    RainbowLearner(
        approximator,
        target_approximator,
        loss_func,
        sampler,
        Vₘₐₓ,
        Vₘᵢₙ,
        n_actions,
        n_atoms,
        support,
        delta_z,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        β_priority,
        rng,
        0.0f0,
    )
end

function (learner::RainbowLearner)(env)
    s = send_to_device(device(learner.approximator), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    logits = learner.approximator(s)
    q = learner.support .* softmax(reshape(logits, :, learner.n_actions))
    vec(sum(q, dims = 1)) |> send_to_host
end

function RLBase.update!(learner::RainbowLearner, batch::NamedTuple)
    Q = learner.approximator
    Qₜ = learner.target_approximator
    γ = learner.sampler.γ
    β = learner.β_priority
    loss_func = learner.loss_func
    n_atoms = learner.n_atoms
    n_actions = learner.n_actions
    support = learner.support
    delta_z = learner.delta_z
    update_horizon = learner.sampler.n
    batch_size = learner.sampler.batch_size
    D = device(Q)
    states = send_to_device(D, batch.state)
    rewards = send_to_device(D, batch.reward)
    terminals = send_to_device(D, batch.terminal)
    next_states = send_to_device(D, batch.next_state)

    actions = CartesianIndex.(batch.action, 1:batch_size)

    target_support =
        reshape(rewards, 1, :) .+
        (reshape(support, :, 1) * reshape((γ^update_horizon) .* (1 .- terminals), 1, :))

    next_logits = Qₜ(next_states)
    next_probs = reshape(softmax(reshape(next_logits, n_atoms, :)), n_atoms, n_actions, :)
    next_q = reshape(sum(support .* next_probs, dims = 1), n_actions, :)
    if haskey(batch, :next_legal_actions_mask)
        l′ = send_to_device(D, batch[:next_legal_actions_mask])
        next_q .+= ifelse.(l′, 0.0f0, typemin(Float32))
    end
    next_prob_select = select_best_probs(next_probs, next_q)

    target_distribution = project_distribution(
        target_support,
        next_prob_select,
        support,
        delta_z,
        learner.Vₘᵢₙ,
        learner.Vₘₐₓ,
    )

    is_use_PER = haskey(batch, :priority)  # is use Prioritized Experience Replay
    if is_use_PER
        updated_priorities = Vector{Float32}(undef, batch_size)
        weights = 1.0f0 ./ ((batch.priority .+ 1f-10) .^ β)
        weights ./= maximum(weights)
        weights = send_to_device(D, weights)
    end

    gs = gradient(Flux.params(Q)) do
        logits = reshape(Q(states), n_atoms, n_actions, :)
        select_logits = logits[:, actions]
        batch_losses = loss_func(select_logits, target_distribution)
        loss =
            is_use_PER ? dot(vec(weights), vec(batch_losses)) * 1 // batch_size :
            mean(batch_losses)
        ignore() do
            if is_use_PER
                updated_priorities .= send_to_host(vec((batch_losses .+ 1f-10) .^ β))
            end
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)

    is_use_PER ? updated_priorities : nothing
end

@inline function select_best_probs(probs, q)
    q_argmax = argmax(q, dims = 1)
    prob_select = @inbounds probs[:, q_argmax] # !!! without @inbounds it would be really slow
    reshape(prob_select, :, length(q_argmax))
end

function project_distribution(supports, weights, target_support, delta_z, vmin, vmax)
    batch_size, n_atoms = size(supports, 2), length(target_support)
    clampped_support = clamp.(supports, vmin, vmax)
    tiled_support = reshape(
        repeat(clampped_support; outer = (n_atoms, 1)),
        n_atoms,
        n_atoms,
        batch_size,
    )

    projection =
        clamp.(
            1 .- abs.(tiled_support .- reshape(target_support, 1, :)) ./ delta_z,
            0,
            1,
        ) .* reshape(weights, n_atoms, 1, batch_size)
    reshape(sum(projection, dims = 1), n_atoms, batch_size)
end


text:
using PhotonTransferCurves
using Test

@testset "PhotonTransferCurves.jl" begin
    # Write your tests here.
end


text:
using CovidData, DataFrames
using Test

@testset "CovidData.jl" begin
  df = CovidData.covidjhudata()
  @test isa(df, DataFrame)
  df = CovidData.statedata()
  @test isa(df, DataFrame)
  df = CovidData.statedata(fillmissingmobility=true)
  @test isa(df, DataFrame)
  df = CovidData.statedata(policies=:indicators)
  @test isa(df, DataFrame)
end


text:
# Autogenerated wrapper script for CImGui_jll for x86_64-w64-mingw32
export compile_commands, libcimgui, libcimgui_helper, libimgui

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Relative path to `compile_commands`
const compile_commands_splitpath = ["share", "compile_commands.json"]

# This will be filled out by __init__() for all products, as it must be done at runtime
compile_commands_path = ""

# compile_commands-specific global declaration
# This will be filled out by __init__()
compile_commands = ""

# Relative path to `libcimgui`
const libcimgui_splitpath = ["bin", "libcimgui.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libcimgui_path = ""

# libcimgui-specific global declaration
# This will be filled out by __init__()
libcimgui_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libcimgui = "libcimgui.dll"

# Relative path to `libcimgui_helper`
const libcimgui_helper_splitpath = ["bin", "libcimgui_helper.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libcimgui_helper_path = ""

# libcimgui_helper-specific global declaration
# This will be filled out by __init__()
libcimgui_helper_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libcimgui_helper = "libcimgui_helper.dll"

# Relative path to `libimgui`
const libimgui_splitpath = ["bin", "libimgui-cpp.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libimgui_path = ""

# libimgui-specific global declaration
# This will be filled out by __init__()
libimgui_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libimgui = "libimgui-cpp.dll"

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    global compile_commands_path = normpath(joinpath(artifact_dir, compile_commands_splitpath...))

    global compile_commands = compile_commands_path
    global libcimgui_path = normpath(joinpath(artifact_dir, libcimgui_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libcimgui_handle = dlopen(libcimgui_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libcimgui_path))

    global libcimgui_helper_path = normpath(joinpath(artifact_dir, libcimgui_helper_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libcimgui_helper_handle = dlopen(libcimgui_helper_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libcimgui_helper_path))

    global libimgui_path = normpath(joinpath(artifact_dir, libimgui_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libimgui_handle = dlopen(libimgui_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libimgui_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
function FindIntersectingPoint(X1,X2,Plane)
    D = length(X1);
    dX = X2-X1;
    Ind = find(dX.!=0);
    Ind = Ind[1];
    t = (-Plane[end]-sum(Plane[1:end-1].*X1))/sum(Plane[1:end-1].*dX);
    X = t.*dX+X1;
end


text:
#Test
function Solver_CentralPlanner(optimizer::Type{<:AbstractOptimizer},resultfile::String="")
#function solver_DeterministicEquivalent(optimizer::Type{<:AbstractOptimizer},resultfile::String="",load_st2::Type{<:AbstractArray{Int64,2}})
#Deterministic Equivalent------------------------------------------------------#
println("Loading inputs...")

(p_D,D,Τ,wind,wind_opt,Ns_H) = loadinputs("data")

(syscost_det,x_w_value,x_e_value)=Investment_OPF_original(optimizer,set_opt_thermalgenerators,set_opt_winds,set_thermalgenerators,set_winds,set_demands,set_nodes,set_nodes_ref,set_nodes_noref,set_scenarios,set_times,P,V,max_demand,R,p_D,D,γ,Τ,wind,wind_opt,Ns_H,links,F_max_dict,B_dict,MapG,MapG_opt,MapD,MapW,MapW_opt,tech_thermal,tech_thermal_opt,tech_wind,tech_wind_opt,capacity_per_unit,varcost,invcost,maxBuilds,ownership,capacity_existingunits,fixedcost,EmissionsRate,HeatRate,fuelprice,life,CRF_thermal_opt,CRF_wind_opt)

resultfile != "" && open(resultfile, truncate = true) do f

    println(f,"*****************OPF Solution*****************")
    println(f,"")
    println(f,"Objective System Cost: ",syscost_det)
    println(f,"")
    println(f,"Wind Investments: ",x_w_value)
    println(f,"")
    println(f,"Thermal Investments: ",x_e_value)
    end

return (syscost_det,x_w_value,x_e_value)
end


text:
module jdplyr

include("Functions.jl")
using DataFrames, Lazy, Statistics

for n in [names(DataFrames); names(Lazy); names(Statistics)]
        @eval export $n
end

end


text:
using Distributed
addprocs(16)
@everywhere include("optimize.jl")
#compute_acoustic_energy()
compute_Eac_along_path()

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
using PrettyTables

struct PolicyValueTestResults
    results::Vector{PolicyValueTestResult}
end

function PolicyValueTestResults(x::MuZeroConfig)
    getpolicytests()
    d = policy_testset(x.game_type)
    tests = []
    for (idx, kv) in enumerate(d)
        game, actions = kv
        push!(tests, PolicyTest(idx, game, actions))
    end
    PolicyValueTestResults(tests)
end

struct PolicyValueTestResultsSVG
    tests::PolicyValueTestResults
end

using Random

svg(x::PolicyValueTestResults) = PolicyValueTestResultsSVG(x)
Base.rand(x::PolicyValueTestResults) = rand(x.results)
Base.getindex(x::PolicyValueTestResults, idx::Integer) = x.results[idx]
Base.getindex(x::PolicyValueTestResults, idx::Vector{<:Integer}) = x.results[idx]
Base.length(x::PolicyValueTestResults) = length(x.results)
Base.lastindex(x::PolicyValueTestResults) = length(x)
Base.iterate(x::PolicyValueTestResults, state=1) = state > length(x) ? nothing : (x[state], state+1)

getdf(x::PolicyValueTestResults; range=1:length(x), limit=3) = getdf(PolicyValueTestResult, x.results[range])

function Base.show(io::IO, m::MIME"text/plain", x::PolicyValueTestResults)
    len = length(x.results)
    println()
    pretty_results(
        io,
        getdf(x);
        title="Policy Test Results (all $len)"
    )
end

function pretty_results(io::IO, df::DataFrame; title)
    pf_hl = Highlighter((data, i, j)->(j==2) && data[i,j] == "Pass", crayon"green");
    pretty_table(
        io,
        df;
        show_row_number=false,
        nosubheader=true,
        title=title,
        highlighters=(pf_hl),
    )
end

function Base.getindex(x::PolicyValueTestResults, r::UnitRange{<:Integer})
    println()
    pretty_results(
        stdout,
        getdf(x; range=r),
        title="Policy Test Results ($(r))"
    )
end


text:
using ClassicalOrthogonalPolynomials, QuasiArrays, InfiniteArrays
using BlockBandedMatrices, BlockArrays, LazyArrays, Test, FillArrays, InfiniteLinearAlgebra, ArrayLayouts, LazyBandedMatrices
import BlockBandedMatrices: _BlockBandedMatrix
import ClassicalOrthogonalPolynomials: Inclusion, ShuffledRFFT
import QuasiArrays: MulQuasiArray

@testset "Fourier" begin
    @testset "ShuffledRFFT" begin
        ret = randn(3,5)
        p = ShuffledRFFT{Float64}(size(ret,1))
        P = ShuffledRFFT{Float64}(size(ret),1)
        @test size(p) == (3,)
        @test size(P) == (3,5)
        for k = 1:size(ret,2)
            @test (P * ret)[:,k] ≈ p * ret[:,k]
        end

        p = ShuffledRFFT{Float64}(size(ret,2))
        P = ShuffledRFFT{Float64}(size(ret),2)
        for k = 1:size(ret,1)
            @test (P * ret)[k,:] ≈ p * ret[k,:]
        end
    end

    @testset "Evaluation" begin
        F = Fourier()

        @test F == F

        @test axes(F,2) isa BlockedUnitRange{InfiniteArrays.InfStepRange{Int,Int}}
        @test axes(F) isa Tuple{<:Inclusion,BlockedUnitRange{InfiniteArrays.InfStepRange{Int,Int}}}
        @test axes(F,2)[Block(2)] == 2:3
        @test F[0.1,1] == 1.0
        @test F[0.1,2] == sin(0.1)
        @test F[0.1,1:4] == [1,sin(0.1),cos(0.1),sin(2*0.1)]
        @test F[0.1,Block(4)] == [sin(3*0.1),cos(3*0.1)]
        @test F[0.1,Block.(1:3)] == [1,sin(0.1),cos(0.1),sin(2*0.1),cos(2*0.1)]

        u = F * PseudoBlockVector([[1,2,3]; zeros(∞)], (axes(F,2),));
        @test u[0.1] == 1 + 2sin(0.1) + 3cos(0.1)
    end

    @testset "Transform" begin
        F = Fourier()
        θ = axes(F,1)
        @test F[:,Base.OneTo(5)] \ cos.(θ) ≈ [0,0,1,0,0]
        @test F[:,Block.(Base.OneTo(5))] \ cos.(θ) ≈ [0,0,1,0,0,0,0,0,0]

        @test (F \ cos.(θ))[Block(2)] ≈ [0,1]
        u = F * (F \ exp.(cos.(θ)))
        @test u[0.1] ≈ exp(cos(0.1))
        @test F[:,Base.OneTo(5)] \ [cos.(θ) sin.(θ)] ≈ [0 0; 0 1; 1 0; 0 0; 0 0]
        @test F[:,Block.(Base.OneTo(3))] \ [cos.(θ) sin.(θ)] ≈ [0 0; 0 1; 1 0; 0 0; 0 0]
        U = F / F \ [exp.(cos.(θ)) cos.(cos.(θ))]
        @test U[0.1,:] ≈ [exp(cos(0.1)),cos(cos(0.1))]
        @test U[[0.1,0.2],:] ≈ [exp(cos(0.1)) cos(cos(0.1)); exp(cos(0.2)) cos(cos(0.2))]
        @test U[0.1,1] ≈ exp(cos(0.1))
        @test U[[0.1,0.2],1] ≈ exp.(cos.([0.1,0.2]))
    end

    @testset "Derivative" begin
        F = Fourier()
        @test F\F === Eye((axes(F,2),))
        @test (F'F)[1:10,1:10] == Diagonal(Vcat(2π,Fill(1.0π,9)))
        @test blockisequal(axes(F'F), (axes(F,2),axes(F,2)))
        D = Derivative(axes(F,1))
        D̃ = (D*F).args[2]
        @test (F\F)*D̃ isa BlockArray
        @test (F \ (D*F))[Block.(1:3),Block.(1:3)] == [0 0 0 0 0; 0 0.0 -1 0 0; 0 1 0 0 0; 0 0 0 0 -2; 0 0 0 2 0]

        u = F * PseudoBlockVector([[1,2,3,4,5]; zeros(∞)], (axes(F,2),));
        @test blockisequal(axes(D̃,2),axes(u.args[2],1))
        @test (D*u)[0.1] ≈ 2cos(0.1) - 3sin(0.1) + 8cos(2*0.1) - 10sin(2*0.1)
    end

    @testset "cos(θ)" begin
        F = Fourier()
        θ = axes(F,1)
        c,s = cos.(θ),sin.(θ)
        @test c .* F isa MulQuasiArray
        @test s .* F isa MulQuasiArray
        @test cos.(θ) .* F isa MulQuasiArray
        @test sin.(θ) .* F isa MulQuasiArray
        X,Y = F \ (c .* F),F \ (s .* F)
        N = 10
        XY = X*Y
        @test XY isa MulMatrix
        @test MemoryLayout(XY) isa LazyBandedMatrices.ApplyBlockBandedLayout{typeof(*)}
        @test MemoryLayout(view(XY,Block.(1:N),Block.(1:N))) isa LazyBandedMatrices.ApplyBlockBandedLayout{typeof(*)}
        @test XY[Block.(1:N),Block.(1:N)] isa BlockSkylineMatrix
        @test blockbandwidths(XY[Block.(1:N),Block.(1:N)]) == (2,2)
        @test XY[Block.(1:N),Block.(1:N)] ≈ ApplyMatrix(*,Y,X)[Block.(1:N),Block.(1:N)]
        @test ApplyMatrix(*,X,X)[Block.(1:N),Block.(1:N)] + ApplyMatrix(*,Y,Y)[Block.(1:N),Block.(1:N)] ≈ Eye(19)
    end
end


text:
mutable struct Image
    pixels::Dict{CartesianIndex{2}, Bool}
    infinite_pixel::Bool
    Image(pixels) = new(pixels, false)
end

function enhance!(image::Image, algo::BitVector)
    moves = CartesianIndex.((
        (-1, -1), (-1, 0), (-1, 1),
        ( 0, -1), ( 0, 0), ( 0, 1),
        ( 1, -1), ( 1, 0), ( 1, 1)
    ))
    min_index = minimum(keys(image.pixels)) + first(moves)
    max_index = maximum(keys(image.pixels)) + last(moves)
    enhanced_pixels = empty(image.pixels)
    for index in min_index:max_index
        bits = [get(image.pixels, index + move, image.infinite_pixel) for move in moves]
        enhanced_pixels[index] = algo[convert(Int, bits) + 1]
    end
    image.pixels = enhanced_pixels
    image.infinite_pixel = image.infinite_pixel ? last(algo) : first(algo)
    image
end

Base.convert(::Type{Int}, bits::Vector{Bool}) = foldl((x, y) -> (x << 1) | y, bits)

function Base.parse(::Type{Image}, lines::Vector{String})
    bits = permutedims(hcat(collect.(lines)...)) .== '#'
    pixels = Dict(zip(CartesianIndices(bits), bits))
    Image(pixels)
end

function part1(input)
    algo = collect(readline(input)) .== '#'
    image = parse(Image, readlines(input)[3:end])
    for _ in 1:2
        enhance!(image, algo)
    end
    sum(values(image.pixels))
end

function part2(input)
    algo = collect(readline(input)) .== '#'
    image = parse(Image, readlines(input)[3:end])
    for _ in 1:50
        enhance!(image, algo)
    end
    sum(values(image.pixels))
end

@show part1("day20/image.txt")

@show part2("day20/image.txt")


text:
@muladd function perform_step!(integrator,cache::SplitEMConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  u = dt*(integrator.f.f1(uprev,p,t) +
           integrator.f.f2(uprev,p,t)) +
           integrator.g(uprev,p,t).*W.dW + uprev
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SplitEMCache,f=integrator.f)
  @unpack rtmp1,rtmp2 = cache
  @unpack t,dt,uprev,u,W,p = integrator

  integrator.g(rtmp2,uprev,p,t)
  if is_diagonal_noise(integrator.sol.prob)
    rmul!(rtmp2,W.dW) # rtmp2 === rtmp3
    @. u = uprev + rtmp2
  else
    mul!(rtmp1,rtmp2,W.dW)
    @. u = uprev + rtmp1
  end

  integrator.f.f1(t,uprev,rtmp1)
  @. u += dt*rtmp1
  integrator.f.f2(t,uprev,rtmp1)
  @. u += dt*rtmp1
end


text:
using Documenter, AVExamples
makedocs(modules=[AVExamples],
        doctest=false, clean=true,
        format =:html,
        authors="Huckleberry Febbo",
        sitename="MAVs",
        pages = Any[
        "Home" => "index.md",
        "Packages"=>Any[
             "packages/nloptcontrol_planner/index.md"
             "packages/obstacle_detector/index.md"
             "packages/vehicle_description/index.md"
             "packages/ros_chrono/index.md"
             "packages/mavs_ros_planner/index.md"
             "packages/contact_sensor/index.md"
             "System"=>Any[
                "packages/system/system_shutdown/index.md"
                "packages/system/logging/index.md"
             ]
                           ],
        "System Demos"=>Any[
              "demos/system/demoA.md"
              "demos/system/demoB.md"
              "demos/system/demoC.md"
              "demos/system/demoD.md"
              "demos/system/demoE.md"
              "demos/system/demoF.md"
                           ],
        "Package Demos"=>Any[
             "demos/vehicle_description/demo.md"
             "demos/nloptcontrol_planner/demo.md"
             "demos/ros_chrono/demo.md"
                            ],
         "Miscellaneous"=>Any[
               "issues/index.md"
                             ]
               ]
               )
deploydocs(
    deps=Deps.pip("mkdocs","python-markdown-math"),
    repo="github.com/JuliaMPC/AVExamples.jl.git",
    target="build",
    osname="linux",
    julia="0.6",
    make=nothing)


text:
using SpecialMatrices
using EllipsisNotation
using DiffEqBase, OrdinaryDiffEq, StochasticDiffEq
using Plots; pyplot()

############ Spatial Parameters

const x_start = -100
const x_end   = 400
const y_start = 0
const y_end   = 100
const dx = 5
const dy = 5
const N = Int((x_end-x_start) / dx + 1)
const M = Int((y_end-y_start) / dy +1)
const X = repmat(collect(x_start:dx:x_end)',M,1)
const Y = repmat(collect(y_start:dy:y_end),1,N)

############ RA Parameters

const DRA = 25.46
const β = .017 # control steepness, higher is steeper
const kA = 0.002
const kmax = 1
const n = 2
const xf = 400

const kp = 0.1
const gamma = 100
const mon = 3.0
const moff = 0.0013
const rdeg1 = 0.0001
const rdeg2 = 0.0001
const jalpha = .85
const jbeta = .850
const bpdeg1 = 0.01
const d = 0.1
const e = 1
const Vr = 2e-2

############ Setup Diffusion Matrices

const Ax = -full(Strang(N))
const Ay = -full(Strang(M))
Ax[2,1] = 2
Ax[end-1,end] = 2
Ax[end,end] = -2*(1+dx*kA)
Ay[1,2] = 2
Ay[end,end-1] = 2
Ax ./= dx^2
Ay ./= dy^2

const Dxtmp = zeros(M,N)
const Dytmp = zeros(M,N)

############# Spatial Constants

const Vmax = 1.0
V(x) = ifelse(x>xf-40,Vmax,zero(typeof(Vmax)))
const VRA = V.(X)
no_cyp(x) = ifelse(x<0 || (xf-40<x && x<=xf),true,false)
const NCYP_kmax = no_cyp.(X)

############# Hox-Krox Parameters

const cₕ = 7.5
const cₖ = 3.0
const κₕ = .4
const κₖ = 4.0
const dₕ = 0.4
const dₖ = 0.4
const nₕ = 2
const nₖ = 2
const m = 2

function ra_gradient(t,u,p,du) # no matrix cache for autodiff, better for stiff
  RAout = view(u,:,:,1)
  RAin = view(u,:,:,2)
  R = view(u,:,:,3)
  RAR = view(u,:,:,4)
  BP = view(u,:,:,5)
  RABP = view(u,:,:,6)
  dRAout = view(du,:,:,1)
  dRAin = view(du,:,:,2)
  dR = view(du,:,:,3)
  dRAR = view(du,:,:,4)
  dBP = view(du,:,:,5)
  dRABP = view(du,:,:,6)
  Vbp, kdeg = p[1],p[2]

  # Doesn't matter because all the time is in the factorization anyways
  Dxtmp = RAout*Ax
  Dytmp = Ay*RAout

  Threads.@threads for i in eachindex(dRAout)
    @inbounds diffRA     = DRA*(Dxtmp[i] + Dytmp[i])
    @inbounds dRAout[i]  = VRA[i] - β*RAout[i] + kp*RAin[i] + diffRA
    @inbounds dRAin[i]   = β*RAout[i] - kp*RAin[i] - kdeg*(RAR[i]/(gamma+RAR[i]))*RAin[i] - mon*RAin[i]*BP[i] + moff*RABP[i] - rdeg1*RAin[i]
    @inbounds dR[i]      = Vr - rdeg2*R[i]- jalpha*RABP[i]*R[i] + jbeta*BP[i]*RAR[i]
    @inbounds dRAR[i]    = jalpha*RABP[i]*R[i] - jbeta*BP[i]*RAR[i]
    @inbounds dBP[i]     = Vbp - bpdeg1*BP[i] - mon*RAin[i]*BP[i] + moff*RABP[i] + jalpha*RABP[i]*R[i] - jbeta*BP[i]*RAR[i] + ((d*RAR[i])/(e+RAR[i]))
    @inbounds dRABP[i]   = mon*RAin[i]*BP[i] - moff*RABP[i] - jalpha*RABP[i]*R[i] + jbeta*BP[i]*RAR[i]
  end
  nothing
end

function ra_gradient2(t,u,p,du) # better for non-stiff
  RAout = view(u,:,:,1)
  RAin = view(u,:,:,2)
  R = view(u,:,:,3)
  RAR = view(u,:,:,4)
  BP = view(u,:,:,5)
  RABP = view(u,:,:,6)
  dRAout = view(du,:,:,1)
  dRAin = view(du,:,:,2)
  dR = view(du,:,:,3)
  dRAR = view(du,:,:,4)
  dBP = view(du,:,:,5)
  dRABP = view(du,:,:,6)
  Vbp, kdeg = p[1],p[2]

  A_mul_B!(Dxtmp,RAout,Ax)
  A_mul_B!(Dytmp,Ay,RAout)

  Threads.@threads for i in eachindex(dRAout)
    @inbounds diffRA     = DRA*(Dxtmp[i] + Dytmp[i])
    @inbounds dRAout[i]  = VRA[i] - β*RAout[i] + kp*RAin[i] + diffRA
    @inbounds dRAin[i]   = β*RAout[i] - kp*RAin[i] - kdeg*(RAR[i]/(gamma+RAR[i]))*RAin[i] - mon*RAin[i]*BP[i] + moff*RABP[i] - rdeg1*RAin[i]
    @inbounds dR[i]      = Vr - rdeg2*R[i]- jalpha*RABP[i]*R[i] + jbeta*BP[i]*RAR[i]
    @inbounds dRAR[i]    = jalpha*RABP[i]*R[i] - jbeta*BP[i]*RAR[i]
    @inbounds dBP[i]     = Vbp - bpdeg1*BP[i] - mon*RAin[i]*BP[i] + moff*RABP[i] + jalpha*RABP[i]*R[i] - jbeta*BP[i]*RAR[i] + ((d*RAR[i])/(e+RAR[i]))
    @inbounds dRABP[i]   = mon*RAin[i]*BP[i] - moff*RABP[i] - jalpha*RABP[i]*R[i] + jbeta*BP[i]*RAR[i]
  end
  nothing
end

function ra_noise(t,u,ϵ,du)
  RAout = view(u,:,:,1)
  RAin = view(u,:,:,2)
  RAR = view(u,:,:,4)
  dRAout = view(du,:,:,1)
  dRAin = view(du,:,:,2)
  dRAR = view(du,:,:,4)

  Threads.@threads for i in eachindex(dRAin)
    @inbounds dRAout[i] = ϵ * RAout[i]
    @inbounds dRAin[i] = ϵ * RAin[i]
    @inbounds dRAR[i] = ϵ * RAR[i]
  end
  nothing
end

function hox_krox(du,u,p,t)
  RAR = view(u,:,:,4)
  gₕ = view(u,:,:,7)
  gₖ = view(u,:,:,8)
  dgₕ = view(du,:,:,7)
  dgₖ = view(du,:,:,8)

  Threads.@threads for i in eachindex(dgₕ)
    @inbounds dgₕ[i] = (cₕ*gₕ[i]^nₕ + (κₕ*RAR[i])^m) / (1 + cₕ*gₕ[i]^nₕ + cₖ*gₖ[i]^nₖ + (κₕ*RAR[i])^m) - dₕ*gₕ[i]
    @inbounds dgₖ[i] = (cₖ*gₖ[i]^nₖ + (κₖ*RAR[i])^m) / (1 + cₕ*gₕ[i]^nₕ + cₖ*gₖ[i]^nₖ + (κₖ*RAR[i])^m) - dₖ*gₖ[i]
  end
  nothing
end

function hox_krox_noise(t,u,ϵₕₖ,du)
  RAout = view(u,:,:,1)
  RAin = view(u,:,:,2)
  RAR = view(u,:,:,4)
  gₕ = view(u,:,:,7)
  gₖ = view(u,:,:,8)
  dRAout = view(du,:,:,1)
  dRAin = view(du,:,:,2)
  dRAR = view(du,:,:,4)
  dgₕ = view(du,:,:,7)
  dgₖ = view(du,:,:,8)

  Threads.@threads for i in eachindex(dRAin)
    @inbounds dgₕ[i] = ϵₕₖ * gₕ[i]
    @inbounds dgₖ[i] = ϵₕₖ * gₖ[i]
  end
  nothing
end

function boundary_sharpening(;Vbp=100.0,kdeg=1e5,ϵ=0.03,ϵₕₖ=0.20)

  ################ Solve ODE

  println("Solve the ODE")
  tspan = (0.0,1e9)
  u0 = zeros(M,N,6)
  f = ParameterizedFunction(ra_gradient,(Vbp,kdeg))
  prob = ODEProblem(f,u0,tspan)
  @time ra_sol = solve(prob,Rosenbrock23(),save_everystep=false,dt=0.5,
                       progress=true,abstol=1e-3,reltol=1e-1,progress_steps=1)

  u = copy(ra_sol[end])
  du = similar(u)
  t = 0
  @time ra_gradient(t,u,(Vbp,kdeg),du)
  println("The maximum du is: $(maximum(du))")

  ################# Add noise to the gradient

  println("Solve the gradient SDE")

  u0 = ra_sol[end]
  tspan2 = (0.0,500.0)
  f2 = ParameterizedFunction(ra_gradient2,(Vbp,kdeg))
  g = ParameterizedFunction(ra_noise,ϵ)
  prob = SDEProblem(f2,g,u0,tspan2)
  @time noisy_ra_sol = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),
                             save_everystep=false,progress_steps=10_000,
                             progress=true,abstol=1e-2,reltol=1e-2)

  ################# Hox-Krox to SS

  u0 = 0.16ones(M,N,8) + 0.20rand(M,N,8) + 0.0005
  u0[:,:,1] .= noisy_ra_sol[end][:,:,1]
  u0[:,:,2] .= noisy_ra_sol[end][:,:,2]
  u0[:,:,3] .= noisy_ra_sol[end][:,:,3]
  u0[:,:,4] .= noisy_ra_sol[end][:,:,4]
  u0[:,:,5] .= noisy_ra_sol[end][:,:,5]
  u0[:,:,6] .= noisy_ra_sol[end][:,:,6]
  u0[:,1:46,7] = ones(X[:,1:46])
  u0[:,:,8] .*= 0.0

  println("Solve the Hox-Krox ODE")

  tspan = (0.0,5000.0)
  prob = ODEProblem(hox_krox,u0,tspan)

  @time sol = solve(prob,Tsit5(),save_everystep=false,progress=true,
              abstol=1e-3,reltol=1e-2,progress_steps=100)

  ########################### Hox Krox Noise Sharpening

  println("Solve the Hox-Krox SDE")

  ghk = ParameterizedFunction(hox_krox_noise,ϵₕₖ)
  prob = SDEProblem(hox_krox,ghk,sol[end],(0.0,10000.0))
  @time sol2 = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),
        progress_steps=1000,saveat=100.0,
        progress=true,abstol=1e-2,reltol=1e-2)

  ############################ Return results

  sol2
end

max_stars_repo_name: ChrisRackauckas/RAMeanVar

max_issues_repo_name: ChrisRackauckas/RAMeanVar

max_forks_repo_name: ChrisRackauckas/RAMeanVar


text:
function getcenter(e::Element)
    dim = typeof(e.nodes[1]).parameters[1]
    T = typeof(e.nodes[1]).parameters[2]
    S = (n -> n.x).(e.nodes)
    v = sum(S) / length(S)
    return Node{dim,T}(v)
end

function substitute_nodes!(newnodes::Array{N}, nodes::Vector{N}) where N <: Node
    for i in 1:(length(newnodes))
        exnodes = findall(n -> isapprox(n.x, newnodes[i].x), nodes)
        if (length(exnodes)) == 1 newnodes[i] = nodes[exnodes[1]] end
        @debug if (length(exnodes)) > 1 error("node in nodelist several times") end
    end
    return newnodes
end

function findlowestancestorwithnodes(e::Element, n::Vector) # nodevector
    if issubset(n, e.nodes) return e end
    if e.parent === nothing return nothing end
    return findlowestancestorwithnodes(e.parent, n)
end

function have_identical_leaves(e::Element, f::Element)
    le = collect_children(e)
    fe = collect_children(f)
    return issetequal(le, fe)
end

function make_master!(e::Element)
    if length(e.attachedto) == 0 return e end
    @debug if (length(e.children) != 0) error("trying to make a non-leaf into a master") end
    oldmaster = get_master(e)
    if oldmaster == e return e end
    @debug if length(oldmaster.children) != 0 error("oldmaster has children") end
    @debug if !issubset(vertices(oldmaster), vertices(e)) error("nodes of old and new masters must coincide") end

    e.attachedto = copy(oldmaster.attachedto) # better with union and filter
    push!(e.attachedto, oldmaster)
    deleteat!(e.attachedto, findall(x -> have_identical_leaves(x, e), e.attachedto))
    e.attachedto = unique!(e.attachedto)
    e.ismaster = true
    slaveleaves = collect_children(e.attachedto)
    (x -> x.attachedto = [e]).(slaveleaves) # collate this?
    (x -> x.ismaster = false).(slaveleaves) # should suffice on oldmaster
    @debug if !e.ismaster error("somehow e seems to be one of its own slaveleaves") end
    return e
end

function descend_until_split(e::Element)
    if length(e.children) != 1 return e end
    return descend_until_split(e.children[1])
end

function attach_children!(newels::Vector{Element{dim,N,M}}, e::Element) where {dim,N,M} # used on face-level
    @debug if length(e.children) != 0 error("already has children") end
    @debug if length(newels) == 0 error("nothing to attach") end

    if (e.ismaster == false)
        (x -> x.attachedto = [get_master(e)]).(newels)
        (x -> x.ismaster = false).(newels)
        e.children = newels
        (x -> x.parent = e).(newels)
        return [get_master(newels[1])]
    end

    altmasters = findall(x -> length(x.children) == 0, descend_until_split.(e.attachedto))
    if length(altmasters) != 0
        newmaster = descend_until_split(e.attachedto[altmasters[1]])
        make_master!(newmaster)
        e.children = newels
        (x -> x.parent = e).(newels)
        (x -> x.attachedto = [newmaster]).(newels)
        (x -> x.ismaster = false).(newels)
        return [get_master(newels[1])]
    end

    e.children = newels
    (x -> x.parent = e).(newels)
    (x -> x.ismaster = true).(newels)

    for m in newels
        for oldslaves in e.attachedto
            for s in descend_until_split(oldslaves).children
                if vertices(s) ⊆ vertices(m)
                    (x -> x.attachedto = [m]).(collect_children(s))
                    push!(m.attachedto, s)
                end
            end
        end
    end

    newelattachments = collect(Iterators.flatten((x -> x.attachedto).(newels)))
    if !(collect_children(e.attachedto) ⊆ collect_children(newelattachments))
        error("not all of the attachments could be propagated the children of e")
    end

    return get_master.(newels)
end

function unattach_children!(e::Element)
    if length(e.children) == 0 return e end
    che = collect_children(e)
    toobigattachments = findall(x -> nothing == findlowestancestorwithnodes(x, e.nodes), get_master.(che))
    if length(toobigattachments) != 0
        e.children = []
        e.ismaster = false
        e.attachedto = [get_master(che[toobigattachments[1]])]
        return e
    end

    attach = Vector{Element}()
    for c in che
        if c.ismaster
            attach = vcat(attach, (c.attachedto))
        else
            attach = vcat(attach, collect_children(get_master(c).attachedto))
            attach = push!(attach, get_master(c))
        end
    end
    attach = (x -> findlowestancestorwithnodes(x, e.nodes)).(attach)
    attach = filter!(x -> length(intersect(collect_children(e), collect_children(x))) == 0, attach)
    e.children = []
    e.ismaster = true
    e.attachedto = unique!(convert(Vector{Element}, attach))
    (x -> x.ismaster = false).(collect_children(e.attachedto)) # elegantly collate this?
    (x -> x.attachedto = [e]).(collect_children(e.attachedto))
    get_master(e)
    return e
end

function derefine!(e::Element)
    @debug if length(e.children) == 0 error("can not derefine leaf") end
    for x in e.children
        if length(x.children) != 0 error("can not derefine cell whose children are not leaves") end
    end
    cs = e.children
    e.children = []
    for b in e.boundaries
        unattach_children!(b)
    end
    return e
end

function isotropic_split!(e::Element{2,4,4})
    return split!(e, [1,2])
end

function split!(e::Element{2,4,4}, axes::Vector{Int})
    if length(e.children) != 0 error("can not split non-leaf") end
    xstep = 1
    ystep = 1
    if 1 ∈ axes xstep = 2 end
    if 2 ∈ axes ystep = 2 end
    if xstep * ystep == 1 error("come on, don't split something into 1 piece, it's embarassing") end
    dim = typeof(e.nodes[1]).parameters[1]
    T = typeof(e.nodes[1]).parameters[2]
    nodemat = Matrix{Node{dim,T}}(undef, 3, 3)
    nodemat[1,1] = e.nodes[1]
    nodemat[3,1] = e.nodes[2]
    nodemat[3,3] = e.nodes[3]
    nodemat[1,3] = e.nodes[4]
    nodemat[2,1] = getcenter(e.boundaries[1])
    nodemat[3,2] = getcenter(e.boundaries[2])
    nodemat[2,3] = getcenter(e.boundaries[3])
    nodemat[1,2] = getcenter(e.boundaries[4])
    nodemat[2,2] = getcenter(e)
    nodesnearby = unique!(collect(Iterators.flatten(getcontainednodes.(get_master.(e.boundaries)))))
    nodemat = substitute_nodes!(nodemat, nodesnearby)
    cellmat = Matrix{Element{2,4,4}}(undef, (xstep, ystep))
    for i in 1:xstep
        for j in 1:ystep
            cellmat[i,j] = Element([nodemat[i,j],nodemat[i + (2 ÷ xstep),j],nodemat[i + (2 ÷ xstep),j + (2 ÷ ystep)],nodemat[i,j + (2 ÷ ystep)]], e)
        end
    end

    for i in 1:xstep
        for j in 1:ystep
            if i < xstep attach_trivially!(cellmat[i,j], cellmat[i + 1,j]) end
            if j < ystep attach_trivially!(cellmat[i,j], cellmat[i,j + 1]) end
        end
    end

    allfaces = collect(Iterators.flatten((x -> x.boundaries).(cellmat)))
    for i in 1:length(e.boundaries)
        side = filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[i] == x), allfaces)
        attach_children!(convert(Array{Element{1,2,2}}, side), e.boundaries[i])
    end

    e.children = reshape(cellmat, length(cellmat)) # [cellmat[1,1],cellmat[2,1],cellmat[2,2],cellmat[1,2]]
    return e
end

function isotropic_split!(e::Element{3,8,6})
    return split!(e, [1,2,3])
end

function split!(e::Element{3,8,6}, axes::Vector{Int})
    if length(e.children) != 0 error("can not split non-leaf") end
    xstep = 1
    ystep = 1
    zstep = 1
    if 1 ∈ axes xstep = 2 end
    if 2 ∈ axes ystep = 2 end
    if 3 ∈ axes zstep = 2 end
    if xstep * ystep * zstep == 1 error("come on, don't split something into 1 piece, it's embarassing") end
    dim = typeof(e.nodes[1]).parameters[1]
    T = typeof(e.nodes[1]).parameters[2]
    nodemat = Array{Node{dim,T}}(undef, (3, 3, 3))
    nodemat[1,1,1] = e.nodes[1]
    nodemat[3,1,1] = e.nodes[2]
    nodemat[3,3,1] = e.nodes[3]
    nodemat[1,3,1] = e.nodes[4]
    nodemat[1,1,3] = e.nodes[5]
    nodemat[3,1,3] = e.nodes[6]
    nodemat[3,3,3] = e.nodes[7]
    nodemat[1,3,3] = e.nodes[8]
    nodemat[2,1,1] = getcenter(e.boundaries[1].boundaries[4])
    nodemat[3,2,1] = getcenter(e.boundaries[1].boundaries[3])
    nodemat[2,3,1] = getcenter(e.boundaries[1].boundaries[2])
    nodemat[1,2,1] = getcenter(e.boundaries[1].boundaries[1])
    nodemat[2,2,1] = getcenter(e.boundaries[1])
    nodemat[2,1,3] = getcenter(e.boundaries[6].boundaries[1])
    nodemat[3,2,3] = getcenter(e.boundaries[6].boundaries[2])
    nodemat[2,3,3] = getcenter(e.boundaries[6].boundaries[3])
    nodemat[1,2,3] = getcenter(e.boundaries[6].boundaries[4])
    nodemat[2,2,3] = getcenter(e.boundaries[6])
    nodemat[3,1,2] = getcenter(e.boundaries[2].boundaries[2])
    nodemat[1,1,2] = getcenter(e.boundaries[2].boundaries[4])
    nodemat[2,1,2] = getcenter(e.boundaries[2])
    nodemat[1,3,2] = getcenter(e.boundaries[4].boundaries[2])
    nodemat[3,3,2] = getcenter(e.boundaries[4].boundaries[4])
    nodemat[2,3,2] = getcenter(e.boundaries[4])
    nodemat[3,2,2] = getcenter(e.boundaries[3])
    nodemat[1,2,2] = getcenter(e.boundaries[5])
    nodemat[2,2,2] = getcenter(e)
    nodesnearby = getnodeneighborhood(e)
    nodemat = substitute_nodes!(nodemat, nodesnearby)
    cellmat = Array{Element{3,8,6}}(undef, (xstep, ystep, zstep))
    for i in 1:xstep
        for j in 1:ystep
            for k in 1:zstep
                cellmat[i,j,k] = Element([nodemat[i,j,k],nodemat[i + (2 ÷ xstep),j,k],nodemat[i + (2 ÷ xstep),j + (2 ÷ ystep),k],nodemat[i,j + (2 ÷ ystep),k],
                nodemat[i,j,k + (2 ÷ zstep)],nodemat[i + (2 ÷ xstep),j,k + (2 ÷ zstep)],nodemat[i + (2 ÷ xstep),j + (2 ÷ ystep),k + (2 ÷ zstep)],nodemat[i,j + (2 ÷ ystep),k + (2 ÷ zstep)]], e)
            end
        end
    end
    for i in 1:xstep
        for j in 1:ystep
            for k in 1:zstep
                if i < xstep attach_trivially!(cellmat[i,j,k], cellmat[i + 1,j,k]) end
                if j < ystep attach_trivially!(cellmat[i,j,k], cellmat[i,j + 1,k]) end
                if k < zstep attach_trivially!(cellmat[i,j,k], cellmat[i,j,k + 1]) end
            end
        end
    end
    allfaces = collect(Iterators.flatten((x -> x.boundaries).(cellmat)))
    for i in 1:length(e.boundaries)
        side = filter(x -> (length(x.attachedto) == 0 && x.boundaryof.boundaries[i] == x), allfaces)
        attach_children!(convert(Array{Element{2,4,4}}, side), e.boundaries[i])
    end

    e.children = reshape(cellmat, length(cellmat))
    # [cellmat[1,1,1],cellmat[2,1,1],cellmat[2,2,1],cellmat[1,2,1],
    # cellmat[1,1,2],cellmat[2,1,2],cellmat[2,2,2],cellmat[1,2,2]]

    return e
end


text:
open("questions-words.txt") do f
    enter = false
    leave = false
    global r = []
    for line in readlines(f)
        words = split(line)
        if leave
            break
        end
        if words[1] != ":" && !enter
            continue
        end
        if words[1] == ":" && words[2] == "family"
            enter = true
            continue
        end
        if enter && !leave && words[1] == ":"
            leave = true
        end
        if enter && !leave
            push!(r, line)
        end
    end

    for v in r
        println(v)
    end
end


text:
function create_system_data(; with_time_series=false, time_series_in_memory=false)
    data = IS.SystemData(; time_series_in_memory=time_series_in_memory)

    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(data, component)

    if with_time_series
        file = joinpath(FORECASTS_DIR, "ComponentsAsColumnsNoTime.json")
        IS.add_time_series_from_file_metadata!(
            data,
            IS.InfrastructureSystemsComponent,
            file,
        )
        time_series = get_all_time_series(data)
        IS.@assert_op length(time_series) > 0
    end

    return data
end

function create_system_data_shared_time_series(; time_series_in_memory=false)
    data = IS.SystemData(; time_series_in_memory=time_series_in_memory)

    name1 = "Component1"
    name2 = "Component2"
    component1 = IS.TestComponent(name1, 5)
    component2 = IS.TestComponent(name2, 6)
    IS.add_component!(data, component1)
    IS.add_component!(data, component2)

    ts = IS.SingleTimeSeries(name="val", data=create_time_array())
    IS.add_time_series!(data, component1, ts)
    IS.add_time_series!(data, component2, ts)

    return data
end

function get_all_time_series(data)
    return collect(IS.get_time_series_multiple(data))
end

function create_time_array()
    dates = collect(
        Dates.DateTime("1/1/2020 00:00:00", "d/m/y H:M:S"):Dates.Hour(1):Dates.DateTime(
            "1/1/2020 23:00:00",
            "d/m/y H:M:S",
        ),
    )
    data = collect(1:24)
    component_name = "gen"
    return TimeSeries.TimeArray(dates, data, [component_name])
end

function create_dates(start_time::Dates.DateTime, resolution, end_time::Dates.DateTime)
    return collect(start_time:resolution:end_time)
end

function create_dates(start_time::String, resolution, end_time::String)
    return create_dates(Dates.DateTime(start_time), resolution, Dates.DateTime(end_time))
end

"""
Verifies that printing an object doesn't crash, which has happened several times.
"""
function verify_show(obj)
    io = IOBuffer()
    show(io, "text/plain", obj)
    val = String(take!(io))
    @test !isempty(val)
    return
end


text:
module TestVI
using GaussianProcesses, Distributions, PDMats, Plots, LinearAlgebra
gr()
using Test, Random
Random.seed!(203617)

@testset "VI" begin
    n = 20
    X = collect(range(-3,stop=3,length=n));
    f = 2*cos.(2*X);
    Y = [rand(Poisson(exp.(f[i]))) for i in 1:n];

    k = Matern(3/2,0.0,0.0)   # Matern 3/2 kernel
    l = PoisLik()             # Poisson likelihood
    gp = GP(X, vec(Y), MeanZero(), k, l)
    set_priors!(gp.kernel,[Normal(-2.0,4.0),Normal(-2.0,4.0)])
    Q, Ω, K = initialise_Q(gp)

    @testset "AutoDiff" begin
        exact = -0.5*exp(Q.m[end] + Q.V[end, end]/2)
        zyg_calc = dv_var_exp(l, Y[end], Q.m[end], Q.V[end, end])
        @test exact ≈ zyg_calc # Test AutoDiff is giving the same results as exact gradient computation
    end

    @testset "Poisson" begin
        Q = vi(gp, nits=10)
        println("Positive Definite: ", isposdef(Q.V))
        println("-------------------")

        xtest = range(minimum(gp.x),stop=maximum(gp.x),length=50)

        nsamps = 500
        ymean = [];
        visamples = Array{Float64}(undef, nsamps, size(xtest, 1))

        for i in 1:nsamps
            visamples[i, :] = rand(gp, xtest, Q)
            push!(ymean, predict_y(gp, xtest)[1])
        end

        q10 = [quantile(visamples[:,i], 0.1) for i in 1:length(xtest)]
        q50 = [quantile(visamples[:,i], 0.5) for i in 1:length(xtest)]
        q90 = [quantile(visamples[:,i], 0.9) for i in 1:length(xtest)]
        # plot(xtest, exp.(q50), ribbon=(exp.(q10), exp.(q90)), leg=true, fmt=:png, label="quantiles")
        plot(xtest, mean(ymean), label="posterior mean", w=2)
        xx = range(-3,stop=3,length=1000);
        f_xx = 2*cos.(2*xx);
        plot!(xx, exp.(f_xx), label="truth")
        scatter!(X,Y, label="data")
        savefig("vi_gp.png")
    end

    # @testset "Gaussian" begin
    #     n=10;                          #number of training points
    #     x = 2π * rand(n);              #predictors
    #     y = sin.(x) + 0.05*randn(n);   #regressors
    #     mZero = MeanZero()                   #Zero mean function
    #     kern = SE(0.0,0.0)                   #Sqaured exponential kernel (note that hyperparameters are on the log scale)
    #     l = GaussLik(0.1)             # Poisson likelihood

    #     logObsNoise = -1.0                        # log standard deviation of observation noise (this is optional)
    #     gp = GP(x, y, mZero, kern, l)       #Fit the GP

    #     Q = vi(gp)

    #     xtest = range(minimum(gp.x),stop=maximum(gp.x),length=50)

    #     nsamps = 500
    #     ymean = [];
    #     visamples = Array{Float64}(undef, nsamps, size(xtest, 1))

    #     for i in 1:nsamps
    #         visamples[i, :] = rand(gp, xtest, Q)
    #         push!(ymean, predict_y(gp, xtest)[1])
    #     end

    #     q10 = [quantile(visamples[i,:], 0.1) for i in 1:length(xtest)]
    #     q50 = [quantile(visamples[:,i], 0.5) for i in 1:length(xtest)]
    #     q90 = [quantile(visamples[:,i], 0.9) for i in 1:length(xtest)]
    #     plot(xtest, exp.(q50), ribbon=(exp.(q10), exp.(q90)), leg=true, fmt=:png, label="quantiles")
    #     plot!(xtest, mean(ymean), label="posterior mean", w=2)
    #     xx = range(-3,stop=3,length=1000);
    #     f_xx = 2*cos.(2*xx);
    #     plot!(xx, exp.(f_xx), label="truth")
    #     scatter!(X,Y, label="data")
    #     savefig("vi_gp.png")
    # end
 end
 end


text:
# Autogenerated wrapper script for Libassuan_jll for i686-w64-mingw32
export libassuan

using Libgpg_error_jll
JLLWrappers.@generate_wrapper_header("Libassuan")
JLLWrappers.@declare_library_product(libassuan, "libassuan-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Libgpg_error_jll)
    JLLWrappers.@init_library_product(
        libassuan,
        "bin\\libassuan-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# BlockBandedMatrix with block range indexes is also block-banded
const Block1 = Block{1,Int}
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
const SubBlockSkylineMatrix{T,LL,UU,R1,R2} =
    SubArray{T,2,BlockSkylineMatrix{T,LL,UU},Tuple{BlockSlice{R1},BlockSlice{R2}}}

const SubBandedBlockBandedMatrix{T,R1,R2} =
    SubArray{T,2,<:BandedBlockBandedMatrix{T},Tuple{BlockSlice{R1},BlockSlice{R2}}}

getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, JR::BlockRange1) = BandedBlockBandedMatrix(view(A, KR, JR))
getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, J::Block1) = BandedBlockBandedMatrix(view(A, KR, J))
getindex(A::BandedBlockBandedMatrix, K::Block1, JR::BlockRange1) = BandedBlockBandedMatrix(view(A, K, JR))

@lazymul AbstractBlockBandedMatrix

MemoryLayout(A::Type{<:PseudoBlockArray{T,N,R}}) where {T,N,R} = MemoryLayout(R)

#############
# BLAS overrides
#############

function materialize!(M::MatMulVecAdd{<:AbstractBlockBandedLayout,<:AbstractStridedLayout,<:AbstractStridedLayout})
    α, A, x_in, β, y_in = M.α, M.A, M.B, M.β, M.C
    if length(x_in) != size(A,2) || length(y_in) != size(A,1)
        throw(DimensionMismatch())
    end

    # impose block structure
    y = PseudoBlockArray(y_in, BlockSizes((cumulsizes(blocksizes(A),1),)))
    x = PseudoBlockArray(x_in, BlockSizes((cumulsizes(blocksizes(A),2),)))

    _fill_lmul!(β, y)

    for J = Block.(1:nblocks(A,2))
        for K = blockcolrange(A,J)
            muladd!(α, view(A,K,J), view(x,J), one(α), view(y,K))
        end
    end
    y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout})
    α, A, X, β, Y = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y)
    for J=Block(1):Block(nblocks(X,2)), N=blockcolrange(X,J), K=blockcolrange(A,N)
        muladd!(α, view(A,K,N), view(X,N,J), one(α), view(Y,K,J))
    end
    Y
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockBandedLayout,<:AbstractColumnMajor,<:AbstractColumnMajor})
    α, A, X_in, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    X = PseudoBlockArray(X_in, BlockSizes((cumulsizes(blocksizes(A),2),[1,size(X_in,2)+1])))
    Y = PseudoBlockArray(Y_in, BlockSizes((cumulsizes(blocksizes(A),1), [1,size(Y_in,2)+1])))
    for N=Block.(1:nblocks(X,1)), K=blockcolrange(A,N)
        muladd!(α, view(A,K,N), view(X,N,Block(1)), one(α), view(Y,K,Block(1)))
    end
    Y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:AbstractBlockBandedLayout,<:AbstractColumnMajor})
    α, A_in, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    A = PseudoBlockArray(A_in, BlockSizes(([1,size(A_in,1)+1],cumulsizes(blocksizes(X),1))))
    Y = PseudoBlockArray(Y_in, BlockSizes(([1,size(Y_in,1)+1],cumulsizes(blocksizes(X),2))))
    for J=Block(1):Block(nblocks(X,2)), N=blockcolrange(X,J)
        muladd!(α, view(A,Block(1),N), view(X,N,J), one(α), view(Y,Block(1),J))
    end
    Y_in
end

#############
# * overrides
#############

*(A::BlockBandedMatrix, B::BlockBandedMatrix) = mul(A,B)
*(A::BlockBandedMatrix, B::BandedBlockBandedMatrix) = mul(A,B)
*(A::BandedBlockBandedMatrix, B::BlockBandedMatrix) = mul(A,B)
*(A::BandedBlockBandedMatrix, B::BandedBlockBandedMatrix) = mul(A,B)
*(A::Matrix, B::BlockBandedMatrix) = mul(A,B)
*(A::BlockBandedMatrix, B::Matrix) = mul(A,B)
*(A::BandedBlockBandedMatrix, B::Matrix) = mul(A,B)
*(A::Matrix, B::BandedBlockBandedMatrix) = mul(A,B)

function add_bandwidths(A::AbstractBlockBandedMatrix,B::AbstractBlockBandedMatrix)
    Al,Au = colblockbandwidths(A)
    Bl,Bu = colblockbandwidths(B)

    l = Vector(Bl)
    u = Vector(Bu)

    for (v,Av) in [(l,Al),(u,Au)]
        n = length(v)
        for i = 1:n
            sel = max(i-Bu[i],1):min(i+Bl[i],length(Av))
            isempty(sel) && continue
            v[i] += maximum(Av[sel])
        end
    end

    l,u
end

function add_bandwidths(A::BlockBandedMatrix,B::BlockBandedMatrix)
    l,u = blockbandwidths(A) .+ blockbandwidths(B)
    Fill(l,nblocks(B,2)), Fill(u,nblocks(B,2))
end

function similar(M::MulAdd{<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout}, ::Type{T}) where T
    A,B = M.A, M.B

    Arows, Acols = A.block_sizes.block_sizes.cumul_sizes
    Brows, Bcols = B.block_sizes.block_sizes.cumul_sizes
    if Acols ≠ Brows
        # diagonal matrices can be converted
        if isdiag(B) && size(A,2) == size(B,1) == size(B,2)
            B = BlockBandedMatrix(B.data, BlockSizes((Acols,Acols)), 0, 0, 0, 0)
        elseif isdiag(A) && size(A,2) == size(B,1) == size(A,1)
            A = BlockBandedMatrix(A.data, BlockSizes((Brows,Brows)), 0, 0, 0, 0)
        else
            throw(DimensionMismatch("*"))
        end
    end
    n,m = size(A,1), size(B,2)

    l,u = add_bandwidths(A,B)
    BlockSkylineMatrix{T}(undef, BlockSkylineSizes(BlockSizes((Arows,Bcols)), l, u))
end

function similar(M::MulAdd{BandedBlockBandedColumnMajor,BandedBlockBandedColumnMajor}, ::Type{T}) where T
    A,B = M.A, M.B

    Arows, Acols = A.block_sizes.block_sizes.cumul_sizes
    Brows, Bcols = B.block_sizes.block_sizes.cumul_sizes
    if Acols ≠ Brows
        # diagonal matrices can be converted
        if isdiag(B) && size(A,2) == size(B,1) == size(B,2)
            # TODO: fix
            B = BandedBlockBandedMatrix(B.data, BlockSizes((Acols,Acols)), 0, 0, 0, 0)
        elseif isdiag(A) && size(A,2) == size(B,1) == size(A,1)
            A = BandedBlockBandedMatrix(A.data, BlockSizes((Brows,Brows)), 0, 0, 0, 0)
        else
            throw(DimensionMismatch("*"))
        end
    end
    n,m = size(A,1), size(B,2)

    bs = BandedBlockBandedSizes(BlockSizes((Arows,Bcols)), A.l+B.l, A.u+B.u, A.λ+B.λ, A.μ+B.μ)

    BandedBlockBandedMatrix{T}(undef, bs)
end

similar(M::MulAdd{<:DiagonalLayout,<:AbstractBlockBandedLayout}, ::Type{T}) where T = similar(M.B,T)
similar(M::MulAdd{<:AbstractBlockBandedLayout,<:DiagonalLayout}, ::Type{T}) where T = similar(M.A,T)

function blocksizes(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockRange1}) where {LL,UU}
    A = parent(V)
    Bs = A.block_sizes.block_sizes

    KR = parentindices(V)[1].block.indices[1]
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(KR[1])-Int(JR[1])

    Bs.cumul_sizes[1]
    @assert KR[1] == JR[1] == 1
    BlockSkylineSizes(BlockSizes((Bs.cumul_sizes[1][KR[1]:KR[end]+1] .- Bs.cumul_sizes[1][KR[1]] .+ 1,
                                       Bs.cumul_sizes[2][JR[1]:JR[end]+1] .- Bs.cumul_sizes[2][JR[1]] .+ 1)),
                           colblockbandwidth(A,1)[1:Int(JR[end])] .- shift, colblockbandwidth(A,2)[1:Int(JR[end])] .+ shift)
end

function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockRange1}) where {LL,UU}
    A = parent(V)
    Bs = A.block_sizes.block_sizes

    KR = parentindices(V)[1].block.indices[1]
    JR = parentindices(V)[2].block.indices[1]

    @assert KR[1] == JR[1] == 1
    blockbandwidths(A)
end

####
# BlockIndexRange subblocks
####

sublayout(::BlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{BlockRange1}}}) = BlockBandedColumnMajor()
sublayout(::BlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{Block1}}}) = ColumnMajor()
sublayout(::BlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
sublayout(::BlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockIndexRange1}, BlockSlice{BlockIndexRange1}}}) = ColumnMajor()

sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{Block1}, BlockSlice{Block1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{Block1}, BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockRange1}, BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{Tuple{BlockSlice{BlockIndexRange1}, BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()

isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()

subblockbandwidths(V::SubBandedBlockBandedMatrix) = subblockbandwidths(parent(V))

function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,Block1})
    A = parent(V)
    Bs = A.block_sizes.block_sizes

    KR = parentindices(V)[1].block.indices[1]
    J = parentindices(V)[2].block
    shift = Int(KR[1])-Int(J)
    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,Block1,BlockRange1})
    A = parent(V)
    Bs = A.block_sizes.block_sizes

    K = parentindices(V)[1].block
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(K)-Int(JR[1])

    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,BlockRange1})
    A = parent(V)
    Bs = A.block_sizes.block_sizes

    KR = parentindices(V)[1].block.indices[1]
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(KR[1])-Int(JR[1])

    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {LL,UU} =
    (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])

function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {T,LL,UU}
    A = parent(V)
    JR = parentindices(V)[2]
    KR = parentindices(V)[1].block
    J = parentindices(V)[2].block
    p = unsafe_convert(Ptr{T}, view(A, first(KR), J))
end

struct ShiftedLayout{T,ML<:MemoryLayout} <: MemoryLayout
    shift::Tuple{Int,Int}  # gives the shift to the start of the memory.
                           # So shift == (0,0) is equivalent to layout
                           # shift == (2,1) has the first two rows and first column all zero
    layout::ML
end

strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockIndexRange1}) where {LL,UU} =
    (1,parent(V).block_sizes.block_strides[Int(Block(parentindices(V)[2]))])

function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockRange1,BlockIndexRange1}) where {T,LL,UU}
    A = parent(V)
    JR = parentindices(V)[2]
    K = first(parentindices(V)[1].block)
    J = Block(JR)
    K ∈ blockcolrange(A, J) || throw(ArgumentError("Pointer is only defined when inside colrange"))
    p = unsafe_convert(Ptr{T}, view(A, K, J))
    p + sizeof(T)*(JR.block.indices[1][1]-1)*stride(V,2)
end

function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU}
    A = parent(V)
    JR = parentindices(V)[2]
    K = parentindices(V)[1].block.block
    kr = parentindices(V)[1].block.indices[1]
    J = parentindices(V)[2].block.block
    jr = parentindices(V)[2].block.indices[1]
    p = unsafe_convert(Ptr{T}, view(A, K, J))
    p + sizeof(T)*(kr[1]-1 + (jr[1]-1)*stride(V,2))
end

strides(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU} =
    (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block.block)])

MemoryLayout(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU} = ColumnMajor()


text:
using JuMP
using Gurobi
using MinimaxCalibratedEBayes
using Test
using Plots
pgfplotsx()
gr()
grid = -3:0.06:3
length(grid)
σ = 0.2
gmix_tmp = GaussianMixturePriorClass(σ, grid, Gurobi.Optimizer)

true_dist = MixtureModel(Normal, [(-2,.2), (+2,.2)])

tupl = ( BarConvTol = 1e-12, BarQCPConvTol = 1e-12, FeasibilityTol = 1e-9, OptimalityTol = 1e-9)

gmix_accurate =  GaussianMixturePriorClass(σ, grid,
                                           Gurobi.Optimizer,
                                           tupl)

Z1 = DiscretizedStandardNormalSample(-6:1.0:6)

mhist_mhist = marginalize(true_dist, Z1)
pdf(mhist_mhist)[1:end] =  pdf(mhist_mhist) .+ quantile(pdf(mhist_mhist), 0.1)
pdf(mhist_mhist) .+ 0.01
Z2 = DiscretizedStandardNormalSample(0.0, mhist_mhist)

plot(Z2)

marginal_target = MarginalDensityTarget(StandardNormalSample(0.0))
prior_target = PriorDensityTarget(0.0)

δs = exp10.(range(-3.0, stop=-1.0, length=70))
δs = range(0.001, 0.2, length=50)

δ = 0.001
mm_1 = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_tmp, prior_target, δ)

mm_2 = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_accurate, prior_target, δ)

mms = [MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_accurate, marginal_target, δ) for δ in δs]

mms[1]
var_mms = sqrt.(var.(mms))
max_bias_squared_mms = worst_case_bias.(mms)
mse_mms = sqrt.( var_mms.^2 .+ max_bias_squared_mms.^2)

using LaTeXStrings
using Plots
gr()
plot(sqrt.(δs), [var_mms max_bias_squared_mms mse_mms],
                     seriestype = :path,
                     xlabel = L"\delta",
                     linewidth=4.0,
                     color=["orange" "purple" "black"],
                     label=["variance" "max_bias_squared" "mse"])

pl = steinminimaxplot(mms[3])
pl1!

δ = 0.0001
δp = 0.021
mm = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_tmp, prior_target, δ)

plot(mm.Q)

var(mm)
worst_case_bias(mm)^2

mycols=["#424395" "#018AC4" "#5EC2DA" "#EBC915" "#EB549A" "#550133"]

plot(marginalize(mms[1].g1, mms[1].Z))

pyplot()
pl =steinminimaxplot(mms[30])
plot!(pl[2])
pl

mm_δp = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_tmp, prior_target, δp)

g1 = mm2[:g1]
g2 = mm2[:g2]

plot(-3:0.01:3, x->pdf(mm2[:g1],x))
plot!(-3:0.01:3, x->pdf(mm2[:g2],x))

#f1_cont = marginalize(g1, StandardNormalSample(0.0))

f1 = marginalize(g1, Z2)
f2 = marginalize(g2, Z2)

plot(f1)
plot!(f2, color="red")
plot!(Z2.mhist, color="black")

ω_δ_prime = -JuMP.dual(mm2[:model][:bound_delta])
ω_δ = mm2[:obj_val]

-JuMP.dual(mm_δp[:model][:bound_delta])

mm2[:obj_val] .- mm_δp[:obj_val]

(mm2[:obj_val] .- mm_δp[:obj_val])/(ω_δ_prime*(δ - δp))

L1 = mm2[:L1]
L2 = mm2[:L2]

Q = ω_δ_prime/δ*(pdf(f2) .- pdf(f1))./pdf(Z2)
Q_0  = (L1+L2)/2 - ω_δ_prime/(2*δ)*sum( (pdf(f2) .- pdf(f1)).* (pdf(f2) .+ pdf(f1)) ./ pdf(Z2))
Q_0_prime  = (L1+L2)/2 - ω_δ_prime/2/δ*sum( (pdf(f2).^2 .- pdf(f1).^2) ./ pdf(Z2))

stein_estimator = DiscretizedAffineEstimator(Z2, Q, Q_0)

max_bias = (ω_δ - δ*ω_δ_prime)/2
lp_worst_bias = worst_case_bias(stein_estimator, Z2, gmix_tmp, prior_target)
max_bias/lp_worst_bias[:max_bias]

plot(Q)

plot(tmp, xlim=(-3.5,3.5))
plot(grid(Z2), Q .+ Q_0)

pdf(f2)[1:5]

plot(pdf(f1))
plot!(pdf(f2), color="red")
plot!(pdf(Z2), color="black")

L1 = mm2[:L1]
L2 = mm2[:L2]

marginal_target(Normal(0,1))
L_mid = (L1 + L2)/2

@test mm2[:obj_val] ≈ L2-L1

gr()
plot(f2)
plot!(f1, color="red")

mm2[:obj_val]

mm2_model = mm2[:model]

normalized_rhs(mm2_model[:bound_delta])
set_normalized_rhs(mm2_model[:bound_delta], 0.1)
optimize!(mm2_model)

mm = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_tmp, marginal_target, 0.1)
m1 = mm[:model]
mmc = m1[:bound_delta]
dual(mmc)
normalized_rhs(mmc)
obj_1 = mm[:obj_val]

mm_plus_delta = MinimaxCalibratedEBayes.modulus_problem(Z2, gmix_tmp, marginal_target, 0.1001)
obj_2 = mm_plus_delta[:obj_val]

obj_2 - obj_1
dual(mmc)*(0.1001 - 0.1)

mmc = m1[:pseudo_chisq_constraint]

dual(m1, ""

mm[:g1]
plot(-3:0.01:3, x->pdf(mm[:g1],x))
plot!(-3:0.01:3, x->pdf(mm[:g2],x), color="red")

plot(-5:0.01:5, x->pdf(mm2[:g1],x))
plot!(-5:0.01:5, x->pdf(mm2[:g2],x), color="red")

# In fact all computations with marginalized object.

# param -> A*x

# finite-dim representation, here: probs
# want to map: probs ->


text:
function construct_ψ(U::Function, dUdx::Function, maxevals::Int=10^4)

    integrand(x::Real,z::Complex) = exp(-2*z*x-U(x))*dUdx(x)
    dz_integrand(x::Real,z::Complex) = -exp(-2*z*x-U(x))*dUdx(x)*2*x
    ψ(z::Complex) = quadgk(x -> integrand(x,z),0,Inf,maxevals=maxevals)[1]
    dψdz(z::Complex) = quadgk(x -> dz_integrand(x,z),0,Inf,maxevals=maxevals)[1]

    return(ψ,dψdz)

end

function construct_eigenfunction(U::Function, dUdx::Function, ψ::Function, γ::Complex, maxevals::Int = 10^4)

    f_plus(x) = (x < 0.0) ? ψ(γ) * exp(γ * x) : exp(γ * x + U(x)) * quadgk(ξ -> dUdx(ξ) * exp(-2 * γ * ξ - U(ξ)), x, Inf, maxevals=maxevals)[1];
    f_minus(x) = (x < 0.0) ? -ψ(γ) * exp(-γ * x + U(x)) * quadgk(ξ -> dUdx(ξ) * exp(2 * γ * ξ - U(ξ)),-Inf,x, maxevals=maxevals)[1] : exp(-γ * x);

    return (f_plus, f_minus);

end

function refresh_perturbation(U::Function, dUdx::Function, psi_plus::Function, γ::Complex, maxevals::Int = 10^4)

    (f_plus, f_minus) = construct_eigenfunction(U, dUdx, psi_plus, γ, maxevals);

    integrand_numerator = x -> exp(-U(x))* (f_plus(x)^2 + f_minus(x)^2);
    integrand_denominator = x -> exp(-U(x))* 2 * f_plus(x)*f_minus(x)

    numerator = quadgk(integrand_numerator, -30,30, maxevals=maxevals)[1];
    denominator = quadgk(integrand_denominator, -30, 30, maxevals=maxevals)[1];

    return numerator/denominator - 1;

end

function construct_eigenfunction_symmetric(U::Function, dUdx::Function, ψ::Function, γ::Complex, sgn::Int, maxevals::Int = 10^4)
    # construct +/- eigenfunction in the symmetric case, depending on value of sgn (-1/+1)

    f(x) = (x < 0.0) ? exp(γ * x) : sgn * exp(γ * x + U(x)) * quadgk(ξ -> dUdx(ξ) * exp(-2 * γ * ξ - U(ξ)), x, Inf, maxevals=maxevals)[1];

    return f

end

function refresh_perturbation_symmetric(U::Function, dUdx::Function, ψ::Function, γ::Complex, sgn::Int, maxevals::Int = 10^4)

    f = construct_eigenfunction_symmetric(U, dUdx, ψ, γ, sgn, maxevals);
    integrand_numerator = x -> exp(-U(x)) * f(x)^2;
    integrand_denominator = x -> exp(-U(x)) * f(x) * f(-x);
    numerator = quadgk(integrand_numerator,-10, 10, maxevals=maxevals)[1];
    denominator = 2 * quadgk(integrand_denominator, 0, 10, maxevals=maxevals)[1];
    return sgn * numerator/denominator - 1;

end

max_stars_repo_name: jbierkens/spectral

max_issues_repo_name: jbierkens/spectral

max_forks_repo_name: jbierkens/spectral