text:
@testset "Sample properties for data frames" begin
    df = DataFrames.DataFrame(:x => [1,2,3,4], :y => ["a", "d", "b", "c"])
    @test Banyan.sample_axes(df) == [1]
    @test Banyan.sample_keys(df) == [:x, :y]
    @test Banyan.sample_min(df, :x) == 1
    @test Banyan.sample_max(df, :x) == 4
    @test Banyan.sample_min(df, :y) == orderinghash("a")
    @test Banyan.sample_max(df, :y) == orderinghash("d")
end

@testset "Sample divisions for data frames" begin
    df1 = DataFrames.DataFrame(:x => [1,2,3,4], :y => ["a", "d", "b", "c"])
    df2 = DataFrames.DataFrame(:x => [1,2,3,4,4], :y => ["a", "d", "b", "b", "c"])
    df3 = DataFrames.DataFrame(:x => [1,2,3,4,5], :y => ["a", "d", "b", "b", "c"])

    # Test df1
    @test sample_divisions(df1, :x) == [
        (1, 2)
        (2, 3)
        (3, 4)
        (4, 4)
    ]
    sample_divisions_df1_y = sample_divisions(df1, :y)
    @test length(sample_divisions_df1_y) == 4
    for (i, division) in enumerate(sample_divisions_df1_y)
        # Each division is a tuple of min/max values representing the range
        # [min, max); i.e., inclusive of minimum and exclusive of the maximum.
        # Minimums and maximums are ignored for the first and final partitions.
        # In other words, the first partition is not lower-bounded and the
        # final partition is not upper-bounded.
        if i < length(sample_divisions_df1_y)
            @test division[2] == sample_divisions_df1_y[i+1][1]
        end

        # The very last range should just be a single "d"
        if i < length(sample_divisions_df1_y)
            @test division[1] < division[2]
        else
            @test division[1] == division[2]
        end
    end

    # Test df2 and df3
    for df in [df2, df3]
        sample_divisions_df_y = sample_divisions(df, :y)
        @test length(sample_divisions_df_y) == 2
        for (i, division) in enumerate(sample_divisions_df_y)
            if i < length(sample_divisions_df_y)
                @test division[2] == sample_divisions_df_y[i+1][1]
            end
            @test division[1] < division[2]
        end
    end

    # Test outliers
    @test sample_divisions(DataFrames.DataFrame(:x=>[5,6,7,8]), :x) == [
        (5, 6)
        (6, 7)
        (7, 8)
        (8, 8)
    ]
    @test sample_divisions(DataFrames.DataFrame(:x=>[5,6,7,8,8,8,8]), :x) == [(5,8)]
    @test sample_divisions(DataFrames.DataFrame(:x=>[5,6,7,8,8,8,20]), :x) == [
        (5, 8)
        (8, 20)
    ]
    @test sample_divisions(DataFrames.DataFrame(:x=>[5,6,7,8,8,8,8,8,20]), :x) == [(5, 20)]
end

@testset "Sample # of groups for data frames" begin
    df1 = DataFrames.DataFrame(:x => [1,2,3,4], :y => ["a", "d", "b", "c"])
    df2 = DataFrames.DataFrame(:x => [1,2,3,4,4], :y => ["a", "d", "b", "b", "c"])
    df3 = DataFrames.DataFrame(:x => [1,2,3,4,5], :y => ["a", "d", "b", "b", "c"])
    @test Banyan.sample_max_ngroups(df1, :x) == 4
    @test Banyan.sample_max_ngroups(df1, :y) == 4
    @test Banyan.sample_max_ngroups(df2, :x) == 2
    @test Banyan.sample_max_ngroups(df2, :y) == 2
    @test Banyan.sample_max_ngroups(df3, :x) == 5
    @test Banyan.sample_max_ngroups(df3, :y) == 2
end

@testset "Sample percentile for data frames" begin
    df1 = DataFrames.DataFrame(:x => [1,2,3,4], :y => ["a", "d", "b", "c"])
    df2 = DataFrames.DataFrame(:x => [1,2,3,4,4], :y => ["a", "d", "b", "b", "c"])
    # df2 = DataFrames.DataFrame(:x => [1,2,3,4,5], :y => ["a", "d", "b", "b", "c"])

    @test sample_percentile(df1, :x, 1, 4) == 1
    @test sample_percentile(df1, :x, 2, 4) == 0.75
    @test sample_percentile(df2, :x, 1, 4) == 1
    @test sample_percentile(df2, :x, 2, 4) == 0.8
end

function test_split_divisions_valid(splits, original=nothing)
    @test length(splits) >= 1
    for (j, split) in enumerate(splits)
        @test length(split) >= 1
        for (i, division) in enumerate(split)
            @test length(division) == 2
            if i < length(split)
                @test division[2] == split[i+1][1]
            end
            @test division[1] <= division[2]
        end

        if j < length(splits)
            @test last(split)[2] == first(splits[j+1])[1]
        end
    end

    if !isnothing(original)
        @test first(first(splits))[1] >= first(original)[1]
        @test last(last(splits))[2] <= last(original)[2]
    end
end

@testset "Splitting sampled divisions into $npieces" for npieces in [2, 1, 5]
    # TOOD: Test `get_partition_idx_from_divisions`

    df = CSV.read(use_data("csv", "Disk", true), DataFrames.DataFrame)
    original_divisions = sample_divisions(df, :species)
    test_split_divisions_valid(Banyan.get_divisions(original_divisions, npieces), original_divisions)
    # @test get_partition_idx_from_divisions("" , original_divisions)
end

@testset "Splitting sampled divisions recursively" begin
    df = CSV.read(use_data("csv", "Disk", true), DataFrames.DataFrame)
    original_divisions = sample_divisions(df, :species)
    test_split_divisions_valid(Banyan.get_divisions(Banyan.get_divisions(original_divisions, 2)[1], 5), original_divisions)
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using Multilane
using MCTS
using POMCP

# dir = "/tmp/sim_data_Aug_22_19_22"
dir = "/scratch/users/zsunberg/sim_data_Aug_23_16_02/"

files = []

N = 3
for i in 1:N
    push!(files, joinpath(dir, "results_$(i)_of_$(N).jld"))
end

tic()
results = @time gather_results(files, save_file=string("results_, ", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld"))
toc()

ext: jl


text:
module SysImageHack

using PackageCompiler: compile_incremental

assetpath(name) = joinpath(@__DIR__, "scripts", name)

function compile_patched_sysimage(sysimage; kwargs...)
    tmp_syso, _curr_syso = compile_incremental(
        assetpath("Project.toml"),
        assetpath("patch.jl");
        kwargs...)
    cp(tmp_syso, sysimage, force=true)
    return
end

end  # module


text:
using SixDOF
using Test

@testset "Forward/Reverse Diff" begin
    include("testderiv.jl")
end


text:
export save_path, load_path

# Write `ceil(log2(max_value))` bits representing `value - 1` into
# `buffer`, where `value` is required to be in the range
# `1:max_value`.
function write_bits!(buffer, value, max_value)
    @assert 1 <= value <= max_value
    max_value -= 1
    value -= 1
    while max_value > 0
        push!(buffer, Bool(value & 1))
        value >>= 1
        max_value >>= 1
    end
end

# Opposite of `write_bits!`.
function read_bits!(buffer, max_value)
    @assert 1 <= max_value
    max_value -= 1
    bit = 1
    value = 0
    while max_value > 0
        if popfirst!(buffer)
            value |= bit
        end
        bit <<= 1
        max_value >>= 1
    end
    return value + 1
end

function save_path(filename, graph, path)
    bits = BitVector(undef, 0)
    write_bits!(bits, length(path), nv(graph) + 1)
    write_bits!(bits, path[1], maximum(vertices(graph)))

    visited = falses(nv(graph))
    for i = 1:length(path) - 1
        v1, v2 = path[i], path[i + 1]
        visited[v1] = true
        neighbors = filter(x -> !visited[x], outneighbors(graph, v1))
        i = findfirst(isequal(v2), sort(neighbors))
        @assert i != nothing
        write_bits!(bits, i, length(neighbors))
    end

    buffer = IOBuffer()
    write(buffer, bits)
    write(filename, reinterpret(UInt8, take!(buffer))[1:ceil(Int, length(bits) / 8)])
end

function load_path(filename, graph)
    bytes = read(filename)
    bits = BitVector(undef, 8 * length(bytes))
    while mod(length(bytes), 8) != 0
        push!(bytes, 0x00)
    end
    read!(IOBuffer(bytes), bits)

    n = read_bits!(bits, nv(graph) + 1)
    path = Vector{eltype(graph)}(undef, 0)
    push!(path, read_bits!(bits, maximum(vertices(graph))))

    visited = falses(nv(graph))
    for i = 1:(n - 1)
        v1 = path[end]
        visited[v1] = true
        neighbors = filter(x -> !visited[x], outneighbors(graph, v1))
        v2 = sort(neighbors)[read_bits!(bits, length(neighbors))]
        push!(path, v2)
    end

    return path
end


text:
export FFTOp
import Base.copy

mutable struct FFTOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: Nothing
  ctprod :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  plan
  iplan
  shift::Bool
  unitary::Bool
end

"""
  FFTOp(T::Type, shape::Tuple, shift=true, unitary=true)

returns an operator which performs an FFT on Arrays of type T

# Arguments:
* `T::Type`       - type of the array to transform
* `shape::Tuple`  - size of the array to transform
* (`shift=true`)  - if true, fftshifts are performed
* (`unitary=true`)  - if true, FFT is normalized such that it is unitary
"""
function FFTOp(T::Type, shape::Tuple, shift=true; unitary=true, cuda::Bool=false)
  if cuda
    @assert CUDA.functional()==true "a functional CUDA setup is required when using the option `cuda=true` in FFTOp"
    plan = plan_fft(CuArray{T}(undef,shape))
    iplan = plan_ifft(CuArray{T}(undef,shape))
  else
    plan = plan_fft(zeros(T, shape);flags=FFTW.MEASURE)
    iplan = plan_ifft(zeros(T, shape);flags=FFTW.MEASURE)
  end
  if unitary
    facF = 1.0/sqrt(prod(shape))
    facB = sqrt(prod(shape))
  else
    facF = 1.0
    facB = prod(shape)
  end

  if shift
    return FFTOp{T}(prod(shape), prod(shape), false, false
              , x->vec(fftshift(plan*fftshift(reshape(x,shape))))*T(facF)
              , nothing
              , y->vec(ifftshift(iplan*ifftshift(reshape(y,shape))))*T(facB)
              , 0, 0, 0
              , plan
              , iplan
              , shift
              , unitary)
  else
    return FFTOp{T}(prod(shape), prod(shape), false, false
            , x->vec(plan*(reshape(x,shape)))*T(facF)
            , nothing
            , y->vec(iplan*(reshape(y,shape)))*T(facB)
            , 0, 0, 0
            , plan
            , iplan
            , shift
            , unitary)
  end
end

function Base.copy(S::FFTOp)
  return FFTOp(eltype(S), size(S.plan), S.shift, unitary=S.unitary)
end


text:
using StatsBase
using Test

n = 5000

# 1D integer counts

x = rand(1:5, n)
w = weights(rand(n))

c = counts(x, 5)
@test size(c) == (5,)
c0 = Int[count(v->v == i, x) for i in 1:5]
@test c == c0
@test counts(x .+ 1, 2:6) == c0
@test proportions(x, 1:5) ≈ (c0 ./ n)

c = counts(x)
@test size(c) == (5,)
c0 = Int[count(v->v == i, x) for i in 1:5]
@test c == c0
@test counts(x .+ 1, 2:6) == c0
@test proportions(x) ≈ (c0 ./ n)

c = counts(x, 5, w)
@test size(c) == (5,)
c0 = Float64[sum(w.values[x .== i]) for i in 1:5]
@test c                      ≈ c0
@test counts(x .+ 1, 2:6, w) ≈ c0
@test proportions(x, 1:5, w) ≈ (c0 ./ sum(w))

c = counts(x, w)
@test size(c) == (5,)
c0 = Float64[sum(w.values[x .== i]) for i in 1:5]
@test c                      ≈ c0
@test counts(x .+ 1, 2:6, w) ≈ c0
@test proportions(x, w)      ≈ (c0 ./ sum(w))

# 2D integer counts

x = rand(1:4, n)
y = rand(1:5, n)
w = weights(rand(n))

c = counts(x, y, (4, 5))
@test size(c) == (4, 5)
c0 = Int[count(t->t != 0,  (x .== i) .& (y .== j)) for i in 1:4, j in 1:5]
@test c == c0
@test counts(x .+ 2, y .+ 3, (3:6, 4:8)) == c0
@test proportions(x, y, (1:4, 1:5)) ≈ (c0 ./ n)

c = counts(x, y)
@test size(c) == (4, 5)
c0 = Int[count(t->t != 0, (x .== i) .& (y .== j)) for i in 1:4, j in 1:5]
@test c == c0
@test counts(x .+ 2, y .+ 3, (3:6, 4:8)) == c0
@test proportions(x, y,) ≈ (c0 ./ n)

c = counts(x, y, (4, 5), w)
@test size(c) == (4, 5)
c0 = Float64[sum(w.values[(x .== i) .& (y .== j)]) for i in 1:4, j in 1:5]
@test c                                     ≈ c0
@test counts(x .+ 2, y .+ 3, (3:6, 4:8), w) ≈ c0
@test proportions(x, y, (1:4, 1:5), w)      ≈ (c0 ./ sum(w))

c = counts(x, y, w)
@test size(c) == (4, 5)
c0 = Float64[sum(w.values[(x .== i) .& (y .== j)]) for i in 1:4, j in 1:5]
@test c                                     ≈ c0
@test counts(x .+ 2, y .+ 3, (3:6, 4:8), w) ≈ c0
@test proportions(x, y, w)                  ≈ (c0 ./ sum(w))

# count map

x = ["a", "b", "a", "a", "b", "c"]
w = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5]

cm = countmap(x)
@test cm["a"] == 3
@test cm["b"] == 2
@test cm["c"] == 1
pm = proportionmap(x)
@test pm["a"] ≈ (1/2)
@test pm["b"] ≈ (1/3)
@test pm["c"] ≈ (1/6)

# testing the radixsort branch of countmap
xx = repeat([6, 1, 3, 1], outer=100_000)
cm = countmap(xx)
@test cm == Dict(1 => 200_000, 3 => 100_000, 6 => 100_000)

# testing the radixsort-based addcounts
xx = repeat([6, 1, 3, 1], outer=100_000)
cm = Dict{Int, Int}()
StatsBase.addcounts_radixsort!(cm,xx)
@test cm == Dict(1 => 200_000, 3 => 100_000, 6 => 100_000)
xx2 = repeat([7, 1, 3, 1], outer=100_000)
StatsBase.addcounts_radixsort!(cm,xx2)
@test cm == Dict(1 => 400_000, 3 => 200_000, 6 => 100_000, 7 => 100_000)

# testing the Dict-based addcounts
cm = Dict{Int, Int}()
StatsBase.addcounts_dict!(cm,xx)
@test cm == Dict(1 => 200_000, 3 => 100_000, 6 => 100_000)

cm = countmap(x, weights(w))
@test cm["a"] == 5.5
@test cm["b"] == 4.5
@test cm["c"] == 3.5

@test cm == countmap(x, w)

pm = proportionmap(x, weights(w))
@test pm["a"] ≈ (5.5 / 13.5)
@test pm["b"] ≈ (4.5 / 13.5)
@test pm["c"] ≈ (3.5 / 13.5)

# testing small bits type
bx = [true, false, true, true, false]
@test countmap(bx) == Dict(true => 3, false => 2)

for T in [UInt8, UInt16, Int8, Int16]
    tx = T[typemin(T), 8, typemax(T), 19, 8]
    @test countmap(tx) == Dict(typemin(T) => 1, typemax(T) => 1, 8 => 2, 19 => 1)
end

@testset "views" begin
    X = view([1,1,1,2,2], 1:5)
    @test countmap(X) == countmap(copy(X))
end


text:
####################
## Canvas Objects ##
####################
abstract type CanvasObject end

include("autogenerated_structs.jl")

# General constructors from dicts and kwargs...
(::Type{CO})(dict) where {CO<:CanvasObject} = Internals.json2canvas(CO, dict)
(::Type{CO})(; kwargs...) where {CO<:CanvasObject} = Internals.json2canvas(CO, kwargs)

# Extract id
Internals.id(x::Int) = x
Internals.id(x::CanvasObject) = x.id
# Get the endpoint path + id
Internals.pid(x::Account) = "/accounts/$(Internals.id(x))"
Internals.pid(x::Assignment) = "/assignments/$(Internals.id(x))"
Internals.pid(x::AssignmentOverride) = "/overrides/$(Internals.id(x))"
Internals.pid(x::Course) = "/courses/$(Internals.id(x))"
Internals.pid(x::File) = "/files/$(Internals.id(x))"
Internals.pid(x::Folder) = "/folders/$(Internals.id(x))"
Internals.pid(x::Group) = "/groups/$(Internals.id(x))"
Internals.pid(x::User) = "/users/$(Internals.id(x))"

# namefield(x) = nothing
function Base.show(io::IO, o::CanvasObject)
    if get(io, :compact, false)
        # uri_id = namefield(o)
        if true # uri_id === nothing
            print(io, typeof(o), "(…)")
        else
            print(io, typeof(o), "($(repr(uri_id)))")
        end
    else
        print(io, "$(typeof(o)) (all fields are Union{T,Nothing}):")
        for field in fieldnames(typeof(o))
            val = getfield(o, field)
            if !(val === nothing)
                println(io)
                print(io, "  $field: ")
                if isa(val, Vector)
                    print(io, typeof(val))
                else
                    show(IOContext(io, :compact => true), val)
                end
            end
        end
    end
end

################
# JSON parsing #
################
function Internals.json2canvas(::Type{CO}, data::AbstractDict{T}) where {CO<:CanvasObject,T}
    # Make sure the struct does not miss any fields
    for k in keys(data)
        if !(Symbol(k) ∈ fieldnames(CO))
            @debug "key `$(k)` missing from `$(CO)` struct definition" CO key=k value=data[k]
        end
    end
    # @assert all(k->Symbol(k) ∈ fieldnames(CO), keys(data))
    CO((begin
            v = get(data, T(x), nothing)
            if v === nothing
                nothing
            else
                # https://github.com/JuliaLang/julia/issues/269
                S = Internals.adjust_fieldtype(fieldtype(CO, x), CO, x)
                Internals.convert′(S, v)
            end
        end
        for x in fieldnames(CO))...)
end

# Work around https://github.com/JuliaLang/julia/issues/269
Internals.adjust_fieldtype(::Type{T}, ::Type, ::Symbol) where T = T

Internals.convert′(::Type{T}, val::T) where T = val
Internals.convert′(::Type{T}, val) where T = error("T=$T, val=$val")
function Internals.convert′(::Type{Union{CO,Nothing}}, val) where CO <: CanvasObject
    return CO(val)
end
function Internals.convert′(::Type{Union{Vector{CO},Nothing}}, val) where CO <: CanvasObject
    return CO.(val)
end
function Internals.convert′(::Type{Union{TimeZones.ZonedDateTime,Nothing}}, val)
    df = Dates.dateformat"yyyy-mm-ddTHH:MM:SSz"
    zdt = parse(TimeZones.ZonedDateTime, val, df)
    zdt′ = TimeZones.astimezone(zdt, TimeZones.localzone())
    return convert(Union{TimeZones.ZonedDateTime,Nothing}, zdt′)
end
function Internals.convert′(::Type{Union{Dates.Date,Nothing}}, val)
    df = Dates.dateformat"yyyy-mm-dd"
    d = parse(Dates.Date, val, df)
    return convert(Union{Dates.Date,Nothing}, d)
end


text:
# load

using Flux
using Flux: onehotbatch, argmax, crossentropy, throttle, @epochs
using BSON: @save, @load
using Base.Iterators: partition

m = Chain(
  Conv((2,2), 1=>16, relu),
  x -> maxpool(x, (2,2)),
  Conv((2,2), 16=>8, relu),
  x -> maxpool(x, (2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(288, 10), softmax) |> gpu

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(argmax(m(x)) .== argmax(y))

evalcb = throttle(() -> @show(accuracy(tX, tY)), 1)
opt = ADAM(params(m))

@load "fashion_mnist_model.bson" m
@load "fashion_mnist_weights.bson" weights
Flux.loadparams!(m, weights)

using MLDatasets # FashionMNIST
using ColorTypes: N0f8, Gray

const Img = Matrix{Gray{N0f8}}

function prepare_test()
    # load full test set
    test_x,  test_y  = FashionMNIST.testdata() # 10_000
    test_imgs = Img.([test_x[:,:,i] for i in 1:10_000])

    # Prepare test set
    tX = cat(4, float.(test_imgs[1:10_000])...) |> gpu
    tY = onehotbatch(test_y[1:10_000], 0:9) |> gpu
    tX, tY
end

tX, tY = prepare_test()
@show(accuracy(tX, tY))
#  accuracy(tX, tY) = 0.3329


text:
# A Classical Unconstrained Test Problem
#
# Rosenbrock, H H, An Automatic Method for finding the Greatest or
# least value of a function. Computer Journal 3 (1960), 175-184.

using GAMS
using JuMP
using Test

function example_rosenbrock(; verbose = true)
   model = Model(GAMS.Optimizer)
   set_optimizer_attribute(model, MOI.Silent(), !verbose)

   @variable(model, -10 <= x1 <= 5, start = -1.2)
   @variable(model, -10 <= x2 <= 10, start = 1.0)

   @NLobjective(model, Min, 100 * (x2 - x1^2)^2 + -(1 - x1)^2)

   if verbose
      print(model)
   end

   JuMP.optimize!(model)

   obj_opt = JuMP.objective_value(model)
   x1_opt = JuMP.value(x1)
   x2_opt = JuMP.value(x2)

   if verbose
      println("Objective value: ", obj_opt)
      println("x1 = ", x1_opt)
      println("x2 = ", x2_opt)
   end
end

example_rosenbrock(verbose = false)


text:
export NeuronType, LIF, Izhikevich, NeuronGroup

abstract type NeuronType end

@with_kw struct LIF <: NeuronType
    dt::Real = 1
    v_thresh::Real = -52.0
    v_rest::Real = -65.0
    v_reset::Real = -65.0
    refractory_period::Int = 5
    voltage_time_constant::Real = 100
    voltage_decay_factor::Real = exp(-dt / voltage_time_constant)
    @assert voltage_decay_factor == exp(-dt / voltage_time_constant)
end

@with_kw struct Izhikevich <: NeuronType
    dt::Real = 1
    a::Union{Real,AbstractMatrix} = 0.02
    b::Union{Real,AbstractMatrix} = 0.2
    c::Union{Real,AbstractMatrix} = -65.0
    d::Union{Real,AbstractMatrix} = 8.0
    v_thresh::Union{Real,AbstractMatrix} = 30.0
    runge_kutta_order::Int = 1
end

function Izhikevich(type::String; dt=1, runge_kutta_order=1)
    if type == "RS" || type == "regular spiking"
        Izhikevich(
            dt=dt,
            a=0.02,
            b=0.2,
            c=-65.0,
            d=8.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "IB" || type == "intrinsically bursting"
        Izhikevich(
            dt=dt,
            a=0.02,
            b=0.2,
            c=-55.0,
            d=4.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "CH" || type == "chattering"
        Izhikevich(
            dt=dt,
            a=0.02,
            b=0.2,
            c=-50.0,
            d=2.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "FS" || type == "fast spiking"
        Izhikevich(
            dt=dt,
            a=0.1,
            b=0.2,
            c=-65.0,
            d=2.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "LTS" || type == "low-threshold spiking"
        Izhikevich(
            dt=dt,
            a=0.02,
            b=0.25,
            c=-65.0,
            d=2.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "TC" || type == "thalamo-cortical"
        Izhikevich(
            dt=dt,
            a=0.02,
            b=0.25,
            c=-65.0,
            d=0.05,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    elseif type == "RZ" || type == "resonator"
        Izhikevich(
            dt=dt,
            a=0.1,
            b=0.26,
            c=-65.0,
            d=2.0,
            v_thresh=30.0,
            runge_kutta_order=runge_kutta_order,
        )
    else
        error("Unknown neuron type. Use one of 'RS' or 'regular spiking', 'IB' or 'intrinsically bursting', 'CH' or 'chattering', 'FS' or 'fast spiking', 'LTS' or 'low-threshold spiking', 'TC' or 'thalamo-cortical', 'RZ' or 'resonator'")
    end
end

struct NeuronGroup
    n::Int
    idx::Union{UnitRange{Int64},Vector{Int64}}
end

Base.:+(A::NeuronGroup, B::NeuronGroup) = NeuronGroup(A.n + B.n, sort([A.idx; B.idx]))


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/16/2014

module AbstractSensorImpl

export AbstractSensor

abstract AbstractSensor

end

module AbstractSensorInterfaces

export updateSensor

function updateSensor()
end

end


text:
import DataFrames
import StatsBase

# d = readtable("./data/people.csv", makefactors = true)

# Convert DataFrames.DataFrame or DataTable object to matrix
function dataframe_to_matrix(dat, factor_indcs, n, p)
    X = zeros(n, p)

    for j = 1:p
        if j ∈ factor_indcs
            X[:, j] = factor_to_float(dat[:, j])
        else
            X[:, j] = convert(Array{Float64, 1}, dat[:, j])
        end
    end
    X
end

# Convert matrix to DataFrames.DataFrame or DataTable object
function matrix_to_dataframe(X_new::Array{Float64, 2}, dat::DataFrames.DataFrame, factor_indcs::Array{Int, 1})
    X_synth = DataFrames.DataFrame()
    p = size(X_new, 2)
    for j = 1:p
        if j ∈ factor_indcs
            X_synth[:, Symbol(:x, j)] = float_to_factor(X_new[:, j],
                                            DataFrames.levels(dat[:, j]))
        else
            X_synth[:, Symbol(:x, j)] = X_new[:, j]
        end
    end
    X_synth
end

function smote_obs(dat::DataFrames.DataFrame, pct = 200, k = 5, column_names = names(dat))
    if pct < 1
        @warn("Percent over-sampling cannot be less than 1. Setting `pct` to 1.")
    end
    if pct < 1
        pct = 1
    end

    if k < 1
        @warn("k cannot be less than 1. Setting `k` to 1.")
    end
    if k < 1
        k = 1
    end

    n, p = size(dat)

    # Calling function has outcome variable in last column
    factor_indcs = factor_columns(dat)
    X = dataframe_to_matrix(dat, factor_indcs, n, p)

    # When pct < 100, only a percentage of cases will be SMOTEd
    if pct < 100
        n_needed = max(k + 1, floor(Int, (pct/100) * n))
    end
    indcs = Colon()
    if pct < 100
        indcs = StatsBase.sample(1:n, n_needed)
    end
    X = X[indcs, :]
    if pct < 100
        # @info(
        #     string("pct is < 100, so only a percentage of cases will be SMOTEd"),
        #     pct,
        #     k,
        #     k+1,
        #     n,
        #     n_needed,
        #     )
    end
    if pct < 100
        pct = 100
    end

    n, p = size(X)

    ranges = column_ranges(X)

    n_obs = round(Int, floor(pct/100))   # num. of artificial ex for each member of X
    X_new = zeros(n_obs * n, p)

    for i = 1:n
        # the k nearest neighbors of case X[i, ]
        xd = rscale(X, X[i, :], ranges)

        for col in factor_indcs
            xd[:, col] = map(x -> x == 0.0 ? 1.0 : 0.0, xd[:, col])
        end

        dd = xd.^2 * ones(p)
        last_idx = (length(dd) ≤ k + 1) ? length(dd) : (k + 1)         # HACK: Find out why `dd` is sometimes less than k+1
        #last_idx = k+1

        if last_idx < k+1
            @warn(
                string("Constraint applied. "),
                i,
                k,
                k+1,
                last_idx,
                )
        end
        k_nns = sortperm(dd)[2:last_idx]

        for l = 1:n_obs
            n_neighbors = (length(k_nns) == k) ? k : length(k_nns)
            neighbor = StatsBase.sample(1:n_neighbors)

            # the attribute values of generated case
            difs = X[k_nns[neighbor], :] - X[i, :]
            X_new[(i - 1) * n_obs + l, :] = X[i, :] + rand() * difs

            # For each of the factor variables, sample at random the original value
            # of Person i or the value that one of Person i's nearest neighbors has.
            for col in factor_indcs
                X_new[(i - 1) * n_obs + l, col] = StatsBase.sample(vcat(X[k_nns[neighbor], col], X[i, col]))
            end
        end
    end
    X_newdf = matrix_to_dataframe(X_new, dat, factor_indcs)
    DataFrames.rename!(X_newdf, column_names)
    X_newdf
end

# This version of the function is to be used when we have no factor
# variables. And it assumes input is simply a numeric matrix, where
# the last column is the outcome (or target) variable.
# NOTE: `pct` is the percent of positive examples relative to total
# sample size to be returned.
function smote_obs(X::Array{S, 2}, pct = 200, k = 5) where {S <: Real}
    if pct < 1
        @warn("Percent over-sampling cannot be less than 1. Setting `pct` to 1.")
    end
    if pct < 1
        pct = 1
    end

    if k < 1
        @warn("k cannot be less than 1. Setting `k` to 1.")
    end
    if k < 1
        k = 1
    end

    n, p = size(X)

    # When pct < 100, only a percentage of cases will be SMOTEd
    if pct < 100
        n_needed = max(k + 1, floor(Int, (pct/100) * n))
    end
    indcs = Colon()
    if pct < 100
        indcs = StatsBase.sample(1:n, n_needed)
    end
    X = X[indcs, :]
    if pct < 100
        # @info(
        #     string("pct is < 100, so only a percentage of cases will be SMOTEd"),
        #     pct,
        #     k,
        #     k+1,
        #     n,
        #     n_needed,
        #     )
    end
    if pct < 100
        pct = 100
    end

    n, p = size(X)

    ranges = column_ranges(X)
    n_obs = floor(Int, pct/100)   # num. of artificial ex for each member of X
    X_new = zeros(n_obs * n, p)

    for i = 1:n
        # The k nearest neighbors of case X[i, ]
        xd = rscale(X, X[i, :], ranges)

        dd = xd.^2 * ones(p)
        last_idx = (length(dd) ≤ k + 1) ? length(dd) : (k + 1)         # HACK: Find out why `dd` is sometimes less than k+1
        #last_idx = k+1

        if last_idx < k+1
            @warn(
                string("Constraint applied. "),
                i,
                k,
                k+1,
                last_idx,
                )
        end

        k_nns = sortperm(dd)[2:last_idx]

        for l = 1:n_obs
            n_neighbors = (length(k_nns) == k) ? k : length(k_nns)
            neighbor = nothing
            neighbor = StatsBase.sample(convert(Array, 1:n_neighbors))
            difs = X[k_nns[neighbor], :] - X[i, :]
            X_new[(i - 1) * n_obs + l, :] = X[i, :] + rand() * difs
        end
    end
    return X_new
end

# m = 150
# X = rand(m, 10)
# y = ones(m)
# X = hcat(X, y)
#
# smote_obs(X, 11)

function cases_needed(y::Array{T, 1}, prop = 0.5) where {T <: Real}
    pos_val = one(T)
    n_minority = count(x -> x == pos_val, y)
    n = length(y)
    res = round(Int, (prop * n) - n_minority)
    res
end

function pct_needed(y::Array{T, 1}, prop = 0.5) where {T <: Real}
    numer = cases_needed(y, prop)
    denom = count(x -> x == 1, y)
    res = 100 * numer/denom
    res
end

# w1 = [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
#
# cases_needed(w1)
# pct_needed(w1)
#
#
#
# X = randn(100, 10)
# y = vcat(zeros(90), ones(10))
# ub_smote(X, y)


text:
mutable struct ncupop_type
    supp # support data
    coe # coefficient data
    partition # the first 'partition' variables commutes with the remaining variables
    constraint # "projection" or "unipotent"
    basis # basis
    obj # "eigen" or "trace"
    ksupp # extending support at the k-th step
    sb # sizes of different blocks
    numb # numbers of different blocks
end

"""
    opt,data = nctssos_first(f::Polynomial{false, T} where T<:Number, x::Vector{PolyVar{false}};
        newton=true, reducebasis=true, TS="block", obj="eigen", merge=false, md=3, solve=true, QUIET=false)

Compute the first step of the NCTSSOS hierarchy for unconstrained noncommutative polynomial optimization.
If `newton=true`, then compute a monomial basis by the Newton chip method.
If `reducebasis=true`, then remove monomials from the monomial basis by diagonal inconsistency.
If `TS="block"`, use maximal chordal extensions; if `TS="MD"`, use approximately smallest chordal
extensions. If obj="eigen", minimize the eigenvalue; if obj="trace", then minimize the trace.
If `merge=true`, perform the PSD block merging. Return the optimum and other auxiliary data.

# Arguments
- `f`: the objective function for unconstrained noncommutative polynomial optimization.
- `x`: the set of noncommuting variables.
- `md`: the tunable parameter for merging blocks.
"""
function nctssos_first(f::Polynomial{false, T} where T<:Number, x::Vector{PolyVar{false}}; d=0, newton=true, reducebasis=true,
    monosquare=true, TS="block", obj="eigen", partition=0, constraint=nothing, merge=false, md=3, solve=true, QUIET=false)
    n,supp,coe = poly_info(f, x)
    opt,data = nctssos_first(supp, coe, n, d=d, newton=newton, reducebasis=reducebasis, monosquare=monosquare, TS=TS, obj=obj,
    partition=partition, constraint=constraint, merge=merge, md=md, solve=solve, QUIET=QUIET)
    return opt,data
end

function nctssos_first(supp::Vector{Vector{UInt16}}, coe, n::Int; d=0, newton=true, reducebasis=true, monosquare=true,
    TS="block", obj="eigen", partition=0, constraint=nothing, merge=false, md=3, solve=true, QUIET=false)
    println("************************NCTSSOS************************")
    println("NCTSSOS is launching...")
    if d == 0
        d = Int(maximum(length.(supp))/2)
    end
    if obj == "trace"
        supp,coe = cyclic_canon(supp, coe)
    else
        supp,coe = sym_canon(supp, coe)
    end
    if newton == true && d == 0
        if obj == "trace"
            basis = newton_cyclic(supp, n, d)
        else
            basis = newton_ncbasis(supp)
        end
    else
        basis = get_ncbasis(n, d)
    end
    if partition > 0
        ind = [_comm(basis[i], partition) == basis[i] for i=1:length(basis)]
        basis = basis[ind]
    end
    if constraint != nothing
        ind = [findfirst(j -> basis[i][j]== basis[i][j+1], 1:length(basis[i])-1) == nothing for i=1:length(basis)]
        basis = basis[ind]
    end
    ksupp = copy(supp)
    if monosquare == true && TS != false
        if obj == "trace"
            append!(ksupp, [_cyclic_canon([basis[i][end:-1:1]; basis[i]]) for i=1:length(basis)])
        else
            append!(ksupp, [[basis[i][end:-1:1]; basis[i]] for i=1:length(basis)])
        end
        if partition > 0
            ksupp = _comm.(ksupp, partition)
        end
        if constraint != nothing
            reduce_cons!.(ksupp, constraint = constraint)
        end
        sort!(ksupp)
        unique!(ksupp)
    end
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    blocks,cl,blocksize,sb,numb,_ = get_ncblocks(ksupp, basis, TS=TS, QUIET=QUIET, merge=merge, md=md, obj=obj, partition=partition, constraint=constraint)
    if reducebasis == true && obj == "eigen" && d == 0
        psupp = copy(supp)
        psupp = psupp[is_sym.(psupp)]
        push!(psupp, UInt16[])
        basis,flag = reducebasis!(psupp, basis, blocks, cl, blocksize)
        if flag == 1
            ksupp = copy(supp)
            if monosquare == true
                if obj == "trace"
                    append!(ksupp, [_cyclic_canon([basis[i][end:-1:1]; basis[i]]) for i=1:length(basis)])
                else
                    append!(ksupp, [[basis[i][end:-1:1]; basis[i]] for i=1:length(basis)])
                end
                sort!(ksupp)
                unique!(ksupp)
            end
            blocks,cl,blocksize,sb,numb,_ = get_ncblocks(ksupp, basis, TS=TS, QUIET=QUIET, merge=merge, md=md, obj=obj)
        end
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure. The maximal size of blocks is $mb.")
    end
    opt,ksupp = ncblockupop(supp, coe, basis, blocks, cl, blocksize, QUIET=QUIET, obj=obj, partition=partition, constraint=constraint, solve=solve)
    data = ncupop_type(supp, coe, partition, constraint, basis, obj, ksupp, sb, numb)
    return opt,data
end

function poly_info(f, x)
    n = length(x)
    mon = monomials(f)
    coe = coefficients(f)
    lm = length(mon)
    supp = [UInt16[] for i=1:lm]
    for i = 1:lm
        ind = mon[i].z .> 0
        vars = mon[i].vars[ind]
        exp = mon[i].z[ind]
        for j = 1:length(vars)
            k = ncbfind(x, n, vars[j], rev=true)
            append!(supp[i], k*ones(UInt16, exp[j]))
        end
    end
    return n,supp,coe
end

"""
    opt,data = nctssos_higher!(data, TS="block", merge=false, md=3, QUIET=false)

Compute higher steps of the NCTSSOS hierarchy.
Return the optimum and other auxiliary data.
"""
function nctssos_higher!(data::ncupop_type; TS="block", merge=false, md=3, solve=true, QUIET=false)
    supp = data.supp
    basis = data.basis
    coe = data.coe
    partition = data.partition
    constraint = data.constraint
    obj = data.obj
    ksupp = data.ksupp
    sb = data.sb
    numb = data.numb
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    blocks,cl,blocksize,sb,numb,status = get_ncblocks(ksupp, basis, sb=sb, numb=numb, TS=TS, QUIET=QUIET, merge=merge, md=md, obj=obj, partition=partition, constraint=constraint)
    opt = nothing
    if status == 1
        if QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure. The maximal size of blocks is $mb.")
        end
        opt,ksupp = ncblockupop(supp, coe, basis, blocks, cl, blocksize, QUIET=QUIET, obj=obj, partition=partition, constraint=constraint, solve=solve)
    end
    data.ksupp = ksupp
    data.sb = sb
    data.numb = numb
    return opt,data
end

function cc(a::Vector{UInt16}, n::Int)
    ua = unique(a)
    ca = zeros(UInt8, n)
    for i = 1:length(ua)
        ca[ua[i]] = count(x->isequal(ua[i], x), a)
    end
    return ca
end

function remove(csupp, dw, n)
    model = Model(optimizer_with_attributes(Mosek.Optimizer))
    set_optimizer_attribute(model, MOI.Silent(), true)
    t = @variable(model)
    alpha = @variable(model, [1:n])
    @constraint(model, [i=1:length(csupp)], alpha'*(csupp[i].-2*dw)<=t)
    @objective(model, Min, t)
    optimize!(model)
    if objective_value(model) >= 0
        return true
    else
        return false
    end
end

function permutation(a)
    b = sparse(a)
    ua = convert(Vector{UInt16}, b.nzind)
    na = convert(Vector{UInt16}, b.nzval)
    return _permutation(ua, na)
end

function _permutation(ua, na)
    if !isempty(ua)
        perm = Vector{UInt16}[]
        for i = 1:length(ua)
            nua = copy(ua)
            nna = copy(na)
            if na[i] == 1
                deleteat!(nua, i)
                deleteat!(nna, i)
            else
                nna[i] -= 1
            end
            temp = _permutation(nua, nna)
            push!.(temp, ua[i])
            append!(perm, temp)
        end
        return perm
    else
        return [UInt16[]]
    end
end

function newton_cyclic(supp, n, d)
    pbasis = get_basis(n,d)
    basis = [UInt16[]]
    csupp = cc.(supp, n)
    pushfirst!(csupp, zeros(UInt8, n))
    sort!(csupp)
    unique!(csupp)
    for i = 2:size(pbasis,2)
        if remove(csupp, pbasis[:,i], n)
            append!(basis, permutation(pbasis[:,i]))
        end
    end
    sort!(basis)
    return basis
end

function get_basis(n, d)
    lb = binomial(n+d, d)
    basis = zeros(UInt8, n, lb)
    i = 0
    t = 1
    while i < d+1
        t += 1
        if basis[n,t-1] == i
           if i < d
              basis[1,t] = i+1
           end
           i += 1
        else
            j = findfirst(x->basis[x,t-1]!=0, 1:n)
            basis[:,t] = basis[:,t-1]
            if j == 1
               basis[1,t] -= 1
               basis[2,t] += 1
            else
               basis[1,t] = basis[j,t] - 1
               basis[j,t] = 0
               basis[j+1,t] += 1
            end
        end
    end
    return basis
end

function _cyclic_canon(a::Vector{UInt16})
    if isempty(a)
        return a
    else
        return minimum([[a[i+1:length(a)]; a[1:i]] for i=0:length(a)-1])
    end
end

function cyclic_canon(supp, coe)
    nsupp = [min(_cyclic_canon(word), _cyclic_canon(reverse(word))) for word in supp]
    sort!(nsupp)
    unique!(nsupp)
    l = length(nsupp)
    ncoe = zeros(l)
    for i = 1:length(supp)
        Locb = ncbfind(nsupp, l, _cyclic_canon(_sym_canon(supp[i])))
        ncoe[Locb] += coe[i]
    end
    return nsupp,ncoe
end

function _sym_canon(a::Vector{UInt16})
    i = 1
    while i <= Int(ceil((length(a)-1)/2))
        if a[i] < a[end+1-i]
            return a
        elseif a[i] > a[end+1-i]
            return reverse(a)
        else
            i += 1
        end
    end
    return a
end

function is_sym(a::Vector{UInt16})
    l = Int(ceil((length(a)-1)/2))
    return isequal(a[1:l], a[end:-1:end-l+1])
end

function sym_canon(supp, coe)
    nsupp = copy(supp)
    nsupp = _sym_canon.(nsupp)
    sort!(nsupp)
    unique!(nsupp)
    l = length(nsupp)
    ncoe = zeros(l)
    for i = 1:length(supp)
        Locb = ncbfind(nsupp, l, _sym_canon(supp[i]))
        ncoe[Locb] += coe[i]
    end
    return nsupp,ncoe
end

function get_ncbasis(n, d; ind=UInt16[i for i=1:n])
    basis = [UInt16[]]
    for i = 1:d
        append!(basis, _get_ncbasis_deg(n, i, ind=ind))
    end
    return basis
end

function _get_ncbasis_deg(n, d; ind=UInt16[i for i=1:n])
    if d > 0
        basis = Vector{UInt16}[]
        for i = 1:n
            temp = _get_ncbasis_deg(n, d-1, ind=ind)
            push!.(temp, ind[i])
            append!(basis, temp)
        end
        return basis
    else
        return [UInt16[]]
    end
end

function newton_ncbasis(supp)
    nbasis = [UInt16[]]
    for bi in supp
        if iseven(length(bi))
            k = Int(length(bi)/2)
            w = bi[end-k+1:end]
            if isequal(w, bi[k:-1:1])
                for j = 1:k
                    push!(nbasis, w[end-j+1:end])
                end
            end
        end
    end
    sort!(nbasis)
    unique!(nbasis)
    return nbasis
end

function ncbfind(A, l, a; rev=false)
    low = 1
    high = l
    while low <= high
        mid = Int(ceil(1/2*(low+high)))
        if isequal(A[mid], a)
           return mid
        elseif isless(A[mid], a)
            if rev == false
                low = mid+1
            else
                high = mid-1
            end
        else
            if rev == false
                high = mid-1
            else
                low = mid+1
            end
        end
    end
    return 0
end

function get_ncgraph(ksupp, basis; obj="eigen", partition=0, constraint=nothing)
    lb = length(basis)
    G = SimpleGraph(lb)
    lksupp = length(ksupp)
    for i = 1:lb, j = i+1:lb
        bi = [basis[i][end:-1:1]; basis[j]]
        bi = reduce!(bi, obj=obj, partition=partition, constraint=constraint)
        if ncbfind(ksupp, lksupp, bi) != 0
           add_edge!(G, i, j)
        end
    end
    return G
end

function get_ncblocks(ksupp, basis; sb=[], numb=[], TS="block", obj="eigen", partition=0, constraint=nothing, QUIET=true, merge=false, md=3)
    if TS == false
        blocksize = [length(basis)]
        blocks = [[i for i=1:length(basis)]]
        cl = 1
    else
        G = get_ncgraph(ksupp, basis, obj=obj, partition=partition, constraint=constraint)
        if TS == "block"
            blocks = connected_components(G)
            blocksize = length.(blocks)
            cl = length(blocksize)
        else
            blocks,cl,blocksize = chordal_cliques!(G, method=TS)
            if merge == true
                blocks,cl,blocksize = clique_merge!(blocks, d=md, QUIET=true)
            end
        end
    end
    nsb = sort(unique(blocksize), rev=true)
    nnumb = [sum(blocksize.== i) for i in nsb]
    if isempty(sb) || nsb!=sb || nnumb!=numb
        status = 1
        if QUIET == false
            println("------------------------------------------------------")
            println("The sizes of PSD blocks:\n$nsb\n$nnumb")
            println("------------------------------------------------------")
        end
    else
        status = 0
        println("No higher NCTSSOS hierarchy!")
    end
    return blocks,cl,blocksize,nsb,nnumb,status
end

function ncblockupop(supp, coe, basis, blocks, cl, blocksize; QUIET=true, obj="eigen", partition=0, constraint=nothing, solve=true)
    ksupp = Vector{Vector{UInt16}}(undef, Int(sum(Int.(blocksize).^2+blocksize)/2))
    k = 1
    for i = 1:cl, j = 1:blocksize[i], r = j:blocksize[i]
        @inbounds bi = [basis[blocks[i][j]][end:-1:1]; basis[blocks[i][r]]]
        @inbounds ksupp[k] = bi
        k += 1
    end
    ksupp = reduce!.(ksupp, obj=obj, partition=partition, constraint=constraint)
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET == false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        cons = [AffExpr(0) for i=1:lksupp]
        for i = 1:cl
            bs = blocksize[i]
            if bs == 1
               @inbounds pos = @variable(model, lower_bound=0)
               @inbounds bi = [basis[blocks[i][1]][end:-1:1]; basis[blocks[i][1]]]
               bi = reduce!(bi, obj=obj, partition=partition, constraint=constraint)
               Locb = ncbfind(ksupp,lksupp,bi)
               @inbounds add_to_expression!(cons[Locb], pos)
            else
               @inbounds pos = @variable(model, [1:bs, 1:bs], PSD)
               for j = 1:blocksize[i], r = j:blocksize[i]
                   @inbounds bi = [basis[blocks[i][j]][end:-1:1]; basis[blocks[i][r]]]
                   bi = reduce!(bi, obj=obj, partition=partition, constraint=constraint)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   if j == r
                       @inbounds add_to_expression!(cons[Locb], pos[j,r])
                   else
                       @inbounds add_to_expression!(cons[Locb], 2, pos[j,r])
                   end
               end
            end
        end
        bc = zeros(lksupp)
        for i = 1:length(supp)
            Locb = ncbfind(ksupp, lksupp, supp[i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] = coe[i]
            end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, cons.==bc)
        @objective(model, Max, lower)
        if QUIET == false
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end

function reduce_cons!(word; constraint="unipotent")
    i = 1
    while i < length(word)
        if word[i] == word[i+1]
            deleteat!(word, i)
            if constraint == "unipotent"
                deleteat!(word, i)
            end
        else
            i += 1
        end
    end
    return word
end

function reduce!(word; obj="eigen", partition=0, constraint=nothing)
    if obj == "trace"
        word = min(_cyclic_canon(word), _cyclic_canon(reverse(word)))
    else
        if partition > 0 && constraint == nothing
            word = min(_comm(word, partition), _comm(reverse(word), partition))
        elseif partition == 0 && constraint != nothing
            cword = copy(word)
            word = min(reduce_cons!(word, constraint = constraint), reduce_cons!(reverse(cword), constraint = constraint))
        elseif partition > 0 && constraint != nothing
            word = min(reduce_cons!(_comm(word, partition), constraint = constraint), reduce_cons!(_comm(reverse(word), partition), constraint = constraint))
        else
            word = _sym_canon(word)
        end
    end
    return word
end

function _comm(word::Vector{UInt16}, partition)
    ind1 = word .<= partition
    ind2 = word .> partition
    return [word[ind1]; word[ind2]]
end


text:
"""
Bridges `CP.Knapsack{KCT, KVT}` to a MILP by adding the corresponding MILP
constraint.
"""
struct Knapsack2MILPBridge{KCT, KVT, T} <: MOIBC.AbstractBridge
    kp::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
    value::Union{
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}},
        Nothing,
    }
end

function MOIBC.bridge_constraint(
    ::Type{Knapsack2MILPBridge{KCT, CP.UNVALUED_KNAPSACK, T}},
    model,
    f::MOI.AbstractVectorFunction,
    s::CP.Knapsack{KCT, CP.UNVALUED_KNAPSACK, T},
) where {KCT, T}
    f_scalars = MOIU.scalarize(f)
    f_items = (KCT == CP.VARIABLE_CAPACITY_KNAPSACK) ? f_scalars[1:(end - 1)] : f_scalars
    f_capacity = (KCT == CP.VARIABLE_CAPACITY_KNAPSACK) ? f_scalars[end] : T(s.capacity)

    # Create the knapsack constraint.
    kp = MOI.add_constraint(
        model,
        dot(f_items, s.weights) - f_capacity,
        MOI.LessThan(zero(T))
    )

    return Knapsack2MILPBridge{KCT, CP.UNVALUED_KNAPSACK, T}(kp, nothing)
end

function MOIBC.bridge_constraint(
    ::Type{Knapsack2MILPBridge{KCT, CP.VALUED_KNAPSACK, T}},
    model,
    f::MOI.AbstractVectorFunction,
    s::CP.Knapsack{KCT, CP.VALUED_KNAPSACK, T},
) where {KCT, T}
    f_scalars = MOIU.scalarize(f)
    f_items = (KCT == CP.VARIABLE_CAPACITY_KNAPSACK) ? f_scalars[1:(end - 2)] : f_scalars[1:(end - 1)]
    f_capacity = (KCT == CP.VARIABLE_CAPACITY_KNAPSACK) ? f_scalars[end - 1] : T(s.capacity)
    f_value = f_scalars[end]

    # Create the knapsack constraint.
    kp = MOI.add_constraint(
        model,
        dot(f_items, s.weights) - f_capacity,
        MOI.LessThan(zero(T))
    )

    # Create the value constraint.
    value = MOI.add_constraint(
        model,
        dot(f_items, s.values) - f_value,
        MOI.EqualTo(zero(T))
    )

    return Knapsack2MILPBridge{KCT, CP.VALUED_KNAPSACK, T}(kp, value)
end

function MOI.supports_constraint(
    ::Type{Knapsack2MILPBridge{KCT, KVT, T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Knapsack{KCT, KVT, T}},
) where {KCT, KVT, T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:Knapsack2MILPBridge{KCT, KVT, T}}) where {KCT, KVT, T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{Knapsack2MILPBridge{KCT, CP.UNVALUED_KNAPSACK, T}}) where {KCT, T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOIB.added_constraint_types(::Type{Knapsack2MILPBridge{KCT, CP.VALUED_KNAPSACK, T}}) where {KCT, T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(
    ::Knapsack2MILPBridge{KCT, KVT, T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {KCT, KVT, T}
    return 1
end

function MOI.get(
    b::Knapsack2MILPBridge{KCT, KVT, T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    },
) where {KCT, KVT, T}
    return (b.value !== nothing) ? 1 : 0
end

function MOI.get(
    b::Knapsack2MILPBridge{KCT, KVT, T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    },
) where {KCT, KVT, T}
    return [b.kp]
end

function MOI.get(
    b::Knapsack2MILPBridge{KCT, KVT, T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    },
) where {KCT, KVT, T}
    return [b.value]
end


text:
using Test
using XEB

@testset "Rule-1" begin
    g = google_layout_53(2, 20)
    add_noise!(g, 1, 7)
    add_noise!(g, 1, 9)
    update_noise!(g)
    @test all(isequal(Noise), gates(g, 1)) && all(isequal(Noise), gates(g, 2))
end

@testset "Rule-2" begin
    g = google_layout_53(2, 20)
    add_noise!(g, 1, 7)
    add_noise!(g, 2, 7)
    update_noise!(g)
    @test all(isequal(Noise), gates(g, 1)) && all(isequal(Noise), gates(g, 2))
end

g = google_layout_53(2, 20)
add_noise!(g, 1, 7)
add_noise!(g, 1, 9)
update_noise!(g)
@test all(isequal(Noise), gates(g, 1)) && all(isequal(Noise), gates(g, 2))
XEB.update_id!(g)
gates(g, 2)


text:
# heat.jl
##########################
using FFTW, DSP, LinearAlgebra, Plots

# initial data
N=256;
n=128;
M=1000;
u=zeros(ComplexF64, N,M);
for p=65:192
    u[p,1]=1;
end

# setting
E=zeros(ComplexF64, N, N);
for k=1:N
    E[k,k]=1
end

t=1/10;

# sqrt(M) times Fourier matrix
omega=exp(im*2*pi/N);
F=zeros(ComplexF64, N, N);
for k=1:N
    for l=1:N
    F[k,l]=omega^(-(k-1)*(l-1));
    end
end

D=zeros(ComplexF64, N, N)
for k=1:n
    D[k,k]=(k-1)^2;
    D[k+n,k+n]=(n-k)^2;
end

G=-(4*pi^2*t*N^(-3))*F'*D*F
A=E+G+G^2/2+G^3/6;

for l=2:M
    u[:,l]=A*u[:,l-1];
end

v=zeros(ComplexF64, N,M)
for l=1:M
    v[:,l]=u[:,M-l+1];
end

p=plot(real.(u'),
   title="heat flow",
    st=:heatmap,
    color=:jet,
    yaxis="time t",
    xaxis="position x")
plot!(yticks = ([0 1000;], [0 100]))
plot!(xticks = ([0 256;], [0 256]))

q=surface(real.(v),
        title="heat flow",
        camera=(80,60),
        color=:jet,
        xaxis="time t",
        yaxis="position x")
plot!(xticks = ([0 1000;], [100 0]))
plot!(yticks = ([0 256;], [0 256]))

plot(p, q,
     layout=(1,2),
     size=(1600,600),
     margin=Plots.Measures.Length(:mm, 10.0))
savefig("heat.png")


text:
module OmegaTest
using Test

export isinferred

function isinferred(f, args...; allow = Union{})
  ret = f(args...)
  inftypes = Base.return_types(f, Base.typesof(args...))
  rettype = ret isa Type ? Type{ret} : typeof(ret)
  rettype <: allow || rettype == Test.typesubtract(inftypes[1], allow)
end

end # module


text:
# ------------------------------------------------------------------------------------------
# # Motivation
#
# Hello, and welcome! We're excited to be your gateway into machine learning. ML is a
# rapidly growing field that's buzzing with opportunity. Why? Because the tools and skills
# employed by ML specialists are extremely powerful and allow them to draw conclusions from
# large data sets quickly and with relative ease.
#
# Take the Celeste project, for example. This is a project that took 178 **tera**bytes of
# data on the visible sky and used it to catalogue 188 millions stars and galaxies.
# "Cataloguing" these stars meant identifying characteristics like their locations, colors,
# sizes, and morphologies. This is an amazing feat, *especially* because this entire
# calculation took under 15 minutes.
#
# <img src="data/Celeste.png" alt="Drawing" style="width: 1000px;"/>
#
# How are Celeste's calculations so fast? To achieve performance on this scale, the Celeste
# team uses the Julia programming language to write their software and supercomputers from
# Lawrence Berkeley National Lab's NERSC as their hardware. In this course, we unfortunately
# won't be able to give you access to a top 10 supercomputer, but we will teach you how to
# use Julia!
#
# We're confident that this course will put you on your way to understanding many of the
# important concepts and "buzz words" in ML. To get you started, we'll teach you how to
# teach a machine to tell the difference between images of apples and bananas, i.e to
# **classify** images as being one or the other type of fruit.
#
# Like Project Celeste, we'll use the [Julia programming language](https://julialang.org/)
# to do this. In particular, we'll be working in [Jupyter notebooks](http://jupyter.org/)
# like this one! (Perhaps you already know that the "ju" in Jupyter comes from Julia.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## What do the images we want to classify look like?
#
# We can use the `Images.jl` package in Julia to load sample images from this dataset. Most
# of the data we will use live in the `data` folder in this repository.
# ------------------------------------------------------------------------------------------

using Images  # To execute hit <shift> + enter

apple = load("data/10_100.jpg")

banana = load("data/104_100.jpg")

# ------------------------------------------------------------------------------------------
# The dataset consists of many images of different fruits, viewed from different positions.
# These images are [available on GitHub here](https://github.com/Horea94/Fruit-Images-
# Dataset).
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## What is the goal?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The ultimate goal is to feed one of these images to the computer and for it to identify
# whether the image represents an apple or a banana!  To do so, we will **train** the
# computer to learn **for itself** how to
# distinguish the two images.
#
# The following notebooks will walk you step by step through the underlying math and machine
# learning concepts you need to know in order to accomplish this classification.
#
# They alternate between two different types of notebooks: those labelled **ML** (Machine
# Learning), which are designed to give a high-level overview of the concepts we need for
# machine learning, but which gloss over some of the technical details; and those labelled
# **Tools**, which dive into the details of coding in Julia that will be key to actually
# implement the machine learning algorithms ourselves.
#
# The notebooks contain many **Exercises**. By doing these exercises in Julia, you will
# learn the basics of machine learning!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Course outline
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The course notebooks are listed below. We recommend that you follow them in order.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# - [01. Representing data in a computer](01. ML - Representing data in a computer.ipynb)
# - [02. Using arrays to store data](02. Tools - Using arrays to store data.ipynb)
#
# - [03. Representing data with models](03. ML - Representing data with models.ipynb)
# - [04. Functions](04. Tools - Functions.ipynb)
# - [05. Building models](05. ML - Building models.ipynb)
# - [06. Adding a function parameter](06. Tools - Adding a function parameter.ipynb)
# - [07. Model complexity](07. ML - Model complexity.ipynb)
# - [08. Multiple function parameters](08. Tools - Multiple function parameters.ipynb)
# - [09. What is learning](09. ML - What is learning.ipynb)
# - [10. Minimizing functions - how a computer learns](10. Tools - Minimizing functions -
# how a computer learns.ipynb)
# - [11. Intro to Neurons](11. ML - Intro to Neurons.ipynb)
# - [12. Learning with a single neuron](12. Tools - Learning with a single neuron.ipynb)
# - [13. Intro to Flux.jl](13. ML - Intro to Flux.jl.ipynb)
# - [14. Learning with a single neuron using Flux.jl](14. Tools - Learning with a single
# neuron using Flux.jl.ipynb)
# - [15. Intro to neural networks](15. ML - Intro to neural networks.ipynb)
# - [16. Using Flux to build a single layer neural net](16. Tools - Using Flux to build a
# single layer neural net.ipynb)
# - [17. Introduction to deep learning](17. ML - Introduction to deep learning.ipynb)
# - [18. Multi-layer neural networks with Flux](18. Tools - Multi-layer neural networks with
# Flux.ipynb)
# - [19. Recognizing handwriting using a neural network](19. Recognizing handwriting using a
# neural network.ipynb)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's get started with [Representing data in a computer](01. ML - Representing data in a
# computer.ipynb
# )!
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/00. Motivation - Classifying fruit.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
import Base: error, convert, show

export  Lexer, Token,
        lex, lex_end, scan_number, scan_string,
        ignore, ignore_whitespace,
        next_char, backup_char, peek_char, current_char,
        accept_char, accept_char_run,
        emit_token, lexeme,
        next_token, drain

include("util.jl")

const EOFChar = Char(0xC0) # Using illegal UTF8 as sentinel

# Because it is practical to use single chars as tokens sometimes
TokenType(ch::Char) = TokenType(Int(ch))
Char(t::TokenType)  = Char(Int(t))

"Keeps track of string of code we want to turn into array of tokens"
mutable struct Lexer
	input	:: String # string being scanned
	start	:: Int    # start position of this item (lexeme)
	pos		:: Int    # current position in the input
    tokens  :: Channel{Token}
	function Lexer(input::String)
		l = new(input, 1, 0, Channel{Token}(32))
		return l
	end
end

function show(io::IO, l::Lexer)
    print(io, "Lexer($(l.start), $(l.pos))")
end

"Advance the lexers position in the input"
function next_char(l::Lexer)
	l.pos = nextind(l.input, l.pos)
    if l.pos > lastindex(l.input)
		return EOFChar
	end
	return l.input[l.pos]
end

"Go one character back in input string"
function backup_char(l::Lexer)
	l.pos = prevind(l.input, l.pos)
	return l.input[l.pos]
end

"Check what the next character will be"
function peek_char(l::Lexer)
	pos = nextind(l.input, l.pos)
    if pos > lastindex(l.input)
		return EOFChar
	end
	return l.input[pos]
end

function current_char(l::Lexer)
	if l.pos < 1
		error("Can't ask for current char before first char has been fetched")
    elseif l.pos > lastindex(l.input)
        return EOFChar
	end
	return l.input[l.pos]
end

"Check if next character is one of among the valid ones"
function accept_char(l::Lexer, valid::AbstractString)
    if peek_char(l) in valid
		next_char(l)
        return true
	end
	return false
end

"Accept a run of characters contained withing array of valid chars"
accept_char_run(l::Lexer, valid::String) = accept_char_run(ch->ch in valid, l)

"Accept characters which `pred` evaluate to true. E.g. `accept_char_run(l, isdigit)`"
function accept_char_run(pred::Function, l::Lexer)
	while pred(peek_char(l))
        next_char(l)
    end
end

"Get lexeme that has been lexed thus far"
function lexeme(l::Lexer)
    if l.start > l.pos
        ""
    else
	    l.input[l.start:l.pos]
    end
end

"Send token of type `t` with lexeme `s` to channel `l.tokens`"
function emit_token(l::Lexer, t::TokenType, s::AbstractString)
    token = Token(t, s)
    put!(l.tokens, token)
    l.start = l.pos + 1
end

emit_token(l::Lexer, t::TokenType) = emit_token(l, t, lexeme(l))

"Skip the current token. E.g. because it is whitespace"
function ignore(l::Lexer)
	l.start = l.pos + 1
end

"Skip whitespace in input"
function ignore_whitespace(l::Lexer)
	while isspace(peek_char(l))
		next_char(l)
	end
	l.start = l.pos + 1
end

"Return this from a lexer state when there is an error"
function error(l::Lexer, error_msg::String)
	token = Token(ERROR, error_msg)
    put!(l.tokens, token)
    return lex_end
end

################### Scan Common Types ###################
"""
Advancing the position in the lexer input passed any number.
It means `lexeme(lexer)` will return a number.
To emit a number token you can just write:

    if scan_number(l)
        emit_token(l, NUMBER)
    end

The reason for doing it this way, is that it retains flexibility
in with which lexer state should be associated with lexing a number.
Several kinds of lexer will need to lex a number. So we want
reusable functions.
"""
function scan_number(l::Lexer)
	# leading sign is optional, but we'll accept it
	accept_char(l, "-+")

	# Could be a hex number, assume it is not first
	digits = "0123456789"
	if accept_char(l, "0") && accept_char(l, "xX")
		digits *= "abcdefABCDEF"
	end
	accept_char_run(l, digits)
	if accept_char(l, ".")
		accept_char_run(l, digits)
	end
	if accept_char(l, "eE")
		accept_char(l, "-+")
		accept_char_run(l, "0123456789")
	end
    return true
end

"""
Advancing the position in the lexer input passed any quoted string.
It means `lexeme(lexer)` will return a quoted string. Example:

    if scan_string(l)
        emit_token(l, STRING)
    end
"""
function scan_string(l::Lexer)
	if !accept_char(l, "\"")
        return false
    end
	while true
		ch = next_char(l)
        if ch == '\\'
            c = next_char(l)
            if c == EOFChar || c == '\n'
                return false
            end
		elseif ch == '"'
            break
		elseif ch == EOFChar
			return false
		end
	end
	return true
end

function scan_identifier(l::Lexer)
    ch = peek_char(l)
    if !isletter(ch)
        return false
    end

    while ch != EOFChar && (isalnum(ch) || ch == '_')
        next_char(l)
        ch = peek_char(l)
    end

    return true
end

################### Lexer Common ###################
"Marker for indicating there is no more input. Since we don't want to use nil in Julia"
function lex_end(l::Lexer)
	return lex_end
end

function lex(input::AbstractString, start::Function)
    l = Lexer(input)
    @async run(l, start)
    return l
end

function run(l::Lexer, start::Function)
    state = start
    while state != lex_end
        state = state(l)
    end
    close(l.tokens)
end

next_token(l::Lexer) = take!(l.tokens)
drain(l::Lexer) = collect(l.tokens)


text:
mutable struct UnfoldAssessState{T}
    weight::Float64
    retvals::Vector{T}
    state::T
end

function process_new!(gen_fn::Unfold{T,U}, params::Tuple, choices::ChoiceMap,
                      key::Int, state::UnfoldAssessState{T}) where {T,U}
    local new_state::T
    kernel_args = (key, state.state, params...)
    submap = get_submap(choices, key)
    (weight, new_state) = assess(gen_fn.kernel, kernel_args, submap)
    state.weight += weight
    state.retvals[key] = new_state
    state.state = new_state
end

function assess(gen_fn::Unfold{T,U}, args::Tuple, choices::ChoiceMap) where {T,U}
    len = args[1]
    init_state = args[2]
    params = args[3:end]
    state = UnfoldAssessState{T}(0., Vector{T}(undef,len), init_state)
    for key=1:len
        process_new!(gen_fn, params, choices, key, state)
    end
    (state.weight, PersistentVector{T}(state.retvals))
end


text:
# This is to make sure the swap operation does not make a mess
function swaptest()
    @testset "swap" begin
        r=Repo{String}(100,3)
        for i in 1:500           # We shuffle the array
            f=mod1(i*12345,100)  # Make some quick random numbers
            t=mod1(i*54321,100)
            Repos._swap(r,f,t)
        end
        for i=1:100          # Now test if the state is still consistent
            @test r.pos[r.ids[i]]==i
        end
    end
end

# Test for setclass and getclass functionality
function classtest()
    @testset "class" begin
        r=Repo{String}(4,4)
        setclass!(r,4,1)
        setclass!(r,3,2)
        setclass!(r,2,3)
        setclass!(r,1,4)
        @test getclass(r,1)==4
        @test getclass(r,2)==3
        @test getclass(r,3)==2
        @test getclass(r,4)==1
    end
end


text:
const brown_dennis = let res_init=zeros(20), jac_init=zeros(20,4), x_init=[25., 5., -5., -1.]

    function res(x, r)
        n = length(r)
        for i = 1:n
            ti = 0.2e0 * Float64(i)
            r[i] = (x[1] + x[2] * ti - exp(ti)) ^ 2 +
                (x[3] + x[4] * sin(ti) - cos(ti))^2
        end
        return r
    end

    function jac(x, j)
        n = size(j)[1]
        for i = 1:n
            ti = 0.2e0 * Float64(i)
            j[i, 1] = 2.0e0 * (x[1] + x[2]*ti - exp(ti))
            j[i, 2] = ti * j[i, 1]
            t = sin(ti)
            j[i, 3] = 2.0e0 *(x[3] + x[4]*t - cos(ti))
            j[i, 4] = t * j[i, 3]
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("brown_dennis", res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end


text:
@testset "Input Data" begin

@testset "Testing function definitions" begin
    @test isdefined(GRASS, :SpecParams)
    @test isdefined(GRASS, :DiskParams)
    @test isdefined(GRASS, :SolarData)
end

@testset "Testing data read-in" begin
    soldat = GRASS.SolarData()
    @test !isempty(soldat.wav)
    @test !isempty(soldat.bis)
    @test length(keys(soldat.wav)) == 41
    @test keys(soldat.wav) == keys(soldat.bis)
    for k in keys(soldat.len)
        @test soldat.len[k] == size(soldat.wav[k], 2)
        @test soldat.len[k] == size(soldat.bis[k], 2)
        @test size(soldat.bis[k], 2) == size(soldat.wav[k], 2)
    end
end

end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using JuMP
using Gurobi
using Suppressor
using CSV
using JLD
using Printf
using Clustering
gurobi_env = @suppress Gurobi.Env()


text:
using TunnelIonization
using Test

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


text:
using Documenter, OpenQuantumTools
using DocumenterTools: Themes

Themes.compile(joinpath(@__DIR__, "src/assets/design.scss"))

makedocs(
    format=Documenter.HTML(
        prettyurls=false,
        assets = ["assets/design.css"],
        mathengine=MathJax(),
    ),
    sitename="OpenQuantumTools",
    pages = Any[
        "Home" => "index.md",
        "Basics" => Any[
            "Math Symbol" => "basics/math_symbol.md",
            "Hamiltonian" => "basics/hamiltonians.md",
            "Coupling" => "basics/couplings.md",
            "Bath" => "basics/bath.md",
            "Annealing" => "basics/annealing.md",
            "Solver" => "basics/solver.md"
        ],
        "Library" => Any[
            "Hamiltonian" => "lib/hamiltonian.md",
            "Matrix Utility" => "lib/matrix_util.md",
            "Coupling" => "lib/coupling.md",
            "Bath"=>"lib/bath.md",
            "Interaction"=>"lib/interactions.md",
            "Annealing"=>"lib/annealing.md",
            "Solver"=>"lib/solvers.md",
            "Callback"=>"lib/callbacks.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/USCqserver/OpenQuantumTools.jl.git",
)

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
using Documenter, Example

makedocs(
    sitename="BinaryTrees",
    repo="github.com:control13/BinaryTrees.jl.git"
)


text:
module test_millboard_functions

using Millboard
using Test

Millboard.set_table_mode(:markdown)

board = [(+) (*); (-) (/)]
@test sprint(show, table(board)) == """
|   | 1 | 2 |
|---|---|---|
| 1 | + | * |
| 2 | - | / |"""

end # module test_millboard_functions


text:
module CommutativeFunctions

export @commutative

"""
	@commutative f(x::T1, y::T2) = ...

Declare a method definition of a binary function to be commutative in the arguments.
If the two arguments are of different types, this will generate a second method with the
arguments swapped that is identical in its definition.

# Example
```jldoctest; setup=:(using CommutativeFunctions)
julia> @commutative f(x::Int, y::Float64) = x + y
f (generic function with 2 methods)

julia> f(1, 2.0)
3.0

julia> f(2.0, 1)
3.0
```
"""
macro commutative(ex)
	if ex.head != :(=)
		throw(ArgumentError("argument to @commutative must be a method definition"))
	end
	lhs = ex.args[1]
	rhs = ex.args[2]
	if lhs.head != :call
		throw(ArgumentError("unable to parse expression"))
	end
	f = lhs.args[1]
	if length(lhs.args) != 3
		throw(ArgumentError("the function $f must have two arguments"))
	end
	arg1 = lhs.args[2]
	arg2 = lhs.args[3]

	quote
		$(esc(f))($arg2, $arg1) = $(rhs.args[2]);
		$(esc(f))($arg1, $arg2) = $(rhs.args[2]);
	end
end

end


text:
using ..CUDASupport
using ..TopOpt: @init_cuda, CPU, GPU
@init_cuda()
import ..TopOpt: whichdevice
using ..GPUUtils, ForwardDiff

abstract type AbstractGPUPenalty{T} <: AbstractPenalty{T} end

whichdevice(::AbstractCPUPenalty) = CPU()
whichdevice(::AbstractGPUPenalty) = GPU()

CUDAnative.pow(d::TD, p::AbstractFloat) where {T, TV, TD <: ForwardDiff.Dual{T, TV, 1}} = ForwardDiff.Dual{T}(CUDAnative.pow(d.value, p), p * d.partials[1] * CUDAnative.pow(d.value, p - 1))

struct GPUPowerPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPUPowerPenalty)(x) = CUDAnative.pow(x, P.p)

struct GPURationalPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPURationalPenalty)(x) = x / (1 + R.p * (1 - x))

struct GPUSinhPenalty{T} <: AbstractGPUPenalty{T}
    p::T
end
@inline (P::GPUSinhPenalty)(x) = CUDAnative.sinh(R.p*x)/CUDAnative.sinh(R.p)

for T in (:PowerPenalty, :RationalPenalty)
    fname = Symbol(:GPU, T)
    @eval @inline CuArrays.cu(p::$T) = $fname(p.p)
end
CuArrays.cu(p::AbstractGPUPenalty) = p

@define_cu(IterativeSolvers.CGStateVariables, :u, :r, :c)
whichdevice(ra::RaggedArray) = whichdevice(ra.offsets)
@define_cu(RaggedArray, :offsets, :values)


text:
module Weight
    export order_weight

    function order_weight(strng)
        modweight = w -> sum( c -> c-'0', w )
        modless = (x,y) -> x[2] == y[2] ? x[1] < y[1] : x[2] < y[2]
        weights_modified = map( w -> (w,modweight(w)), split(strng))
        sort!( weights_modified, lt=modless )
        join( map( wm -> wm[1], weights_modified ), ' ' )
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Inverting for the mass-balance coefficients: b_max
error("This is not currently running...")
using  JLD, Plots, Printf, LinearAlgebra

@views av(A)    = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
@views av_xa(A) = 0.5.*(A[1:end-1,:].+A[2:end,:])
@views av_ya(A) = 0.5.*(A[:,1:end-1].+A[:,2:end])
@views inn(A)   = A[2:end-1,2:end-1]

@views function iceflow(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max)
    println("Initialising ice flow model ... ")
    # physics
    s2y      = 3600*24*365.25  # seconds to years
    rho_i    = 910.0           # ice density
    g        = 9.81            # gravity acceleration
    npow     = 3.0             # Glen's power law exponent
    a0       = 1.5e-24         # Glen's law enhancement term
    # numerics
    @assert (dx>0 && dy>0) "dx and dy need to be positive"
    nx, ny   = size(Zbed,1), size(Zbed,2) # numerical grid resolution
    @assert (nx, ny) == size(Zbed) == size(Hice) == size(Mask) "Sizes don't match"
    itMax    = 1e5             # number of iteration (max)
    nout     = 200             # error check frequency
    tolnl    = 1e-6            # nonlinear tolerance
    epsi     = 1e-4            # small number
    damp     = 0.85            # convergence accelerator (this is a tuning parameter, dependent on e.g. grid resolution)
    dtausc   = 1.0/3.0         # iterative dtau scaling
    # inversion
    tol_inv  = 1e-4            # inversion tolerance
    nt_inv   = 30              # number of inversion steps
    nout_inv = 1               # inversion plot frequency
    nsm      = 300             # number of smoothing steps
    tau1     = 5e-5            # inversion update step B_max
    tau2     = 0.1             # inversion update step z_ELA
    nsmG     = 1               # number of Gam smoothing steps
    # derived physics
    a        = 2.0*a0/(npow+2)*(rho_i*g)^npow*s2y
    lx, ly   = nx*dx, ny*dy
    # derived numerics
    xc, yc   = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny)
    xv, yv   = 0.5*(xc[1:end-1].+xc[2:end]), 0.5*(yc[1:end-1].+yc[2:end])
    (Xc,Yc)  = ([x for x=xc,y=yc], [y for x=xc,y=yc])
    cfl      = max(dx^2,dy^2)/4.1
    # array initialization
    Err      = zeros(nx  , ny  )
    dSdx     = zeros(nx-1, ny  )
    dSdy     = zeros(nx  , ny-1)
    gradS    = zeros(nx-1, ny-1)
    D        = zeros(nx-1, ny-1)
    qHx      = zeros(nx-1, ny-2)
    qHy      = zeros(nx-2, ny-1)
    dtau     = zeros(nx-2, ny-2)
    ResH     = zeros(nx-2, ny-2)
    dHdt     = zeros(nx-2, ny-2)
    Vx       = zeros(nx-1, ny-1)
    Vy       = zeros(nx-1, ny-1)
    M        = zeros(nx  , ny  )
    B        = zeros(nx  , ny  )
    H        = zeros(nx  , ny  )
    Gam      = zeros(nx  , ny  )
    # initial condition
    S        = zeros(nx  , ny  )
    # initial conditions
    B       .= Zbed
    H       .= Hice
    Yc2      = Yc .- minimum(Yc); Yc2 .= Yc2./maximum(Yc2)
    grad_b   = (1.3517 .- 0.014158.*(60.0.+Yc2*20.0))./100.0.*0.91 # From doi: 10.1017/jog.2016.75
    z_ELA    = 1300.0 .- Yc2*300.0
    B_max    = b_max.*ones(nx, ny)
    S       .= B .+ H
    if do_visu
        FS = 7
        H_v     = fill(NaN, nx, ny)
        B_max_v = fill(NaN, nx, ny)
    end
    err_inv0 = 1.0
    err_evo1=[]; err_evo2=[]
    # inversion loop
    for it_inv = 1:nt_inv
        H     .= Hice
        S     .= B .+ H
        println("forward solver:")
        # forward solver
        it = 1; err = 2*tolnl
        while err>tolnl && it<itMax
            Err   .= H
            # mass balance
            M     .= min.(grad_b.*(S .- z_ELA), B_max)
            # compute diffusivity
            dSdx  .= diff(S, dims=1)/dx
            dSdy  .= diff(S, dims=2)/dy
            gradS .= sqrt.(av_ya(dSdx).^2 .+ av_xa(dSdy).^2)
            D     .= a*av(H).^(npow+2) .* gradS.^(npow-1)
            # compute flux
            qHx   .= .-av_ya(D).*diff(S[:,2:end-1], dims=1)/dx
            qHy   .= .-av_xa(D).*diff(S[2:end-1,:], dims=2)/dy
            # update ice thickness
            dtau  .= dtausc*min.(10.0, cfl./(epsi .+ av(D)))
            ResH  .= .-(diff(qHx, dims=1)/dx .+ diff(qHy, dims=2)/dy) .+ inn(M)
            dHdt  .= dHdt.*damp .+ ResH
            H[2:end-1,2:end-1] .= max.(0.0, inn(H) .+ dtau.*dHdt)
            # apply mask
            H[Mask.==0] .= 0.0
            # update surface
            S     .= B .+ H
            if mod(it, nout)==0
                # error check
                Err .= Err .- H
                err = norm(Err)/length(Err)
                @printf(" it = %d, error = %1.2e \n", it, err)
                if isnan(err) error("NaNs") end # safeguard
            end
            it += 1
        end
        # modified from inversion (Visnjevic 2018)
        if it_inv==1 rel = 1.0 else rel = 0.2 end
        Gam .= (1.0-rel)*Gam .+ rel*(Hice .- H)
        for ism = 1:nsmG
            smooth!(Gam)
        end
        # z_ELA  .= z_ELA .- tau2*Gam
        B_max  .= B_max .+ tau1*Gam
        for ism = 1:nsm
            # smooth!(z_ELA) # assumes tau2=min(dx^2,dy^2)/4.1
            smooth!(B_max) # assumes tau2=min(dx^2,dy^2)/4.1
        end
        B_max[B_max.<1e-3].=1e-3
        if it_inv==1 err_inv0 = sum(abs.(Gam)) end
        err_inv = sum(abs.(Gam))/err_inv0
        @printf("--> it_inv = %d, error = %1.2e \n\n", it_inv, err_inv)
        if (err_inv<tol_inv) break end

        if mod(it_inv, nout_inv)==0 && do_visu
            H_v.=H;         H_v[Mask.==0].=NaN
            B_max_v.=B_max; B_max_v[Mask.==0].=NaN
            push!(err_evo1, it_inv); push!(err_evo2, err_inv)
            p1 = plot(err_evo1, err_evo2, xlims=(1, nt_inv), ylims=(1e-1, 1.0), legend=false, xlabel="# iterations", ylabel="error", xaxis=font(FS, "Courier"), yaxis=font(FS, "Courier"), linewidth=2, markershape=:circle, markersize=3, framestyle=:box, labels="max(error)")
            p2 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(Gam, dims=2)'    , c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="Gam", titlefontsize=FS, titlefont="Courier")
            p3 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(B_max_v, dims=2)', c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="B_max", titlefontsize=FS, titlefont="Courier")
            p4 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(H_v, dims=2)'    , c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="Ice thickness", titlefontsize=FS, titlefont="Courier")
            plot(p1, p2, p3, p4, size=(400,400), dpi=200); frame(anim) #background_color=:transparent, foreground_color=:white
            # savefig("../output_inv/iceflow_it_inv$(it_inv).png")
        end
    end
    # compute velocities
    Vx .= -D./(av(H) .+ epsi).*av_ya(dSdx)
    Vy .= -D./(av(H) .+ epsi).*av_xa(dSdy)
    # return as GeoArrays
    return  as_geoarray(H,  Zbed, name=:thickness),
            as_geoarray(S,  Zbed, name=:surface),
            as_geoarray(M,  Zbed, name=:smb),
            as_geoarray(Vx, Zbed, name=:vel_x, staggerd=true),
            as_geoarray(Vy, Zbed, name=:vel_y, staggerd=true)
end
# ------------------------------------------------------------------------------
include("../scripts/helpers.jl")

# load the data
print("Loading the data ... ")
Zbed, Hice, Mask, dx, dy, xc, yc = load_data(; nx=96) # nx=96,160 are included in the repo
                                                                      # other numbers will trigger a 2GB download
println("done.")

# apply some smoothing
print("Applying some smoothing ... ")
for is=1:2 # two smoothing steps
    smooth!(Zbed)
    smooth!(Hice)
end
println("done.")

# calculate mass balance coefficients for given spatial grid
grad_b, z_ELA, b_max = mass_balance_constants(xc, yc)

# handle output
do_visu = true
do_save = true

# visu and save
nx, ny = size(Hice)
if do_visu
    # ENV["GKSwstype"]="nul"
    !ispath("../output_inv") && mkdir("../output_inv")
    !ispath("../output_inv/gif_$(nx)x$(ny)") && mkdir("../output_inv/gif_$(nx)x$(ny)"); loadpath = "../output_inv/gif_$(nx)x$(ny)"; anim = Animation(loadpath,String[])
    println("Animation directory: $(anim.dir)")
end

# run the inversion SIA flow model
H, S, M, Vx, Vy = iceflow_inverse(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max; do_visu)

# visualisation
if do_visu
    gif(anim, "../output_inv/iceflow_inv_$(nx)x$(ny).gif", fps = 5)

    H_v = fill(NaN, nx, ny)
    S_v = fill(NaN, nx, ny)
    M_v = fill(NaN, nx, ny)
    V_v = fill(NaN, nx-2, ny-2)

    # outputs
    FS  = 7
    H_v.=H; H_v[Mask.==0].=NaN
    S_v.=S; S_v[Mask.==0].=NaN
    M_v.=M; M_v[Mask.==0].=NaN
    V_v.=sqrt.(av(Vx).^2 .+ av(Vy).^2); V_v[inn(H).==0].=NaN
    p1 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(S_v, dims=2)', c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="Surface elev. [m]", titlefontsize=FS, titlefont="Courier")
    p2 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(H_v, dims=2)', c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="Ice thickness [m]", titlefontsize=FS, titlefont="Courier")
    p3 = heatmap(xc[2:end-1]./1e3, reverse(yc[2:end-1])./1e3, reverse(log10.(V_v), dims=2)', c=:batlow, aspect_ratio=1, xlims=(xc[2], xc[end-1])./1e3, ylims=(yc[end-1], yc[2])./1e3, clims=(0.1, 2.0), yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="log10(vel) [m/yr]", titlefontsize=FS, titlefont="Courier")
    p4 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(M_v, dims=2)', c=:devon, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="Mass Bal. rate [m/yr]", titlefontsize=FS, titlefont="Courier")
    # display(plot(p1, p2, p3, p4, size=(400,400)))
    plot(p1, p2, p3, p4, size=(400,400), dpi=200) #background_color=:transparent, foreground_color=:white
    savefig("../output_inv/iceflow_inv_out1_$(nx)x$(ny).png")

    # error
    H_diff = Hice.-H; H_diff[Mask.==0] .= NaN
    Hice[Mask.==0] .= NaN
    H[Mask.==0]    .= NaN
    FS = 7
    p1 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(Hice, dims=2)'  , c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="H data [m]", titlefontsize=FS, titlefont="Courier")
    p2 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(H, dims=2)'     , c=:davos, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="H model [m]", titlefontsize=FS, titlefont="Courier")
    p3 = heatmap(xc./1e3, reverse(yc)./1e3, reverse(H_diff, dims=2)', c=:viridis, aspect_ratio=1, xlims=(xc[1], xc[end])./1e3, ylims=(yc[end], yc[1])./1e3, yaxis=font(FS, "Courier"), ticks=nothing, framestyle=:box, title="H (data-model) [m]", titlefontsize=FS, titlefont="Courier")
    # display(plot(p1, p2, p3, layout=(1, 3), size=(500,160)))
    plot(p1, p2, p3, layout=(1, 3), size=(500,160), dpi=200) #background_color=:transparent, foreground_color=:white
    savefig("../output_inv/iceflow_inv_out2_$(nx)x$(ny).png")
end

if do_save
    save("../output_inv/iceflow_inv_$(nx)x$(ny).jld", "Hice", convert(Matrix{Float32}, Hice),
                                                      "Mask", convert(Matrix{Float32}, Mask),
                                                      "H"   , convert(Matrix{Float32}, H),
                                                      "S"   , convert(Matrix{Float32}, S),
                                                      "M"   , convert(Matrix{Float32}, M),
                                                      "Vx"  , convert(Matrix{Float32}, Vx),
                                                      "Vy"  , convert(Matrix{Float32}, Vy),
                                                      "xc", xc, "yc", yc)
end

println("... done.")


text:
# julia -i osc.jl

# (v0.7) pkg> add https://github.com/fundamental/OSC.jl#master

#=
* Requirements
  - [BlenderPlot.jl](https://github.com/wookay/BlenderPlot.jl)

* Run OSC Server in Blender Julia Console

@pyimport osc
serv = osc.Server("bpy")
serv.start(5005)
# serv.stop()      # stop it to quit Blender gracefully

=#

using Sockets, OSC
function send_osc(oscmsg, port=5005)
    udp = UDPSocket()
    send(udp, ip"127.0.0.1", port, oscmsg.data)
end

function handle_knob(msg)
    println("handle_knob ", msg)
    if length(msg) >= 3
        oscmsg = OscMsg("/knob", "ii", Int32(msg[2]), Int32(msg[3]))
        send_osc(oscmsg)
    end
end

import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        wait(cond)
        msg = codeunits(unsafe_string(r_ecb[].message))
        handle_knob(msg)
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)


text:
module SparseMatrixLILs

import Base: getindex, size, setindex!, length, findnz, nnz, +, -, map, map!, /, \, *
export SparseMatrixLIL

"""
Store sparse matrix as as vector of sparse vectors (list of lists)
This is to enable more efficient updates (i.e. setindex!)
Each vector represents a COLUMN on the matrix
"""
immutable SparseMatrixLIL{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    data :: Vector{SparseVector{Tv,Ti}}
    m :: Int
    n :: Int
end

"""
        SparseMatrixLIL(A :: SparseMatrixLIL)

Returns A without copying.

        SparseMatrixLIL(A :: SparseMatrixCSC{Tv,Ti})

Converts to `SparseMatrixLIL` format (copying).

        SparseMatrixLIL(A :: AbstractMatrix{Tv})

Converts to `SparseMatrixLIL` format (copying).
"""
SparseMatrixLIL(A :: SparseMatrixLIL) = A
SparseMatrixLIL(A :: SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer} =
    SparseMatrixLIL{Tv,Ti}(map(v -> A[:,v], 1:size(A,2)), size(A)...)
SparseMatrixLIL(A :: AbstractMatrix{Tv}) where {Tv} =
    SparseMatrixLIL{Tv,Int}(map(v -> sparsevec(A[:,v]), 1:size(A,2)), size(A)...)

size(A::SparseMatrixLIL) = (A.m,A.n)
size(A::SparseMatrixLIL,d::Integer) = (A.m,A.n)[d]
length(A::SparseMatrixLIL) = A.m * A.n

@inline getindex(A::SparseMatrixLIL, i::Integer, j::Integer) = A.data[j][i]
getindex(A::SparseMatrixLIL, I::AbstractVector, j::Integer) = SparseMatrixLIL([A.data[j][I]])

function getindex(A::SparseMatrixLIL{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti<:Integer}
    data = Vector{SparseVector{Tv,Ti}}(length(J))
    for ji in 1:length(J)
        data[ji] = A.data[J[ji]][I]
    end
    SparseMatrixLIL{Tv,Ti}(data, length(I), length(J))
end

@inline setindex!(A::SparseMatrixLIL, x, i::Integer, j::Integer) = A.data[j][i] = x
@inline setindex!(A::SparseMatrixLIL, x::AbstractVector, I::AbstractVector, j::Integer) = A.data[j][I] = x
@inline function setindex!(A::SparseMatrixLIL{Tv,Ti}, x::AbstractMatrix, I::AbstractVector, J::AbstractVector) where {Tv,Ti<:Integer}
    for ji in 1:length(J)
        A.data[J[ji]][I] = view(x, :, ji)
    end
    x
end

function nnz(A::SparseMatrixLIL{Tv,Ti}) where {Tv,Ti<:Integer}
    n = 0
    for j = 1:A.n
        n += nnz(A.data[j])
    end
    n
end

function findnz(S::SparseMatrixLIL{Tv,Ti}) where {Tv,Ti}
    numnz = nnz(S)
    I = Vector{Ti}(numnz)
    J = Vector{Ti}(numnz)
    V = Vector{Tv}(numnz)

    count = 1
    for col = 1:S.n
        nzval = S.data[col].nzval
        nzind = S.data[col].nzind
        for u = 1:length(nzind)
            if nzval[u] != 0
                I[count] = nzind[u]
                J[count] = col
                V[count] = nzval[u]
                count += 1
            end
        end
    end

    count -= 1
    if numnz != count
        deleteat!(I, (count+1):numnz)
        deleteat!(J, (count+1):numnz)
        deleteat!(V, (count+1):numnz)
    end

    return (I, J, V)
end

map(f::Tf, A::SparseMatrixLIL) where {Tf} =
    SparseMatrixLIL(map(v -> map(f,v), A.data), A.m, A.n)

map(f::Tf, A::SparseMatrixLIL, Bs::Vararg{SparseMatrixLIL,N}) where {Tf,N} =
    SparseMatrixLIL(map((v...) -> map(f,v...), A.data, map(B -> B.data, Bs)...), A.m, A.n)

function map!(f::Tf, C::SparseMatrixLIL, A::SparseMatrixLIL, Bs::Vararg{SparseMatrixLIL,N}) where {Tf,N}
    for i = 1:C.n
        map!(f, C.data[i], A.data[i], map(B -> B.data[i], Bs)...)
    end
    C
end

broadcast(f::Tf, A::SparseMatrixLIL) where {Tf} =
    SparseMatrixLIL(map(v -> map(f,v), A.data), A.m, A.n)

function broadcast(f::Tf, A::SparseMatrixLIL, Bs::Vararg{SparseMatrixLIL,N}) where {Tf,N}
    #shape = Base.Broadcast.broadcast_indices(A.data, map(B -> B.data, Bs)...)
    #n = Base.Broadcast.broadcast_shape(size(A.data), map(B -> size(B.data), Bs)...)[1]
    #keeps, Idefaults = Base.Broadcast.map_newindexer(shape, A, Bs)
    for B in Bs
        if A.n != B.n
            error("not implemented")
        end
    end
    data = Vector{SparseVector{Tv,Ti}}(A.n)
    for i = 1:A.n
        data[i] = broadcast(f, A.data[i], map(B -> B.data[i], Bs)...)
    end
    # SparseMatrixLIL(map((v...) -> map(f,v...), A.data, map(B -> B.data, Bs)...), A.m, A.n)
    SparseMatrixLIL(data, length(data[1]), A.n)
end

function broadcast!(f::Tf, C::SparseMatrixLIL, A::SparseMatrixLIL, Bs::Vararg{SparseMatrixLIL,N}) where {Tf,N}
    for i = 1:C.n
        broadcast!(f, C.data[i], A.data[i], map(B -> B.data[i], Bs)...)
    end
    C
end

## Binary arithmetic and boolean operators
(+)(A::SparseMatrixLIL, B::SparseMatrixLIL) = map(+, A, B)
(-)(A::SparseMatrixLIL, B::SparseMatrixLIL) = map(-, A, B)
(+)(A::SparseMatrixLIL, B::Array) = Array(A) + B
(+)(A::Array, B::SparseMatrixLIL) = A + Array(B)
(-)(A::SparseMatrixLIL, B::Array) = Array(A) - B
(-)(A::Array, B::SparseMatrixLIL) = A - Array(B)

for f in (:/, :\, :*)
    if f != :/
        @eval ($f)(A::Number, B::SparseMatrixLIL) = SparseMatrixLIL(map(v -> ($f)(A,v), B.data), B.m, B.n)
    end
    if f != :\
        @eval ($f)(A::SparseMatrixLIL, B::Number) = SparseMatrixLIL(map(v -> ($f)(v,B), A.data), A.m, A.n)
    end
end

end


text:
type GuidType
    guidance_mode::UInt
    num_waypts::UInt
    current_wypt::UInt
    WaypointList::Array{WGS_Pos,1}
    tgo::Float64
    miss::Float64
    miss_vec_ned::Array{Float64,1}
    eta_dive_cmd::Float64
    eta_turn_cmd::Float64
    alt_cmd::Float64
    gamma_cmd::Float64
    Mach_cmd::Float64
end
################################################################################
## Cooperative Formation
function GuidInit(wyptlist::Array{WGS_Pos,1}, num_waypoints::Float64)
    guidance_mode = 0;
    current_wypt = 1;
    tgo = 99999.0;
    miss = 9999.0;
    miss_vec_ned = zeros(3);
    eta_dive_cmd = 0.0;
    eta_turn_cmd = 0.0;
    alt_cmd = 0.0;
    gamma_cmd = 0.0;
    Mach_cmd = 0.0;

    Guidance_Data = GuidType(guidance_mode, num_waypoints, current_wypt, wyptlist,
                             tgo, miss, miss_vec_ned, eta_turn_cmd, eta_dive_cmd,
                             alt_cmd, gamma_cmd, Mach_cmd );

    return Guidance_Data;
end

################################################################################
## Cooperative Formation
function CoopFormation(nav_state::NavStates, guidance_data::GuidType, SS::Array{Float64,2},
    missile_id::UInt,  salvo_size::UInt, message_list::Array{CoopFormation,1})

    # Computing NED to Body DCM
    Tr_BN = nav_state.euler.Tr_bn;
    Tr_NE = nav_state.position.Tr_ne;
    rm_ecef = nav_state.position.rm_ecef;
    vm_ecef = nav_state.velocity.vm_ecef;
    nn = length(message_list);

    ## Clearing the summation variables
    sumvec = zeros(3);
    sum = 0;
    ii = missile_id;
    icoop_flag = true;

    ## Gravity Vector in NED frame
    tmpvec = zeros;
    tmpvec[3] = 1.0;
    gvec_LL[1] = 0.0;
    gvec_LL[2] = 0.0;
    gvec_LL[3] = go;

    rm_ned = Tr_NE*rm_ecef;
    vm_ned = Tr_NE*vm_ecef;

    ## Computing the Guidance Command
    jj = 0;

    for jj = 1:NN

        ## Determine Communication Topology
        if ( message_list[jj].time >= 0 )
            SS[jj] = 1;
        else
            SS[jj] = 0;
        end

        ## For right now, Missile 0 is the leader.
        if( jj == 1 )
            rm_lead_ecef = message_list[jj].rm_ecef;
            vm_lead_ecef = message_list[jj].vm_ecef;
            rm_lead_ned = Tr_NE*rm_lead_ecef;
            vm_lead_ned = Tr_NE*vm_lead_ecef;

            ## Unit Velocity Vector
            unit_vmlead_ned = vm_lead_ned/norm(vm_lead_ned);

            ## Gamma and Azimuth of Lead Vector
            gamma_lead = atan2(-unit_vmlead_ned[3], sqrt( unit_vmlead_ned[1]^2 + unit_vmlead_ned[2]^2 ));
            azimuth_lead = atan2( unit_vmlead_ned[1], unit_vmlead_ned[0] );

            ## NED to Velocity Frame Transformation Matrix
            Roll_lead = 0.0;
            Tr_VN_LEAD = TR_BN(Roll_lead, gamma_lead, azimuth_lead );
            Tr_NV_LEAD = Tr_VN_LEAD';

            #delta_vel[0] = 0.0;
            #delta_vel[1] = YY(jj) - YY(missile_number);
            #delta_vel[2] = 0.0;

            delta_vel = zeros(3);
            delta_vel[1] = 0.0;
            delta_vel[2] = 200.0*sin(0.20*guidtime + (ii-1)*2.0*MC::PI/(salvo_size - 1));
            delta_vel[3] = 200.0*cos(0.20*guidtime + (ii-1)*2.0*MC::PI/(salvo_size - 1));
            delta_ned = Tr_NV_LEAD*delta_vel;

            wn = 0.5;
            zz = 0.70;
            Kpf = wn*wn;
            Kvf = 2.0*wn*zz;
            sumvec = sumvec + SS[jj]*( Kpf*(rm_lead_ned - rm_ned - delta_ned ) + Kvf*(vm_lead_ned - vm_ned) );

        end

        jj = jj + 1;
    end

    ## Guidance Command
    acc_cmd_ned = sumvec - gvec_LL;
    acc_cmd_vel = Tr_VN*acc_cmd_ned;

    speed_accel_cmd = acc_cmd_vel[0];
    eta_turn_cmd = acc_cmd_vel[1];
    eta_dive_cmd = acc_cmd_vel[2];

    turn_lim = 3.0*go;
    speed_accel_cmd = Bound(speed_accel_cmd, -turn_lim, turn_lim );

    vcmd_coop = vcmd_coop + dt50*speed_accel_cmd;
    eta_turn_cmd = Bound(eta_turn_cmd, -turn_lim, turn_lim );
    eta_dive_cmd = Bound(eta_dive_cmd, -turn_lim, turn_lim );

end

################################################################################
## Altitude Command
function AltitudeControl(nav_state::NavStates, alt_cmd::Float64)

    ## Parameter Computation
    # Computing NED to Body DCM
    Tr_BN = nav_state.euler.Tr_bn;
    Tr_NE = nav_state.wgs_pos.Tr_ne;
    rm_ecef = nav_state.wgs_pos.rm_ecef;
    vm_ecef = nav_state.velocity.vm_ecef;
    vm_ned = Tr_NE*vm_ecef;

    gamma = atan2( -vm_ned[3], sqrt( vm_ned[1]^2 + vm_ned[2]^2 ) );
    alt = nav_state.wgs_pos.altitude;
    Vmag = nav_state.velocity.vmag;

    ## Altitude Control Gains
    Kalt = 0.210;

    ## Altitude Error
    alt_err = alt_cmd - alt;
    println("alt_err = ", alt_err)

    ## Gamma Command
    gammac = Kalt/Vmag*alt_err;

    gamlim = 10.0*d2r;
    gammac = Bound(gammac, -gamlim, gamlim);

    ## Gamma Control Gains
    wn = 1.2;
    zz = 0.70;
    Kvf = 2*zz*wn;
    Kpf = wn*wn;

    ## Eta Dive Command
    eta_dive_cmd = -Kpf*Vmag*(gammac - gamma);
    eta_dive_cmd = Bound(eta_dive_cmd, -1.5*go, 1.5*go );

    return eta_dive_cmd;
end