text:
### Abstract variable IOStreams

abstract VariableIOStream

### BasicVariableIOStream

type BasicVariableIOStream <: VariableIOStream
  stream::IOStream
  size::Tuple
  n::Integer
end

BasicVariableIOStream(size::Tuple, n::Integer, filename::AbstractString, mode::AbstractString="w") =
  BasicVariableIOStream(open(filename, mode), size, n)

BasicVariableIOStream(
  size::Tuple,
  n::Integer;
  filepath::AbstractString="",
  filesuffix::AbstractString="csv",
  mode::AbstractString="w"
) =
  BasicVariableIOStream(size, n, joinpath(filepath, "value."*filesuffix), mode)

Base.close(iostream::BasicVariableIOStream) = close(iostream.stream)

Base.write(iostream::BasicVariableIOStream, state::VariableState) = write(iostream.stream, join(state.value, ','), "\n")

Base.write(iostream::BasicVariableIOStream, nstate::BasicUnvVariableNState) = writedlm(iostream.stream, nstate.value)

Base.write(iostream::BasicVariableIOStream, nstate::BasicMuvVariableNState) = writedlm(iostream.stream, nstate.value', ',')

function Base.write(iostream::BasicVariableIOStream, nstate::BasicMavVariableNState)
  statelen = prod(nstate.size)
  for i in 1:nstate.n
    write(iostream.stream, join(nstate.value[1+(i-1)*statelen:i*statelen], ','), "\n")
  end
end

Base.read!{N<:Number}(iostream::BasicVariableIOStream, nstate::BasicUnvVariableNState{N}) =
  nstate.value = vec(readdlm(iostream.stream, ',', N))

Base.read!{N<:Number}(iostream::BasicVariableIOStream, nstate::BasicMuvVariableNState{N}) =
  nstate.value = readdlm(iostream.stream, ',', N)'

function Base.read!{N<:Number}(iostream::BasicVariableIOStream, nstate::BasicMavVariableNState{N})
  statelen = prod(iostream.size)
  line = 1
  while !eof(iostream.stream)
    nstate.value[1+(line-1)*statelen:line*statelen] = [parse(N, c) for c in split(chomp(readline(iostream.stream)), ',')]
    line += 1
  end
end

function Base.read{N<:Number}(iostream::BasicVariableIOStream, T::Type{N})
  nstate::VariableNState
  l = length(iostream.size)

  if l == 0
    nstate = BasicUnvVariableNState(iostream.n, T)
  elseif l == 1
    nstate = BasicMuvVariableNState(iostream.size[1], iostream.n, T)
  elseif l == 2
    nstate = BasicMavVariableNState(iostream.size, iostream.n, T)
  else
    error("BasicVariableIOStream.size must be a tuple of length 0 or 1 or 2, got $(iostream.size) length")
  end

  read!(iostream, nstate)

  nstate
end

Base.show(io::IO, iostream::BasicVariableIOStream) =
  print(io, "BasicVariableIOStream: state size = $(iostream.size), number of states = $(iostream.n)")

Base.writemime(io::IO, ::MIME"text/plain", iostream::BasicVariableIOStream) = show(io, iostream)


text:
"""
    visualize(<arguments>; <keyword arguments>) -> Plot

Make a plot from an output collected by running necessary simulations. A convenient function to run both `simulate` and `plot` together.

See also: [`visualize!`](@ref), [`simulate`](@ref), [`plot`](@ref), [`manipulate`](@ref)

# Examples
```julia-repl
julia> @system S(Controller) begin
           a(a) => a ~ accumulate(init=1)
       end;

julia> visualize(S, :time, :a; stop=5, kind=:line)
       ┌────────────────────────────────────────┐
    32 │                                       :│
       │                                      : │
       │                                     :  │
       │                                    :   │
       │                                   :    │
       │                                  :     │
       │                                 :      │
  a    │                                :       │
       │                              .'        │
       │                            .'          │
       │                          .'            │
       │                       ..'              │
       │                   ..''                 │
       │             ....''                     │
     1 │.........''''                           │
       └────────────────────────────────────────┘
       0                                        5
                      time (hr)
```
"""
visualize(a...; kw...) = plot(a...; kw...)
"""
    visualize!(p, <arguments>; <keyword arguments>) -> Plot

Update an existing `Plot` object `p` by appending a new graph made with `visualize`.

See also: [`visualize`](@ref)

# Arguments
- `p::Union{Plot,Nothing}`: plot object to be updated; `nothing` creates a new plot.
"""
visualize!(a...; kw...) = plot!(a...; kw...)

visualize(S::Type{<:System}, x, y; kw...) = visualize!(nothing, S, x, y; kw...)
visualize!(p, S::Type{<:System}, x, y;
    config=(), group=(), xstep=(),
    stop=nothing, snap=nothing,
    ylab=nothing, legend=nothing, names=nothing, colors=nothing, plotopts...
) = begin
    G = configure(group)
    C = @config config + !G
    n = length(C)

    legend!(k, S) = begin
        isnothing(legend) && (legend = string(k))
        #TODO: support custom unit for legend?
        u = vartype(S, k) |> unittype
        legend = label(legend, u)
    end

    #HACK: support indirect referencing of label by variable name
    names = if names isa Symbol
        k = names
        legend!(k, S)
        repeat([k], n)
    elseif !isnothing(names)
        names
    elseif isempty(G)
        [""]
    elseif G isa Vector
        # temporary numeric labels
        string.(1:n)
    else
        K, V = only(G)
        k, v = only(V)
        T = K == Symbol(0) ? S : typefor(K, scopeof(S))
        legend!(k, T)
        string.(v)
    end
    isnothing(colors) && (colors = repeat([nothing], n))
    isnothing(ylab) && (ylab = y)

    s(c) = simulate(S; target=[x, y], configs=@config(c + !xstep), stop, snap, verbose=false)
    r = s(C[1])
    p = plot!(p, r, x, y; ylab, legend, name=names[1], color=colors[1], plotopts...)
    for i in 2:n
        r = s(C[i])
        p = plot!(p, r, x, y; name=names[i], color=colors[i], plotopts...)
    end
    p
end

visualize(SS::Vector, x, y; kw...) = visualize!(nothing, SS, x, y; kw...)
visualize!(p, SS::Vector, x, y; configs=[], names=nothing, colors=nothing, kw...) = begin
    n = length(SS)
    isempty(configs) && (configs = repeat([()], n))
    @assert length(configs) == n
    isnothing(names) && (names = string.(nameof.(SS)))
    isnothing(colors) && (colors = repeat([nothing], n))

    for (S, config, name, color) in zip(SS, configs, names, colors)
        p = visualize!(p, S, x, y; config, name, color, kw...)
    end
    p
end

visualize(S::Type{<:System}, x, y::Vector; kw...) = visualize!(nothing, S, x, y; kw...)
visualize!(p, S::Type{<:System}, x, y::Vector;
    config=(), xstep=(),
    stop=nothing, snap=nothing,
    plotopts...
) = begin
    r = simulate(S; target=[x, y], configs=@config(config + !xstep), stop, snap, verbose=false)
    plot!(p, r, x, y; plotopts...)
end

visualize(df::DataFrame, S::Type{<:System}, x, y; kw...) = visualize!(nothing, df, S, x, y; kw...)
visualize!(p, df::DataFrame, S::Type{<:System}, x, y; config=(), kw...) = visualize!(p, df, [S], x, y; configs=[config], kw...)
visualize(df::DataFrame, SS::Vector, x, y; kw...) = visualize!(nothing, df, SS, x, y; kw...)
visualize!(p, df::DataFrame, SS::Vector, x, y;
    configs=[], xstep=(),
    stop=nothing, snap=nothing,
    xlab=nothing, ylab=nothing, name=nothing, names=nothing, colors=nothing, xunit=nothing, yunit=nothing, plotopts...
) = begin
    x = x isa Pair ? x : x => x
    y = y isa Pair ? y : y => y
    xo, xe = x
    yo, ye = y
    xlab = isnothing(xlab) ? xe : xlab
    ylab = isnothing(ylab) ? ye : ylab

    u(n) = extractunit(df, n)
    isnothing(xunit) && (xunit = u(xo))
    isnothing(yunit) && (yunit = u(yo))

    n = length(SS)
    isempty(configs) && (configs = repeat([()], n))
    @assert length(configs) == n
    isnothing(names) && (names = string.(nameof.(SS)))
    isnothing(colors) && (colors = repeat([nothing], n))

    p = plot!(p, df, xo, yo; kind=:scatter, name, xlab, ylab, xunit, yunit, plotopts...)
    for (S, c, name, color) in zip(SS, configs, names, colors)
        cs = isnothing(xstep) ? c : @config(c + !xstep)
        r = simulate(S; target=[xe, ye], configs=cs, stop, snap, verbose=false)
        p = plot!(p, r, xe, ye; kind=:line, name, color, xunit, yunit, plotopts...)
    end
    p
end

visualize(obs::DataFrame, S::Type{<:System}, y::Vector; kw...) = visualize!(nothing, obs, S, y; kw...)
visualize!(p, obs::DataFrame, S::Type{<:System}, y::Vector;
    index,
    config=(), configs=[],
    stop=nothing, snap=nothing,
    names=nothing, plotopts...
) = begin
    #HACK: use copy due to normalize!
    obs = copy(obs)
    I = parseindex(index, S) |> keys |> collect
    T = parsetarget(y, S)
    Yo = T |> keys |> collect
    Ye = T |> values |> collect

    est = simulate(S; config, configs, index, target=Ye, stop, snap, verbose=false)
    normalize!(obs, est, on=I)
    df = DataFrames.innerjoin(obs, est, on=I)
    Xs = extractarray.(Ref(df), Yo)
    Ys = extractarray.(Ref(df), Ye)

    isnothing(names) && (names = [o == e ? "$o" : "$o ⇒ $e" for (o, e) in T])

    _visualize_obs_vs_est!(p, Xs, Ys; names, plotopts...)
end

visualize(obs::DataFrame, S::Type{<:System}, y; kw...) = visualize!(nothing, obs, S, y; kw...)
visualize!(p, obs::DataFrame, S::Type{<:System}, y; config=(), configs=[], name="", kw...) = visualize!(p, obs, [(; system=S, config, configs)], y; names=[name], kw...)
visualize(obs::DataFrame, maps::Vector{<:NamedTuple}, y; kw...) = visualize!(nothing, obs, maps, y; kw...)
visualize!(p, obs::DataFrame, maps::Vector{<:NamedTuple}, y;
    index,
    config=(), configs=[],
    stop=nothing, snap=nothing,
    names=nothing, plotopts...
) = begin
    #HACK: use copy due to normalize!
    obs = copy(obs)
    I = parseindex(index, obs) |> keys |> collect
    yo, ye = parsetarget(y, obs) |> collect |> only

    ests = map(m -> simulate(; m..., index, target=ye, stop, snap, verbose=false), maps)
    normalize!(obs, ests..., on=I)
    dfs = map(est -> DataFrames.innerjoin(obs, est, on=I), ests)
    Xs = extractarray.(dfs, yo)
    Ys = extractarray.(dfs, ye)

    isnothing(names) && (names = ["$(m.system)" for m in maps])

    _visualize_obs_vs_est!(p, Xs, Ys; names, plotopts...)
end

_visualize_obs_vs_est!(p, Xs, Ys; xlab=nothing, ylab=nothing, names=nothing, colors=nothing, lim=nothing, plotopts...) = begin
    isnothing(xlab) && (xlab = "Observation")
    isnothing(ylab) && (ylab = "Model")

    n = length(Ys)
    isnothing(names) && (names = repeat([""], n))
    isnothing(colors) && (colors = repeat([nothing], n))

    if isnothing(lim)
        l = findlim.(deunitfy.([Xs..., Ys...]))
        lim = (minimum(minimum.(l)), maximum(maximum.(l)))
    end
    L = [lim[1], lim[2]]

    for (X, Y, name) in zip(Xs, Ys, names)
        p = plot!(p, X, Y; kind=:scatter, name, xlab, ylab, xlim=lim, ylim=lim, aspect=1, plotopts...)
    end
    !isnothing(lim) && plot!(p, L, L, kind=:line, color=:lightgray, name="")
    p
end

visualize(S::Type{<:System}, x, y, z;
    config=(), xstep=(), ystep=(),
    stop=nothing, snap=nothing,
    plotopts...
) = begin
    configs = @config config + xstep * ystep
    r = simulate(S; index=[x, y], target=z, configs, stop, snap, verbose=false)
    plot(r, x, y, z; plotopts...)
end

export visualize, visualize!

ext: jl


text:
## convert.jl Convenience function for type conversion
## (c) 2014 David A. van Leeuwen

## just give me the array!

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

import Base: convert

## from array to namedarray, the fallback constructor
convert(::Type{NamedArray}, a::AbstractArray) = NamedArray(a)
convert(::Type{NamedVector}, a::AbstractArray) = NamedArray(a)

## to array
## array(n::NamedArray) = n.array
## convert{T,N}(::Type{AbstractArray{T,N}}, a::NamedArray{T,N}) = a.array
## convert{T}(::Type{AbstractArray{T}}, a::NamedArray{T}) = a.array
## convert{T,N}(::Type{AbstractArray}, a::NamedArray{T,N}) = a.array

## to other type
convert(::Type{NamedArray{T}}, n::NamedArray) where {T} = NamedArray(T.(n.array), n.dicts, n.dimnames)

@inline convert(::Type{NamedArray}, n::NamedArray) = n


text:
using Yao, Zygote, Flux.Optimise
include("patch.jl")

# make a one qubit GHZ state as learning target
# NOTE: this can be an arbitrary one qubit state
t = ArrayReg(bit"0") + ArrayReg(bit"1")
normalize!(t)

# calculate the fidelity with GHZ state
function fid(xs)
    r = zero_state(1)
    U = mat(chain(Rx(xs[1]), Rz(xs[2])))
    return abs(statevec(t)' * U * statevec(r))
end

# simply tell Zygote to get the gradient and start training
function train()
    opt = ADAM()
    xs = rand(2)
    for _ in 1:1000
        println(fid(xs))
        Optimise.update!(opt, xs, -fid'(xs))
    end
    return xs
end

train()

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
module TestLazyArrays

include("preamble.jl")
using LazyArrays: Vcat, Hcat

@testset "Vcat" begin
    @testset for xs in @expressions begin
            Vcat(1:5, 2:3)
            Vcat(zeros(1, 2), ones(2, 2))
        end
        xs = @eval $xs
        @test collect(Map(identity), xs) == collect(xs)[:]
    end

end

@testset "Hcat" begin
    @testset for xs in @expressions begin
            Hcat(1:2, 3:4)
            Hcat(zeros(2, 1), ones(2, 2))
        end
        xs = @eval $xs
        @test collect(Map(identity), xs) == collect(xs)[:]
    end

end

end  # module


text:
include("MatrixMarket.jl")
using MatrixMarket
using LinearLeastSquares

# read in the data
A = full(MatrixMarket.mmread("largeCorpus.mtx"));

# extract the classes of each document
classes = A[:,1];
# TODO: modify classes so that 4 5 6 are 1 2 3
classes[classes .> 3] = classes[classes .> 3] - 3;
A = A[:, 2:end];

# split into train/test
numData = size(A, 1);
data = randperm(numData);
ind = floor(numData*0.7);
training = data[1:ind];
test = data[ind+1:end];
trainDocuments = A[training,:];
trainClasses = classes[training,:];
testDocuments = A[test,:];
testClasses = classes[test,:];

# change all other than sports to -1 (sports is 1)
holdClass = 1;
trainClasses[trainClasses .!= holdClass] = -1;
trainClasses[trainClasses .== holdClass] = 1;
testClasses[testClasses .!= holdClass] = -1;
testClasses[testClasses .== holdClass] = 1;

# build the problem and solve with LLS
lambda = 100;
w = Variable(size(A, 2));
v = Variable();
objective = sum_squares(trainDocuments * w + v - trainClasses) + lambda * sum_squares(w);
optval = minimize!(objective);

# print out the 5 words most indicative of sports and nonsports
words = String[];
f = open("largeCorpusfeatures.txt");
for i = 1:length(evaluate(w))
  push!(words, readline(f))
end
indices = sortperm(vec(evaluate(w)));
for i = 1:5
  print(words[indices[i]])
end
for i = 0:4
  print(words[indices[length(words) - i]])
end

# calculate training error
yhat = sign(trainDocuments * evaluate(w) + evaluate(v));
trainCE =  1/size(trainClasses,1)*sum(trainClasses .!= yhat)

# calculate performance of our classifier on the test set
yhat2 = sign(testDocuments * evaluate(w) + evaluate(v));
testCE = 1/size(testClasses,1)*sum(testClasses .!= yhat2)


text:
# # Transfer Learning with Flux

# This article is intended to be a general guide to how transfer learning works in the Flux ecosystem.
# We assume a certain familiarity of the reader with the concept of transfer learning. Having said that,
# we will start off with a basic definition of the setup and what we are trying to achieve. There are many
# resources online that go in depth as to why transfer learning is an effective tool to solve many ML
# problems, and we recommend checking some of those out.

# Machine Learning today has evolved to use many highly trained models in a general task,
# where they are tuned to perform especially well on a subset of the problem.

# This is one of the key ways in which larger (or smaller) models are used in practice. They are trained on
# a general problem, achieving good results on the test set, and then subsequently tuned on specialised datasets.

# In this process, our model is already pretty well trained on the problem, so we don't need to train it
# all over again as if from scratch. In fact, as it so happens, we don't need to do that at all! We only
# need to tune the last couple of layers to get the most performance from our models. The exact last number of layers
# is dependant on the problem setup and the expected outcome, but a common tip is to train the last few `Dense`
# layers in a more complicated model.

# So let's try to simulate the problem in Flux.

# We'll tune a pretrained ResNet from Metalhead as a proxy. We will tune the `Dense` layers in there on a new set of images.

using Flux, Metalhead
using Flux: @epochs
resnet = ResNet().layers

# If we intended to add a new class of objects in there, we need only `reshape` the output from the previous layers accordingly.
# Our model would look something like so:

# ```julia
# model = Chain(
#   resnet[1:end-2],               # We only need to pull out the dense layer in here
#   x -> reshape(x, size_we_want), # / global_avg_pooling layer
#   Dense(reshaped_input_features, n_classes)
# )
# ```

# We will use the [Dogs vs. Cats](https://www.kaggle.com/c/dogs-vs-cats/data) dataset from Kaggle for our use here.
# Make sure to extract the images in a `train` folder.

# The `datatloader.jl` script contains some functions that will help us load batches of images, shuffled between
# dogs and cats along with their correct labels.

include("dataloader.jl")

# Finally, the model looks something like:

model = Chain(
  resnet[1:end-2],
  Dense(2048, 1000),
  Dense(1000, 256),
  Dense(256, 2),        # we get 2048 features out, and we have 2 classes
)

# To speed up training, let’s move everything over to the GPU

model = model |> gpu
dataset = [gpu.(load_batch(10)) for i in 1:10]

# After this, we only need to define the other parts of the training pipeline like we usually do.

opt = ADAM()
loss(x,y) = Flux.Losses.logitcrossentropy(model(x), y)

# Now to train
# As discussed earlier, we don’t need to pass all the parameters to our training loop. Only the ones we need to
# fine-tune. Note that we could have picked and chosen the layers we want to train individually as well, but this
# is sufficient for our use as of now.

ps = Flux.params(model[2:end])  # ignore the already trained layers of the ResNet

# And now, let's train!

@epochs 2 Flux.train!(loss, ps, dataset, opt)

# And there you have it, a pretrained model, fine tuned to tell the the dogs from the cats.

# We can verify this too.

imgs, labels = gpu.(load_batch(10))
display(model(imgs))

labels

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
module L_
    using MultiInterface

    abstract type SFunc{I,O} end
    struct MulVecs{U<:AbstractFloat,V<:AbstractFloat,W<:AbstractFloat} <: SFunc{Tuple{Vector{U}, Vector{V}}, Vector{W}, Nothing} end

    @interface sample(sf::SFunc{I,O}, i::I)::O where {I,O}

    @impl begin
        struct SampleMulVecs end
        function sample(sf::MulVecs{U,V,W}, x::Tuple{Vector{U}, Vector{V}})::Vector{W} where {U,V,W}
            return x[1] .* x[2]
        end
    end
end


text:
"""
   IO functions for the NEID spectrograph
   https://neid.psu.edu/
Author: Eric Ford and collaborators
Created: August 2020
"""

"""Create Dataframe containing filenames and key data for all files neid*.fits in directory"""
function make_manifest(data_path::String ; max_spectra_to_use::Int = 1000 )
    df_filenames = EchelleInstruments.make_manifest(data_path,r"^neidL1_\d+[T\.]\d+\.fits$")
    #=
    dir_filelist = readdir(data_path,join=true)
    idx_spectra = map(fn->occursin(r"^neidL1_\d+[T\.]\d+\.fits$", last(split(fn,'/')) ),dir_filelist)
    spectra_filelist = dir_filelist[idx_spectra]
    =#
    #=
    df_files = DataFrame(Filename = String[], target = String[], bjd = Float64[], ssbz=Float64[] )
    map(fn->add_metadata_from_fits!(df_files,fn),spectra_filelist)
    df_files
    =#
    #@assert length(spectra_filelist) >= 1
    #df_files = DataFrame(read_metadata(spectra_filelist[1]))
    df_files = DataFrame(read_metadata(df_filenames.Filename[1]))
    keys = propertynames(df_files)
    allowmissing!(df_files, keys[map(k->k∉[:Filename, :bjd, :target, :airmass],keys)] )
    #println("df_files = ", df_files)

    #if length(spectra_filelist) >= 2
    if length(df_filenames.Filename) >= 2
        map(fn->add_metadata_from_fits!(df_files,fn),df_filenames.Filename[2:end])
    end
    #=
    for i in 2:length(spectra_filelist)
        try
            add_metadata_from_fits!(df_files,spectra_filelist[i])

        catch
            println("# Problem reading ", spectra_filelist[i])
            fields_to_save = metadata_symbols_default(NEID2D())
            fields_str_to_save = metadata_strings_default(NEID2D())

            metadata_tmp = read_metadata_from_fits(spectra_filelist[i],fields=fields_to_save,fields_str=fields_str_to_save)
            println("Before: ", metadata_tmp)
            #println("After : ", metadata_tmp)
            #push!(df_files,metadata_tmp)
        end
    end
    =#
    df_files

end

"""Create Dict containing filename and default metadata from file."""
function read_metadata(fn::Union{String,FITS})
    fields_to_save = metadata_symbols_default(NEID2D())
    fields_str_to_save = metadata_strings_default(NEID2D())
    dict = read_metadata_from_fits(fn,fields=fields_to_save,fields_str=fields_str_to_save)
    check_metadata_fields_expected_present(dict,fields_to_save)
    dict[:Filename] = fn
    if haskey(dict,:airmass) && typeof(dict[:airmass]) == String
       dict[:airmass] = parse(Float64,dict[:airmass])
    end
    return dict
end

function add_metadata_from_fits!(df::DataFrame, fn::String)
    println("# Reading metadata from ", fn)
    metadata_keep = read_metadata(fn)
    #if metadata_keep[:target] != "Solar" return df end
    calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]

   if any(map(ss->contains( metadata_keep[:target],ss),calibration_target_substrings))
      return df
   end

    for k in propertynames(df)
        if typeof(metadata_keep[k]) == Missing continue end
        if typeof(metadata_keep[k]) == Nothing
            metadata_keep[k] = missing
            continue
        end
        if (eltype(df[!,k]) <: Union{Real,Union{<:Real,Missing}}) && !(typeof(metadata_keep[k]) <: Union{Real,Union{<:Real,Missing}})
            metadata_keep[k] = missing
        end
    end
    #println("metadata_keep = ", metadata_keep)
    if haskey(metadata_keep,:airmass) && typeof(metadata_keep[:airmass]) == String
       metadata_keep[:airmass] = parse(Float64,metadata_keep[:airmass])
    end
    push!(df, metadata_keep)
    return df
end

""" Read NEID data from FITS file, and return in a Spectra2DBasic object."""
function read_data   end

function read_data(f::FITS, metadata::Dict{Symbol,Any} )
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["SCIWAVE"]]), FITSIO.read(f[img_idx["SCIFLUX"]]), FITSIO.read(f[img_idx["SCIVAR"]])
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR ) where  { AR<:AbstractRange }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["SCIWAVE"]],:,orders_to_read), FITSIO.read(f[img_idx["SCIFLUX"]],:,orders_to_read), FITSIO.read(f[img_idx["SCIVAR"]],:,orders_to_read)
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    read_data(f, metadata)
end

function read_data(fn::String)
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_data(f, metadata)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index})
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata)
end

function read_data(fn::String, orders_to_read::AR ) where  { AR<:AbstractRange }
    f = FITS(fn)
    metadata = read_metadata(f)
    read_data(f,metadata, orders_to_read)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index}, orders_to_read::AR ) where  { AR<:AbstractRange }
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata, orders_to_read)
end

function apply_barycentric_correction!(λ::AbstractArray{T,1}, z::Real) where { T<:Real }

end

""" Read space delimited file with differential extinction corrections, interpolate to bjd's in df and insert into df[:,diff_ext_rv]. """
function read_differential_extinctions!(fn::String, df::DataFrame, df_time_col::Symbol = :bjd)
    df_diff_ext = CSV.read(fn, DataFrame, delim='\t')
    @assert any(isequal(:JD),propertynames(df_diff_ext))
    @assert any(isequal(:delta_vr),propertynames(df_diff_ext))
    diff_ext = LinearInterpolation(df_diff_ext.JD, df_diff_ext.delta_vr)
    df[!,:diff_ext_rv] = -diff_ext(df[!,df_time_col])/100 # cm/s -> m/s
    return df
end

""" Read normalization function for approximate solar continuum+blaze correction
TODO: Figure out mapping between these orders and what's stored in the FITS file before using. """
function read_solar_normalization(fn::String = joinpath(pkgdir(EchelleInstruments),"data","neid","NEID_normalization_function.csv") )
    CSV.read(fn, DataFrame)
end

""" Read calibrator data """
function read_cal_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR )  where { AR<:AbstractRange  }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: FITSIO.ImageHDU,2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["CALWAVE"]],:,orders_to_read), FITSIO.read(f[img_idx["CALFLUX"]],:,orders_to_read), FITSIO.read(f[img_idx["CALVAR"]],:,orders_to_read)
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_cal_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    read_cal_data(f, metadata)
end

function read_cal_data(fn::String)
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_cal_data(f, metadata)
end


text:
## im2col

a = [1 4
     2 5
     3 6]
b = [7 10
     8 11
     9 12]
x = [a b]
cols = im2col(x,size(a),"distinct")
@test cols == [vec(a) vec(b)]
@test cols == im2col!(similar(cols),x,size(a),size(a))
@test_throws DimensionMismatch im2col!(Array(Int,1,2),x,size(a),size(a))

x = randn(4,1,7)
sz = (2,1,3)
f = rand(sz)
@test_approx_eq  vec(f)' * im2col(x,sz,"sliding")  imfilter(x,f,"inner")

x = reshape(1:(3*7),3,7)
@test x[:] == im2col(x,(1,1))[:]
@test im2col(x,(3,1),(1,3)) == x[:,1:3:end]
@test_throws DimensionMismatch im2col(x,(3,1),(1,4))
@test_throws DimensionMismatch im2col(x,(4,1))
@test_throws ArgumentError im2col(x,(1,1),"foo")

## col2im

szI = (4,3,2)
szW  = (2,1,2)
N = length(szI)
onetuple = ntuple(d->1,N)
I = randn(szI)

# identity
@test I == col2im(reshape(I,1,length(I)),onetuple,szI,"distinct")
@test I == col2im(reshape(I,1,length(I)),onetuple,szI,"sliding")
@test I == col2im(reshape(I,1,length(I)),onetuple,szI,onetuple)

# distinct
cols = im2col(I,szW,"distinct")
@test I == col2im(cols,szW,szI,"distinct")
@test I == col2im(cols,szW,szI,szW)
@test I == col2im!(zeros(I),cols,szW,szW)

# sliding with averaging
cols = im2col(I,szW,"sliding")
@test I == col2im(cols,szW,szI)
@test I == col2im(cols,szW,szI,"sliding")
@test I == col2im(cols,szW,szI,"sliding","average")
@test I == col2im(cols,szW,szI,onetuple)
@test I == col2im(cols,szW,szI,ntuple(d->1,N),+,0,/)

# sliding
I = rand(3,3)
cols = im2col(I,(2,2),"sliding")
w = [1 2 1; 2 4 2; 1 2 1]
ref = I .* w
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),"sliding","sum")
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),(1,1),+,0,nothing)
@test_approx_eq  ref  col2im!(zeros(I),cols,(2,2),(1,1),+)
ref = I .^ w
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),(1,1),*,1,nothing)
@test_approx_eq  ref  col2im!(ones(I),cols,(2,2),(1,1),*)
ref = I
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),"sliding","average")
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),(1,1),+,0,/)
ref = I .* w - w
@test_approx_eq  ref  col2im(cols,(2,2),(3,3),(1,1),+,0,-)

# errors
I = rand(Int,4,4); cols = im2col(I,(2,2),"sliding")
@test_throws ArgumentError col2im(cols,(2,2),(4,4),"sliding","average")
I = rand(4,4); cols = im2col(I,(2,2),"sliding")
@test_throws DimensionMismatch col2im(cols,(2,2),(3,4))
@test_throws DimensionMismatch col2im(cols,(2,2),(4,4),"distinct")
@test_throws ArgumentError col2im(cols,(2,2),(4,4),"sliding","foo")
@test_throws ArgumentError col2im(cols,(2,2),(4,4),"foo")


text:
#######################################################################
# KetModel
#######################################################################
struct KetModel{O<:Val,M<:TightbindingModel,R<:Union{Missing,Real}}
    model::M
    normalization::R
    maporbitals::O       # Val(false) or Val(true) to aid in type-stability
    singlesitekets::Bool
end
"""
    ketmodel(a; region = missing, sublats = missing, normalization = 1, maporbitals = false, singlesitekets = false)

Create an `KetModel` of amplitude `a` on any site in the specified `region` and `sublats`.
For single-column kets, the amplitude `a` can be a `Number`, an `AbstractVector`, or for a
position-dependent amplitude a function of the form `r -> ...` returning either. For
multi-column kets, make `a` an `AbstractMatrix{<:Number}` or a function returning one, which
will be sliced into each ket column as appropriate. An error will be thrown if the slicing
is impossible, due e.g. to a mismatch of `size(a, 1)` and the number of orbitals in an
applicable sublattice.

# Keyword arguments

If keyword `normalization` is not `missing` or `false`, each column of the ket will be
rescaled to have norm `normalization` when the `KetModel` is applied to a specific
Hamiltonian. If a ket column `iszero`, however, it will not be normalized.

If keyword `maporbitals = true` and amplitude `a` is a scalar or a scalar function, `a`
will be applied to each orbital independently. This is particularly useful in multiorbital
systems with random amplitudes, e.g. `a = r -> randn()`. If `a` is not a scalar and
`maporbitals == true`, an error will be thrown.

If keyword `singlesitekets = true`, then the model represents a multicolumn ket, where each
column (or block of columns for `a::AbstractMatrix`) has amplitude `a` on a single site of
those selected by `region` and `sublats` (as opposed to having the same amplitude `a` on all
said sites if `singlesitekets = false`). This is useful e.g. to build a basis for the selected
sites.

Keywords `region` and `sublats` are the same as for `siteselector`. Only sites at position
`r` in sublattice with name `s::NameType` will be selected if `region(r) && s in sublats` is
true. Any missing `region` or `sublat` will not be used to constraint the selection.

The keyword `sublats` allows the following formats:

    sublats = :A           # Onsite on sublat :A only
    sublats = (:A,)        # Same as above
    sublats = (:A, :B)     # Onsite on sublat :A and :B

# Ket algebra

`KetModel`s created with `ket` can added or substracted together or be multiplied by scalars
to build more elaborate `KetModel`s, e.g. `ket(1) - 3 * ket(2, region = r -> norm(r) < 10)`.
Only models with the same `maporbitals` can be combined. When combining two models with
different `singlesitekets`, the result has `singlesitekets = true`.

# Examples

```jldoctest
julia> k = ketmodel(1, sublats=:A) - ketmodel(1, sublats=:B)
KetModel{2}: model with 2 terms
  Normalization : 1
  Map orbitals  : Val{false}()
  OnsiteTerm{Int64}:
    Sublattices      : A
    Coefficient      : 1
  OnsiteTerm{Int64}:
    Sublattices      : B
    Coefficient      : -1
```
# See also
    `ket`, `onsite`, `orbitalstructure`

"""
ketmodel(f; normalization = 1, maporbitals::Bool = false, singlesitekets::Bool = false, kw...) =
    KetModel(onsite(f; kw...), sanitize_normalization(normalization), Val(maporbitals), singlesitekets)

maporbitals(m::KetModel{Val{true}}) = true
maporbitals(m::KetModel{Val{false}}) = false

sanitize_normalization(::Missing) = missing
sanitize_normalization(b::Bool) = b ? 1 : missing
sanitize_normalization(b) = b

function Base.show(io::IO, ::MIME"text/plain", k::KetModel{<:Any,M}) where {N,M<:TightbindingModel{N}}
    i = get(io, :indent, "")
    ioindent = IOContext(io, :indent => "$i  ")
    print(io, "$(i)KetModel{$N}: model with $N terms
$i  Normalization    : $(k.normalization)
$i  Map orbitals     : $(k.maporbitals)
$i  Single-site kets : $(k.singlesitekets)")
    foreach(t -> print(ioindent, "\n", t), k.model.terms)
end

Base.:*(x::Number, k::KetModel) =
    KetModel(k.model * x, k.normalization, k.maporbitals, k.singlesitekets)
Base.:*(k::KetModel, x::Number) =
    KetModel(x * k.model, k.normalization, k.maporbitals, k.singlesitekets)
Base.:-(k::KetModel) =
    KetModel(-k.model, k.normalization, k.maporbitals, k.singlesitekets)
Base.:-(k1::KetModel, k2::KetModel) = k1 + (-k2)

function Base.:+(k1::KetModel, k2::KetModel)
    newnormalization = sanitize_norm(k1.normalization, k2.normalization)
    if k1.maporbitals == k2.maporbitals
        return KetModel(k1.model + k2.model, newnormalization,
            k1.maporbitals, k1.singlesitekets || k2.singlesitekets)
    else
        throw(ArgumentError("Cannot combine ket models with different `maporbitals`"))
    end
end

sanitize_norm(n1::Number, n2::Number) = n1 ≈ n2 ? n1 : _normwarn()
sanitize_norm(::Missing, ::Missing) = missing
sanitize_norm(n1, n2) = _normwarn()

function _normwarn()
    @warn "Combining `KetModel`s with different normalizations, choosing `normalization = missing`"
    return missing
end

resolve(k::KetModel, lat::AbstractLattice) =
    KetModel(resolve(k.model, lat), k.normalization, k.maporbitals, k.singlesitekets)

#######################################################################
# Ket
#######################################################################
# The eltype T of any ket must be equal to orbitaltype(orbstruct), i.e. Number or SVector
struct Ket{T,M<:AbstractMatrix{T},O<:OrbitalStructure{T}} <: AbstractMatrix{T}
    amplitudes::M
    orbstruct::O
    function Ket{T,M,O}(amplitudes, orbstruct) where {T,M<:AbstractMatrix{T},O<:OrbitalStructure{T}}
        check_compatible_ket(amplitudes, orbstruct)
        new(amplitudes, orbstruct)
    end
end

Ket(amplitudes::M, orbstruct::O) where {T,M<:AbstractMatrix{T},O<:OrbitalStructure{T}} =
    Ket{T,M,O}(amplitudes, orbstruct)

Ket(amplitudes::M, orbstruct::O) where {T,M<:AbstractMatrix,O<:OrbitalStructure{T}} =
    Ket{T,M,O}(T.(amplitudes), orbstruct)

check_compatible_ket(a::AbstractMatrix{T}, o::OrbitalStructure) where {T} =
    T == orbitaltype(o) && size(a, 1) == dimh(o) ||
        throw(ArgumentError("Ket is incompatible with OrbitalStructure"))

check_compatible_ket(ket::Ket, h::Hamiltonian) =
    orbitalstructure(ket) == orbitalstructure(h) ||
        throw(ArgumentError("Ket is incompatible with Hamiltonian"))

function Base.show(io::IO, ::MIME"text/plain", k::Ket{T}) where {T}
    ioindent = IOContext(io, :indent => "  ")
    print(io, "Ket{$T}: ket with a $(size(k.amplitudes, 1)) × $(size(k.amplitudes, 2)) amplitude matrix\n")
  show(ioindent, k.orbstruct)
end

orbitalstructure(k::Ket) = k.orbstruct

# AbstractArray interface
Base.getindex(k::Ket, args...) = getindex(k.amplitudes, args...)
Base.setindex!(k::Ket, args...) = setindex!(k.amplitudes, args...)
Base.firstindex(k::Ket) = firstindex(k.amplitudes)
Base.lastindex(k::Ket) = lastindex(k.amplitudes)
Base.length(k::Ket) = length(k.amplitudes)
Base.size(k::Ket, args...) = size(k.amplitudes)
Base.IndexStyle(k::Ket) = IndexStyle(k.amplitudes)
Base.similar(k::Ket) = Ket(similar(k.amplitudes), k.orbstruct)
Base.parent(k::Ket) = k.amplitudes

"""
    ket(m::AbstractArray, o::OrbitalStructure)
    ket(m::AbstractArray, h::Hamiltonian)

Construct a `Ket` `|k⟩` with amplitudes `⟨i|k⟩ = m[i]`, which can be scalars or `SVector`s
depending on the number of orbitals on site `i`. If `m` is an `AbstractMatrix` instead of an
`AbstractVector`, the `Ket` represents a multi-column ket (i.e. a collection of kets `|kⱼ⟩`,
one per column), such that `⟨i|kⱼ⟩ = m[i,j]`. The orbitals per sublattice are encoded in `o
= orbitalstructure(h)`.

    ket(km::KetModel, h::Hamiltonian)

Construct a `Ket` by applying model `km` to Hamiltonian `h` (see also `ketmodel` for
details).

    ket(h::Hamiltonian)

Construct a zero ket, equivalent to `ket(0, h; maporbitals = true, normalization = missing)`

# See also
    `ketmodel`, `ket!`, `onsite`
"""
ket(m::AbstractVector, o::OrbitalStructure) = Ket(hcat(m), o)
ket(m::AbstractMatrix, o::OrbitalStructure) = Ket(m, o)
ket(m::AbstractArray, h::Hamiltonian) = ket(m, orbitalstructure(h))
ket(h::Hamiltonian) = ket(ketmodel(0; maporbitals = true, normalization = missing), h)

#######################################################################
# flatten and unflatten
#######################################################################

flatten(k::Ket) = ket(flatten(parent(k), orbitalstructure(k)), flatten(orbitalstructure(k)))

unflatten(k::Ket, o::OrbitalStructure) = Ket(unflatten_orbitals(parent(k), o), o)

#######################################################################
# ket(::KetModel, ::Hamiltonian)
#######################################################################
function ket(kmodel::KetModel, h::Hamiltonian)
    kmodel´ = resolve(kmodel, h.lattice)
    ncols = guess_ket_columns(kmodel´, h)
    T = orbitaltype(h)
    kmat = Matrix{T}(undef, size(h, 1), ncols)
    k = ket(kmat, h)
    return ket!(k, kmodel´, h)
end

function guess_ket_columns(km, h)
    cols = 0
    r = first(allsitepositions(h.lattice))
    if km.singlesitekets
        cols = 0
        for term in km.model.terms
            t = term(r,r)
            for (sublat, orbs) in enumerate(orbitals(h)), i in siterange(h.lattice, sublat)
                if i in term.selector
                    cols += _term_columns(t, orbs)
                end
            end
        end
    else
        # This guess assumes that all model terms have the same size(amplitude, 2)
        term = first(km.model.terms)
        cols = _term_columns(term(r, r))
    end
    return cols
end

_term_columns(::Number, orbs...) = 1
_term_columns(::AbstractVector, orbs...) = 1
_term_columns(t::AbstractMatrix, orbs...) = size(t, 2)
_term_columns(::UniformScaling, orbs) = length(orbs)
_term_columns(x, orbs...) = throw(ArgumentError("Ket model amplitude should be of type `T` or a function `a(r)::T`, where `T is a `Number`, an `AbstractVector`, an `AbstractMatrix` or, in the case of `singlesitekets = true`, a `Uniformscaling`"))

# Model application, possibly flattening if target requires it (like bloch!)
# The type instability in ket! (due to orbs in multi-orbital h's) is harmless
function ket!(k::Ket, kmodel, h)
    kmat = parent(k)
    T = eltype(kmat)
    fill!(kmat, zero(T))
    kmodel´ = resolve(kmodel, h.lattice)    # resolve sublat names into sublat indices
    coloffset = Ref(0)                      # column counter for singlesitekets = true
    for (sublat, orbs) in enumerate(orbitals(h))
        ket_applyterms_sublat!(kmat, sublat, orbs, h, kmodel´, coloffset)
    end
    @assert !kmodel.singlesitekets || coloffset[] == size(kmat, 2) "Internal error, bug in ket!"
    kmodel.normalization === missing || normalize_columns!(kmat, kmodel.normalization)
    return ket(kmat, k.orbstruct)
end

# function barrier for orbs type-stability
function ket_applyterms_sublat!(kmat::AbstractArray{T}, sublat, orbs, h, kmodel, coloffset) where {T}
    allpos = allsitepositions(h.lattice)
    orbstruct = orbitalstructure(h)
    for i in siterange(h.lattice, sublat), term in kmodel.model.terms
        if i in term.selector
            r = allpos[i]
            orbsmat = ket_orbs_matrix(kmodel.maporbitals, term, r, orbs)
            di = size(orbsmat, 2)
            if kmodel.singlesitekets
                kmatview = view(kmat, :, coloffset[]+1:coloffset[]+di)
                coloffset[] += di
            else
                di == size(kmat, 2) ||
                    throw(ArgumentError("All ket model terms should have amplitudes with the same number of columns"))
                kmatview = view(kmat, :, 1:di)  # for type stability in branch
            end
            copy_rows!(kmatview, i, orbstruct, orbsmat, orbs)
        end
    end
    return kmat
end

# should return only orbs matrix, so that to_orbtype pads with zeros correctly
# must evaluate term once per orbital, in case it is a random function
ket_orbs_matrix(::Val{true}, term, r, orbs) = rows_to_matrix((orb -> ensure_singlerow(term(r, r))).(orbs))
ket_orbs_matrix(::Val{false}, term, r, orbs) = ensure_orbs_matrix(term(r, r), orbs)

ensure_singlerow(row::Number) = row
ensure_singlerow(row) = size(row, 1) == 1 ? row :
    throw(ArgumentError("Expected a scalar or a single-row matrix in ket model with `maporbitals = true`, got $(size(row, 1)) rows"))

rows_to_matrix(rows::NTuple{<:Any,Number}) = hcat(vcat(rows...))
rows_to_matrix(rows::NTuple{<:Any,AbstractVector}) = hcat(vcat(rows...))
rows_to_matrix(rows) = vcat(rows...)

ensure_orbs_matrix(x::Number, ::NTuple{N}) where {N} = N == 1 ? SMatrix{1,1}(x) :
    throw(ArgumentError("Expected an array with $N rows in ket model with `maporbitals = false`, got a scalar"))
ensure_orbs_matrix(v::AbstractVector, ::NTuple{N}) where {N} = length(v) == N ? hcat(v) :
    throw(ArgumentError("Expected an array with $N rows in ket model with `maporbitals = false`, got a vector of length $(length(v))"))
ensure_orbs_matrix(mat::AbstractMatrix, ::NTuple{N}) where {N} = size(mat, 1) == N ? mat :
    throw(ArgumentError("Expected an array with $N rows in ket model with `maporbitals = false`, got $(size(mat, 1)) rows"))
ensure_orbs_matrix(mat::UniformScaling, ::NTuple{N}) where {N} = SMatrix{N,N}(mat)

# copy SVectors
function copy_rows!(kmat::AbstractMatrix{T}, i, orbstruct, orbsmat, orbs) where {T<:SVector}
    kmat[i, :] .+= to_orbtype.(eachcol(orbsmat), T, Ref(orbs))
    return kmat
end

# copy Scalars
function copy_rows!(kmat::AbstractMatrix{T}, i, orbstruct, orbsmat, orbs) where {T<:Number}
    row = flatoffset_site(i, orbstruct) + 1
    dr = size(orbsmat, 1)
    kmat[row:row+dr-1, :] .+= T.(orbsmat)
    return kmat
end

to_orbtype(t::Number, ::Type{S}, t1::NTuple{1}) where {S<:SVector} = padtotype(t, S)
to_orbtype(t::AbstractVector, ::Type{S}, t1::NTuple{N}) where {N,S} = padtotype(SVector{N}(t), S)

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using Documenter
using ResumableFunctions
using SimJulia

makedocs(
  sitename = "SimJulia",
  authors = "Ben Lauwens",
  pages    = [
    "Home" => "index.md",
    "Tutorial" => "tutorial.md",
    "Topical Guides" => ["Basics" => "guides/basics.md",
                         "Environments" => "guides/environments.md",
                         "Events" => "guides/events.md",],
    "Examples" => ["Ross" => "examples/ross.md", "Latency" =>
                   "examples/Latency.md"],
    "API" => "api.md"
  ]
)

deploydocs(
  repo = "github.com/BenLauwens/SimJulia.jl.git"
)


text:
# ---
# title: 424. Longest Repeating Character Replacement
# id: problem424
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/longest-repeating-character-replacement/description/>
# hidden: true
# ---
#
# Given a string `s` that consists of only uppercase English letters, you can
# perform at most `k` operations on that string.
#
# In one operation, you can choose **any** character of the string and change it
# to any other uppercase English character.
#
# Find the length of the longest sub-string containing all repeating letters you
# can get after performing the above operations.
#
# **Note:**
# Both the string's length and _k_ will not exceed 104.
#
# **Example 1:**
#
#
#
#     Input:
#     s = "ABAB", k = 2
#
#     Output:
#     4
#
#     Explanation:
#     Replace the two 'A's with two 'B's or vice versa.
#
#
#
#
# **Example 2:**
#
#
#
#     Input:
#     s = "AABABBA", k = 1
#
#     Output:
#     4
#
#     Explanation:
#     Replace the one 'A' in the middle with 'B' and form "AABBBBA".
#     The substring "BBBB" has the longest repeating letters, which is 4.
#
#
#
#
#
## @lc code=start
using LeetCode

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


text:
# wrapers for old functionality

## removed in CUDA 11

@checked function cusparseScsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_api()
    @runtime_ccall((:cusparseScsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{Cfloat}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{Cfloat}, CuPtr{Cint}, CuPtr{Cint}, cusparseAction_t,
                    cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseDcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_api()
    @runtime_ccall((:cusparseDcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{Cdouble}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{Cdouble}, CuPtr{Cint}, CuPtr{Cint},
                    cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseCcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_api()
    @runtime_ccall((:cusparseCcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{cuComplex}, CuPtr{Cint},
                    CuPtr{Cint}, CuPtr{cuComplex}, CuPtr{Cint}, CuPtr{Cint},
                    cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end

@checked function cusparseZcsr2csc(handle, m, n, nnz, csrSortedVal, csrSortedRowPtr,
                                   csrSortedColInd, cscSortedVal, cscSortedRowInd,
                                   cscSortedColPtr, copyValues, idxBase)
    initialize_api()
    @runtime_ccall((:cusparseZcsr2csc, libcusparse()), cusparseStatus_t,
                   (cusparseHandle_t, Cint, Cint, Cint, CuPtr{cuDoubleComplex},
                    CuPtr{Cint}, CuPtr{Cint}, CuPtr{cuDoubleComplex}, CuPtr{Cint},
                    CuPtr{Cint}, cusparseAction_t, cusparseIndexBase_t),
                   handle, m, n, nnz, csrSortedVal, csrSortedRowPtr, csrSortedColInd,
                   cscSortedVal, cscSortedRowInd, cscSortedColPtr, copyValues, idxBase)
end


text:
using Distributions

function test_gaussian_kl_loss_layer(backend::Backend, T, eps)
    println("-- Testing GaussianKLLossLayer on $(typeof(backend)){$T}...")

    ############################################################
    # Prepare Data for Testing
    ############################################################
    tensor_dim = abs(rand(Int)) % 4 + 2
    dims = tuple((abs(rand(Int,tensor_dim)) % 6 + 6)...)
    println("    > $dims")
    mus = rand(T, dims)
    sigmas = sqrt(rand(T, dims).^2)

    ############################################################
    # Setup
    ############################################################
    weight = 1.1
    layer  = GaussianKLLossLayer(; bottoms=[:predictions, :labels], weight=weight)
    mu_blob  = make_blob(backend, T, dims)
    sigma_blob = make_blob(backend, T, dims)

    mu_diff_blob  = make_blob(backend, T, dims)
    sigma_diff_blob  = make_blob(backend, T, dims)

    copy!(mu_blob, mus)
    copy!(sigma_blob, sigmas)
    inputs = Blob[mu_blob, sigma_blob]
    diffs = Blob[mu_diff_blob, sigma_diff_blob]

    state = setup(backend, layer, inputs, diffs)

    forward(backend, state, inputs)

    n = length(mu_blob)
    loss = 0.5(sum(mus.^2 + sigmas.^2 - 2log(sigmas)) - n)
    loss *= weight/get_num(mu_blob)
    @test -eps < loss-state.loss < eps

    backward(backend, state, inputs, diffs)
    grad = mus
    grad *= weight/get_num(mu_blob)
    diff = similar(grad)
    copy!(diff, diffs[1])
    @test all(-eps .< grad - diff .< eps)

    grad = sigmas - 1./sigmas
    grad *= weight/get_num(mu_blob)
    diff = similar(grad)
    copy!(diff, diffs[2])
    @test all(-eps .< grad - diff .< eps)

    shutdown(backend, state)
end

function verify_kl(q_dist, p_dist, claimed_kl)
    sample_size = 2_000_000
    q_samples = rand(q_dist, sample_size)
    empirical_kl_samples = logpdf(q_dist, q_samples) - logpdf(p_dist, q_samples)
    empirical_kl = mean(empirical_kl_samples)
    tol = 4 * std(empirical_kl_samples) / sqrt(sample_size)
    min_diff = 1e-2 * std(empirical_kl_samples) / sqrt(sample_size)
    @test_approx_eq_eps empirical_kl claimed_kl tol
end

function test_gaussian_kl_loss_layer_2(backend::Backend, T, eps)
    println("-- Testing GaussianKLLossLayer again on $(typeof(backend)){$T}...")

    ############################################################
    # Prepare Data for Testing
    ############################################################
    tensor_dim = 2
    dims = (p, n) = (5, 6)
    println("        > $dims")

    ############################################################
    # Setup
    ############################################################
    layer = GaussianKLLossLayer(bottoms=[:mean_dummy, :sd_dummy])

    mu = 5rand(dims) - 2
    sigma = 2rand(dims) + 0.01

    inputs = Blob[
        make_blob(backend, convert(Array{T}, mu)),    # mean blob
        make_blob(backend, convert(Array{T}, sigma))]    # sd blob

    diffs = Blob[
        make_blob(backend, rand(T, dims) + 4), # mean diff
        make_blob(backend, 4rand(T, dims) - 7)]    # sd diff

    state = setup(backend, layer, inputs, diffs)

    ############################################################
    # Forward Propagation
    ############################################################
    forward(backend, state, inputs)

    loss = 0.
    for i in 1:n, j in 1:p
        loss -= 1 + 2log(sigma[j, i]) - mu[j, i]^2 - sigma[j, i]^2
    end
    loss /= 2dims[end]
    @test loss >= 0  # it's a KL divergence
    @test -eps < loss - state.loss < eps

    q_dist = MvNormal(mu[:], sigma[:])
    p_dist = MvNormal(zeros(n * p), 1)
    verify_kl(q_dist, p_dist, n * state.loss)

    ############################################################
    # Backward Propagation
    ############################################################
    backward(backend, state, inputs, diffs)
    grad_mu = mu / dims[end]
    grad_sigma = (sigma - 1./sigma) ./ dims[end]

    got_diff_mu = Array(T, dims...)
    copy!(got_diff_mu, diffs[1])
    @test all(-eps .< grad_mu - got_diff_mu .< eps)

    got_diff_sigma = Array(T, dims...)
    copy!(got_diff_sigma, diffs[2])
    @test all(-eps .< grad_sigma - got_diff_sigma .< eps)

    for i in 1:2
        input_i = Array(T, size(inputs[i]))
        copy!(input_i, inputs[i])
        diff_i = Array(T, size(diffs[i]))
        copy!(diff_i, diffs[i])

        for k in 1:p, j in 1:n
            x0 = input_i[k, j]
            function f(x::FloatingPoint)
                input_i[k, j] = x
                copy!(inputs[i], input_i)
                forward(backend, state, inputs)
                input_i[k, j] = x0
                copy!(inputs[i], input_i)
                state.loss
            end
            test_deriv(f, x0, diff_i[k, j])
        end
    end

    ############################################################
    # Shutdown
    ############################################################
    shutdown(backend, state)
end

function test_gaussian_kl_loss_layer(backend::Backend)
  test_gaussian_kl_loss_layer(backend, Float32, 5e-2)
  test_gaussian_kl_loss_layer(backend, Float64, 1e-8)
  test_gaussian_kl_loss_layer_2(backend, Float64, 1e-8)
end

if test_gpu
  test_gaussian_kl_loss_layer(backend_gpu)
end

if test_cpu
  test_gaussian_kl_loss_layer(backend_cpu)
end


text:
module Extenders

module GastonExtender

    import Gaston
    export gplot

    function Gaston.gnuplot_init(args::String)

        if Gaston.gnuplot_state.running
            Gaston.gnuplot_exit()
        end

        pr = (0,0,0,0) # stdin, stdout, stderr, pid
        try
            pr = Gaston.popen3(`gnuplot $args`)
        catch
            error("There was a problem starting up gnuplot.")
        end
        # It's possible that `popen3` runs successfully, but gnuplot exits
        # immediately. Double-check that gnuplot is running at this point.
        if Base.process_running(pr[4])
            Gaston.gnuplot_state.running = true
            Gaston.gnuplot_state.fid = pr
            # Start tasks to read and write gnuplot's pipes
            yield()  # get async tasks started (code blocks without this line)
            notify(Gaston.StartPipes)
        else
            error("There was a problem starting up gnuplot.")
        end
    end

    function gplot(args...)
        Gaston.gnuplot_init("--persist")
        Gaston.plot(args...)
        Gaston.gnuplot_exit()
    end
end

end

max_stars_repo_name: Ashymad/JuliaModules

max_issues_repo_name: Ashymad/JuliaModules

max_forks_repo_name: Ashymad/JuliaModules


text:
"""
    tryparsenext_internal(::Type{<:TimeType}, str, pos, len, df::DateFormat, raise=false)

Parses the string according to the directives within the DateFormat. The specified TimeType
type determines the type of and order of tokens returned. If the given DateFormat or string
does not provide a required token a default value will be used. When the string cannot be
parsed the returned value tuple will be null if `raise` is false otherwise an exception will
be thrown.

Returns a 2-element tuple `(values, pos)`:
* `values::Nullable{Tuple}`: A tuple which contains a value for each token as specified by
  the passed in type.
* `pos::Int`: The character index at which parsing stopped.
"""
@generated function tryparsenext_internal{T<:TimeType}(
                                                       ::Type{T}, str::AbstractString, pos::Int, len::Int, df::DateFormat, endchar='\0', raise::Bool=false,
)
    letters = character_codes(df)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]

    output_tokens = CONVERSION_TRANSLATIONS[T]
    output_names = Symbol[genvar(t) for t in output_tokens]
    output_defaults = ([CONVERSION_DEFAULTS[t] for t in output_tokens]...)
    R = typeof(output_defaults)

    # Pre-assign output variables to defaults. Ensures that all output variables are
    # assigned as the value tuple returned from `tryparsenext_core` may not include all
    # of the required variables.
    assign_defaults = Expr[
        quote
            $name = $default
        end
        for (name, default) in zip(output_names, output_defaults)
    ]

    # Unpacks the value tuple returned by `tryparsenext_core` into separate variables.
    value_tuple = Expr(:tuple, value_names...)

    assign_value_till = Expr[begin
        quote
        ($i <= num_parsed) && ($name = unsafe_val[$i])
        end
    end for (i,name) in enumerate(value_names)]

    quote
        values, pos, num_parsed = tryparsenext_core(str, pos, len, df, raise)
        $(assign_defaults...)
        unsafe_val = unsafe_get(values)
        $(assign_value_till...)
        if isnull(values)
            if (pos <= len && str[pos] == endchar) ||
                num_parsed == $(length(value_names))
                # finished parsing and found an extra char,
                # or parsing was terminated by a delimiter
                return Nullable{$R}($(Expr(:tuple, output_names...))), pos
            end
            return Nullable{$R}(), pos
        end
        return Nullable{$R}($(Expr(:tuple, output_names...))), pos
    end
end

"""
    tryparsenext_core(str::AbstractString, pos::Int, len::Int, df::DateFormat, raise=false)

Parses the string according to the directives within the DateFormat. Parsing will start at
character index `pos` and will stop when all directives are used or we have parsed up to
the end of the string, `len`. When a directive cannot be parsed the returned value tuple
will be null if `raise` is false otherwise an exception will be thrown.

Returns a 3-element tuple `(values, pos, num_parsed)`:
* `values::Nullable{Tuple}`: A tuple which contains a value for each `DatePart` within the
  `DateFormat` in the order in which they occur. If the string ends before we finish parsing
  all the directives the missing values will be filled in with default values.
* `pos::Int`: The character index at which parsing stopped.
* `num_parsed::Int`: The number of values which were parsed and stored within `values`.
  Useful for distinguishing parsed values from default values.
"""
@generated function tryparsenext_core(
    str::AbstractString, pos::Int, len::Int, df::DateFormat, raise::Bool=false,
)
    directives = _directives(df)
    letters = character_codes(directives)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]
    value_defaults = ([CONVERSION_DEFAULTS[t] for t in tokens]...)
    R = typeof(value_defaults)

    # Pre-assign variables to defaults. Allows us to use `@goto done` without worrying about
    # unassigned variables.
    assign_defaults = Expr[
        quote
            $name = $default
        end
        for (name, default) in zip(value_names, value_defaults)
    ]

    vi = 1
    parsers = Expr[
        begin
            if directives[i] <: DatePart
                name = value_names[vi]
                nullable = Symbol(:nullable_, name)
                vi += 1
                quote
                    pos > len && @goto done
                    $nullable, next_pos = tryparsenext(directives[$i], str, pos, len, locale)
                    isnull($nullable) && @goto error
                    $name = unsafe_get($nullable)
                    pos = next_pos
                    num_parsed += 1
                    directive_index += 1
                end
            else
                quote
                    pos > len && @goto done
                    nullable_delim, next_pos = tryparsenext(directives[$i], str, pos, len, locale)
                    isnull(nullable_delim) && @goto error
                    pos = next_pos
                    directive_index += 1
                end
            end
        end
        for i in 1:length(directives)
    ]

    quote
        directives = df.tokens
        locale::DateLocale = df.locale

        num_parsed = 0
        directive_index = 1

        $(assign_defaults...)
        $(parsers...)

        pos > len || @goto error

        @label done
        return Nullable{$R}($(Expr(:tuple, value_names...))), pos, num_parsed

        @label error
        if raise
            if directive_index > length(directives)
                throw(ArgumentError("Found extra characters at the end of date time string"))
            else
                d = directives[directive_index]
                throw(ArgumentError("Unable to parse date time. Expected directive $d at char $pos"))
            end
        end
        return Nullable{$R}($(Expr(:tuple, value_names...)), $(VERSION < v"0.6.0-dev")), pos, num_parsed
    end
end


text:
using Pkg.TOML

const SEMAGRAMS_VERSION = let
  project_file = normpath(joinpath(@__DIR__, "..", "Project.toml"))
  project = TOML.parsefile(project_file)
  VersionNumber(project["version"])
end

const BUNDLE_URL = "https://unpkg.com/semagrams@$SEMAGRAMS_VERSION/dist/app.bundle.js"
const BUNDLE_DIR = normpath(joinpath(@__DIR__, "bundles"))
const BUNDLE_PATH = normpath(joinpath(BUNDLE_DIR, "app.bundle.js"))

function download_js_bundle()
  mkpath(BUNDLE_DIR)

  @info "Downloading Semagrams bundle from unpkg"
  download(BUNDLE_URL, BUNDLE_PATH)
end

function build_js_bundle()
  mkpath(BUNDLE_DIR)

  cd(joinpath(@__DIR__, "..", "javascript"))

  @info "Building Semagrams bundle"
  run(`npm run devbuild`)

  run(`cp dist/app.bundle.js ../deps/bundles/`)

  run(`rm -rf dist`)
end

if get(ENV, "SEMAGRAMS_DEV", "false") == "true"
  build_js_bundle()
else
  download_js_bundle()
end


text:
@testset "Test NLP Stopping unconstrained" begin

  # We create a simple function to test
  A = rand(5, 5);
  Q = A' * A

  f(x) = x' * Q * x
  nlp = ADNLPModel(f, zeros(5))
  nlp_at_x = NLPAtX(zeros(5))
  meta = StoppingMeta(optimality0 = 0.0,
                      max_cntrs   = init_max_counters(),
                      optimality_check = (x,y) -> unconstrained_check(x,y))
  stop_nlp = NLPStopping(nlp, meta, nlp_at_x)

  src = StopRemoteControl()
  stop_nlp_src = NLPStopping(nlp, meta, src, nlp_at_x)

  a = zeros(5)
  fill_in!(stop_nlp, a)

  # we make sure that the fill_in! function works properly
  @test obj(nlp, a) == stop_nlp.current_state.fx
  @test grad(nlp, a) == stop_nlp.current_state.gx
  @test stop_nlp.meta.optimality0 == 0.0

  # we make sure the optimality check works properly
  @test stop!(stop_nlp)
  @test stop_nlp.current_state.current_score == unconstrained_check(stop_nlp.pb, stop_nlp.current_state)
  # we make sure the counter of stop works properly
  @test stop_nlp.meta.nb_of_stop == 1

  reinit!(stop_nlp, rstate = true, x = ones(5))
  @test stop_nlp.current_state.x == ones(5)
  @test isnan(stop_nlp.current_state.fx)
  @test stop_nlp.meta.nb_of_stop == 0

  #We know test how to initialize the counter:
  test_max_cntrs = init_max_counters(obj = 2)
  stop_nlp_cntrs = NLPStopping(nlp, max_cntrs = test_max_cntrs)
  @test stop_nlp_cntrs.meta.max_cntrs[:neval_obj] == 2
  @test stop_nlp_cntrs.meta.max_cntrs[:neval_grad] == typemax(Int)
  @test stop_nlp_cntrs.meta.max_cntrs[:neval_sum] == typemax(Int)

  reinit!(stop_nlp.current_state)
  @test unconstrained_check(stop_nlp.pb, stop_nlp.current_state) >= 0.0
  reinit!(stop_nlp.current_state)
  @test unconstrained2nd_check(stop_nlp.pb, stop_nlp.current_state) >= 0.0
  @test stop_nlp.current_state.Hx != nothing

  #We now test the _unbounded_problem_check:
  @test stop_nlp.pb.meta.minimize #we are minimizing
  stop_nlp.current_state.fx = - 1.0e50 #default meta.unbounded_threshold
  stop!(stop_nlp)
  @test :UnboundedPb in status(stop_nlp, list = true) # the problem is unbounded as fx <= - 1.0e50
  stop_nlp.meta.unbounded_pb = false #reinitialize
  #Let us now consider a maximization problem:
  nlp_max = ADNLPModel(NLPModelMeta(5, minimize = false), Counters(), ADNLPModels.ForwardDiffAD(5), f, x->[])
  stop_nlp.pb = nlp_max
  @test !stop_nlp.pb.meta.minimize
  stop!(stop_nlp)
  @test !(:UnboundedPb in status(stop_nlp, list = true)) # the problem is NOT unbounded as fx <= 1.0e50

  #Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
  #for the proper way of defining an ADNLPModel
  nlp_bnd = ADNLPModel(NLPModelMeta(5, x0=zeros(5), lvar=zeros(5), uvar=zeros(5)),
                       Counters(), ADNLPModels.ForwardDiffAD(5), f, x->[])
  stop_bnd = NLPStopping(nlp_bnd)
  fill_in!(stop_bnd, zeros(5))
  @test KKT(stop_bnd.pb, stop_bnd.current_state) == 0.0
  reinit!(stop_bnd.current_state)
  @test optim_check_bounded(stop_bnd.pb, stop_bnd.current_state) == 0.0

  stop_bnd.meta.optimality_check = (x,y) -> NaN
  start!(stop_bnd)
  @test stop_bnd.meta.domainerror == true
  reinit!(stop_bnd, rcounters = true)
  @test neval_grad(stop_bnd.pb) == 0
  @test stop_bnd.meta.domainerror == false
  stop!(stop_bnd)
  @test stop_bnd.meta.domainerror == true

  stop_bnd.meta.optimality_check = (x,y) -> 0.0
  reinit!(stop_bnd, rcounters = true)
  @test neval_grad(stop_bnd.pb) == 0
  fill_in!(stop_bnd, zeros(5), mu = ones(5), lambda = zeros(0))
  @test stop_bnd.current_state.mu == ones(5)
  @test stop_bnd.current_state.lambda == zeros(0)

  update!(stop_bnd.current_state, fx = NaN, current_score = 0.0)
  stop!(stop_bnd)
  @test !isnan(stop_bnd.current_state.fx)

end


text:
using Dates
using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

const path = joinpath(DEPOT_PATH[1], "registries", "General")
RegistryCI.test(path)

# Test the BioJuliaRegistry
RegistryCI.with_temp_depot() do
    Pkg.Registry.add(Pkg.RegistrySpec(url = "https://github.com/BioJulia/BioJuliaRegistry.git"))
    # Test this will validate the BioJuliaRegistry, when providing General as an
    # optional dependency. BJW.
    biopath = joinpath(DEPOT_PATH[1], "registries", "BioJuliaRegistry")
    RegistryCI.test(biopath, registry_deps = ["https://github.com/JuliaRegistries/General.git"])
end


text:
__precompile__()
module MINLPLib

using JuMP
using JSON
using Glob

include("features.jl")

METAATTRS = ["LIBRARY", "NAME",
            "NVARS",
            "NBINVARS", "NINTVARS", "NNLVARS",
            "NCONS",
            "NLINCONS", "NNLCONS",
            "OBJSENSE", "OBJTYPE", "NLOPERANDS",
            "STATUS",
            "OBJVAL", "OBJBOUND",
            "SOURCE"]

PROTECTED_LIBS = []

minlplib_dir = joinpath(dirname(pathof(MINLPLib)), "..")

special_instances = String[]

export fetch_model

end #module

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function RK4(SLPDE::SemiLinearPDE{1, S}, T::Number, n::Int) where S
    U0 = SLPDE.U0[1]
    U = deepcopy(U0) # solution spherical harmonic coefficients
    V = zero(U) # solution Fourier coefficients
    F = zero(U) # solution function values
    NU = zero(U) # nonlinearity spherical harmonic coefficients
    NV = zero(U) # nonlinearity Fourier coefficients
    NF = zero(U) # nonlinearity function values

    # RK4 stages
    K1 = zero(U)
    K2 = zero(U)
    K3 = zero(U)
    K4 = zero(U)
    NK1 = zero(U)
    NK2 = zero(U)
    NK3 = zero(U)

    P = plan_sph2fourier(U; sketch = :none)
    Ps = FastTransforms.plan_synthesis(U)
    Pa = FastTransforms.plan_analysis(U)

    L = create_linear_operator(SLPDE.L, size(U, 1), size(U, 2)÷2 + 1)
    h = S(T)/n

    for k=1:n
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, U, V, F, NF, NV, NU)
        K1 .= L.*U .+ NU
        NK1 .= U .+ (0.5h).*K1
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK1, V, F, NF, NV, NU)
        K2 .= L.*NK1 .+ NU
        NK2 .= U .+ (0.5h).*K2
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK2, V, F, NF, NV, NU)
        K3 .= L.*NK2 .+ NU
        NK3 .= U .+ h.*K3
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK3, V, F, NF, NV, NU)
        K4 .= L.*NK3 .+ NU
        U .= U + (h/6.0).*(K1 .+ 2.0.*(K2 .+ K3) .+ K4)
    end

    U
end

function RK4(SLPDE::SemiLinearPDE{2, S}, T::Number, n::Int) where S
    U0 = SLPDE.U0
    U = deepcopy(U0) # solution spherical harmonic coefficients
    V = zero(U[1]), zero(U[2]) # solution Fourier coefficients
    F = zero(U[1]), zero(U[2]) # solution function values
    NU = zero(U[1]), zero(U[2]) # nonlinearity spherical harmonic coefficients
    NV = zero(U[1]), zero(U[2]) # nonlinearity Fourier coefficients
    NF = zero(U[1]), zero(U[2]) # nonlinearity function values

    # RK4 stages
    K1 = zero(U[1]), zero(U[2])
    K2 = zero(U[1]), zero(U[2])
    K3 = zero(U[1]), zero(U[2])
    K4 = zero(U[1]), zero(U[2])
    NK1 = zero(U[1]), zero(U[2])
    NK2 = zero(U[1]), zero(U[2])
    NK3 = zero(U[1]), zero(U[2])

    P = plan_sph2fourier(U[1]; sketch = :none), plan_sph2fourier(U[2]; sketch = :none)
    Ps = FastTransforms.plan_synthesis(U[1]), FastTransforms.plan_synthesis(U[2])
    Pa = FastTransforms.plan_analysis(U[1]), FastTransforms.plan_analysis(U[2])

    L = create_linear_operator(SLPDE.L[1], size(U[1], 1), size(U[1], 2)÷2 + 1), create_linear_operator(SLPDE.L[2], size(U[2], 1), size(U[2], 2)÷2 + 1)
    h = S(T)/n

    for k=1:n
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, U, V, F, NF, NV, NU)
        K1[1] .= L[1].*U[1] .+ NU[1]
        K1[2] .= L[2].*U[2] .+ NU[2]
        NK1[1] .= U[1] .+ (0.5h).*K1[1]
        NK1[2] .= U[2] .+ (0.5h).*K1[2]
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK1, V, F, NF, NV, NU)
        K2[1] .= L[1].*NK1[1] .+ NU[1]
        K2[2] .= L[2].*NK1[2] .+ NU[2]
        NK2[1] .= U[1] .+ (0.5h).*K2[1]
        NK2[2] .= U[2] .+ (0.5h).*K2[2]
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK2, V, F, NF, NV, NU)
        K3[1] .= L[1].*NK2[1] .+ NU[1]
        K3[2] .= L[2].*NK2[2] .+ NU[2]
        NK3[1] .= U[1] .+ h.*K3[1]
        NK3[2] .= U[2] .+ h.*K3[2]
        evaluate_nonlinear_operator!(SLPDE.N, P, Ps, Pa, NK3, V, F, NF, NV, NU)
        K4[1] .= L[1].*NK3[1] .+ NU[1]
        K4[2] .= L[2].*NK3[2] .+ NU[2]
        U[1] .= U[1] + (h/6.0).*(K1[1] .+ 2.0.*(K2[1] .+ K3[1]) .+ K4[1])
        U[2] .= U[2] + (h/6.0).*(K1[2] .+ 2.0.*(K2[2] .+ K3[2]) .+ K4[2])
    end

    U
end


text:
module ParentSelection
    using ..Representation
    # Random parent selection
    function random_parent_selection(population, λ, μ_parents)
        groups_indices = split_indices(length(population.members), λ, μ_parents)
        groups = Array{Representation.Organism, 1}[]
        n_groups = ceil(Int, (μ_parents * λ) / length(population.members))
        for i in 1:n_groups
            group = population.members[groups_indices[i, :]]
            push!(groups, group)
        end
        return groups
    end
    # Aux function
    function split_indices(population_length, λ, μ_parents)
        multiplier = ceil(Int, (λ * μ_parents) / population_length)
        total_size = floor(Int, population_length * multiplier)
        group_size = div(total_size, population_length)
        all_indices = ceil.(Int, rand(total_size) * population_length)
        groups_indices = reshape(all_indices, (group_size, div(total_size, group_size)))
        return groups_indices

    end

end # module ParentSelection


text:
# Functions associated with HamiltonianCPU, which maintains the actual internal
# interaction types and orchestrates energy/field calculations.

function validate_and_clean_interactions(ints::Vector{<:AbstractInteraction}, crystal::Crystal, latsize::Vector{Int64})
    # Validate all interactions
    for int in ints
        if isa(int, QuadraticInteraction)
            b = int.bond

            # Verify that both basis sites indexed actually exist
            if !(1 <= b.i <= nbasis(crystal)) || !(1 <= b.j <= nbasis(crystal))
                error("Provided interaction $(repr(MIME("text/plain"), int)) indexes a non-existent basis site.")
            end

            # Verify that the interactions are symmetry-consistent
            if !is_coupling_valid(crystal, b, int.J)
                println("Symmetry-violating interaction: $(repr(MIME("text/plain"), int)).")
                if b.i == b.j && iszero(b.n)
                    println("Allowed single-ion anisotropy for this atom:")
                else
                    println("Allowed exchange for this bond:")
                end
                print_allowed_coupling(crystal, b; prefix="    ")
                println("Use `print_bond(crystal, bond)` for more information.")
                error("Interaction violates symmetry.")
            end

            # Verify that no bond wraps the entire system
            bs = all_symmetry_related_bonds(crystal, b)
            wraps = any(bs) do b
                any(abs.(b.n) .>= latsize)
            end
            if wraps
                println("Distance-violating interaction: $int.")
                error("Interaction wraps system.")
            end
        end
    end

    return ints
end

"""
    HamiltonianCPU

Stores and orchestrates the types that perform the actual implementations
of all interactions internally.
"""
struct HamiltonianCPU
    ext_field   :: Union{Nothing, ExternalFieldCPU}
    heisenbergs :: Vector{HeisenbergCPU}
    diag_coups  :: Vector{DiagonalCouplingCPU}
    gen_coups   :: Vector{GeneralCouplingCPU}
    dipole_int  :: Union{Nothing, DipoleRealCPU, DipoleFourierCPU}
    spin_mags   :: Vector{Float64}
end

"""
    HamiltonianCPU(ints, crystal, latsize, sites_info::Vector{SiteInfo})

Construct a `HamiltonianCPU` from a list of interactions, converting
each of the interactions into the proper backend type specialized
for the given `crystal` and `latsize`.

Note that `sites_info` must be complete when passed to this constructor.
"""
function HamiltonianCPU(ints::Vector{<:AbstractInteraction}, crystal::Crystal,
                        latsize::Vector{Int64}, sites_info::Vector{SiteInfo};
                        μB=BOHR_MAGNETON::Float64, μ0=VACUUM_PERM::Float64)
    ext_field   = nothing
    heisenbergs = Vector{HeisenbergCPU}()
    diag_coups  = Vector{DiagonalCouplingCPU}()
    gen_coups   = Vector{GeneralCouplingCPU}()
    dipole_int  = nothing
    spin_mags   = [site.S for site in sites_info]

    ints = validate_and_clean_interactions(ints, crystal, latsize)

    for int in ints
        if isa(int, ExternalField)
            if isnothing(ext_field)
                ext_field = ExternalFieldCPU(int, sites_info; μB=μB)
            else
                ext_field.Bgs .+= ExternalFieldCPU(int, sites_info; μB=μB).Bgs
            end
        elseif isa(int, QuadraticInteraction)
            int_impl = convert_quadratic(int, crystal, sites_info)
            if isa(int_impl, HeisenbergCPU)
                push!(heisenbergs, int_impl)
            elseif isa(int_impl, DiagonalCouplingCPU)
                push!(diag_coups, int_impl)
            elseif isa(int_impl, GeneralCouplingCPU)
                push!(gen_coups, int_impl)
            else
                error("Quadratic interaction failed to convert to known backend type.")
            end
        elseif isa(int, DipoleDipole)
            if !isnothing(dipole_int)
                @warn "Provided multiple dipole interactions. Only using last one."
            end
            dipole_int = DipoleFourierCPU(int, crystal, latsize, sites_info; μB=μB, μ0=μ0)
        else
            error("$(int) failed to convert to known backend type.")
        end
    end

    return HamiltonianCPU(
        ext_field, heisenbergs, diag_coups, gen_coups, dipole_int, spin_mags
    )
end

function energy(spins::Array{Vec3}, ℋ::HamiltonianCPU) :: Float64
    E = 0.0
    if !isnothing(ℋ.ext_field)
        E += energy(spins, ℋ.ext_field)
    end
    for heisen in ℋ.heisenbergs
        E += energy(spins, heisen)
    end
    for diag_coup in ℋ.diag_coups
        E += energy(spins, diag_coup)
    end
    for gen_coup in ℋ.gen_coups
        E += energy(spins, gen_coup)
    end
    if !isnothing(ℋ.dipole_int)
        E += energy(spins, ℋ.dipole_int)
    end
    return E
end

"""
Updates `B` in-place to hold the local field on `spins` under `ℋ`,
defined as:

``𝐁_i = -∇_{𝐬_i} ℋ / S_i``

with ``𝐬_i`` the unit-vector variable at site i, and ``S_i`` is
the magnitude of the associated spin.

Note that all `_accum_neggrad!` functions should return _just_ the
``-∇_{𝐬_i} ℋ`` term, as the scaling by spin magnitude happens in
this function. Likewise, all code which utilizes local fields should
be calling _this_ function, not the `_accum_neggrad!`'s directly.
"""
function field!(B::Array{Vec3}, spins::Array{Vec3}, ℋ::HamiltonianCPU)
    fill!(B, SA[0.0, 0.0, 0.0])
    if !isnothing(ℋ.ext_field)
        _accum_neggrad!(B, ℋ.ext_field)
    end
    for heisen in ℋ.heisenbergs
        _accum_neggrad!(B, spins, heisen)
    end
    for diag_coup in ℋ.diag_coups
        _accum_neggrad!(B, spins, diag_coup)
    end
    for gen_coup in ℋ.gen_coups
        _accum_neggrad!(B, spins, gen_coup)
    end
    if !isnothing(ℋ.dipole_int)
        _accum_neggrad!(B, spins, ℋ.dipole_int)
    end

    # Normalize each gradient by the spin magnitude on that sublattice
    for idx in CartesianIndices(B)
        S = ℋ.spin_mags[idx[1]]
        B[idx] /= S
    end
end


text: const DBL_MAX = Float64(0x1.fffffffffffffp+1023)


text:
# code from https://github.com/JuliaGraphics/FreeType.jl/tree/master/gen

using Clang.Generators
using H3_jll

const h3api_dir = normpath(H3_jll.find_artifact_dir(), "include/h3")
const h3api_header = normpath(h3api_dir, "h3api.h")

const H3_VERSION = v"3.7.2"
const h3_dir = normpath(@__DIR__, "h3-$H3_VERSION")
!isdir(h3_dir) && throw(string("need ", h3_dir))
const h3_include_dir = normpath(h3_dir, "src/h3lib/include")
((root, dirs, files),) = walkdir(h3_include_dir)
const skips = ["constants.h", "polygonAlgos.h"]
const h3lib_headers = normpath.(root, filter(x -> endswith(x, ".h") && !in(x, skips), files))

const headers = vcat([h3api_header], h3lib_headers)
const args = vcat(get_default_args(), "-I$h3api_dir")
const options = load_options(joinpath(@__DIR__, "generator.toml"))

ctx = create_context(headers, args, options)

build!(ctx)


text:
# a bunch of functions usefull for working with
# Microchip's RTDM (Real Time Data Monitor) interface
module MicrochipRTDM

include("ExceptionHandler.jl")
include("ParseMap.jl")
include("RTDM.jl")

end # module


text:
export phase_shift, phase_shift!

"""
  phase_shift!(C::SeisChannel)

Phase shift SeisChannel if starttime is not aligned with sampling period.

For example, a channel with starttime of 2016-07-04T00:00:00.008 and a sampling
rate of 100 Hz, the starttime is off the sampling period (0.01) by 2 millisecond
and should be at 2016-07-04T00:00:00.01. The waveform will be phase shifted by
2 milliseconds in the frequency domain. If ϕshift = false, will only correct
starttime of channel.

# Arguments
- `C::SeisChannel`: SeisChannel.
- `ϕshift::Bool`: Boolean to allow/disable phase-shifting.
"""
function phase_shift!(C::SeisChannel; ϕshift::Bool=true)

    # get time offset from sampling rate
    dt = Float32(1. / C.fs)
    off = mod(mod(C.t[1,2],1000000)*Float32(1e-6),dt)
    n = length(C.x)

    if dt - off <= eps(eltype(off))
        off = 0.
    end

    if off != 0.
        if off <= dt / 2.
            off = -off
        else
            off = dt - off
        end

        # phase shift data
        if ϕshift
            freq = FFTW.rfftfreq(length(C.x),C.fs)
            fftdata = rfft(C.x)
            fftdata .= fftdata .* exp.(1im .* 2 .* pi .* freq .* off)
            C.x[:] = irfft(fftdata,length(C.x))
        end

        # reset time of signal
        C.t[1,2] += off * 1e6
    end
    return nothing
end
phase_shift(C::SeisChannel;ϕshift::Bool=true) = (U = deepcopy(C);
                phase_shift!(U,ϕshift=ϕshift);return U)

"""

  phase_shift!(S::SeisData)

For example, a channel with starttime of 2016-07-04T00:00:00.008 and a sampling
rate of 100 Hz, the starttime is off the sampling period (0.01) by 2 millisecond
and should be at 2016-07-04T00:00:00.01. The waveform will be phase shifted by
2 milliseconds in the frequency domain. If ϕshift = false, will only correct
starttime of channel.

This function will not work properly for sampling rates >= 1000 Hz due to Julia
DateTime resolving time only to milliseconds!

# Arguments
- `S::SeisData`: SeisData.
- `ϕshift::Bool`: Boolean to allow/disable phase-shifting.
"""
function phase_shift!(S::SeisData;ϕshift::Bool=true)
    @inbounds for i = 1:S.n
        phase_shift!(S[i],ϕshift=ϕshift)
    end
    return nothing
end
phase_shift(S::SeisData;ϕshift::Bool=true) = (U = deepcopy(S);
            phase_shift!(U,ϕshift=ϕshift);return U)


text:
export ParticleType
export photon, electron, positron, neutron, proton

@enum ParticleType photon=1 electron=2 positron=3 neutron=4 proton=5

for pt in instances(ParticleType)
    fname = Symbol("is", pt)
    @eval $fname(p) = p.typ == $pt
    eval(Expr(:export, fname))
end

function load(path::AbstractString, T)
    open(path) do io
        load(io, T)
    end
end

function compute_u_v_w(u::Float32, v::Float32, sign_w::Float32)
    tmp = Float64(u)^2 + Float64(v)^2
    if tmp <= 1
        w = Float32(sign_w) * Float32(√(1 - tmp))
    else
        w = Float32(0)
        tmp = √(tmp)
        u = Float32(u/tmp)
        v = Float32(v/tmp)
    end
    u,v,w
end

@noinline function call_fenced(f::F, arg::A) where {F,A}
    f(arg)
end

function kwshow(io::IO, o; calle=nameof(typeof(o)))
    print(io, calle, "(")
    for pname in propertynames(o)
        pval = getproperty(o, pname)
        print(io, string(pname), "=")
        show(io, pval)
        print(io, ", ")
    end
    print(io, ")")
end


text:
# This file is a part of SimilaritySearch.jl

export CommonPrefixDissimilarity, GenericLevenshteinDistance, StringHammingDistance, LevenshteinDistance, LcsDistance

"""
    CommonPrefixDissimilarity()

Uses the common prefix as a measure of dissimilarity between two strings
"""
struct CommonPrefixDissimilarity <: SemiMetric
end

"""
    StringHammingDistance()

The hamming distance counts the differences between two equally sized strings
"""
struct StringHammingDistance <: SemiMetric
end

"""
    GenericLevenshteinDistance(;icost, dcost, rcost)

The levenshtein distance measures the minimum number of edit operations to convert one string into another.
The costs insertion `icost`, deletion cost `dcost`, and replace cost `rcost`. Not thread safe, use a copy of for each thread.
"""
struct GenericLevenshteinDistance <: SemiMetric
    icost::Int32 # insertion cost
    dcost::Int32 # deletion cost
    rcost::Int32 # replace cost

    Cpool::Vector{Vector{Int16}}
end

GenericLevenshteinDistance(; icost=1, dcost=1, rcost=1) =
    GenericLevenshteinDistance(icost, dcost, rcost, [Vector{Int16}(undef, 64) for i in 1:Threads.nthreads()])

"""
    LevenshteinDistance()

Instantiates a GenericLevenshteinDistance object to perform traditional levenshtein distance
"""
LevenshteinDistance() = GenericLevenshteinDistance()

"""
    LcsDistance()

Instantiates a GenericLevenshteinDistance object to perform LCS distance
"""
LcsDistance() = GenericLevenshteinDistance(rcost=2)

"""
    common_prefix(a, b)

Computes the length of the common prefix among two strings represented as arrays
"""
function common_prefix(a, b)
    len_a = length(a)
    len_b = length(b)
    i = 1
    min_len = min(len_a, len_b)
    @inbounds while i <= min_len && a[i] == b[i]
    	i += 1
    end

    i - 1
end

"""
    evaluate(::CommonPrefixDissimilarity, a, b)

Computes a dissimilarity based on the common prefix between two strings
"""
evaluate(::CommonPrefixDissimilarity, a, b) = 1.0 - common_prefix(a, b) / min(length(a), length(b))

"""
    evaluate(::GenericLevenshteinDistance, a, b)

Computes the edit distance between two strings, this is a low level function
"""
function evaluate(lev::GenericLevenshteinDistance, a, b)
    if length(a) < length(b)
        a, b = b, a
    end

    alen = length(a)
    blen = length(b)

    alen == 0 && return blen
    blen == 0 && return alen

    C = lev.Cpool[Threads.threadid()]
    resize!(C, blen+1)
    @inbounds for i in 0:blen
        C[i+1] = i
    end

    prevA = 0
    @inbounds for i in 1:alen
        prevA = i
        prevC = C[1]
        j = 1

        while j <= blen
            cost = a[i] == b[j] ? 0 : lev.rcost
            C[j] = prevA
            j += 1
            prevA = min(C[j]+lev.dcost, prevA+lev.icost, prevC+cost)
            prevC = C[j]
        end

        C[j] = prevA
    end

    prevA
end

"""
     evaluate(::StringHammingDistance, a, b)

Computes the hamming distance between two sequences of the same length
"""
function evaluate(::StringHammingDistance, a, b)
    d = 0

    @inbounds for i = 1:length(a)
        d += Int(a[i] != b[i])
    end

    d
end

# function kerrormatch(a::T1, b::T2, errors::Integer)::Bool where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return alen == blen
#     blen == 0 && return true

#     C::Vector{Int} = Vector{Int}(0:blen)

#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
# 	    end

#         C[j] = prevA
#         if prevA <= errors
#             return true
#         end
#     end

#     return false
# end

# function best_match_levenshtein(a::T1, b::T2)::Int where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return blen
#     blen == 0 && return alen

#     C::Vector{Int} = 1:blen |> collect

#     mindist = alen
#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
#         end

#         C[j] = prevA
#         if prevA < mindist
#             mindist = prevA
#         end
#     end

#     return mindist
# end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
# using Revise, Plots
using OrdinaryDiffEq, ForwardDiff, Test
	using BifurcationKit, LinearAlgebra, Parameters, Setfield
	const BK = BifurcationKit
	const FD = ForwardDiff

norminf = x -> norm(x, Inf)

function Fsl!(f, u, p, t = 0)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2
	return f
end

Fsl(x, p) = Fsl!(similar(x), x, p)
dFsl(x, dx, p) = FD.derivative(t -> Fsl(x .+ t .* dx, p), 0.)
jet = BK.getJet(Fsl; matrixfree = false)

par_sl = (r = 0.5, μ = 0., ν = 1.0, c3 = 1.0, c5 = 0.0,)
par_hopf = (@set par_sl.r = 0.1)
u0 = [.001, .001]

function FslMono!(f, x, p, t)
	u = x[1:2]
	du = x[3:4]
	Fsl!(f[1:2], u, p, t)
	f[3:4] .= dFsl(u, du, p)
end
####################################################################################################
# continuation
optconteq = ContinuationPar(ds = -0.01, detectBifurcation = 3, pMin = -0.5, nInversion = 4)
br, = continuation(Fsl, u0, par_hopf, (@lens _.r), optconteq)
####################################################################################################
prob = ODEProblem(Fsl!, u0, (0., 100.), par_hopf)
probMono = ODEProblem(FslMono!, vcat(u0, u0), (0., 100.), par_hopf)
####################################################################################################
sol = solve(probMono, KenCarp4(), abstol =1e-9, reltol=1e-6)
sol = solve(prob, KenCarp4(), abstol =1e-9, reltol=1e-6)
# plot(sol[1,:], sol[2,:])
####################################################################################################
section(x, T) = x[1] #* x[end]
section(x, T, dx, dT) = dx[1] #* x[end]
# standard simple shooting
M = 1
dM = 1
_pb = ShootingProblem(prob, KenCarp4(), 1, section; abstol =1e-10, reltol=1e-9)

initpo = [0.13, 0., 6.]
res = _pb(initpo, par_hopf)

# test of the differential of the shooting method

_dx = rand(3)
resAD = FD.derivative(z -> _pb(initpo .+ z .* _dx, par_hopf), 0.)
resFD = (_pb(initpo .+ 1e-8 .* _dx, par_hopf) - _pb(initpo, par_hopf)) .* 1e8
resAN = _pb(initpo, par_hopf, _dx; δ = 1e-8)
@test norm(resAN - resFD, Inf) < 4e-6
@test norm(resAN - resAD, Inf) < 4e-6
####################################################################################################
# test shooting interface M = 1
_pb = ShootingProblem(prob, Rodas4(), [initpo[1:end-1]]; abstol =1e-10, reltol=1e-9)
res = _pb(initpo, par_hopf)
res = _pb(initpo, par_hopf, initpo)

# test the jacobian of the functional in the case M=1
_Jad = FD.jacobian( x -> _pb(x, par_hopf), initpo)
_Jana = _pb(Val(:JacobianMatrix), initpo, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-7

_pb2 = ShootingProblem(prob, Rodas4(), probMono, Rodas4(autodiff=false), [initpo[1:end-1]]; abstol = 1e-10, reltol = 1e-9)
res = _pb2(initpo, par_hopf)
res = _pb2(initpo, par_hopf, initpo)
BK.isSimple(_pb2)

# we test this using Newton - Continuation
optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 20)
# deflationOp = BK.DeflationOperator(2, (x,y) -> dot(x[1:end-1], y[1:end-1]), 1.0, [zeros(3)])
outpo, _, flag,_ = newton(_pb,
	initpo, par_hopf,
	optn;
	jacobianPO = :autodiffDense,
	normN = norminf)
	@test flag

BK.getPeriod(_pb, outpo, par_hopf)
BK.getAmplitude(_pb, outpo, par_hopf)
BK.getMaximum(_pb, outpo, par_hopf)
BK.getPeriodicOrbit(_pb, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds= -0.01, pMax = 4.0, maxSteps = 30, detectBifurcation = 2, nev = 2, newtonOptions = (@set optn.tol = 1e-7), precisionStability = 1e-5)
	br_pok2, = continuation(_pb,
		outpo, par_hopf, (@lens _.r),
		opts_po_cont;
		tangentAlgo = BorderedPred(),
		verbosity = 0, plot = false,
		jacobianPO = :autodiffDense,
		recordFromSolution = (u, p) -> norm(u[1:2]),
		normC = norminf)
# plot(br_pok2)
####################################################################################################
# test automatic branch switching
br_pok2, = continuation(jet..., br, 1, opts_po_cont, ShootingProblem(1, prob, KenCarp4();  abstol = 1e-10, reltol = 1e-9); normC = norminf, jacobianPO = :autodiffDense)

# test matrix-free computation of floquet coefficients
eil = EigKrylovKit(dim = 2, x₀=rand(2))
opts_po_contMF = @set opts_po_cont.newtonOptions.eigsolver = eil
opts_po_contMF = @set opts_po_cont.detectBifurcation = 0
br_pok2, = continuation(jet...,br,1, opts_po_contMF, ShootingProblem(1, prob, Rodas4(); abstol = 1e-10, reltol = 1e-9); jacobianPO = :autodiffDense, normC = norminf, plot=false)

# case with 2 sections
br_pok2_s2, = continuation(jet..., br, 1, opts_po_cont, ShootingProblem(2, prob, KenCarp4();  abstol = 1e-10, reltol = 1e-9); normC = norminf, jacobianPO = :autodiffDense)

####################################################################################################
# test shooting interface M > 1
_pb = ShootingProblem(prob, KenCarp4(), [initpo[1:end-1],initpo[1:end-1],initpo[1:end-1]]; abstol =1e-10, reltol=1e-9)
initpo = [0.13, 0, 0, 0.13, 0, 0.13 , 6.3]
res = _pb(initpo, par_hopf)
res = _pb(initpo, par_hopf, initpo)
# test the jacobian of the functional in the case M=1
_Jad = FD.jacobian( x -> _pb(x, par_hopf), initpo)
_Jana = _pb(Val(:JacobianMatrix), initpo, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-7
####################################################################################################
# Single Poincaré Shooting with hyperplane parametrization
normals = [[-1., 0.]]
centers = [zeros(2)]

probPsh = PoincareShootingProblem(2, prob, Rodas4(), probMono, Rodas4(); abstol=1e-10, reltol=1e-9)
probPsh = PoincareShootingProblem(2, prob, Rodas4(); rtol = abstol=1e-10, reltol=1e-9)

probPsh = PoincareShootingProblem(prob, Rodas4(),
		probMono, Rodas4(),
		normals, centers; abstol = 1e-10, reltol = 1e-9)

initpo_bar = BK.R(probPsh, [0, 0.4], 1)

BK.E(probPsh, [1.0,], 1)
initpo_bar = [0.4]

probPsh(initpo_bar, par_hopf)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x-> probPsh(x, par_hopf), initpo_bar)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-4

ls = DefaultLS()
	eil = EigKrylovKit(dim = 1, x₀ = rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-8,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = BK.DeflationOperator(2, dot, 1.0, [zero(initpo_bar)])
	outpo, _, flag, = newton(probPsh,
			initpo_bar, par_hopf,
			optn; normN = norminf, jacobianPO = :autodiffDenseAnalytical)
	@test flag

BK.getPeriod(probPsh, outpo, par_hopf)
BK.getAmplitude(probPsh, outpo, par_hopf)
BK.getMaximum(probPsh, outpo, par_hopf)
BK.getPeriodicOrbit(probPsh, outpo, par_hopf)

probPsh = PoincareShootingProblem(prob, Rodas4(),
		# probMono, Rodas4(autodiff=false),
		normals, centers; abstol = 1e-10, reltol = 1e-9)

probPsh(outpo, par_hopf)
probPsh(outpo, par_hopf, outpo)
probPsh([0.30429879744900434], par_hopf)
probPsh([0.30429879744900434], (r = 0.09243096156871472, μ = 0.0, ν = 1.0, c3 = 1.0, c5 = 0.0))
BK.evolve(probPsh.flow,[0.0, 0.30429879744900434], (r = 0.094243096156871472, μ = 0.0, ν = 1.0, c3 = 1.0, c5 = 0.0), Inf64) # this gives an error in DiffEqBase

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.015, ds= 0.01, pMax = 4.0, maxSteps = 30, newtonOptions = setproperties(optn; tol = 1e-7, eigsolver = eil), detectBifurcation = 0)
	# br_pok2, = continuation(
	# 	probPsh, outpo, par_hopf, (@lens _.r),
	# 	opts_po_cont; verbosity = 2,
	# 	tangentAlgo = BorderedPred(),
	# 	jacobianPO = :autodiffDenseAnalytical,
	# 	plot = false, normC = norminf)
# plot(br_pok2)
####################################################################################################
# normals = [[-1., 0.], [1, -1]]
# centers = [zeros(2), zeros(2)]
# initpo_bar = [1.04, -1.04/√2]

normals = [[-1., 0.], [1, 0]]
centers = [zeros(2), zeros(2)]
initpo_bar = [0.2, -0.2]

probPsh = PoincareShootingProblem(prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9)
# version with analytical jacobian
probPsh2 = PoincareShootingProblem(prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9, δ = 0)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x-> probPsh(x, par_hopf), initpo_bar)
# _Jphifd = BifurcationKit.finiteDifferences(x->probPsh.flow(x, par_hopf, Inf64), [0,0.4]; δ=1e-8)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-5

ls = DefaultLS()
	eil = EigKrylovKit(dim = 1, x₀=rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = DeflationOperator(2.0, dot, 1.0, [zero(initpo_bar)])
	outpo, = newton(probPsh2, initpo_bar, par_hopf, optn; normN = norminf, jacobianPO = :autodiffDenseAnalytical)
	outpo, = newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf, jacobianPO = :autodiffDenseAnalytical)

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.01, pMax = 4.0, maxSteps = 50, newtonOptions = (@set optn.tol = 1e-9), detectBifurcation = 3, nev = 2)
	br_pok2, = continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0,
		tangentAlgo = BorderedPred(),
		jacobianPO = :autodiffDenseAnalytical,
		plot = false, normC = norminf)

####################################################################################################
normals = [[-1., 0.], [1, 0], [0, 1]]
centers = [zeros(2), zeros(2), zeros(2)]
initpo = [[0., 0.4], [0, -.3], [0.3, 0]]

probPsh = PoincareShootingProblem(prob, KenCarp4(), normals, centers; abstol=1e-10, reltol=1e-9)

initpo_bar = reduce(vcat, [BK.R(probPsh, initpo[ii], ii) for ii in eachindex(centers)])
# same with projection function
initpo_bar = reduce(vcat, BK.projection(probPsh, initpo))

# test of the other projection function
BK.projection(probPsh, reduce(hcat, initpo)')

probPsh(initpo_bar, par_hopf; verbose = true)

# test of the analytical formula for jacobian of the functional
_Jad = BifurcationKit.finiteDifferences( x -> probPsh(x, par_hopf), initpo_bar)
_Jana = probPsh(Val(:JacobianMatrix), initpo_bar, par_hopf)
@test norm(_Jad - _Jana, Inf) < 1e-5

outpo, = newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf, jacobianPO = :autodiffDenseAnalytical)

for ii=1:length(normals)
	@show BK.E(probPsh	, [outpo[ii]], ii)
end

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.005, pMax = 4.0, maxSteps = 10, newtonOptions = setproperties(optn; tol = 1e-8), detectBifurcation = 3)
	br_hpsh, upo , _= continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; normC = norminf, jacobianPO = :autodiffDenseAnalytical)
# plot(br_hpsh)
####################################################################################################
# test automatic branch switching with most possible options
# calls with analytical jacobians
br_psh, = continuation(jet..., br, 1, (@set opts_po_cont.ds = 0.005), PoincareShootingProblem(2, prob, KenCarp4(); abstol=1e-10, reltol=1e-9); normC = norminf, jacobianPO = :autodiffDenseAnalytical)

# test Iterative Floquet eigen solver
@set! opts_po_cont.newtonOptions.eigsolver.dim = 20
@set! opts_po_cont.newtonOptions.eigsolver.x₀ = rand(2)
br_sh, = continuation(jet..., br, 1, ContinuationPar(opts_po_cont; ds = 0.005, saveSolEveryStep = 1), ShootingProblem(2, prob, KenCarp4(); abstol=1e-10, reltol=1e-9); normC = norminf, jacobianPO = :autodiffDenseAnalytical)

# test MonodromyQaD
# BK.MonodromyQaD(br_sh.functional, br.sol)

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo_bar), maxiter = 500, verbose = false)
@set! opts_po_cont.detectBifurcation = 0
@set! opts_po_cont.newtonOptions.linsolver = ls

for M in (1,2), jacobianPO in (:autodiffMF, :MatrixFree, :autodiffDenseAnalytical, :FiniteDifferencesDense)
	@info M, jacobianPO, "PS"

	# specific to Poincaré Shooting
	jacPO = jacobianPO == :autodiffMF ? :FiniteDifferencesDense : jacobianPO

	br_psh, = continuation(jet..., br, 1, (@set opts_po_cont.ds = 0.005), PoincareShootingProblem(M, prob, Rodas4P(); abstol=1e-10, reltol=1e-9, parallel = true); normC = norminf, updateSectionEveryStep = 2, linearAlgo = BorderingBLS(solver = (@set ls.N = M), checkPrecision = false), jacobianPO = jacPO, verbosity = 0)

	br_ssh, = continuation(jet..., br, 1, (@set opts_po_cont.ds = 0.005),
	ShootingProblem(M, prob, Rodas4P(); abstol=1e-10, reltol=1e-9, parallel = true); normC = norminf, updateSectionEveryStep = 2, jacobianPO = jacobianPO,
	linearAlgo = BorderingBLS(solver = (@set ls.N = 2M + 1), checkPrecision = false), verbosity = 0)
end