text:
@test isregular(rand(4))
@test isregular(rand(4,4,4,4))
@test !isregular(rand(4,5))

@test mirror([1]) == [1]
@test mirror([2,3]) == [2,-3]
@test mirror([2,3,4]) == [2,-3,4]
@test mirror(@SVector([2,3,4])) == @SVector [2,-3,4]
@test mirror([4.9,5,6,7]) == [4.9,-5,6,-7]

@test upsample([1]) == [1,0]
@test upsample([1]; skipfirst=true) == [0,1]
@test upsample([1,2]) == [1,0,2,0]
@test upsample([1,2]; skipfirst=true) == [0,1,0,2]
@test downsample([1,2]) == [1]
@test downsample([1,2]; skipfirst=true) == [2]
@test downsample([1,2,3,4]) == [1,3]
@test downsample([1,2,3,4]; skipfirst=true) == [2,4]

signalsize = 64

@test_throws ArgumentError coefficientcount(ones(signalsize); level=0)
@test coefficientcount(ones(signalsize)) == signalsize
@test coefficientcount(rand(signalsize); threshold=1.01) == 0
@test coefficientcount(rand(signalsize).-2.5; threshold=1.01) == signalsize
@test coefficientcount([0,2,6,7,8]; threshold=6.6) == 2
@test coefficientcount([10,-11,6,7,8.0,-5,-8,0]; threshold=3, level=4) == 3
@test coefficientcount(ones(signalsize,signalsize)) == signalsize * signalsize
@test coefficientcount(rand(signalsize,signalsize); threshold=1.01) == 0
@test coefficientcount([-1 2;3 4]; threshold=2.5) == 2

signal = [1,2]
output = [1,2]
upto = 2
@test partition(signal) == output
@test partition(signal; upto) == output
@test combine(signal) == output
@test combine(signal; upto) == output

signal = [-1,2,3.3,4]
output = [-1,3.3,2,4]
upto = 4
@test partition(signal) == output
@test partition(signal; upto) == output

signal = [-1,3.3,2,4]
output = [-1,2,3.3,4]
@test combine(signal) == output
@test combine(signal; upto) == output

signal = randn(128)
@test signal == combine(partition(signal))

# with strides
signal = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]

stride(s, i, l) = s:i:(s + (l - 1) * i)

@test partition(signal, stride(3, 2, 2))[1:2] == [3,5]
@test partition(signal, stride(3, 2, 4))[1:4] == [3,7,5,9]
@test partition(signal, stride(1, 3, 4))[1:4] == [1,7,4,10]
@test partition(signal, stride(1, 1, 8)) ==  partition(signal; upto=8)

range = stride(3, 2, 2)
@test combine(signal, range)[range] == [1,2]

range = stride(3, 2, 4)
@test combine(signal, range)[range] == [1,3,2,4]

range = stride(1, 3, 4)
@test combine(signal, range)[range] == [1,3,2,4]
@test combine(signal, stride(1, 1, 8)) ==  combine(signal; upto=8)


text:
# This file is for functions that are new functions specifically for working with
# NamedDimsArrays, rather than overloads of existing functions

"""
    rename(nda::NamedDimsArray, names)
    rename(nda::NamedDimsArray, pairs)

Returns a new `NamedDimsArray` with the given dimension `names` or `pairs`
of `old=>new` names.
`rename` outright replaces the names; while still wrapping the same backing array.
Unlike the constructor, it does not require that new names are compatible
with the old names (though you do still need to match the number of dimensions).
"""
rename(nda::NamedDimsArray, names) = NamedDimsArray(parent(nda), names)

function rename(nda::NamedDimsArray, pairs::Vararg{Pair{Symbol,Symbol}})
    names = dimnames(nda)
    new_names = ntuple(i -> _rename(names[i], pairs...), length(names))
    return NamedDimsArray(parent(nda), new_names)
end

"""
    refine_names(x, names)

Refine the names of the dimensions of `x` to match `names`.
This is like [`rename`](ref), but it only affects unnamed dimensions.
I.e. dimensions of a `NamedDimsArray` called `:_`, or any dimensions of an
`AbstractArray` in general.
"""
@inline function refine_names(orig::NamedDimsArray{L}, names::Tuple) where {L}
    new_names = unify_names(names, L)
    return NamedDimsArray{new_names}(parent(orig))
end
@inline refine_names(orig::AbstractArray, names::Tuple) = NamedDimsArray{names}(orig)

"""
    unname(A::NamedDimsArray) -> AbstractArray
    unname(A) === A

Return the input array `A` without any dimension names.

For `NamedDimsArray`s this returns the parent array, equivalent to calling `parent`, but for
anything else it simply returns the input.
"""
unname(x::NamedDimsArray) = parent(x)
unname(x) = x

"""
    dimnames(A) -> Tuple
    dimnames(A, d) -> Symbol

Return the names of all the dimensions of the array `A`,
or just the one for the `d`-th dimension.

Gives wildcards `:_` if this is not a `NamedDimsArray`.
Like `size(A, d)`, it allows `d > ndims(A)`, in this case all the trailing dimension are given the wildcard name (`:_`).
"""
dimnames(::Type{<:NamedDimsArray{L}}) where {L} = L
dimnames(::Type{<:AbstractArray{T,N}}) where {T,N} = ntuple(_ -> :_, N)
dimnames(x::T) where {T<:AbstractArray} = dimnames(T)

function dimnames(AT::Type{<:AbstractArray{T,N}}, d::Integer) where {T,N}
    if 1 <= d <= N
        return dimnames(AT)[d]
    elseif d > N
        return :_
    else
        throw(DimensionMismatch("dimnames: dimension out of range"))
    end
end
dimnames(x::T, d::Integer) where {T<:AbstractArray} = dimnames(T, d)

ext: jl


text:
let
    pomdp = BabyPOMDP()
    up = updater(pomdp)

    bmdp = GenerativeBeliefMDP(pomdp, up)
    b = initialstate(bmdp, Random.GLOBAL_RNG)
end


text:
"""
EMEWS simple noop example using Julia code.
Please refer to the Python example, for more details:
https://github.com/emews/EQ-Py/tree/master/examples/noop
Authors: Carmine Spagnuolo and Giuseppe D'Ambrosio

"""

include("../ext/EQ-JL/eqjl.jl")
using .eqjl
function algo()
    eqjl.OUT_put("1;2;3;4;5;6")
    result = eqjl.IN_get()
    eqjl.OUT_put("33;45;46")
    result = eqjl.IN_get()
    eqjl.OUT_put("FINAL")
end
@async algo()

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
using DimensionalData, Test

@test name(nothing) == ""
@test name(Nothing) == ""
@test dims(1) == nothing
@test_throws ErrorException dims(nothing)
@test refdims(1) == ()


text:
@inline function Phi_e(Cell,s,z,ϕ_tf,D,res0)
   """
   Electrolyte Potential Transfer Function

   Phi_e(CellData,s,z)

   """

κ_eff_Neg = Cell.Const.κ*Cell.Neg.ϵ_e^Cell.Neg.κ_brug
κ_eff_Sep = Cell.Const.κ*Cell.Sep.ϵ_e^Cell.Sep.κ_brug
κ_eff_Pos = Cell.Const.κ*Cell.Pos.ϵ_e^Cell.Pos.κ_brug
σ_eff_Neg = Cell.Neg.σ*Cell.Neg.ϵ_s^Cell.Neg.σ_brug #Effective Conductivity Neg
σ_eff_Pos = Cell.Pos.σ*Cell.Pos.ϵ_s^Cell.Pos.σ_brug #Effective Conductivity Pos

#Defining SOC
θ_neg = Cell.Const.SOC * (Cell.Neg.θ_100-Cell.Neg.θ_0) + Cell.Neg.θ_0
θ_pos = Cell.Const.SOC * (Cell.Pos.θ_100-Cell.Pos.θ_0) + Cell.Pos.θ_0

#Prepare for j0
ce0 = Cell.Const.ce0
cs_max_neg = Cell.Neg.cs_max
cs0_neg = cs_max_neg * θ_neg
cs_max_pos = Cell.Pos.cs_max
cs0_pos = cs_max_pos * θ_pos
α_neg = Cell.Neg.α
α_pos = Cell.Pos.α

#Current Flux Density
if Cell.Const.CellTyp == "Doyle_94"
   κ_neg = Cell.Neg.k_norm/Cell.Neg.cs_max/ce0^(1-α_neg)
   κ_pos = Cell.Pos.k_norm/Cell.Pos.cs_max/ce0^(1-α_pos)
   j0_neg = κ_neg*(ce0*(cs_max_neg-cs0_neg))^(1-α_neg)*cs0_neg^α_neg
   j0_pos = κ_pos*(ce0*(cs_max_pos-cs0_pos))^(1-α_pos)*cs0_pos^α_pos
else
   j0_neg = Cell.Neg.k_norm*(Cell.Const.ce0*cs0_neg*(Cell.Neg.cs_max-cs0_neg))^(1-Cell.Neg.α)
   j0_pos = Cell.Pos.k_norm*(Cell.Const.ce0*cs0_pos*(Cell.Pos.cs_max-cs0_pos))^(1-Cell.Pos.α)
end

#Resistances
Rtot_neg = R*Cell.Const.T/(j0_neg*F^2) + Cell.Neg.RFilm
Rtot_pos = R*Cell.Const.T/(j0_pos*F^2) + Cell.Pos.RFilm
#Rtot_neg = R*Cell.Const.T/(j0_neg*Cell.Const.CC_A*F) + Cell.Neg.RFilm # Ωm²
#Rtot_pos = R*Cell.Const.T/(j0_pos*Cell.Const.CC_A*F) + Cell.Pos.RFilm # Ωm²

#OCP derivative
∂Uocp_neg = Cell.Const.∂Uocp("Neg",θ_neg)/cs_max_neg
∂Uocp_pos = Cell.Const.∂Uocp("Pos",θ_pos)/cs_max_pos

ν_neg = @. Cell.Neg.L*sqrt((Cell.Neg.as/σ_eff_Neg+Cell.Neg.as/κ_eff_Neg)/(Rtot_neg+∂Uocp_neg*(Cell.Neg.Rs/(F*Cell.Neg.Ds))*(tanh(Cell.Neg.β)/(tanh(Cell.Neg.β)-Cell.Neg.β)))) #Condensing Variable - eq. 4.13
ν_neg_∞ = @. Cell.Neg.L*sqrt(Cell.Neg.as*((1/κ_eff_Neg)+(1/σ_eff_Neg))/(Rtot_neg))
ν_pos = @. Cell.Pos.L*sqrt((Cell.Pos.as/σ_eff_Pos+Cell.Pos.as/κ_eff_Pos)/(Rtot_pos+∂Uocp_pos*(Cell.Pos.Rs/(F*Cell.Pos.Ds ))*(tanh(Cell.Pos.β)/(tanh(Cell.Pos.β)-Cell.Pos.β)))) #Condensing Variable - eq. 4.13
ν_pos_∞ = @. Cell.Pos.L*sqrt(Cell.Pos.as*((1/κ_eff_Pos)+(1/σ_eff_Pos))/(Rtot_pos))

i=Int64(1)
# Loop Tf's
   for pt in z
         if pt <= Cell.Neg.L+eps()
            ϕ_tf[i,:] = @. (Cell.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg*pt/Cell.Neg.L)) - pt*ν_neg*sinh(ν_neg) + Cell.Neg.L*(cosh(ν_neg)-cosh(ν_neg*(Cell.Neg.L-pt)/Cell.Neg.L)))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg)*ν_neg) #Lee. Eqn. 4.22
            zero_tf = @. -(pt^2)/(2*Cell.Const.CC_A*κ_eff_Neg*Cell.Neg.L)
            D[i]  =  @. (Cell.Neg.L*(σ_eff_Neg/κ_eff_Neg)*(1-cosh(ν_neg_∞*pt/Cell.Neg.L)) - pt*ν_neg_∞*sinh(ν_neg_∞) + Cell.Neg.L*(cosh(ν_neg_∞)-cosh(ν_neg_∞*(Cell.Neg.L-pt)/Cell.Neg.L)))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*sinh(ν_neg_∞)*ν_neg_∞) #Lee. Eqn. 4.22 @ ∞
            ϕ_tf[i,findall(s.==0)] .= zero_tf

         elseif pt <= Cell.Neg.L + Cell.Sep.L + eps()
            ϕ_tf[i,:] = @. (Cell.Neg.L - pt)/(Cell.Const.CC_A*κ_eff_Sep) + (Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg) #Lee. Eqn. 4.23
            zero_tf = @. (2*κ_eff_Neg*Cell.Neg.L-κ_eff_Sep*Cell.Neg.L-2*κ_eff_Neg*pt)/(2*Cell.Const.CC_A*κ_eff_Neg*κ_eff_Sep)
            D[i] = @. (Cell.Neg.L - pt)/(Cell.Const.CC_A*κ_eff_Sep) + (Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg_∞/2)-ν_neg_∞))/(Cell.Const.CC_A*(κ_eff_Neg+σ_eff_Neg)*ν_neg_∞) #Lee. Eqn. 4.23 @ ∞
            ϕ_tf[i,findall(s.==0)] .= zero_tf

         else
            ϕ_tf[i,:] = @. -Cell.Sep.L/(Cell.Const.CC_A*κ_eff_Sep) + Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg/2)-ν_neg)/(Cell.Const.CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg) + (Cell.Pos.L*(-σ_eff_Pos*cosh(ν_pos) + σ_eff_Pos*cosh((Cell.Const.Ltot-pt)*ν_pos/Cell.Pos.L) +  κ_eff_Pos*(cosh((pt-Cell.Neg.L-Cell.Sep.L)*ν_pos/Cell.Pos.L)-1)) - (pt-Cell.Neg.L-Cell.Sep.L)*κ_eff_Pos*sinh(ν_pos)*ν_pos)/(Cell.Const.CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos*sinh(ν_pos))
            zero_tf = @. - Cell.Neg.L/(2*Cell.Const.CC_A*κ_eff_Neg) - Cell.Sep.L/(Cell.Const.CC_A*κ_eff_Sep) - (Cell.Pos.L-(Cell.Const.Ltot-pt)^2/Cell.Pos.L)/(2*Cell.Const.CC_A*κ_eff_Pos)
            D[i] = @. -Cell.Sep.L/(Cell.Const.CC_A*κ_eff_Sep) + Cell.Neg.L*((1-σ_eff_Neg/κ_eff_Neg)*tanh(ν_neg_∞/2)-ν_neg_∞)/(Cell.Const.CC_A*(σ_eff_Neg+ κ_eff_Neg)*ν_neg_∞) + (Cell.Pos.L*(-σ_eff_Pos*cosh(ν_pos_∞) + σ_eff_Pos*cosh((Cell.Const.Ltot-pt)*ν_pos_∞/Cell.Pos.L) +  κ_eff_Pos*(cosh((pt-Cell.Neg.L-Cell.Sep.L)*ν_pos_∞/Cell.Pos.L)-1)) - (pt-Cell.Neg.L-Cell.Sep.L)*κ_eff_Pos*sinh(ν_pos_∞)*ν_pos_∞)/(Cell.Const.CC_A*κ_eff_Pos*(κ_eff_Pos+σ_eff_Pos)*ν_pos_∞*sinh(ν_pos_∞))
            ϕ_tf[i,findall(s.==0)] .= zero_tf
            end
      i+=1
 end
res0 .= zeros(length(z))
end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
# Test the JSON converter functionality

## Import packages

import PowerModels; const _PM = PowerModels
import PowerModelsACDC; const _PMACDC = PowerModelsACDC
import FlexPlan; const _FP = FlexPlan
import Cbc
optimizer = _FP.optimizer_with_attributes(Cbc.Optimizer, # Options: <https://github.com/jump-dev/Cbc.jl#options>
    "logLevel" => 0,
)

## Script parameters

file = "test/data/json_converter/case6_input_file_2018-2019.json"

## Parse the JSON file to easily check the input

#import JSON
#d = JSON.parsefile(file)

## Convert JSON file

mn_data = _FP.convert_JSON(file) # Conversion caveats and function parameters: see function documentation

## Instantiate model and solve network expansion problem

# Transmission network

setting = Dict("conv_losses_mp" => true)
result = _FP.simple_stoch_flex_tnep(mn_data, _PM.DCPPowerModel, optimizer; setting)
# Two-step alternative: exposes `pm`
#pm = _PM.instantiate_model(mn_data, _PM.DCPPowerModel, _FP.post_simple_stoch_flex_tnep; setting, ref_extensions=[_FP.ref_add_gen!, _PMACDC.add_ref_dcgrid!, _PMACDC.add_candidate_dcgrid!, _FP.ref_add_storage!, _FP.ref_add_ne_storage!, _FP.ref_add_flex_load!, _PM.ref_add_on_off_va_bounds!, _PM.ref_add_ne_branch!])
#result = _PM.optimize_model!(pm; optimizer=optimizer)

# Distribution network

#result = _FP.simple_stoch_flex_tnep(mn_data, _FP.BFARadPowerModel, optimizer)
# Two-step alternative: exposes `pm`
#pm = _PM.instantiate_model(mn_data, _FP.BFARadPowerModel, _FP.post_simple_stoch_flex_tnep; ref_extensions=[_FP.ref_add_gen!, _FP.ref_add_storage!, _FP.ref_add_ne_storage!, _FP.ref_add_flex_load!, _PM.ref_add_on_off_va_bounds!, _FP.ref_add_ne_branch_allbranches!, _FP.ref_add_frb_branch!, _FP.ref_add_oltc_branch!])
#result = _PM.optimize_model!(pm; optimizer=optimizer, solution_processors=[_PM.sol_data_model!])

println("Test completed")


text:
@testset "Deterministic" begin
    img = generate_surface(SynthSphere(50), 1.0, [0,0,1])
    p,q = synthetic_gradient(SynthSphere(40))
    Z,p,q = Deterministic(p = p, q = q, β=10.0)(img)
    @test p[76,76] ≈ 5.390522989786133e-7
    @test p[1,1] ≈ 0.0
    @test p[151,151] ≈ 0.0
    @test q[90,90] ≈ 0.40305773574048137
    @test q[101,36] ≈ 0.0
    @test q[22, 136] ≈ 0.0
    @test q[42,42] ≈ 0.0
end


text:
@block SimonDanisch ["2d"] begin

    @cell "Test heatmap + image overlap" [image, heatmap, transparency] begin
        heatmap(rand(32, 32))
        image!(map(x->RGBAf0(x,0.5, 0.5, 0.8), rand(32,32)))
    end

    @cell "poly and colormap" [poly, colormap, colorrang] begin
        # example by @Paulms from JuliaPlots/Makie.jl#310
        points = Point2f0[[0.0, 0.0], [0.1, 0.0], [0.1, 0.1], [0.0, 0.1]]
        colors = [0.0 ,0.0, 0.5, 0.0]
        scene = poly(points, color = colors, colorrange = (0.0,1.0))
        points = Point2f0[[0.1, 0.1], [0.2, 0.1], [0.2, 0.2], [0.1, 0.2]]
        colors = [0.5,0.5,1.0,0.3]
        poly!(scene, points, color = colors, colorrange = (0.0,1.0))
        scene
    end
    @cell "quiver" [quiver, arrows, vectorfield, gradient] begin
        x = range(-2, stop = 2, length = 21)
        arrows(x, x, rand(21, 21), rand(21, 21), arrowsize = 0.05)
    end
    @cell "Arrows on hemisphere" [arrows, quiver, mesh] begin
        using GeometryBasics
        s = Sphere(Point3f0(0), 0.9f0)
        scene = mesh(s, transparency=true, alpha=0.05)
        pos = decompose(Point3f0, s)
        dirs = decompose_normals(s)
        arrows!(scene, pos, dirs, arrowcolor=:red, arrowsize=0.1, linecolor=:red)
    end
    @cell "image" [image] begin
        vbox(
            image(AbstractPlotting.logo(), scale_plot = false),
            image(rand(100, 500), scale_plot = false),
        )
    end
    @cell "FEM polygon 2D" [fem, poly] begin
        coordinates = [
            0.0 0.0;
            0.5 0.0;
            1.0 0.0;
            0.0 0.5;
            0.5 0.5;
            1.0 0.5;
            0.0 1.0;
            0.5 1.0;
            1.0 1.0;
        ]
        connectivity = [
            1 2 5;
            1 4 5;
            2 3 6;
            2 5 6;
            4 5 8;
            4 7 8;
            5 6 9;
            5 8 9;
        ]
        color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
        poly(coordinates, connectivity, color = color, strokecolor = (:black, 0.6), strokewidth = 4)
    end
    @cell "FEM mesh 2D" [fem, mesh] begin
        coordinates = [
            0.0 0.0;
            0.5 0.0;
            1.0 0.0;
            0.0 0.5;
            0.5 0.5;
            1.0 0.5;
            0.0 1.0;
            0.5 1.0;
            1.0 1.0;
        ]
        connectivity = [
            1 2 5;
            1 4 5;
            2 3 6;
            2 5 6;
            4 5 8;
            4 7 8;
            5 6 9;
            5 8 9;
        ]
        color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
        scene = mesh(coordinates, connectivity, color = color, shading = false)
        wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 3)
    end
    @cell "colored triangle" [mesh, polygon] begin
        mesh(
            [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)], color = [:red, :green, :blue],
            shading = false
        )
    end
    @cell "colored triangle" [polygon] begin
        poly(
            [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)],
            color = [:red, :green, :blue],
            strokecolor = :black, strokewidth = 2
        )
    end
    @cell "Subscenes" [image, scatter, subscene] begin
        img = rand(RGBAf0, 100, 100)
        scene = image(img, show_axis = false)
        subscene = Scene(scene, IRect(100, 100, 300, 300))
        scatter!(subscene, rand(100) * 200, rand(100) * 200, markersize = 4)
        scene
    end

    @cell "scale_plot" [scale] begin
        t = range(0, stop=1, length=500) # time steps
        θ = (6π) .* t    # angles
        x =  # x coords of spiral
        y =  # y coords of spiral
        p1 = lines(t .* cos.(θ), t .* sin.(θ);
            color = t, colormap = :algae, linewidth = 8,
            scale_plot = false)
    end

    @cell "Polygons" [poly, polygon, linesegments] begin
        using GeometryBasics
        scene = Scene(resolution = (500, 500))
        points = decompose(Point2f0, Circle(Point2f0(50), 50f0))
        pol = poly!(scene, points, color = :gray, strokewidth = 10, strokecolor = :red)
        # Optimized forms
        poly!(scene, [Circle(Point2f0(50+300), 50f0)], color = :gray, strokewidth = 10, strokecolor = :red)
        poly!(scene, [Circle(Point2f0(50+i, 50+i), 10f0) for i = 1:100:400], color = :red)
        poly!(scene, [FRect2D(50+i, 50+i, 20, 20) for i = 1:100:400], strokewidth = 2, strokecolor = :green)
        linesegments!(scene,
            [Point2f0(50 + i, 50 + i) => Point2f0(i + 70, i + 70) for i = 1:100:400], linewidth = 8, color = :purple
        )
    end

    @cell "Contour Function" [contour] begin
        r = range(-10, stop = 10, length = 512)
        z = ((x, y)-> sin(x) + cos(y)).(r, r')
        contour(r, r, z, levels = 5, colormap = :viridis, linewidth = 3)
    end
    @cell "Hbox" [lines, scatter, hbox] begin
        t = range(-122277.9, stop=-14798.0, length=29542)
        x = -42 .- randn(length(t))
        sc1 = scatter(t, x, color=:black, markersize=sqrt(length(t)/20))
        sc2 = lines(t[1:end-1], diff(x), color = :blue)
        hbox(sc2, sc1)
    end
    @cell "Customize Axes" [lines, axis] begin
        x = LinRange(0,3pi,200); y = sin.(x)
        lin = lines(x, y, padding = (0.0, 0.0), axis = (
            names = (axisnames = ("", ""),),
            grid = (linewidth = (0, 0),),
        ))
    end
    @cell "contour" [contour] begin
        y = range(-0.997669, stop = 0.997669, length = 23)
        contour(range(-0.99, stop = 0.99, length = 23), y, rand(23, 23), levels = 10)
    end

    @cell "Text Annotation" [text, align, annotation] begin
        text(
            ". This is an annotation!",
            position = (300, 200),
            align = (:center,  :center),
            textsize = 60,
            font = "Blackchancery"
        )
    end

    @cell "Text rotation" [text, rotation] begin
        scene = Scene()
        pos = (500, 500)
        posis = Point2f0[]
        for r in range(0, stop = 2pi, length = 20)
            global pos, posis
            p = pos .+ (sin(r)*100.0, cos(r) * 100)
            push!(posis, p)
            t = text!(
                scene, "test",
                position = p,
                textsize = 50,
                rotation = 1.5pi - r,
                align = (:center, :center)
            )
        end
        scatter!(scene, posis, markersize = 10)
    end
    @cell "Standard deviation band" [band, lines, statistics] begin
        # Sample 100 Brownian motion path and plot the mean trajectory together
        # with a ±1σ band (visualizing uncertainty as marginal standard deviation).
        using Statistics
        n, m = 100, 101
        t = range(0, 1, length=m)
        X = cumsum(randn(n, m), dims = 2)
        X = X .- X[:, 1]
        μ = vec(mean(X, dims=1)) # mean
        lines(t, μ)              # plot mean line
        σ = vec(std(X, dims=1))  # stddev
        band!(t, μ + σ, μ - σ)   # plot stddev band
    end

    @cell "Streamplot animation" ["streamplot", "animation"] begin
        v(x::Point2{T}, t) where T = Point2{T}(one(T) * x[2] * t, 4 * x[1])
        sf = Node(Base.Fix2(v, 0e0))
        title_str = Node("t = 0.00")
        sp = streamplot(sf, -2..2, -2..2;
                        linewidth = 2, padding = (0, 0),
                        arrow_size = 0.09, colormap =:magma)
        sc = title(sp, title_str)
        record(sc, @replace_with_a_path(mp4), LinRange(0, 20, 5)) do i
            sf[] = Base.Fix2(v, i)
            title_str[] = "t = $(round(i; sigdigits = 2))"
        end
    end
end

@block AnshulSinghvi ["colors"] begin

    @cell "Line changing colour" [colors, lines, animation] begin
        scene = lines(rand(10); linewidth=10)

        record(scene, @replace_with_a_path(mp4), 1:255; framerate = 60) do i
               scene.plots[2][:color] = RGBf0(i/255, (255 - i)/255, 0) # animate scene
        end
    end

    @cell "streamplot" [arrows, lines, streamplot] begin
        struct FitzhughNagumo{T}
            ϵ::T
            s::T
            γ::T
            β::T
        end

        P = FitzhughNagumo(0.1, 0.0, 1.5, 0.8)
        f(x, P::FitzhughNagumo) = Point2f0(
            (x[1]-x[2]-x[1]^3+P.s)/P.ϵ,
            P.γ*x[1]-x[2] + P.β
        )
        f(x) = f(x, P)
        streamplot(f, -1.5..1.5, -1.5..1.5, colormap = :magma)
    end
end

@block AnshulSinghvi ["Transformations"] begin
    @cell "Transforming lines" [transformation, lines] begin
        N = 7 # number of colours in default palette
        sc = Scene()
        st = Stepper(sc, @replace_with_a_path)

        xs = 0:9        # data
        ys = zeros(10)

        for i in 1:N    # plot lines
            lines!(sc,
                xs, ys;
                color = AbstractPlotting.default_palettes.color[][i],
                limits = FRect((0, 0), (10, 10)),
                linewidth = 5
            ) # plot lines with colors
        end

        center!(sc)

        step!(st)

        for (i, rot) in enumerate(LinRange(0, π/2, N))
            rotate!(sc.plots[i+1], rot)
            arc!(sc,
                Point2f0(0),
                (8 - i),
                pi/2,
                (pi/2-rot);
                color = sc.plots[i+1].color,
                linewidth = 5,
                linestyle = :dash
            )
        end

        step!(st)
    end
end


text:
@testset "context" begin

@test ctx == CuCurrentContext()
@test ctx === CuCurrentContext()

@test_throws ErrorException deepcopy(ctx)

let ctx2 = CuContext(dev)
    @test ctx2 == CuCurrentContext()    # ctor implicitly pushes
    activate(ctx)
    @test ctx == CuCurrentContext()

    @test device(ctx2) == dev

    destroy!(ctx2)
end

instances = length(CUDAdrv.context_instances)
CuContext(dev) do ctx2
    @test length(CUDAdrv.context_instances) == instances+1
    @test ctx2 == CuCurrentContext()
    @test ctx != ctx2
end
@test length(CUDAdrv.context_instances) == instances
@test ctx == CuCurrentContext()

@test device(ctx) == dev
@test device() == dev
synchronize(ctx)
synchronize()

end

@testset "primary context" begin

pctx = CuPrimaryContext(dev)

@test !isactive(pctx)
unsafe_reset!(pctx)
@test !isactive(pctx)

@test flags(pctx) == CUDAdrv.SCHED_AUTO
setflags!(pctx, CUDAdrv.SCHED_BLOCKING_SYNC)
@test flags(pctx) == CUDAdrv.SCHED_BLOCKING_SYNC

CuContext(pctx) do ctx
    @test CUDAdrv.isvalid(ctx)
    @test isactive(pctx)
end
GC.gc()
@test !isactive(pctx)

CuContext(pctx) do ctx
    @test CUDAdrv.isvalid(ctx)
    @test isactive(pctx)

    unsafe_reset!(pctx)

    @test !isactive(pctx)
    @test !CUDAdrv.isvalid(ctx)
end

end

@testset "cache config" begin

cache_config!(CUDAdrv.FUNC_CACHE_PREFER_L1)
@test cache_config() == CUDAdrv.FUNC_CACHE_PREFER_L1

end

@testset "shmem config" begin

shmem_config!(CUDAdrv.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE)
@test shmem_config() == CUDAdrv.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE

end


text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)

struct DataLoader{D}
    data::D
    batchsize::Int
    nobs::Int
    partial::Bool
    imax::Int
    indices::Vector{Int}
    shuffle::Bool
end

"""
    DataLoader(data; batchsize=1, shuffle=false, partial=true)

An object that iterates over mini-batches of `data`, each mini-batch containing `batchsize` observations
(except possibly the last one).

Takes as input a single data tensor, or a tuple (or a named tuple) of tensors.
The last dimension in each tensor is considered to be the observation dimension.

If `shuffle=true`, shuffles the observations each time iterations are re-started.
If `partial=false`, drops the last mini-batch if it is smaller than the batchsize.

The original data is preserved in the `data` field of the DataLoader.

Usage example:

    Xtrain = rand(10, 100)
    train_loader = DataLoader(Xtrain, batchsize=2)
    # iterate over 50 mini-batches of size 2
    for x in train_loader
        @assert size(x) == (10, 2)
        ...
    end

    train_loader.data   # original dataset

    # similar, but yielding tuples
    train_loader = DataLoader((Xtrain,), batchsize=2)
    for (x,) in train_loader
        @assert size(x) == (10, 2)
        ...
    end

    Xtrain = rand(10, 100)
    Ytrain = rand(100)
    train_loader = DataLoader((Xtrain, Ytrain), batchsize=2, shuffle=true)
    for epoch in 1:100
        for (x, y) in train_loader
            @assert size(x) == (10, 2)
            @assert size(y) == (2,)
            ...
        end
    end

    # train for 10 epochs
    using IterTools: ncycle
    Flux.train!(loss, ps, ncycle(train_loader, 10), opt)

    # can use NamedTuple to name tensors
    train_loader = DataLoader((images=Xtrain, labels=Ytrain), batchsize=2, shuffle=true)
    for datum in train_loader
        @assert size(datum.images) == (10, 2)
        @assert size(datum.labels) == (2,)
    end
"""
function DataLoader(data; batchsize=1, shuffle=false, partial=true)
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))

    n = _nobs(data)
    if n < batchsize
        @warn "Number of observations less than batchsize, decreasing the batchsize to $n"
        batchsize = n
    end
    imax = partial ? n : n - batchsize + 1
    DataLoader(data, batchsize, n, partial, imax, [1:n;], shuffle)
end

@propagate_inbounds function Base.iterate(d::DataLoader, i=0)     # returns data in d.indices[i+1:i+batchsize]
    i >= d.imax && return nothing
    if d.shuffle && i == 0
        shuffle!(d.indices)
    end
    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i+1:nexti]
    batch = _getobs(d.data, ids)
    return (batch, nexti)
end

function Base.length(d::DataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

_nobs(data::AbstractArray) = size(data)[end]

function _nobs(data::Union{Tuple, NamedTuple})
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))
    n = _nobs(data[1])
    if !all(x -> _nobs(x) == n, Base.tail(data))
        throw(DimensionMismatch("All data should contain same number of observations"))
    end
    return n
end

_getobs(data::AbstractArray, i) = data[ntuple(i -> Colon(), Val(ndims(data) - 1))..., i]
_getobs(data::Union{Tuple, NamedTuple}, i) = map(Base.Fix2(_getobs, i), data)

Base.eltype(::DataLoader{D}) where D = D


text:
using TomoForward
using XfromProjections

# using MKLSparse # uncomment if you've installed MKLSparse, which will boost the performance

# generate synthetic data
img = zeros(128, 128, 128)
img[40:70, 40:60, 50:70] .= 1.0
H, W, nslice = size(img)

nangles = 90
detcount = Int(floor(size(img,1)*1.4))
proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])

isdefined_A = @isdefined A
if isdefined_A == false
    A = fp_op_parallel2d_strip(proj_geom, size(img, 1), size(img, 2))
end

p = zeros(nangles, nslice, detcount)
for i=1:nslice
    p[:,i,:] .= reshape(A * vec(img[:,:,i]), nangles, detcount)
end

# fbp slice by slice
q = filter_proj(p)
img3d = bp_slices(q, A, H, W)

# using PyPlot
# imshow(img3d[:,:,60])

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset verbose=true "test_seq_launch" begin

    @testset "execution_mode_1" begin
        println("\n\n\n")
        network = PSCOPF.Networks.Network()
        TS = PSCOPF.create_target_timepoints(DateTime("2015-01-01T11:00:00"))
        mode = PSCOPF.PSCOPF_MODE_1
        ECH = PSCOPF.generate_ech(network, TS, mode)

        sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

        exec_context = PSCOPF.PSCOPFContext(network, TS, mode)
        PSCOPF.run!(exec_context, sequence)
    end

    @testset "execution_mode_2" begin
        println("\n\n\n")
        network = PSCOPF.Networks.Network()
        TS = PSCOPF.create_target_timepoints(DateTime("2015-01-01T11:00:00"))
        mode = PSCOPF.PSCOPF_MODE_2
        ECH = PSCOPF.generate_ech(network, TS, mode)

        sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

        exec_context = PSCOPF.PSCOPFContext(network, TS, mode)
        PSCOPF.run!(exec_context, sequence)
    end

    @testset "execution_mode_3" begin
        println("\n\n\n")
        network = PSCOPF.Networks.Network()
        TS = PSCOPF.create_target_timepoints(DateTime("2015-01-01T11:00:00"))
        mode = PSCOPF.PSCOPF_MODE_3
        ECH = PSCOPF.generate_ech(network, TS, mode)

        sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

        exec_context = PSCOPF.PSCOPFContext(network, TS, mode)
        PSCOPF.run!(exec_context, sequence)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
## The Earth System Data Cube
![](http://earthsystemdatacube.net/wp-content/uploads/2015/07/EarthDataCube3.png "The DataCube")

Some info on the project...
"""
module ESDL
import Dates: Date
export Date
export ESDLdir, getAxis
export Cube, getCubeData,readcubedata,CubeMem,CubeAxis, TimeAxis, VariableAxis, LonAxis, LatAxis, SpatialPointAxis, saveCube, loadCube,
        RangeAxis, CategoricalAxis, MSCAxis, ScaleAxis, QuantileAxis, MethodAxis, cubeinfo, @caxis_str,splitdim,
        axVal2Index, mapCubeSimple, concatenateCubes, NetCDFCube, mergeAxes, caxes, subsetcube, CubeMask, renameaxis! #From Cube module
export registerDATFunction, mapCube, reduceCube, getAxis, InDims, OutDims, (..), Dataset, ESDLDataset,S3Cube,
        CubeTable, AsArray,AsAxisArray,AsDataFrame, cubefittable, TableAggregator, fittable #From DAT module
export cubeAnomalies, removeMSC, gapFillMSC, normalizeTS,DATfitOnline,
  sampleLandPoints, toPointAxis, getMSC, filterTSFFT, getNpY,
  getMedSC, extractLonLats,simpleAnomalies,spatialinterp,cubefromshape, exportcube #From Proc module
export rmCube # From CachedArrays
export @loadOrGenerate # from ESDL Tools
import Zarr
global const ESDLDefaults = (
  workdir = Ref("./"),
  recal   = Ref(false),
  compressor = Ref{Zarr.Compressor}(Zarr.NoCompressor()),
  chunksize  = Ref{Any}(:input),
  max_cache  = Ref(1e8),
  cubedir    = Ref(""),
)
global const workdir=ESDLDefaults.workdir
global const recal=ESDLDefaults.recal
function __init__()
  ESDLDefaults.workdir[]   = get(ENV,"ESDL_WORKDIR","./")
  ESDLDefaults.max_cache[] = parse(Float64,get(ENV,"ESDL_MAX_CACHE","100")) * 1e6
  ESDLDefaults.cubedir[]   = if isdir("/home/jovyan/work/datacube/ESDCv2.0.0/esdc-8d-0.25deg-184x90x90-2.0.0.zarr/")
    "/home/jovyan/work/datacube/ESDCv2.0.0/esdc-8d-0.25deg-184x90x90-2.0.0.zarr/"
  else
    get(ENV,"ESDL_CUBEDIR","")
  end
end
ESDLdir(x::String)=ESDLDefaults.workdir[]=x
recalculate(x::Bool)=ESDLDefaults.recal[]=x
recalculate()=ESDLDefaults.recal[]
ESDLdir()=ESDLDefaults.workdir[]
export ESDLdir

include("ESDLTools.jl")
include("Cubes/Cubes.jl")
#include("CubeAPI/CubeAPI.jl")
include("DAT/DAT.jl")
include("Proc/Proc.jl")

using .Cubes, .DAT, .Proc, .ESDLTools

#include("precompile.jl")
#_precompile_()

end # module


text:
#####
##### Evaluating the mistake rate of an AlphaZero agent
##### on the standard benchmark from Pascal Pons
#####

const SOLVER_ERROR =
  "To use this script, you should install Pascal Pons' Connect Four Solver.\n" *
  "See games/connect-four/solver/README.md for instructions."

const SESSION_DIR = "sessions/connect-four"
const SAVE_FILE = "pons-benchmark-results.json"
const PLOT_FILE = "pons-benchmark-results.png"
const ITC_STRIDE = 5
const DEBUG_MODE = false # Launches a quick run on a tiny dataset to help debug

# A benchmark to evaluate connect-four agents is available at:
#   http://blog.gamesolver.org/solving-connect-four/02-test-protocol/
# It features 6000 positions along with their expected scores.

# Let `ne` the number of elapsed moves. A game stage is either
# `begin` (ne <= 14), `middle` (14 < ne <= 28) or `end` (ne > 28)
const STAGES = [:beginning, :middle, :end]

# Let `nr` the number of remaining moves. A difficulty level is either
# `easy` (nr < 14), `medium` (14 <= nr < 28) or `hard` (nr > 28).
const DIFFICULTIES = [:easy, :medium, :hard]

# Return either `nothing` or a `(stage, level)` tuple
function parse_test_filename(name)
  try
    m = match(r"^Test_L(\d)_R(\d)$", name)
    L = parse(Int, m.captures[1])
    R = parse(Int, m.captures[2])
    return (STAGES[L], DIFFICULTIES[R])
  catch e
    rethrow(e)
    return nothing
  end
end

struct Bench
  stage :: Symbol
  difficulty :: Symbol
  entries :: Vector{Tuple{String, Int}}
end

function load_benchmarks(dir)
  @assert isdir(dir)
  benchmarks = Bench[]
  files = readdir(dir)
  for bf in files
    meta = parse_test_filename(bf)
    f = joinpath(dir, bf)
    if !isnothing(meta) && isfile(f)
      stage, difficulty = meta
      entries = []
      for L in readlines(f)
        try
          L = split(L)
          push!(entries, (L[1], parse(Int, L[2])))
        catch e
          rethrow(e)
        end
      end
    push!(benchmarks, Bench(stage, difficulty, entries))
    end
  end
  rank(b) = (
    findfirst(==(b.difficulty), DIFFICULTIES),
    findfirst(==(b.stage), STAGES))
  sort!(benchmarks, by=rank)
  return benchmarks
end

const BENCHMARKS_DIR = joinpath(@__DIR__, "..", "benchmark")
const BENCHMARKS = load_benchmarks(BENCHMARKS_DIR)

####
#### Testing code
####

using AlphaZero
include("../main.jl")
using .ConnectFour: Game, Solver, Training
using ProgressMeter
using Formatting

function state_of_string(str)
  g = Game()
  for c in str
    a = parse(Int, c)
    GI.play!(g, a)
  end
  return g
end

function optimal_on(solver, player, e)
  s = state_of_string(e[1])
  as, π = AlphaZero.think(player, s)
  api = argmax(π)
  qs = [Solver.qvalue(solver, s, a) for a in as]
  return sign(qs[api]) == maximum(sign, qs)
end

function test_player_on(player, benchmark, progress)
  solver = Solver.Player()
  nerr = 0
  entries = benchmark.entries
  DEBUG_MODE && (entries = Iterators.take(entries, 10))
  for e in entries
    optimal_on(solver, player, e) || (nerr += 1)
    next!(progress)
  end
  return nerr / length(entries)
end

function test_player(player)
  errs = Float64[]
  for bench in BENCHMARKS
    p = Progress(length(bench.entries))
    err = test_player_on(player, bench, p)
    push!(errs, err)
    err_str =
    println("($(bench.stage), $(bench.difficulty)): $(fmt(".2f", 100 * err))%")
  end
  return errs
end

function load_alphazero_player(dir)
  env = AlphaZero.UserInterface.load_env(
    Game, Training.Network{Game},
    AlphaZero.Log.Logger(devnull), dir, params=Training.params)
  mcts_params = MctsParams(env.params.arena.mcts,
    temperature=ConstSchedule(0.),
    dirichlet_noise_ϵ=0.)
  return MctsPlayer(env.bestnn, mcts_params)
end

####
#### Running all the tests
####

const Errors = Vector{Float64} # One error rate for each benchmark

# Type gathering all produced data.
struct Results
  minmax :: Errors # Minmax baseline
  alphazero :: Errors #
  alphazero_training :: Vector{Tuple{Int, Errors}}
end

function generate_data(session_dir)
  println("Testing the minmax baseline")
  player = MinMax.Player{Game}(depth=5, τ=0)
  minmax_errs = test_player(player)
  println("")

  println("Testing AlphaZero")
  player = load_alphazero_player(session_dir)
  az_errs = test_player(player)
  println("")

  println("Testing AlphaZero during training")
  az_training_errs = []
  itc = 0
  while true
    itdir = joinpath(session_dir, "iterations", "$itc")
    AlphaZero.UserInterface.valid_session_dir(itdir) || break
    DEBUG_MODE && itc > 10 && break
    println("Iteration $itc")
    player = load_alphazero_player(itdir)
    push!(az_training_errs, (itc, test_player(player)))
    itc += ITC_STRIDE
  end
  return Results(minmax_errs, az_errs, az_training_errs)
end

#####
##### Main
#####

if !isdir(Solver.DEFAULT_SOLVER_DIR)
  println(stderr, SOLVER_ERROR)
  exit()
end

using JSON2
using JSON3

JSON3.StructType(::Type{Results}) = JSON3.Struct()

# Regenerate data only if needed as it can take a while.
if isfile(SAVE_FILE)
  data = open(SAVE_FILE, "r") do io
    JSON3.read(io, Results)
  end
else
  data = generate_data(SESSION_DIR)
  open(SAVE_FILE, "w") do io
    JSON2.pretty(io, JSON3.write(data))
  end
end

import Plots

# i is the benchmark number
function sub_plot(benchs, results, n)
  bench = benchs[n]
  sname(s) = uppercase(string(s)[1]) * string(s)[2:end]
  title = "$(sname(bench.stage)) - $(sname(bench.difficulty))"
  azt = results.alphazero_training
  xs = [itc for (itc, _) in azt]
  ys = [100 * errs[n] for (itc, errs) in azt]
  y_minmax = 100 * results.minmax[n]
  #ymax = max(5, maximum(ys))
  ymax = max(maximum(ys), y_minmax)
  plt = Plots.plot(
    xs, ys, ylims=(0, ymax),
    # xlabel="Training Iteration",
    ylabel="Error Rate (in %)",
    legend=nothing,
    title=title,
    titlefontsize=10,
    xguidefontsize=8,
    yguidefontsize=8)
  Plots.hline!(plt, [y_minmax])
  return plt
end

function plot_results(results)
  plts = [sub_plot(BENCHMARKS, results, n) for n in eachindex(BENCHMARKS)]
  return Plots.plot(plts..., layout=(3, 2))
end

plot = plot_results(data)
Plots.savefig(plot, PLOT_FILE)


text:
function example()
    println("Hello ")
    @goto world
    println("Never printed")
    @label world
    println("world")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Before running this, please make sure to activate and instantiate the
# environment with [this `Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-crabs-xgb/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-crabs-xgb/Manifest.toml).
# For instance, copy these files to a folder 'EX-crabs-xgb', `cd` to it and
#
# ```julia
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# This example is inspired from [this post](https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/) showing how to use XGBoost.
#
# ## First steps
#
# Again, the crabs dataset is so common that there is a  simple load function for it:

using MLJ
using StatsBase
using Random
using PyPlot
using CategoricalArrays
using PrettyPrinting
import DataFrames

X, y = @load_crabs
X = DataFrames.DataFrame(X)
@show size(X)
@show y[1:3]
first(X, 3) |> pretty

# It's a classification problem with the following classes:

levels(y) |> pprint

# Note that the dataset is currently sorted by target, let's shuffle it to avoid the obvious issues this may cause

Random.seed!(523)
perm = randperm(length(y))
X = X[perm,:]
y = y[perm];

# It's not a very big dataset so we will likely overfit it badly using something as sophisticated as XGBoost but it will do for a demonstration.

train, test = partition(collect(eachindex(y)), 0.70, shuffle=true, rng=52)
XGBC = @load XGBoostClassifier
xgb_model = XGBC()

# Let's check whether the training and  is balanced, `StatsBase.countmap` is useful for that:

countmap(y[train]) |> pprint

# which is pretty balanced. You could check the same on the test set and full set and the same comment would still hold.
#
# ## XGBoost machine
#
# Wrap a machine around an XGBoost model (XGB) and the data:

xgb  = XGBC()
xgbm = machine(xgb, X, y)

# We will tune it varying the number of rounds used and generate a learning curve

r = range(xgb, :num_round, lower=50, upper=500)
curve = learning_curve(xgbm, range=r, resolution=50,
                        measure=L1HingeLoss())

# Let's have a look

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("HingeLoss", fontsize=14)
xticks([10, 100, 200, 500], fontsize=12)

# \figalt{Cross entropy vs Num Round}{EX-crabs-xgb-curve1.svg}
#
# So, in short, using more rounds helps. Let's arbitrarily fix it to 200.

xgb.num_round = 200;

# ### More tuning (1)
#
# Let's now tune the maximum depth of each tree and the minimum child weight in the boosting.

r1 = range(xgb, :max_depth, lower=3, upper=10)
r2 = range(xgb, :min_child_weight, lower=0, upper=5)

tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=11), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

# Great, as always we can investigate the tuning by using `report` and can, for instance, plot a heatmap of the measurements:

r = report(mtm)

res = r.plotting

md = res.parameter_values[:,1]
mcw = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, res.measurements)

xlabel("Maximum tree depth", fontsize=14)
ylabel("Minimum child weight", fontsize=14)
xticks(3:2:10, fontsize=12)
yticks(fontsize=12)

# \figalt{Hyperparameter heatmap}{EX-crabs-xgb-heatmap.svg}
#
# Let's extract the optimal model and inspect its parameters:

xgb = fitted_params(mtm).best_model
@show xgb.max_depth
@show xgb.min_child_weight

# ### More tuning (2)
#
# Let's examine the effect of `gamma`:

xgbm = machine(xgb, X, y)
r = range(xgb, :gamma, lower=0, upper=10)
curve = learning_curve!(xgbm, range=r, resolution=30,
                        measure=cross_entropy);

# it looks like the `gamma` parameter substantially affects model performance:

@show round(minimum(curve.measurements), sigdigits=3)
@show round(maximum(curve.measurements), sigdigits=3)

# ### More tuning (3)
#
# Let's examine the effect of `subsample` and `colsample_bytree`:

r1 = range(xgb, :subsample, lower=0.6, upper=1.0)
r2 = range(xgb, :colsample_bytree, lower=0.6, upper=1.0)
tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=234), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

# and the usual procedure to visualise it:

r = report(mtm)

res = r.plotting

ss = res.parameter_values[:,1]
cbt = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(ss, cbt, res.measurements)

xlabel("Sub sample", fontsize=14)
ylabel("Col sample by tree", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

# \figalt{Hyperparameter heatmap}{EX-crabs-xgb-heatmap2.svg}
#
# Let's retrieve the best models:

xgb = fitted_params(mtm).best_model
@show xgb.subsample
@show xgb.colsample_bytree

# We could continue with more fine tuning but given how small the dataset is, it doesn't make much sense.
# How does it fare on the test set?

ŷ = predict_mode(mtm, rows=test)
round(accuracy(ŷ, y[test]), sigdigits=3)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# deps.jl is created at the end of a successful build, so rm
# to ensure that failed builds are missing this file.
if isfile("deps.jl")
  rm("deps.jl")
end
include("build_petscs.jl")


text:
#
# start - covered by test cases
#

bv(v::Float64, bins::Int64, min=-1.0, max=1.0) = minimum([int64(floor(bins * (v - min) / (max - min))), bins-1])

bin_vector(vec::Vector{Float64}, min::Float64, max::Float64, bins::Int64) = map(v->bv(v, bins, min, max), vec)

function unbin_value(v::Int64, bins::Int64, min=-1.0, max=1.0; mode="centre")
  known_mode = (mode == "centre" ||
                mode == "lower"  ||
                mode == "upper")
  @assert known_mode "Mode may be any of the following: [\"centre\", \"lower\", \"upper\"]"

  delta = (max - min) / float64(bins)
  u = (v - 1) * delta + min

  if     mode == "centre"
    return u + 0.5 * delta
  elseif mode == "upper"
    return u + delta
  end
  return u
end

function bin_matrix(m::Matrix{Float64}, min::Float64, max::Float64, bins::Int64)
  r = zeros(size(m))
  for i=1:size(m)[1]
    r[i,:] = bin_vector(squeeze(m[i,:],1), min, max, bins)
  end
  convert(Matrix{Int64}, r)
end

function unbin_matrix(m::Matrix{Float64}, min::Float64, max::Float64, bins::Int64; mode="centre")
  r = zeros(size(m))
  for i=1:size(m)[1]
    r[i,:] = unbin_vector(squeeze(m[i,:],1), min, max, bins, mode=mode)
  end
  convert(Matrix{Int64}, r)
end

unbin_vector(vec::Vector{Float64}, min::Float64, max::Float64, bins::Int64; mode="centre") = map(v->unbin_value(v, bins, min, max, mode=mode), vec)

i2b(v::Int64, n::Int64) = [(((1 << (n-i)) & v)>0)?1.0:0.0 for i=1:n]
b2i(v::Vector{Float64}) = int64(foldl(+, [(v[i]>0)?(1<<(length(v)-i)):0 for i=1:length(v)]))

iv2b(v::Vector{Int64},   n::Int64) = foldl(vcat, [i2b(u ,n) for u in v])
b2iv(v::Vector{Float64}, n::Int64) = foldl(vcat, [b2i(v[i:i+(n-1)]) for i=1:n:length(v)])

bv2dv(v::Vector{Float64}, bins::Int64; mode="centre") = map(x->unbin_value(int64(x+1), bins, -1.0, 1.0, mode=mode), b2iv(v, int(ceil(log2(bins)))))

function binarise_matrix(A::Matrix{Float64}, bins::Int64)
  N = int(ceil(log2(bins)))
  B = zeros(size(A)[1], size(A)[2]* N)
  C = bin_matrix(A, -1.0, 1.0, bins)
  for row_index = 1:size(A)[1]
    B[row_index,:] = iv2b(squeeze(C[row_index,:],1), N)
  end
  B
end

# .- because bin_vector in [1, bins]
binarise_vector(v::Vector{Float64}, bins::Int64) = iv2b(bin_vector(v, -1.0, 1.0, bins), int(ceil(log2(bins))))

function unbinarise_matrix(A::Matrix{Float64}, bins::Int64; mode="centre")
  N = int(ceil(log2(bins)))
  w = size(A)[1]
  v = int(size(A)[2] / N)
  B = zeros(w, v)
  for row_index = 1:w
    B[row_index,:] = b2iv(squeeze(A[row_index,:],1), N)
  end
  map(x->unbin_value(int64(x+1), bins, -1.0, 1.0, mode=mode), B)
end

#
# end   - covered by test cases
#


text:
const ColKey = Union{Symbol,AbstractString}

"""
    coerce(A, specs...; tight=false, verbosity=1)

Given a table `A`, return a copy of `A`, ensuring that the element
scitypes of the columns match the new specification, `specs`. There
are three valid specifiations:

(i) one or more `column_name=>Scitype` pairs:

    coerce(X, col1=>Sciyype1, col2=>Scitype2, ... ; verbosity=1)

(ii) one or more `OldScitype=>NewScitype` pairs (`OldScitype` covering
both the `OldScitype` and `Union{Missing,OldScitype}` cases):

    coerce(X, OldScitype1=>NewSciyype1, OldScitype2=>NewScitype2, ... ; verbosity=1)

(iii) a dictionary of scientific types keyed on column names:

    coerce(X, d::AbstractDict{<:ColKey, <:Type}; verbosity=1)

where `ColKey = Union{Symbol,AbstractString}`.

### Examples

Specifiying  `column_name=>Scitype` pairs:

```
using CategoricalArrays, DataFrames, Tables
X = DataFrame(name=["Siri", "Robo", "Alexa", "Cortana"],
              height=[152, missing, 148, 163],
              rating=[1, 5, 2, 1])
Xc = coerce(X, :name=>Multiclass, :height=>Continuous, :rating=>OrderedFactor)
schema(Xc).scitypes # (Multiclass, Continuous, OrderedFactor)
```

Specifying `OldScitype=>NewScitype` pairs:

```
X  = (x = [1, 2, 3],
      y = rand(3),
      z = [10, 20, 30])
Xc = coerce(X, Count=>Continuous)
schema(Xfixed).scitypes # (Continuous, Continuous, Continuous)
```
"""
coerce(X, a...; kw...) = coerce(Val(ST.trait(X)), X, a...; kw...)

# Non tabular data is not supported
coerce(::Val{:other}, X, a...; kw...) =
    throw(CoercionError("`coerce` is undefined for non-tabular data."))

_bad_dictionary() = throw(ArgumentError(
    "A dictionary specifying a scitype conversion "*
    "must have type `AbstractDict{<:ColKey, <:Type}`. It's keys must "*
    "be column names and its values be scientific types. "*
    "E.g., `Dict(:cats=>Continuous, :dogs=>Textual`. "))
coerce(::Val{:table}, X, types_dict::AbstractDict; kw...) =
    _bad_dictionary()

_bad_specs() =
        throw(ArgumentError(
        "Invalid `specs` in `coerce(X, specs...;  kwargs...)`. "*
        "Valid `specs` are: (i) one or more pairs of "*
        "the form `column_name=>Scitype`; (ii) one or more pairs "*
        "of the from `OldScitype=>NewScitype`; or (iii) a "*
        "dictionary of scientific "*
        "types keyed on column names. "))
coerce(::Val{:table}, X, specs...; kw...) = _bad_specs()

function coerce(::Val{:table},
                X,
                types_dict::AbstractDict{<:ColKey, <:Type};
                kw...)
    isempty(types_dict) && return X
    names  = schema(X).names
    X_ct   = Tables.columntable(X)
    ct_new = (_coerce_col(X_ct, col, types_dict; kw...) for col in names)
    return Tables.materializer(X)(NamedTuple{names}(ct_new))
end

# -------------------------------------------------------------
# utilities for coerce

struct CoercionError <: Exception
    m::String
end

function _coerce_col(Xcol,
                     name,
                     types_dict::AbstractDict{Symbol, <:Type};
                     kw...)
    y = Tables.getcolumn(Xcol, name)
    if haskey(types_dict, name)
        coerce_type = types_dict[name]
        return coerce(y, coerce_type; kw...)
    end
    return y
end

# -------------------------------------------------------------
# alternative ways to do coercion, both for coerce and coerce!

# The following extends the two methods so that a mixture of
# symbol=>type and type=>type pairs can be specified in place of a
# dictionary:

feature_scitype_pairs(p::Pair{<:ColKey,<:Type}, X) = [Symbol(first(p)) => last(p), ]
function feature_scitype_pairs(p::Pair{<:Type,<:Type}, X)
    from_scitype = first(p)
    to_scitype = last(p)
    sch = schema(X)
    ret = Pair{Symbol,Type}[]
    for j in eachindex(sch.names)
        if sch.scitypes[j] <: Union{Missing,from_scitype}
            push!(ret, Pair(sch.names[j], to_scitype))
        end
    end
    return ret
end

for c in (:coerce, :coerce!)
    ex = quote
        function $c(::Val{:table},
                    X,
                    mixed_pairs::Pair{<:Union{<:ColKey,<:Type},<:Type}...;
                    kw...)
            components = map(p -> feature_scitype_pairs(p, X), mixed_pairs)
            pairs = vcat(components...)

            # must construct dictionary by hand to check no conflicts:
            scitype_given_feature = Dict{Symbol,Type}()
            for p in pairs
                feature = first(p)
                if haskey(scitype_given_feature, feature)
                    throw(ArgumentError("`coerce` argments cannot be "*
                                        "resolved to determined a "*
                                        "*unique* scitype for each "*
                                        "feature. "))
                else
                    scitype_given_feature[feature] = last(p)
                end
            end

            return $c(X, scitype_given_feature; kw...)
        end
    end
    eval(ex)
end

# -------------------------------------------------------------
# In place coercion

"""
coerce!(X, ...)

Same as [`ScientificTypes.coerce`](@ref) except it does the modification in
place provided `X` supports in-place modification (at the moment, only the
DataFrame! does). An error is thrown otherwise. The arguments are the same as
`coerce`.

"""
coerce!(X, a...;  kw...) = begin
    coerce!(Val(ST.trait(X)), X, a...; kw...)
end

coerce!(::Val{:other}, X, a...; kw...) =
    throw(CoercionError("`coerce!` is undefined for non-tabular data."))

coerce!(::Val{:table}, X, types_dict::AbstractDict; kw...) =
    _bad_dictionary()

coerce!(::Val{:table}, X, specs...; kw...) = _bad_specs()

function coerce!(::Val{:table},
                 X,
                 types_dict::AbstractDict{<:ColKey, <:Type};
                 kw...)
    # DataFrame --> coerce_df!
    if is_type(X, :DataFrames, :DataFrame)
        return coerce_df!(X, types_dict; kw...)
    end
    # Everything else
    throw(ArgumentError("In place coercion not supported for $(typeof(X))." *
                        "Try `coerce` instead."))
end

# -------------------------------------------------------------
# utilities for coerce!

"""
    coerce_df!(df, pairs...; kw...)

In place coercion for a dataframe.(Unexported method)
"""
function coerce_df!(df, tdict::AbstractDict{<:ColKey, <:Type}; kw...)
    names = schema(df).names
    for name in names
        name in keys(tdict) || continue
        coerce_type = tdict[name]
        df[!, name] = coerce(df[!, name], coerce_type; kw...)
    end
    return df
end

"""
    is_type(X, spkg, stype)

Check that an object `X` is of a given type that may be defined in a package
that is not loaded in the current environment.
As an example say `DataFrames` is not loaded in the current environment, a
function from some package could still return a DataFrame in which case it
can be checked with

```
is_type(X, :DataFrames, :DataFrame)
```
"""
function is_type(X, spkg::Symbol, stype::Symbol)
    # If the package is loaded, then it will just be `stype`
    # otherwise it will be `spkg.stype`
    rx = Regex("^($spkg\\.)?$stype")
    return ifelse(match(rx, "$(typeof(X))") === nothing, false, true)
end


text:
module Provider

using ZipFile
using Shapefile
using DBaseReader

export polygons, records, projection, names, describenames

abstract AbstractProvider
abstract ShapefileProvider <: AbstractProvider

###############################################################################
#### All Provider's must implement these three functions
# to return the appropriate data
"""
`polygons(provider::AbstractProvider)`

Returns an array of polygons (shapes) associated with a provider
"""
function polygons(provider::AbstractProvider)
  return []
end

"""
`records(provider::AbstractProvider)`

Returns a `Dict` of records (without polygons) associated with a provider
"""
function records(provider::AbstractProvider)
  return Dict()
end

"""
`projection(provider::AbstractProvider)`

Returns the wkid of the Provider's crs/srs/projection as an `AbstractString`.
"""
function projection(provider::AbstractProvider)
  return ""
end

# And may optionally implement these two functions for returning record column
# names and descriptions
"""
`names(provider::AbstractProvider)`

Returns an array containing the names of the columns contained in the
Provider's records
"""
function names(provider::AbstractProvider)
  return []
end

"""
`describenames(provider::AbstractProvider)`

Returns a `Dict` associating a Provider's columns names with a description
of the column.
"""
function describenames(provider::AbstractProvider)
  return Dict()
end
###############################################################################

include("shapefile_provider.jl")
include("us_tiger_provider.jl")

end # module


text:
struct KaHyParPartitioner
    configuration::Union{Nothing, Symbol, String}
    imbalance::Union{Nothing, Float64}
    weight::Any
    verbose::Bool
end

@enum KaHyParConfiguration begin
    kahypar_configuration_communication_volume = 1
    kahypar_configuration_hyperedge_cut = 2
end

function KaHyParPartitioner(;
    configuration::Union{Nothing, KaHyParConfiguration, String}=nothing,
    imbalance::Union{Nothing, Float64} = nothing,
    weight::Any = VertexCount(),
    verbose = false
)
    if isa(configuration, KaHyParConfiguration)
        if configuration == kahypar_configuration_communication_volume
            configuration = :connectivity
        elseif configuration == kahypar_configuration_hyperedge_cut
            configuration = :edgecut
        end
    end

    return KaHyParPartitioner(
        configuration,
        imbalance,
        weight,
        verbose,
    )
end

function partition_stripe(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, K, method::KaHyParPartitioner; kwargs...) where {Tv, Ti}
    return partition_stripe(SparseMatrixCSC(A), K, method)
    #TODO gotta add the column permutation stuff
end

function partition_stripe(A::SparseMatrixCSC, K, method::KaHyParPartitioner; kwargs...)
    m, n = size(A)
    if method.imbalance === nothing
        imbalance = ()
    else
        imbalance = (:imbalance=>method.imbalance,)
    end
    if method.verbose
        asg = KaHyPar.partition(
            KaHyPar.HyperGraph(
                adjointpattern(A),
                compute_weight(A, method.weight),
                ones(Int64, m)
            ),
            K;
            configuration=method.configuration,
            imbalance...
        )
    else
        asg = @suppress KaHyPar.partition(
            KaHyPar.HyperGraph(
                adjointpattern(A),
                compute_weight(A, method.weight),
                ones(Int64, m)
            ),
            K;
            configuration=method.configuration,
            imbalance...
        )
    end
    return MapPartition(K, asg .+ true)
end


text:
# Statistics-related math functions

# scalar functions

xlogx(x::FloatingPoint) = x > zero(x) ? x * log(x) : zero(x)
xlogx(x::Real) = xlogx(float(x))

xlogy{T<:FloatingPoint}(x::T, y::T) = x > zero(T) ? x * log(y) : zero(x)
xlogy{T<:Real}(x::T, y::T) = xlogy(float(x), float(y))
xlogy(x::Real, y::Real) = xlogy(promote(x, y)...)

logistic(x::FloatingPoint) = one(x) / (one(x) + exp(-x))
logistic(x::Real) = logistic(float(x))

logit(x::FloatingPoint) = log(x / (one(x) - x))
logit(x::Real) = logit(float(x))

softplus(x::FloatingPoint) = x <= 0 ? log1p(exp(x)) : x + log1p(exp(-x))
softplus(x::Real) = softplus(float(x))

invsoftplus(x::FloatingPoint) = log(exp(x) - one(x))
invsoftplus(x::Real) = invsoftplus(float(x))

@vectorize_1arg Real xlogx
@vectorize_2arg Real xlogy
@vectorize_1arg Real logistic
@vectorize_1arg Real logit
@vectorize_1arg Real softplus
@vectorize_1arg Real invsoftplus

## logsumexp

logsumexp{T<:FloatingPoint}(x::T, y::T) = x > y ? x + log1p(exp(y - x)) : y + log1p(exp(x - y))
logsumexp{T<:Real}(x::T, y::T) = logsumexp(float(x), float(y))
logsumexp(x::Real, y::Real) = logsumexp(promote(x, y)...)

function logsumexp{T<:Real}(x::AbstractArray{T})
    isempty(x) && return -Inf
    u = maximum(x)
    s = 0.
    for i = 1:length(x)
    	@inbounds s += exp(x[i] - u)
    end
    log(s) + u
end

## softmax

function softmax!{T<:Real}(r::AbstractArray, x::AbstractArray{T})
	n = length(x)
	length(r) == n || throw(DimensionMismatch("Inconsistent array lengths."))
	u = maximum(x)
	s = 0.
	for i = 1:n
		@inbounds s += (r[i] = exp(x[i] - u))
	end
	for i = 1:n
		r[i] /= s
	end
	r
end

softmax!{T<:FloatingPoint}(x::AbstractArray{T}) = softmax!(x, x)
softmax{T<:Real}(x::AbstractArray{T}) = softmax!(Array(Float64, size(x)), x)


text:
cd(@__DIR__); allNGWPs = false; include("setups/toronto.jl")
gr(dpi = 200)

## pedestrian signal 16 most important PC-NGWP vectors
PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal)
fp = load("../datasets/new_toronto.jld", "fp")
G_Sig.f = reshape(fp, (N, 1))
dmatrix_PC = ngwp_analysis(G_Sig, PC_NGWP)
dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GP_dual)
important_idx = sortperm(dvec_pc_ngwp[:].^2; rev = true)
println("================fp-PC-NGWP-top-basis-vectors=================")
for i in 1:16
    dr, dc = BS_pc_ngwp.levlist[important_idx[i]]
    w = PC_NGWP[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual, dr, dc)
    print("(j, k, l) = ($(j), $(k), $(l))  ")
    if j == jmax
        print("φ_{$(GP_dual.ind[dr]-1)}")
    end
    println()
    sgn = (maximum(w) > -minimum(w)) * 2 - 1
    gplot(A, X; width=1)
    scatter_gplot!(X; marker = sgn .* w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    savefig(plt,
        "../figs/Toronto_fp_DAG_PC_NGWP_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
"""
    PermuteDims <: Augmentor.Operation

Description
--------------

Permute the dimensions of the given array with the predefined
permutation `perm`. This operation is particularly useful if the
order of the dimensions needs to be different than the default
"julian" layout (described below).

Augmentor expects the given images to be in vertical-major layout
for which the colors are encoded in the element type itself. Many
deep learning frameworks however require their input in a
different order. For example it is not untypical that separate
color channels are expected to be encoded in the third dimension.

Usage
--------------

    PermuteDims(perm)

    PermuteDims(perm...)

Arguments
--------------

- **`perm`** : The concrete dimension permutation that should be
    used. Has to be specified as a `Vararg{Int}` or as a `NTuple`
    of `Int`. The length of `perm` has to match the number of
    dimensions of the expected input image to that operation.

See also
--------------

[`SplitChannels`](@ref), [`CombineChannels`](@ref), [`augment`](@ref)

Examples
--------------

```julia-repl
julia> using Augmentor, Colors

julia> A = rand(10, 5, 3) # width=10, height=5, and 3 color channels
10×5×3 Array{Float64,3}:
[...]

julia> img = augment(A, PermuteDims(3,2,1) |> CombineChannels(RGB))
5×10 Array{RGB{Float64},2}:
[...]

julia> img2 = testpattern()
300×400 Array{RGBA{N0f8},2}:
[...]

julia> B = augment(img2, SplitChannels() |> PermuteDims(3,2,1))
400×300×4 Array{N0f8,3}:
[...]
```
"""
struct PermuteDims{N,perm,iperm} <: Operation end
PermuteDims() = throw(MethodError(PermuteDims, ()))
PermuteDims(perm::Tuple{}) = throw(MethodError(PermuteDims, (perm,)))
PermuteDims(perm::NTuple{N,Int}) where {N} = PermuteDims{N,perm,invperm(perm)}()
PermuteDims(perm::Vararg{Int,N}) where {N} = PermuteDims{N,perm,invperm(perm)}()

@inline supports_eager(::Type{<:PermuteDims}) = true
@inline supports_lazy(::Type{<:PermuteDims}) = true

function applyeager(op::PermuteDims{N,perm}, img::AbstractArray{T,N}, param) where {T,N,perm}
    plain_array(permutedims(img, perm))
end

function applylazy(op::PermuteDims{N,perm,iperm}, img::AbstractArray{T,N}, param) where {T,N,perm,iperm}
    PermutedDimsArray{T,N,perm,iperm,typeof(img)}(img)
end

function showconstruction(io::IO, op::PermuteDims{N,perm}) where {N,perm}
    print(io, nameof(typeof(op)), '(', join(map(string, perm),", "), ')')
end

function Base.show(io::IO, op::PermuteDims{N,perm}) where {N,perm}
    if get(io, :compact, false)
        print(io, "Permute dimension order to ", perm)
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end

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

"""
    Reshape <: Augmentor.Operation

Description
--------------

Reinterpret the shape of the given array of numbers or colorants.
This is useful for example to create singleton-dimensions that
deep learning frameworks may need for colorless images, or for
converting an image array to a feature vector (and vice versa).

Usage
--------------

    Reshape(dims)

    Reshape(dims...)

Arguments
--------------

- **`dims`** : The new sizes for each dimension of the output
    image. Has to be specified as a `Vararg{Int}` or as a
    `NTuple` of `Int`.

See also
--------------

[`CombineChannels`](@ref), [`augment`](@ref)

Examples
--------------

```julia-repl
julia> using Augmentor, Colors

julia> A = rand(10,10)
10×10 Array{Float64,2}:
[...]

julia> augment(A, Reshape(10,10,1)) # add trailing singleton dimension
10×10×1 Array{Float64,3}:
[...]
```
"""
struct Reshape{N} <: Operation
    dims::NTuple{N,Int}
end
Reshape() = throw(MethodError(Reshape, ()))
Reshape(dims::Tuple{}) = throw(MethodError(Reshape, (dims,)))
Reshape(dims::Int...) = Reshape(dims)

@inline supports_eager(::Type{<:Reshape}) = false
@inline supports_lazy(::Type{<:Reshape}) = true

applylazy(op::Reshape, img::AbstractArray, param) = reshape(plain_axes(img), op.dims)

function showconstruction(io::IO, op::Reshape)
    print(io, nameof(typeof(op)), '(', join(map(string, op.dims),", "), ')')
end

function Base.show(io::IO, op::Reshape{N}) where N
    if get(io, :compact, false)
        if N == 1
            print(io, "Reshape array to ", first(op.dims), "-element vector")
        else
            print(io, "Reshape array to ", join(op.dims,"×"))
        end
    else
        print(io, "Augmentor.")
        print(io, nameof(typeof(op)), '{', N, '}', '(', join(map(string, op.dims),", "), ')')
    end
end


text:
# This file contains functions for the Cost Efficiency Hölder DEA model
"""
    CostHolderDEAModel
An data structure representing a cost Hölder DEA model.
"""
struct CostHolderDEAModel <: AbstractCostDEAModel
    n::Int64
    m::Int64
    s::Int64
    l::Union{Int64,Float64}
    isweighted::Bool
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    rts::Symbol
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deacostholder(X, Y, W; l)
Compute cost efficiency using data envelopment analysis for
inputs `X`, outputs `Y` and price of inputs `W`.

# Hölder norm `l` specification
- `1`.
- `2`.
- `Inf`.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];

julia> Y = [1; 1; 1; 1; 1; 1; 1; 1];

julia> W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> deacostholder(X, Y, W, l = 1)
Cost Hölder L1 DEA Model
DMUs = 8; Inputs = 2; Outputs = 1
Orientation = Input; Returns to Scale = VRS
──────────────────────────────
   Cost  Technical  Allocative
──────────────────────────────
1   0.0        0.0         0.0
2   1.0        0.0         1.0
3   1.0        0.0         1.0
4   3.0        2.0         1.0
5   6.0        4.0         2.0
6   3.0        0.0         3.0
7   3.0        1.0         2.0
8   5.6        0.6         5.0
──────────────────────────────
```
"""
function deacostholder(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector};
    l::Union{Int64,Float64}, weight::Bool = false,
    rts::Symbol = :VRS, monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::CostHolderDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end

    if l != 1 && l != 2 && l != Inf
        throw(ArgumentError("l must by :1, :2, or :Inf"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get minimum cost targets and lambdas
    n = nx

    Xtarget, clambdaeff = deamincost(X, Y, W, rts = rts, dispos = :Strong, optimizer = optimizer)
    Ytarget = Y[:,:]

    # Cost, technical and allocative efficiency
    cefficiency  = vec(sum(W .* X, dims = 2) .- sum(W .* Xtarget, dims = 2))

    if weight
        if l == 1
            normalization = maximum(W .* X, dims = 2)
        elseif l == 2
            normalization = sqrt.(sum((W .* X).^2, dims = 2))
        elseif l == Inf
            normalization = sum(W .* X, dims = 2)
        end
    else
        if l == 1
            normalization = maximum(W, dims = 2)
        elseif l == 2
            normalization = sqrt.(sum((W ).^2, dims = 2))
        elseif l == Inf
            normalization = sum(W, dims = 2)
        end
    end
    normalization = vec(normalization)

    techefficiency = efficiency(deaholder(X, Y, l = l, weight = weight, orient = :Input, rts = rts, slack = false, optimizer = optimizer))

    if monetary
        techefficiency = techefficiency .* normalization
    else
        cefficiency = cefficiency ./ normalization
    end

    allocefficiency = cefficiency - techefficiency

    return CostHolderDEAModel(n, m, s, l, weight, monetary, names, rts, cefficiency, clambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)
end

function Base.show(io::IO, x::CostHolderDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    l = x.l
    isweighted = x.isweighted
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Cost Hölder L", string(x.l), " DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Input")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if isweighted
            print(io, "Weighted (weakly) Hölder distance function \n")
        end
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Cost", "Technical", "Allocative"], dmunames))
    end

end


text:
struct MTGRU_Single end
struct MTGRU_Double end

"""
    MTGRU_variational(enc_tform, init_enc, init_post, mt_enc, mt_post,
        x0decoder, mtbias_deocder, gen_model, d_y, mtbias_only)

A sequence-to-sequence model 𝒴→𝒴 which imbibes `T_steps` steps of the sequence,
and predicts `T_steps` (including the original `T_steps`), but unlike the
`Seq2SeqModel`, the `MTSeqModel_Pred` can *modulate* its prediction based on
hierarchical 'multi-task' variables. The `Pred` suffix refers to the function
of the encoders. I had previously made the `E3` suffix version (see docstring),
but this was optimized for learning a latent representation of the entire
sequence, which is effective for disentangled-like representation and likelihood
but less so for prediction, since optimizing the evidence of the **observed**
sequence appears to perform worse on future observations than an RNN trained
only to predict. For this model, the encoder learns a latent code that is
optimized for **predictive** performance, not **reconstructive** performance.
See the paper, and similar observations have previously been made -- see e.g.
Karl et al. 2016, Deep VB Filter.

As for the `Seq2SeqModel`, the posterior over the initial state of the generator
is learned variationally via the `init_enc` RNN with the final state passed to
the `init_post` layer which generates the mean and variance of the latent
posterior. Since the posterior is typically lower dimensional than the state of
the generator network, the latent variable is expanded via the `x0decoder`
network into the state of the `generator`. The MT variable is encoded with a
full-length sequence encoding via the `mt_enc` network, and the final state is
passed to the `mt_post` layer to generate the posterior mean and variance.
The MT variable enters the generator network either:

1. by predicting all parameters of the network via a `MTGRU_NoU` construction
(see docstring).
2. by forming a latent (constant) input at each time step, which functions as a
variable bias of the generator recurrent cell. This may be directly from the
latent variables, or via a decoder (`mtbias_decoder`).

Once these variables are sampled, the `generator` is iterated in an open-loop
fashion for `T_steps`, similar to the `Seq2SeqModel`. The penultimate slot `d`
is a record of the dimension of 𝒴, used currently just for `show`; and the
`mtbias_only` variable dictates whether either approach (1) or (2) is used as
above.

--------------------
This type can be called via
```julia
    m = MTSeqModel_Pred(...)
    m(y::Vector)
```
where `y` is a vector of slices of the sequence for the encoder. This means each
`y[i]` is an `AbstractMatrix` of size ``d \\times n_{\\text{batch}}`.

Options include `T_enc`, `T_steps` (already discussed), `reset` (which performs
a state reset before execution; true by default), and `stoch`, whether to
randomly sample the latent state, or use the state mean (default true).
"""
struct MTGRU_variational{U,A,B,V,W,S,N,M,T} <: MTDSModel
    enc_tform::U
    init_enc::A
    init_post::B
    mt_enc::A
    mt_post::B
    x0decoder::V
    mtbias_decoder::W
    generator::S
    dec_tform::N
    model_type::M
    logstd::T
    d::Int
    d_mt::Int
    mtbias_only::Bool
end

# Model utils
Flux.@treelike MTGRU_variational
model_type(m::MTGRU_variational{U,A,B,V,W,S,N,M,T}) where {U,A,B,V,W,S,N,M,T} = M

unpack(m::MTGRU_variational) = Flux.children(m)[1:9] # unpack struct; exclude `d`, `mtbias_only`.
unpack_inference(m::MTGRU_variational) = Flux.children(m)[1:5]
unpack_generative(m::MTGRU_variational) = (m.x0decoder, m.mtbias_decoder, m.generator, m.dec_tform)

load!(m::MTGRU_variational, fname::String) = load!(Flux.params(m), fname)

function Base.show(io::IO, l::MTGRU_variational)
    d_x0, d_mt = size(l.init_post.Dense1.W, 1), l.d_mt
    out_type = Dict(MTGRU_Single=>"Deterministic",
                    MTGRU_Double=>"Probabilistic")[model_type(l)]
    mttype = l.mtbias_only ? "Bias-only" : "Full-MT Recurrent Cell (" *
                    get_strtype_wo_params(l.generator) * ")"
    print(io, "MTGRU_variational(", mttype, ", d_x0=", d_x0, ", d_mt=", d_mt, ", ", out_type, ")")
end

################################################################################
##                                                                            ##
##                          Encoding the prior sequence                       ##
##                                                                            ##
################################################################################

function encode(m::MTGRU_variational, Y::AbstractArray{T}, U::AbstractArray{T}; T_enc=size(Y, 2),
                T_x0=T_enc, stochastic=true) where T
    encode(m, vcat(Y, U); T_enc=T_enc, T_x0=T_x0, stochastic=stochastic)
end

encode(m::MTGRU_variational, Y::MaskedArray{T}, U::AbstractArray{T}; T_enc=size(Y, 2), T_x0=T_enc,
       stochastic=true) where T = encode(m, vcat(stacked(Y), U); T_enc=T_enc, T_x0=T_x0,
                                         stochastic=stochastic)

function encode(m::MTGRU_variational, data_enc::AbstractArray; T_enc=size(data_enc,2), T_x0=T_enc,
                           stochastic=true, x0=nothing, z=nothing)
    enc_tform, init_enc, init_post, mt_enc, mt_post = unpack_inference(m)

    # RNN 1: Amortized inference for initial state, h0, of generator
    x0, μ_x0, σ_x0 = posterior_sample(init_enc, init_post, data_enc, T_enc, stochastic)

    # RNN 2: Amortized inference for z
    # Technically this should be conditioned on the sample x₀, nevertheless, the posterior
    # of x0 is usually very tight, and the `mt_enc` has access to the same information.
    z, μ_z, σ_z = posterior_sample(mt_enc, mt_post, data_enc, T_x0, stochastic)

    return (z, μ_z, σ_z), (x0, μ_x0, σ_x0)
end

# common superclass properties
is_amortized(m::MTGRU_variational) = true
is_amortized(m::MTGRU_variational{U,A,B,V,W,S,N,M,T}) where {A <: LookupTable, U,B,V,W,S,N,M,T} = false
has_x0_encoding(m::MTGRU_variational) = true

################################################################################
##                                                                            ##
##                                Forward Model                               ##
##                                                                            ##
################################################################################

function forward(m::MTGRU_variational, z, U; T_steps=size(U, 2))
    n_u, n_t, n_b = size(U)
    forward(m, zeros(eltype(U), size(m.x0decoder.W, 2), n_b), z, U; T_steps=T_steps)
end

function forward(m::MTGRU_variational, x0::AbstractVecOrMat{T}, z::AbstractVecOrMat{T},
        U::AbstractArray{T}; T_steps=size(U, 2)) where T <: AbstractFloat
    # note that T_steps is currently ignored (this is no great issue).
    @argcheck size(z, 2) == size(U, 3)

    # currently don't modulate the initial state.
    return _forward_given_latent(m, x0, z, U, T_steps)
end

    function _forward_given_latent(m::MTGRU_variational, x0::AbstractVector, z::AbstractVector,
        U::AbstractArray{T,2}, T_steps::Int) where T
        x0, z, U = ensure_matrix(x0), ensure_matrix(z), ensure_tensor(U)
        return dropdims(_forward_given_latent(m, x0, z, U, T_steps), dims=3)
    end

    function _forward_given_latent(m::MTGRU_variational, x0::AbstractMatrix, z::AbstractMatrix,
                                U::AbstractArray{T,3}, T_steps::Int) where T
        x0decoder, mtbias_decoder, gen_model, dec_tform = unpack_generative(m)

        if dec_tform.layers[2] isa MTDenseGenerator
            _mtdec = dec_tform.layers[2](z)
            decoder(x) = let l1=dec_tform.layers[1](x); _mtdec(l1) + dec_tform.layers[3](l1); end
        else
            decoder = dec_tform
        end

        if m.mtbias_only
            gen_model.state = x0decoder(x0)
            mtbias = mtbias_decoder(z)   # if linear, this is just the identity ∵ GRUCell.Wi
            yhat_ts = [unsqueeze(decoder(gen_model(mtbias)),2) for t in 1:T_steps]
        else
            # Get GRU models from samples from the MT GRU model
            posterior_grus = gen_model(z) # output: BatchedGRU (def. above)
            # Set state (batch-wise) to x0 sample
            posterior_grus.state = x0decoder(x0)  # 2×linear d_x0 × nbatch → d_x × nbatch
            # Run generative model

            yhat_ts = [unsqueeze(decoder(posterior_grus(U[:,tt,:])),2) for tt in 1:T_steps]
        end
        return hcat(yhat_ts...)
    end

function (m::MTGRU_variational)(x0::AbstractVecOrMat, z::AbstractVecOrMat, U::AbstractArray;
          T_steps=size(U,2))
    _forward_given_latent(m, x0, z, U, T_steps)
end

# function _posterior_sample(enc, dec, input::Vector, T_max, stochastic=true, input_cat=nothing)
#     Flux.reset!(enc)
#     if input_cat === nothing
#         for tt = 1:T_max; enc(input[tt]); end
#     else
#         d, nbatch = size(input[1])
#         if nbatch == 1 && nbatch != size(input_cat,2)
#             gmove = Flux.has_cuarrays() && input_cat isa Flux.CuArray ? gpu : identity
#             input_expander = gmove(ones(Float32, 1, size(input_cat, 2)))
#         else
#             input_expander = 1
#         end
#         for tt = 1:T_max; enc(vcat(input[tt] * input_expander, input_cat)); end
#     end
#     μ_, σ_ = dec(enc.state)
#     n, d = size(μ_)
#     smp = μ_ + randn_repar(σ_, n, d, stochastic)
#     return smp, μ_, σ_
# end

################################################################################
##                                                                            ##
##                             Model "Constructor"                            ##
##                                                                            ##
################################################################################

function create_mtgru(d_u, d_x, d_y, d_enc_state, d_enc_x0, d_mt; encoder=:GRU,
    out_heads=1, d_hidden=d_x, mtbias_only=false, d_hidden_mt=32, mt_is_linear=true,
    fixb=false, fixb_version=nothing, spherical_var=false, lkp_init_std_logit=-2f0)

    @argcheck out_heads in 1:2
    @assert !(fixb && mtbias_only) "cannot do mtbias AND fix 'b'"
    @assert !fixb || something(fixb_version, 1) ∈ [1,2] "fixb_version must be in [1,2]."

    if fixb && fixb_version === nothing
        @warn "fixed bias chosen, but `fixb_version` not specified. Defaulting to v1."
        fixb_version = 1
    end

    # ENCODER TRANSFORM FROM OBS SPACE => ENC SPACE
    tform_enc = identity
    d_rnn_in = d_y+d_u
    model_type = out_heads == 1 ? MTGRU_Single() : MTGRU_Double()

    # 2× ENCODERS for inference of (a) initial state (init_enc); (b) sequence level MT var (mt_enc)
    @argcheck encoder in [:LSTM, :GRU, :Bidirectional, :LookupTable]

    if encoder == :LookupTable
        # args 2/3: initial posterior noise for mu/logstd
        init_enc = LookupTable(d_enc_x0, 0.01f0, lkp_init_std_logit)
        mt_enc = LookupTable(d_mt, 0.01f0, lkp_init_std_logit)
    else
        rnn_constructor = Dict(:LSTM=>LSTM, :GRU=>GRU, :Bidirectional=>BRNNenc)[encoder]
        init_enc = rnn_constructor(d_rnn_in, d_enc_state)
        mt_enc = rnn_constructor(d_rnn_in, d_enc_state)
        (encoder == :Bidirectional) && (d_enc_state *= 2)
    end

    # POSTERIOR OVER LVM CORR. TO (a), (b). Ignored if using LookupTable
    init_post = MultiDense(Dense(d_enc_state, d_enc_x0, identity), Dense(d_enc_state, d_enc_x0, σ))
    init_post.Dense2.b.data .= -2   # initialize posteriors to be low variance
    mt_post = MultiDense(Dense(d_enc_state, d_mt, identity), Dense(d_enc_state, d_mt, σ))
    mt_post.Dense2.b.data .= -2

    # decode from LV (a) --> size of generative hidden state
    x0decoder = Dense(d_enc_x0, d_x, identity)

    # -------------- GENERATIVE MODEL -----------------
    ###################################################

    if !mtbias_only
        d_out = !fixb ? 3*d_x*(d_x+d_u+1) : (fixb_version == 1 ? 3*d_x*(d_x+d_u) :
            3*d_x*(d_x+d_u) + 2*d_x)
        par_gen_net = mt_is_linear ? mlp(d_mt, d_out) : mlp(d_mt, d_hidden_mt,
                                                                    d_out; activation=tanh)
        if !fixb
            gen_rnn = MTGRU(d_u, d_x, par_gen_net)
        else
            if fixb_version == 1
                MTGRU_b = Flux.param(randn(Float32, d_x*3))
                gen_rnn = MTGRU_fixb(d_u, d_x, par_gen_net, MTGRU_b)
            elseif fixb_version == 2
                MTGRU_b = Flux.param(randn(Float32, d_x))
                gen_rnn = MTGRU_fixb_v2(d_u, d_x, par_gen_net, MTGRU_b)
            else
                error("fixb_version: Unreachable error.")
            end
        end
        mtbias_decoder = identity
    else
        mtbias_decoder = mt_is_linear ? identity : Dense(d_mt, d_hidden_mt, swish)
        _d_in = mt_is_linear ? d_mt : d_hidden_mt
        gen_rnn = GRU(_d_in, d_x)
    end

    if out_heads == 1
        # decoder = Chain(Dense(d_x, d_hidden, relu), Dense(d_hidden, d_y, identity))
        decoder = Chain(Dense(d_x, d_hidden, relu),
                        MTDenseGenerator(d_mt, d_hidden, d_y, identity),
                        Dense(d_hidden, d_y, identity))   # this is super hacky -- see forward pass
    elseif out_heads == 2
        @error "out_heads=2 reqs a little re-engineering of the forward pass due to the MTDense"*
        " layer. Really we need to make a generator for the whole Chain of the decoder, not simply"*
        " the MTDense part."
        decoder = Chain(Dense(d_x, d_hidden, relu), MultiDense(
            MTDenseGenerator(d_mt, d_hidden, d_y, identity),
            MTDenseGenerator(d_mt, d_hidden, d_y, identity)))
    end

    emission_std = spherical_var ? Tracker.param([0.0f0]) : Tracker.param(zeros(Float32, d_y))

    m = MTGRU_variational(tform_enc, init_enc, init_post, mt_enc, mt_post,
        x0decoder, mtbias_decoder, gen_rnn, decoder, model_type, emission_std,
        d_y, d_mt, mtbias_only)

    return m
end

"""
    create_model_opt_dict(...)
Takes same arguments as `create_model` but just wraps up arguments in a Dict for
saving out as JSON/YAML/BSON etc.
"""
function create_model_opt_dict(d_x, d_x0, d_y, d_enc_state, d_mt, d_chaos=0; encoder=:GRU,
    cnn=false, out_heads=1, d_hidden=d_x, mtbias_only=false, d_hidden_mt=32,
    mt_is_linear=true, decoder_fudge_layer=false, model_purpose=:llh, fixb=false,
    fixb_version=nothing)

    Dict("d_x"=>d_x, "d_x0"=>d_x0, "d_y"=>d_y, "d_enc_state"=>d_enc_state, "d_mt"=>d_mt, "d_chaos"=>d_chaos,
    "encoder"=>encoder, "cnn"=>cnn, "out_heads"=>out_heads, "d_hidden"=>d_hidden, "mtbias_only"=>mtbias_only,
    "d_hidden_mt"=>d_hidden_mt, "mt_is_linear"=>mt_is_linear, "decoder_fudge_layer"=>decoder_fudge_layer,
    "model_purpose"=>model_purpose, "fixb"=>fixb, "fixb_version"=>fixb_version)
end

"""
    load_model_from_def(ymlfile)
We have YAML files saved with metadata associated with the model parameters
in the /data folder. These YAML files point to the relevant parameters too, and
so suffices to construct the model and load the parameters. That is the purpose
of this function.
"""
function load_model_from_def(ymlfile::String)
    model_details = YAML.load_file(ymlfile)
    filename = model_details["filename"]
    constructor = model_details["constructor"]
    @argcheck constructor == "mtmodel.create_model"
    D = model_details["model_def"]
    D["encoder"] = Symbol(D["encoder"])
    D["model_purpose"] = Symbol(D["model_purpose"])
    m = create_model(D["d_x"], D["d_x0"], D["d_y"], D["d_enc_state"], D["d_mt"], D["d_chaos"];
        encoder=D["encoder"], cnn=D["cnn"], out_heads=D["out_heads"], d_hidden=D["d_hidden"],
        mtbias_only=D["mtbias_only"], d_hidden_mt=D["d_hidden_mt"], mt_is_linear=D["mt_is_linear"],
        decoder_fudge_layer=D["decoder_fudge_layer"], model_purpose=D["model_purpose"], fixb=D["fixb"],
        fixb_version=D["fixb_version"])
    modelutils.load!(m, filename)
    return m
end

# # Example saving model info
# YAML.write_file("saved_models/mtgru_video_fixbgrp_pred_450.yml", Dict(
#     "filename"=>"data/mtgru_video_fixbgrp_pred_450.bson",
#     "description"=>"MT GRU model with posterior mean GROUPED BY SEQ IDENTITY fixed (v2) bias. 450x2 epochs. T=80, Tenc=20.",
#     "constructor"=>"mtmodel.create_model",
#     "model_def"=>mtmodel.create_model_opt_dict(64, 6, 4, 40, 2; encoder=:GRU, cnn=true, d_hidden=128,
#     model_purpose=:pred, mtbias_only=false, decoder_fudge_layer=true, fixb=true, fixb_version=2)
# ))

####################################################################################################
#                                                                                #
#     Overloading some objective subroutines in order to accommodate both z and x0 latent vars     #
#                                                                                #
####################################################################################################

function elbo_w_kl(m::MTGRU_variational, y::AbstractArray, u::AbstractArray,
                   data_enc::AbstractArray; kl_coeff=1.0f0, stochastic=true, logstd_prior=nothing,
                   T_enc=size(u, 2))
    nbatch = size(y,3)
    zstuff, x0stuff = encode(m, data_enc; stochastic=stochastic, T_enc=T_enc)
    zsmp, zμ, zσ = zstuff
    x0smp, x0μ, x0σ = x0stuff
    ŷ = m(x0smp, zsmp, u)
    nllh = _gauss_nllh_batchavg(y, ŷ, m.logstd)
    kl = stochastic * kl_coeff * (-0.5f0 * sum(1 .+ 2*log.(zσ) - zμ.*zμ - zσ.*zσ)) / nbatch
    kl += stochastic * kl_coeff * (-0.5f0 * sum(1 .+ 2*log.(x0σ) - x0μ.*x0μ - x0σ.*x0σ)) / nbatch
    # Prior regularization of emission stdev.
    if !(logstd_prior === nothing)
        nllh += sum(x->x^2, (m.logstd .- logstd_prior[1])./logstd_prior[2])/2
    end
    nllh + kl, kl
end

"""
    forward_multiple_z(m::MTDSModel, x0, z, u::AbstractMatrix; maxbatch=200)
Performs multiple forward passes, partitioned according to the `maxbatch` argument, with a single
`x0` and `u`, but different `z` values. Returns a matrix of size d_y × T × n_b.
"""
function forward_multiple_z(m::MTGRU_variational, x0, z, u::AbstractMatrix; maxbatch=200)
    M = size(z, 2)
    ŷs = map(Iterators.partition(1:M, maxbatch)) do cbatch
        u_rep = repeat(u, outer=(1, 1, length(cbatch)));
        x0_rep = repeat(x0, outer=(1, length(cbatch)));
        m(x0_rep, z[:, cbatch], u_rep)
    end
    cat(ŷs..., dims=3)   # problematic if length(ŷs) is too large, but there's no nice `reduce` version in Julia yet.
end


text:
using Robotlib
using Robotlib: Rt2T, T2R, T2t
function leds_to_tf(L::AbstractMatrix, usemean = false)
    size(L) == (3, 3) || throw(ArgumentError("Expected a size 3×3 matrix"))

    if usemean
        m = mean(L, dims = 2)[:]
        v1 = normalize(L[:,1] - m)
        v2 = normalize(L[:,2] - m)
    else
        m = L[:,1]
        v1  = normalize(L[:,2]-m)
        v2  = normalize(L[:,3]-m)
    end
    v3 = normalize(cross(v1, v2))
    v2 = normalize(cross(v3, v1))
    R = [v1 v2 v3]
    det(R) < 0 && (R = [v1 v2 -v3])
    Rt2T(R, m)
end

function leds_to_tf(L::AbstractVector{<:AbstractMatrix})
    length(L) == 3 || throw(ArgumentError("Expected three matrices"))
    NL, N = size(L[1])
    NL > N && @warn "More LEDs than number of data points, are you sure you should not transpose the L matrices?"
    T = Matrix{Float64}[]
    for i = 1:N
        D = [L[1][:, i] L[2][:, i] L[3][:, i]]
        all(isfinite, D) || continue
        push!(T, leds_to_tf(D))
    end
    T
end

"""
    T_RB_T, T_TF_TCP, normDist = calibAXYB(POSES, MEASURED; ortho=true, estimator = \\, refine = true)

This functions solves the problem AX = YB

- `POSES` is always the position of the tool frame from the robot FK
- `MEASURED` is the position of the object mounted on the tool frame, for
    instance if measuring with the Nikon it is the diod frame. If measuring
    with a laser scanner mounted on the TF, it is the position of the SCANNER
    in the measured object frame, not the object in the scanner frame!
- `ortho` is a flag indicating if found matrices are to be orthogonalized
- `refine`: run a second, local optimization routine to improve the accuracy (recommended)
"""
function calibAXYB(POSES, MEASURED; ortho = true, estimator = \, refine = true)
    length(POSES) == length(MEASURED) ||
        throw(ArgumentError("POSES and MEASURED must have the same length"))
    T = promote_type(eltype(POSES[1]), eltype(MEASURED[1]))
    A = Matrix{T}[]
    B = Vector{T}[]
    I = Robotlib.I3
    Z = zeros(3, 3)

    for i = 1:length(POSES)
        RA = T2R(POSES[i])
        BB = MEASURED[i][1:3, 1:4]
        RB = T2R(BB)
        TB = T2t(BB)

        At = [
            RA Z Z Z -I*RB[1, 1] -I*RB[2, 1] -I*RB[3, 1] Z
            Z RA Z Z -I*RB[1, 2] -I*RB[2, 2] -I*RB[3, 2] Z
            Z Z RA Z -I*RB[1, 3] -I*RB[2, 3] -I*RB[3, 3] Z
            Z Z Z RA -I*TB[1] -I*TB[2] -I*TB[3] -I
        ]

        b = [zeros(9); -T2t(POSES[i])]
        push!(A, At)
        push!(B, b)
    end
    A = reduce(vcat, A)
    B = reduce(vcat, B)

    w = estimator(A, B)
    T_RB_T = reshape(w[13:24], 3, 4)
    T_RB_T = [T_RB_T; [0 0 0 1]]
    T_TF_TCP = reshape(w[1:12], 3, 4)
    T_TF_TCP = [T_TF_TCP; [0 0 0 1]]
    if ortho
        T_RB_T = Robotlib.orthonormal(T_RB_T)
        T_TF_TCP = Robotlib.orthonormal(T_TF_TCP)
        # Second optimization
        B2 = B - A[:, 1:9] * vec(T_TF_TCP[1:3, 1:3]) - A[:, 13:21] * vec(T_RB_T[1:3, 1:3])
        A2 = A[:, [10:12; 22:24]]
        w2 = estimator(A2, B2)
        T_TF_TCP[1:3, 4] = w2[1:3]
        T_RB_T[1:3, 4] = w2[4:6]
    end

    # Eval
    YN = similar(MEASURED)
    MX = similar(MEASURED)
    normDist = zeros(length(MEASURED))
    function cost(T_RB_T, T_TF_TCP)
        for i = 1:length(POSES)
            YN[i] = Robotlib.orthonormal(T_RB_T * MEASURED[i])
            MX[i] = Robotlib.orthonormal(POSES[i] * T_TF_TCP)
            dist = (YN[i][1:3, 4] - MX[i][1:3, 4])'
            normDist[i] = norm(dist)
        end
        sum(normDist)
    end
    cost(T_RB_T, T_TF_TCP)
    if refine
        p0 = zeros(12)
        adapter = function(p)
            T_RB_Ti = T_RB_T*Rt2T(rpy2R(p[4:6]), p[1:3])
            T_TF_TCPi = T_TF_TCP*Rt2T(rpy2R(p[10:12]), p[7:9])
            T_RB_Ti, T_TF_TCPi
        end
        cost(p) = cost(adapter(p)...)
        res = Optim.optimize(
            cost,
            p0,
            BFGS(),
            Optim.Options(
                store_trace       = true,
                show_trace        = true,
                show_every        = 2,
                iterations        = 100,
                allow_f_increases = false,
                time_limit        = 100,
                g_tol             = 1e-12,
            ),
        )
        T_RB_T, T_TF_TCP = adapter(res.minimizer)
        cost(T_RB_T, T_TF_TCP)
    end

    T_RB_T, T_TF_TCP, normDist
end

random_SE3(σt=1) = Rt2T(orthonormal(randn(3,3)), σt*randn(3))

"""
    T_RB_T, T_TF_S, POSES, MEASURED, LEDs = simulateCalibration_AXYB(N; σt=0.1, σT=1, σy=0)

σt denots the std of the translation of the sensor frame. σT is the std for the translation of the camera to the robot. σy is measurement noise
"""
function simulateCalibration_AXYB(N=100; σt=0.1, σT=1, σy=0)
    POSES = [random_SE3(0.5*σT) for _ = 1:N]
    T_RB_T = random_SE3(σT)
    T_TF_S = random_SE3(σt)
    MEASURED = [trinv(T_RB_T)*P*T_TF_S for P in POSES]
    LEDs = map(1:3) do l
        L = map(1:N) do i
            t = T2t(MEASURED[i])
            R = T2R(MEASURED[i])
            t .+ (l == 1 ? 0 : σt*R[:,l-1]) # origin + one of the basis vectors
        end
        reduce(hcat, L) .+ σy .* randn.()
    end
    T_RB_T, T_TF_S, POSES, MEASURED, LEDs
end

max_forks_repo_name: baggepinnen/Robotlib


text:
include("random.jl")
include("matrix.jl")

HEIGHT = 20
WIDTH = 80
playground = matrix::make_matrix(2, [HEIGHT, WIDTH], 0)
playground_temp = matrix::make_matrix(2, [HEIGHT, WIDTH], 0)

"""Prints the playground"""
function print_playground()
    for i = 0 : (HEIGHT - 1)
        for j = 0 : (WIDTH - 1)
            if (playground[i][j] == 1)
                print("0")
            else
                print(".")
            end
        end
        println()
    end
end

"""Initializes the playground randomly."""
function initialize_playground()
    somme = 0
    for i = 0 : (HEIGHT - 1)
        for j = 0 : (WIDTH - 1)
            playground[i][j] = random::randint(0, 2)
        end
    end
end

"""returns the number of alive neighbours of point (h,w)."""
function get_neigh(h :: Int64, w :: Int64) :: Int64
    voisins = 0
    if (h >= 1 && playground[h - 1][w] == 1)
        voisins = voisins + 1
    end
    if (w >= 1 && playground[h][w - 1] == 1)
        voisins = voisins + 1
    end
    if (h >= 1 && w >= 1 && playground[h - 1][w - 1] == 1)
        voisins = voisins + 1
    end
    if (h < HEIGHT - 1 && playground[h + 1][w] == 1)
        voisins = voisins + 1
    end
    if (w < WIDTH - 1 && playground[h][w + 1] == 1)
        voisins = voisins + 1
    end
    if (h < HEIGHT - 1 && w < WIDTH - 1 && playground[h + 1][w + 1] == 1)
        voisins = voisins + 1
    end
    if (h < HEIGHT - 1 && w >= 1 && playground[h + 1][w - 1] == 1)
        voisins = voisins + 1
    end
    if (h >= 1 && w < WIDTH - 1 && playground[h - 1][w + 1] == 1)
        voisins = voisins + 1
    end
    return voisins
end

"""Computes a single step for a single pixel."""
function handle_pixel(h :: Int64, w :: Int64)
    voisins = get_neigh(h, w)
    if (playground[h][w] == 0) && (voisins == 3)
        playground_temp[h][w] = 1
    elseif (playground[h][w] == 1) && (voisins < 2 || voisins > 3)
        playground_temp[h][w] = 0
    end
end

"""copies back from playground_temp to playground."""
function copy_back()
    for i = 0 : (HEIGHT - 1)
        for j = 0 : (WIDTH - 1)
            playground[i][j] = playground_temp[i][j]
        end
    end
end

function step()
    for i = 0 : (HEIGHT - 1)
        for j = 0 : (WIDTH - 1)
            handle_pixel(i, j)
        end
    end
    copy_back()
    print_playground()
end

function run()
    initialize_playground()
    while true
        a = input_int() #sert à temporiser
        step()
    end
end


text:
push!(LOAD_PATH, "./src")

using Test
using YaDistances: cosine_similarity

@testset "Basic tests" begin
  @test cosine_similarity([1., 1, 1], [2., 2, 2]) ≈ 1. ## "same direction"
	@test cosine_similarity([-1., -1], [2., 2]) ≈ -1.    ## "opposite direction"
	@test cosine_similarity([1., 0], [0., 1]) ≈ 0.       ## "orthogonal"
end

max_stars_repo_path: Data Science From Scratch/test/ya_distances_test.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
function initialize!(
                status,
                parameters::AbstractParameters,
                problem,
                information,
                options,
                args...;
                kargs...
        )

    m = string(typeof(parameters))
    println("Method $m is calling the template function.")
    println("Redefine your own function using template")
    error("Method $m is calling the template function.")
    return State(0.0, zeros(0))
end

function update_state!(
        status,
        parameters::AbstractParameters,
        problem,
        information,
        options,
        args...;
        kargs...
)

    return
end

function final_stage(
        status,
        parameters::AbstractParameters,
        problem,
        information,
        options,
        args...;
        kargs...
)
    return
end

function stop_criteria!(
        status,
        parameters::AbstractParameters,
        problem,
        information,
        options,
        args...;
        kargs...
    )

    status.stop = status.stop || call_limit_stop_check(status, information, options) ||
                  iteration_stop_check(status, information, options)  ||
                  time_stop_check(status, information, options)       ||
                  accuracy_stop_check(status, information, options)
    return
end


text:
using IHT, PLINK, RegressionTools
(xbed, ybed) = read_plink_data("test")
kevin_result = L0_reg(xbed, ybed, 2)