text:
"""
    KernelTensorProduct

Tensor product of kernels.

# Definition

For inputs ``x = (x_1, \\ldots, x_n)`` and ``x' = (x'_1, \\ldots, x'_n)``, the tensor
product of kernels ``k_1, \\ldots, k_n`` is defined as
```math
k(x, x'; k_1, \\ldots, k_n) = \\Big(\\bigotimes_{i=1}^n k_i\\Big)(x, x') = \\prod_{i=1}^n k_i(x_i, x'_i).
```

# Construction

The simplest way to specify a `KernelTensorProduct` is to use the overloaded `tensor`
operator or its alias `⊗` (can be typed by `\\otimes<tab>`).
```jldoctest tensorproduct
julia> k1 = SqExponentialKernel(); k2 = LinearKernel(); X = rand(5, 2);

julia> kernelmatrix(k1 ⊗ k2, RowVecs(X)) == kernelmatrix(k1, X[:, 1]) .* kernelmatrix(k2, X[:, 2])
true
```

You can also specify a `KernelTensorProduct` by providing kernels as individual arguments
or as an iterable data structure such as a `Tuple` or a `Vector`. Using a tuple or
individual arguments guarantees that `KernelTensorProduct` is concretely typed but might
lead to large compilation times if the number of kernels is large.
```jldoctest tensorproduct
julia> KernelTensorProduct(k1, k2) == k1 ⊗ k2
true

julia> KernelTensorProduct((k1, k2)) == k1 ⊗ k2
true

julia> KernelTensorProduct([k1, k2]) == k1 ⊗ k2
true
```
"""
struct KernelTensorProduct{K} <: Kernel
    kernels::K
end

function KernelTensorProduct(kernel::Kernel, kernels::Kernel...)
    return KernelTensorProduct((kernel, kernels...))
end

@functor KernelTensorProduct

Base.length(kernel::KernelTensorProduct) = length(kernel.kernels)

function (kernel::KernelTensorProduct)(x, y)
    if !(length(x) == length(y) == length(kernel))
        throw(DimensionMismatch("number of kernels and number of features
are not consistent"))
    end
    return prod(k(xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
end

function validate_domain(k::KernelTensorProduct, x::AbstractVector)
    return dim(x) == length(k) ||
           error("number of kernels and groups of features are not consistent")
end

# Utility for slicing up inputs.
slices(x::AbstractVector{<:Real}) = (x,)
slices(x::ColVecs) = eachrow(x.X)
slices(x::RowVecs) = eachcol(x.X)

function kernelmatrix!(K::AbstractMatrix, k::KernelTensorProduct, x::AbstractVector)
    validate_inplace_dims(K, x)
    validate_domain(k, x)

    first_kernels, tail_kernels = Iterators.peel(k.kernels)
    first_x, tail_x = Iterators.peel(slices(x))

    # handle first kernel and input
    kernelmatrix!(K, first_kernels, first_x)

    # handle remaining kernels and inputs
    Ktmp = similar(K)
    for (ki, xi) in zip(tail_kernels, tail_x)
        kernelmatrix!(Ktmp, ki, xi)
        hadamard!(K, K, Ktmp)
    end

    return K
end

function kernelmatrix!(
    K::AbstractMatrix, k::KernelTensorProduct, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    validate_domain(k, x)

    first_kernels, tail_kernels = Iterators.peel(k.kernels)
    first_x, tail_x = Iterators.peel(slices(x))
    first_y, tail_y = Iterators.peel(slices(y))

    # handle first kernel and inputs
    kernelmatrix!(K, first_kernels, first_x, first_y)

    # handle remaining kernels and inputs
    Ktmp = similar(K)
    for (ki, xi, yi) in zip(tail_kernels, tail_x, tail_y)
        kernelmatrix!(Ktmp, ki, xi, yi)
        hadamard!(K, K, Ktmp)
    end

    return K
end

function kernelmatrix_diag!(K::AbstractVector, k::KernelTensorProduct, x::AbstractVector)
    validate_inplace_dims(K, x)
    validate_domain(k, x)

    first_kernels, tail_kernels = Iterators.peel(k.kernels)
    first_x, tail_x = Iterators.peel(slices(x))

    # handle first kernel and input
    kernelmatrix_diag!(K, first_kernels, first_x)

    # handle remaining kernels and inputs
    Ktmp = similar(K)
    for (ki, xi) in zip(tail_kernels, tail_x)
        kernelmatrix_diag!(Ktmp, ki, xi)
        hadamard!(K, K, Ktmp)
    end

    return K
end

function kernelmatrix_diag!(
    K::AbstractVector, k::KernelTensorProduct, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    validate_domain(k, x)

    first_kernels, tail_kernels = Iterators.peel(k.kernels)
    first_x, tail_x = Iterators.peel(slices(x))
    first_y, tail_y = Iterators.peel(slices(y))

    # handle first kernel and inputs
    kernelmatrix_diag!(K, first_kernels, first_x, first_y)

    # handle remaining kernels and inputs
    Ktmp = similar(K)
    for (ki, xi, yi) in zip(tail_kernels, tail_x, tail_y)
        kernelmatrix_diag!(Ktmp, ki, xi, yi)
        hadamard!(K, K, Ktmp)
    end

    return K
end

function kernelmatrix(k::KernelTensorProduct, x::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix, hadamard, k.kernels, slices(x))
end

function kernelmatrix(k::KernelTensorProduct, x::AbstractVector, y::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix, hadamard, k.kernels, slices(x), slices(y))
end

function kernelmatrix_diag(k::KernelTensorProduct, x::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix_diag, hadamard, k.kernels, slices(x))
end

function kernelmatrix_diag(k::KernelTensorProduct, x::AbstractVector, y::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix_diag, hadamard, k.kernels, slices(x), slices(y))
end

Base.show(io::IO, kernel::KernelTensorProduct) = printshifted(io, kernel, 0)

function Base.:(==)(x::KernelTensorProduct, y::KernelTensorProduct)
    return (
        length(x.kernels) == length(y.kernels) &&
        all(kx == ky for (kx, ky) in zip(x.kernels, y.kernels))
    )
end

function printshifted(io::IO, kernel::KernelTensorProduct, shift::Int)
    print(io, "Tensor product of ", length(kernel), " kernels:")
    for k in kernel.kernels
        print(io, "\n")
        for _ in 1:(shift + 1)
            print(io, "\t")
        end
        printshifted(io, k, shift + 2)
    end
end


text:
"""
$(SIGNATURES)

This function turns an initial guess for a Fold/Hopf point into a solution to the Fold/Hopf problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `J  = (x, p) -> d_xF(x, p)` associated jacobian
- `br` results returned after a call to [`continuation`](@ref)
- `ind_bif` bifurcation index in `br`
- `par` parameters used for the vector field
- `lens` parameter axis used to locate the Fold/Hopf point.
- `options::NewtonPar`

# Optional arguments:
- `Jt = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jt` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jt = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1,v2]`.
- `normN = norm`
- `options` You can pass newton parameters different from the ones stored in `br` by using this argument `options`.
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver
"""
function newton(F, J, br::ContResult, ind_bif::Int64, par, lens::Lens; Jt = nothing, d2F = nothing, normN = norm, options = br.contparams.newtonOptions, kwargs...)
	if length(br.bifpoint) > 0 && br.bifpoint[ind_bif].type == :hopf
		return newtonHopf(F, J, br, ind_bif, par, lens; Jt = Jt, d2F = d2F, normN = normN, options = options, kwargs...)
	elseif br.foldpoint[ind_bif].type == :fold
		return newtonFold(F, J, br, ind_bif, par, lens; Jt = Jt, d2F = d2F, normN = normN, options = options, kwargs...)
	end
	@error "Bifurcation type $(br[ind_bif].type) not yet handle for codim2 newton / continuation"
end

"""
$(SIGNATURES)

codim 2 continuation of Fold / Hopf points. This function turns an initial guess for a Fold/Hopf point into a curve of Fold/Hopf points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `br` results returned after a call to [`continuation`](@ref)
- `ind_bif` bifurcation index in `br`
- `par` set of parameters
- `lens1` parameter axis for parameter 1
- `lens2` parameter axis for parameter 2
- `options_cont` arguments to be passed to the regular [`continuation`](@ref)

# Optional arguments:

- `Jt = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose
- `d2F = p -> ((x, p, v1, v2) -> d2F(x, p, v1, v2))` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`.

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Hopf point in `br` you want to refine.

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""

function continuation(F, J, br::ContResult, ind_bif::Int64, par, lens1::Lens, lens2::Lens, options_cont::ContinuationPar ; Jt = nothing, d2F = nothing, kwargs...)
	if length(br.bifpoint) > 0 && br.bifpoint[ind_bif].type == :hopf
		return continuationHopf(F, J, br, ind_bif, par, lens1, lens2, options_cont; Jt = Jt, d2F = d2F, kwargs...)
	elseif br.foldpoint[ind_bif].type == :fold
		return continuationFold(F, J, br, ind_bif, par, lens1, lens2, options_cont; Jt = Jt, d2F = d2F, kwargs...)
	end
	@error "Bifurcation type $(br[ind_bif].type) not yet handle for codim2 newton / continuation"
end

ext: jl


text:
#=
hello:
- Julia version: 1.0.0
- Author: Baur Urazalinov
- Date: 2018-09-10
=#

println("Hello World!")


text:
module TestCheckboxes

using FLoops:
    _box_detection_works,
    _make_closure_with_a_box,
    _make_closure_without_a_box,
    has_boxed_variables,
    verify_no_boxes,
    HasBoxedVariableError
using Test

function _make_closure_with_two_boxes()
    local a
    local b = 1
    set(v) = (z = a; a = b; b = v; z)
    return set
end

function test__box_detection_works()
    @test _box_detection_works()
end

function test_verify_no_boxes()
    @test (verify_no_boxes(_make_closure_without_a_box()); true)
    f = _make_closure_with_a_box()
    @test_throws HasBoxedVariableError(f) verify_no_boxes(f)
end

function test_HasBoxedVariableError()
    @testset "one box" begin
        err = HasBoxedVariableError(_make_closure_with_a_box())
        @debug "Example `HasBoxedVariableError`" err
        msg = sprint(showerror, err)
        @test occursin("Consider adding declarations such as `local a`", msg)
        @test occursin("1 boxed variable:", msg)
        @test occursin(r"\b a \b"x, msg)
    end
    @testset "two boxes" begin
        err = HasBoxedVariableError(_make_closure_with_two_boxes())
        @debug "Example `HasBoxedVariableError`" err
        msg = sprint(showerror, err)
        @test occursin(r"Consider adding declarations such as `local (a|b)`", msg)
        @test occursin("2 boxed variables:", msg)
        @test occursin(r"\b a \b"x, msg)
        @test occursin(r"\b b \b"x, msg)
    end
end

function with_box_1()
    a = 0
    function closure()
        a += 1
    end
    Val(has_boxed_variables(closure))
end

function with_box_3()
    (a, b, c) = (1, 2, 3)
    function closure()
        c = a + b + c
    end
    Val(has_boxed_variables(closure))
end

function with_box_7()
    (a, b, c, d, e, f, g) = (1, 2, 3, 4, 5, 6, 7)
    function closure()
        g = a + b + c + d + e + f + g
    end
    Val(has_boxed_variables(closure))
end

function with_box_10()
    (a, b, c, d, e, f, g, h, i, j) = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
    function closure()
        j = a + b + c + d + e + f + g + h + i + j
    end
    Val(has_boxed_variables(closure))
end

function test_inferrability()
    @test @inferred(with_box_1()) == Val(true)
    @test @inferred(with_box_3()) == Val(true)
    @test @inferred(with_box_7()) == Val(true)
    @test @inferred(with_box_10()) == Val(true)
end

end  # module


text:
module Spinner

using ..Ahorn, Maple

# Spinners in the base game without texture, faster than looking up sprite and trying to load it
const spinnerWithoutTexture = Set{String}([
    "Core",
    "Rainbow",
    "core",
    "rainbow"
])

const defaultSpinnerTexture = "danger/crystal/fg_blue03"

function rotatingSpinnerFinalizer(entity::Maple.RotateSpinner)
    x, y = Int(entity.data["x"]), Int(entity.data["y"])
    entity.data["x"], entity.data["y"] = x + 32, y
    entity.data["nodes"] = [(x, y)]
end

const movingSpinner = Union{Maple.TrackSpinner, Maple.RotateSpinner}

const placements = Ahorn.PlacementDict(
    "Dust Sprite" => Ahorn.EntityPlacement(
        Maple.Spinner,
        "point",
        Dict{String, Any}(
            "dust" => true
        )
    ),
    "Dust Sprite (Attached)" => Ahorn.EntityPlacement(
        Maple.Spinner,
        "point",
        Dict{String, Any}(
            "dust" => true,
            "attachToSolid" => true,
        )
    ),

    "Dust Sprite (Rotating, Clockwise)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => true,
            "clockwise" => true
        ),
        rotatingSpinnerFinalizer
    ),
    "Dust Sprite (Rotating)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => true,
            "clockwise" => false
        ),
        rotatingSpinnerFinalizer
    ),
    "Blade (Rotating, Clockwise)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => false,
            "clockwise" => true
        ),
        rotatingSpinnerFinalizer
    ),
    "Blade (Rotating)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => false,
            "clockwise" => false
        ),
        rotatingSpinnerFinalizer
    ),
    "Star (Rotating, Clockwise)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => false,
            "star" => true,
            "clockwise" => true
        ),
        rotatingSpinnerFinalizer
    ),
    "Star (Rotating)" => Ahorn.EntityPlacement(
        Maple.RotateSpinner,
        "point",
        Dict{String, Any}(
            "dust" => false,
            "star" => true,
            "clockwise" => false
        ),
        rotatingSpinnerFinalizer
    ),
)

crystalSpinnerColors = Maple.crystal_colors
for color in crystalSpinnerColors
    for attached in false:true
        key = "Crystal Spinner ($(uppercasefirst(color))$(attached ? ", Attached" : ""))"
        placements[key] = Ahorn.EntityPlacement(
            Maple.Spinner,
            "point",
            Dict{String, Any}(
                "color" => color,
                "attachToSolid" => attached,
            )
        )
    end
end

speeds = Maple.track_spinner_speeds
for speed in speeds, dusty in false:true, star in false:true
    if star && dusty
        continue
    end

    key = (star ? "Star" : (dusty ? "Dust Sprite" : "Blade")) * " (Track, $(uppercasefirst(speed)))"
    placements[key] = Ahorn.EntityPlacement(
        Maple.TrackSpinner,
        "line",
        Dict{String, Any}(
            "dust" => dusty,
            "star" => star,
            "speed" => speed
        )
    )
end

Ahorn.editingOptions(entity::Maple.Spinner) = Dict{String, Any}(
    "color" => crystalSpinnerColors
)

Ahorn.editingOptions(entity::Maple.TrackSpinner) = Dict{String, Any}(
    "color" => crystalSpinnerColors,
    "speed" => speeds
)

Ahorn.nodeLimits(entity::Maple.TrackSpinner) = 1, 1
Ahorn.nodeLimits(entity::Maple.RotateSpinner) = 1, 1

function Ahorn.selection(entity::Maple.Spinner)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 8, y - 8, 16, 16)
end

function Ahorn.selection(entity::Maple.RotateSpinner)
    nx, ny = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    return [Ahorn.Rectangle(x - 8, y - 8, 16, 16), Ahorn.Rectangle(nx - 8, ny - 8, 16, 16)]
end

function Ahorn.selection(entity::Maple.TrackSpinner)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    return [Ahorn.Rectangle(startX - 8, startY - 8, 16, 16), Ahorn.Rectangle(stopX - 8, stopY - 8, 16, 16)]
end

function renderMovingSpinner(ctx::Ahorn.Cairo.CairoContext, entity::movingSpinner, x::Number, y::Number)
    dusty = get(entity.data, "dust", false)
    star = get(entity.data, "star", false)

    if star
        Ahorn.drawSprite(ctx, "danger/starfish13", x, y)

    elseif dusty
        Ahorn.drawSprite(ctx, "danger/dustcreature/base00", x, y)
        Ahorn.drawSprite(ctx, "danger/dustcreature/center00", x, y)

    else
        Ahorn.drawSprite(ctx, "danger/blade00", x, y)
    end
end

# TODO - Add support for background
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Spinner, room::Maple.Room)
    # Custom attributes from Everest
    dusty = get(entity.data, "dust", false)
    color = lowercase(get(entity.data, "color", "blue"))

    if dusty
        Ahorn.drawSprite(ctx, "danger/dustcreature/base00", 0, 0)
        Ahorn.drawSprite(ctx, "danger/dustcreature/center00", 0, 0)

    else
        resource = "danger/crystal/fg_$(color)03"

        if color in spinnerWithoutTexture
            Ahorn.drawSprite(ctx, defaultSpinnerTexture, 0, 0)

        else
            fgSprite = Ahorn.getSprite(resource, "Gameplay")
            texture = fgSprite.width == 0 || fgSprite.height == 0 ? defaultSpinnerTexture : resource

            Ahorn.drawSprite(ctx, texture, 0, 0)
        end
    end
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.RotateSpinner, room::Maple.Room)
    clockwise = get(entity.data, "clockwise", false)
    dir = clockwise ? 1 : -1

    centerX, centerY = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    radius = sqrt((centerX - x)^2 + (centerY - y)^2)

    Ahorn.drawCircle(ctx, centerX, centerY, radius, Ahorn.colors.selection_selected_fc)
    Ahorn.drawArrow(ctx, centerX + radius, centerY, centerX + radius, centerY + 0.001 * dir, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawArrow(ctx, centerX - radius, centerY, centerX - radius, centerY + 0.001 * -dir, Ahorn.colors.selection_selected_fc, headLength=6)

    renderMovingSpinner(ctx, entity, x, y)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = entity.data["nodes"][1]

    renderMovingSpinner(ctx, entity, stopX, stopY)
    Ahorn.drawArrow(ctx, startX, startY, stopX, stopY, Ahorn.colors.selection_selected_fc, headLength=10)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.RotateSpinner, room::Maple.Room)
    centerX, centerY = Int.(entity.data["nodes"][1])

    renderMovingSpinner(ctx, entity, centerX, centerY)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)

    renderMovingSpinner(ctx, entity, startX, startY)
end

end


text:
@defstruct CropLayer Layer (
  name :: AbstractString = "crop",
  random_crop :: Bool = false,
  random_mirror :: Bool = false,
  (crop_size :: NTuple{2, Int} = (0,0), crop_size[1] > 0 && crop_size[2] > 0),
  (bottoms :: Vector{Symbol} = [], length(bottoms) > 0),
  (tops :: Vector{Symbol} = [], length(tops) == length(bottoms))
)

type CropLayerState <: LayerState
  layer      :: CropLayer
  blobs      :: Vector{Blob}
end

function setup(backend::Backend, layer::CropLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  for i = 1:length(inputs)
    @assert ndims(inputs[i]) == 4 "Input blob $(layer.bottoms[i]) should be a 4D tensor"
    @assert isa(diffs[i], NullBlob) # Back-propagation for crop-layer is not implemented
  end

  blobs = Array{Blob}(length(inputs))
  for i = 1:length(inputs)
    width, height, channels, num = size(inputs[i])
    @assert layer.crop_size[1] <= width && layer.crop_size[2] <= height
    blobs[i] = make_blob(backend, eltype(inputs[i]),
        (layer.crop_size[1], layer.crop_size[2],channels,num))
  end

  return CropLayerState(layer, blobs)
end

function shutdown(backend::Backend, state::CropLayerState)
  map(destroy, state.blobs)
end

function crop_blob{T}(input::Array{T}, output::Array{T}, crop_size::NTuple{2,Int}, offsets::NTuple{2,Int})
  crop_w = crop_size[1]; w_off = offsets[1]
  crop_h = crop_size[2]; h_off = offsets[2]
  num = size(input, 4); channels = size(input, 3)

  for n = 1:num
    for c = 1:channels
      for h = 1:crop_h
        @simd for w = 1:crop_w
          @inbounds output[w,h,c,n] = input[w+w_off,h+h_off,c,n]
        end
      end
    end
  end
end
function mirror_crop_blob{T}(input::Array{T}, output::Array{T}, crop_size::NTuple{2,Int}, offsets::NTuple{2,Int})
  crop_w = crop_size[1]; w_off = offsets[1]
  crop_h = crop_size[2]; h_off = offsets[2]
  num = size(input, 4); channels = size(input, 3)

  for n = 1:num
    for c = 1:channels
      for h = 1:crop_h
        @simd for w = 1:crop_w
          @inbounds output[crop_w-w+1,h,c,n] = input[w+w_off,h+h_off,c,n]
        end
      end
    end
  end
end

function forward(backend::CPUBackend, state::CropLayerState, inputs::Vector{Blob})
  crop_size = state.layer.crop_size
  for i = 1:length(inputs)
    input = inputs[i].data
    output = state.blobs[i].data
    if state.layer.random_crop
      w_off = abs(rand(Int)) % (size(input,1)-crop_size[1]+1)
      h_off = abs(rand(Int)) % (size(input,2)-crop_size[2]+1)
    else
      w_off = div(size(input,1)-crop_size[1], 2)
      h_off = div(size(input,2)-crop_size[2], 2)
    end
    if state.layer.random_mirror && rand(UInt)%2 == 0
      mirror_crop_blob(input, output, crop_size, (w_off, h_off))
    else
      crop_blob(input, output, crop_size, (w_off, h_off))
    end
  end
end

function backward(backend::Backend, state::CropLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  # backward for a crop layer could be implemented, but since crop layer is typically used
  # directly on top of a data layer, which does not need back propagation, we don't implement
  # backward for CropLayer
end


text:
## Casette Powered Tracking
## ========================
Cassette.@context TrackErrorCtx

function Cassette.execute(ctx::TrackErrorCtx, ::typeof(condf), ω, x, y)
  conjoinerror!(ctx.metadata, Cassette.overdub(ctx, y, ω))
  Cassette.overdub(ctx, x, ω)
end

function casetteapplytrackerr(f, args...)
  sbw = SoftBoolWrapper(SoftBool(Val{true}))
  fx = Cassette.overdub(TrackErrorCtx(metadata=sbw), f, args...)
  (fx, sbw.sb)
end


text:
#ДАНО: Робот - рядом с горизонтальной перегородкой (под ней), бесконечно продолжающейся в обе стороны,
#в которой имеется проход шириной в одну клетку.
#РЕЗУЛЬТАТ: Робот - в клетке под проходом

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function check_border(r, side, num) #двигает робота на num шагов и проверяет наличие прохода сверху
    for i in 1:num
        move!(r,side)
        if isborder(r,Nord) == false
            break
        end
    end
end

function nord_space(r)
    num = 1
    side = Ost
    while isborder(r,Nord) == true
        check_border(r,side,num)
        side = reverse(side)
        num += 1
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
fname = :ST_TRsparse
c = Combi(hessian_sparse,PDataST,solve_modelST_TR,preprocessST,decreaseGen,TparamsST())
include("Template.jl")


text:
transD_GP.SkyTEM1DInversion.splitlineconvandlast(soundings, 50, 2; zstart=zstart, extendfrac=extendfrac, dz=dz, nlayers=nlayers,
    figsize=(14,7), vmin=lo, vmax=hi, postfix="_$(regtype)_β²_$(β²)_bg_$(10. ^σ0)Spm", preferEright=true,
    showplot=true, logscale=true, prefix=zipsaveprefix)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
module Tests

include("matching.jl")
include("destructuring.jl")
include("dispatch.jl")
include("metamodule.jl")
include("reflection.jl")
include("slurp-optimizations.jl")
include("extra.jl")

end


text:
##

using Flux
using CUDA
using Transformers
using Transformers.Basic
using Flux: onehot, onecold
using Transformers: PositionEmbedding
#flux function for update parameters
using Flux: gradient,  @functor
using Flux.Optimise: update!
using Flux:onehotbatch, logitcrossentropy
using Statistics: norm, mean, std
using Flux: Momentum
## ---  Model and Datasoruce

struct GPTPositionEmbedding{F, W <: AbstractArray{F}}
	embedding::W
end

@functor GPTPositionEmbedding

function GPTPositionEmbedding(size::Int, max_len::Int = 1024)
	GPTPositionEmbedding( zeros(Float32,size, max_len))
end

function (pe::GPTPositionEmbedding{F})(x::AbstractArray{F}) where F
	pe.embedding[:, 1:size(x,2)]  # just return the needed position embedding
end

######
labels = collect(1:10)
startsym = 11
endsym = 12
unksym = 0
labels = [unksym, startsym, endsym, labels...]
vocab = Vocabulary(labels, unksym)

#function for generate training datas
sample_data() = rand(1:10, 10)
#function for adding start & end symbol
preprocess(x) = [startsym, x..., endsym]
#
# sample = preprocess.(sample_data())
# encoded_sample = vocab(sample[1]) #use Vocabulary to encode the training data
#
struct  SimpleTransfomerEncoder
	embd::Embed
	pe::GPTPositionEmbedding
	#pe::PositionEmbedding
	encode1::Transformer
	encode2::Transformer
end

@functor SimpleTransfomerEncoder

struct  SimpleTransfomerDecoder
	embd::Embed
	pe::GPTPositionEmbedding
	#pe::PositionEmbedding
	decode1::TransformerDecoder
	decode2::TransformerDecoder
	linear::Positionwise
end

@functor SimpleTransfomerDecoder

function SimpleTransfomer()
	#512 is the size of internal embedding representation
	embd=Embed(512, length(vocab))
	pe = GPTPositionEmbedding(512)
	#pe = PositionEmbedding(512, maxInputLen; trainable=true)   # max length of input string

	embd_dec=Embed(512, length(vocab))
	pe_dec = GPTPositionEmbedding(512)
	#pe_dec = PositionEmbedding(512, maxInputLen; trainable=true)

	encode_t1 = Transformer(512, 8, 64, 2048)
	encode_t2 = Transformer(512, 8, 64, 2048)
	decode_t1 = TransformerDecoder(512, 8, 64, 2048)
	decode_t2 = TransformerDecoder(512, 8, 64, 2048)
	linear = Positionwise(Dense(512, length(vocab)))
    encoder=SimpleTransfomerEncoder(embd, pe, encode_t1, encode_t2)
	decoder=SimpleTransfomerDecoder(embd_dec, pe_dec, decode_t1,decode_t2, linear)
	encoder, decoder
end

function (t::SimpleTransfomerEncoder)(x)
	we = t.embd(x, inv(sqrt(512)))
	e = we .+ t.pe(we)
	e=t.encode1(e)
	t.encode2(e)
end
function (t::SimpleTransfomerDecoder)(x, m)
	we = t.embd(x, inv(sqrt(512)))
	e = we .+ t.pe(we)
	e = t.decode1(e, m)
	e = t.decode2(e, m)
	t.linear(e)
end

## ---

opt = ADAM(1e-4)
#opt=Momentum(0.001)
model = SimpleTransfomer()

## ---

using Flux: onecold
# assume x is list of numbers, This method will add start and end symbol,
# encode it as vocab, then feed it to the network and decode the output
# resulting output should have begin and end symbol
#list of numbers to a list of numbers with start and finish.
function translate(x, encoder, decoder)
    seq = [startsym]
	ix=vocab.(preprocess(x))  #add start and end, and vocab encoded
	encoderBatchOutput=encoder(ix)

	len = length(ix)
	local modelOut
    for i = 1:2len     # go twice as long, stop when there is end symbol
		modelOut = decoder(vocab(seq), encoderBatchOutput)
		dec = softmax(modelOut)  # distribution over vocab
		ntok = onecold(vocab, dec)  # decode from vocab to labels
        push!(seq, ntok[end])
        if ntok[end] == endsym
		 	break
		 end
    end
	seq   # note seq contains decoded output from itereation of the loop, one symbol at the time
  end

## ---

function train!(enc, dec, cnt=2000)
	@info "start training"
	local l
	local y, yhat
	ps=params(enc, dec)
	for i = 0:cnt
		trainingSize = 1:100
		batch=hcat((vocab.(preprocess(sample_data())) for i = trainingSize)...)
		#view the data as one long stream, then onehot encode it

		#assume output should be like input
		expected=batch[2:end,:]
		expectedStream=onehot(vocab, view(expected, :))  # using vocab one hot the array

		gs = gradient(ps) do
	        input = batch[1:end-1, :]
			encoderBatchOutput= enc(input)
			yhat = dec(input, encoderBatchOutput)   # output should be vocabs
			generatedStream = reshape(yhat, size(yhat)[1],:)

			l=logitcrossentropy(generatedStream, expectedStream)
		end
		if i % 10 == 0
			validateModel( enc, dec)
			@show i, l
		end
		update!(opt, ps, gs)
	end
end

## ---

function lossMeasure(b, yhat)
	expected = view(b, :)  # one dim list of output, vocab encoded
    y=onehot(vocab, expected)
	generatedStream = reshape(yhat, size(yhat)[1],:)
	logitcrossentropy(generatedStream, y)
end

# assuming batch is vocab encoded,
function doAModelRun(batch, encoder, decoder)
	encoderBatchOutput=encoder(batch)
	yhat = decoder(batch, encoderBatchOutput)
	@show loss=lossMeasure(batch, yhat)
end

function validateModel(encoder, decoder)
    validationData = [5,5,6,4,1,2,3,4,7,10]
	validationResult = translate(validationData, encoder, decoder)
	println(validationData, "->", validationResult)
end

## ---

train!(model...)
## ---  In case we need to do a step by step training

before = validateModel(model...)
@show before
#Get a batch, try run the model for a batch before and after a single training loop
batch=hcat(vocab.(preprocess(sample_data()) for i = 1:2)...)
@show doAModelRun(batch, model...)
train!(model..., 1)
@show doAModelRun(batch, model...)

#Translate the fix data again
after = validateModel(model ...)
@show after

max_stars_repo_name: daryoush/CudaAndFluxExperiments

max_issues_repo_name: daryoush/CudaAndFluxExperiments

max_forks_repo_name: daryoush/CudaAndFluxExperiments


text:
using JuMP, Base.Test
include("../src/OnePhase.jl")
include("problems.jl")
include("kkt_system_solvers.jl")
include("linear_system_solvers.jl")

function unit_tests()
    test_compare_columns()
    test_compute_indicies()
    test_linear_solvers()
    test_kkt_solvers()
end

function basic_tests(solver)
    @testset "rosenbrook" begin
        test_rosenbrook1(solver)
        test_rosenbrook2(solver)
        test_rosenbrook3(solver)
        test_rosenbrook4(solver)
    end

    @testset "LP" begin
        test_toy_lp1(solver)
        test_toy_lp2(solver)
        test_toy_lp3(solver)
        test_toy_lp4(solver)
        test_toy_lp5(solver)
        test_toy_lp6(solver)
        test_toy_lp7(solver)
    end

    @testset "infeasible" begin
        model = toy_lp_inf1()
        setsolver(model,solver)
        @test :Infeasible == solve(model)

        model = toy_lp_inf2()
        setsolver(model,solver)
        @test :Infeasible == solve(model)

        model = circle_nc_inf1()
        setsolver(model,solver)
        @test :Infeasible == solve(model)
    end

    @testset "convex_nlp" begin
        @testset "circle1" begin
            model = circle1()
            setsolver(model,solver)
            @test :Optimal == solve(model)
            check_circle1(model)
        end
        @testset "circle2" begin
            model = circle2()
            setsolver(model,solver)
            @test :Optimal == solve(model)
            check_circle2(model)
        end

        @testset "quad_opt" begin
            model = quad_opt()
            setsolver(model,solver)
            if @test :Optimal == solve(model)
                check_quad_opt(model)
            end
        end
    end

    @testset "nonconvex" begin
        model = circle_nc1()
        setsolver(model,solver)
        @test :Optimal == solve(model)
        check_circle_nc1(model)

        model = circle_nc2()
        setsolver(model,solver)
        @test :Optimal == solve(model)
        check_circle_nc2(model)
    end

    @testset "unbounded_opt_val" begin
        model = lp_unbd()
        setsolver(model,solver)
        @test :Unbounded == solve(model)

        model = circle_nc_unbd()
        setsolver(model,solver)
        status = solve(model)
        @test status == :Unbounded

        model = quad_unbd()
        setsolver(model,solver)
        status = solve(model)
        @test status == :Unbounded
    end

    @testset "unbounded_feasible_region" begin
        # problems with unbounded feasible region
        test_unbd_feas(solver)
    end

    @testset "starting point" begin
        test_starting_point(solver,0.5)
        test_starting_point(solver,-0.5)
    end
end

function basic_tests()
    max_it = 100
    output_level = 3
    a_norm_penalty = 1e-4
    @testset "basic_tests" begin
        @testset "cholseky linear system solve" begin
            solver = OnePhase.OnePhaseSolver(term!max_it=max_it,
            a_norm_penalty = a_norm_penalty,
            output_level=output_level,
            kkt!kkt_solver_type=:schur)
            basic_tests(solver)
        end
        @testset "Ma57 linear system solve" begin
            solver = OnePhase.OnePhaseSolver(term!max_it=max_it,
            a_norm_penalty = a_norm_penalty,
            output_level=output_level,
            kkt!kkt_solver_type=:symmetric,
            kkt!linear_solver_type=:HSL)

            basic_tests(solver)
        end
        @testset "Ma57 linear system solve with clever elimination" begin
            solver = OnePhase.OnePhaseSolver(term!max_it=max_it,
            a_norm_penalty = a_norm_penalty,
            output_level=output_level,
            kkt!kkt_solver_type=:clever_symmetric,
            kkt!linear_solver_type=:HSL)
            basic_tests(solver)
        end

        @testset "LDLT julia linear system solve" begin
            solver = OnePhase.OnePhaseSolver(term!max_it=max_it,
            a_norm_penalty = a_norm_penalty,
            output_level=output_level,
            kkt!kkt_solver_type=:clever_symmetric,
            kkt!linear_solver_type=:julia)
            basic_tests(solver)
        end
    end
end

# lets run the tests!
unit_tests()
basic_tests()


text:
"""
This file initializes the IscaTools module by setting and determining the
ECMWF reanalysis parameters to be analyzed and the regions upon which the data
are to be extracted from.  Functionalities include:
    - Setting up of reanalysis module type
    - Setting up of reanalysis parameters to be analyzed
    - Setting up of time steps upon which data are to be downloaded
    - Setting up of region of analysis based on ClimateEasy

"""

function iscananmean(data)
    dataii = @view data[data .!= NaN32]
    if dataii != []; return mean(dataii); else; return NaN32; end
end

function iscaanalysis(
    init::AbstractDict, iroot::AbstractDict;
    modID::AbstractString, parID::AbstractString,
    plvls::Union{AbstractString,Real,Vector{<:Real}}
)

    if typeof(plvls) <: Array
        for p in plvls
            imod,ipar,itime = iscainitialize(init;modID=modID,parID=parID,pressure=p);
            iscaanalysis(imod,ipar,itime,iroot); putinfo(imod,ipar,itime,iroot)
        end
    else
        if plvls == "sfc"
              imod,ipar,itime = iscainitialize(init;modID=modID,parID=parID)
        else; imod,ipar,itime = iscainitialize(init;modID=modID,parID=parID,pressure=plvls);
        end
        iscaanalysis(imod,ipar,itime,iroot); putinfo(imod,ipar,itime,iroot)
    end

end

function iscaanalysis(
    imod::AbstractDict, ipar::AbstractDict, itime::AbstractDict, iroot::AbstractDict
)

    if itime["nhr"] > 1;

        if any(uppercase(itime["calendar"]) .== ["NOLEAP","GREGORIAN"])
            for irun = 1 : itime["nruns"]
                iscaanalysis365diurnal(imod,ipar,itime,iroot,irun);
            end
        elseif itime["ndy"] == 360
            for irun = 1 : itime["nruns"];
                iscaanalysis360diurnal(imod,ipar,itime,iroot,irun);
            end
        elseif itime["ndy"] == 30 && mod(12,itime["nruns"]) == 0
            nruns = round(Int,itime["nruns"]/12)
            for irun = 1 : nruns;
                iscaanalysis30diurnal(imod,ipar,itime,iroot,irun);
            end
        else
            for irun = 1 : itime["nruns"];
                iscaanalysisgendiurnal(imod,ipar,itime,iroot,irun);
            end
        end

    else

        if any(uppercase(itime["calendar"]) .== ["NOLEAP","GREGORIAN"])
            for irun = 1 : itime["nruns"]
                iscaanalysis365daily(imod,ipar,itime,iroot,irun);
            end
        elseif itime["ndy"] == 360
            for irun = 1 : itime["nruns"];
                iscaanalysis360daily(imod,ipar,itime,iroot,irun);
            end
        elseif itime["ndy"] == 30 && mod(12,itime["nruns"]) == 0
            nruns = round(Int,itime["nruns"]/12)
            for irun = 1 : nruns;
                iscaanalysis30daily(imod,ipar,itime,iroot,irun);
            end
        else
            for irun = 1 : itime["nruns"];
                iscaanalysisgendaily(imod,ipar,itime,iroot,irun);
            end
        end

    end

end

function iscaanalysis360daily(
    imod::AbstractDict, ipar::AbstractDict, itime::AbstractDict, iroot::AbstractDict,
    irun::Integer
)

    rds,rvar = iscarawread(ipar,iroot,irun=irun);
    attr = Dict{AbstractString,Any}();
    lon = rds["lon"][:]*1; nlon = length(lon)
    lat = rds["lat"][:]*1; nlat = length(lat)
    lvl = ipar["level"];

    if lvl != "sfc"; raw = rvar[:,:,lvl,:]*1; else; raw = rvar[:]*1; end; close(rds);
    raw = reshape(raw,nlon,nlat,:,12)

    davg = zeros(Float32,nlon,nlat,13); dstd = zeros(Float32,nlon,nlat,13);
    dmax = zeros(Float32,nlon,nlat,13); dmin = zeros(Float32,nlon,nlat,13);

    zavg = zeros(Float32,nlat,13); zstd = zeros(Float32,nlat,13);
    zmax = zeros(Float32,nlat,13); zmin = zeros(Float32,nlat,13);

    mavg = zeros(Float32,nlon,13); mstd = zeros(Float32,nlon,13);
    mmax = zeros(Float32,nlon,13); mmin = zeros(Float32,nlon,13);

    for mo = 1 : 12

        @info "$(Dates.now()) - Analyzing $(uppercase(ipar["name"])) data for MONTH $mo of RUN $irun ..."
        rawii = @view raw[:,:,:,mo]

        @debug "$(Dates.now()) - Extracting information on monthly climatology ..."
        davg[:,:,mo] = mean(rawii,dims=3);
        dstd[:,:,mo] = std(rawii,dims=3);
        dmax[:,:,mo] = maximum(rawii,dims=3);
        dmin[:,:,mo] = minimum(rawii,dims=3);

    end

    @info "$(Dates.now()) - Calculating yearly climatology for $(uppercase(ipar["name"])) for RUN $irun ..."
    davg[:,:,end] = mean(davg[:,:,1:12],dims=3);
    dstd[:,:,end] = mean(dstd[:,:,1:12],dims=3);
    dmax[:,:,end] = maximum(dmax[:,:,1:12],dims=3);
    dmin[:,:,end] = minimum(dmin[:,:,1:12],dims=3);

    @info "$(Dates.now()) - Calculating zonal-averaged climatology for $(uppercase(ipar["name"])) for RUN $irun ..."
    for imo = 1 : 12, ilat = 1 : nlat
        zavg[ilat,imo] = iscananmean(@view davg[:,ilat,imo]);
        zstd[ilat,imo] = iscananmean(@view dstd[:,ilat,imo]);
        zmax[ilat,imo] = iscananmean(@view dmax[:,ilat,imo]);
        zmin[ilat,imo] = iscananmean(@view dmin[:,ilat,imo]);
    end

    @info "$(Dates.now()) - Calculating meridional-averaged climatology for $(uppercase(ipar["name"])) for RUN $irun ..."
    for imo = 1 : 12, ilon = 1 : nlon
        mavg[ilon,imo] = iscananmean(@view davg[ilon,:,imo]);
        mstd[ilon,imo] = iscananmean(@view dstd[ilon,:,imo]);
        mmax[ilon,imo] = iscananmean(@view dmax[ilon,:,imo]);
        mmin[ilon,imo] = iscananmean(@view dmin[ilon,:,imo]);
    end

    iscaanasavedaily(
        [davg,dstd,dmax,dmin],
        [zavg,zstd,zmax,zmin],
        [mavg,mstd,mmax,mmin],
        imod,ipar,iroot,irun=irun
    )

end

function iscaanasavedaily(
    data::Array{Array{Float32,3},1},
    zdata::Array{Array{Float32,2},1},
    mdata::Array{Array{Float32,2},1},
    imod::AbstractDict, ipar::AbstractDict, iroot::AbstractDict;
    irun::Integer
)

    @info "$(Dates.now()) - Saving analysed $(uppercase(ipar["name"])) for RUN $irun ..."

    afol = iscaanafolder(ipar,iroot);
    fana = iscaananame(ipar,irun=irun);
    afnc = joinpath(afol,fana);

    if isfile(afnc)
        @info "$(Dates.now()) - Stale NetCDF file $(afnc) detected.  Overwriting ..."
        rm(afnc);
    end

    @debug "$(Dates.now()) - Creating NetCDF file $(afnc) for analyzed $(uppercase(ipar["name"])) for RUN $irun ..."

    ds = Dataset(afnc,"c");
    ds.dim["longitude"] = length(imod["lon"])
    ds.dim["latitude"] = length(imod["lat"])
    ds.dim["month"] = 12;

    nclon = defVar(ds,"lon", Float64, ("longitude",), attrib = Dict(
        "long_name"      => "longitude",
        "units"          => "degrees_E",
        "cartesian_axis" => "X",
        "edges"          => "lonb",
    ))

    nclat = defVar(ds,"lat", Float64, ("latitude",), attrib = Dict(
        "long_name"      => "latitude",
        "units"          => "degrees_N",
        "cartesian_axis" => "Y",
        "edges"          => "latb",
    ))

    nclon[:] = imod["lon"]
    nclat[:] = imod["lat"]

    attr_var = Dict(
        "long_name" => ipar["name"],
        "units"     => ipar["unit"],
    );

    @debug "$(Dates.now()) - Saving analyzed $(uppercase(ipar["name"])) data for RUN $irun to NetCDF file $(afnc) ..."

    v = defVar(ds,"domain_yearly_mean_climatology",Float32,
               ("longitude","latitude"),attrib=attr_var);
    v[:] = data[1][:,:,end];

    v = defVar(ds,"domain_yearly_std_climatology",Float32,
               ("longitude","latitude"),attrib=attr_var);
    v[:] = data[2][:,:,end];

    v = defVar(ds,"domain_yearly_maximum_climatology",Float32,
               ("longitude","latitude"),attrib=attr_var);
    v[:] = data[3][:,:,end];

    v = defVar(ds,"domain_yearly_minimum_climatology",Float32,
               ("longitude","latitude"),attrib=attr_var);
    v[:] = data[4][:,:,end];

    v = defVar(ds,"domain_monthly_mean_climatology",Float32,
               ("longitude","latitude","month"),attrib=attr_var);
    v[:] = data[1][:,:,1:12];

    v = defVar(ds,"domain_monthly_std_climatology",Float32,
               ("longitude","latitude","month"),attrib=attr_var);
    v[:] = data[2][:,:,1:12];

    v = defVar(ds,"domain_monthly_maximum_climatology",Float32,
               ("longitude","latitude","month"),attrib=attr_var);
    v[:] = data[3][:,:,1:12];

    v = defVar(ds,"domain_monthly_minimum_climatology",Float32,
               ("longitude","latitude","month"),attrib=attr_var);
    v[:] = data[4][:,:,1:12];

    v = defVar(ds,"zonalavg_yearly_mean_climatology",Float32,
               ("latitude",),attrib=attr_var);
    v[:] = zdata[1][:,end];

    v = defVar(ds,"zonalavg_yearly_std_climatology",Float32,
               ("latitude",),attrib=attr_var);
    v[:] = zdata[2][:,end];

    v = defVar(ds,"zonalavg_yearly_maximum_climatology",Float32,
               ("latitude",),attrib=attr_var);
    v[:] = zdata[3][:,end];

    v = defVar(ds,"zonalavg_yearly_minimum_climatology",Float32,
               ("latitude",),attrib=attr_var);
    v[:] = zdata[4][:,end];

    v = defVar(ds,"zonalavg_monthly_mean_climatology",Float32,
               ("latitude","month"),attrib=attr_var);
    v[:] = zdata[1][:,1:12];

    v = defVar(ds,"zonalavg_monthly_std_climatology",Float32,
               ("latitude","month"),attrib=attr_var);
    v[:] = zdata[2][:,1:12];

    v = defVar(ds,"zonalavg_monthly_maximum_climatology",Float32,
               ("latitude","month"),attrib=attr_var);
    v[:] = zdata[3][:,1:12];

    v = defVar(ds,"zonalavg_monthly_minimum_climatology",Float32,
               ("latitude","month"),attrib=attr_var);
    v[:] = zdata[4][:,1:12];

    v = defVar(ds,"meridionalavg_yearly_mean_climatology",Float32,
               ("longitude",),attrib=attr_var);
    v[:] = mdata[1][:,end];

    v = defVar(ds,"meridionalavg_yearly_std_climatology",Float32,
               ("longitude",),attrib=attr_var);
    v[:] = mdata[2][:,end];

    v = defVar(ds,"meridionalavg_yearly_maximum_climatology",Float32,
               ("longitude",),attrib=attr_var);
    v[:] = mdata[3][:,end];

    v = defVar(ds,"meridionalavg_yearly_minimum_climatology",Float32,
               ("longitude",),attrib=attr_var);
    v[:] = mdata[4][:,end];

    v = defVar(ds,"meridionalavg_monthly_mean_climatology",Float32,
               ("longitude","month"),attrib=attr_var);
    v[:] = mdata[1][:,1:12];

    v = defVar(ds,"meridionalavg_monthly_std_climatology",Float32,
               ("longitude","month"),attrib=attr_var);
    v[:] = mdata[2][:,1:12];

    v = defVar(ds,"meridionalavg_monthly_maximum_climatology",Float32,
               ("longitude","month"),attrib=attr_var);
    v[:] = mdata[3][:,1:12];

    v = defVar(ds,"meridionalavg_monthly_minimum_climatology",Float32,
               ("longitude","month"),attrib=attr_var);
    v[:] = mdata[4][:,1:12];

    close(ds);

    @info "$(Dates.now()) - Analysed $(uppercase(ipar["name"])) data for RUN $irun has been saved into file $(afnc) and moved to the data directory $(afol)."

end


text:
"""
ap_test_time(re, factors; cov_method = "time iid", lag = 0, test_method = "asymptotic")

TIME-SERIES ASSET PRICING TESTS

Asset pricing models predict a restriction on the intercepts of the time-series regressions.
The asset pricing model is valid if the intercepts of time-series regressions are jointly zero.

Regress asset excess returns on factors to obtain betas

INPUTS
`re': T x N matrix of excess returns, where T is the number of periods for each test assets and N is the number of test assets
`factors': T x K matrix of factors. Factors must be excess returns for time-series tests
`cov_method': method for calculating the variance-covariance matrix of the error terms:
            - time iid: error terms are iid across time but cross-sectionally correlated
            - iid: error terms are iiid both across time and assets
            - NW: Newey-West standard errors with m lags
            - HH: Hansen-Hodrick standard errors with m lag
`lag': number of lags used in the Newey-West or Hansen-Hodrick standard errors
`test_method': method of the time-series asset pricing tests:
            - asymptotic: OLS asymptotic distribution
            - finite: OLS finite-sample distribution - Gibbons Ross and Shanken (1989) test statistics
            - GMM asymptotic: asymptotic GMM exploiting the moment conditions in time-series regressions
            - GMM finite: finite-distribution GMM exploiting the moment conditions in time-series regressions

OUTPUTS
`alpha': N x 1 vector of alphas of the test assets
`beta': N x K matrix of the betas of the test assets. There are N test assets and K factors
`resid': T x N matrix of the pricing errors
`cov_resid': N x N matrix of covariance matrix of the pricing errors

`test_statistics': Wald-statistics of the asset pricing tests
`pvalue': pvalue of the test statistics. For asymptotics, the test statistics has a Chi-square distribution. For finite sample, it has a F distribution

"""

function ap_test_time(re, factors, cov_method = "time iid", lag = 0, test_method = "asymptotic")
    (T, K) = size(factors);
    N = size(re, 2);
    Ef = mean(factors, dims = 1)'; # mean values of factors K x 1
    cov_f = cov(factors); # coviance matrix for the factors K x K

    X = [ones(T, 1) factors]; # T x (K+1)
    (alpha, beta, resid, cov_resid) = OLSFn(re, X, cov_method, lag); # time-series regressions

    # Methods of time series test
    if test_method == "asymptotic" # OLS Asymptotic
        test_statistics = T * ((1.0 .+ Ef' * (cov_f \ Ef)) \ (alpha')) * (cov_resid \ alpha);
        pvalue = 1.0 .- cdf.(Chisq(N), test_statistics);

    elseif test_method == "finite" # OLS Finite (Gibbons Ross and Shanken (1989))
        test_statistics = ((T - N - K) / N) * ((1.0 .+ Ef' * (cov_f \ Ef)) \ (alpha')) * (cov_resid \ alpha);
        pvalue = 1.0 .- cdf.(FDist(N, T - N - K), test_statistics);

    elseif test_method == "GMM asymptotic" || test_method == "GMM finite" # GMM Asymptotic
        h = [ones(T, 1) factors];
        f = kron(ones(1, K+1), resid) .* kron(h, ones(1, N)); # Moment Conditions - p.218 in Cochrane (2000)
        d = - kron(h' * h / T, I(N)); # Jacobian of the moment conditions
        dinv = d \ I; # Inverse of the Jacobian

        Ef2 = CovMoment(f, cov_method, lag); # spectral density matrix
        vp = (1.0 / T) * dinv * Ef2 * dinv'; # asymptotic covariance of the regression coefficiets - p.218 equation (172) in Cochrane
        wt = vp[1:N, 1:N] \ I; # covariance matrix of the intercepts (alphas)
        test_statistics = alpha' * wt * alpha;
        pvalue = 1.0 .- cdf.(Chisq(N), test_statistics);

        if test_method == "GMM finite" # GMM Finite
            test_statistics = test_statistics * ((T - N - K) / N) * ((1.0 .+ Ef' * (cov_f \ Ef)) \ I) / (T * ((1.0 .+ Ef' * (cov_f \ Ef)) \ I));
            pvalue = 1.0 .- cdf.(FDist(N, T - N - K), test_statistics);
        end
    end

    return (alpha, beta, resid, cov_resid, test_statistics, pvalue)
end


text:
"""
    runserver(pipe_in=stdin, pipe_out=stdout[, env_path])

Run a `LanguageServerInstance` reading from `pipe_in` and writing to `pipe_out`.

The same options can be passed to `runserver` as to
[`LanguageServerInstance`](@ref). If `env_path` is not specified,
attempt to pick an environment by considering in order of priority:

1. [`ARGS`[1]](@ref): the first command-line argument passed to the
   invocation of `julia`.
2. The Julia project containing [`pwd()`](@ref).
3. The default Julia environment withing `.julia/environments/v#.#`.

# Examples

The following invocation of Julia would set `env_path` to
`/home/example/repos/Example.jl`:

```sh
julia --project=/path/to/LanguageServer.jl \\
  -e "using LanguageServer; runserver()" \\
  /home/example/repos/Example.jl
```

If there was a `Project.toml` or `JuliaProject.toml` in
`/home/example/repos/Example.jl/`, the following invocation would set
`env_path` to `/home/example/repos/Example.jl/`; otherwise it would be
set to `.julia/environments/v#.#` where `v#.#` is the major/minor
version of Julia being invoked.

```sh
julia --project=/path/to/LanguageServer.jl \\
  -e "using LanguageServer; runserver()"
```
"""
function runserver(pipe_in=stdin, pipe_out=stdout, env_path=choose_env(),
                   depot_path="", err_handler=nothing, symserver_store_path=nothing)
    server = LanguageServerInstance(pipe_in, pipe_out, env_path, depot_path,
                                    err_handler, symserver_store_path)
    run(server)
end

function choose_env()
    maybe_dirname = x -> x !== nothing ? dirname(x) : nothing
    something(get(ARGS, 1, nothing),                            # 1. path passed explicitly
              maybe_dirname(Base.current_project(pwd())),       # 2. parent project of pwd()
              maybe_dirname(Base.load_path_expand("@v#.#")))    # 3. default "global" env
end


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

mutable struct UnprepareNetworkPoliciesRequest <: SwaggerModel
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName

    function UnprepareNetworkPoliciesRequest(;serviceName=nothing)
        o = new()
        validate_property(UnprepareNetworkPoliciesRequest, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        o
    end
end # type UnprepareNetworkPoliciesRequest

const _property_map_UnprepareNetworkPoliciesRequest = Dict{Symbol,Symbol}(Symbol("serviceName")=>Symbol("serviceName"))
const _property_types_UnprepareNetworkPoliciesRequest = Dict{Symbol,String}(Symbol("serviceName")=>"String")
Base.propertynames(::Type{ UnprepareNetworkPoliciesRequest }) = collect(keys(_property_map_UnprepareNetworkPoliciesRequest))
Swagger.property_type(::Type{ UnprepareNetworkPoliciesRequest }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UnprepareNetworkPoliciesRequest[name]))}
Swagger.field_name(::Type{ UnprepareNetworkPoliciesRequest }, property_name::Symbol) =  _property_map_UnprepareNetworkPoliciesRequest[property_name]

function check_required(o::UnprepareNetworkPoliciesRequest)
    true
end

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


text:
@inline function initialize!(integrator,cache::DiscreteConstantCache,f=integrator.f)
  integrator.kshortsize = 0
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function initialize!(integrator,cache::DiscreteCache,f=integrator.f)
  integrator.kshortsize = 0
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::DiscreteConstantCache,f=integrator.f)
  if discrete_apply_map(integrator.alg)
    if discrete_scale_by_time(integrator.alg)
      integrator.u = integrator.uprev + integrator.dt*f(integrator.t,integrator.uprev)
    else
      integrator.u = f(integrator.t,integrator.uprev)
    end
  end
end

@inline function perform_step!(integrator,cache::DiscreteCache,f=integrator.f)
  @unpack u,uprev,dt,t = integrator
  @unpack du = cache
  if discrete_apply_map(integrator.alg)
    if discrete_scale_by_time(integrator.alg)
      f(t,uprev,du)
      for i in eachindex(integrator.u)
        u[i] = @muladd uprev[i] + dt*du[i]
      end
    else
      f(t,uprev,u)
    end
  end
end

@inline function initialize!(integrator,cache::EulerConstantCache,f=integrator.f)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::EulerConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  k = integrator.fsalfirst
  u = muladd(dt,k,uprev)
  k = f(t+dt,u) # For the interpolation, needs k at the updated point
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::EulerCache,f=integrator.f)
  integrator.kshortsize = 2
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # For the interpolation, needs k at the updated point
end

@inline function perform_step!(integrator,cache::EulerCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  for i in uidx
    u[i] = muladd(dt,integrator.fsalfirst[i],uprev[i])
  end
  f(t+dt,u,integrator.fsallast) # For the interpolation, needs k at the updated point
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::MidpointConstantCache,f=integrator.f)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::MidpointConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  halfdt = dt/2
  k = integrator.fsalfirst
  k = f(t+halfdt,uprev+halfdt*k)
  u = uprev + dt*k
  integrator.fsallast = f(t+dt,u) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::MidpointCache,f=integrator.f)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # FSAL for interpolation
end

@inline function perform_step!(integrator,cache::MidpointCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k,du,tmp,fsalfirst = cache
  halfdt = dt/2
  for i in uidx
    tmp[i] = muladd(halfdt,integrator.fsalfirst[i],uprev[i])
  end
  f(t+halfdt,tmp,du)
  for i in uidx
    u[i] = muladd(dt,du[i],uprev[i])
  end
  f(t+dt,u,k)
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::RK4ConstantCache,f=integrator.f)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::RK4ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  halfdt = dt/2
  k₁ =integrator.fsalfirst
  ttmp = t+halfdt
  k₂ = f(ttmp,muladd(halfdt,k₁,uprev))
  k₃ = f(ttmp,muladd(halfdt,k₂,uprev))
  k₄ = f(t+dt,muladd(dt,k₃,uprev))
  u = muladd(dt/6,muladd(2,(k₂ + k₃),k₁+k₄),uprev)
  k = f(t+dt,u)
  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::RK4Cache,f=integrator.f)
  @unpack tmp,fsalfirst,k₂,k₃,k₄,k = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # pre-start FSAL
end

@inline function perform_step!(integrator,cache::RK4Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack tmp,fsalfirst,k₂,k₃,k₄,k = cache
  k₁ = fsalfirst
  halfdt = dt/2
  ttmp = t+halfdt
  for i in uidx
    tmp[i] = muladd(halfdt,k₁[i],uprev[i])
  end
  f(ttmp,tmp,k₂)
  for i in uidx
    tmp[i] = muladd(halfdt,k₂[i],uprev[i])
  end
  f(ttmp,tmp,k₃)
  for i in uidx
    tmp[i] = muladd(dt,k₃[i],uprev[i])
  end
  f(t+dt,tmp,k₄)
  for i in uidx
    u[i] = muladd(dt/6,muladd(2,(k₂[i] + k₃[i]),k₁[i] + k₄[i]),uprev[i])
  end
  f(t+dt,u,k)
  @pack integrator = t,dt,u
end


text:
"""
    soft_theta(x::T, k=T(1)) where T

This is a version of the theta function that uses a soft transition and is differentiable.
It is implemented as `1 / (1+exp(-2 * k * x))`.
A larger `k` makes the edge harder.
"""
function soft_theta(x::T, k=1) where T
    one(T) / (one(T) + exp(-T(2)*k*x))
end

"""
    soft_theta_pw(x::T, ϵ=T(0.01)) where T

A different version of the soft theta. Uses a `cos` instead of `exp` functions.
Furthermore, is piecewise defined.
A smaller `e` makes the edge harder.
"""
function soft_theta_pw(x::T, ϵ=T(0.01)) where T
    if x .> ϵ
        one(T)
    elseif x .< -ϵ
        zero(T)
    else
        (one(T) - cos((x + ϵ) * (T(π) / (T(2)*ϵ)))) / T(2)
    end
end

"""
    soft_delta(x::T, k=T(1)) where T

This is a smooth version of the delta function that uses a soft peak and is differentiable.
Based on a Gauss function.
A larger `k` makes the edge harder.
"""
function soft_delta(x::T, k=T(1)) where T
    one(T) / (k * √(T(π))) * exp(-abs2(x / k))
end

"""
    soft_delta_pw(x::T, ϵ=T(1)) where T

This is a smooth version of the delta function that uses a soft peak and is differentiable.
Based on a Gauss function.
A larger `k` makes the edge harder.
"""
function soft_delta_pw(x::T, ϵ=T(0.01)) where T
    if abs(x) > abs(ϵ)
        zero(T)
    else
        (one(T) + cos(x * (T(π) / ϵ))) / 2
    end
end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
using DataFrames, Dash, DashHtmlComponents, DashCoreComponents, UrlDownload, PlotlyJS

df6 = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/country_indicators.csv"))

rename!(df6, Dict(:"Year" => "year"))

dropmissing!(df6)

available_indicators = unique(df6[:, "Indicator Name"])
years = unique(df6[:, "year"])

app = dash()

app.layout = html_div() do
    html_div(
        children = [
            html_div(
                children = [
                    dcc_dropdown(
                        id = "crossfilter-xaxis-column",
                        options = [
                            (label = i, value = i)
                            for i in available_indicators
                        ],
                        value = "Fertility rate, total (births per woman)",
                    ),
                    dcc_radioitems(
                        id = "crossfilter-xaxis-type",
                        options = [
                            (label = i, value = i) for i in ["linear", "log"]
                        ],
                        value = "linear",
                    ),
                ],
                style = (width = "49%", display = "inline-block"),
            ),
            html_div(
                children = [
                    dcc_dropdown(
                        id = "crossfilter-yaxis-column",
                        options = [
                            (label = i, value = i)
                            for i in available_indicators
                        ],
                        value = "Life expectancy at birth, total (years)",
                    ),
                    dcc_radioitems(
                        id = "crossfilter-yaxis-type",
                        options = [
                            (label = i, value = i) for i in ["linear", "log"]
                        ],
                        value = "linear",
                    ),
                ],
                style = (
                    width = "49%",
                    float = "right",
                    display = "inline-block",
                ),
            ),
        ],
        style = (
            borderBottom = "thin lightgrey solid",
            backgroundColor = "rgb(250, 250, 250)",
            padding = "10px 5px",
        ),
    ),
    html_div(
        children = [
            dcc_graph(id = "crossfilter-indicator-scatter"),
            dcc_slider(
                id = "crossfilter-year-slider",
                min = minimum(years),
                max = maximum(years),
                marks = Dict([Symbol(v) => Symbol(v) for v in years]),
                value = minimum(years),
                step = nothing,
            ),
        ],
        style = (
            width = "49%",
            display = "inline-block"
        ),
    ),
    html_div(
        children = [
            dcc_graph(id = "x-time-series"),
            dcc_graph(id = "y-time-series"),
        ],
        style = (width = "49%", display = "inline-block"),
    )
end

callback!(
    app,
    Output("crossfilter-indicator-scatter", "figure"),
    Input("crossfilter-xaxis-column", "value"),
    Input("crossfilter-yaxis-column", "value"),
    Input("crossfilter-xaxis-type", "value"),
    Input("crossfilter-yaxis-type", "value"),
    Input("crossfilter-year-slider", "value"),
) do xaxis_column_name, yaxis_column_name, xaxis_type, yaxis_type, year_slider_value

    df6f = df6[df6.year .== year_slider_value, :]

    return Plot(
        df6f[df6f[Symbol("Indicator Name")] .== xaxis_column_name, :Value],
        df6f[df6f[Symbol("Indicator Name")] .== yaxis_column_name, :Value],
        Layout(
            xaxis_type = xaxis_type == "Linear" ? "linear" : "log",
            xaxis_title = xaxis_column_name,
            yaxis_title = yaxis_column_name,
            yaxis_type = yaxis_type == "Linear" ? "linear" : "log",
            hovermode = "closest",
            height = 450,
        ),
        kind = "scatter",
        text = df6f[
            df6f[Symbol("Indicator Name")] .== yaxis_column_name,
            Symbol("Country Name"),
        ],
        customdata = df6f[
            df6f[Symbol("Indicator Name")] .== yaxis_column_name,
            Symbol("Country Name"),
        ],
        mode = "markers",
        marker_size = 15,
        marker_opacity = 0.5,
        marker_line_width = 0.5,
        marker_line_color = "white",
    )
end

function create_time_series(df6f, axis_type, title)
    Plot(
        df6f[:, :year],
        df6f[:, :Value],
        Layout(
            yaxis_type = axis_type == "Linear" ? "linear" : "log",
            xaxis_showgrid = false,
            annotations = [attr(
                x = 0,
                y = 0.85,
                xanchor = "left",
                yanchor = "bottom",
                xref = "paper",
                yref = "paper",
                showarrow = false,
                align = "left",
                bgcolor = "rgba(255, 255, 255, 0.5)",
                text = title,
            )],
            height = 225,
        ),
        mode = "lines+markers",
    )
end

callback!(
    app,
    Output("x-time-series", "figure"),
    Input("crossfilter-indicator-scatter", "hoverData"),
    Input("crossfilter-xaxis-column", "value"),
    Input("crossfilter-xaxis-type", "value"),
) do hover_data, xaxis_column_name, axis_type
    country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
    df6f = df6[df6[:, Symbol("Country Name")].==country_name, :]
    df6f = df6f[df6f[:, Symbol("Indicator Name")].==xaxis_column_name, :]
    title = "$(country_name)
$(xaxis_column_name)"
    return create_time_series(df6f, axis_type, title)

end

callback!(
    app,
    Output("y-time-series", "figure"),
    Input("crossfilter-indicator-scatter", "hoverData"),
    Input("crossfilter-yaxis-column", "value"),
    Input("crossfilter-yaxis-type", "value"),
) do hover_data, yaxis_column_name, axis_type
    country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
    df6f = df6[df6[:, Symbol("Country Name")].==country_name, :]
    df6f = df6f[df6f[:, Symbol("Indicator Name")].==yaxis_column_name, :]
    return create_time_series(df6f, axis_type, yaxis_column_name)
end

run_server(app, "0.0.0.0", debug = true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
"""
	Utility functions for MRForecaster.jl
"""

"""
	split_hyper_linear(hyper)

Splits a set of hyperparameters into an array of parameters
for a broken power law.

Returns an array of intercepts, slopes, dispersions and transition points.

"""
function split_hyper_linear(hyper::AbstractVector)
	c0, slope, sigma, trans =
	hyper[1], hyper[2:1+N_POPS], hyper[2+N_POPS:1+2*N_POPS],
	hyper[2+2*N_POPS:end]

	# Compute individual intercepts
	c = Vector{Float64}(undef, N_POPS)
	c[1] = c0
	for i = 2:N_POPS
		c[i] = c[i-1] + trans[i-1]*(slope[i-1]-slope[i])
	end

	return c, slope, sigma, trans
end

"""
	piece_linear(hyper, M, prob_R)

Generate a radius from the given mass M from the piecewise model
specified by the hyperparameters in hyper, with a dispersion
parameter of prob_R.
"""
function piece_linear(hyper, M, prob_R)
	c, slope, sigma, trans = split_hyper_linear(hyper)
	if !isa(M, AbstractVector)
		mval = M
		M = Vector{Float64}(undef, 1)
		M[1] = mval
		prob_rval = prob_R
		prob_R = Vector{Float64}(undef, 1)
		prob_R[1] = prob_rval
	end
	R = similar(M)
	for i in 1:N_POPS
		# Determine the masses in population i
		ind = indicate(M, trans, i)
		mu = c[i] .+ (M[ind] .* slope[i])
		R[ind] = quantile.(Normal.(mu, sigma[i]), prob_R[ind])
	end
	return R
end

# Preallocate vectors
const trans_compare = Vector{Float64}(undef, N_POPS+1)
trans_compare[1] = -Inf
trans_compare[end] = Inf
"""
	indicate(M, trans, i)

Determine the masses in the array M that belong to
population i, as specified by the transition parameter trans.
"""
function indicate(M, trans, i)
	# Create array for comparison
	trans_compare[2:N_POPS] = trans
	ind = (M .> trans_compare[i]) .& (M .<= trans_compare[i+1])
	return ind
end

"""
	prob_R_given_M(radii, M, hyper)

Given an array of masses and set of hyperparameters,
return the probability of the radius.
"""
function prob_R_given_M(radii, M, hyper)
	# Get hyperparameteers
	c, slope, sigma, trans = split_hyper_linear(hyper)

	prob = similar(M)
	for i in 1:N_POPS
		# Determine the masses in population i
		ind = indicate(M, trans, i)
		mu = c[i] .+ (M[ind] .* slope[i])
		prob[ind] = pdf.(Normal.(mu, sigma[i]), radii)
	end

	prob /= sum(prob)
	return prob
end


text:
const W_per_cm2 = typeof(1.0u"W/cm^2")
const kg_per_m3 = typeof(1.0u"kg/m^3")
const m_per_s = typeof(1.0u"m" / 1.0u"s")
const hertz = typeof(1.0u"Hz")
const seconds = typeof(1.0u"s")

const PressureArray{x} = Array{T, x} where T <: Union{Unitful.Pressure, Number, W_per_cm2}

squeeze(A::PressureArray{1}) = A[:]
squeeze(A::PressureArray{2}) = A[1, :]
squeeze(A::PressureArray{3}) = A[1, :, :]
squeeze(A::PressureArray{4}) = A[1, :, :, :]

"""
struct Medium
    density: Density [kg/m³]
          c: Speed of sound [m/s]
end
"""
struct Medium
    density::kg_per_m3          # [kg/m³]
    c::m_per_s                  # [m/s]
end

Medium() = Medium(1000.0u"kg/m^3", 1480.0u"m/s")
#TODO: make a selector for 20C water and standard human tissue.

"""
struct Excitation
                f0: Center frequency [Hz]
    pulse_duration: On-time / duration of the excitation pulse [s]
        duty_cycle: [0, 1]
    total_duration: Total duration of waveform [s]
                fs: Sampling frequency of the received signal [Hz]
end
"""
struct Excitation
    f0::hertz                   # [Hz]
    pulse_duration::seconds     # [s]
    duty_cycle::Float64
    total_duration::seconds     # [s]
    fs::hertz                   # [Hz]
end

Excitation() = Excitation(5e6u"Hz", 0.1u"s", 1.0, 10.0u"s", 1.0u"Hz")
Excitation(f0, pulse_duration::Number, duty_cycle) = Excitation(f0, Float64(pulse_duration), Float64(duty_cycle), 10u"s", 1.0u"Hz")


text:
""" Useful tools for analysing and comparing between multiple PEtriNet models
"""

module ModelComparison

using AlgebraicPetri
using Catlab
using Catlab.CategoricalAlgebra
import Catlab.CategoricalAlgebra: Subobject

export petri_homomorphisms, compare, PetriSubobject

""" petri_homomorphisms
Calculates all homomorphisms between two PetriNet models and returns an array
of ACSetTransformations that describe these homomorphisms.

NOTE:
This function does restrict to homomorphisms which preserve the transition
signatures (number of input/output wires).
"""
function petri_homomorphisms(p1::AbstractPetriNet, p2::AbstractPetriNet; kw...)
  results = ACSetTransformation[]
  sigsTS = Set{Vector{Int64}}()
  homs = homomorphisms(PetriNet(p1), PetriNet(p2); kw...)
  for transform in homs
    if all([length(inputs(p1, t)) == length(inputs(p2, transform.components[:T](t))) &&
            length(outputs(p1, t)) == length(outputs(p2, transform.components[:T](t)))
            for t in 1:nt(p1)]) &&
        vcat(transform.components[:T].func, transform.components[:S].func) ∉ sigsTS

      push!(results, ACSetTransformation(deepcopy(transform.components), p1, p2))
      push!(sigsTS, vcat(transform.components[:T].func,
                         transform.components[:S].func))
    end
  end
  results
end

""" compare
Calculates all homomorphisms from a single Petri net (`apex`) to other Petri
nets. This is returned as a Multispan with `apex` as the apex and the
homomorphisms as the legs.
"""
function compare(apex::T, pn::T...; monic=[:S, :T, :I, :O], kw...) where {T <: AbstractPetriNet}
  legs = map(enumerate(pn)) do (i, p)
    hp = petri_homomorphisms(apex, p; monic=monic, kw...)
    isempty(hp) && error("No homomorphism from apex to leg $i which is monic in $monic")
    hp
  end
  Multispan(apex, vcat(legs...))
end

end

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
facts("Frequency Adaptation") do

  context("returns all indices once in the first block") do
    for(reps in 1:20)
      n = rand(1:50)
      fa = BlackBoxOptim.FrequencyAdapter(n)
      block = Int64[]
      for(i in 1:n)
        mi = BlackBoxOptim.next(fa)
        @fact (1  <= mi <= n) => true
        push!(block, mi)
        BlackBoxOptim.update!(fa, mi, rand())
      end
      @fact sort(block) => collect(1:n)
    end
  end

  context("increases the frequency of a method that has higher progress values") do
    for(reps in 1:20)
      n = rand(2:50)
      fa = BlackBoxOptim.FrequencyAdapter(n)
      up1(mi) = BlackBoxOptim.update!(fa, mi, rand() + ((mi == 1) ? 0.4 : 0.0))
      counts = zeros(Int64, n)
      uprep(rs) = begin
        for(i in 1:rs)
          mi = BlackBoxOptim.next(fa)
          counts[mi] += 1
          up1(mi)
        end
      end
      uprep(n)
      @fact counts == ones(Int64, n) => true
      uprep(20*n)
      @fact counts[1] > (sum(counts[2:end])/(n-1)) => true
    end
  end
end


text:
function (S::NLAnderson{false,<:NLSolverCache})(integrator)
  nlcache = S.cache
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,tmp,W,κ,tol,c,γ,max_iter,min_iter,zs,gs = nlcache
  mass_matrix = integrator.f.mass_matrix
  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
  else
    f = integrator.f
  end
  # precalculations
  κtol = κ*tol

  # zs = zeros(S.n+1)
  # gs = zeros(S.n+1)
  residuals = fill(zero(eltype(z)), length(z), S.n+1)
  alphas = fill(zero(eltype(z)), S.n)
  # initial step of NLAnderson iteration
  zs[1] = z
  iter = 1
  tstep = t + c*dt
  u = @. tmp + γ*z
  if mass_matrix == I
    z₊ = dt .* f(u, p, tstep)
  else
    mz = mass_matrix * z
    z₊ = dt .* f(u, p, tstep) .- mz .+ z
  end
  gs[1] = z₊
  dz = z₊ .- z
  ndz = integrator.opts.internalnorm(dz)
  zptr = zs[S.n+1]
  gptr = gs[S.n+1]
  for t = (S.n+1):-1:2
    zs[t] = zs[t-1]
    gs[t] = gs[t-1]
  end
  zs[1] = zptr
  gs[1] = gptr
  # zs = circshift(zs, 1)
  # gs = circshift(gs, 1)
  z = z₊
  zs[1] = z
  η = nlcache.ηold
  do_anderson = true

  # anderson acceleration for fixed point iteration
  fail_convergence = false
  while (do_anderson || iter < min_iter) && iter < max_iter
    iter += 1
    u = @. tmp + γ*z
    if mass_matrix == I
      z₊ = dt .* f(u, p, tstep)
    else
      mz = mass_matrix * z
      z₊ = dt .* f(u, p, tstep) .- mz .+ z
    end
    gs[1] = z₊

    mk = min(S.n, iter-1)
    ztmp = zs[1] .- gs[1]
    for i in 2:mk+1
      @. residuals[:,i-1] = gs[i] - zs[i] + ztmp
    end
    resqr = qr(residuals[:,1:mk], Val(true))
    if typeof(ztmp) <: AbstractArray
      alphas[1:mk] = resqr \ ztmp
    else
      alphas[1:mk] = resqr \ [ztmp]
    end
    for i = 1:mk
        z₊ += alphas[i].*(gs[i+1] .- gs[1])
    end
    zptr = zs[S.n+1]
    gptr = gs[S.n+1]
    for t = (S.n+1):-1:2
      zs[t] = zs[t-1]
      gs[t] = gs[t-1]
    end
    zs[1] = zptr
    gs[1] = gptr
    # zs = circshift(zs, 1)
    # gs = circshift(gs, 1)
    zs[1] = z₊
    ndzprev = ndz
    dz = z₊ .- z
    ndz = integrator.opts.internalnorm(dz)
    θ = ndz/ndzprev
    if θ > 1 || ndz*(θ^(max_iter - iter)/(1-θ)) > κtol
      fail_convergence = true
      break
    end
    η = θ/(1-θ)
    do_anderson = (η*ndz > κtol)
    z = z₊
  end

  if (iter >= max_iter && do_anderson) || fail_convergence
    integrator.force_stepfail = true
    return (z, η, iter, true)
  end

  return (z, η, iter, false)
end

function (S::NLAnderson{true})(integrator)
  nlcache = S.cache
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,z₊,dz,b,tmp,κ,tol,k,c,γ,min_iter,max_iter,zs,gs = nlcache
  ztmp = b
  mass_matrix = integrator.f.mass_matrix
  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
  else
    f = integrator.f
  end
  # precalculations
  κtol = κ*tol

  # zs = [zero(vec(z)) for i in 1:S.n+1]
  # gs = [zero(vec(z)) for i in 1:S.n+1]
  residuals = zeros(length(vec(z)), S.n+1)
  alphas = zeros(S.n)
  # initial step of NLAnderson iteration
  zs[1] .= vec(z)
  iter = 1
  tstep = t + c*dt
  @. u = tmp + γ*z
  # z₊ = dt*f(u, p, tstep)
  f(k, u, p, tstep)
  if mass_matrix == I
    @. z₊ = dt*k
  else
    @. z₊ = dt*k + z
    mul!(ztmp, mass_matrix, z)
    @. z₊ -= ztmp
  end
  gs[1] .= vec(z₊)
  @. dz = z₊ - z
  ndz = integrator.opts.internalnorm(dz)
  zptr = zs[S.n+1]
  gptr = gs[S.n+1]
  for t = (S.n+1):-1:2
    zs[t] = zs[t-1]
    gs[t] = gs[t-1]
  end
  zs[1] = zptr
  gs[1] = gptr
  # zs = circshift(zs, 1)
  # gs = circshift(gs, 1)
  z .= z₊
  zs[1] .= vec(z)
  η = nlcache.ηold
  do_anderson = true
  # anderson acceleration for fixed point iteration
  fail_convergence = false
  while (do_anderson || iter < min_iter) && iter < max_iter
    iter += 1
    @. u = tmp + γ*z
    f(k, u, p, tstep)
    if mass_matrix == I
      @. z₊ = dt*k
    else
      ztmp = b
      @. z₊ = dt*k + z
      mul!(ztmp, mass_matrix, z)
      @. z₊ -= ztmp
    end
    gs[1] .= vec(z₊)

    mk = min(S.n, iter-1)
    ztmp = vec(b)
    @. ztmp = zs[1] - gs[1]
    for i in 2:mk+1
      @. residuals[:,i-1] = gs[i] - zs[i] + ztmp
    end
    resqr = qr!(residuals[:,1:mk], Val(true))
    ldiv!(@view(alphas[1:mk]), resqr, ztmp)
    vecz₊ = vec(z₊)
    for i = 1:mk
        vecz₊ .+= alphas[i].*(gs[i+1] .- gs[1])
    end
    zptr = zs[S.n+1]
    gptr = gs[S.n+1]
    for t = (S.n+1):-1:2
      zs[t] = zs[t-1]
      gs[t] = gs[t-1]
    end
    zs[1] = zptr
    gs[1] = gptr
    # zs = circshift(zs, 1)
    # gs = circshift(gs, 1)
    zs[1] .= vec(z₊)
    ndzprev = ndz
    @. dz = z₊ - z
    ndz = integrator.opts.internalnorm(dz)
    θ = ndz/ndzprev
    if θ > 1 || ndz*(θ^(max_iter - iter)/(1-θ)) > κtol
      fail_convergence = true
      break
    end
    η = θ/(1-θ)
    do_anderson = (η*ndz > κtol)
    @. z = z₊
  end

  if (iter >= max_iter && do_anderson) || fail_convergence
    integrator.force_stepfail = true
    return (z, η, iter, true)
  end

  return (z, η, iter, false)
end