text:
using FixedPointNumbers, Statistics, Random, Test
using FixedPointNumbers: bitwidth, rawtype, nbitsfrac
using Base.Checked

SP = VERSION >= v"1.6.0-DEV.771" ? " " : "" # JuliaLang/julia #37085

"""
    target(X::Type, Ss...; ex = :default)

Return a generator which enumerates the target types for testing.

# Arguments
- `X`: target base type
- `Ss`: symbols for specifying the target raw types
  - `:i*` : a `Signed` type if `X === Fixed`, or an `Unsigned` type if `X === Normed`
  - `:s*` : a `Signed` type (not yet supported)
  - `:u*` : an `Unsigned` type (not yet supported)
- `ex`: exhaustivity of `f`s (see also the [`target_f`](@ref) function)
  - `:heavy`: all supported `f`s
  - `:default`: same as `:heavy` for 8-/16-bit types, and same as `:light` otherwise
  - `:light`: important `f`s for byte boundaries and floating point types
  - `:thin`: maximum and half `f`s per type

# Example
```julia
julia> collect(target(Normed, :i8, :i32; ex = :default))
21-element Array{DataType,1}:
 Normed{UInt8,1}
 Normed{UInt8,2}
 Normed{UInt8,3}
 Normed{UInt8,4}
 Normed{UInt8,5}
 Normed{UInt8,6}
 Normed{UInt8,7}
 Normed{UInt8,8}
 Normed{UInt32,1}
 Normed{UInt32,7}
 Normed{UInt32,8}
 Normed{UInt32,9}
 Normed{UInt32,10}
 Normed{UInt32,11}
 Normed{UInt32,15}
 Normed{UInt32,16}
 Normed{UInt32,17}
 Normed{UInt32,23}
 Normed{UInt32,24}
 Normed{UInt32,31}
 Normed{UInt32,32}
```
"""
function target(X::Type, Ss...; ex = :default)
    Ts = symbol_to_inttype.(X, Ss)
    (X{T,f} for T in Ts for f in target_f(X, T; ex = ex))
end
target(X::Type; ex = :default) = target(X, :i8, :i16, :i32, :i64, :i128; ex = ex)

"""
    target_f(X::Type, T::Type; ex = :default)

Return a tuple or range of the number of fractional bits `f` to be tested.

# Arguments
The `X` specifies the target base type, i.e. `Fixed` or `Normed`, and the `T`
specifies the target raw type.

## `ex` keyword
The `ex` specifies the exhaustivity of `f`s.
The following are examples of `target_f(Normed, T)`. The marker `x` means the
target and the marker `-` means not the target.

### `:heavy` -- all supported `f`s
```
            |    3                   2                   1                  |
          f |2 1 0 9 8 7 6 5:4 3 2 1 0 9 8 7:6 5 4 3 2 1 0 9:8 7 6 5 4 3 2 1|
T == UInt8  |               :               :               :x x x x x x x x|
T == UInt16 |               :               :x x x x x x x x:x x x x x x x x|
T == UInt32 |x x x x x x x x:x x x x x x x x:x x x x x x x x:x x x x x x x x|
```
### `:default` -- same as `:heavy` for 8-/16-bit types, and same as `:light` otherwise
```
            |    3                   2                   1                  |
          f |2 1 0 9 8 7 6 5:4 3 2 1 0 9 8 7:6 5 4 3 2 1 0 9:8 7 6 5 4 3 2 1|
T == UInt8  |               :               :               :x x x x x x x x|
T == UInt16 |               :               :x x x x x x x x:x x x x x x x x|
T == UInt32 |x x - - - - - -:x x - - - - - x:x x - - - x x x:x x - - - - - x|
```

### `:light` -- important `f`s for byte boundaries and floating point types
```
            |    3                   2                   1                  |
          f |2 1 0 9 8 7 6 5:4 3 2 1 0 9 8 7:6 5 4 3 2 1 0 9:8 7 6 5 4 3 2 1|
T == UInt8  |               :               :               :x x - - - - - x|
T == UInt16 |               :               :x x - - - x x x:x x - - - - - x|
T == UInt32 |x x - - - - - -:x x - - - - - x:x x - - - x x x:x x - - - - - x|
                             |                         |
                             +--precision(Float32)     +--precision(Float16)
```

### `:thin` -- maximum and half `f`s per type
```
            |    3                   2                   1                  |
          f |2 1 0 9 8 7 6 5:4 3 2 1 0 9 8 7:6 5 4 3 2 1 0 9:8 7 6 5 4 3 2 1|
T == UInt8  |               :               :               :x - - - x - - -|
T == UInt16 |               :               :x - - - - - - -:x - - - - - - -|
T == UInt32 |x - - - - - - -:- - - - - - - -:x - - - - - - -:- - - - - - - -|
```
"""
function target_f(X::Type, T::Type{<:Integer}; ex = :default)
    f_min = X === Fixed ? 0 : 1
    f_max = bitwidth(T) - (T <: Signed)
    ex === :heavy && return f_min:f_max
    ex === :default && bitwidth(T) <= 16 && return f_min:f_max
    ex === :thin && return ((f_max + 1) ÷ 2, f_max)
    if ex === :light || ex === :default
        itr = Iterators.filter(x -> x <= f_max, target_f_series(X, T))
        return (itr...,)
    end
    error()
end

target_f_series(::Type{Fixed}, T::Type{<:Integer}) =
    (0, 1, 7, 8, 9,
     10, 11, 15, 16, 17,
     23, 24, 31, 32, 33,
     52, 53, 63, 64, 65,
     112, 113, 127)

target_f_series(::Type{Normed}, T::Type{<:Integer}) =
    (1, 7, 8, 9,
     10, 11, 15, 16, 17,
     23, 24, 31, 32, 33,
     52, 53, 63, 64, 65,
     112, 113, 127, 128)

# generator for cartesian product
function xypairs(::Type{X}) where X
    xs = typemin(X):eps(X):typemax(X)
    ((x, y) for x in xs, y in xs)
end

function test_floattype(TX::Type)
    @testset "floattype($X)" for X in target(TX, :i8, :i16, :i32, :i64; ex = :heavy)
        @test typemax(X) <= maxintfloat(floattype(X))
    end
end

function test_convert_from_nan(TX::Type)
    @testset "$X(nan)" for X in target(TX; ex = :thin)
        @test_throws ArgumentError X(Inf)
        @test_throws ArgumentError X(-Inf32)
        @test_throws ArgumentError X(NaN)
    end
end

function test_rem_type(TX::Type)
    @testset "% $X" for X in target(TX, :i8, :i16; ex = :thin)
        xs = typemin(X):0.1:typemax(X)
        @test all(x -> x % X === X(x), xs)
    end
end

function test_rem_nan(TX::Type)
    # TODO: avoid undefined behavior
    @testset "nan % $X" for X in target(TX, :i8, :i16, :i32, :i64; ex = :thin)
        @test NaN % X === NaN32 % X === NaN16 % X === zero(X)
    end
end

function test_neg(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xs = typemin(X):eps(X):typemax(X)
        fneg(x) = -float(x)
        @test all(x -> wrapping_neg(wrapping_neg(x)) === x, xs)
        @test all(x -> saturating_neg(x) === clamp(fneg(x), X), xs)
        @test all(x -> !(typemin(X) <= fneg(x) <= typemax(X)) ||
                       wrapping_neg(x) === checked_neg(x) === fneg(x) % X, xs)
    end
end

function test_abs(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xs = typemin(X):eps(X):typemax(X)
        fabs(x) = abs(float(x))
        @test all(x -> wrapping_abs(x) === (x > 0 ? x : wrapping_neg(x)), xs)
        @test all(x -> saturating_abs(x) === clamp(fabs(x), X), xs)
        @test all(x -> !(typemin(X) <= fabs(x) <= typemax(X)) ||
                       wrapping_abs(x) === checked_abs(x) === fabs(x) % X, xs)
    end
end

function test_add(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xys = xypairs(X)
        fadd(x, y) = float(x) + float(y)
        @test all(((x, y),) -> wrapping_sub(wrapping_add(x, y), y) === x, xys)
        @test all(((x, y),) -> saturating_add(x, y) === clamp(fadd(x, y), X), xys)
        @test all(((x, y),) -> !(typemin(X) <= fadd(x, y) <= typemax(X)) ||
                               wrapping_add(x, y) === checked_add(x, y) === fadd(x, y) % X, xys)
    end
end

function test_sub(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xys = xypairs(X)
        fsub(x, y) = float(x) - float(y)
        @test all(((x, y),) -> wrapping_add(wrapping_sub(x, y), y) === x, xys)
        @test all(((x, y),) -> saturating_sub(x, y) === clamp(fsub(x, y), X), xys)
        @test all(((x, y),) -> !(typemin(X) <= fsub(x, y) <= typemax(X)) ||
                               wrapping_sub(x, y) === checked_sub(x, y) === fsub(x, y) % X, xys)
    end
end

function test_mul(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xys = xypairs(X)
        fmul(x, y) = float(x) * float(y) # note that precision(Float32) < 32
        @test all(((x, y),) -> wrapping_mul(x, y) === fmul(x, y) % X, xys)
        @test all(((x, y),) -> saturating_mul(x, y) === clamp(fmul(x, y), X), xys)
        @test all(((x, y),) -> !(typemin(X) <= fmul(x, y) <= typemax(X)) ||
                               wrapping_mul(x, y) === checked_mul(x, y), xys)
    end
end

function test_fdiv(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        xys = xypairs(X)
        fdiv(x, y) = oftype(float(x), big(x) / big(y))
        fdivz(x, y) = y === zero(y) ? float(y) : fdiv(x, y)
        @test all(((x, y),) -> wrapping_fdiv(x, y) === fdivz(x, y) % X, xys)
        @test all(((x, y),) -> saturating_fdiv(x, y) === clamp(fdiv(x, y), X), xys)
        @test all(((x, y),) -> !(typemin(X) <= fdiv(x, y) <= typemax(X)) ||
                               wrapping_fdiv(x, y) === checked_fdiv(x, y), xys)
    end
end

function test_div(TX::Type)
    for X in target(TX, :i8; ex = :thin)
        T = rawtype(X)
        xys = xypairs(X)
        fdiv(x, y) = oftype(float(x), big(x) / big(y))
        @test all(xys) do (x, y)
            rem_t(x) = x > typemax(T) ? typemin(T) : unsafe_trunc(T, x)
            z = y === zero(y) ? float(y) : fdiv(x, y)
            return (wrapping_div(x, y) === rem_t(trunc(z))) &
                   (wrapping_fld(x, y) === rem_t(floor(z))) &
                   (wrapping_cld(x, y) === rem_t( ceil(z)))
        end
        @test all(xys) do (x, y)
            clamp_t(x) = isnan(x) ? zero(T) : trunc(T, clamp(x, typemin(T), typemax(T)))
            z = fdiv(x, y)
            return (saturating_div(x, y) === clamp_t(trunc(z))) &
                   (saturating_fld(x, y) === clamp_t(floor(z))) &
                   (saturating_cld(x, y) === clamp_t( ceil(z)))
        end
        @test all(xys) do (x, y)
            z = fdiv(x, y)
            t = !(typemin(T) <= trunc(z) <= typemax(T)) || wrapping_div(x, y) === checked_div(x, y)
            f = !(typemin(T) <= floor(z) <= typemax(T)) || wrapping_fld(x, y) === checked_fld(x, y)
            c = !(typemin(T) <=  ceil(z) <= typemax(T)) || wrapping_cld(x, y) === checked_cld(x, y)
            return t & f & c
        end
    end
end

function test_div_3arg(TX::Type)
    for X in target(TX; ex = :thin)
        @test div(eps(X), typemax(X), RoundToZero) === div(eps(X), typemax(X))
        @test div(eps(X), typemax(X), RoundDown)   === fld(eps(X), typemax(X))
        @test div(eps(X), typemax(X), RoundUp)     === cld(eps(X), typemax(X))
    end
end

function test_fld1_mod1(TX::Type)
    for X in target(TX, :i8, :i16; ex = :thin)
        T = rawtype(X)
        eps2 = eps(X) + eps(X)
        xs = reinterpret.(X, T.((17, 16, 15, 14)))
        @test all(fld1.(xs, eps2) .=== T.((9, 8, 8, 7)))
        @test_throws DivideError fld1(eps(X), zero(X))

        @test all(mod1.(xs, eps2) .=== reinterpret.(X, T.((1, 2, 1, 2))))
        @test_throws DivideError mod1(eps(X), zero(X))

        d, r = fldmod1(typemin(X), eps2)
        @test d isa T && r isa X && ((d - 1) * eps2 + r) % X === typemin(X)
        d, r = fldmod1(typemax(X), eps2)
        @test d isa T && r isa X && ((d - 1) * eps2 + r) % X === typemax(X)
    end
end

function test_isapprox(TX::Type)
    @testset "approx $X" for X in target(TX, :i8, :i16; ex = :light)
        xs = typemin(X):eps(X):typemax(X)-eps(X)
        @test all(x -> x ≈ x + eps(X), xs)
        @test all(x -> x + eps(X) ≈ x, xs)
        @test !any(x -> x - eps(X) ≈ x + eps(X), xs)
    end

end

function test_clamp_nan(TX::Type)
    @testset "clamp(nan, $X)" for X in target(TX; ex = :thin)
        @test clamp( Inf, X) === clamp( Inf32, X) === typemax(X)
        @test clamp(-Inf, X) === clamp(-Inf32, X) === typemin(X)
        @test clamp( NaN, X) === clamp( NaN32, X) === zero(X)
    end
end

function test_isinteger(TX::Type)
    @testset "isinteger(::$X)" for X in target(TX, :i8, :i16)
        xs = typemin(X):eps(X):typemax(X)
        @test all(x -> isinteger(x) == isinteger(float(x)), xs)
    end
end

function test_rand(TX::Type)
    @testset "rand(::$X)" for X in target(TX; ex = :thin)
        @test isa(rand(X), X)
        a = rand(X, (3, 5))
        @test ndims(a) == 2 && eltype(a) === X
        @test size(a) == (3, 5)
    end
end


text:
# simple methods to process audio signals

import DSP.resample, DSP.arraysplit

"""
    mono(audio::SampleBuf{T, N})

Converts a multichannel audio to single channel

# Output

Return the mono channel audio as a SampleBuf{T, 1} with initial sample rate.

# Arguments

## `audio`

An audio of type SampleBuf{T, N} which is to be converted to single channel audio.

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(1000), 10)
audio_two_channel = SampleBuf(rand(1000, 2), 10)
audio_multi_channel = SampleBuf(rand(1000, 4), 10)

mono(audio_one_channel)
mono(audio_two_channel)
mono(audio_multi_channel)
```

"""
function mono(audio::SampleBuf{T, N}) where {T <: AbstractFloat, N}
    SampleBuf{T, 1}(
        vec(monoch(audio).data),
        audio.samplerate
    )
end

# special code for fixed-point samples, to avoid overflow
function mono(audio::SampleBuf{T, 2}) where {T <: Fixed}
    nchannels = nchannels(audio)
    if nchannels == 1
        SampleBuf{T, 1}(vec(audio.data), audio.samplerate)
    elseif nchannels == 2
        nsamples = nframes(audio)
        buffer = Array{T}(undef, nsamples)
        for i = 1:nsamples
            @inbounds a = audio.data[i, 1].i
            @inbounds b = audio.data[i, 2].i
            m = (a >> 1) + (b >> 1) + (a & b & 1)
            @inbounds buffer[i] = T(m, 0)
        end
        SampleBuf{T, 1}(buffer, audio.samplerate)
    else
        SampleBuf{T, 1}(
            map(T, vec(mean(map(Float32, audio.data)))),
            audio.samplerate
        )
    end
end

"""
    resample(audio::SampleBuf{T, 2}, samplerate::Real)

Resamples input audio with a provided sample rate

# Output

Return the resampled audio as a SampleBuf{T, 2} with given sample rate.

# Arguments

## `audio`

An audio of type SampleBuf{T, 2} which is to be resampled.

## `samplerate`

A parameter to specify the samplerate for final output's sample rate

# Example

```julia
using MusicProcessing

audio_two_channel = SampleBuf(rand(1000, 2), 10)

resample(audio_two_channel, 5)
```
"""
function resample(audio::SampleBuf{T, 2}, samplerate::Real) where T
    rate = samplerate / audio.samplerate
    filter = DSP.resample_filter(rate, 512, 1.0, 140)
    SampleBuf{T, 2}(
        mapslices(audio.data, dims=1) do data
            DSP.resample(data, rate, filter)
        end,
        samplerate
    )
end

"""
    resample(audio::SampleBuf{T, 2}, samplerate::Real)

Resamples input audio with a provided sample rate

"""
function resample(audio::SampleBuf{T, 1}, samplerate::Real) where {T, F}
    SampleBuf{T, 1}(
        DSP.resample(audio.data, samplerate / audio.samplerate),
        samplerate
    )
end

"""
    duration(audio::SampleBuf)

Returns the duration of given audio, in seconds

# Output

Return the duration of audio(in seconds).

# Arguments

## `audio`

An audio of type SampleBuf whose duration is to be found.

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(1000), 10)
audio_two_channel = SampleBuf(rand(1000, 2), 10)
audio_multi_channel = SampleBuf(rand(1000, 4), 10)

duration(audio_one_channel)
duration(audio_two_channel)
duration(audio_multi_channel)
```

"""
function duration(audio::SampleBuf)
    nframes(audio) / samplerate(audio)
end

"""
    play(audio::SampleBuf{Float32})

play the audio on local computer using PortAudio
"""
function play(audio::SampleBuf{Float32})
    # import PortAudio on-demand
    @eval import PortAudio
    nchannels = nchannels(audio)
    stream = PortAudio.PortAudioStream(2, nchannels)
    try
        write(stream, audio)
    finally
        close(stream)
    end
end
play(audio::SampleBuf{T}) where T = play(map(Float32, audio))

"""
    pitchshift(audio::SampleBuf{T, N}, semitones::Real = 8)

Returns pitchshifted audio waveform

# Output

Return the pitchshifted audio waveform as a SampleBuf{T, N}.

# Arguments

## `audio`

An audio of type SampleBuf{T, N} whose pitch is to be shifted.

## `semitones`

A parameter used to define the new sample rate for the audio.\n
`new samplerate` = 2.0 ^ (semitones / 12.0)

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(10000), 10)
audio_two_channel = SampleBuf(rand(10000, 2), 10)
audio_multi_channel = SampleBuf(rand(10000, 4), 10)

pitchshift(audio_one_channel, 8)
pitchshift(audio_two_channel, 8)
pitchshift(audio_multi_channel, 8)

```

"""
function pitchshift(audio::SampleBuf{T, N}, semitones::Real = 0) where {T, N}
    rate = 2.0 ^ (semitones / 12.0)
    shifted = resample(slowdown(audio, rate), audio.samplerate / rate)
    SampleBuf{T, N}(
        shifted.data,
        audio.samplerate
    )
end

"""
    speedup(audio::SampleBuf, speed::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)

Returns a speedup audio waveform using various parameters

# Output

Return the speedup audio waveform as a SampleBuf{T, N} using windowsize and hopsize.

# Arguments

## `audio`

An audio of type SampleBuf{T, N} whose pitch is to be shifted.

## `windowsize`

A parameter used to define the window size for stft and istft functions.

## `hopsize`

A parameter used to define the number audio of frames between stft columns.

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(10000), 10)
audio_two_channel = SampleBuf(rand(10000, 2), 10)
audio_multi_channel = SampleBuf(rand(10000, 4), 10)

speedup(audio_one_channel, 2048, 512)
speedup(audio_one_channel, 2048, 512)
speedup(audio_one_channel, 2048, 512)

```

"""
function speedup(audio::SampleBuf, speed::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)
    S = stft(audio, windowsize, hopsize; kwargs...)
    S = phase_vocoder(S, speed, hopsize)
    istft(S, audio.samplerate, windowsize, hopsize; kwargs...)
end

"""
    slowdown(audio::SampleBuf, ratio::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)

Returns a slowdowned audio using various parameters

# Output

Return the slowdowned audio waveform as a SampleBuf{T, N} using windowsize and hopsize.

# Arguments

## `audio`

An audio of type SampleBuf{T, N} whose speed is to be slowdowned.

## `windowsize`

A parameter used to define the window size for stft and istft functions.

## `hopsize`

A parameter used to define the number audio of frames between stft columns.

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(10000), 10)
audio_two_channel = SampleBuf(rand(10000, 2), 10)
audio_multi_channel = SampleBuf(rand(10000, 4), 10)

slowdown(audio_one_channel, 2048, 512)
slowdown(audio_two_channel, 2048, 512)
slowdown(audio_multi_channel, 2048, 512)

```

"""
function slowdown(audio::SampleBuf, ratio::Real, windowsize::Int = 1024, hopsize::Int = windowsize >> 2; kwargs...)
    speedup(audio, 1.0 / ratio, windowsize, hopsize; kwargs...)
end

"""
    zero_crossing_rate(audio::SampleBuf{T, 1}, framesize::Int = 1024, hopsize::Int = framesize >> 2)

Returns zero crossing rate in a audio waveform

# Output

Returns a array of zero-crossings in y along the selected axis of audio using various parameters

# Arguments

## `audio`

An audio of type SampleBuf{T, N} whose speed is to be slowdowned.

## `framesize`

A parameter used to define the framesize to calculate no of zero crossings.

## `hopsize`

A parameter used to define the number audio of frames between frames size columns.

# Example

```julia
using MusicProcessing

audio_one_channel = SampleBuf(rand(10000), 10)

zero_crossing_rate(audio_one_channel, 2048, 512)
```

"""
function zero_crossing_rate(audio::SampleBuf{T, 1}, framesize::Int = 1024, hopsize::Int = framesize >> 2) where T
    nframes = MusicProcessing.nframes(length(audio.data), framesize, hopsize)
    result = Array{Float32}(undef,nframes)

    offset = 0
    for i = 1:nframes
        result[i] = zero_crossings(audio.data, framesize, offset) / framesize
        offset += hopsize
    end

    result
end

ext: jl


text:
# REAL SPACE
#-------------------------------------------------------------------------------------------------#

"""
    RealSpaceDataGrid{D,T} <: AbstractRealSpaceData{D}

A data grid defined in real space, containing data of type T.
"""
struct RealSpaceDataGrid{D,T} <: AbstractRealSpaceData{D}
    # Basis vectors defining the lattice
    latt::BasisVectors{D}
    # Shift of the origin from the lattice
    orig::SVector{D,Float64}
    # The actual data grid
    grid::Array{T,D}
    # Inner constructor
    function RealSpaceDataGrid(
        latt::BasisVectors{D},
        orig::AbstractVector{<:Real},
        grid::Array{T,D}
    ) where {D,T}
        @assert length(orig) == D "Origin vector does not have the lattice dimensionality"
        # Make sure the shift values lie in (-0.5, 0.5]
        orig = orig - floor.(orig)
        return new{D,T}(latt, orig, grid)
    end
end

"""
    RealSpaceDataGrid(
        latt::AbstractLattice{D},
        orig::AbstractVector{<:Real},
        grid::AbstractArray{T,D};
        prim=true
    )

Creates a real space data grid using lattice information from an `AbstractLattice`. By default,
data is assumed to be given in terms of the primitive lattice (as is usually the case for
computational data).
"""
function RealSpaceDataGrid(
    latt::AbstractLattice{D},
    orig::AbstractVector{<:Real},
    grid::AbstractArray{T,D};
    prim=true
) where {D,T}
    # Conversion for safety
    l = RealLattice(latt)
    if prim
        return RealSpaceDataGrid(prim(l), orig, grid)
    else
        return RealSpaceDataGrid(conv(l), orig, grid)
    end
end

# getindex supports arbitrary integer indices for RealSpaceDataGrid
function Base.getindex(g::RealSpaceDataGrid, inds...)
    # Perform modulo math to get the indices
    # WARNING: Julia % is the remainder function, not modulo!
    imod = mod.(inds .- 1,  gridsize(g)) .+ 1
    return getindex(grid(g), imod...)
end

# Iterator definitions: pass through matrix iteration
Base.iterate(g::RealSpaceDataGrid) = iterate(grid(g))
Base.iterate(g::RealSpaceDataGrid, state) = iterate(grid(g), state)

"""
    basis(g::RealSpaceDataGrid{D,T}) -> BasisVectors{D}

Gets the basis vectors of a `RealSpaceDataGrid`.
"""
basis(g::RealSpaceDataGrid) = g.latt
shift(g::RealSpaceDataGrid) = g.orig

"""
    grid(g::RealSpaceDataGrid{D,T}) -> Array{T,D}

Gets the array that backs a `RealSpaceDataGrid{D,T}`, which is an `Array{T,D}`.
"""
grid(g::RealSpaceDataGrid) = g.grid
# Size of the data grid in entries per dimension
# TODO: should we overload Base.size() as well?

"""
    gridsize(g::RealSpaceDataGrid{D,T}) -> NTuple{D,Int}

Gets the dimensions of the backing array corresponding to a `RealSpaceDataGrid`.
"""
gridsize(g::RealSpaceDataGrid) = size(g.grid)
#Base.size(g::RealSpaceDataGrid) = gridsize(g)

"""
    volume(g::RealSpaceDataGrid) -> Float64

Gets the crystal volume associated with a `RealSpaceDataGrid`.

By default, units are assumed to be cubic angstroms.
"""
volume(g::RealSpaceDataGrid) = volume(basis(g))

"""
    voxelsize(g::RealSpaceDataGrid) -> Float64

Gets the size of a single voxel of a `RealSpaceDataGrid`.

By default, units are assumed to be cubic angstroms.
"""
voxelsize(g::RealSpaceDataGrid) = volume(g) / prod(gridsize(g))

"""
    Xtal.grid_check(g1::RealSpaceDataGrid, g2::RealSpaceDataGrid)

Performs a check on two `RealSpaceDataGrid`s to ensure that the basis, origin shift, and grid
dimensions are the same before performing mathematical operations.
"""
function grid_check(g1::RealSpaceDataGrid, g2::RealSpaceDataGrid)
    @assert basis(g1) === basis(g2) "Grid basis vectors for each grid are not identical."
    @assert shift(g1) === shift(g2) "Grid shifts from origin are not identical."
    @assert size(grid(g1)) === size(grid(g2)) "Grid sizes are different."
    return nothing
end

# TODO: ensure that type inference works
# We might be able to remove some type parameter tests as well
function Base.:+(g1::RealSpaceDataGrid{D,T1}, g2::RealSpaceDataGrid{D,T2}) where {D,T1,T2}
    # Check that the grids are identical
    grid_check(g1, g2)
    # Add the two datagrids elementwise
    newgrid = grid(g1) + grid(g2)
    return RealSpaceDataGrid(basis(g1), shift(g1), newgrid)
end

function Base.:*(g1::RealSpaceDataGrid{D,T1}, g2::RealSpaceDataGrid{D,T2}) where {D,T1,T2}
    # Check that the grids are identical
    grid_check(g1, g2)
    # Add the two datagrids elementwise
    newgrid = grid(g1) .* grid(g2)
    return RealSpaceDataGrid(basis(g1), shift(g1), newgrid)
end

function Base.:*(s::Number, g::RealSpaceDataGrid{D,T}) where {D,T}
    newgrid = s * grid(g)
    return RealSpaceDataGrid(basis(g), shift(g), newgrid)
end

Base.:*(g::RealSpaceDataGrid, s::Number) = s * g

# Multiply everything in the datagrid by -1
Base.:-(g::RealSpaceDataGrid) = RealSpaceDataGrid(basis(g), shift(g), -grid(g))
# Subtract two datagrids
Base.:-(g1::RealSpaceDataGrid, g2::RealSpaceDataGrid) = +(g1, -g2)

"""
    integrate(g::RealSpaceDataGrid{D,T<:Number}) -> <:Number

Performs an integration across all voxels, returning a scalar value.
"""
function integrate(g::RealSpaceDataGrid{D,T}) where {D,T<:Number}
    return sum(grid(g)) * voxelsize(g)
end

"""
    fft(g::RealSpaceDataGrid{D,<:Number}; maxhkl=zeros(Int,D)) -> HKLData{D,<:Complex}

Performs a fast Fourier transform on the data in a `RealSpaceDataGrid{D,<:Number}` and
generates an `HKLData{D,<:Complex}`.
"""
function FFTW.fft(
    g::RealSpaceDataGrid{D,<:Number};
    maxhkl::AbstractVector{<:Integer} = zeros(Int, D)
) where D
    # Generate the bounds needed for the HKLdata
    bounds = [range(div(sz, 2, RoundUp) - sz, length = sz) for sz in gridsize(g)]
    @debug string("Bounds:\t", bounds)
    # Calculate the shift factors to put the values in the right places
    shifts = [last(b) for b in bounds]
    # Take the grid fft
    # Permute the elements to get the indexing working
    f = circshift(fft(grid(g)), shifts .+ 1)
    # With defined bounds, truncate the data
    if all(!iszero, maxhkl)
        ind = [(-x:x) .- first(b) .+ 1 for (x,b) in zip(abs.(maxhkl), bounds)]
        f = f[ind...]
        bounds = [-x:x for x in abs.(maxhkl)]
    end
    return HKLData(f, bounds)
end

#=
"""
    interpolate(g::RealSpaceDataGrid{D,T}, inds...)

Evaluates a datagrid at a fractional index by linearly interpolating the values between nearest
set of points.
"""
function interpolate(g::RealSpaceDataGrid{D,T}, inds::Vararg{<:Real,D}) where {D,T}
    # Get the floor
    f = floor.(inds)
    c = ceil.(inds)
    # If they're the same, it's probably because the indices were integers
    # Just index the array
    if f == c
        return g[inds...]
    end
    # Get the set of points needed to perform the interpolation
end

"""
    interpolate(g::RealSpaceDataGrid{D,T}, r::AbstractVector{<:Real})

Evaluates a datagrid at the reduced coordinate `r` using linear interpolation.
"""
function interpolate(g::RealSpaceDataGrid, r::AbstractVector{<:Real})
    inds = reduced .* size
end
=#

# RECIPROCAL SPACE
#-------------------------------------------------------------------------------------------------#

"""
    KPointGrid{D} <: AbstractKPoints{D}

Contains a grid used to generate k-points during a calculation.

The grid itself is given as an `SMatrix{D,D,Int}`, and can be interpreted as a set of `D` vectors
given in terms of the primitive basis. These vectors can alternatively be used to construct a
supercell.

The shift of the k-point mesh off Γ is given as an `SVector{D,Float64}`.
"""
struct KPointGrid{D} <: AbstractKPoints{D}
    grid::SMatrix{D,D,Int}
    orig::SVector{D,Float64}
    function KPointGrid{D}(grid::AbstractMatrix{<:Integer}, orig::AbstractVector{<:Real}) where D
        # only allow positive values in the grid matrix
        @assert all(x -> x > 0, grid) "negative values are disallowed in the grid matrix"
        # Keep shift inside the Brillouin zone
        orig = orig -  round.(orig)
        return new(grid, orig)
    end
end

"""
    KPointList{D} <: AbstractKPoints{D}

Contains a list of k-points. This is useful for describing an ordered list of k-points that are
not associated with a mesh - for instance, those used in band structure calculations. This can
also be used to store lists of k-points that are generated by a `KPointGrid{D}`.
"""
struct KPointList{D} <: AbstractKPoints{D}
    list::Vector{SVector{D,Float64}}
end

# Index like the vector it is internally
Base.getindex(k::KPointList, i) = k.list[i]

function Base.setindex!(k::KPointList, v::AbstractVector, i)
    k.list[i] = v
end

"""
    nkpt(k::KPointList{D}) -> Int

Gets the number of k-points in a `KPointList`.
"""
nkpt(k::KPointList{D}) where D = length(k.list)

#= TODO: figure out how to get a k-point list from a grid
function KPointList{D}(k::KPointGrid)

end
=#

"""
    BandAtKPoint

Stores information about a band's energy and its occupancy at a specific k-point.
"""
struct BandAtKPoint
    # Energies
    e::Vector{Float64}
    # Occupancy
    occ::Vector{Float64}
    function BandAtKPoint(e::AbstractVector{<:Real}, occ::AbstractVector{<:Real})
        @assert length(e) == length(occ) "Size of energy and occupancy arrays do not match."
        return new(e, occ)
    end
end

"""
    BandAtKPoint(eocc::AbstractVector{NTuple{2,<:Real}})

Constructs a new `BandAtKPoint` from a vector containing tuples of energy and occupancy data
(in that order).
"""
function BandAtKPoint(eocc::AbstractVector{<:NTuple{2,<:Real}})
    return BandAtKPoint([x[1] for x in eocc], [x[2] for x in eocc])
end

# Access any pair of energy and occupancy with indexing
function Base.getindex(b::BandAtKPoint, inds...)
    return (b.e[inds...], b.occ[inds...])
end

"""
    nband(b::BandAtKPoint) -> Int

Returns the number of bands associated with a k-point.
"""
nband(b::BandAtKPoint) = length(b.e)

"""
    BandStructure{D}

Stores information about an electronic band structure, including the list of k-points used to
generate the data (as a `KPointList{D}`)and the band information at every k-point (as a
`Vector{BandAtKPoint}`).
"""
struct BandStructure{D} <: AbstractReciprocalSpaceData{D}
    # k-points for which band data is defined
    kpts::KPointList{D}
    # Set of energy and occupancy data
    bands::Vector{BandAtKPoint}
    function BandStructure{D}(kpts::KPointList{D}, bands::Vector{BandAtKPoint}) where D
        @assert nkpt(kpts) == length(bands) "Incorrect number of k-points or band datasets."
        @assert _allsame(length(bands)) "Number of bands is inconsistent."
        return new(kpts, bands)
    end
end

"""
    BandStructure{D}(kpts::AbstractKPoints{D}, bands::AbstractVector{<:BandAtKPoint}) where D

Generates a new band structure from k-point information and a vector containing band information
at each k-point.
"""
function BandStructure{D}(kpts::AbstractKPoints{D}, bands::AbstractVector{<:BandAtKPoint}) where D
    return BandStructure{D}(kpts, bands)
end

# Get the pair of a k-point and associated band data
function Base.getindex(b::BandStructure{D}, inds...) where D
    return (b.kpts[inds...], b.bands[inds...])
end

nkpt(b::BandStructure{D}) where D = nkpt(b.kpts)
nband(b::BandStructure{D}) where D = nband(b.bands[1])

"""
    FatBands{D} <: AbstractReciprocalSpaceData{D}

Stores information relevant to plotting fatbands.

- FatBands.bands: matrix of energies at each [kpt, band].
- FatBands.projband: array of lm-decomposed band structure. [orbital, ion, band, kpt].
- FatBands.cband: array of complex-valued contributions to band structure.
"""
struct FatBands{D} <: AbstractReciprocalSpaceData{D}
    bands::Matrix{Float64}
    projband::Array{Float64,4}
    cband::Array{Complex{Float64},4}
end

"""
    HKLData{D,T} <: AbstractReciprocalSpaceData{D}

Stores information associated with specific sets of reciprocal lattice vectors. Data can be
accessed and modified using regular indexing, where indices may be negative.
"""
struct HKLData{D,T} <: AbstractReciprocalSpaceData{D}
    # the actual data
    data::Array{T,D}
    # the bounds in each dimension
    # mutable since the dimensions of Array{D,T} can be changed, in principle
    bounds::MVector{D,UnitRange{Int}}
    function HKLData(
        data::AbstractArray{T,D},
        bounds::AbstractVector{<:AbstractRange{<:Integer}}
    ) where {D,T}
        # The size of the array should match the bounds given
        # For instance, a HKLData with bounds  [-10:10, -10:10, -10:10] should Be
        # [21, 21, 21]
        @assert [s for s in size(data)] == [length(r) for r in bounds] string(
            "Array size incompatible with bounds."
        )
        return new{D,T}(data, bounds)
    end
end

# Needed because HKLData will nearly always have unexpected indices
Base.has_offset_axes(g::HKLData) = true

"""
    Xtal.shiftbounds(g::HKLData{D,T}, inds) -> NTuple{D,<:Integer}

Checks that integer array indices used to access data in an `HKLData` are valid and shifts them to
access the correct portions of the backing array.
"""
function shiftbounds(hkl::HKLData{D,T}, inds) where {D,T}
    # Check that all the indices are in bounds
    if !mapreduce((i,r) -> i in r, &, inds, bounds(hkl))
        # TODO: does this produce a reasonable error message with correct bounds?
        throw(BoundsError(hkl, inds))
    end
    # Adjust the indices to match the array
    # Subtract the minimum index then add 1
    # So if the range is -10:10, an index of 0 should be 0 - -10 + 1 = 11
    i = inds .- minimum.(bounds(hkl)) .+ 1
    return i
end

bounds(hkl::HKLData) = hkl.bounds

# HKLData now supports indexing by Miller index
function Base.getindex(g::HKLData{D,T}, inds...) where {D,T}
    i = shiftbounds(g, inds)
    return g.data[i...]
end

function Base.setindex!(g::HKLData{D,T}, x::T, inds...) where {D,T}
    i = shiftbounds(g, inds)
    g.data[i...] = x
end

function HKLData(
    a::AbstractArray{T,D},
    bounds::Vararg{AbstractUnitRange{<:Integer},D}
) where {D,T}
    return HKLData(a, bounds)
end

function Base.zeros(
    ::Type{HKLData{D,T}},
    bounds::Vararg{AbstractUnitRange{<:Integer}, D}
) where {D,T}
    data = zeros(T, length.(bounds))
    return HKLData(data, MVector{D,UnitRange{Int}}(bounds))
end

"""
    ReciprocalWavefunction{D,T<:Real} <: AbstractReciprocalSpaceData{D}

Contains a wavefunction stored by k-points and bands in a planewave basis. Used to store data in
VASP WAVECAR files. Each k-point is expected to have the same number of bands.

Every band has associated data containing coefficients of the constituent planewaves stored in a
`HKLData{D,Complex{T}}`. Unlike most data structures provided by this package, the type of
complex number used does not default to `Float64`: wavefunction data is often supplied as a
`Complex{Float32}` to reduce the size of the data.
"""
struct ReciprocalWavefunction{D,T<:Real} <: AbstractReciprocalSpaceData{D}
    # Reciprocal lattice on which the k-points are defined
    rlatt::BasisVectors{D}
    # Planewave coefficients: a Matrix (size nkpt*maxnband) of HKLData
    waves::Matrix{HKLData{D,Complex{T}}}
    function ReciprocalWavefunction{D,T}(
        rlatt::BasisVectors{D},
        waves::AbstractMatrix{HKLData{D,Complex{T}}}
    ) where {D,T<:Real}
        # Checks were removed here
        return new(rlatt, waves)
    end
end

function ReciprocalWavefunction{D,T}(
    latt::AbstractLattice{D},
    waves::AbstractMatrix{HKLData{D,Complex{T}}};
    prim = true
) where {D,T<:Real}
    M = prim ? prim(latt) : conv(latt)
    return ReciprocalWavefunction{D,T}(M, waves)
end

# Getting indices should pull from the waves struct: wf[kpt, band]
Base.getindex(wf::ReciprocalWavefunction{D,T}, inds...) where {D,T} = wf.waves[inds...]

nkpt(wf::ReciprocalWavefunction{D,T}) where {D,T} = size(wf.waves, 1)
nband(wf::ReciprocalWavefunction{D,T}) where {D,T} = size(wf.waves, 2)

# DENSITY OF STATES
#-------------------------------------------------------------------------------------------------#

"""
    DensityOfStates

Contains the total density of states information.
"""
struct DensityOfStates <: AbstractDensityOfStates
    # Fermi energy
    fermi::Float64
    # Energy at each point
    energy::Vector{Float64}
    # Number of states at each energy
    dos::Vector{Float64}
    # Integrated density of states (electron count up to that energy)
    int::Vector{Float64}
    function DensityOfStates(
        fermi::Real,
        energy::AbstractVector{<:Real},
        dos::AbstractVector{<:Real},
        int::AbstractVector{<:Real}
    )
        @assert size(energy) == size(dos) == size(int) string(
            "Number of energy and DOS entries do not match."
        )
        return new(fermi, energy, dos, int)
    end
end

# Create a DOS curve without integrated DOS data
"""
    DensityOfStates(fermi::Real, energy::AbstractVector{<:Real}, dos::AbstractVector{<:Real})

Creates a `DensityOfStates` object without integrated DOS information by integrating the state
density.
"""
function DensityOfStates(
    fermi::Real,
    energy::AbstractVector{<:Real},
    dos::AbstractVector{<:Real}
)
    # Generate an integrated DOS vector
    # TODO: perhaps we can improve this with a trapezoidal approximation?
    int = [sum(dos[1:n]) for n in 1:length(dos)]
    return DensityOfStates(fermi, energy, dos, int)
end

# Getting an index produces a tuple of the energy,
Base.getindex(d::DensityOfStates, ind) = (d.energy[ind], d.dos[ind], d.int[ind])

# Setting up for gaussian smearing of DOS curve
function gaussian(
    dos::DensityOfStates,
    sigma::Real,
)
    # Center gaussian within energy range
    # Define gaussian
    # Permute gaussian to prevent shifting of smear output
    ctr = dos.energy[div(length(dos.energy), 2)]
    g = exp.(-((dos.energy .- ctr)/sigma).^2)
    return circshift(g, 1 - findmax(g)[2])
end

"""
    smear(dos::DensityOfStates, sigma::Real)

Smears the DOS function by convoluting it with a gaussian, sigma will determine how
much the DOS is smeared.
"""
function smear(
    dos::DensityOfStates,
    sigma::Real
)
    # Convolute DOS with gaussian using Fourier transforms
    gsmear = real(ifft(fft(dos.dos) .* fft(gaussian(dos, sigma))))
    # Normalize gaussian smearing so that sum(gaussian smear) = sum(DOS)
    return DensityOfStates(fermi(dos), energies(dos), (gsmear)*(sum(dos.dos)/sum(gsmear)))
end

"""
    ProjectedDensityOfStates

Contains projected density of states information.
"""
struct ProjectedDensityOfStates
    # Fermi energy
    fermi::Float64
    # Range of energies
    energy::Vector{Float64}
    # Matrix containing projected DOS data
    # Columns are the components, rows are the energies
    dos::Matrix{Float64}
    # Integrated DOS in the same format
#    int::Matrix{Float64}
    function ProjectedDensityOfStates(
        fermi::Real,
        energy::AbstractVector{<:Real},
        dos::AbstractMatrix{<:Real}
        #int::AbstractMatrix{<:Real}
    )
        @assert length(energy) == size(dos,2) "Number of energy and DOS entries do not match."
        return new(fermi, energy, dos) #, int)
    end
end

#="""
    ProjectedDensityOfStates(
        fermi::Real,
        energy::AbstractVector{<:Real},
        dos::AbstractVector{<:Real}
    )

Creates a `ProjectedDensityOfStates` object without integrated DOS information by integrating
the state density.
"""
function ProjectedDensityOfStates(
    fermi::Real,
    energy::AbstractVector{<:Real},
    dos::AbstractMatrix{<:Real}
)
    # Generate an integrated DOS matrix
    # TODO: perhaps we can improve this with a trapezoidal approximation?
    int = hcat(vec(sum(M[:,1:n], dims=2)) for n in 1:size(M,2))
    return ProjectedDensityOfStates(fermi, energy, dos, int)
end=#

"""
    fermi(d::AbstractDensityOfStates) -> Float64

Gets the Fermi energy from DOS data. There are no guarantees on the unit of energy used!
"""
fermi(d::AbstractDensityOfStates) = d.fermi

"""
    energies(d::AbstractDensityOfStates; usefermi=false) -> Vector{Float64}

Gets the range of energies in the dataset. If `usefermi` is set to true, the energies returned will
be adjusted such that the Fermi energy is set to zero.

There are no guarantees on the unit of energy used!
"""
energies(d::AbstractDensityOfStates; usefermi=false) = d.energy .- (usefermi * d.fermi)

"""
    nelectrons(d::DensityOfStates)

Gets the approximate number of electrons that are needed to reach the Fermi level.
"""
function nelectrons(d::DensityOfStates)
    # Get the nearest entries to the Fermi energy
    E = energies(d, usefermi=true)
    # Find the closest pair of energies to the Fermi energy
    inds = partialsortperm(abs.(E), 1:2)
    # Interpolate linearly between the nearest points
    # TODO: there's probably a nicer way to do this.
    (e1, e2) = [d[i][1] for i in inds]
    (i1, i2) = [d[i][3] for i in inds]
    m = (i2 - i1)/(e2 - e1)
    b = i1 - m*e1
    return m*fermi(d) + b
end

# ATOMIC DATA
#-------------------------------------------------------------------------------------------------#

"""
    AtomicData{D,T}

Data associated with individual atoms in a structure.

This is a type alias for `Dict{AtomPosition{D},T}`. Keys are `AtomPosition` entries, and the values
may be of any type.
"""
const AtomicData{D,T} = Dict{AtomPosition{D},T} where {D,T}

# Note: @computed structs cannot be documented normally
# Use an @doc after the struct, like such
@computed struct SphericalComponents{Lmax}
    v::NTuple{(Lmax+1)^2,Float64}
    # Default constructor without parameters takes numbers directly
    function SphericalComponents(x::Vararg{<:Real,N}) where N
        L = sqrt(length(x)) - 1
        if isinteger(L)
            Lmax = Int(L)
        else
            throw(ArgumentError(string("Cannot determine Lmax from number of arguments.")))
        end
        return new{Lmax}(x)
    end
    # Default constructor with parameters takes any iterator
    function SphericalComponents{Lmax}(v) where Lmax
        @assert length(v) == (Lmax+1)^2 "For Lmax == $Lmax, iterator have length $((Lmax+1)^2)"
        return new{Lmax}(Tuple(v))
    end
end

@doc """
    SphericalComponents{Lmax}

Real spherical harmonic components up to `Lmax`. This can be used to describe atomic orbitals or
projections of data onto atomic sites.
""" SphericalComponents

SphericalComponents(v::SVector{N,<:Real}) where N = SphericalComponents(v...)
SphericalComponents(t::NTuple{N,<:Real}) where N = SphericalComponents(t...)

"""
    Xtal.sc_ind(l::Integer, m::Integer) -> Int

Gets the associated linear index for a pair of (l,m) values used in `SphericalComponents`.
"""
sc_ind(l::Integer, m::Integer) = l^2 + l + 1 + m

# TODO: finish the inverse of the above function
# sc_ind(x) =

function Base.getindex(s::SphericalComponents{Lmax}, l::Integer, m::Integer) where Lmax
    abs(m) <= l || error("|m| must be less than l")
    l <= Lmax || error("l exceeds lmax ($Lmax)")
    return s.v[sc_ind(l, m)]
end


text:
mutable struct ExplorativeMarginQss <: HybridQss
    occ::SVDD.SVDDClassifier
    occ_eps::SVDD.SVDDnegEps
    occ_eps_init_strat::SVDD.InitializationStrategyCombined
    solver::JuMP.MOI.OptimizerWithAttributes
    optimizer::QuerySynthesisOptimizer
    boundary_shift_agg_func::Symbol
    lambda::Float64
    use_penalty::Bool
    eps
    function ExplorativeMarginQss(occ, data; solver::JuMP.MOI.OptimizerWithAttributes, optimizer::QuerySynthesisOptimizer,
                                  boundary_shift_agg_func=:maximum, lambda=1.0, use_penalty=true)
        !isa(occ.kernel_fct, SquaredExponentialKernel) && throw(ArgumentError("Invalid kernel type $(typeof(occ.kernel_fct)). Expected type is a SquaredExponentialKernel."))
        occ_eps = SVDD.SVDDnegEps(data, fill(:U, size(data, 2)))
        occ_params = SVDD.get_model_params(occ)
        !haskey(occ_params, :C) && !haskey(occ_params, :C1) && throw(ArgumentError("Invalid base learner type $(typeof(occ)). Cannot retrieve C parameter for SVDDnegEps."))
        occ_eps_init_strat = SVDD.SimpleCombinedStrategy(SVDD.FixedGammaStrategy(MLKernels.GaussianKernel(occ.kernel_fct.alpha.value.x)),
                                                         SVDD.FixedCStrategy(haskey(occ_params, :C) ? occ_params[:C] : occ_params[:C1]))
        new(occ, occ_eps, occ_eps_init_strat, solver, optimizer, boundary_shift_agg_func, lambda, use_penalty, nothing)
    end
end

function qs_score_function(qs::ExplorativeMarginQss, data::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Function where T <: Real
    # Init epsilon if not yet set
    if qs.eps === nothing
        data_target, data_outliers = SVDD.generate_binary_data_for_tuning(data)
        qs.eps = max(0.0, qs.lambda * estimate_boundary_shift_epsilon(qs.occ, data_outliers, agg_func=qs.boundary_shift_agg_func))
        SVDD.set_eps!(qs.occ_eps, qs.eps)
    end
    if (qs.occ_eps.state === SVDD.model_created)
        SVDD.initialize!(qs.occ_eps, qs.occ_eps_init_strat)
    end
    SVDD.set_data!(qs.occ_eps, data)
    SVDD.set_pools!(qs.occ_eps, labels)
    # Workaround: redirect solver output
    stdout_orig, stderr_orig = stdout, stderr
    redirect_stdout(); redirect_stderr()
    status = SVDD.fit!(qs.occ_eps, qs.solver)
    redirect_stdout(stdout_orig); redirect_stderr(stderr_orig)
    if status === JuMP.MathOptInterface.ALMOST_LOCALLY_SOLVED ||
        status === JuMP.MathOptInterface.LOCALLY_SOLVED
        warn(getlogger(@__MODULE__), "Qss occ_eps not solved to optimality. Solver status: $status.")
    elseif status !== JuMP.MathOptInterface.OPTIMAL
        error(getlogger(@__MODULE__), "Qss occ_eps not solved to optimality. Solver status: $status.")
    end
    oc_scoring(x) = -abs.(SVDD.predict(qs.occ_eps, x))

    # Train penalty binary SVM if outlier labels are available
    if qs.use_penalty && haskey(labels, :Lout)
        labels_binary = fill(:inlier, size(data, 2))
        labels_binary[labels[:Lout]] .= :outlier
        model_binary = LIBSVM.svmtrain(data, labels_binary, gamma=qs.occ_eps.kernel_fct.alpha.value.x, cost=100_000.0)
        penalty(x) = LIBSVM.svmpredict(model_binary, x)[2][1, :]
        return x -> oc_scoring(x) .- max.(zeros(size(x, 2)), -penalty(x))
    end
    return oc_scoring
end


text:
using JLD2
import DataStructures: OrderedDict
import Test: @test, @testset
import DSGE: @test_matrix_approx_eq, update!
using DSGE

### Model
m = KrusellSmithCT()

### Model indices

# Endogenous states
endo = m.endogenous_states
@testset "Model indices" begin
    @test n_states(m) == length(endo[:value_function]) .+ length(endo[:distribution]) .+ length(endo[:log_aggregate_tfp]) + 6
    @test get_setting(m, :n_jump_vars) == length(endo[:value_function])
    @test get_setting(m, :n_state_vars) == length(endo[:distribution]) .+ length(endo[:log_aggregate_tfp]) - get_setting(m, :n_state_vars_unreduce)

    # Exogenous shocks
    @test n_shocks_exogenous(m) == 1

    # Expectation shocks
    @test n_shocks_expectational(m) == length(endo[:value_function])
end

### Steady state
#=
# Load test values
# TODO: These file paths are out of date
params = load("../../../test_outputs/krusell_smith/saved_params.jld2", "params")
grids = load("../../../test_outputs/krusell_smith/saved_params.jld2", "grids")
init_params = load("../../../test_outputs/krusell_smith/saved_params.jld2", "init_params")
approx_params = load("../../../test_outputs/krusell_smith/saved_params.jld2", "approx_params")
vars_ss_mat = load("../../../test_outputs/krusell_smith/vars_ss.jld2")
vars_ss_mat = vars_ss_mat["varsSS"]

# Update parameters
m[:γ] = params[:γ]
m[:ρ] = params[:ρ]
m[:δ] = params[:δ]
m[:α] = params[:α]
m[:σ_tfp] = params[:σ_tfp]
m[:ρ_tfp] = params[:ρ_tfp]
m[:λ1] = params[:λ1]
m[:λ2] = params[:λ2]
m[:μ] = params[:μ]
m[:τ] = params[:τ]
=#
# Update settings
update!(m.settings[:I], Setting(:I, 100))
update!(m.settings[:amin], Setting(:amin, 0.))
update!(m.settings[:amax], Setting(:amax, 100.))
update!(m.settings[:a], Setting(:a, collect(range(0., stop=100., length=100))))
update!(m.settings[:da], Setting(:da, 100 ./ (100 - 1)))
update!(m.settings[:aaa], Setting(:aaa, reshape(get_setting(m, :aa), 2 * 100, 1)))
update!(m.settings[:J], Setting(:J, 2))
update!(m.settings[:z], Setting(:z, [0.; 1.]))
update!(m.settings[:dz], Setting(:dz, 1.))
update!(m.settings[:zz], Setting(:zz, ones(100, 1) * get_setting(m, :z)'))
update!(m.settings[:zzz], Setting(:zzz, reshape(get_setting(m, :zz), 2 * 100, 1)))
endo = m.endogenous_states
update!(m.settings[:n_jump_vars], Setting(:n_jump_vars, length(endo[:value_function])))
update!(m.settings[:n_state_vars], Setting(:n_state_vars, length(endo[:distribution]) .+ n_shocks_exogenous(m)))
update!(m.settings[:n_state_vars_unreduce], Setting(:n_state_vars_unreduce, 0))
update!(m.settings[:rmax], Setting(:rmin, .0001))
update!(m.settings[:rmax], Setting(:rmax, m[:ρ].value))
update!(m.settings[:r0], Setting(:r0, .005))
update!(m.settings[:maxit_HJB], Setting(:maxit_HJB, 100))
update!(m.settings[:crit_HJB], Setting(:crit_HJB, 1e-6))
update!(m.settings[:Δ_HJB], Setting(:Δ_HJB, 1e4))
update!(m.settings[:Ir], Setting(:Ir, 100))
update!(m.settings[:crit_S], Setting(:crit_S, 1e-5))
update!(m.settings[:reduce_state_vars], Setting(:reduce_state_vars, true))
update!(m.settings[:reduce_v], Setting(:reduce_v, true))
update!(m.settings[:krylov_dim], Setting(:krylov_dim, 5))
update!(m.settings[:n_knots], Setting(:n_knots, 12))
update!(m.settings[:c_power], Setting(:c_power, 7))
update!(m.settings[:knots_dict], Setting(:knots_dict, Dict(1 => collect(range(get_setting(m, :amin),
                                                                              stop=get_setting(m, :amax),
                                                                              length=(12 - 1))))))
update!(m.settings[:spline_grid], Setting(:spline_grid, get_setting(m, :a)))
update!(m.settings[:n_prior], Setting(:n_prior, 1))
update!(m.settings[:n_post], Setting(:n_post, 2))
#=
# compare steady state values
steadystate!(m)
@testset "Steady State" begin
    @test @test_matrix_approx_eq vec(m[:V_ss].value) vec(vars_ss_mat[:VSS])
    @test @test_matrix_approx_eq vec(m[:gg_ss].value) vec(vars_ss_mat[:ggSS])
    @test m[:K_ss].value ≈ vars_ss_mat[:KSS] # allow for machine error
    @test m[:r_ss].value ≈ vars_ss_mat[:rSS]
    @test m[:w_ss].value ≈ vars_ss_mat[:wSS]
    @test m[:Y_ss].value ≈ vars_ss_mat[:YSS]
    @test m[:C_ss].value ≈ vars_ss_mat[:CSS]
    @test m[:I_ss].value ≈ vars_ss_mat[:ISS]
end
=#
### Equilibrium conditions
Γ0, Γ1, Ψ, Π, C = eqcond(m)

# Transition and measurement equations
TTT, RRR, CCC, ~ = solve(m)
# measurement equation here

# Matrices have expected dimensions
@testset "Check that outputted matrices from solve have the right dimensions" begin
    n_vars = n_states(m)
    @test size(Γ0) == (n_vars, n_vars)
    @test size(Γ1) == (n_vars, n_vars)
    @test size(vec(C)) == (n_vars,)
    @test size(Ψ) == (n_vars, 1)
    @test size(Π) == (n_vars, n_shocks_expectational(m))
    if get_setting(m, :reduce_v)
        @test size(TTT) == (get_setting(m, :n_splined) + get_setting(m, :n_state_vars_red), get_setting(m, :n_splined) + get_setting(m, :n_state_vars_red))
        @test size(RRR) == (get_setting(m, :n_splined) + get_setting(m, :n_state_vars_red), 1)
        @test size(vec(CCC)) == (get_setting(m, :n_splined) + get_setting(m, :n_state_vars_red),)
    elseif get_setting(m, :reduce_state_vars)
        @test size(TTT) == (get_setting(m, :n_jump_vars) + get_setting(m, :n_state_vars_red), get_setting(m, :n_state_vars_red))
        @test size(RRR) == (get_setting(m, :n_jump_vars) + get_setting(m, :n_state_vars_red), 1)
        @test size(vec(CCC)) == (get_setting(m, :n_jump_vars) + get_setting(m, :n_state_vars_red),)
    else
        @test size(TTT) == (n_vars, n_vars)
        @test size(RRR) == (n_vars, 1)
        @test size(vec(CCC)) == (n_vars,)
    end
end


text:
export optcontract, treecontract, optimaltree

"""
find the optimal contraction tree.
"""
function TensorOperations.optimaltree(network, size_dict::Dict{Int, Int}=Dict{Int, Int}())
    unique_tokens = union(network...)
    optimaltree(network, Dict(token=>Power{:χ}(get(size_dict, token, 1),1) for token in unique_tokens))
end

function leg_analysis(IVS...)
    IALL = union(IVS...)
    II = intersect(IVS...)
    IC = setdiff(IALL, II)
    IALL, II, IC
end

_treecontract(tree::Int, ixs, xs, iy::Nothing) = xs[tree], ixs[tree]
function _treecontract(tree::Int, ixs, xs, iy)
    iy0, y = ixs[tree], xs[tree]
    einsum(EinCode{(iy0,), iy}(), (C,)), iy
end

function _treecontract(tree, ixs, xs, iy)
    i, j = tree
    A, IA = _treecontract(i, ixs, xs, nothing)
    B, IB = _treecontract(j, ixs, xs, nothing)
    _iy = iy == nothing ? Tuple(leg_analysis(IA, IB)[3]) : iy
    ixs = (IA, IB)
    code = EinCode{ixs, _iy}()
    @show size(A), size(B)
    einsum(code, (A, B)), _iy
end

function treecontract(tree, ixs, xs, iy)
    _treecontract(tree, ixs, xs, iy) |> first
end

function get_size_dict(ixs, xs)
    nt = length(ixs)
    size_dict = Dict{Int, Int}()
    @inbounds for i = 1:nt
        for (N, leg) in zip(size(xs[i]), ixs[i])
            if haskey(size_dict, leg)
                size_dict[leg] == N || throw(DimensionMismatch("size of contraction leg $leg not match."))
            else
                size_dict[leg] = N
            end
        end
    end
    return size_dict
end

function optcontract(ixs, xs, iy)
    size_dict = get_size_dict(ixs, xs)

    print("finding optimal tree ...")
    tree, cost = optimaltree(ixs, size_dict)
    @show tree, cost
    treecontract(tree, ixs, xs, iy)
end


text:
"""
    Shapes

A module encapsulating types and methods for bidimensional shapes.

Exported types: `PolygonXor`

Useful functions:
"""
module Shapes

using LinearAlgebra
using StaticArrays
using FixedPointNumbers
using FastClosures
using Triangle
using DataStructures # SortedSet
import Clipper

norm²(x) = x[1]*x[1]+x[2]*x[2]
distance²(x,y) = norm²(x-y)
Path{D,T} = Vector{SVector{D,T}}
# @inline det(u::StaticVector{2}, v::StaticVector{2}) =
# 	@inbounds u[1]*v[2]-u[2]*v[1]
# @inline det(p::StaticVector{2}, q::StaticVector{2}, r::StaticVector{2}) =
# 	det(q-p, r-p)

# Clipper types and constants««1
# default number of bits for Clipper types
# this *must* be a 64-bit type:
const _CLIPPER_FIXED = Fixed{Int64,16}
const _CLIPPER_ONE = FixedPointNumbers.rawone(_CLIPPER_FIXED)
# constants
const _CLIPPER_ENUM = (#««
	clip=(
		union       =Clipper.ClipTypeUnion,
		intersection=Clipper.ClipTypeIntersection,
		difference  =Clipper.ClipTypeDifference,
		xor         =Clipper.ClipTypeXor,
	),
	# the enum order in `Clipper.jl` is wrong, we fix it here:
	ends=(
		fill  = Clipper.EndType(0), #Clipper.EndTypeClosedPolygon,
		loop  = Clipper.EndType(1), #Clipper.EndTypeClosedLine,
		butt  = Clipper.EndType(2), #Clipper.EndTypeOpenButt,
		square= Clipper.EndType(3), #Clipper.EndTypeOpenSquare,
		round = Clipper.EndType(4), #Clipper.EndTypeOpenRound,
	),
	join=(
		square=Clipper.JoinTypeSquare,
		round =Clipper.JoinTypeRound,
		miter =Clipper.JoinTypeMiter,
	),
	fill=(
		nonzero=Clipper.PolyFillTypeNonZero,
		evenodd=Clipper.PolyFillTypeEvenOdd,
		positive=Clipper.PolyFillTypePositive,
	),
)#»»
# conversions««2

@inline toc(x::Real) = reinterpret(convert(_CLIPPER_FIXED, x))
@inline toc(x::StaticVector{2,<:Real}) = Clipper.IntPoint(toc.(x)...)
@inline toc(x) = toc.(x) # recursively...

@inline fromc(x::Int64) = Float64(reinterpret(_CLIPPER_FIXED,x))
@inline fromc(x::Clipper.IntPoint) = SA[fromc(x.X), fromc(x.Y)]
@inline fromc(x::StaticVector) = SA[fromc.(x)...]
@inline fromc(x) = fromc.(x)

# this is a special case; the `delta` parameter wants a Float64,
# and we *might* use a different type for `delta` than for coordinates:
# @inline to_clipper_float(T::Type{<:Real}, x, n=1)::Float64 = x*_CLIPPER_ONE^n

# numbers...
# Wrappers for Clipper calls««1
# We wrap all Clipper objects in a NamedTuple with the original type
# struct Marked{T,X}
# 	data::X
# 	@inline Marked{T}(x::X) where{T,X} = new{T,X}(x)
# end

# @inline ClipperClip(T::Type) = (T, Clipper.Clip())
# @inline ClipperClip(T::Type) = Marked{T}(Clipper.Clip())
# @inline ClipperOffset(T::Type, miterLimit::Real, roundPrecision::Real) =
# 	Marked{T}(Clipper.ClipperOffset(Float64(miterLimit),
# 		to_clipper_float(T, roundPrecision)))

# @inline add_path!(Type{<:Real}, c, path, args...) =
# 	Clipper.add_path!(c, to_clipper(T, path), args...)
# @inline add_paths!(T::Type{<:Real}, c, paths, args...) =
# 	Clipper.add_paths!(c, to_clipper(T, paths), args...)

# @inline execute(T::Type{<:Real}, c::Clipper.Clip, args...) =
# 	from_clipper(T, Clipper.execute(c, args...)[2])
# @inline execute(T::Type{<:Real}, c::Clipper.ClipperOffset, delta::Real) =
# 	from_clipper(T, Clipper.execute(c, to_clipper_float(T, delta)))

# Clipper calls on Path values««1
@inline function clip(op::Symbol,
		v1::AbstractVector{Path{2,T}},
		v2::AbstractVector{Path{2,T}};
		fill = :nonzero)::Vector{Path{2,T}} where {T}
	c = Clipper.Clip()
	Clipper.add_paths!(c, toc(v1), Clipper.PolyTypeSubject, true) # closed=true
	Clipper.add_paths!(c, toc(v2), Clipper.PolyTypeClip, true)

	f = _CLIPPER_ENUM.fill[fill]
	return fromc(Clipper.execute(c, _CLIPPER_ENUM.clip[op], f, f)[2])
end
@inline function offset(v::AbstractVector{<:Path{2}}, r::Real;
		join = :round,
		ends = :fill,
		miter_limit = 2.,
		precision = 0.2
		)::Vector{Path{2,Float64}}
	c = Clipper.ClipperOffset(miter_limit, precision*_CLIPPER_ONE)
	Clipper.add_paths!(c, toc(v),
		_CLIPPER_ENUM.join[join], _CLIPPER_ENUM.ends[ends])
	return fromc(Clipper.execute(c, r*_CLIPPER_ONE))
end
@inline function offset(v::AbstractVector{<:Path{2}}, r::AbstractVector{<:Real};
		join = :round,
		ends = :fill,
		miter_limit = 2.,
		precision = 0.2
		)::Vector{Vector{Path{2,Float64}}} where{T}
	# “Simultaneously” computes offset for several offset values.
	# Used by path_extrude().
	c = Clipper.ClipperOffset(miter_limit, precision*_CLIPPER_ONE)
	Clipper.add_paths!(c, toc(v),
		_CLIPPER_ENUM.join[join], _CLIPPER_ENUM.ends[ends])
	return [ fromc(Clipper.execute(c, ρ*_CLIPPER_ONE)) for ρ in r]
end
@inline simplify(p::AbstractVector{<:Path{2}}; fill=:nonzero) =
	fromc(Clipper.simplify_polygons(toc(p), _CLIPPER_ENUM.fill[fill]))
"""
    orientation(p::Path{2})

Returns `true` iff p is a direct loop (i.e. if area >= 0).
"""
@inline orientation(p::Path{2}) = Clipper.orientation(toc(p))
@inline area(p::Path{2}) = Clipper.area(toc(p)) / _CLIPPER_ONE^2

"""
    point_in_polygon(pt, poly::Path{2})

Returns 1 if point is in the interior, -1 on boundary, and 0 outside the
given polygon.

Polygon is assumed not self-intersecting.
"""
@inline point_in_polygon(point::StaticVector{2}, path) =
	Clipper.pointinpolygon(toc(point), toc(path))

# # convex hull of 2d points (monotone chain)««1
# # this is the version using MiniQhull: #««
# # convex_hull(points::AbstractVector{<:AnyVec(2)}) =
# #		# Delaunay triangulation:
# #		let T = MiniQhull.delaunay([p[i] for i in 1:2, p in points]),
# #				N = length(points),
# #				M = zeros(Bool, N, N) # ⚠ memory O(N²)
# #		# mark all edges:
# #		for (a,b,c) in eachcol(T)
# #			b1 = points[b] - points[a]
# #			c1 = points[c] - points[a]
# #			d = b1[1]*c1[2] - b1[2]*c1[1] # determinant === orientation of triangle
# #			if d < 0
# #				M[a,b] = M[b,c] = M[c,a] = true
# #			else
# #				M[b,a] = M[c,b] = M[a,c] = true
# #			end
# #		end
# #		# list of remaining edges (retrograde oriented)
# #		L= sort([(i,j) for i in 1:N, j in 1:N if M[i,j] && ! M[j,i]], by=v->v[1])
# #		next(i) = L[searchsorted(L, i, by=y->y[1])][1][2]
# #		R = zeros(Int, length(L))
# #		R[1:2] .= L[1] # initialize with first edge
# #		for i in 3:length(R)
# #			R[i] = next(R[i-1])
# #		end
# #		# returns in retrograde ordering (OpenSCAD convention):
# #		points[R]
# # end#»»
# """
#     convex_hull_list(points)
#
# Returns the convex hull of the points, as a list of indexes (in direct
# order, starting at a reproducible index in the list of points).
# """
# function convex_hull_list(points::AbstractVector{<:StaticVector{2}})
# 	list = Int[]
# 	length(points) ≤ 2 && return points
# 	p = sortperm(points) # lex by x then y
# 	k = 1
# 	pdet = @closure (i,j,k) -> begin
# 		qi = points[i]; qj = points[j]; qk = points[k]
# 		vij = qj-qi; vik = qk-qi
# 		return vij[2]*vik[1]-vij[1]*vik[2]
# 	end
# 	for i in p # lower hull
# 		while length(list) ≥ 2 && pdet(list[end-1], list[end], i) ≥ 0
# 			pop!(list)
# 		end
# 		push!(list, i)
# 	end
# 	t = length(list)
# 	for i in reverse(p) # upper hull
# 		while length(list) ≥ t+1 && pdet(list[end-1], list[end], i) ≥ 0
# 			pop!(list)
# 		end
# 		push!(list, i)
# 	end
# 	pop!(list) # remove loop
# 	return list
# end
# """
#     convex_hull([vector of 2d points])
#
# Returns the convex hull (as a vector of 2d points, ordered in direct
# order).
# """
# @inline convex_hull(points::AbstractVector{<:StaticVector{2}}) =
# 	points[convex_hull_list(points)]
# Draw path««1
# # """
# #     draw(path, width; kwargs...)
# #
# #     ends=:round|:square|:butt|:closed
# #     join=:round|:miter|:square
# # """
# # function draw(path::Path{2,T}, width::Real;
# # 		ends::Symbol = :round, join::Symbol = :round,
# # 		miter_limit::Float64 = 2.0, precision::Real = 0.2) where{T}
# # 	CT = clipper_type(T)
# # 	RT = clipper_rettype(T)
# # 	c = ClipperOffset(miter_limit, clipper_float(CT, precision))
# # 	println("join=$join, round=$round")
# # 	Clipper.add_path!(c, clipper_path(path),
# # 		JoinType(Val(join)), EndType(Val(ends)))
# # 	println("$(clipper_type(T)) $(CT(1.)); prec=$(Float64(CT(precision)))")
# # 	ret = clipper_unpath.(RT, Clipper.execute(c, clipper_float(CT, width)/2))
# # 	return PolyUnion(ret)
# # end
# #
# PolygonXor««1
"""
    PolygonXor{T}(polygon...)
    PolygonXor{T}([points], [points]...)

Exclusive union of several polygons.
(This may be used to represent either the disjoint union of several polygons,
or a polygon with holes, etc.)
"""
struct PolygonXor{T}
	paths::Vector{Vector{SVector{2,T}}}
# 	@inline PolygonXor{T}(paths::AbstractVector{<:AbstractVector{<:StaticVector{2,<:Real}}}) where{T<:Real} =
# 		new{T}(Vector{SVector{2,T}}.(paths))
	@inline PolygonXor{T}(paths::AbstractVector) where{T} = new{T}(paths)
	@inline PolygonXor(paths::AbstractVector) =
		PolygonXor{eltype(eltype(eltype(paths)))}(paths)
end
@inline PolygonXor{T}(paths::AbstractVector{<:StaticVector{2,<:Real}}...
	) where{T} = PolygonXor{T}([paths...])

@inline coordtype(::Type{PolygonXor{T}}) where{T} = T
@inline coordtype(s::PolygonXor) = coordtype(typeof(s))
@inline paths(p::PolygonXor) = p.paths
@inline vertices(p::PolygonXor) = reduce(vcat, paths(p))
@inline simplify(p::P; fill=:nonzero) where{P<:PolygonXor} =
	P(simplify(paths(p); fill))

"""
    loops(::PolygonXor)

Returns the lists of indexes for each loop.
Perimeters are oriented ↺ and holes ↻.
"""
function loops(p::PolygonXor)
	firstindex = zeros(Int, length(p.paths))
	firstindex[1] = 0
	for i in 1:length(p.paths)-1
		firstindex[i+1] = firstindex[i] + length(p.paths[i])
	end
	return [firstindex[i]+1 : firstindex[i] + length(p.paths[i])
		for i in eachindex(p.paths)]
end

"""
    identify_polygons(s::PolygonXor)

Given a `PolygonXor` defined as the exclusive union of polygons and holes,
returns a vector mapping individual polygons of `s` to:
 - positive indices 1, 2, 3... for the exterior polygons of `s`;
 - negative indices for the holes, indicating in which perimeter polygon is the corresponding hole.
"""
function identify_polygons(s::PolygonXor)
	n = length(paths(s))
	m = [ point_in_polygon(first(p), q) for p in paths(s), q in paths(s) ]
	index = zeros(Int, n)
	hole = zeros(Int, n)
	c = 0 # number of connected components
	for i in 1:n
		# if this is a hole, then it will have m[i,j] > 0 for some i
		k = findfirst(m[i,:] .> 0)
		if k == nothing
			index[i] = (c+= 1)
		else
			hole[i] = k
		end
	end
	for i in 1:n
		if hole[i] > 0
			index[i] = -index[hole[i]]
		end
	end
	return index
end

@inline clip(op, p::P, q::P) where{P<:PolygonXor} =
	P(clip(op, paths(p), paths(q), fill=:nonzero))
@inline clip(op, s::P...) where{P<:PolygonXor} = reduce((p,q)->clip(op, p,q), s)

@inline offset(s::PolygonXor{T}, r; kwargs...) where{T} =
	PolygonXor{T}(offset(paths(s), r; kwargs...))

@inline area(shape::PolygonXor) = sum(area.(paths(shape)))
@inline convex_hull(shape::PolygonXor{T}...) where{T} =
	PolygonXor{T}(convex_hull([vertices.(shape)...;]))

# Minkowski sum of polygons and their unions ««1
"""
    segment_xray(p, a, b)

Returns intersection abscissa for segments (p, p+(∞,0)) and (a,b),
or `nothing`.
"""
function segment_xray(p,a,b)#««
	(a[2] > b[2]) && ((a,b) = (b,a)) # guarantees yb ≥ ya
	(x,y) = p; (xa,ya) = a; (xb,yb) = b
	(ya > y) && return nothing
	(yb < y) && return nothing # now ya ≤ y ≤ yb
	if ya == yb # both are == y
		(u,v) = minmax(xa,xb)
		return (v < x) ? nothing : max(u,x)
	end
	# ya < y < yb
	xc =  (xa*(yb-y)+xb*(y-ya))/(yb-ya)
	return (xc ≥ x) ? xc : nothing
end#»»
"""
    polygon_xray(point, loop)

Returns the segment in `loop` intersecting the ray originating
from `point` in the ``(+∞,0)`` direction. The segment is returned as
a pair (destination point, x).
"""
function polygon_xray(point, loop)#««
	x = nothing; k = nothing
	b = last(loop)
	for (i, a) in pairs(loop)
		v = segment_xray(point, a, b)
		v == nothing && @goto next
		x == nothing && @goto found
		x < v && @goto next
		@label found
		x = v; k = i
		@label next
		b = a
	end
	return (k, x)
end#»»

"""
    connect_holes(s::PolygonXor)

Returns a representation of `s` as a disjoint union of pseudo-simple loops
(the holes of each polygon being joined to its perimeter via a short “cut”).
"""
function connect_holes(s::PolygonXor)#««
	poly_id = identify_polygons(s); n = length(poly_id)
	loops = [ [i] for (i,a) in pairs(poly_id) if a > 0 ]
	for (i, a) in pairs(poly_id); a < 0 && push!(loops[-a], i); end
	# loops[i] is (perimeter, hole1, hole2, ...)

	newpaths = [ copy(paths(s)[first(l)]) for l in loops ]
	for ((_, holes...), newpath) in zip(loops, newpaths)
		xmax = [ findmax(paths(s)[h]) for h in holes ]
		for i in sort(eachindex(holes); by=i->xmax[i][1], rev=true)
			h = holes[i]; hole = paths(s)[h]
			(p, j) = xmax[i] # j = index, p = xmax point
			(k, x) = polygon_xray(p, newpath)
			q = (x, p[2])
			(q ≠ newpath[k]) && insert!(newpath, k, q)
			# the hole path starts at point #j, so j:end followed by start:j
			splice!(newpath, k+1:k, [hole[j:end];hole[1:j];[q]])
		end
# 		println("# newpath ($(length(newpath))):")
# 		for p in newpath; println("$(p[1]+.2rand())\t$(p[2]+.2rand())"); end
	end
	# this is *not* a PolygonXor, because of all the degeneracies
	# rather, it must be converted back to that format with a Clipper call
	return newpaths
end#»»

function minkowski_sum(m1::PolygonXor, m2::PolygonXor)#««
	m1c = connect_holes(m1) # vector of Paths
	m2c = connect_holes(m2) # idem
	result = nothing
	for p1 in m1c, p2 in m2c
		# vector of paths
		s = fromc(Clipper.minkowski_sum(toc(p1), toc(p2)))
		s = PolygonXor(simplify(s; fill=:nonzero))
		result = (result == nothing) ? s : clip(:union, result, s);

		# fill holes:
		a1 = first(p1); a2 = first(p2)
		result = clip(:union, result,
			PolygonXor([b1+a2 for b1 in p1]),
			PolygonXor([b2+a1 for b2 in p2]))
	end
	return result
end#»»

# Triangulation and reconstruction from triangle««1
"""
    point_in_polygon(poly; orientation=orientation(poly))
Returns one point in the given `poly` (given as a list of points).
"""
function point_in_polygon(poly; orientation=orientation(poly))#««
	@assert length(poly) ≥ 3
	a,b,c = poly # first three points
	ac = c-a
	# compute min (squared) distance from b to other sides of the polygon:
	q = poly[end]; qa = a-q
	h2 = cross(qa,(a-b))^2 / dot(qa,qa)
	for i in 3:length(poly)-1
		p = poly[i]; q = poly[i+1]; pq = q-p
		h2 = min(h2, cross(pq, (p-b))^2 / dot(pq, pq))
	end

	# vector pointing inside the polygon:
	v = SA[-ac[2], ac[1]]
	t = √(h2/4/dot(ac,ac))
	return orientation ? b + t*v : b - t*v
end#»»
# triangulates a simple loop
function triangulate(v::AbstractVector{<:SVector{2,<:Real}})
	tri = basic_triangulation(
		Matrix{Float64}([transpose.(v)...;]),
		collect(1:length(v)))
	return orientation(v) ? [ (t[1], t[2], t[3]) for t in tri ] :
	 [(t[1], t[3], t[2]) for t in tri ]
end
function triangulate(m::PolygonXor)#««
# 	v = vertices(m)
	id = identify_polygons(m)
# 	edges = Matrix{Int}(undef, sum(length.(paths(m))), 2)
	tri = NTuple{3,Int}[]
	peri = loops(m)
	# triangulate each loop separately
	for (l, k) in pairs(id)
		k > 0 || continue
		# l is outer loop, h is list of holes
		h = filter((@closure j->id[j] == -k), eachindex(id))
		holes = [ point_in_polygon(p; orientation=false) for p in paths(m)[h] ]
		push!(h, k)
		plist = paths(m)[h]
		edges = Matrix{Int}(undef, sum(length.(plist)), 2)
		c = 0
		labels = [peri[h]...;]
		for p in plist
			n = length(p)
			for i in 1:n-1
				edges[c+i,:] = labels[[c+i,c+i+1]]
			end
			edges[c+n,:] = labels[[c+n,c+1]]
			c+= n
		end
		v = [plist...;]
		triangles = constrained_triangulation(
			Matrix{Float64}([transpose.(v)...;]),
			labels, edges, fill(true, size(edges,1)),
			[ p[i] for p in holes, i in 1:2 ])
		tri = [tri; [ (t[1], t[2], t[3]) for t in triangles ]]
	end
	return tri
end#»»
# reconstruction from triangles
# this is used by 3d->2d projection:
function PolygonXor{T}(points::AbstractVector{<:SVector{2}},
	faces::AbstractVector{<:Tuple{<:Integer,<:Integer,<:Integer}}) where{T}
	v=[PolygonXor{T}([points[f[1]], points[f[2]], points[f[3]]]) for f in faces ]
	return clip(:union, v...)
end

# this is used by 3d->2d slicing:
function glue_segments(points::AbstractVector, segments)
	n = maximum(x[2] for x in segments)
	# compute the set of all edges
	edges = [ SortedSet{eltype(eltype(segments))}() for _ in 1:n ]
	for (a, b) in segments
		push!(edges[a], b)
		push!(edges[b], a)
	end
	# decomposes the graph as a set of loops
	loops = Vector{eltype(points)}[]
	for (a0, e) in pairs(edges)
		isempty(e) && continue
		a = a0; l = [a]
		while !isempty(e)
			b = first(e)
			delete!(edges[a], b)
			delete!(edges[b], a)
			(b == a0) && (length(l) > 1) && break
			push!(l, b)
			e = edges[b]
			a = b
		end
		push!(loops, points[l])
	end
	return PolygonXor{eltype(eltype(points))}(loops)
end

# Intersections (2d)««1

struct HalfPlane{T<:Real}
	# equation a*x+b = 0
	a::SVector{2,T}
	b::T
end
@inline (h::HalfPlane)(v::AbstractVector) = dot(h.a, v) + h.b

function Base.intersect(halfplane::HalfPlane,
		path::AbstractVector{<:StaticVector{2}})
	newpath = similar(path, 0)
	n = length(path)
	s = halfplane.(path) # vector of signs
	for i in 1:n
		j = i+1; (j > n) && (j = 1)
		(si, sj) = (s[i], s[j])
		si >= 0 && push!(newpath, path[i])
		if si*sj < 0
			newpoint = (s[j]*path[i]-s[i]*path[j])/(s[j]-s[i])
			(isempty(newpath) || last(newpath) ≠ newpoint) && push!(newpath, newpoint)
		end
	end
	return newpath
end

@inline function bound(s::PolygonXor)
	r = maximum(norm(v, Inf) for v in vertices(s))
	return r + one(r)
end

@inline function bounding_square(s::PolygonXor)
	r = bound(s)
	return [SA[r,r], SA[-r,r], SA[-r,-r], SA[r,-r]]
end

@inline function Base.intersect(halfplane::HalfPlane, s::PolygonXor)
	c = intersect(halfplane, bounding_square(s))
	return clip(:intersection, PolygonXor{coordtype(s)}(c), s)
end

# Path extrusion««1
# triangulate_between: triangulate between two parallel paths««2
"""
		triangulate_between(poly1, poly2, start1, start2)

Given two polygons `poly1` and `poly2`, both of them represented as a
vector of paths, and produced as offsets from a common path,
find a triangulation for the region between the two polygons.

This functions returns a pair `(triangulation, edge)`, where:

 - the triangulation is a vector of triples of integers,
where each point is represented by its index. Indices in `poly1` start at
value `start1`, and in `poly2` at `start2`.

 - the edge is a pair `(lastidx1, lastidx2)` corresponding to the last
	 points visited on each polygon. (this will be useful for closing the
	 extrusion).

"""
function triangulate_between(poly1, poly2, start1=1, start2=1)
	Big = typemax(eltype(eltype(eltype(poly1))))
	triangles = NTuple{3,Int}[]
	# head is the marker of current leading edge
	# headpoint[i] is the point marked to by head[i]
	# headidx is the new index for this marked point
	# status[i][j] is the number of last used point in j-th path of i-th poly
	head = [(1,1), (1,1)]
	headpoint = [poly1[1][1], poly2[1][1]]
	headidx = [start1, start2]
	status = zeros.(Int,length.((poly1, poly2)))
	# so far we used exactly one point on each side:
	status[1][1] = status[2][1] = 1

	# we need a way to convert (poly, path, index) to integer index««
	function first_indices(start::Int, l::Vector{Int})::Vector{Int}
		f = zeros.(Int, length(l))
		f[1] = start
		for i in 1:length(l)-1
			@inbounds f[i+1] = f[i] + l[i]
		end
		f
	end
	# firstindex[poly][path] is the first index for this path
	# firstindex[1][1] = start1
	# firstindex[1][2] = start1 + len(poly1[1]) etc.
	firstindex = (first_indices(start1, length.(poly1)),
								first_indices(start2, length.(poly2)))
	newindex(poly::Int, path::Int, index::Int)::Int =
		firstindex[poly][path] + index - 1
#»»
	# computing diagonal distances to find the smallest one:««
	distance(pt, path, i) =
		i > length(path) ? Big : distance²(pt, path[i])

	closest(pt, poly, status) =
		findmin([distance(pt, poly[i], status[i]+1) for i in eachindex(poly)])
#»»

	while true
		d1, i1 = closest(headpoint[2], poly1, status[1])
		d2, i2 = closest(headpoint[1], poly2, status[2])
		# if no more points are left, we return:
		(d1 == d2 == Big) && break

		if d1 < d2 # we append a point from poly1
			# add the triangle: head1, head2, newpoint
			s = status[1][i1] += 1
			newidx = newindex(1, i1, s)
			push!(triangles, (headidx[1], headidx[2], newidx))
			# update head1 to point to new point
			headidx[1] = newidx
			head[1] = (i1, s)
			headpoint[1] = poly1[i1][s]
		else
			# add the triangle: head1, head2, newpoint
			s = status[2][i2] += 1
			newidx = newindex(2, i2, s)
			push!(triangles, (headidx[1], headidx[2], newidx))
			# update head1 to point to new point
			headidx[2] = newidx
			head[2] = (i2, s)
			headpoint[2] = poly2[i2][s]
		end
	end
	(triangles, (headidx[1], headidx[2]))
end#»»
# path_extrude««
"""
		path_extrude(path, poly, options...)

Extrudes the given polygon (a path of points forming a simple loop)
along the given path. Both arguments are provided as a
`Vector{SVector{2}}`.

Returns a (points, triangulation) pair.
"""
function path_extrude(path::Path{2}, poly::Path{2};
		join=:round, closed=true, miter_limit=2.0, precision=0.2)
	# in kwargs: ends = closed ? :fill : :butt
	N = length(poly)
	# offset_path is a vector of vector of paths
	@assert closed == true "Open-path extrusion is currently not implemented by ClipperLib"
	offset_path = offset([path], [pt[1] for pt in poly ];
		join, miter_limit, ends=closed ? :fill : :single)
	new_points = SVector{3,eltype(eltype(eltype(offset_path)))}[]
	first_face = Int[]; last_ff = 1
	for (a, paths) in zip(poly, offset_path)
		push!(first_face, last_ff)
		for p in paths
			pxyz = [SA[b[1], b[2], a[2]] for b in p]
			new_points = [new_points; pxyz]
			last_ff+= length(p)
		end
	end
	# new_points is a flat list of all 3d points produced
# 	new_points = [[
# 		[ SA[[pt[1], pt[2], poly[i][2]]] for pt in [p...;] ]
# 		for (i, p) in pairs(offset_path)
# 	]...;]

	# first index for each path
# 	first_face = cumsum([1; # initial
# 		map(p->sum(length.(p)), offset_path)])
# # 	println("first_face=$first_face")

	triangles = map(1:N) do i
		i1 = mod1(i+1, N)
		triangulate_between(offset_path[i], offset_path[i1],
			first_face[i], first_face[i1])
		# XXX keep the last edge for closing the poly
	end
	# this completes the set of triangles for the tube:
	tube_triangles = vcat([ t[1] for t in triangles ]...)
	last_face = [ t[2][1] for t in triangles ]
# 	println("last_face=$last_face")
	# here we decide if it is closed or open
	# if open, triangulate the two facets
	# if closed, join them together
	if closed
		more_triangles = vcat(map(1:N) do i
			j = (i%N)+1
			[ SA[first_face[i], last_face[i], first_face[j]],
				SA[first_face[j], last_face[i], last_face[j]] ]
		end...)
# 		println("more_triangles=$more_triangles")
		tube_triangles = [ tube_triangles; more_triangles ]
	else
	# TODO: triangulate the surface
	# or, for now, close with two non-triangular facets...
		more_triangles = [ reverse(first_face), last_face ]
# 		println("more_triangles=$more_triangles")
	end
	return (new_points, reverse.(tube_triangles))
end#»»
# Exports««1

# export Path
# export PolygonXor
end


text:
function axis_angle_to_quaternion(x)
    @assert length(x) == 3
    θ = norm(x)
    if θ > 0.0
        r = x ./ θ
        q = Quaternion(cos(0.5 * θ), sin(0.5 * θ) * r)
    else
        q = Quaternion(1.0, 0.0, 0.0, 0.0)
    end
    return q
end

function daxis_angle_to_quaterniondx(x)
    θ = norm(x)
    if θ > 0.0
        r = x ./ θ

        ∂qw∂x = -0.5 * sin(0.5 * θ) * transpose(x) ./ θ
        ∂qx∂x = 0.5 * cos(0.5 * θ) * transpose(x) ./ θ * r[1] + [sin(0.5 * θ) / θ 0.0 0.0] - sin(0.5 * θ) * x[1] / θ^2 * transpose(x) ./ θ
        ∂qy∂x = 0.5 * cos(0.5 * θ) * transpose(x) ./ θ * r[2] + [0.0 sin(0.5 * θ) / θ 0.0] - sin(0.5 * θ) * x[2] / θ^2 * transpose(x) ./ θ
        ∂qz∂x = 0.5 * cos(0.5 * θ) * transpose(x) ./ θ * r[3] + [0.0 0.0 sin(0.5 * θ) / θ] - sin(0.5 * θ) * x[3] / θ^2 * transpose(x) ./ θ

        return [
                ∂qw∂x;
                ∂qx∂x;
                ∂qy∂x;
                ∂qz∂x;
               ]
    else
        return [
                    0.0  0.0  0.0;
                    0.5  0.0  0.0;
                    0.0  0.5  0.0;
                    0.0  0.0  0.5;
                ]
    end
end

function axes_pair_to_quaternion(n1, n2)
	if norm(n1 + n2, Inf) < 1e-5
		n2 = n2 + 1e-5ones(3)
	end

	reg(x) = 1e-20 * (x == 0) + x
	# provides the quaternion that rotates n1 into n2, assuming n1 and n2 are normalized
	n1 ./= reg(norm(n1))
	n2 ./= reg(norm(n2))
	n3 = skew(n1)*n2
	cθ = n1' * n2 # cosine
	sθ = norm(n3) # sine
	axis = n3 ./ reg(sθ)
	tanθhalf = sθ / reg(1 + cθ)
	q = [1; tanθhalf * axis]
	q /= norm(q)
	return Quaternion(q...)
end


text:
module TestTwiddle

using Twiddle, Base.Test

@testset "Enumerating nibbles" begin
    @test Twiddle.enumerate_nibbles(0x4216CEDF4216CEDF) == 0x1112233411122334
end

@testset "Counting nibbles" begin
    @test Twiddle.count_nonzero_nibbles(0x0F11F111F11111F1) == 15
    @test Twiddle.count_zero_nibbles(0x0F11F111F11111F1) == 1
    @test Twiddle.count_one_nibbles(0x0F11F111F11111F1) == 4
end

@testset "Counting bitpairs" begin
    @test Twiddle.count_zero_bitpairs(0x185DF69C185DF69C) == 6
    @test Twiddle.count_nonzero_bitpairs(0x185DF69C185DF69C) == 26
    @test Twiddle.count_one_bitpairs(0x185DF69C185DF69C) == 8
end

@testset "Masking nibbles" begin
    @test Twiddle.nibble_mask(0x185DF69C185DF69C, 0xDDDDDDDDDDDDDDDD) == 0x000F0000000F0000
end

end


text:
#=

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

=#
using Distributed
addprocs(3)

max_num = Int(2e6)

@everywhere is_prime(n) = !any(n%i == 0 for i in 2:BigInt(floor(sqrt(n))))
partial_sum(start, stop) = sum(i for i in start:stop if is_prime(i))

res = @distributed (+) for i = 2:max_num
    if is_prime(i)
        i
    else
        0
    end
end

# sieve of erathostenes
# soooo much faster
max_num =  Int(2e6)
prime_candidates = trues(max_num)
for i in 2:max_num
    if prime_candidates[i]
        prime_candidates[i^2:i:max_num] .= false
    end
end

sum(findall(prime_candidates))-1

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DensityRatioValidation(k; [parameters])

Desntity ratio validation where weights are first obtained with density
ratio estimation, and then used in `k`-fold weighted cross-validation.

## Parameters

* `shuffle`   - Shuffle the data before folding (default to `true`)
* `estimator` - Density ratio estimator (default to `LSIF()`)
* `optlib`    - Optimization library (default to `default_optlib(estimator)`)
* `lambda`    - Power of density ratios (default to `1.0`)

Please see [DensityRatioEstimation.jl]
(https://github.com/JuliaEarth/DensityRatioEstimation.jl)
for a list of supported estimators.

## References

* Hoffimann et al. 2020. [Geostatistical Learning: Challenges and Opportunities]
  (https://arxiv.org/abs/2102.08791)
"""
struct DensityRatioValidation{T,E,O} <: ErrorEstimationMethod
  k::Int
  shuffle::Bool
  lambda::T
  dre::E
  optlib::O
  loss::Dict{Symbol,SupervisedLoss}
end

function DensityRatioValidation(k::Int; shuffle=true, lambda=1.0, loss=Dict(),
                                estimator=LSIF(), optlib=default_optlib(estimator))
  @assert k > 0 "number of folds must be positive"
  @assert 0 ≤ lambda ≤ 1 "lambda must lie in [0,1]"
  T = typeof(lambda)
  E = typeof(estimator)
  O = typeof(optlib)
  DensityRatioValidation{T,E,O}(k, shuffle, lambda, estimator, optlib, loss)
end

function Base.error(solver::LearningSolver, problem::LearningProblem,
                    method::DensityRatioValidation)
  tdata = targetdata(problem)
  vars = collect(inputvars(task(problem)))

  # density-ratio weights
  weighting = DensityRatioWeighting(tdata, vars,
                                    estimator=method.dre,
                                    optlib=method.optlib)

  # random folds
  folding = UniformFolding(method.k, method.shuffle)

  wcv = WeightedCrossValidation(weighting, folding,
                                lambda=method.lambda,
                                loss=method.loss)

  error(solver, problem, wcv)
end


text:
const soh     = Libdl.dlopen(Pkg.dir("IScatterSpectrum")*"/libiscatspe")
const plascb  = Libdl.dlsym(soh, :plas_)
const plasmaf = Libdl.dlsym(soh, :plasma_)
const pf32    = cglobal(plascb, Float32)
const pi32    = cglobal(plascb, Int32)

zr     = unsafe_wrap(Vector{Float32}, pf32, 100)
zi     = unsafe_wrap(Vector{Float32}, pf32+4*100, 100)
x      = unsafe_wrap(Vector{Float32}, pf32+4*200, 100)
y      = unsafe_wrap(Vector{Float32}, pf32+4*300, 1)
f      = unsafe_wrap(Vector{Float32}, pf32+4*301, 120)
scalef = unsafe_wrap(Vector{Float32}, pf32+4*421, 1)
nx     = unsafe_wrap(Vector{Int32}, pi32+4*422, 1)

function plasma!(freq::Vector{Float32})
    for (k, f) in enumerate(freq)
        unsafe_store!(pf32, f, 301+k)
    end
    unsafe_store!(pi32, Int32(length(freq)), 423)
    return
end

function plasma!(scalef::Float32=1f0, y::Float32=0f0)
    unsafe_store!(pf32,      y, 301)
    unsafe_store!(pf32, scalef, 422)
    ccall(plasmaf, Void, ())
    nx = unsafe_load(pi32, 423)
    hcat([unsafe_load(pf32, k) for k in 1:nx],
         [unsafe_load(pf32, k+100) for k in 1:nx])
end


text:
int_rules_7_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.1.1*(a+b*arcsinh(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~(n'), ~x) => ~x * (~a + ~b * asinh(~c * ~x)) ^ ~n - ~b * ~c * ~n * Antiderivative((~x * (~a + ~b * asinh(~c * ~x)) ^ (~n - 1)) / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => (sqrt(1 + (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (~c / (~b * (~n + 1))) * Antiderivative((~x * (~a + ~b * asinh(~c * ~x)) ^ (~n + 1)) / sqrt(1 + (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * asinh(~(c') * ~x)) ^ ~n, ~x) => (1 / (~b * ~c)) * Subst(Antiderivative((~x) ^ ~n * cosh(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * asinh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/7.1.1 (a+b arcsinh(c x))^n.jl


text:
# A number of data tables, mainly used for coordinate conversion

const x2pix = [
       0,     1,     4,     5,    16,    17,    20,    21,    64,    65,
      68,    69,    80,    81,    84,    85,   256,   257,   260,   261,
     272,   273,   276,   277,   320,   321,   324,   325,   336,   337,
     340,   341,  1024,  1025,  1028,  1029,  1040,  1041,  1044,  1045,
    1088,  1089,  1092,  1093,  1104,  1105,  1108,  1109,  1280,  1281,
    1284,  1285,  1296,  1297,  1300,  1301,  1344,  1345,  1348,  1349,
    1360,  1361,  1364,  1365,  4096,  4097,  4100,  4101,  4112,  4113,
    4116,  4117,  4160,  4161,  4164,  4165,  4176,  4177,  4180,  4181,
    4352,  4353,  4356,  4357,  4368,  4369,  4372,  4373,  4416,  4417,
    4420,  4421,  4432,  4433,  4436,  4437,  5120,  5121,  5124,  5125,
    5136,  5137,  5140,  5141,  5184,  5185,  5188,  5189,  5200,  5201,
    5204,  5205,  5376,  5377,  5380,  5381,  5392,  5393,  5396,  5397,
    5440,  5441,  5444,  5445,  5456,  5457,  5460,  5461 ]

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

const y2pix = [
        0,     2,     8,    10,    32,    34,    40,    42,   128,   130,
      136,   138,   160,   162,   168,   170,   512,   514,   520,   522,
      544,   546,   552,   554,   640,   642,   648,   650,   672,   674,
      680,   682,  2048,  2050,  2056,  2058,  2080,  2082,  2088,  2090,
     2176,  2178,  2184,  2186,  2208,  2210,  2216,  2218,  2560,  2562,
     2568,  2570,  2592,  2594,  2600,  2602,  2688,  2690,  2696,  2698,
     2720,  2722,  2728,  2730,  8192,  8194,  8200,  8202,  8224,  8226,
     8232,  8234,  8320,  8322,  8328,  8330,  8352,  8354,  8360,  8362,
     8704,  8706,  8712,  8714,  8736,  8738,  8744,  8746,  8832,  8834,
     8840,  8842,  8864,  8866,  8872,  8874, 10240, 10242, 10248, 10250,
    10272, 10274, 10280, 10282, 10368, 10370, 10376, 10378, 10400, 10402,
    10408, 10410, 10752, 10754, 10760, 10762, 10784, 10786, 10792, 10794,
    10880, 10882, 10888, 10890, 10912, 10914, 10920, 10922 ]

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

const pix2x = [
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  0,  1,  0,  1,  2,  3,  2,  3,  0,  1,  0,  1,  2,  3,  2,  3,
  4,  5,  4,  5,  6,  7,  6,  7,  4,  5,  4,  5,  6,  7,  6,  7,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
  8,  9,  8,  9, 10, 11, 10, 11,  8,  9,  8,  9, 10, 11, 10, 11,
 12, 13, 12, 13, 14, 15, 14, 15, 12, 13, 12, 13, 14, 15, 14, 15,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 16, 17, 16, 17, 18, 19, 18, 19, 16, 17, 16, 17, 18, 19, 18, 19,
 20, 21, 20, 21, 22, 23, 22, 23, 20, 21, 20, 21, 22, 23, 22, 23,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31,
 24, 25, 24, 25, 26, 27, 26, 27, 24, 25, 24, 25, 26, 27, 26, 27,
 28, 29, 28, 29, 30, 31, 30, 31, 28, 29, 28, 29, 30, 31, 30, 31 ]

const pix2y = [
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  0,  0,  1,  1,  0,  0,  1,  1,  2,  2,  3,  3,  2,  2,  3,  3,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  4,  4,  5,  5,  4,  4,  5,  5,  6,  6,  7,  7,  6,  6,  7,  7,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
  8,  8,  9,  9,  8,  8,  9,  9, 10, 10, 11, 11, 10, 10, 11, 11,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
 12, 12, 13, 13, 12, 12, 13, 13, 14, 14, 15, 15, 14, 14, 15, 15,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 16, 16, 17, 17, 16, 16, 17, 17, 18, 18, 19, 19, 18, 18, 19, 19,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 20, 20, 21, 21, 20, 20, 21, 21, 22, 22, 23, 23, 22, 22, 23, 23,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 24, 24, 25, 25, 24, 24, 25, 25, 26, 26, 27, 27, 26, 26, 27, 27,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31,
 28, 28, 29, 29, 28, 28, 29, 29, 30, 30, 31, 31, 30, 30, 31, 31 ]


text:
# include all my files, and setup the environment -
include("Include.jl")

# main -
function main()

    # generate the problem dictionary -
    problem_dictionary = generate_problem_dictionary()

    # pass the problem dictionary to the GLPK wrapper -
    return calculate_optimal_flux_distribution(problem_dictionary)
end

# execute main -
results = main()


text:
#quickly export dataframe columns:

using JLD2, DataFrames, CSV, Dates
using RCall
include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

#using Distributions

## Need to have oligo_df_plus loaded! This has zeros replaced with small amount:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"

##
#y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)]))
y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)]))
sbp0=[1 -1 1 1 1 1; 1 0 1 -1 -1 1; 0 0 0 1 -1 0; 1 0 1 0 0 -1; 1 0 -1 0 0 0]
ilr_y, bal=ilr(y,sbp0)

wisp=findall(dayofyear.(oligo_df_plus[:,:date]) .< 166)
su=findall((dayofyear.(oligo_df_plus[:,:date]) .>= 166) .& (dayofyear.(oligo_df_plus[:,:date]) .< 258))
fa=findall(dayofyear.(oligo_df_plus[:,:date]) .>= 258)

test_fa = y[fa,:]
test_wisp = y[wisp,:]
test_su = y[su,:]
## Ternary Plots, courtesy of the lovely Ternary package in R:

R"require(compositions)"
#R"require(png)"

R"""
W_tr=t($sbp0)
V=gsi.buildilrBase(W_tr)
y=acomp($test_fa)
colnames(y) <- c('O1','O2','O3','O4','O5','O6')
png(file="/home/kristen/Documents/V6V8_analysis/paper_V6V8/EM_submission/test_dendro_fa.png")
CoDaDendrogram(y,V)
dev.off()
"""
##

R"""
W_tr=t($sbp0)
V=gsi.buildilrBase(W_tr)
y=acomp($test_wisp)
colnames(y) <- c('O1','O2','O3','O4','O5','O6')
png(file="/home/kristen/Documents/V6V8_analysis/paper_V6V8/EM_submission/test_dendro_wisp.png")
CoDaDendrogram(y,V)
dev.off()
"""
##

##
R"""
W_tr=t($sbp0)
V=gsi.buildilrBase(W_tr)
y=acomp($test_su)
colnames(y) <- c('O1','O2','O3','O4','O5','O6')
png(file="/home/kristen/Documents/V6V8_analysis/paper_V6V8/EM_submission/test_dendro_su.png")
CoDaDendrogram(y,V)
dev.off()
"""

#break up year into seasons:

# X_wisp = convert(Array{Union{Float64,Missing}},oligo_df_plus[wisp,vars])
# Y_wisp = ilr_y[wisp,:]
# #jj=findall(ismissing.(X_wisp)) #-> this should be zero

# X_su = convert(Array{Union{Float64,Missing}},oligo_df_plus[su,vars])
# Y_su = ilr_y[su,:]
# jj=findall(ismissing.(X_su))
# qq=unique(map(x->x[1], jj)) #no nutrient values for one of the samples, exclude
# ii=setdiff(collect(1:size(X_su,1)),qq)
# X_su=X_su[ii,:]
# Y_su=Y_su[ii,:]