text:
module PLX
using ..Struct, ..Common, PiecewiseIncreasingRanges
import Compat.String

export PLXSpikeChannel, PLXEventChannel, PLXContinuousChannel, PLXFile, readplx

const MagicNumber = 0x58454c50

function struct_string(fn::Function, bytes::Vector{UInt8})
    last_byte = findfirst(bytes, 0)
    fn(last_byte == 0 ? bytes : bytes[1:last_byte-1])
end

@struct PL_FileHeader begin
    MagicNumber::UInt32
    Version::Int32
    Comment::String(128)
    ADFrequency::Int32
    NumDSPChannels::Int32
    NumEventChannels::Int32
    NumSlowChannels::Int32
    NumPointsWave::Int32
    NumPointsPreThr::Int32

    Year::Int32
    Month::Int32
    Day::Int32
    Hour::Int32
    Minute::Int32
    Second::Int32

    FastRead::Int32
    WaveformFreq::Int32
    LastTimestamp::Float64

    Trodalness::UInt8
    DataTrodalness::UInt8
    BitsPerSpikeSample::UInt8
    BitsPerSlowSample::UInt8
    SpikeMaxMagnitudeMV::UInt16
    SlowMaxMagnitudeMV::UInt16
    SpikePreAmpGain::UInt16

    Padding::UInt8(46)

    TSCounts::Int32(5, 130)
    WFCounts::Int32(5, 130)

    EVCounts::Int32(512)
end

@struct PL_ChanHeader begin
    Name::String(32)
    SIGName::String(32)
    Channel::Int32
    WFRate::Int32
    SIG::Int32
    Ref::Int32
    Gain::Int32
    Filter::Int32
    Threshold::Int32
    Method::Int32
    NUnits::Int32
    Template::Int16(64, 5)
    Fit::Int32(5)
    SortWidth::Int32
    Boxes::Int16(4, 2, 5)
    SortBeg::Int32
    Comment::String(128)
    Padding::Int32(11)
end

@struct PL_EventHeader begin
    Name::String(32)
    Channel::Int32
    Comment::String(128)
    Padding::Int32(33)
end

@struct PL_SlowChannelHeader begin
    Name::String(32)
    Channel::Int32
    ADFreq::Int32
    Gain::Int32
    Enabled::Int32
    PreAmpGain::Int32

    SpikeChannel::Int32

    Comment::String(128)
    Padding::Int32(28)
end

type PLXSpikeChannel <: SpikeChannel
    header::PL_ChanHeader
    sample_rate::Int32
    voltage_multiplier::Float64
    bits_per_sample::Int8
    times::Vector{Float64}
    unit_numbers::Vector{Int16}
    data::Matrix{Int16}

    function PLXSpikeChannel(header::PL_ChanHeader, fheader::PL_FileHeader)
        voltage_multiplier = 1e-3*(if fheader.Version >= 105
                fheader.SpikeMaxMagnitudeMV/
                    ((1 << (fheader.BitsPerSpikeSample-1))*header.Gain*fheader.SpikePreAmpGain)
            elseif fheader.Version >= 103
                fheader.SpikeMaxMagnitudeMV/
                    ((1 << (fheader.BitsPerSpikeSample-1))*header.Gain*1000)
            else
                3000/(2048*header.Gain*1000)
            end)
        return new(header, fheader.ADFrequency, voltage_multiplier, fheader.BitsPerSpikeSample, Float64[], Int16[])
    end
end

type PLXEventChannel <: EventChannel
    header::PL_EventHeader
    times::Vector{Float64}
    data::Vector{Int16}

    PLXEventChannel(header::PL_EventHeader) = new(header, Float64[])
end

type PLXContinuousChannel <: ContinuousChannel
    header::PL_SlowChannelHeader
    voltage_multiplier::Float64
    data::Vector{Int16}
    times::PiecewiseIncreasingRange{Float64,StepRange{Int64,Int64},Int32}

    function PLXContinuousChannel(header::PL_SlowChannelHeader, fheader::PL_FileHeader)
        voltage_multiplier = 1e-3*(if fheader.Version >= 103
                fheader.SlowMaxMagnitudeMV/
                    ((1 << (fheader.BitsPerSpikeSample-1))*header.Gain*header.PreAmpGain)
            elseif fheader.Version == 102
                5000/(2048*header.Gain*header.PreAmpGain)
            else
                5000/(2048*1000*header.Gain)
            end)

        return new(header, voltage_multiplier)
    end
end

Common.label(ch::Union{PLXSpikeChannel,PLXEventChannel,PLXContinuousChannel}) =
    ch.header.Name

immutable PLXFile <: File
    header::PL_FileHeader
    spike_channels::Channels{PLXSpikeChannel}
    event_channels::Channels{PLXEventChannel}
    continuous_channels::Channels{PLXContinuousChannel}
end
PLXFile(header::PL_FileHeader) =
    PLXFile(header, Channels{PLXSpikeChannel}(), Channels{PLXEventChannel}(),
            Channels{PLXContinuousChannel}())

function readplx(ios::IOStream; lfps::Bool=true, waveforms::Bool=true)
    const maxChannels = 32*1024

    fheader = read(ios, PL_FileHeader)
    if fheader.MagicNumber != MagicNumber
        error("$file_name does not appear to be a PLX file")
    elseif fheader.NumDSPChannels < 0 || fheader.NumDSPChannels > maxChannels ||
            fheader.NumEventChannels < 0 || fheader.NumEventChannels > maxChannels ||
            fheader.NumSlowChannels < 0 || fheader.NumSlowChannels > maxChannels
        error("PLX file header specifies an invalid number of channels")
    elseif fheader.ADFrequency <= 0
        error("PLX file header specifies ADFrequency <= 0")
    end
    x = PLXFile(fheader)

    for i = 1:x.header.NumDSPChannels
        header = read(ios, PL_ChanHeader)
        x.spike_channels[header.Channel] = PLXSpikeChannel(header, fheader)
    end
    if x.header.NumDSPChannels > 0 && fheader.NumPointsWave <= 0
        error("PLX file header specifies no points in waveform")
    end

    for i=1:x.header.NumEventChannels
        header = read(ios, PL_EventHeader)
        x.event_channels[header.Channel] = PLXEventChannel(header)
    end

    for i=1:x.header.NumSlowChannels
        header = read(ios, PL_SlowChannelHeader)
        lfps && (x.continuous_channels[header.Channel+1] = PLXContinuousChannel(header, fheader))
    end

    data_offset = position(ios)

    # Read through the file once to determine how much memory to allocate
    nspikes = zeros(Int, endof(x.spike_channels))
    nevents = zeros(Int, endof(x.event_channels))
    ntimestamps = zeros(Int, endof(x.continuous_channels))
    nsamples = zeros(Int, endof(x.continuous_channels))

    nblocks = 0
    max_offset = div(filesize(ios)-data_offset, 2)
    contents = Mmap.mmap(ios, Vector{Int16}, max_offset)::Vector{Int16}
    cur_offset = 1
    while cur_offset < max_offset
        block_type = contents[cur_offset]
        ch = contents[cur_offset+4]

        if block_type == 1          # spike
            nspikes[ch] += 1
        elseif block_type == 4      # event
            nevents[ch] += 1
        elseif block_type == 5      # continuous
            if lfps
                ntimestamps[ch+1] += 1
                nsamples[ch+1] += contents[cur_offset+7]
            end
        else
            error("invalid data block type ", block_type)
        end

        cur_offset += contents[cur_offset+7]+8
        nblocks += 1
    end

    # Allocate
    for (i, channel) = enumerate(x.spike_channels)
        resize!(channel.times, nspikes[i])
        resize!(channel.unit_numbers, nspikes[i])
        if waveforms
            channel.data = Array(Int16, fheader.NumPointsWave, nspikes[i])
        end
    end
    for (i, channel) = enumerate(x.event_channels)
        resize!(channel.times, nevents[i])
        if i == 257
            channel.data = Array(Int16, nevents[i])
        end
    end

    if lfps
        sample_dts = Array(Int64, length(ntimestamps))
        timestamp_ranges = Array(Vector{StepRange{Int64,Int64}}, length(ntimestamps))
        for (i, channel) = enumerate(x.continuous_channels)
            sample_dt = div(fheader.ADFrequency, channel.header.ADFreq)
            if sample_dt * channel.header.ADFreq != fheader.ADFrequency
                warn("""
                    channel $i frequency $(channel.header.ADFreq) is non-integer
                    multiple of AD frequency $(x.header.ADFrequency)
                """)
            end
            sample_dts[i] = sample_dt
            timestamp_ranges[i] = Array(StepRange{Int64,Int64}, ntimestamps[i])
            channel.data = Array(Int16, nsamples[i])
        end
    else
        sample_dts = Array(Int64, 0)
        timestamp_ranges = Array(Vector{StepRange{Int64,Int64}}, 0)
    end

    # Read through the file again
    cur_spike = zeros(Int, size(nspikes))
    cur_event = zeros(Int, size(nevents))
    cur_timestamp = zeros(Int, size(ntimestamps))
    cur_sample = zeros(Int, size(nsamples))
    cur_offset = 1
    while cur_offset < max_offset
        block_type = contents[cur_offset]
        timestamp = (convert(Int64, reinterpret(UInt16, contents[cur_offset+1])) << 32) +
            reinterpret(UInt16, contents[cur_offset+2]) +
            (convert(Int64, reinterpret(UInt16, contents[cur_offset+3])) << 16)
        ch = convert(Int, contents[cur_offset+4])

        if block_type == 1          # spike
            let channel = x.spike_channels[ch]
                num = (cur_spike[ch] += 1)

                channel.unit_numbers[num] = contents[cur_offset+5]
                channel.times[num] = timestamp/x.header.ADFrequency

                block_samples = contents[cur_offset+6]*contents[cur_offset+7]
                if waveforms
                    wf = channel.data
                    if block_samples == 0
                        for i = 1:size(wf, 1)
                            wf[i, num] = 0
                        end
                    else
                        for i = 1:block_samples
                            wf[i, num] = contents[cur_offset+7+i]
                        end
                    end
                end
                cur_offset += block_samples
            end
        elseif block_type == 4      # event
            let channel = x.event_channels[ch]
                num = (cur_event[ch] += 1)

                channel.times[num] = timestamp/x.header.ADFrequency
                if ch == 257
                    channel.data[num] = contents[cur_offset+5]
                end
            end
        elseif block_type == 5      # continuous
            block_samples = contents[cur_offset+7]
            if block_samples > 0
                if lfps
                    t = (cur_timestamp[ch+1] += 1)
                    channel = x.continuous_channels[ch+1]
                    num = cur_sample[ch+1]

                    sample_dt = sample_dts[ch+1]
                    timestamp_ranges[ch+1][t] = timestamp:sample_dt:timestamp+sample_dt*(block_samples-1)

                    start_offset = cur_offset + 7
                    samples = channel.data
                    for i = 1:block_samples
                        samples[num+i] = contents[start_offset+i]
                    end
                    cur_sample[ch+1] += block_samples
                end
                cur_offset += block_samples
            end
        end
        cur_offset += 8
    end

    @assert cur_spike == nspikes
    @assert cur_event == nevents
    @assert cur_timestamp == ntimestamps
    @assert cur_sample == nsamples

    if lfps
        for (i, channel) = enumerate(x.continuous_channels)
            channel.times = PiecewiseIncreasingRange(timestamp_ranges[i], x.header.ADFrequency)
        end
    end

    close(ios)

    return x
end

chmax(ch) = 2^(bitspersample(ch)-1)*voltagemultiplier(ch)*1e3

function PLXFile(;spike_channels=Channels{PLXSpikeChannel}(),
                  event_channels=Channels{PLXEventChannel}(),
                  continuous_channels=Channels{PLXContinuousChannel}())
    ADFrequency = 0
    NumPointsWave = 0
    LastTimestamp = 0.
    BitsPerSpikeSample = 16
    BitsPerSlowSample = 16
    SpikeMaxMagnitudeMV::UInt16 = 1
    SlowMaxMagnitudeMV::UInt16 = 1
    SpikePreAmpGain::UInt16 = 1000
    if !isempty(spike_channels)
        ADFrequency = round(Int32, samplerate(first(spike_channels)))
        for ch in spike_channels
            samplerate(ch) == ADFrequency || throw(ArgumentError("spike channels do not have matching sample rates"))
            if hasdata(ch)
                NumPointsWave = size(data(ch), 1)
            end
            LastTimestamp = max(LastTimestamp, last(times(ch)))
            mag = round(UInt16, chmax(ch)*SpikePreAmpGain)
            SpikeMaxMagnitudeMV = max(SpikeMaxMagnitudeMV, mag)
        end
    elseif !isempty(continuous_channels)
        ADFrequency = round(Int32, samplerate(first(continuous_channels)))
        for ch in continuous_channels
            ADFrequency = max(ADFrequency, samplerate(ch))
            LastTimestamp = max(LastTimestamp, last(times(ch)))
            mag = round(UInt16, chmax(ch))
            SlowMaxMagnitudeMV = max(SlowMaxMagnitudeMV, mag)
        end
    end
    LastTimestamp *= ADFrequency

    header = PL_FileHeader(MagicNumber, 105, "Created with NeuroIO.jl", ADFrequency,
                           length(spike_channels), length(event_channels), length(continuous_channels),
                           NumPointsWave, div(NumPointsWave, 2), 0, 0, 0, 0, 0, 0, 0, ADFrequency,
                           LastTimestamp, 1, 1, BitsPerSpikeSample, BitsPerSlowSample,
                           SpikeMaxMagnitudeMV, SlowMaxMagnitudeMV, SpikePreAmpGain, zeros(UInt8, 46),
                           zeros(Int32, 5, 130), zeros(Int32, 5, 130), zeros(Int32, 512))

    new_spike_channels = Channels{PLXSpikeChannel}()
    for (i, ch) in enumerate(spike_channels)
        if isa(ch, PLXSpikeChannel)
            new_spike_channels[i] = ch
            continue
        end
        Gain = round(Int32, SpikeMaxMagnitudeMV/(SpikePreAmpGain*chmax(ch)))
        NUnits = maximum(unitnumbers(ch))
        chheader = PL_ChanHeader(label(ch), "", i, 10, i, 0, Gain, 0, 0, 2, NUnits, zeros(Int16, 64, 5),
                                 zeros(Int32, 5), NumPointsWave, zeros(Int16, 4, 2, 5), 0, "", zeros(Int32, 11))
        newch = new_spike_channels[i] = PLXSpikeChannel(chheader, header)
        newch.times = times(ch)
        newch.unit_numbers = unitnumbers(ch)
        hasdata(ch) && (newch.data = data(ch))
    end

    new_event_channels = Channels{PLXEventChannel}()
    have_codes = false
    for (i, ch) in enumerate(event_channels)
        if isa(ch, PLXEventChannel)
            new_event_channels[i] = ch
            continue
        end

        if hasdata(ch) && !isempty(times(ch))
            if !have_codes
                i = 257
                have_codes = true
            else
                warn("multiple event channels have codes; ignoring codes for channel $i")
            end
        end
        chheader = PL_EventHeader(label(ch), i, "", zeros(Int332, 33))
        newch = new_event_channels[i] = PLXEventChannel(chheader)
        newch.times = times(ch)
        hasdata(ch) && i == 257 && (newch.data = data(ch))
    end

    new_continuous_channels = Channels{PLXContinuousChannel}()
    chs = Int[]
    for (i, ch) in enumerate(continuous_channels)
        hasdata(ch) && push!(chs, i)
    end
    dvecs = datavecs(continuous_channels, chs)

    idatavec = 0
    for (i, ch) in enumerate(continuous_channels)
        if isa(ch, PLXContinuousChannel)
            new_continuous_channels[i] = ch
            continue
        end

        Gain = round(Int32, SlowMaxMagnitudeMV/chmax(ch))
        chheader = PL_SlowChannelHeader(label(ch), i, samplerate(ch), Gain, hasdata(ch), 1, 0, "", zeros(Int32, 28))
        newch = new_continuous_channels[i] = PLXContinuousChannel(chheader, header)
        if hasdata(ch)
            newch.times = times(ch)
            newch.data = dvecs[idatavec += 1]
        end
    end

    validate_headers(PLXFile(header, new_spike_channels, new_event_channels, new_continuous_channels))
end

function validate_headers(x::PLXFile)
    x.header.NumDSPChannels = length(x.spike_channels)
    x.header.NumEventChannels = length(x.event_channels)
    x.header.NumSlowChannels = length(x.continuous_channels)

    nwfpoints = 0
    for (i, ch) in enumerate(x.spike_channels)
        length(ch.times) == length(ch.unit_numbers) ||
            throw(ArgumentError("number of spike times and unit numbers do not match for spike channel $i"))
        if isdefined(ch, :data)
            length(ch.times) == size(ch.data, 2) ||
                throw(ArgumentError("number of spike times and waveforms do not match for spike channel $i"))
            if nwfpoints == 0
                nwfpoints = size(ch.data, 1)
            elseif nwfpoints != size(ch.data, 1)
                throw(ArgumentError("all channels must have the same number of points per waveform"))
            end
        end

        ch.header.NUnits = maximum(ch.unit_numbers)
        ch.header.Channel = i
        x.header.NumPointsWave = nwfpoints

        if i <= 129
            counts = x.header.WFCounts
            for unumber in ch.unit_numbers
                unumber > 4 && continue
                counts[unumber+1, i] += 1
            end
        end
    end

    for (i, ch) in enumerate(x.event_channels)
        if isdefined(ch, :data)
            if i == 257
                length(ch.times) == length(ch.data) ||
                    throw(ArgumentError("number of event times and codes do not match"))
            else
                throw(ArgumentError("event codes were present on channel $i; only channel 257 may have associated event codes"))
            end
        end
        ch.header.EVCounts[300+i] = length(ch.times)
        ch.header.Channel = i
    end

    for (i, ch) in enumerate(x.continuous_channels)
        length(ch.times) == length(ch.data) ||
            throw(ArgumentError("number of samples and times do not match for spike $i"))
        ch.header.Channel = i-1
    end
    x
end

function integer_timestamp(time::Float64, ADFrequency::Int32)
    timestamp = round(Int64, time*ADFrequency)
    (Int16(timestamp >>> 32), timestamp % Int16, (timestamp >>> 16) % Int16)
end

const EMPTY_WAVEFORMS = Array(Int16, 0, 0)
const EMPTY_CODES = Int16[]
function Base.write(io::IO, x::PLXFile)
    validate_headers(x)

    write(io, x.header)
    for ch in x.spike_channels
        write(io, ch.header)
    end
    for ch in x.event_channels
        write(io, ch.header)
    end
    for ch in x.continuous_channels
        write(io, ch.header)
    end

    ADFrequency = x.header.ADFrequency
    for (ich, ch) = enumerate(x.spike_channels)
        has_waveforms = isdefined(ch, :data)
        waveforms = has_waveforms ? ch.data : EMPTY_WAVEFORMS
        waveform_size = has_waveforms ? size(waveforms, 1) : 0
        blockhead = Array(Int16, 8+waveform_size)
        blockhead[1] = 1
        blockhead[5] = ich
        blockhead[7] = has_waveforms
        blockhead[8] = waveform_size

        times = ch.times
        unit_numbers = ch.unit_numbers
        for i = 1:length(times)
            (blockhead[2], blockhead[3], blockhead[4]) = integer_timestamp(times[i], ADFrequency)
            blockhead[6] = unit_numbers[i]
            if has_waveforms
                for iwf = 1:waveform_size
                    blockhead[8+iwf] = waveforms[iwf, i]
                end
            end
            write(io, blockhead)
        end
    end

    blockhead = zeros(Int16, 8)
    for (ich, ch) = enumerate(x.event_channels)
        has_codes = isdefined(ch, :data)
        codes = has_codes ? ch.data : EMPTY_CODES
        blockhead[1] = 4
        blockhead[5] = ich
        if !has_codes
            blockhead[6] = 0
        end

        times = ch.times
        for i = 1:length(times)
            (blockhead[2], blockhead[3], blockhead[4]) = integer_timestamp(times[i], ADFrequency)
            has_codes && (blockhead[6] = codes[i])
            write(io, blockhead)
        end
    end

    blockhead[1] = 5
    blockhead[6] = 0
    blockhead[7] = 1
    for (ich, ch) = enumerate(x.continuous_channels)
        times = ch.times
        blockhead[5] = ich-1
        fac = ADFrequency/times.divisor
        n = 0
        for irg = 1:length(times.ranges)
            rg = times.ranges[irg]
            for irgstart = 1:typemax(Int16):length(rg)
                dur = min(typemax(Int16), length(rg)-irgstart+1)
                timestamp = round(Int, rg[irgstart]*fac)
                blockhead[2] = Int16(timestamp >>> 32)
                blockhead[3] = timestamp % Int16
                blockhead[4] = (timestamp >>> 16) % Int16
                blockhead[8] = dur
                write(io, blockhead)
                write(io, sub(ch.data, n+1:n+dur))
                n += dur
            end
        end
    end
end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export is_visited
"""
  `is_visited(fst::WFST, i = get_initial(fst;single=true))`

Return an array of booleans indicating if the i-th state of the fst is visted starting form `i`.
"""
function is_visited(fst::WFST, i = get_initial(fst;single=true))
  if isempty(fst)
    return Bool[]
  else
    d = DFS(fst,i)
    for i in d
      nothing
    end
    return d.is_visited
  end
end

ext: jl


text:
using NumCME
using Test
using LinearAlgebra: norm

## Bursting gene model
𝕊 = [[-1, 1, 0] [1, -1, 0] [0, 0, 1] [0, 0, -1]]
x₀ = [1, 0, 0]
k₀₁ = 0.05
k₁₀ = 0.1
λ = 5.0
γ = 1.0

θ = [k₀₁, k₁₀, λ, γ]

α₁ = propensity() do x, p
    p[1]*x[1]
end
α₂ = propensity() do x, p
    p[2]*x[2]
end
α₂tv = propensity((t,p) -> max(0.0, 1.0 - sin(π * t / 2))) do x, p
    p[2]*x[2]
end
α₂tvj = propensity() do t, x, p
    max(0.0, 1.0 - sin(π * t / 2))*p[2]*x[2]
end
α₃ = propensity() do x, p
    p[3]*x[2]
end
α₄ = propensity() do x, p
    p[4]*x[3]
end

propensities_ti = [α₁, α₂, α₃, α₄]
propensities_tv = [α₁, α₂tv, α₃, α₄]
propensities_tvj = [α₁, α₂tvj, α₃, α₄]

𝔛 = StateSpaceSparse(𝕊, x₀)
expand!(𝔛, 2)
𝐀 = FspMatrixSparse(𝔛, propensities_ti, parameters=θ)
@test size(𝐀, 1) == get_state_count(𝔛) + get_sink_count(𝔛)
@test size(𝐀, 2) == get_state_count(𝔛) + get_sink_count(𝔛)
𝐯 = ones(Float64, size(𝐀, 1))
𝐰 = 𝐀(1.0) * 𝐯
@test sum(𝐰) ≈ 0.0 atol = 1.0e-14
𝐰 = 𝐀 * 𝐯
@test sum(𝐰) ≈ 0.0 atol = 1.0e-14

# Test mat-vec for time-varying matrix
𝔛 = StateSpaceSparse(𝕊, x₀)
expand!(𝔛, 2)
A1 = FspMatrixSparse(𝔛, propensities_tv, parameters=θ)
@test size(A1, 1) == get_state_count(𝔛) + get_sink_count(𝔛)
@test size(A1, 2) == get_state_count(𝔛) + get_sink_count(𝔛)
𝐯 = ones(Float64, size(A1, 1))
w1 = A1(1.0) * 𝐯
@test sum(w1) ≈ 0.0 atol = 1.0e-14
w1 = A1 * 𝐯
@test sum(w1) ≈ 0.0 atol = 1.0e-14

A2 = FspMatrixSparse(𝔛, propensities_tvj, parameters=θ)
w2 = A2(1.0) * 𝐯
@test sum(w2) ≈ 0.0 atol = 1.0e-14
w2 = A2 * 𝐯

@test sum(w2) ≈ 0.0 atol = 1.0e-14
@test norm(w1 -w2) ≈ 0


text:
export OptM

abstract type PDFEstimator <: MLJBase.Unsupervised end

@with_kw struct OptM
    method::Union{Optim.AbstractOptimizer, Flux.Optimise.AbstractOptimiser} = NelderMead()
    maxiters::Int64 = Int64(typemax(Int8))
    atol::Float64 = Float64(eps(Float32))
end

function loss_4obj(logpx)
    lin_sca([
        -mean(logpx),
        -std(logpx),
        maximum(logpx),
        -minimum(logpx),
    ])
end

function loss_1obj(logpx)
    -mean(logpx)
end

long_optms = OptM[
    # - Flux
    OptM(
        method=AMSGrad(),
    ),
    # - Optim
    # -- Zeroth order
    # OptM(
    #     method=ParticleSwarm(),
    # ),
    OptM(
        method=SimulatedAnnealing(),
    ),
    OptM(
        method=NelderMead(),
    ),
    # -- First order
    # OptM(
    #     method=AcceleratedGradientDescent(),
    # ),
    # OptM(
    #     method=MomentumGradientDescent(),
    # ),
    # OptM(
    #     method=OACCEL(),
    # ),
    # OptM(
    #     method=NGMRES(),
    # ),
    OptM(
        method=ConjugateGradient(
            alphaguess=LineSearches.InitialHagerZhang(),
            linesearch=LineSearches.HagerZhang(),
            eta=1/2,
            manifold=Flat(),
        ),
    ),
    # OptM(
    #     method=GradientDescent(
    #         alphaguess=LineSearches.InitialHagerZhang(),
    #         linesearch=LineSearches.HagerZhang(),
    #         manifold=Flat(),
    #     ),
    # ),
    # OptM(
    #     method=LBFGS(
    #         alphaguess=LineSearches.InitialHagerZhang(),
    #         linesearch=LineSearches.HagerZhang(),
    #         manifold=Flat(),
    #     ),
    # ),
    OptM(
        method=BFGS(
            alphaguess=LineSearches.InitialHagerZhang(),
            linesearch=LineSearches.HagerZhang(),
            manifold=Flat(),
        ),
    ),
    # -- Second order
    OptM(
        method=Optim.KrylovTrustRegion(
            initial_radius=1.0,
            max_radius=128.0,
            eta=1/8,
            rho_lower=1/4,
            rho_upper=3/4,
            cg_tol=Float64(eps(Float32)),
        ),
    ),
    OptM(
        method=NewtonTrustRegion(
            initial_delta=1.0,
            delta_hat=128.0,
            eta=1/8,
            rho_lower=1/4,
            rho_upper=3/4,
            use_fg=true,
        ),
    ),
    OptM(
        method=Newton(
            alphaguess=LineSearches.InitialHagerZhang(),
            linesearch=LineSearches.HagerZhang(),
        ),
    ),
]
short_optms = OptM[
    OptM(
        method=ConjugateGradient(
            alphaguess=LineSearches.InitialHagerZhang(),
            linesearch=LineSearches.HagerZhang(),
            eta=1/2,
            manifold=Flat(),
        ),
    ),
]


text:
function pixelshuffle(x::AbstractArray{T, 4}, ratios = (2, 2)) where T
    iH, iW, C, B = size(x)
    pH, pW = ratios
    oC, oH, oW = C ÷ (pH*pW), iH * pH, iW * pW

    x = reshape(x, iH, iW, pH, pW, oC, B)
    x = permutedims(x, (4, 1, 3, 2, 5, 6))  # pH, iH, pW, iW, oC, B
    x = reshape(x, oH, oW, oC, B)

    return x
end

"""
    PixelShuffle(scale, kernels_in[, kernels_out])

Pixel shuffle layer that upscales height and width of `x` by `scale`. Has reduced
checkerboard artifacts compared to `ConvTranspose`

Introduced in [Real-Time Single Image and Video Super-Resolution Using
an EfficientSub-Pixel Convolutional Neural Network](https://arxiv.org/pdf/1609.05158.pdf).

"""
struct PixelShuffle
    conv::Conv
    scales::Tuple{Int, Int}
end

Flux.@functor PixelShuffle

function PixelShuffle(scales::Tuple{Int, Int}, k_in, k_out = k_in)
    return PixelShuffle(
        Conv((1, 1), k_in => k_out * scales[1] * scales[2]),
        scales,
    )
end

PixelShuffle(scale::Int, k_in, k_out = k_in) = PixelShuffle((scale, scale), k_in, k_out)

function (ps::PixelShuffle)(x)
    pixelshuffle(ps.conv(x), ps.scales)
end


text:
"""
Based on
https://hyperopt.github.io/hyperopt/
"""
module TestLogUniform

using StatsBase
using Test
using TreeParzen

# we do this because current version of StatsBase getting pulled in
# doesn't support edges keyword argument
# NOTE: IF YOU DONT HAVE -INF IN THE FIRST VALUE OF EDGES
# YOU MIGHT MISS SOME POINTS -- USE AT OWN RISK
# Each invocation is relatively memory intensive,
# don't use with millions of points (or use at own risk)
# When statsbase is available:
#   fit(StatsBase.Histogram, samples; edges=bin_edges).weights
function edgeshist(values, edges::Vector{})
       x = values .< edges'
       y = sum(x; dims=1)
       return diff(dropdims(y; dims=1))
end

logspace_edges(emin, emax, nbins) = exp.(emin : ((emax - emin) / nbins) : emax)

# test that distributions come out right
@testset "loguniform" begin

    N = 10_000
    COUNTMAX = 1200
    COUNTMIN = 800
    N_BINS = 10

    samples = [TreeParzen.Resolve.node(HP.LogUniform(:lu, -2.0, 2.0), TreeParzen.Trials.ValsDict()) for i in 1:N]
    bin_edges  = logspace_edges(-2.0, 2.0, N_BINS)

    @test length(samples) == N
    @test exp(-2) < minimum(samples)
    @test maximum(samples) < exp(2)
    h = edgeshist(samples, bin_edges)
    @test all(COUNTMIN .< h)
    @test all(h .< COUNTMAX)
end

end
true


text:
function make_simulation(model, dir, filename; Δt=0.002,
                                               stop_time=40.0,
                                               iteration_interval=1000,
                                               schedule_interval=0.5,
                                               output_fields = :b)

    progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time: $(round(Int, sim.model.clock.time))"

    simulation = Simulation(model, Δt=Δt,
                                   stop_time=stop_time,
                                   iteration_interval=iteration_interval,
                                   progress=progress)

    scheduler = TimeInterval(schedule_interval)

    #NOTE: Nan-checker currently not working
    #simulation.diagnostics[:nc] = nc

    output_writer!(simulation, output_fields, model, scheduler, dir, filename)

    return simulation
end
function get_fields(model)
    u, v, w = model.velocities
    b = model.tracers.b

    s = sqrt(v^2 + w^2)
    ω = ∂y(w) - ∂z(v)

    s_field = ComputedField(s)
    ω_field = ComputedField(ω)

    return u, v, w, b, s_field, ω_field
end
function output_writer!(simulation, output_fields, model, scheduler, dir, filename)
    u, v, w, b, s_field, ω_field = get_fields(model)
    if output_fields == :b
        simulation.output_writers[:fields] =
            JLD2OutputWriter(model, b,
                             dir = dir,
                             schedule = scheduler,
                             prefix = filename,
                             force = true)
    elseif output_fields == :all
        simulation.output_writers[:fields] =
            JLD2OutputWriter(model, (b=b, v=v, w=w, ω=ω_field, s=s_field),
                             dir = dir,
                             schedule = scheduler,
                             prefix = filename,
                             force = true)
    else
        throw(ArgumentError("Unknown output field(s)"))
    end
    return nothing
end


text:
module MMKP

using JuMP

export load_folder

struct Item
    profit::Int64
    resources::Vector{Int64}
end
function Item(v::Vector{Int64})
    Item(first(v), v[2:end])
end

const Group = Vector{Item}

struct MMKP_Data
    resource_bounds::Vector{Int}
    groups::Vector{Group}
end

mutable struct MMKP_Prob
    id::Dict
    model::Union{Model,Nothing}
end

include("loaders.jl")
include("moi_model.jl")
include("load_mps.jl")
problems = load_problems()

function get_id(path)
    class, instance = splitpath(path)[end-1:end]
    Dict("class"=>class, "instance"=>first(split(instance, ".")))
end

function load_file(path, loader; n_threads=4)
    id = get_id(path)
    problem = loader(path)
    model = jump_model(problem, n_threads)
    MMKP_Prob(id, model)
end

_get_loader(path) = first(filter(i-> !isnothing(i),
        map(fl->occursin(fl[1], path) ? fl[2] : nothing, folder_keys)))
function load_subfolder(path)
    l = _get_loader(path)
    map(f->load_file(f, l), readdir(path, join=true))
end

function load_folder(path=joinpath(@__DIR__, "benchmark_problems"))
    vcat(load_subfolder.(readdir(path, join=true))...)
end

end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
# using GeometryPrimitives: orthoaxes
using FFTW
using IterativeSolvers: LOBPCGIterator
export MaxwellGrid, MaxwellData, t2c, c2t, kx_c2t, kx_t2c, zx_t2c, kxinv_c2t, kxinv_t2c, ε⁻¹_dot, ε_dot_approx, M, M̂, P, P̂, Mₖ, M̂ₖ, t2c!, c2t!, kcross_c2t!, kcross_t2c!, zcross_t2c!, kcrossinv_c2t!, kcrossinv_t2c!, ε⁻¹_dot!, ε_dot_approx!, M!, M̂!, P!, P̂!, calc_kpg, H_Mₖ_H

struct MaxwellGrid{T<:Real}
    Δx::T
    Δy::T
    Δz::T
    Nx::Int
    Ny::Int
    Nz::Int
    δx::T
    δy::T
    δz::T
    x::StepRangeLen{T}
    y::StepRangeLen{T}
    z::StepRangeLen{T}
    g⃗::Array{SVector{3,T},3}
end

MaxwellGrid(Δx::T,Δy::T,Δz::T,Nx::Int,Ny::Int,Nz::Int) where {T<:Real} = MaxwellGrid{T}(
    Δx,
    Δy,
    Δz,
    Nx,
    Ny,
    Nz,
    Δx / Nx,    # δx
    Δy / Ny,    # δy
    Δz / Nz,    # δz
    ( ( Δx / Nx ) .* (0:(Nx-1))) .- Δx/2.,  # x
    ( ( Δy / Ny ) .* (0:(Ny-1))) .- Δy/2.,  # y
    ( ( Δz / Nz ) .* (0:(Nz-1))) .- Δz/2.,  # z
    # [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)], # g⃗
	[ SVector{T,3}(gx,gy,gz) for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)], # g⃗
)

MaxwellGrid(Δx::T,Δy::T,Nx::Int,Ny::Int) where {T<:Real} = MaxwellGrid{T}(
    Δx,
    Δy,
    1.,
    Nx,
    Ny,
    1,
    Δx / Nx,    # δx
    Δy / Ny,    # δy
    1.,    # δz
    ( ( Δx / Nx ) .* (0:(Nx-1))) .- Δx/2.,  # x
    ( ( Δy / Ny ) .* (0:(Ny-1))) .- Δy/2.,  # y
    0.0:1.0:0.0,  # z
    [ SVector{T,3}(gx,gy,gz) for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(1,1.0)], # g⃗
)

mutable struct MaxwellData{T<:Real}
	grid::MaxwellGrid			# static grid data
	Neigs::Int					# number of eigenvalues/vectors found
	N::Int						# number of real/reciprocal space grid points (N = Nx * Ny * Nz)
	k::T
    ω²::T
    ω²ₖ::T
    ω::T
    ωₖ::T
    H⃗::Array{Complex{T},2}
    H::Array{Complex{T},4}
    e::Array{Complex{T},4}
    d::Array{Complex{T},4}
    mn::Array{T,5}
	kpg_mag::Array{T,3}
    𝓕::FFTW.cFFTWPlan
	𝓕⁻¹::AbstractFFTs.ScaledPlan
    𝓕!::FFTW.cFFTWPlan
	𝓕⁻¹!::AbstractFFTs.ScaledPlan
	solver::LOBPCGIterator
end

MaxwellData(k::T,g::MaxwellGrid{T},Neigs::Int) where {T<:Real} = MaxwellData(
	g,
	Neigs,
	N=g.Nx*g.Ny*g.Nz,
	k,
    0.0,
    0.0,
    0.0,
    0.0,
    H⃗ = randn(Complex{T},(2*g.Nx*g.Ny*g.Nz,Neigs)),
    H = randn(Complex{T},(2,g.Nx,g.Ny,g.Nz)),
    e = randn(Complex{T},(3,g.Nx,g.Ny,g.Nz)),
    d = randn(Complex{T},(3,g.Nx,g.Ny,g.Nz)),
    ( (kpg_mag, mn) = calc_kpg(k,g.Δx,g.Δy,g.Δz,g.Nx,g.Ny,g.Nz); mn),  # mn
	kpg_mag,
	plan_fft(randn(Complex{T}, (3,g.Nx,g.Ny,g.Nz))),  # planned FFT operator 𝓕
	plan_ifft(randn(Complex{T}, (3,g.Nx,g.Ny,g.Nz))),
	𝓕! = plan_fft!(randn(Complex{T}, (3,g.Nx,g.Ny,g.Nz))),
	𝓕⁻¹! = plan_ifft!(randn(Complex{T}, (3,g.Nx,g.Ny,g.Nz))),  # planned in-place FFT operator 𝓕!
	LOBPCGIterator(M̂!(ε⁻¹,H,e,d,kpg_mag,mn,𝓕!,𝓕⁻¹!), 							# Helmholtz Operator
				   false,									   					# "largest", true: find largest eigenvals, false: find smallest
				   H⃗,															 # pre-allocated eigenvectors
				   Neigs,														# number of eigenval/vec pairs to find
				   P=P̂!(ε⁻¹,H,e,d,kpg_mag,mn,𝓕!,𝓕⁻¹!),						   # (Right-)Preconditioner, approximate inverse Helmoltz Operator
				   nothing),													# Constraints
)

MaxwellData(k,g::MaxwellGrid) = MaxwellData(k,g,1)
MaxwellData(k::T,Δx::T,Δy::T,Δz::T,Nx::Int,Ny::Int,Nz::Int) where {T<:Real} = MaxwellData{T}(k,MaxwellGrid{T}(Δx,Δy,Δz,Nx,Ny,Nz))
MaxwellData(k::T,Δx::T,Δy::T,Nx::Int,Ny::Int) where {T<:Real} = MaxwellData{T}(k,MaxwellGrid{T}(Δx,Δy,Nx,Ny))

# non-Mutating Operators

function calc_kpg(kz::T,g⃗::Array{Array{T,1},3})::Tuple{Array{T,3},Array{T,5}} where {T<:Real}
	g⃗ₜ_zero_mask = Zygote.@ignore [ sum(abs2.(gg[1:2])) for gg in g⃗ ] .> 0.
	g⃗ₜ_zero_mask! = Zygote.@ignore .!(g⃗ₜ_zero_mask)
	ŷ = [0.; 1. ;0.]
	k⃗ = [0.;0.;kz]
	@tullio kpg[a,i,j,k] := k⃗[a] - g⃗[i,j,k][a] nograd=g⃗ fastmath=false
	@tullio kpg_mag[i,j,k] := sqrt <| kpg[a,i,j,k]^2 fastmath=false
	zxinds = [2; 1; 3]
	zxscales = [-1; 1. ;0.] #[[0. -1. 0.]; [-1. 0. 0.]; [0. 0. 0.]]
	@tullio kpg_nt[a,i,j,k] := zxscales[a] * kpg[zxinds[a],i,j,k] * g⃗ₜ_zero_mask[i,j,k] + ŷ[a] * g⃗ₜ_zero_mask![i,j,k]  nograd=(zxscales,zxinds,ŷ,g⃗ₜ_zero_mask,g⃗ₜ_zero_mask!) fastmath=false
	@tullio kpg_nmag[i,j,k] := sqrt <| kpg_nt[a,i,j,k]^2 fastmath=false
	@tullio kpg_n[a,i,j,k] := kpg_nt[a,i,j,k] / kpg_nmag[i,j,k] fastmath=false
	xinds1 = [2; 3; 1]
	xinds2 = [3; 1; 2]
	@tullio kpg_mt[a,i,j,k] := kpg_n[xinds1[a],i,j,k] * kpg[xinds2[a],i,j,k] - kpg[xinds1[a],i,j,k] * kpg_n[xinds2[a],i,j,k] nograd=(xinds1,xinds2) fastmath=false
	@tullio kpg_mmag[i,j,k] := sqrt <| kpg_mt[a,i,j,k]^2 fastmath=false
	@tullio kpg_m[a,i,j,k] := kpg_mt[a,i,j,k] / kpg_mmag[i,j,k] fastmath=false
	kpg_mn_basis = [[1. 0.] ; [0. 1.]]
	@tullio kpg_mn[a,b,i,j,k] := kpg_mn_basis[b,1] * kpg_m[a,i,j,k] + kpg_mn_basis[b,2] * kpg_n[a,i,j,k] nograd=kpg_mn_basis fastmath=false
	return kpg_mag, kpg_mn
end

function calc_kpg(kz::T,Δx::T,Δy::T,Δz::T,Nx::Int,Ny::Int,Nz::Int)::Tuple{Array{T,3},Array{T,5}} where {T<:Real}
	g⃗ = Zygote.@ignore [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)]
	calc_kpg(kz,g⃗)
end
# function calc_kpg(kz::T,Δx::T,Δy::T,Δz::T,Nx::Int,Ny::Int,Nz::Int)::Tuple{Array{T,3},Array{T,5}} where T <: Real
# 	g⃗ = Zygote.@ignore [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)]
# 	calc_kpg(kz,g⃗)
# end

# Non-Mutating Operators

"""
    t2c: v⃗ (transverse vector) → a⃗ (cartesian vector)
"""
function t2c(H,mn)
    @tullio h[a,i,j,k] := H[b,i,j,k] * mn[a,b,i,j,k] fastmath=false
end

"""
    c2t: a⃗ (cartesian vector) → v⃗ (transverse vector)
"""
function c2t(h,mn)
    @tullio H[a,i,j,k] := h[b,i,j,k] * mn[b,a,i,j,k] fastmath=false
end

"""
    kx_t2c: a⃗ (cartesian vector) = k⃗ × v⃗ (transverse vector)
"""
function kx_t2c(H,mn,kpg_mag)
	kxscales = [-1.; 1.]
	kxinds = [2; 1]
    @tullio d[a,i,j,k] := kxscales[b] * H[kxinds[b],i,j,k] * mn[a,b,i,j,k] * kpg_mag[i,j,k] nograd=(kxscales,kxinds) fastmath=false
end

"""
    kx_c2t: v⃗ (transverse vector) = k⃗ × a⃗ (cartesian vector)
"""
function kx_c2t(e⃗,mn,kpg_mag)
	kxscales = [-1.; 1.]
    kxinds = [2; 1]
    @tullio H[b,i,j,k] := kxscales[b] * e⃗[a,i,j,k] * mn[a,kxinds[b],i,j,k] * kpg_mag[i,j,k] nograd=(kxinds,kxscales) fastmath=false
end

"""
    kxinv_t2c: compute a⃗ (cartestion vector) st. v⃗ (cartesian vector from two trans. vector components) ≈ k⃗ × a⃗
    This neglects the component of a⃗ parallel to k⃗ (not available by inverting this cross product)
"""
function kxinv_t2c(H,mn,kpg_mag)
	kxinvscales = [1.; -1.]
	kxinds = [2; 1]
    @tullio e⃗[a,i,j,k] := kxscales[b] * H[kxinds[b],i,j,k] * mn[a,b,i,j,k] / kpg_mag[i,j,k] nograd=(kxscales,kxinds) fastmath=false
end

"""
    kxinv_c2t: compute  v⃗ (transverse 2-vector) st. a⃗ (cartestion 3-vector) = k⃗ × v⃗
    This cross product inversion is exact because v⃗ is transverse (perp.) to k⃗
"""
function kxinv_c2t(d⃗,mn,kpg_mag)
	kxscales = [1.; -1.]
    kxinds = [2; 1]
    @tullio H[b,i,j,k] := kxscales[b] * d⃗[a,i,j,k] * mn[a,kxinds[b],i,j,k] / kpg_mag[i,j,k] nograd=(kxinds,kxscales) fastmath=false
end

"""
    zx_t2c: a⃗ (cartesian vector) = ẑ × v⃗ (transverse vector)
"""
function zx_t2c(H,mn)
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxH[a,i,j,k] := zxscales[a] * H[b,i,j,k] * mn[zxinds[a],b,i,j,k] nograd=(zxscales,zxinds) fastmath=false
end

"""
    ε⁻¹_dot_t: e⃗  = ε⁻¹ ⋅ d⃗ (transverse vectors)
"""
function ε⁻¹_dot_t(d⃗,ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * fft(d⃗,(2:4))[b,i,j,k] fastmath=false
	return ifft(e⃗,(2:4))
end

"""
    ε⁻¹_dot: e⃗  = ε⁻¹ ⋅ d⃗ (cartesian vectors)
"""
function ε⁻¹_dot(d⃗,ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * d⃗[b,i,j,k] fastmath=false
end

"""
    ε_dot_approx: approximate     d⃗  = ε ⋅ e⃗
                    using         d⃗  ≈  e⃗ * ( 3 / Tr(ε⁻¹) )
    (all cartesian vectors)
"""
function ε_dot_approx(e⃗,ε⁻¹)
    @tullio d⃗[b,i,j,k] := e⃗[b,i,j,k] * 3 / ε⁻¹[a,a,i,j,k] fastmath=false
end

function M(H,ε⁻¹,mn,kpg_mag)
    kx_c2t(ε⁻¹_dot_t(kx_t2c(H,mn,kpg_mag),ε⁻¹),mn,kpg_mag)
end

function M(H,ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)
    kx_c2t( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * kx_t2c(H,mn,kpg_mag), ε⁻¹), mn,kpg_mag)
end

function M(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag)::Array{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(M(HinA,ε⁻¹,mn,kpg_mag))
end

function M(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Array{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(M(HinA,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹))
end

M̂(ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> M(H,ε⁻¹,mn,kpg_mag,𝓕,𝓕⁻¹)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

function P(Hin::AbstractArray{ComplexF64,4},ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::AbstractArray{ComplexF64,4}
    kxinv_c2t( 𝓕 * ε_dot_approx( 𝓕⁻¹ * kxinv_t2c(H,mn,kpg_mag),ε⁻¹),mn,kpg_mag)
end

function P(Hin::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag,𝓕)::AbstractArray{ComplexF64,1}
    HinA = reshape(Hin,(2,size(ε⁻¹)[end-2:end]...))
    return vec(P(HinA,ε⁻¹,mn,kpg_mag,𝓕))
end

P̂(ε⁻¹,mn,kpg_mag,𝓕) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> P(H,ε⁻¹,mn,kpg_mag,𝓕)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

function Mₖ(H,ε⁻¹,mn,kpg_mag)
    kx_c2t(ε⁻¹_dot_t(zx_t2c(H,mn),ε⁻¹),mn,kpg_mag)
end

function Mₖ(H,ε⁻¹,mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)
    kx_c2t( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * zx_t2c(H,mn), ε⁻¹), mn,kpg_mag)
end

function Mₖ(H::AbstractArray{ComplexF64,1},ε⁻¹,mn,kpg_mag)::AbstractArray{ComplexF64,1}
    Ha = reshape(H,(2,size(ε⁻¹)[end-2:end]...))
    return vec(Mₖ(Ha,ε⁻¹,mn,kpg_mag))
end

function Mₖ(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{ComplexF64,5},mn,kpg_mag,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Array{ComplexF64,1}
    Ha = reshape(H,(2,size(ε⁻¹)[end-2:end]...))
    return vec(M(Ha,ε⁻¹,mn,kpg_mag,𝓕))
end

M̂ₖ(ε⁻¹,mn,kpg_mag,𝓕) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> Mₖ(H,ε⁻¹,mn,kpg_mag,𝓕)::AbstractArray{ComplexF64,1},*(2,size(ε⁻¹)[end-2:end]...),ishermitian=true,ismutating=false)

function H_Mₖ_H(H::AbstractArray{ComplexF64,4},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn)
	kxinds = [2; 1]
	kxscales = [-1.; 1.]
	@tullio out[_] := conj.(H)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_dot_t(zx_t2c(H,mn),ε⁻¹)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) fastmath=false
	return abs(out[1])
end

# function H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{ComplexF64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)::Float64
# 	kxinds = [2; 1]
# 	kxscales = [-1.; 1.]
# 	@tullio out[_] := conj.(H)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ( 𝓕⁻¹ * ε⁻¹_dot( 𝓕 * zx_t2c(H,mn), ε⁻¹) )[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) fastmath=false
# 	return abs(out[1])
# end

# H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{ComplexF64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)  = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn,𝓕,𝓕⁻¹)
H_Mₖ_H(H::AbstractArray{ComplexF64,1},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn) = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn)
# H_Mₖ_H(H::AbstractArray{ComplexF64,2},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn,𝓕::FFTW.cFFTWPlan,𝓕⁻¹)  = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn,𝓕,𝓕⁻¹)
H_Mₖ_H(H::AbstractArray{ComplexF64,2},ε⁻¹::AbstractArray{Float64,5},kpg_mag,mn) = H_Mₖ_H(reshape(H,(2,size(ε⁻¹)[end-2:end]...)),ε⁻¹,kpg_mag,mn)

###### Mutating Operators #######

function kx_tc!(d,H,m,n,mag)::Array{ComplexF64,4}
    # @assert size(Y) === size(X)
    # @assert size(d,4) == 3
    # @assert size(H,4) === 2
    @avx for k ∈ axes(d,3), j ∈ axes(d,2), i ∈ axes(d,1), l in 0:0
        d[i,j,k,1+l] = ( H[i,j,k,1] * n[i,j,k,1+l] - H[i,j,k,2] * m[i,j,k,1+l] ) * -mag[i,j,k]
        d[i,j,k,2+l] = ( H[i,j,k,1] * n[i,j,k,2+l] - H[i,j,k,2] * m[i,j,k,2+l] ) * -mag[i,j,k]
        d[i,j,k,3+l] = ( H[i,j,k,1] * n[i,j,k,3+l] - H[i,j,k,2] * m[i,j,k,3+l] ) * -mag[i,j,k]
    end
    return d
end

function kx_ct!(H,e,m,n,mag,Ninv)::Array{ComplexF64,4}
    # @assert size(Y) === size(X)
    # @assert size(e,4) == 3
    # @assert size(H,4) === 2
    @avx for k ∈ axes(H,3), j ∈ axes(H,2), i ∈ axes(H,1), l in 0:0
        scale = mag[i,j,k] * Ninv
        H[i,j,k,1+l] =  (	e[i,j,k,1+l] * n[i,j,k,1+l] + e[i,j,k,2+l] * n[i,j,k,2+l] + e[i,j,k,3+l] * n[i,j,k,3+l]	) * -scale  # -mag[i,j,k] * Ninv
		H[i,j,k,2+l] =  (	e[i,j,k,1+l] * m[i,j,k,1+l] + e[i,j,k,2+l] * m[i,j,k,2+l] + e[i,j,k,3+l] * m[i,j,k,3+l]	) * scale   # mag[i,j,k] * Ninv
    end
    return H
end

function eid!(e,ei,d)::Array{ComplexF64,4}
    # @assert size(e,4) === 3
    # @assert size(d,4) === 3
    # @assert size(ei,4) === 3
    # @assert size(ei,5) === 3
    @avx for k ∈ axes(e,3), j ∈ axes(e,2), i ∈ axes(e,1), l in 0:0, h in 0:0
        e[i,j,k,1+h] =  ei[i,j,k,1+h,1+l]*d[i,j,k,1+l] + ei[i,j,k,2+h,1+l]*d[i,j,k,2+l] + ei[i,j,k,3+h,1+l]*d[i,j,k,3+l]
        e[i,j,k,2+h] =  ei[i,j,k,1+h,2+l]*d[i,j,k,1+l] + ei[i,j,k,2+h,2+l]*d[i,j,k,2+l] + ei[i,j,k,3+h,2+l]*d[i,j,k,3+l]
        e[i,j,k,3+h] =  ei[i,j,k,1+h,3+l]*d[i,j,k,1+l] + ei[i,j,k,2+h,3+l]*d[i,j,k,2+l] + ei[i,j,k,3+h,3+l]*d[i,j,k,3+l]
    end
    return e
end

# function M!(Hout::AbstractArray{T,4},Hin::AbstractArray{T,4},e,d,ε⁻¹,m,n,mag,𝓕!,𝓕⁻¹!,Ninv) where {T<:Real}
function M!(Hout,Hin,e,d,ε⁻¹,m,n,mag,𝓕!,𝓕⁻¹!,Ninv)
    kx_tc!(d,Hin,m,n,mag);
    mul!(d,𝓕!,d);
    eid!(e,ε⁻¹,d);
    mul!(d,𝓕⁻¹!,d);
    kx_ct!(Hout,e,m,n,mag,Ninv);
end

#
# function t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @fastmath @inbounds ds.e[1,i,j,k] = ( Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[1,2,i,j,k] )
#         @fastmath @inbounds ds.e[2,i,j,k] = ( Hin[1,i,j,k] * ds.mn[2,1,i,j,k] + Hin[2,i,j,k] * ds.mn[2,2,i,j,k] )
#     	@fastmath @inbounds ds.e[3,i,j,k] = ( Hin[1,i,j,k] * ds.mn[3,1,i,j,k] + Hin[2,i,j,k] * ds.mn[3,2,i,j,k] )
# 	end
#     return ds.e
# end
#
# function c2t!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @fastmath @inbounds ds.e[1,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[2,1,i,j,k] + Hin[3,i,j,k] * ds.mn[3,1,i,j,k]
#         @fastmath @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,2,i,j,k] + Hin[2,i,j,k] * ds.mn[2,2,i,j,k] + Hin[3,i,j,k] * ds.mn[3,2,i,j,k]
#     end
#     return ds.e
# end
#
# function zcross_t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @fastmath @inbounds ds.e[1,i,j,k] = -Hin[1,i,j,k] * ds.mn[2,1,i,j,k] - Hin[2,i,j,k] * ds.mn[2,2,i,j,k]
#         @fastmath @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.mn[1,1,i,j,k] + Hin[2,i,j,k] * ds.mn[1,2,i,j,k]
#     end
#     return ds.e
# end
#
# function kcross_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @fastmath @inbounds ds.d[1,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[1,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[1,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
#         @fastmath @inbounds ds.d[2,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[2,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[2,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
#         @fastmath @inbounds ds.d[3,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[3,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[3,1,i,j,k] ) * -ds.kpg_mag[i,j,k]
#     end
#     return ds.d
# end
#
# function kcross_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# 		@fastmath @inbounds  ds.H[1,i,j,k] =  (	ds.e[1,i,j,k] * ds.mn[1,2,i,j,k] + ds.e[2,i,j,k] * ds.mn[2,2,i,j,k] + ds.e[3,i,j,k] * ds.mn[3,2,i,j,k]	) * -ds.kpg_mag[i,j,k]
# 		@fastmath @inbounds  ds.H[2,i,j,k] =  (	ds.e[1,i,j,k] * ds.mn[1,1,i,j,k] + ds.e[2,i,j,k] * ds.mn[2,1,i,j,k] + ds.e[3,i,j,k] * ds.mn[3,1,i,j,k]	) * ds.kpg_mag[i,j,k]
#     end
#     return ds.H
# end
#
# function ε⁻¹_dot!(ε⁻¹::Array{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @fastmath @inbounds ds.e[1,i,j,k] =  ε⁻¹[1,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,1,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,1,i,j,k]*ds.d[3,i,j,k]
#         @fastmath @inbounds ds.e[2,i,j,k] =  ε⁻¹[1,2,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,2,i,j,k]*ds.d[3,i,j,k]
#         @fastmath @inbounds ds.e[3,i,j,k] =  ε⁻¹[1,3,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,3,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,3,i,j,k]*ds.d[3,i,j,k]
#         # ds.e[1,i,j,k] =  ε⁻¹[1,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[1,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[1,3,i,j,k]*ds.d[3,i,j,k]
#         # ds.e[2,i,j,k] =  ε⁻¹[2,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[2,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[2,3,i,j,k]*ds.d[3,i,j,k]
#         # ds.e[3,i,j,k] =  ε⁻¹[3,1,i,j,k]*ds.d[1,i,j,k] + ε⁻¹[3,2,i,j,k]*ds.d[2,i,j,k] + ε⁻¹[3,3,i,j,k]*ds.d[3,i,j,k]
#     end
#     return ds.e
# end
#
# function M!(ε⁻¹::Array{Float64,5},ds::MaxwellData)::Array{ComplexF64,4}
#     kcross_t2c!(ds);
#     mul!(ds.d,ds.𝓕!,ds.d);
#     ε⁻¹_dot!(ε⁻¹,ds);
# 	mul!(ds.e,ds.𝓕⁻¹!::AbstractFFTs.ScaledPlan,ds.e);
#     kcross_c2t!(ds)
# end

function M!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{Float64,5},ds::MaxwellData)::Array{ComplexF64,1}
    @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
    M!(ε⁻¹,ds);
    @inbounds Hout .= vec(ds.H)
end

M̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData) = LinearMap{ComplexF64}((2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true) do y::AbstractVector{ComplexF64},x::AbstractVector{ComplexF64}
    M!(y,x,ε⁻¹,ds)::AbstractArray{ComplexF64,1}
    end

# function M̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData)
#     function f!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         M!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(f!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end

### Preconditioner P̂ & Component Operators (approximate inverse operations of M̂) ###

function kcrossinv_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds scale::Float64 = inv(ds.kpg_mag[i,j,k])
        @fastmath @inbounds ds.e[1,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[1,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[1,1,i,j,k] ) * scale
        @fastmath @inbounds ds.e[2,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[2,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[2,1,i,j,k] ) * scale
        @fastmath @inbounds ds.e[3,i,j,k] = ( ds.H[1,i,j,k] * ds.mn[3,2,i,j,k] - ds.H[2,i,j,k] * ds.mn[3,1,i,j,k] ) * scale
    end
    return ds.e
end

function kcrossinv_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds scale = -inv(ds.kpg_mag[i,j,k])
        @fastmath @inbounds ds.H[1,i,j,k] =  (	ds.d[1,i,j,k] * ds.mn[1,2,i,j,k] + ds.d[2,i,j,k] * ds.mn[2,2,i,j,k] + ds.d[3,i,j,k] * ds.mn[3,2,i,j,k]	) * -scale
        @fastmath @inbounds ds.H[2,i,j,k] =  (	ds.d[1,i,j,k] * ds.mn[1,1,i,j,k] + ds.d[2,i,j,k] * ds.mn[2,1,i,j,k] + ds.d[3,i,j,k] * ds.mn[3,1,i,j,k]	) * scale
    end
    return ds.H
end

function ε_dot_approx!(ε⁻¹::AbstractArray{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    @fastmath @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
        @fastmath @inbounds ε_ave = 3. * inv( ε⁻¹[1,1,i,j,k] + ε⁻¹[2,2,i,j,k] + ε⁻¹[3,3,i,j,k] ) # tr(ε⁻¹[:,:,i,j,k])
        @fastmath @inbounds ds.d[1,i,j,k] =  ε_ave * ds.e[1,i,j,k]
        @fastmath @inbounds ds.d[2,i,j,k] =  ε_ave * ds.e[2,i,j,k]
        @fastmath @inbounds ds.d[3,i,j,k] =  ε_ave * ds.e[3,i,j,k]
    end
    return ds.d
end

function P!(ε⁻¹::AbstractArray{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,4}
    kcrossinv_t2c!(ds);
    # ds.𝓕⁻¹! * ds.e;
    # ldiv!(ds.e,ds.𝓕!,ds.e)
	mul!(ds.e,ds.𝓕⁻¹!,ds.e);
    ε_dot_approx!(ε⁻¹,ds);
    # ds.𝓕! * ds.d;
    mul!(ds.d,ds.𝓕!,ds.d);
    kcrossinv_c2t!(ds)
end

function P!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{Float64,5},ds::MaxwellData)::AbstractArray{ComplexF64,1}
    @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
    P!(ε⁻¹,ds);
    @inbounds Hout .= vec(ds.H)
end

P̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData) = LinearMap{ComplexF64}((2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true) do y::AbstractVector{ComplexF64},x::AbstractVector{ComplexF64}
	P!(y,x,ε⁻¹,ds)::AbstractArray{ComplexF64,1}
    end

# function P̂!(ε⁻¹::Array{Float64,5},ds::MaxwellData)
#     function fp!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         P!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(fp!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end

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

# function zcross_t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds ds.e[1,i,j,k] = -Hin[1,i,j,k] * ds.kpG[i,j,k].m[2] - Hin[2,i,j,k] * ds.kpG[i,j,k].n[2]
#         @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[1]
#     end
#     return ds.e
# end
#
# function kcross_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = -ds.kpg_mag #-ds.kpG[i,j,k].mag
#         ds.d[1,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         ds.d[2,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         ds.d[3,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return ds.d
# end
#
# function kcross_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = ds.kpG[i,j,k].mag
#         at1 = ds.e[1,i,j,k] * ds.kpG[i,j,k].m[1] + ds.e[2,i,j,k] * ds.kpG[i,j,k].m[2] + ds.e[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         at2 = ds.e[1,i,j,k] * ds.kpG[i,j,k].n[1] + ds.e[2,i,j,k] * ds.kpG[i,j,k].n[2] + ds.e[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         ds.H[1,i,j,k] =  -at2 * scale
#         ds.H[2,i,j,k] =  at1 * scale
#     end
#     return ds.H
# end
#
# function kcrossinv_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = 1 / ds.kpG[i,j,k].mag
#         ds.e[1,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         ds.e[2,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         ds.e[3,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return ds.e
# end
#
# function kcrossinv_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = -1 / ds.kpG[i,j,k].mag
#         at1 = ds.d[1,i,j,k] * ds.kpG[i,j,k].m[1] + ds.d[2,i,j,k] * ds.kpG[i,j,k].m[2] + ds.d[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         at2 = ds.d[1,i,j,k] * ds.kpG[i,j,k].n[1] + ds.d[2,i,j,k] * ds.kpG[i,j,k].n[2] + ds.d[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         ds.H[1,i,j,k] =  -at2 * scale
#         ds.H[2,i,j,k] =  at1 * scale
#     end
#     return ds.H
# end
#
# function ε⁻¹_dot!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         ds.e[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,1]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,1]*ds.d[3,i,j,k]
#         ds.e[2,i,j,k] =  ε⁻¹[i,j,k][1,2]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,2]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,2]*ds.d[3,i,j,k]
#         ds.e[3,i,j,k] =  ε⁻¹[i,j,k][1,3]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,3]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,3]*ds.d[3,i,j,k]
#         # ds.e[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][1,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][1,3]*Hin[3,i,j,k]
#         # ds.e[2,i,j,k] =  ε⁻¹[i,j,k][2,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][2,3]*Hin[3,i,j,k]
#         # ds.e[3,i,j,k] =  ε⁻¹[i,j,k][3,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][3,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,3]*Hin[3,i,j,k]
#     end
#     return ds.e
# end
#
# function ε_dot_approx!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         ε_ave = 3 / tr(ε⁻¹[i,j,k])
#         ds.d[1,i,j,k] =  ε_ave * ds.e[1,i,j,k]
#         ds.d[2,i,j,k] =  ε_ave * ds.e[2,i,j,k]
#         ds.d[3,i,j,k] =  ε_ave * ds.e[3,i,j,k]
#     end
#     return ds.d
# end
#
# function M!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     kcross_t2c!(ds);
#     # ds.𝓕! * ds.d;
#     mul!(ds.d,ds.𝓕!,ds.d);
#     ε⁻¹_dot!(ε⁻¹,ds);
#     # ds.𝓕⁻¹! * ds.e;
#     mul!(ds.e,ds.𝓕⁻¹!,ds.e)
#     kcross_c2t!(ds)
# end
#
# function M!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     # copyto!(ds.H,reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz)))
#     @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     M!(ε⁻¹,ds);
#     # copyto!(Hout,vec(ds.H))
#     @inbounds Hout .= vec(ds.H)
# end
#
# function M̂!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)
#     function f!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         M!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(f!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end
#
# function P!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     kcrossinv_t2c!(ds);
#     # ds.𝓕⁻¹! * ds.e;
#     mul!(ds.e,ds.𝓕⁻¹!,ds.e)
#     ε_dot_approx!(ε⁻¹,ds);
#     # ds.𝓕! * ds.d;
#     mul!(ds.d,ds.𝓕!,ds.d);
#     kcrossinv_c2t!(ds)
# end
#
# function P!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     # copyto!(ds.H,reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz)))
#     @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     P!(ε⁻¹,ds);
#     # copyto!(Hout,vec(ds.H))
#     @inbounds Hout .= vec(ds.H)
# end
#
# function P̂!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)
#     function fp!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         P!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(fp!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end


text:
using Distributions, CmdStan, Random, Parameters
ProjDir = @__DIR__
cd(ProjDir)
include("LNR.jl")

stream = open("LNR_Model.stan","r")
Model = read(stream,String)
close(stream)

stanmodel = Stanmodel(Sample(save_warmup=false, num_warmup=1000,
  num_samples=1000, thin=1), nchains=1, name="LNR", model=Model,
  printsummary=false, output_format=:mcmcchains)

Random.seed!(343)
Nreps = 50
Nr = 3
Nobs = 10

for i in 1:Nreps
    μ = rand(Normal(0 ,3), Nr)
    σ = rand(Uniform(.2, 1))
    dist = LNR(μ=μ, σ=σ, ϕ=0.0)
    data = rand(dist, Nobs)
    choice = map(x->x[1], data)
    rt = map(x->x[2], data)
    stanInput = Dict("choice"=>choice, "rt"=>rt, "Nr"=>Nr, "N"=>Nobs)
    rc, chain, cnames = stan(stanmodel, stanInput, ProjDir)
end

max_stars_repo_name: itsdfish/Miscellaneous

max_issues_repo_name: itsdfish/Miscellaneous

max_forks_repo_name: itsdfish/Miscellaneous


text:
"""
# Bucket Sort

Bucket Sort is a divide and conquer sorting algorithm.
It distributes elements into buckets and sort each bucket before merging all buckets
Buckets are arrays that each stores a different range of elements
Insertion sort is used to sort the buckets because its fast in small-sized arrays

## Bucket Calculation

number of buckets = length of array
size of buckets = (largest element - smallest element + 1) / length of array

## Bucket Sort Steps

1. Create an array of buckets
2. Add each elements into its respective bucket
3. Sort each non-empty bucket
4. Concatenate elements of non-empty buckets

## Example

1. Array of [3, 5, 2, 9] has 4 buckets
2. The size of buckets = (9 - 2 + 1) / 4 = 2
3. The range of elements can be assign to each bucket, x
    bucket 1: 2 <= x < 4
    bucket 2: 4 <= x < 6
    bucket 3: 6 <= x < 8
    bucket 4: 8 <= x <= 9
4. The distribution of elements in the buckets = [3, 2][5][][9]
5. After sorting the elements each bucket = [2, 3][5][][9]
6. Concatenate all elements = [2, 3, 5, 9]

# Example
```julia
BucketSort!([3, 5, 2, 9]) # returns [2, 3, 5, 9]
```

Contributed By:- [Ming Liang](https://github.com/DrakonDarquesse)
"""

function BucketSort!(arr::Vector{T}, l::Int=length(arr)) where T
    if l==0
        return
    end

    max = maximum(arr)
    min = minimum(arr)
    r = (max-min+1)/l
    buckets = Array{T, 1}[[] for a in 1:l]

    for i in arr
        push!(buckets[floor(Int, (i-min)/r)+1], i)
    end

    for bucket in buckets
        for i in 1:length(bucket)
            val = bucket[i]
            idx = i
            while idx > 1 && val < bucket[idx-1]
                bucket[idx] = bucket[idx-1]
                idx = idx - 1
            end

            bucket[idx] = val
        end
    end
    copy!(arr, vcat(buckets...))
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
@testset "NOCB" begin
    @testset "Default" begin
        tester = ImputorTester(NOCB)

        test_hashing(tester)
        test_equality(tester)
        test_vector(tester)
        test_matrix(tester)
        # test_cube(tester)
        test_dataframe(tester)
        test_groupby(tester)
        test_axisarray(tester)
        test_nameddimsarray(tester)
        test_keyedarray(tester)
        test_columntable(tester)
        test_rowtable(tester)

        test_limited(tester)

        @testset "Cube" begin
            a = allowmissing(1.0:1.0:60.0)
            a[[2, 7, 18, 23, 34, 41, 55, 59, 60]] .= missing
            C = collect(reshape(a, 5, 4, 3))

            # Cube tests are expected to fail
            @test_throws MethodError impute(C, tester.imp(; tester.kwargs...); dims=3)
        end
    end

    @testset "Floats" begin
        a = allowmissing(1.0:1.0:20.0)
        a[[2, 3, 7]] .= missing

        result = impute(a, NOCB())
        expected = copy(a)
        expected[2] = 4.0
        expected[3] = 4.0
        expected[7] = 8.0

        @test result == expected
        @test result == Impute.nocb(a)

        a2 = copy(a)
        Impute.nocb!(a2)
        @test a2 == result

        # Test LOCF at endpoints
        b = ones(Union{Float64, Missing}, 20)
        b[[1, 3, 20]] .= missing
        result = Impute.nocb(b)
        @test result[1] == 1.0
        @test ismissing(result[20])

        # Test limiting
        a[11:15] .= missing

        expected = copy(a)
        @test isequal(impute(a, NOCB(; limit=0)), expected)

        expected[2] = 4.0
        expected[3] = 4.0
        expected[7] = 8.0
        expected[13:15] .= 16.0

        @test isequal(impute(a, NOCB(; limit=3)), expected)
    end

    @testset "Ints" begin
        a = allowmissing(1:1:20)
        a[[2, 3, 7]] .= missing

        result = impute(a, NOCB())
        expected = copy(a)
        expected[2] = 4
        expected[3] = 4
        expected[7] = 8

        @test result == expected
    end

    @testset "Strings" begin
        a = allowmissing([randstring(4) for i in 1:20])
        a[[2, 3, 7]] .= missing

        result = impute(a, NOCB())
        expected = copy(a)
        expected[2] = expected[4]
        expected[3] = expected[4]
        expected[7] = expected[8]

        @test result == expected
    end
end


text:
println("""MIT License

Copyright (c) 2021 Wenlu Shi

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

	""")


text:
# SPDX-License-Identifier: MIT

"""
    struct ProbabilitySemifield{T<:AbstractFloat} <: Semfield
        val::T
    end

Log-semifield is defined as :
  * ``x \\oplus y \\triangleq x + y``
  * ``x \\otimes y \\triangleq x \\cdot y``
  * ``x \\oslash y \\triangleq \\frac{x}{y}``
``\\forall x, y \\in [0, 1]``.
"""
struct ProbabilitySemifield{T<:AbstractFloat} <: Semifield
    val::T
end

Base.:+(x::ProbabilitySemifield{T}, y::ProbabilitySemifield{T}) where T =
    ProbabilitySemifield{T}(x.val + y.val)
Base.:*(x::ProbabilitySemifield, y::ProbabilitySemifield) =
    ProbabilitySemifield(x.val * y.val)
Base.:/(x::ProbabilitySemifield, y::ProbabilitySemifield) =
    ProbabilitySemifield(x.val / y.val)
Base.zero(::Type{ProbabilitySemifield{T}}) where T =
    ProbabilitySemifield{T}(T(0))
Base.one(::Type{ProbabilitySemifield{T}}) where T =
    ProbabilitySemifield{T}(T(1))
Base.isless(x::ProbabilitySemifield, y::ProbabilitySemifield) =
    isless(x.val, y.val)
Base.typemin(x::Type{ProbabilitySemifield{T}}) where T = zero(T)
Base.typemax(x::Type{ProbabilitySemifield{T}}) where T = one(T)


text:
matrix = zeros(Float64, (1000,1000,1000))
# use matrix, then when done set variable to 0 to garbage collect the matrix:
matrix = 0 # large memory pool will now be collected when needed

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
push!(LOAD_PATH, "./models")
push!(LOAD_PATH, ".")
push!(LOAD_PATH, "./SparseInverseProblems/src")
using SparseInverseProblemsMod
using SuperResModels
using TestCases
using Utils
using PyCall
@pyimport numpy as np
include("1d_parameters.jl")

# To see the progress in pmap
# To use PmapProgressMeter you need to clone it manually
# do in Julia: Pkg.clone("https://github.com/slundberg/PmapProgressMeter.jl")
# if this is no longer existent, uncoment the alternative pmap in this code.
@everywhere using ProgressMeter
@everywhere using PmapProgressMeter

model_static = SuperResModels.Fourier1d(f_c, x_max, filter_x, num_x)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, v_max, tau, K, num_v*K)
if K*tau != 1.0
    warning(" The function match_points < Utils considers K*tau = 1 to set distance in Omega, change this if it's not the case")
end

@everywhere begin
    model_static = $model_static
    model_dynamic = $model_dynamic
    x_max = $x_max
    v_max = $v_max
    minWeights = $minWeights
    maxWeights = $maxWeights
    K = $K
    tau = $tau
    dataFolder = $dataFolder
    test_case = $test_case
    noises_data = $noises_data
    noises_position = $noises_position
    bins = $bins
    density = $density
    cases = $cases
end

results = pmap(x -> begin sleep(1) ; Utils.generate_and_reconstruct_all(model_static, model_dynamic, bins, density, algIter, test_case, noises_data, noises_position, cases) end, Progress(num_trials), 1:num_trials)
# results = pmap(x -> Utils.generate_and_reconstruct_all(model_static, model_dynamic, bins, density, test_case, noises_data, noises_position, cases), 1:num_trials)
results_array = vcat([result[3] for result in results]...)
separations_array = vcat([result[1] for result in results]...)
separation_dyn_array = vcat([result[2] for result in results]...)

#
now_str = string(now())
now_str = replace(now_str, ":", "-")
now_str = replace(now_str, ".", "-")

# Save results
saveFolder = dataFolder*"/"*now_str
mkdir(saveFolder)
cp("1d_parameters.jl", string(saveFolder, "/1d_parameters.jl"))
cd(saveFolder)

println(" Simulations finished, we save the files into folder :/" *saveFolder )

np.save("datanoise", noises_data)
np.save("positionnoise.npy", noises_position)
np.save("separationDynamic.npy", separation_dyn_array)
np.save("separations.npy", separations_array)
np.save("results.npy", results_array)