text:
#=--------------------------------------------------------------------

syntheticimages - Functions for creating various synthetic test images
                  for evaluating edge detectors. Most of these images
                  cause considerable grief for gradient based
                  operators.

Copyright (c) 2015-2017 Peter Kovesi
peterkovesi.com

MIT License:

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

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

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

PK August 2015   Original porting from MATLAB to Julia
   November 2017 Julia 0.6
   October 2018  Julia 0.7/1.0

---------------------------------------------------------------------=#

export step2line, circsine, starsine, noiseonf
export nophase, quantizephase, swapphase

#----------------------------------------------------------------------

# step2line
"""
A phase congruent test image that interpolates from a step to a line.

Generates a test image where the feature type changes from a step edge to a line
feature from top to bottom.  Gradient based edge detectors will only correctly
mark the step-like feature towards the top of the image and incorrectly mark two
features towards the bottom of the image whereas phase congruency will correctly
mark a single feature from top to bottom.  In general, natural images contain a
roughly uniform distribution of the full continuum of feature types from step to
line.

```
Usage:
    img = step2line(sze; nscales=50, ampexponent=-1, ncycles=1.5, phasecycles=0.25)

Arguments:
      sze::Integer - Number of rows in test image, defaults to 512.

Keyword Arguments:
  nscales::Integer - No of Fourier components used to construct the signal.
                     Defaults to 50.
 ampexponent::Real - Decay exponent of amplitude with frequency.
                     A value of -1 will produce amplitude inversely
                     proportional to frequency (corresponds to step feature).
                     A value of -2 will result in the line feature
                     appearing as a triangular waveform. Defaults to -1.
     ncycles::Real - Number of wave cycles across the width of the image.
                     Defaults to 1.5
 phasecycles::Real - Number of feature type phase cycles going vertically
                     down the image. Defaults to 0.25 giving a sequence of feature
                     phase congruency angle varying from 0 to pi/2.
Returns:
   img::Array{Float64,2} - The test image.

Examples of use:
  > img = step2line()                              # Default pattern
  > img = step2line(ncycles=3, ampexponent=-1.5);  # 3 cycles, 'soft' step to line
  > img = step2line(ncycles=3, ampexponent=-1.5, phasecycles = 3);

```
See also:  [`circsine`](@ref), [`starsine`](@ref)
"""
function step2line(sze::Integer = 512; nscales::Integer=50, ampexponent::Real = -1,
                   ncycles::Real = 1.5, phasecycles::Real = 0.25)

    # Construct vector of angles over desired number of cycles
    x = (0:(sze-1))/(sze-1)*ncycles*2*pi

    img = zeros(sze,sze)
    phaseoffset = 0.0

    for row = 1:sze
        signal = zeros(1,sze)
        for scale = 1:2:(nscales*2-1)
            img[row,:] .+= scale^float(ampexponent).*sin.(scale.*x .+ phaseoffset)
        end
        phaseoffset += phasecycles*2*pi/sze
    end

#=
    figure
    colormap(gray)
    imagesc(im), axis('off') , title('step to line feature interpolation')

    range = 3.2
    s = 'Profiles having phase congruency at 0/180, 30/210, 60/240 and 90/270 degrees'

    figure
    subplot(4,1,1), plot(im(1,:)) , title(s), axis([0,sze,-range,range]), axis('off')
    subplot(4,1,2), plot(im(fix(sze/3),:)), axis([0,sze,-range,range]), axis('off')
    subplot(4,1,3), plot(im(fix(2*sze/3),:)), axis([0,sze,-range,range]), axis('off')
    subplot(4,1,4), plot(im(sze,:)), axis([0,sze,-range,range]), axis('off')
=#
    return img
end

#----------------------------------------------------------------------
# circsine
"""
Generate a phase congruent circular sine wave grating.

Useful for testing the isotropy of response of a feature dectector.

```
Usage:    img = circsine(sze; wavelength = 40, nscales = 50, ampexponent = -1,
                          offset = 0, p = 2, trim = false)

Arguments:
   sze::Integer  - The size of the square image to be produced. Defaults to 512.

Keyword arguments:
 wavelength::Real - The wavelength in pixels of the sine wave. Defaults to 40.
 nscales::Integer - No of Fourier components used to construct the
                    signal. This is typically 1, if you want a simple sine
                    wave, or >50 if you want to build a phase congruent
                    waveform. Defaults to 50.
ampexponent::Real - Decay exponent of amplitude with frequency.
                    A value of -1 will produce amplitude inversely
                    proportional to frequency (this will produce a step
                    feature if offset is 0)
                    A value of -2 with an offset of pi/2 will result in a
                    triangular waveform.  Defaults to -1;
     offset::Real - Angle of phase congruency at which the features of the
                    star pattern will be generated at. This controls the feature type.
                    0 for a step-like feature, pi/2 for a line/triangular-like feature.
                    Defaults to 0. If nscales = 1 use pi/2 to get continuity
                    at the centre.
      p::Integer  - Optional parameter specifying the norm to use in
                    calculating the radius from the centre. This defaults to
                    2, resulting in a circular pattern.  Large values gives
                    a square pattern
      trim::Bool  - Optional boolean flag indicating whether you want the
                    circular pattern trimmed from the corners leaving
                    only complete circles. Defaults to false.
Returns:
   img::Array{Float64,2} - The test image.

Examples:
> circsine(nscales = 1) - A simple circular sine wave pattern
> circsine(nscales = 50, ampexponent = -1, offset =  0)     - Square waveform
> circsine(nscales = 50, ampexponent = -2, offset = pi/2)   - Triangular waveform
> circsine(nscales = 50, ampexponent = -1.5, offset = pi/4) - Something in between
                                                              square and triangular
> circsine(nscales = 50, ampexponent = -1.5, offset = 0)    - Looks like a square but is not.
```
See also: [`starsine`](@ref), [`step2line`](@ref)
"""
function circsine(sze::Integer=512; wavelength::Real=40, nscales::Integer=50,
                  ampexponent::Real=-1, offset::Real=0, p::Integer=2,
                  trim::Bool=false)

    if isodd(p)
	error("p should be an even number")
    end

    # Place origin at centre for odd sized image, and below and the the
    # right of centre for an even sized image
    if iseven(sze)
	l = -sze/2
	u = sze/2-1
    else
	l = -(sze-1)/2
	u = (sze-1)/2
    end

    # Grid of radius values
    r = [(x.^p + y.^p).^(1.0/p) for x = l:u, y = l:u]

    img = zeros(size(r))

    for scale = 1:2:(2*nscales-1)
        @. img += scale^float(ampexponent) * sin(scale * r * 2*pi/wavelength + offset)
    end

    if trim     # Remove circular pattern from the 'corners'
	cycles = floor(sze/2/wavelength) # No of complete cycles within sze/2
	@. img *= (r < cycles*wavelength) + (r >= cycles*wavelength)
    end

    return img
end

#----------------------------------------------------------------------
# starsine
"""
Generate a phase congruent star shaped sine wave grating.

Useful for testing the behaviour of feature detectors at line junctions.

```
Usage:    img = starsine(sze; ncycles=10, nscales=50, ampexponent=-1, offset=0)

Argument:
     sze::Integer - The size of the square image to be produced. Defaults to 512.

Keyword arguments:
    ncycles::Real - The number of sine wave cycles around centre point.
                    Typically an integer, but any value can be used.
 nscales::Integer - No of fourier components used to construct the
                    signal. This is typically 1, if you want a simple sine
                    wave, or >50 if you want to build a phase congruent
                    waveform.  Defaults to 50.
ampexponent::Real - Decay exponent of amplitude with frequency.
                    A value of -1 will produce amplitude inversely
                    proportional to frequency (this will produce a step
                    feature if offset is 0)
                    A value of -2 with an offset of pi/2 will result in a
                    triangular waveform.
     offset::Real - Angle of phase congruency at which the features of the
                    star pattern will be generated at. This controls the feature type.
                    0 for a step-like feature, pi/2 for a line/triangular-like feature.
Returns:
   img::Array{Float64,2} - The test image.

Examples:
> starsine(nscales = 1) - A simple sine wave pattern radiating out
                          from the centre. Use 'offset' if you wish to
                          rotate it a bit.
> starsine(nscales = 50, ampexponent = -1, offset =  0)     - Square waveform
> starsine(nscales = 50, ampexponent = -2, offset = pi/2)   - Triangular waveform
> starsine(nscales = 50, ampexponent = -1.5, offset = pi/4) - Something in between
                                                              square and triangular
> starsine(nscales = 50, ampexponent = -1.5, offset = 0)    - Looks like a square but is not.
```
See also: [`circsine`](@ref), [`step2line`](@ref)
"""
function starsine(sze::Integer=512; ncycles::Real=10, nscales::Integer=50,
                  ampexponent::Real=-1, offset::Real=0)

    # Place origin at centre for odd sized image, and below and the the
    # right of centre for an even sized image
    if iseven(sze)
	l = -sze/2
	u = sze/2-1
    else
	l = -(sze-1)/2
	u = (sze-1)/2
    end

    # Grid of angular values
    theta = [atan(y,x) for x = l:u, y = l:u]

    img = zeros(size(theta))

    for scale = 1:2:(nscales*2 - 1)
        @. img += scale^float(ampexponent)*sin(scale*ncycles*theta + offset)
    end

    return img
end

#----------------------------------------------------------------------
# noiseonf
"""
Create \$1/f^p\$ spectrum noise images.

When displayed as a surface these images also generate great landscape
terrain.
```
Usage: img = noiseonf(sze, p)

Arguments:
    sze::Tuple{Integer, Integer} or ::Integer
              - A tuple (rows, cols) or single value specifying size of
                image to produce.
      p::Real - Exponent of spectrum decay = 1/(f^p)

Returns:
   img::Array{Float64,2} - The noise image with specified spectrum.

Reference values for p:
            p = 0   - raw Gaussian noise image.
              = 1   - gives the supposedly 1/f 'standard' drop-off for
                      'natural' images.
              = 1.5 - seems to give the most interesting 'cloud patterns'.
              = > 2 - produces 'blobby' images.
```
"""
function noiseonf(sze::Tuple{Integer, Integer}, p::Real)

    (rows,cols) = sze

    # Generate an image of random Gaussian noise, mean 0, std dev 1.
    img = randn(rows,cols)

    imgfft = fft(img)
    mag = abs.(imgfft)               # Get magnitude
    phase = imgfft./mag              # and phase

    # Construct the amplitude spectrum filter
    # Add 1 to avoid divide by 0 problems later
    radius = filtergrid(rows,cols) * max(rows,cols) .+ 1
    filter = 1 ./ (radius.^p)

    # Reconstruct fft of noise image, but now with the specified amplitude
    # spectrum
    newfft =  filter .* phase
    img .= real.(ifft(newfft))

    return img
end

function noiseonf(sze::Integer, p::Real)
    return noiseonf((sze, sze), p)
end

#--------------------------------------------------------------------
# nophase
"""
Randomize image phase leaving amplitude spectrum unchanged.
```
Usage:   newimg = nophase(img)

Argument:       img::AbstractArray{T,2} where T <: Real - Input image

Returns:     newimg::Array{Float64,2} - Image with randomized phase
```
In general most images will be destroyed by this transform.  However, some
textures are reproduced in an 'amplitude only' image quite well.  Typically
these are textures which have an amplitude spectrum that have a limited number
of isolated peaks. That is, a texture made up from a limited number of strong
harmonics.

See also: [`noiseonf`](@ref), [`quantizephase`](@ref), [`swapphase`](@ref)
"""
function nophase(img::AbstractArray{T,2}) where T <: Real

    # Take FFT, get magnitude.
    IMG = fft(img)
    mag = abs.(IMG)

    # Generate random phase values
    # ** ? Should I just randomize the phase for half the spectrum and
    # generate the other half as the complex conjugate ? **
    phaseAng = rand(Float64, size(img))*2*pi
    phase = cos.(phaseAng) .+ im.*sin.(phaseAng)

    # Reconstruct fft of image using original amplitude and randomized phase and
    # invert.  Note that we must take the real part because the phase
    # randomization applied above did not respect the complex conjugate symmetry
    # that we should have in a real valued signal
    newfft =  mag .* phase
    return real.(ifft(newfft))
end

#--------------------------------------------------------------------
# quantizephase
"""
Quantize phase values in an image.
```
Usage:  qimg = quantizephase(img, N)

Arguments: img::Array{T,2} where T <: Real - Image to be processed
                     N::Integer - Desired number of quantized phase values

Returns:  qimg::Array{Float64,2} - Phase quantized image
```
Phase values in an image are important.  However, despite this, they can be
quantized very heavily with little perceptual loss.  The value of N can be
as low as 4, or even 3!  Using N = 2 is also worth a look.

See also: [`swapphase`](@ref)
"""
function quantizephase(img::AbstractArray{T,2}, N::Integer) where T <: Real

    IMG = fft(img)
    amp = abs.(IMG)
    phase = angle.(IMG)

    # Quantize the phase values as follows:
    # Add pi - .001 so that values range [0 - 2pi)
    # Divide by 2pi so that values range [0 - 1)
    # Scale by N so that values range [0 - N)
    # Round twoards 0 using floor giving integers [0 - N-1]
    # Scale by 2*pi/N to give N discrete phase values  [0 - 2*pi)
    # Subtract pi so that discrete values range [-pi - pi)
    # Add pi/N to counteract the phase shift induced by rounding towards 0
    # using floor
    phase = floor.( (phase .+ pi .- 0.001)/(2*pi) * N) * (2*pi)/N .- pi  .+ pi/N

    # Reconstruct Fourier transform with quantized phase values and take inverse
    # to obtain the new image.
    QIMG = amp.*(cos.(phase) .+ im*sin.(phase))
    return real.(ifft(QIMG))
end

#--------------------------------------------------------------------
# swapphase
"""
Demonstrates phase - amplitude swapping between images.
```
Usage:   (newimg1, newimg2) = swapphase(img1, img2)

Arguments:
 img1, img2::Array{<:Real,2} - Two images of same size to be used as input

Returns:
    newimg1::Array{Float64,2} - Image obtained from the phase of img1
                                and the magnitude of img2.
    newimg2::Array{Float64,2} - Phase of img2, magnitude of img1.
```

See also: [`quantizephase`](@ref), [`nophase`](@ref)
"""
function swapphase(img1::AbstractArray{T,2}, img2::AbstractArray{T,2}) where T <: Real

    if size(img1) != size(img2)
	error("Images must be the same size")
    end

    # Take FFTs, get magnitude and phase
    IMG1 = fft(img1);         IMG2 = fft(img2)
    mag1 = abs.(IMG1);        mag2 = abs.(IMG2)
    phase1 = IMG1 ./ mag1;    phase2 = IMG2./mag2

    # Now swap amplitude and phase between images.
    NEWIMG1 = mag2.*phase1;   NEWIMG2 = mag1.*phase2
    newimg1 = real.(ifft(NEWIMG1))
    newimg2 = real.(ifft(NEWIMG2))

    # Scale image values 0-1
    newimg1 = imgnormalize(newimg1);
    newimg2 = imgnormalize(newimg2);

    return newimg1, newimg2
end

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
module LocalSearch
using ..JOH

function greedy_flip(sol::JOH.Solution)
    sol(sol.bitlist, sol.problem)

end

"""Performs an exhaustive search of the local bit flip neighborhood, over and
over, until no improvement is found. """
function greedy_flip(sol::Solution, problem::Problem)
    sol = CompleteSolution(sol.bitlist, problem)
    improved = true
    while improved
        improved = greedy_flip_internal!(sol, problem)
    end
    Solution(sol)
end
#
# """Loop over every bit in a bitarry, and calculate the score if the bit were
# flipped. Then, if an improved score were found, flip the bit that lead to the
# greatest improvement, and return true. Else, return false. """
# function greedy_flip_internal!(sol::CompleteSolution, problem::Problem)::Bool
#     index_to_change = 0
#     best_found_score = sol.score
#     # println("best found score is $best_found_score")
#     feas = best_found_score > 0
#     # println("feas is $feas")
#     for i in 1:length(sol.bitlist)
#         # println("starting score is $(sol.score)")
#         if flip_bit!(sol, problem, i, feas=feas)
#             # println("resulting flip scores $(sol.score)")
#             if sol.score > best_found_score
#                 # println("new high found")
#                 best_found_score = sol.score
#                 index_to_change = i
#             else
#                 # println("feas short circuit")
#             end
#
#             flip_bit!(sol, problem, i) #flip the bit back
#         end
#         # println("ending score is $(sol.score)")
#     end
#     if index_to_change > 0
#         # println("changing an index $index_to_change")
#         flip_bit!(sol, problem, index_to_change)
#         # println("score is now $(sol.score)")
#         return true
#     end
#     return false
# end
#
end

ext: jl

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
export ItemKNN

"""
    ItemKNN(
        data::DataAccessor,
        k::Int
    )

[Item-based CF](https://dl.acm.org/citation.cfm?id=963776) that provides a way to model item-item concepts by utilizing the similarities of items in the CF paradigm. `k` represents number of neighbors.

Item properties are relatively stable compared to the users' tastes, and the number of items is generally smaller than the number of users. Hence, while user-based CF successfully captures the similarities of users' complex tastes, modeling item-item concepts could be much more promising in terms of both scalability and overall accuracy.

Item-based CF defines a similarity between an item ``i`` and ``j`` as:

```math
s_{i,j} = \\frac{ \\sum_{u \\in \\mathcal{U}_{i \\cap j}}  (r_{u, i} - \\overline{r}_i) (r_{u, j} - \\overline{r}_j)}
{ \\sqrt{\\sum_{u \\in \\mathcal{U}_{i \\cap j}} (r_{u,i} - \\overline{r}_i)^2} \\sqrt{\\sum_{u \\in \\mathcal{U}_{i \\cap j}} (r_{u, j} - \\overline{r}_j)^2} },
```

where ``\\mathcal{U}_{i \\cap j}`` is a set of users that both of ``r_{u,i}`` and ``r_{u, j}`` are not missing, and ``\\overline{r}_i, \\overline{r}_j`` are mean values of ``i``-th and ``j``-th column in ``R``. Similarly to the user-based algorithm, for the ``t``-th nearest-neighborhood item ``\\tau(t)``, prediction can be done by top-``k`` weighted sum of target user's feedbacks:

```math
r_{u,i} = \\frac{\\sum^k_{t=1} s_{i,\\tau(t)} \\cdot r_{u,\\tau(t)} }{ \\sum^k_{t=1} s_{i,\\tau(t)} }.
```

In case that the number of items is smaller than users, item-based CF could be a more reasonable choice than the user-based approach.
"""
struct ItemKNN <: Recommender
    data::DataAccessor
    k::Int
    sim::AbstractMatrix

    function ItemKNN(data::DataAccessor, k::Int)
        n_item = size(data.R, 2)
        new(data, k, matrix(n_item, n_item))
    end
end

ItemKNN(data::DataAccessor) = ItemKNN(data, 5)

isbuilt(recommender::ItemKNN) = isfilled(recommender.sim)

function build!(recommender::ItemKNN; adjusted_cosine::Bool=false)
    # cosine similarity

    R = copy(recommender.data.R)
    n_row, n_col = size(R)

    if adjusted_cosine
        # subtract mean
        for ri in 1:n_row
            indices = broadcast(!isnan, R[ri, :])
            vmean = mean(R[ri, indices])
            R[ri, indices] .-= vmean
        end
    end

    # unlike pearson correlation, matrix can be filled by zeros for cosine similarity
    R[isnan.(R)] .= 0

    # compute L2 nrom of each column
    norms = sqrt.(sum(R.^2, dims=1))

    for ci in 1:n_col
        for cj in ci:n_col
            numer = dot(R[:, ci], R[:, cj])
            denom = norms[ci] * norms[cj]
            s = numer / denom

            recommender.sim[ci, cj] = s
            if (ci != cj); recommender.sim[cj, ci] = s; end
        end
    end

    # NaN similarities are converted into zeros
    recommender.sim[isnan.(recommender.sim)] .= 0
end

function predict(recommender::ItemKNN, u::Int, i::Int)
    check_build_status(recommender)

    numer = denom = 0

    # negative similarities are filtered
    pairs = collect(zip(1:size(recommender.data.R)[2], max.(recommender.sim[i, :], 0)))
    ordered_pairs = sort(pairs, by=tuple->last(tuple), rev=true)[1:recommender.k]

    for (j, s) in ordered_pairs
        r = recommender.data.R[u, j]
        if isnan(r); continue; end

        numer += s * r
        denom += s
    end

    (denom == 0) ? 0 : numer / denom
end


text:
"""
Control-based continuation: Bifurcation and stability analysis for physical experiments
David A.W. Barton
Mechanical Systems and Signal Processing, Volume 84, Part B, 1 February 2017, Pages 54-64
"""
module _NTMD_MSSP_2017

export NTMD_MSSP_2017

using ..NonlinearSciMLData: NonlinearSciMLData, describe
using DataDeps: register, DataDep, unpack

struct NTMD_MSSP_2017 end

function NonlinearSciMLData.paper_doi(::Type{NTMD_MSSP_2017})
    return "https://doi.org/10.1016/j.ymssp.2015.12.039"
end
function NonlinearSciMLData.data_doi(::Type{NTMD_MSSP_2017})
    return "https://doi.org/10.5523/bris.1esk6h8klq0bg1gdcr78a2ypyz"
end

function NonlinearSciMLData.describe(::Type{NTMD_MSSP_2017})
    return "Experimental data associated with the paper \"Control-based continuation: bifurcation and stability analysis for physical experiments\" by David A.W. Barton"
end

function NonlinearSciMLData.experiments(::Type{NTMD_MSSP_2017})
    return Dict(
        :E1ai =>
            "Stable orbit. Uncontrolled with random inputs. Variety of random amplitudes [0.1:0.1:0.5] x 20 sec x 10 replicates.",
        :E1aii =>
            "Stable orbit. Uncontrolled with random inputs. A long time capture (2 minutes) x2.",
        :E1bi =>
            "Stable orbit. Controlled with random inputs. Variety of random amplitudes [0.1:0.1:0.5] x 20 sec x 10 replicates.",
        :E1bii =>
            "Stable orbit. Controlled with random inputs. A long time capture (2 minutes) x2.",
        :E2bi =>
            "Unstable orbit. Controlled with random inputs. Variety of random amplitudes [0.1:0.1:0.5] x 20 sec x 10 replicates.",
        :E2bii =>
            "Unstable orbit. Controlled with random inputs. A long time capture (2 minutes) x2.",
        :E2ci =>
            "Unstable orbit. Uncontrolled release followed by return to control. 40 replicates.",
        :E3a =>
            "Continuation of a single branch. Variety of random amplitudes [0.1, 0.3, 0.5] x 20 sec.",
        :E3b =>
            "Continuation of a single branch. Several replicates for fixed random amplitude (0.5).",
        :E4 => "Continuation of the full surface.",
    )
end

function get_meta(data)
    timestamp = DateTime(floor.(Int, data["timestamp_created"])...)
    fourier = data["fourier"]
end

function get_observations(data)
    # base is the input
    # x is the output

    # unfortunately it would have been better to have an accelerometer on the base rather than a displacement sensor...
end

function NonlinearSciMLData.experiment(::Type{NTMD_MSSP_2017}, label::Symbol)
    if label === :E1ai
    elseif label === :E1aii
    elseif label === :E1bi
    elseif label === :E1bii
    elseif label === :E2bi
    elseif label === :E2bii
    elseif label === :E2ci
    elseif label === :E2cii
    end
end

function NonlinearSciMLData.register_data(::Type{NTMD_MSSP_2017})
    return register(DataDep(
        string(nameof(NonlinearSciMLData)) * "." * string(nameof(NTMD_MSSP_2017)),
        describe(NTMD_MSSP_2017),
        "http://data.bris.ac.uk/datasets/1esk6h8klq0bg1gdcr78a2ypyz/2015_MSSP_Barton.zip",
        "711bb90e42c4cdbf3c98829108b84489e39801f567868d6096a24361672a40f5";
        post_fetch_method=unpack,
    ))
end

end


text:
nk = (5,5,5)
R  = DFW.Vec3(2, 0, 0)
fermi = 11.4394
n_ωh = 300
n_ωv = 50
ωh = -30.0
ωv = 0.5

job = DFW.DFControl.Job(assetfile(""))
set_magnetization!(job, :Ni1 => [0.0,0.0,1.0], :Ni => [0.0,0.0,-1.0])
hami = DFW.readhami(job)

ω_grid = DFW.setup_ω_grid(ωh, ωv, n_ωh, n_ωv)
kpoints = DFW.ExchangeKGrid(hami, DFW.uniform_kgrid(nk...), R)

@test isapprox(sum(sum.(kpoints.hamiltonian_kgrid.eigvals)), 28310.156014291606)
@test isapprox(sum(kpoints.D),0.07393119748780791 - 1.0312022655727505e-15im)

g_caches = [fill!(similar(kpoints.hamiltonian_kgrid.eigvecs[1]), zero(ComplexF64)) for i=1:3]
G        =fill!(similar(kpoints.hamiltonian_kgrid.eigvecs[1]), zero(ComplexF64))
fill!(G, zero(ComplexF64))
DFW.integrate_Gk!(G, ω_grid[1], fermi, kpoints, g_caches);
exch     = calc_exchanges(hami, atoms(job), fermi; R=R, site_diagonal=false, nk=nk, n_ωh = n_ωh, n_ωv = n_ωv, ωh = ωh, ωv = ωv )
maxJ = maximum([tr(e.J) for e in exch])
@test isapprox(maxJ, 26.044428709929104)

exch1    = calc_exchanges(hami, atoms(job), fermi; R=R, site_diagonal=true, nk=nk, n_ωh = n_ωh, n_ωv = n_ωv, ωh = ωh, ωv = ωv )
maxJ1 = maximum([sum(e.J) for e in exch1])
@test isapprox(maxJ, maxJ1)

hami_soc = [DFW.TbBlock(b.R_cryst, b.R_cart, DFW.NonColinMatrix(DFW.up(b.block), DFW.down(b.block)), DFW.NonColinMatrix(DFW.up(b.tb_block), DFW.down(b.tb_block))) for b in hami]

set_projections!(job.structure, :Ni1 => [:d], :Ni2 => [:d], :O => [:p], soc=true)
exch     = calc_exchanges(hami_soc, atoms(job), fermi; R=R, site_diagonal=false, nk=nk, n_ωh = n_ωh, n_ωv = n_ωv, ωh = ωh, ωv = ωv )
maxJ1 = maximum([tr(e.J) for e in exch])
@test isapprox(maxJ1, maxJ)


text:
"""
    initialize!(cb::CallbackSet,u,t,integrator::DEIntegrator)

Recursively apply `initialize!` and return whether any modified u
"""
function initialize!(cb::CallbackSet,u,t,integrator::DEIntegrator)
  initialize!(u,t,integrator,false,cb.continuous_callbacks...,cb.discrete_callbacks...)
end
initialize!(cb::CallbackSet{Tuple{},Tuple{}},u,t,integrator::DEIntegrator) = false
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback,cs::DECallback...)
  c.initialize(c,u,t,integrator)
  initialize!(u,t,integrator,any_modified || integrator.u_modified,cs...)
end
function initialize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback)
  c.initialize(c,u,t,integrator)
  any_modified || integrator.u_modified
end

"""
    finalize!(cb::CallbackSet,u,t,integrator::DEIntegrator)

Recursively apply `finalize!` and return whether any modified u
"""
function finalize!(cb::CallbackSet,u,t,integrator::DEIntegrator)
  finalize!(u,t,integrator,false,cb.continuous_callbacks...,cb.discrete_callbacks...)
end
finalize!(cb::CallbackSet{Tuple{},Tuple{}},u,t,integrator::DEIntegrator) = false
function finalize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback,cs::DECallback...)
  c.finalize(c,u,t,integrator)
  finalize!(u,t,integrator,any_modified || integrator.u_modified,cs...)
end
function finalize!(u,t,integrator::DEIntegrator,any_modified::Bool,
                     c::DECallback)
  c.finalize(c,u,t,integrator)
  any_modified || integrator.u_modified
end

# Helpers
Base.isempty(cb::CallbackSet) = isempty(cb.continuous_callbacks) && isempty(cb.discrete_callbacks)
Base.isempty(cb::AbstractContinuousCallback) = false
Base.isempty(cb::AbstractDiscreteCallback) = false

has_continuous_callback(cb::DiscreteCallback) = false
has_continuous_callback(cb::ContinuousCallback) = true
has_continuous_callback(cb::VectorContinuousCallback) = true
has_continuous_callback(cb::CallbackSet) = !isempty(cb.continuous_callbacks)
has_continuous_callback(cb::Nothing) = false

#======================================================#
# Callback handling
#======================================================#

function get_tmp(integrator::DEIntegrator, callback)
  _tmp = get_tmp_cache(integrator)
  _tmp === nothing && return nothing
  _cache = first(_tmp)
  if callback.idxs === nothing
    tmp = _cache
  elseif !(callback.idxs isa Number)
    tmp = @view _cache[callback.idxs]
  else
    tmp = nothing
  end
  return tmp
end

function get_condition(integrator::DEIntegrator, callback, abst)
  tmp = get_tmp(integrator, callback)
  ismutable = !(tmp === nothing)
  if abst == integrator.t
    if callback.idxs === nothing
      tmp = integrator.u
    elseif callback.idxs isa Number
      tmp = integrator.u[callback.idxs]
    else
      tmp = @view integrator.u[callback.idxs]
    end
  elseif abst == integrator.tprev
    if callback.idxs === nothing
      tmp = integrator.uprev
    elseif callback.idxs isa Number
      tmp = integrator.uprev[callback.idxs]
    else
      tmp = @view integrator.uprev[callback.idxs]
    end
  else
    if ismutable
      if callback.idxs === nothing
        integrator(tmp,abst,Val{0})
      else
        integrator(tmp,abst,Val{0},idxs=callback.idxs)
      end
    else
      if callback.idxs === nothing
        tmp = integrator(abst,Val{0})
      else
        tmp = integrator(abst,Val{0},idxs=callback.idxs)
      end
    end
    # ismutable && !(callback.idxs isa Number) ? integrator(tmp,abst,Val{0},idxs=callback.idxs) :
    #                                                 tmp = integrator(abst,Val{0},idxs=callback.idxs)
  end
  integrator.sol.destats.ncondition += 1
  if callback isa VectorContinuousCallback
    callback.condition(@view(integrator.callback_cache.tmp_condition[1:callback.len]),tmp,abst,integrator)
    return @view(integrator.callback_cache.tmp_condition[1:callback.len])
  else
    return callback.condition(tmp,abst,integrator)
  end
end

# Use Recursion to find the first callback for type-stability

# Base Case: Only one callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback)
  (find_callback_time(integrator,callback,1)...,1,1)
end

# Starting Case: Compute on the first callback
function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback, args...)
  find_first_continuous_callback(integrator,find_callback_time(integrator,callback,1)...,1,1,args...)
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,
                                        event_occured::Bool,event_idx::Int,idx::Int,counter::Int,
                                        callback2)
  counter += 1 # counter is idx for callback2.
  tmin2,upcrossing2,event_occurred2,event_idx2 = find_callback_time(integrator,callback2,counter)

  if event_occurred2 && (tmin2 < tmin || !event_occured)
    return tmin2,upcrossing2,true,event_idx2,counter,counter
  else
    return tmin,upcrossing,event_occured,event_idx,idx,counter
  end
end

function find_first_continuous_callback(integrator,tmin::Number,upcrossing::Number,event_occured::Bool,event_idx::Int,idx::Int,counter::Int,callback2,args...)
  find_first_continuous_callback(integrator,find_first_continuous_callback(integrator,tmin,upcrossing,event_occured,event_idx,idx,counter,callback2)...,args...)
end

@inline function determine_event_occurance(integrator,callback::VectorContinuousCallback,counter)
  event_occurred = false
  if callback.interp_points!=0
    addsteps!(integrator)
  end

  ts = range(integrator.tprev, stop=integrator.t, length=callback.interp_points)

  #=
  # Faster but can be inaccurate
  if callback.interp_points > 1
    dt = (integrator.t - integrator.tprev) / (callback.interp_points-1)
  else
    dt = integrator.dt
  end
  ts = integrator.tprev:dt:integrator.t
  =#

  interp_index = 0
  # Check if the event occured
  previous_condition = @views(integrator.callback_cache.previous_condition[1:callback.len])

  if callback.idxs === nothing
    callback.condition(previous_condition,integrator.uprev,integrator.tprev,integrator)
  else
    callback.condition(previous_condition,integrator.uprev[callback.idxs],integrator.tprev,integrator)
  end
  integrator.sol.destats.ncondition += 1

  ivec = integrator.vector_event_last_time
  prev_sign = @view(integrator.callback_cache.prev_sign[1:callback.len])
  next_sign = @view(integrator.callback_cache.next_sign[1:callback.len])

  if integrator.event_last_time == counter && minimum(ODE_DEFAULT_NORM(ArrayInterfaceCore.allowed_getindex(previous_condition,ivec),integrator.t)) <= 100ODE_DEFAULT_NORM(integrator.last_event_error,integrator.t)

    # If there was a previous event, utilize the derivative at the start to
    # chose the previous sign. If the derivative is positive at tprev, then
    # we treat `prev_sign` as negetive, and if the derivative is negative then we
    # treat `prev_sign` as positive, regardless of the postiivity/negativity
    # of the true value due to it being =0 sans floating point issues.

    # Only due this if the discontinuity did not move it far away from an event
    # Since near even we use direction instead of location to reset

    if callback.interp_points==0
      addsteps!(integrator)
    end

    # Evaluate condition slightly in future
    abst = integrator.tprev+integrator.dt * callback.repeat_nudge
    tmp_condition = get_condition(integrator, callback, abst)
    @. prev_sign = sign(previous_condition)
    prev_sign[ivec] = sign(tmp_condition[ivec])
  else
      @. prev_sign = sign(previous_condition)
  end

  prev_sign_index = 1
  abst = integrator.t
  next_condition = get_condition(integrator, callback, abst)
  @. next_sign = sign(next_condition)

  event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
  if sum(event_idx) != 0
    event_occurred = true
    interp_index = callback.interp_points
  end

  if callback.interp_points!=0 && !isdiscrete(integrator.alg) && sum(event_idx) != length(event_idx) # Use the interpolants for safety checking
    fallback = true
    for i in 2:length(ts)
      abst = ts[i]
      copyto!(next_sign,get_condition(integrator, callback, abst))
      _event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
      if sum(_event_idx) != 0
        event_occurred = true
        event_idx = _event_idx
        interp_index = i
        fallback = false
        break
      else
        prev_sign_index = i
      end
    end

    if fallback
      # If you get here, then you need to reset the event_idx to the
      # non-interpolated version

      abst = integrator.t
      next_condition = get_condition(integrator, callback, abst)
      @. next_sign = sign(next_condition)
      event_idx = findall_events!(next_sign,callback.affect!,callback.affect_neg!,prev_sign)
      interp_index = callback.interp_points
    end
  end

  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx
end

@inline function determine_event_occurance(integrator,callback::ContinuousCallback,counter)
  event_occurred = false
  if callback.interp_points!=0
    addsteps!(integrator)
  end

  ts = range(integrator.tprev, stop=integrator.t, length=callback.interp_points)

  #=
  # Faster but can be inaccurate
  if callback.interp_points > 1
    dt = (integrator.t - integrator.tprev) / (callback.interp_points-1)
  else
    dt = integrator.dt
  end
  ts = integrator.tprev:dt:integrator.t
  =#

  interp_index = 0

  # Check if the event occured
  if callback.idxs === nothing
    previous_condition = callback.condition(integrator.uprev,integrator.tprev,integrator)
  else
    @views previous_condition = callback.condition(integrator.uprev[callback.idxs],integrator.tprev,integrator)
  end
  integrator.sol.destats.ncondition += 1

  prev_sign = 0.0
  next_sign = 0.0

  if integrator.event_last_time == counter && minimum(ODE_DEFAULT_NORM(previous_condition,integrator.t)) <= 100ODE_DEFAULT_NORM(integrator.last_event_error,integrator.t)

    # If there was a previous event, utilize the derivative at the start to
    # chose the previous sign. If the derivative is positive at tprev, then
    # we treat `prev_sign` as negetive, and if the derivative is negative then we
    # treat `prev_sign` as positive, regardless of the postiivity/negativity
    # of the true value due to it being =0 sans floating point issues.

    # Only due this if the discontinuity did not move it far away from an event
    # Since near even we use direction instead of location to reset

    if callback.interp_points==0
      addsteps!(integrator)
    end

    # Evaluate condition slightly in future
    abst = integrator.tprev+integrator.dt * callback.repeat_nudge
    tmp_condition = get_condition(integrator, callback, abst)
    prev_sign = sign(tmp_condition)
  else
    prev_sign = sign(previous_condition)
  end

  prev_sign_index = 1
  abst = integrator.t
  next_condition = get_condition(integrator, callback, abst)
  next_sign = sign(next_condition)

  if ((prev_sign < 0 && callback.affect! !== nothing) || (prev_sign > 0 && callback.affect_neg! !== nothing)) && prev_sign*next_sign<=0
    event_occurred = true
    interp_index = callback.interp_points
  elseif callback.interp_points!=0 && !isdiscrete(integrator.alg) # Use the interpolants for safety checking
    for i in 2:length(ts)
      abst = ts[i]
      new_sign = get_condition(integrator, callback, abst)
      if ((prev_sign < 0 && callback.affect! !== nothing) || (prev_sign > 0 && callback.affect_neg! !== nothing)) && prev_sign*new_sign<0
        event_occurred = true
        interp_index = i
        break
      else
        prev_sign_index = i
      end
    end
  end
  event_idx = 1

  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx
end

# rough implementation, needs multiple type handling
# always ensures that if r = bisection(f, (x0, x1))
# then either f(nextfloat(r)) == 0 or f(nextfloat(r)) * f(r) < 0
function bisection(f, tup, t_forward::Bool, rootfind::SciMLBase.RootfindOpt, abstol, reltol; maxiters=1000)
  if rootfind == SciMLBase.LeftRootFind
    NonlinearSolve.solve(NonlinearSolve.NonlinearProblem{false}(f, tup), NonlinearSolve.Falsi(), abstol=abstol, reltol=reltol).left
  else
    NonlinearSolve.solve(NonlinearSolve.NonlinearProblem{false}(f, tup), NonlinearSolve.Falsi(), abstol=abstol, reltol=reltol).right
  end
end

"""
findall_events!(next_sign,affect!,affect_neg!,prev_sign)

Modifies `next_sign` to be an array of booleans for if there is a sign change
in the interval between prev_sign and next_sign
"""
function findall_events!(next_sign::Union{Array,SubArray},affect!::F1,affect_neg!::F2,prev_sign::Union{Array,SubArray}) where {F1,F2}
  @inbounds for i in 1:length(prev_sign)
    next_sign[i] = ((prev_sign[i] < 0 && affect! !== nothing) || (prev_sign[i] > 0 && affect_neg! !== nothing)) && prev_sign[i]*next_sign[i]<=0
  end
  next_sign
end

function findall_events!(next_sign,affect!::F1,affect_neg!::F2,prev_sign) where {F1,F2}
  hasaffect::Bool = affect! !== nothing
  hasaffectneg::Bool = affect_neg! !== nothing
  f = (n,p)-> ((p < 0 && hasaffect) || (p > 0 && hasaffectneg)) && p*n<=0
  A = map!(f,next_sign,next_sign,prev_sign)
  next_sign
end

function find_callback_time(integrator,callback::ContinuousCallback,counter)
  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx = determine_event_occurance(integrator,callback,counter)
  if event_occurred
    if callback.condition === nothing
      new_t = zero(typeof(integrator.t))
    else
      if callback.interp_points!=0
        top_t = ts[interp_index] # Top at the smallest
        bottom_t = ts[prev_sign_index]
      else
        top_t = integrator.t
        bottom_t = integrator.tprev
      end
      if callback.rootfind != SciMLBase.NoRootFind && !isdiscrete(integrator.alg)
        zero_func(abst, p=nothing) = get_condition(integrator, callback, abst)
        if zero_func(top_t) == 0
          Θ = top_t
        else
          if integrator.event_last_time == counter &&
            abs(zero_func(bottom_t)) <= 100abs(integrator.last_event_error) &&
            prev_sign_index == 1

            # Determined that there is an event by derivative
            # But floating point error may make the end point negative

            bottom_t += integrator.dt * callback.repeat_nudge
            sign_top = sign(zero_func(top_t))
            sign(zero_func(bottom_t)) * sign_top >= zero(sign_top) && error("Double callback crossing floating pointer reducer errored. Report this issue.")
          end
          Θ = bisection(zero_func, (bottom_t, top_t), isone(integrator.tdir), callback.rootfind, callback.abstol, callback.reltol)
          integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ), Θ)
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = Θ - integrator.tprev
      elseif interp_index != callback.interp_points && !isdiscrete(integrator.alg)
        new_t = ts[interp_index] - integrator.tprev
      else
        # If no solve and no interpolants, just use endpoint
        new_t = integrator.dt
      end
    end
  else
    new_t = zero(typeof(integrator.t))
  end

  new_t,prev_sign,event_occurred,event_idx
end

function find_callback_time(integrator,callback::VectorContinuousCallback,counter)
  event_occurred,interp_index,ts,prev_sign,prev_sign_index,event_idx = determine_event_occurance(integrator,callback,counter)
  if event_occurred
    if callback.condition === nothing
      new_t = zero(typeof(integrator.t))
      min_event_idx = findfirst(isequal(1),event_idx)
    else
      if callback.interp_points!=0
        top_t = ts[interp_index] # Top at the smallest
        bottom_t = ts[prev_sign_index]
      else
        top_t = integrator.t
        bottom_t = integrator.tprev
      end
      if callback.rootfind != SciMLBase.NoRootFind && !isdiscrete(integrator.alg)
        min_t = nextfloat(top_t)
        min_event_idx = -1
        for idx in 1:length(event_idx)
          if ArrayInterfaceCore.allowed_getindex(event_idx,idx) != 0
            zero_func(abst, p=nothing) = ArrayInterfaceCore.allowed_getindex(get_condition(integrator, callback, abst),idx)
            if zero_func(top_t) == 0
              Θ = top_t
            else
              if integrator.event_last_time == counter &&
                integrator.vector_event_last_time == idx &&
                abs(zero_func(bottom_t)) <= 100abs(integrator.last_event_error) &&
                prev_sign_index == 1

                # Determined that there is an event by derivative
                # But floating point error may make the end point negative

                bottom_t += integrator.dt * callback.repeat_nudge
                sign_top = sign(zero_func(top_t))
                sign(zero_func(bottom_t)) * sign_top >= zero(sign_top) && error("Double callback crossing floating pointer reducer errored. Report this issue.")
              end

              Θ = bisection(zero_func, (bottom_t, top_t), isone(integrator.tdir), callback.rootfind, callback.abstol, callback.reltol)
              if integrator.tdir * Θ < integrator.tdir * min_t
                integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ), Θ)
              end
            end
            if integrator.tdir * Θ < integrator.tdir * min_t
              min_event_idx = idx
              min_t = Θ
            end
          end
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = min_t -integrator.tprev
      elseif interp_index != callback.interp_points && !isdiscrete(integrator.alg)
        new_t = ts[interp_index] - integrator.tprev
        min_event_idx = findfirst(isequal(1),event_idx)
      else
        # If no solve and no interpolants, just use endpoint
        new_t = integrator.dt
        min_event_idx = findfirst(isequal(1),event_idx)
      end
    end
  else
    new_t = zero(typeof(integrator.t))
    min_event_idx = 1
  end

  if event_occurred && min_event_idx < 0
    error("Callback handling failed. Please file an issue with code to reproduce.")
  end

  new_t,ArrayInterfaceCore.allowed_getindex(prev_sign,min_event_idx),event_occurred::Bool,min_event_idx::Int
end

function apply_callback!(integrator,callback::Union{ContinuousCallback,VectorContinuousCallback},cb_time,prev_sign,event_idx)

  if isadaptive(integrator)
    set_proposed_dt!(integrator, integrator.tdir * max(nextfloat(integrator.opts.dtmin), integrator.tdir * callback.dtrelax * integrator.dt))
  end

  change_t_via_interpolation!(integrator,integrator.tprev+cb_time)

  # handle saveat
  _, savedexactly = savevalues!(integrator)
  saved_in_cb = true

  @inbounds if callback.save_positions[1]
    # if already saved then skip saving
    savedexactly || savevalues!(integrator,true)
  end

  integrator.u_modified = true

  if prev_sign < 0
    if callback.affect! === nothing
      integrator.u_modified = false
    else
      callback isa VectorContinuousCallback ? callback.affect!(integrator,event_idx) : callback.affect!(integrator)
    end
  elseif prev_sign > 0
    if callback.affect_neg! === nothing
      integrator.u_modified = false
    else
      callback isa VectorContinuousCallback ? callback.affect_neg!(integrator,event_idx) : callback.affect_neg!(integrator)
    end
  end

  if integrator.u_modified
    reeval_internals_due_to_modification!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
    return true,saved_in_cb
  end
  false,saved_in_cb
end

#Base Case: Just one
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback)
  saved_in_cb = false
  if callback.condition(integrator.u,integrator.t,integrator)
    # handle saveat
    _, savedexactly = savevalues!(integrator)
    saved_in_cb = true
    @inbounds if callback.save_positions[1]
      # if already saved then skip saving
      savedexactly || savevalues!(integrator,true)
    end
    integrator.u_modified = true
    callback.affect!(integrator)
    @inbounds if callback.save_positions[2]
      savevalues!(integrator,true)
      saved_in_cb = true
    end
  end
  integrator.sol.destats.ncondition += 1
  integrator.u_modified,saved_in_cb
end

#Starting: Get bool from first and do next
@inline function apply_discrete_callback!(integrator,callback::DiscreteCallback,args...)
  apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback,args...)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,apply_discrete_callback!(integrator,callback)...,args...)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

@inline function apply_discrete_callback!(integrator,discrete_modified::Bool,saved_in_cb::Bool,callback::DiscreteCallback)
  bool,saved_in_cb2 = apply_discrete_callback!(integrator,callback)
  discrete_modified || bool, saved_in_cb || saved_in_cb2
end

max_vector_callback_length_int(cs::CallbackSet) = max_vector_callback_length_int(cs.continuous_callbacks...)
max_vector_callback_length_int() = nothing
function max_vector_callback_length_int(continuous_callbacks...)
  all(cb->cb isa ContinuousCallback,continuous_callbacks) && return nothing
  maxlen = -1
  for cb in continuous_callbacks
    if cb isa VectorContinuousCallback && cb.len > maxlen
      maxlen = cb.len
    end
  end
  maxlen
end

function max_vector_callback_length(cs::CallbackSet)
  continuous_callbacks = cs.continuous_callbacks
  maxlen_cb = nothing
  maxlen = -1
  for cb in continuous_callbacks
    if cb isa VectorContinuousCallback && cb.len > maxlen
      maxlen = cb.len
      maxlen_cb = cb
    end
  end
  maxlen_cb
end

"""
$(TYPEDEF)
"""
mutable struct CallbackCache{conditionType,signType}
  tmp_condition::conditionType
  previous_condition::conditionType
  next_sign::signType
  prev_sign::signType
end

function CallbackCache(u,max_len,::Type{conditionType},::Type{signType}) where {conditionType,signType}
    tmp_condition = similar(u, conditionType, max_len)
    previous_condition = similar(u, conditionType, max_len)
    next_sign = similar(u, signType, max_len)
    prev_sign = similar(u, signType, max_len)
    CallbackCache(tmp_condition,previous_condition,next_sign,prev_sign)
end

function CallbackCache(max_len,::Type{conditionType},::Type{signType}) where {conditionType,signType}
    tmp_condition = zeros(conditionType, max_len)
    previous_condition = zeros(conditionType, max_len)
    next_sign = zeros(signType, max_len)
    prev_sign = zeros(signType, max_len)
    CallbackCache(tmp_condition,previous_condition,next_sign,prev_sign)
end


text:
# Interface methods

rasterlayers(::Type{EarthEnv{HabitatHeterogeneity}}) =
    (:cv, :evenness, :range, :shannon, :simpson, :std, :Contrast, :Correlation,
     :Dissimilarity, :Entropy, :Homogeneity, :Maximum, :Uniformity, :Variance)

rastername(::Type{EarthEnv{HabitatHeterogeneity}}, layer::Symbol, resolution::Integer) =
    "$(layer)_$(resolution)km.tif"

rasterpath(T::Type{<:EarthEnv{HabitatHeterogeneity}}, layer::Symbol, resolution::Integer) =
    joinpath(rasterpath(T), string(resolution) * "km", rastername(T, layer, resolution))

rasterurl(T::Type{EarthEnv{HabitatHeterogeneity}}, layer, resolution::Integer) =
    joinpath(rasterurl(T), "$(layer)_01_05_$(string(resolution))km_$(_getprecision(layer, resolution)).tif")

function download_raster(T::Type{EarthEnv{HabitatHeterogeneity}}; layer::Symbol=:cv, resolution::Integer=25)
    resolution in [1, 5, 25] || throw(ArgumentError("The resolution must be 1km, 5km, or 25km"))
    layer in rasterlayers(T) || throw(ArgumentError("The layer must be one of $(rasterlayers(T))"))

    path = rasterpath(T, layer, resolution)
    url = rasterurl(T, layer, resolution)
    return _maybe_download(url, path)
end

# Utility methods

_pathsegment(::Type{HabitatHeterogeneity}) = "habitat_heterogeneity"

function _getprecision(layer, resolution)
    precision = "uint16"
    # WELLCOME TO HELL
    precision = ((resolution == 10) && (layer == :cv)) ? "uint32" : precision
    precision = (layer == :Contrast) ? "uint32" : precision
    precision = (layer == :Dissimilarity) ? "uint32" : precision
    precision = (layer == :Variance) ? "uint32" : precision
    precision = ((resolution == 25) && (layer == :Entropy)) ? "uint32" : precision
    # /HELL
end


text:
# This file was generated, do not modify it. # hide
X_oos = MLJ.table(hcat(y1_oos, y2_oos))
m_judge = machine(judge, X_oos, y)


text:
"""
The `Modular` type is used by a lot of functions internally.

It has two fields:
- `A` is the adjacency matrix, of size *R* x *C*
- `S` is the community partition, of size (*R+C*) x *c*, where *c* is the number of modules
"""
type Modular
   A::Array{Int64, 2}
   S::Array{Int64, 2}
end


text:
function trace!(tr, state, iteration, method::NelderMead, options)
    dt = Dict()
    if options.extended_trace
        dt["centroid"] = state.x_centroid
        dt["step_type"] = state.step_type
    end
    update!(tr,
    iteration,
    state.f_lowest,
    state.nm_x,
    dt,
    options.store_trace,
    options.show_trace,
    options.show_every,
    options.callback)
end

function trace!(tr, state, iteration, method::Union{ParticleSwarm, SimulatedAnnealing}, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
    end
    g_norm = NaN
    update!(tr,
            state.iteration,
            state.f_x,
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

function trace!(tr, state, iteration, method::BFGS, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(state.g)
        dt["~inv(H)"] = copy(state.invH)
        dt["Current step size"] = state.alpha
    end
    g_norm = vecnorm(state.g, Inf)
    update!(tr,
    iteration,
    state.f_x,
    g_norm,
    dt,
    options.store_trace,
    options.show_trace,
    options.show_every,
    options.callback)
end

function trace!(tr, state, iteration, method::Union{LBFGS, AcceleratedGradientDescent, GradientDescent, MomentumGradientDescent, ConjugateGradient}, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(state.g)
        dt["Current step size"] = state.alpha
    end
    g_norm = vecnorm(state.g, Inf)
    update!(tr,
            iteration,
            state.f_x,
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

function trace!(tr, state, iteration, method::Newton, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(state.g)
        dt["h(x)"] = copy(state.H)
    end
    g_norm = vecnorm(state.g, Inf)
    update!(tr,
            iteration,
            state.f_x,
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

function trace!(tr, state, iteration, method::NewtonTrustRegion, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(state.g)
        dt["h(x)"] = copy(state.H)
        dt["delta"] = copy(state.delta)
        dt["interior"] = state.interior
        dt["hard case"] = state.hard_case
        dt["reached_subproblem_solution"] = state.reached_subproblem_solution
        dt["lambda"] = state.lambda
    end
    g_norm = norm(state.g, Inf)
    update!(tr,
            iteration,
            state.f_x,
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end


text:
const settings_file = joinpath((@__FILE__) |> dirname |> dirname, "online_package.json")

export set_up
"""
    set_up(username, github_token; get_travis_token(github_token), ssh_keygen_file = "ssh-keygen")

set up `OnlinePackage`.

get a `github_token` [here](https://github.com/settings/tokens/new). make
sure to check the `"public_repo"` scope.

the default `ssh_keygen_file` assumes ssh-keygen is in your path. if not,
it often comes prepacked with git; check `PATH_TO_GIT/usr/bin/ssh-keygen"`.
"""
set_up(username, github_token; ssh_keygen_file = "ssh-keygen", travis_token = get_travis_token(github_token)) =
    open(settings_file, "w") do io
        JSON.print(io, Dict(
            "username" => username,
            "github_token" => github_token,
            "travis_token" => travis_token,
            "ssh_keygen_file" => ssh_keygen_file))
    end

function settings(name; settings_file = settings_file)
    if !ispath(settings_file)
        error("Cannot find settings. Please `set_up`")
    end
    dict = JSON.parsefile(settings_file)
    if !haskey(dict, name)
        error("Missing setting $name")
    end
    dict[name]
end


text:
"""
    Before calling this function one should check that `startpos` is inside `det`: `in(startpos, det`

"""
function drift_charge!(
                            drift_path::Vector{CartesianPoint{T}},
                            det::SolidStateDetector{T, :cartesian},
                            point_types::PointTypes{T, 3, :cartesian},
                            grid::Grid{T, 3, :cartesian},
                            startpos::CartesianPoint{T},
                            Δt::T,
                            velocity_field::Interpolations.Extrapolation{<:StaticVector{3}, 3};
                            verbose::Bool = true
                        )::Nothing where {T <: SSDFloat}
    done::Bool = false
    drift_path[1] = startpos
    null_step::CartesianVector{T} = CartesianVector{T}(0, 0, 0)
    @inbounds for istep in eachindex(drift_path)[2:end] #end]
        if done == false
            current_pos::CartesianPoint{T} = drift_path[istep - 1]
            stepvector::CartesianVector{T} = get_velocity_vector(velocity_field, current_pos) * Δt
            if stepvector == null_step done = true end
            next_pos::CartesianPoint{T} = current_pos + stepvector
            if next_pos in point_types
                drift_path[istep] = next_pos
            elseif next_pos in det
                drift_path[istep] = next_pos
            else
                crossing_pos::CartesianPoint{T}, cd_point_type::UInt8, boundary_index::Int, surface_normal::CartesianVector{T} = get_crossing_pos(det, grid, current_pos, next_pos)
                if cd_point_type == CD_ELECTRODE
                    done = true
                    drift_path[istep] = crossing_pos
                elseif cd_point_type == CD_FLOATING_BOUNDARY
                    projected_vector = project_to_plane(stepvector, surface_normal )
                    next_pos = current_pos + projected_vector
                    # ToDo: We actually need a time array as well to do this properly...
                    small_projected_vector = projected_vector * T(0.001)
                    i::Int = 0
                    while i < 1000 && !(next_pos in det)
                        next_pos -= small_projected_vector
                        i += 1
                    end
                    if i == 1000 && verbose @warn("Handling of charge at floating boundary did not work as intended. Start Position (Cart): $startpos") end
                    drift_path[istep] = next_pos
                elseif cd_point_type == CD_BULK
                    if verbose @warn ("Internal error for charge starting at $startpos") end
                    drift_path[istep] = current_pos
                    done = true
                else # elseif cd_point_type == CD_OUTSIDE
                    if verbose @warn ("Internal error for charge starting at $startpos") end
                    drift_path[istep] = current_pos
                    done = true
                end
            end
        elseif done == true
            drift_path[istep] = drift_path[istep-1]
        end
    end
    return nothing
end

# Point types for charge drift: Defined in DetectorGeometries/DetectorGeometries.jl
# const CD_ELECTRODE = 0x00
# const CD_OUTSIDE = 0x01
# const CD_BULK = 0x02
# const CD_FLOATING_BOUNDARY = 0x04 # not 0x03, so that one could use bit operations here...

function get_crossing_pos(  detector::SolidStateDetector{T, :cartesian}, grid::Grid{T, 3}, point_in::CartesianPoint{T}, point_out::CartesianPoint{T};
                            max_n_iter::Int = 2000)::Tuple{CartesianPoint{T}, UInt8, Int, CartesianVector{T}} where {T <: SSDFloat}
    point_mid::CartesianPoint{T} = T(0.5) * (point_in + point_out)
    cd_point_type::UInt8, contact_idx::Int, surface_normal::CartesianVector{T} = point_type(detector, grid, point_mid)
    for i in 1:max_n_iter
        if cd_point_type == CD_BULK
            point_in = point_mid
        elseif cd_point_type == CD_OUTSIDE
            point_out = point_mid
        elseif cd_point_type == CD_ELECTRODE
            break
        else #elseif cd_point_type == CD_FLOATING_BOUNDARY
            break
        end
        point_mid = T(0.5) * (point_in + point_out)
        cd_point_type, contact_idx, surface_normal = point_type(detector, grid, point_mid)
    end
    return point_mid, cd_point_type, contact_idx, surface_normal
end


text:
# The transformer below, not intended for the general user, wraps a
# single data-manipulating function without parameters, and is used by
# the @pipeline macro.  Generally, the only time the user will need to
# implement such a tranformer is when it also depends on parameters,
# and the recommendation is to instead subtype `Static` with a new
# custom type. In the context of a learning network, the user can use
# the `node` function to overload data-manipulating operations for
# nodes (if they are not already overloaded out-of-the-box).

const STATIC_TRANSFORMER_DESCR = "Applies a given data transformation
`f` (either a function or callable)."

"""
$STATIC_TRANSFORMER_DESCR

## Field

* `f=identity`: function or callable object to use for the data transformation.
"""
mutable struct WrappedFunction <: MLJBase.Static
    f
end
WrappedFunction(;f=identity) = WrappedFunction(f)

fitted_params(::WrappedFunction) = NamedTuple()
fit(::WrappedFunction, ::Integer, _) = nothing, nothing, NamedTuple()
transform(model::WrappedFunction, fitresult, Xnew) = (model.f)(Xnew)

MMI.is_wrapper(::Type{<:WrappedFunction}) = true
MMI.is_pure_julia(::Type{<:WrappedFunction}) = true
MMI.load_path(::Type{<:WrappedFunction}) = "MLJBase.WrappedFunction"
MMI.package_name(::Type{<:WrappedFunction}) = "MLJBase"


text:
using Random

function selectrandom(array, fitness)
    # fitness is not necessarly
    # 引数を合わせるため
    array = shuffle(array)
    return array
end

function selectroulette(array, fitness)
    populationSize = size(array, 1)
    geneSize = size(array, 2)
    prob = fitness./ sum(fitness)

    selectedIndex = pselection(prob, populationSize)
    # Intに限定
    selected = reshape(zeros(Int, length(array)), (populationSize, geneSize))
    for (i, j) in enumerate(selectedIndex)
        selected[i, :] = array[j, :]
    end
    return selected
end

# probability selection
function pselection(prob::Vector{Float64}, N::Int)
    cp = cumsum(prob)
    selected = Array{Int}(undef, N)
    for i in 1:N
        j = 1
        r = rand()
        while cp[j] < r
            j += 1
        end
        selected[i] = j
    end
    return selected
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
module NumericFuns

    using Compat

    export

    # mathfuns
    sqr, rcp, rsqrt, rcbrt, xlogx, xlogy,
    sigmoid, logistic, logit, softplus, invsoftplus, logsumexp,

    # functors (Note: export of specific functors are in functors.jl)
    Functor, UnaryFunctor, BinaryFunctor, TernaryFunctor,
    evaluate, @functor1, @functor2,

    # rtypes
    fptype, arithtype, result_type

    # source files

    include("mathfuns.jl")
    include("functors.jl")
    include("rtypes.jl")

end # module


text:
using LinearAlgebraicRepresentation
using ViewerGL, LinearAlgebra, SparseArrays, DataStructures, NearestNeighbors
L = Lar = LinearAlgebraicRepresentation
GL = ViewerGL

# initialization
n = 3
v,(_,ev,fv) = L.cuboidGrid([n,n],true)
v = [v; zeros(1,(n+1)*(n+1))]
g = (v,ev,fv)

# range di numeri casuali piccoli 0e10-2  x  1e10-2
interval =  0e10-2 : 1e10-2

# matrice 4x4 di numeri casuali piccoli (dipende dagli zeri)
A(interval) = (1_000_000 ./ rand(interval,4,4))

# Esempio di complesso 2D immerso in 3D
G1 = L.Struct([  g, L.r(pi/2,0,0) + A(interval),  g , L.r(0,pi/2,0), g , L.r(0,0,pi/2), L.s(-1,-1,-1), g ,
		L.t(2,0,0), L.r(pi/2,0,0),  g , L.t(0,0,2), L.r(0,pi/2,0), L.s(-1,-1,-1), g ]);
V1,EV1,FV1 = Lar.struct2lar(G1);
VV1 = [[k] for  k=1:size(V1,2)];
#GL.VIEW(push!(GL.numbering(.5)((V1,[VV1,EV1,FV1]), GL.COLORS[1], 0.1),GL.GLFrame2));

# Esempio di complesso 2D immerso in 3D
G2 = L.Struct([  g, L.r(pi/2,0,0) + A(interval)./10, g , L.r(0,pi/2,0), g , L.r(0,0,pi/2), g ,
		L.t(2,0,0), L.r(pi/2,0,0),  g , L.t(0,0,2), L.r(0,pi/2,0), g ]);
V2,EV2,FV2 = Lar.struct2lar(G2);
VV2 = [[k] for  k=1:size(V2,2)];
#GL.VIEW(push!(GL.numbering(.5)((V2,[VV2,EV2,FV2]), GL.COLORS[1], 0.1),GL.GLFrame2));

# complesso unione combinatoria di G1 e G2
G = L.Struct([ G1, G2 ])
V,EV,FV = Lar.struct2lar(G)
VV = [[k] for  k=1:size(V,2)]
GL.VIEW(push!(GL.numbering(.5)((V,[VV,EV,FV]), GL.COLORS[1], 0.1),GL.GLFrame2));
Delta_0 = Lar.coboundary_0(EV::Lar.Cells)
Delta_1 = Lar.coboundary_1(FV::Lar.Cells, EV::Lar.Cells)

W,EW,FW,FE = Lar.chaincongruence(V, Delta_0, Delta_1, epsilon=0.1);

WW = [[k] for  k=1:size(W,2)]
GL.VIEW(push!(GL.numbering(.5)((Lar.Points(W),Lar.Cells[WW,EW,FW]), GL.COLORS[1], 0.1),GL.GLFrame2));


text:
################################################################################
# Residual
################################################################################
function residual!(prob::GameProblem{KN,n,m,T,SVd,SVx}) where {KN,n,m,T,SVd,SVx}
	residual!(prob, prob.pdtraj)
	return nothing
end

function residual!(prob::GameProblem{KN,n,m,T,SVd,SVx}, pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}
	N = prob.probsize.N
	p = prob.probsize.p
    pu = prob.probsize.pu
	core = prob.core
	model = prob.model
	game_obj = prob.game_obj

	# Initialization
	prob.core.res .= 0.0
    stamp = VStamp()
	∇dyn = zeros(MMatrix{n,(n+m),T,n*(n+m)})

	# Cost
	cost_gradient!(game_obj, pdtraj)
	for i = 1:p
		# State cost
		for k = 1:N
			stampify!(stamp, :opt, i, :x, 1, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].q) : nothing
			end
		end
		# Control Cost
		for k = 1:N-1
			stampify!(stamp, :opt, i, :u, i, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].r[pu[i]]) : nothing
			end
		end
	end
	# Dynamics penalty
	for k = 1:N-1
		∇dynamics!(∇dyn, model, pdtraj, k)
		for i = 1:p
			λik = pdtraj.du[i][k]
			stampify!(stamp, :opt, i, :x, 1, k)
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:x][1]]'*λik) : nothing
			stampify!(stamp, :opt, i, :u, i, k)
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:u][i]]'*λik) : nothing
			stampify!(stamp, :opt, i, :x, 1, k+1)
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], -λik) : nothing
		end
	end

	# Constraints
	constraint_residual!(prob, pdtraj)

	# Dynamics
    for k = 1:N-1
        stampify!(stamp, :dyn, 1, :x, 1, k)
        valid(stamp, N, p) ? add2sub(core.res_sub[stamp], dynamics_residual(model, pdtraj, k)) : nothing # shouldn't allocate
    end
    return nothing
end

function regularize_residual!(core::NewtonCore, opts::Options, pdtraj::PrimalDualTraj, pdtraj_ref::PrimalDualTraj)
	N = core.probsize.N
	p = core.probsize.p
	pu = core.probsize.pu

	stamp = VStamp()
	for k = 1:N-1
		x = state(pdtraj.pr[k+1])
		x_ref = state(pdtraj_ref.pr[k+1])
		u = control(pdtraj.pr[k])
		u_ref = control(pdtraj_ref.pr[k])
		for i = 1:p
			stampify!(stamp, :opt, i, :x, 1, k+1)
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], opts.reg.x*(x - x_ref)) : nothing
			stampify!(stamp, :opt, i, :u, i, k)
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], opts.reg.u*(u[pu[i]] - u_ref[pu[i]])) : nothing
		end
	end
	return nothing
end

function residual_norm(prob::GameProblem{KN,n,m,T,SVd,SVx}) where {KN,n,m,T,SVd,SVx}
	residual_norm(prob, prob.pdtraj)
end

function residual_norm(prob::GameProblem{KN,n,m,T,SVd,SVx},
		pdtraj::PrimalDualTraj{KN,n,m,T,SVd}; update_residual::Bool=true) where {KN,n,m,T,SVd,SVx}
	update_residual && residual!(prob, pdtraj)
	res_norm = norm(prob.core.res, 1)/length(prob.core.res)
	return res_norm
end

################################################################################
# Residual Jacobian
################################################################################

function residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx}) where {KN,n,m,T,SVd,SVx}
	residual_jacobian!(prob, prob.pdtraj)
	return nothing
end

function residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}
    N = prob.probsize.N
	p = prob.probsize.p
    pu = prob.probsize.pu
	model = prob.model
	core = prob.core
    game_obj = prob.game_obj

	# Reset!
	sparse_zero!(prob.core.jac)

    # Allocations
    stamp = Stamp()
    ∇dyn = zeros(MMatrix{n,(n+m),T,n*(n+m)})

	# Cost function
	cost_hessian!(game_obj, pdtraj)
	# Cost
	for i = 1:p
		# State cost
		for k = 1:N
			stampify!(stamp, :opt, i, :x, 1, k, :x, 1, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], game_obj.E[i][j].cost[k].Q) : nothing
			end
		end
		# Control cost
		for k = 1:N-1
			stampify!(stamp, :opt, i, :u, i, k, :u, i, k)
			n_obj = length(game_obj.E[i])
			for j = 1:n_obj
				valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], game_obj.E[i][j].cost[k].R[pu[i],pu[i]]) : nothing
			end
		end
	end

	# Constraints
	constraint_jacobian_residual!(prob, pdtraj)

	# Dynamics
	for k = 1:N-1
        ∇dynamics!(∇dyn, model, pdtraj, k)
        # Bottom Left
        stampify!(stamp, :dyn, 1, :x, 1, k, :x, 1, k)
		∇dyn_x = ∇dyn[:,core.dyn[:x][1]]
		valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:x][1]]) : nothing
        for i = 1:p
            stampify!(stamp, :dyn, 1, :x, 1, k, :u, i, k)
            valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:u][i]]) : nothing
        end
        stampify!(stamp, :dyn, 1, :x, 1, k, :x, 1, k+1)
        valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], -1.0) : nothing
        # Top Right
        for i = 1:p
            stampify!(stamp, :opt, i, :x, 1, k, :λ, i, k)
            valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:x][1]]') : nothing
            stampify!(stamp, :opt, i, :u, i, k, :λ, i, k)
            valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:u][i]]') : nothing
            stampify!(stamp, :opt, i, :x, 1, k+1, :λ, i, k)
            valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], -1.0) : nothing
        end
    end
    return nothing
end

function regularize_residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx}) where {KN,n,m,T,SVd,SVx}
	N = prob.probsize.N
	p = prob.probsize.p
	pu = prob.probsize.pu
	core = prob.core
	opts = prob.opts

	stamp = Stamp()
	for k = 1:N-1
		for i = 1:p
			stampify!(stamp, :opt, i, :x, 1, k+1, :x, 1, k+1)
			valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], opts.reg.x) : nothing
			stampify!(stamp, :opt, i, :u, i, k, :u, i, k)
			valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], opts.reg.u) : nothing
		end
	end
	return nothing
end

################################################################################
# Iterative Best Response Residual
################################################################################
function ibr_residual!(prob::GameProblem{KN,n,m,T,SVd,SVx}, i::Int) where {KN,n,m,T,SVd,SVx}
	ibr_residual!(prob, prob.pdtraj, i)
	return nothing
end

function ibr_residual!(prob::GameProblem{KN,n,m,T,SVd,SVx}, pdtraj::PrimalDualTraj{KN,n,m,T,SVd}, i::Int) where {KN,n,m,T,SVd,SVx}
	N = prob.probsize.N
	p = prob.probsize.p
    pu = prob.probsize.pu
	core = prob.core
	model = prob.model
	game_obj = prob.game_obj

	# Initialization
	prob.core.res .= 0.0
    stamp = VStamp()
	∇dyn = zeros(MMatrix{n,(n+m),T,n*(n+m)})

	# Cost
	cost_gradient!(game_obj, pdtraj, i)
	# State cost
	for k = 1:N
		stampify!(stamp, :opt, i, :x, 1, k)
		n_obj = length(game_obj.E[i])
		for j = 1:n_obj
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].q) : nothing
		end
	end
	# Control Cost
	for k = 1:N-1
		stampify!(stamp, :opt, i, :u, i, k)
		n_obj = length(game_obj.E[i])
		for j = 1:n_obj
			valid(stamp, N, p) ? add2sub(core.res_sub[stamp], game_obj.E[i][j].cost[k].r[pu[i]]) : nothing
		end
	end
	# Dynamics penalty
	for k = 1:N-1
		∇dynamics!(∇dyn, model, pdtraj, k)
		λik = pdtraj.du[i][k]
		stampify!(stamp, :opt, i, :x, 1, k)
		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:x][1]]'*λik) : nothing
		stampify!(stamp, :opt, i, :u, i, k)
		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], ∇dyn[:,core.dyn[:u][i]]'*λik) : nothing
		stampify!(stamp, :opt, i, :x, 1, k+1)
		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], -λik) : nothing
	end

	# Constraints
	constraint_residual!(prob, pdtraj)

	# Dynamics
    for k = 1:N-1
        stampify!(stamp, :dyn, 1, :x, 1, k)
        valid(stamp, N, p) ? add2sub(core.res_sub[stamp], dynamics_residual(model, pdtraj, k)) : nothing # shouldn't allocate
    end
    return nothing
end

function regularize_ibr_residual!(core::NewtonCore, opts::Options, pdtraj::PrimalDualTraj, pdtraj_ref::PrimalDualTraj, i::Int)
	N = core.probsize.N
	p = core.probsize.p
	pu = core.probsize.pu

	stamp = VStamp()
	for k = 1:N-1
		x = state(pdtraj.pr[k+1])
		x_ref = state(pdtraj_ref.pr[k+1])
		u = control(pdtraj.pr[k])
		u_ref = control(pdtraj_ref.pr[k])
		stampify!(stamp, :opt, i, :x, 1, k+1)
		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], opts.reg.x*(x - x_ref)) : nothing
		stampify!(stamp, :opt, i, :u, i, k)
		valid(stamp, N, p) ? add2sub(core.res_sub[stamp], opts.reg.u*(u[pu[i]] - u_ref[pu[i]])) : nothing
	end
	return nothing
end

################################################################################
# Iterative Best Response Residual Jacobian
################################################################################

function ibr_residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx}, i::Int) where {KN,n,m,T,SVd,SVx}
	ibr_residual_jacobian!(prob, prob.pdtraj, i)
	return nothing
end

function ibr_residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}, i::Int) where {KN,n,m,T,SVd,SVx}
    N = prob.probsize.N
	p = prob.probsize.p
    pu = prob.probsize.pu
	model = prob.model
	core = prob.core
    game_obj = prob.game_obj

	# Reset!
	sparse_zero!(prob.core.jac)

    # Allocations
    stamp = Stamp()
    ∇dyn = zeros(MMatrix{n,(n+m),T,n*(n+m)})

	# Cost function
	cost_hessian!(game_obj, pdtraj, i)
	# Cost
	# State cost
	for k = 1:N
		stampify!(stamp, :opt, i, :x, 1, k, :x, 1, k)
		n_obj = length(game_obj.E[i])
		for j = 1:n_obj
			valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], game_obj.E[i][j].cost[k].Q) : nothing
		end
	end
	# Control cost
	for k = 1:N-1
		stampify!(stamp, :opt, i, :u, i, k, :u, i, k)
		n_obj = length(game_obj.E[i])
		for j = 1:n_obj
			valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], game_obj.E[i][j].cost[k].R[pu[i],pu[i]]) : nothing
		end
	end

	# Constraints
	constraint_jacobian_residual!(prob, pdtraj)

	# Dynamics
	for k = 1:N-1
        ∇dynamics!(∇dyn, model, pdtraj, k)
        # Bottom Left
        stampify!(stamp, :dyn, 1, :x, 1, k, :x, 1, k)
		∇dyn_x = ∇dyn[:,core.dyn[:x][1]]
		valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:x][1]]) : nothing

        stampify!(stamp, :dyn, 1, :x, 1, k, :u, i, k)
        valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:u][i]]) : nothing

        stampify!(stamp, :dyn, 1, :x, 1, k, :x, 1, k+1)
        valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], -1.0) : nothing
        # Top Right
        stampify!(stamp, :opt, i, :x, 1, k, :λ, i, k)
        valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:x][1]]') : nothing
        stampify!(stamp, :opt, i, :u, i, k, :λ, i, k)
        valid(stamp, N, p) ? add2sub(core.jac_sub[stamp], ∇dyn[:,core.dyn[:u][i]]') : nothing
        stampify!(stamp, :opt, i, :x, 1, k+1, :λ, i, k)
        valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], -1.0) : nothing
    end
    return nothing
end

function regularize_ibr_residual_jacobian!(prob::GameProblem{KN,n,m,T,SVd,SVx}, i::Int) where {KN,n,m,T,SVd,SVx}
	N = prob.probsize.N
	p = prob.probsize.p
	pu = prob.probsize.pu
	core = prob.core
	opts = prob.opts

	stamp = Stamp()
	for k = 1:N-1
		stampify!(stamp, :opt, i, :x, 1, k+1, :x, 1, k+1)
		valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], opts.reg.x) : nothing
		stampify!(stamp, :opt, i, :u, i, k, :u, i, k)
		valid(stamp, N, p) ? addI2sub(core.jac_sub[stamp], opts.reg.u) : nothing
	end
	return nothing
end


text:
global OpenGLver="3.2"

module Processing

using OpenGL
using SDL

#include("constants.jl")

export @setup, @draw
export displayHeight, displayWidth, height, winSize, width
# export cursor focused, frameCount, frameRate, noCursor
# export createShape, loadShape
export arc, ellipse, line, point, quad, rect, triangle
# export bezier, bezierDetail, bezierPoint, bezierTangent, curve, curveDetail, curvePoint, curveTangent, curveTightness
export sphere, sphereDetail
# export box
# export ellipseMode, noSmooth, rectMode, smooth, strokeCap, strokeJoin, strokeWeight
export beginShape, endShape, vertex
# export bezierVertex, curveVertex, quadraticVertex
# export shape, shapeMode
export background, fill, noFill
# export colorMode, noStroke, Stroke
# export alpha, blue, brightness, color, green, hue, lerpColor, red, saturation
# export createImage
# export image, imageMode, loadImage, noTint, requestImage, tint
# export texture, textureMode, textureWrap
# export blend, copy, filter, get, loadPixels, set, updatePixels
export blendMode
# export createGraphics, hint
# export loadShader, resetShader, shader
# export createFont, loadFont, text, textFont
# export textAlign, textLeading, textMode, textSize, textWidth,
# export textAscent, textDescent

# state structures

type stateStruct
	Sdet::Real #sphere detail
	noFill::Bool #should drawing elements be filled or not?
	h::Integer #height of display window
	w::Integer #width of display window
	dispH::Integer #height of display area
	dispW::Integer #width of display area
	winTitle::String #title of display window
	iconTitle::String #title of docked icon
	bpp::Integer #Bits Per Pixel (BPP)

	function stateStruct()
		new(
			5, #give spheres decent detail, but be kind to low-end systems
			1, #don't fill by default
			480, #Height - default = 480
			640, #Width - default = 640
			480, #Display height - default = 480
			640, #Display width - default = 640
			"Processing.jl",
			"Processing.jl",
			16 #BPP - default = 16bpp
			)
	end
end

# macros to simulate Processing environment

macro setup(body)
	processingState = stateStruct()
	println(processingState)

	SDL_Init(SDL_INIT_VIDEO)
	videoFlags = (SDL_OPENGL | SDL_GL_DOUBLEBUFFER | SDL_HWPALETTE | SDL_RESIZABLE | SDL_HWSURFACE | SDL_HWACCEL)
	SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)
	SDL_SetVideoMode(processingState.w, processingState.h, processingState.bpp, videoFlags)
	SDL_WM_SetCaption(processingState.winTitle, processingState.iconTitle)

	glViewport(0, 0, processingState.dispW, processingState.dispH)
	glClearColor(0.0, 0.0, 0.0, 0.0)
	glClearDepth(1.0)
	glDepthFunc(GL_LESS)
	glEnable(GL_DEPTH_TEST)
	glShadeModel(GL_SMOOTH)

	glMatrixMode(GL_PROJECTION)
	glLoadIdentity()

	gluPerspective(45.0,processingState.w/processingState.h,0.1,100.0)

	glMatrixMode(GL_MODELVIEW)

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
	glLoadIdentity()

	body

	SDL_GL_SwapBuffers()
end

macro draw(body)
	while true
		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
		glLoadIdentity()

		body

		SDL_GL_SwapBuffers()
	end
end

# Environment

#cursor

function displayHeight(howHigh)
	processingState.dispH = howHigh
end

function displayWidth(howWide)
	processingState.dispW = howWide
end

#focused
#frameCount
#frameRate

function height(howHigh)
	processingState.h = howHigh
end

#noCursor

function winSize(howWide, howHigh)
	width(howWide)
	height(howHigh)
end

function width(howWide)
	processingState.w = howWide
end

# Shape

#createShape
#loadShape

## 2D Primitives

function arc()
    gluBeginCurve()
    gluEndCurve()
end

function ellipse(xcent,ycent,Mrad,mrad)
    shapeList = glGenLists(1)
    glNewList(shapeList, GL_COMPILE)
        quad = gluNewQuadric()
        gluQuadricDrawStyle(quad, GLU_FILL)
        glTranslate(xcent,ycent)
        gluDisk(quad,mrad,Mrad,10,10)
        glTranslate(-xcent,-ycent)
        gluDeleteQuadric(quad)
    glEndList()
end

function line(x1,y1,x2,y2)
    glBegin(GL_LINE)
        glVertex(x1,y1)
        glVertex(x2,y2)
    glEnd()
end

function point(x,y)
    glBegin(GL_POINTS)
        glVertex(x,y)
    glEnd()
end

function quad(x1,y1,x2,y2,x3,y3,x4,y4)
    glBegin(GL_QUAD)
        glVertex(x1,y1)
        glVertex(x2,y2)
        glVertex(x3,y3)
        glVertex(x4,y4)
    glEnd()
end

function rect(xcent,ycent,width,height)
    glBegin(GL_QUAD)
        glVertex(xcent+width,ycent-height)
        glVertex(xcent+width,ycent+height)
        glVertex(xcent-width,ycent+height)
        glVertex(xcent-widht,ycent-height)
    glEnd()
end

function triangle(x1,y1,x2,y2,x3,y3)
    glBegin(GL_POLYGON)
        glVertex(x1,y1)
        glVertex(x2,y2)
        glVertex(x3,y3)
    glEnd()
end

## Curves

#bezier
#bezierDetail
#bezierPoint
#bezierTangent
#curve
#curveDetail
#curvePoint
#curveTangent
#curveTightness

## 3D Primitives

#box

function sphere(xcent,ycent,radius)
    shapeList = glGenLists(1)
    glNewList(shapeList, GL_COMPILE)
        quad = gluNewQuadric()
        gluQuadricDrawStyle(quad, GLU_FILL)
        glTranslate(xcent,ycent)
        gluSphere(quad,radius,processingState.Sdet,processingState.Sdet)
        glTranslate(-xcent,-ycent)
        gluDeleteQuadric(quad)
    glEndList()
end

function sphereDetail(detail)
    processingState.Sdet = detail
end

## Attributes

#ellipseMode
#noSmooth
#rectMode
#smooth
#strokeCap
#strokeJoin
#strokeWeight

## Vertex

function beginShape()
    glBegin()
end

#bezierVertex
#curveVertex

function endShape()
    glEnd()
end

#quadraticVertex

function vertex(x,y)
    glVertex(x,y)
end

## Loading & Displaying

#shape
#shapeMode

# Color

## Setting

function background(r, g, b, a)
    glClearColor(r, g, b, a)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()
end

#colorMode

function fill(r, g, b, a)
    if processingState.noFill == false
        glColor(r, g, b, a)
    end
end

function noFill()
    processingState.noFill = (processingState.noFill ? false : true)
end

#noStroke
#stroke

## Creating & Reading

#alpha
#blue
#brightness
#color
#green
#hue
#lerpColor
#red
#saturation

# Image

#createImage

## Loading & Displaying

#image
#imageMode
#loadImage
#noTint
#requestImage
#tint

## Textures

#texture
#textureMode
#textureWrap

## Pixels

#blend
#copy
#filter
#get
#loadPixels
#set
#updatePixels

# Rendering

function blendMode()
    glEnable(GL_SRC_ALPHA, GL_ONE)
end

#createGraphics
#hint

## Shaders

#loadShader
#resetShader
#shader

# Typography

## Loading & Displaying

#createFont
#loadFont
#text
#textFont

## Attributes

#textAlign
#textLeading
#textMode
#textSize
#textWidth

## Metrics

#textAscent
#textDescent

end # module Processing


text:
using Base.Test
include("splinesSerial.jl")

@testset "splines" begin
	point = [0.5, 0.7, 1.2]
	V = [1.0,3.4,5.5,6.9,2.4,8.0];
	FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
	@testset "sel" begin
 		@test typeof(sel(1)(V))==Float64
		@test typeof(sel(1)(FV))==Array{Int64,1}
		@test sel(4)(V)==6.9
		@test sel(5)(FV)==[5,7,8]
	end
	@testset "larBernsteinBasis" begin
		@test typeof(larBernsteinBasis(sel(1),3)([0.5,2,7.3]))==Array{Float64,1}
		@test typeof(larBernsteinBasis(sel(5), 6)(V))==
		Array{Float64,1}
		H=[1.2,3.4,5.6,7.8,9.0]
		F=larBernsteinBasis(sel(4),5)(H)
		@test length(F)==6
	end

	@testset "larTensorProdSurface" begin
		controlpoints = [[[0,0,0],[2,-4,2]],[[0,3,1],[4,0,0]]]
		basis = larBernsteinBasis(sel(1), 1)
		point = [0.5, 2, 7.3]
		T=larTensorProdSurface([basis, basis], controlpoints)(point)
		@test typeof(T)==Array{Float64,1}
		@test length(T)==3
	end

	@testset "larBilinearSurface" begin
		controlpoints = [[[0,0,0],[2,-4,2]],[[0,3,1],[4,0,0]]]
		Bil=larBilinearSurface(controlpoints)([0.5,0.7])
		@test typeof(Bil)==Array{Float64,1}
		@test length(Bil)==3
	end

	@testset "larBiquadraticSurface" begin
		controlpoints = [[[0,0,0],[2,0,1],[3,1,1]],[[1,3,-1],[2,2,0],[3,2,0]],[[-2,4,0],[2,5,1],[1,3,2]]]
		mapping = larBiquadraticSurface(controlpoints)([0.7, 0.2])
		@test typeof(mapping)==Array{Float64,1}
		@test length(mapping)==3
	end

	@testset "larBicubicSurface" begin
		controlpoints = [[[0,0,0], [0,3,4], [0,6,3], [0,10,0]],[[3,0,2], [2,2.5,5], [3,6,5], [4,8,2]],[[6,0,2], [8,3,5], [7,6,4.5], [6,10,2.5]],[[10,0,0], [11,3,4], [11,6,3], [10,9,0]]]
		mapping = larBicubicSurface(controlpoints)([0.7,0.2])
		@test typeof(mapping)==Array{Float64,1}
		@test length(mapping)==3
	end

	@testset "larBezier"begin
		controldata = [[0,1],[0,0],[1,1],[1,0]]
		U = sel(2)
		@test typeof(larBezier(U, controldata)(point))==
		Array{Float64,1}
		G=larBezier(U, controldata)(point)
		@test length(G)==2
		@test typeof(larBezier(sel(3), controldata)(point))==
		Array{Float64,1}
		F=larBezier(sel(3), controldata)(point)
		@test length(F)==2
	end

	@testset "larBezierCurve" begin
		controldata = [[1,2,-5],[2,-4,2]]
		@test typeof(larBezierCurve(controldata)(point))==Array{Float64,1}
		@test length(larBezierCurve(controldata)(point))==3
	end

	@testset "bsplineBasis" begin
		knots = [0,0,0,1,1,2,2,3,3,4,4,4]
		ncontrols = 9
		degree = 2
		H=bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0])
		@test typeof(bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0]))==Array{Float64,1}
		@test length(bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0]))==9
	end

	@testset "larCoonsPatch" begin
		Su0 = larBezier(sel(1), [[0,0,0],[10,0,0]])
		Su1 = larBezier(sel(1), [[0,10,0],[2.5,10,3],[5,10,-3],[7.5,10,3],[10,10,0]])
		Sv0 = larBezier(sel(1), [[0,0,0],[0,0,3],[0,10,3],[0,10,0]])
		Sv1 = larBezier(sel(2), [[10,0,0],[10,5,3],[10,10,0]])
		point = [0.7, 0.2]
		out = larCoonsPatch([Su0,Su1,Sv0,Sv1])(point)
		@test typeof(out)==Array{Float64,1}
	end

	@testset "t_bspline"begin
		b1 = larBezier(S1, [[0,1,0],[0,1,5]])
		b2 = larBezier(S1, [[0,0,0],[0,0,5]])
		b3 = larBezier(S1, [[1,0,0],[2,-1,2.5],[1,0,5]])
		b4 = larBezier(S1, [[1,1,0],[1,1,5]])
		b5 = larBezier(S1, [[0,1,0],[0,1,5]])
		controls = [b1, b2, b3, b4, b5]
		knots = [0,1,2,3,4,5,6,7]           # periodic B-spline
		knots = [0,0,0,1,2,3,3,3]           # non-periodic B-spline
		tbspline = t_bspline(S2, 2, knots, controls)
		B=tbspline([0.7, 0.2])
		@test typeof(B)==Array{Float64,1}
		@test length(B)==3
	end

	@testset "t_nurbs" begin
		knots = [0,0,0,1,1,2,2,3,3,4,4,4]
		p=sqrt(2)/2.0
		controls = [[-1,0,1], [-p,p,p], [0,1,1], [p,p,p],[1,0,1], [p,-p,p],[0,-1,1], [-p,-p,p], [-1,0,1]]
		c1 = larBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
		c2 = larBezier(S1, [[-p,p,0,p],[-p,p,p,p]])
		c3 = larBezier(S1, [[0,1,0,1],[0,1,1,1]])
		c4 = larBezier(S1, [[p,p,0,p],[p,p,p,p]])
		c5 = larBezier(S1, [[1,0,0,1],[1,0,1,1]])
		c6 = larBezier(S1, [[p,-p,0,p],[p,-p,p,p]])
		c7 = larBezier(S1, [[0,-1,0,1],[0,-1,1,1]])
		c8 = larBezier(S1, [[-p,-p,0,p],[-p,-p,p,p]])
		c9 = larBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
		controls = [c1,c2,c3,c4,c5,c6,c7,c8,c9]
		tnurbs = t_nurbs(S2,2,knots,controls)
		N=tnurbs([0.7, 0.2])
		@test typeof(N)==Array{Float64,1}
		@test length(N)==3
	end

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export link

"""
    link(;kwargs...)
    link(children::Any;kwargs...)
    link(children_maker::Function;kwargs...)

A Link component.
Link component
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): Provide the content for the Link
- `id` (String; optional): The id of the link
- `style` (Dict; optional): Inline styles
- `className` (String; optional): Provide a custom className to be applied to the containing `<a>` node
- `disabled` (Bool; optional): Specify if the control should be disabled, or not
- `inline` (Bool; optional): Specify whether you want the inline version of this control
- `size` (a value equal to: 'sm', 'md', 'lg'; optional): Specify the size of the Link. Currently supports either `sm`, 'md' (default) or 'lg` as an option.
- `n_clicks` (Real; optional): Number of times this link has been clicked
- `visited` (Bool; optional): Specify whether you want the link to receive visited styles after the link has been clicked
- `href` (String; optional): Provide the `href` attribute for the `<a>` node
"""
function link(; kwargs...)
        available_props = Symbol[:children, :id, :style, :className, :disabled, :inline, :size, :n_clicks, :visited, :href]
        wild_props = Symbol[]
        return Component("link", "Link", "dash_carbon_components", available_props, wild_props; kwargs...)
end

link(children::Any; kwargs...) = link(;kwargs..., children = children)
link(children_maker::Function; kwargs...) = link(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


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

mutable struct RollbackStatusInfo <: SwaggerModel
    successfullyRolledbackInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: successfullyRolledbackInstanceCount
    failedRolledbackInstanceCount::Any # spec type: Union{ Nothing, Int32 } # spec name: failedRolledbackInstanceCount
    rollbackError::Any # spec type: Union{ Nothing, ApiError } # spec name: rollbackError

    function RollbackStatusInfo(;successfullyRolledbackInstanceCount=nothing, failedRolledbackInstanceCount=nothing, rollbackError=nothing)
        o = new()
        validate_property(RollbackStatusInfo, Symbol("successfullyRolledbackInstanceCount"), successfullyRolledbackInstanceCount)
        setfield!(o, Symbol("successfullyRolledbackInstanceCount"), successfullyRolledbackInstanceCount)
        validate_property(RollbackStatusInfo, Symbol("failedRolledbackInstanceCount"), failedRolledbackInstanceCount)
        setfield!(o, Symbol("failedRolledbackInstanceCount"), failedRolledbackInstanceCount)
        validate_property(RollbackStatusInfo, Symbol("rollbackError"), rollbackError)
        setfield!(o, Symbol("rollbackError"), rollbackError)
        o
    end
end # type RollbackStatusInfo

const _property_map_RollbackStatusInfo = Dict{Symbol,Symbol}(Symbol("successfullyRolledbackInstanceCount")=>Symbol("successfullyRolledbackInstanceCount"), Symbol("failedRolledbackInstanceCount")=>Symbol("failedRolledbackInstanceCount"), Symbol("rollbackError")=>Symbol("rollbackError"))
const _property_types_RollbackStatusInfo = Dict{Symbol,String}(Symbol("successfullyRolledbackInstanceCount")=>"Int32", Symbol("failedRolledbackInstanceCount")=>"Int32", Symbol("rollbackError")=>"ApiError")
Base.propertynames(::Type{ RollbackStatusInfo }) = collect(keys(_property_map_RollbackStatusInfo))
Swagger.property_type(::Type{ RollbackStatusInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RollbackStatusInfo[name]))}
Swagger.field_name(::Type{ RollbackStatusInfo }, property_name::Symbol) =  _property_map_RollbackStatusInfo[property_name]

function check_required(o::RollbackStatusInfo)
    true
end

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


text:
partial = false
include("lorenz.jl")
include("../extra/makie.jl")
using Makie, GLVisualize
using GeometryTypes, Colors
using Bridge: _viridis
skippoints = 1;
scene = Scene(resolution = (800, 800))
eyeposition = Float32[-37.2262, 65.8173, 37.4411]
sphere = Sphere(Point3f0(0,0,0), 1.0f0)
circle = Sphere(Point2f0(0,0), 1.0f0)
#RGB(Bridge._viridis[20]...)
maxviri = 200
viri = _viridis[round.(Int,range(1, stop=maxviri, length=length(Xtrue)))]
lines(( X2.yy), linewidth = 2.0, color = map(x->RGBA(Float32.(x)..., 0.9f0), viri) )
#scatter([X2.yy[1].+ SVector(-0.2,-0.2,1)], marker = '↘',  markersize = 0.8, color=:black)

ma = Float32[10,10,10]
axis(map(x -> UnitRange(0, round.(Int,x)), ma)..., showaxis = (true,true,true), showgrid = (false,false,false))
ma2 = ma .+ 2
for (x, c) in zip(ℝ{3}[(ma2[1], 0, 0), (0, ma2[2], 0),(0, 0, ma2[3]), (0,0,0)], ['x', 'y', 'z', '•'])
    scatter([x], marker = c, markersize = 1., color = :black)
end

scatter(collect(Models.foci(P)), marker=circle, markersize = .5,  strokewidth = 0.1, strokecolor = :white, color = :black)

scatter([Models.foci(P)[1] + SVector(-0.2,-0.2,1.0)], marker= 'a', markersize = .8, color = :black)
scatter([Models.foci(P)[2] + SVector(-0.2,-0.2,1.0)], marker='b', markersize = .8, color = :black)

center!(scene)
println(Makie.getscreen(scene).cameras[:perspective].view.value)

push!(Makie.getscreen(scene).cameras[:perspective].view,
    @SArray Float32[0.433338 0.901231 -3.7998f-7 1.95197; -0.237208 0.114057 0.96474 -20.4296; 0.869454 -0.418059 0.263205 -90.2712; 0.0 0.0 0.0 1.0]
)

Makie.save("output/lorenz1.png", scene);