text:
function has_subgraph(tait::Tait, m::Match)
    all(has_vertex(tait, v) for v in m.vertices) || return false
    all(has_half_edge(tait, he) for he in m.half_edges) || return false
    return true
end

check(tait::Tait, m::Match) = false

function check(tait::Tait, m::Match{:string_genus})
    v = m.vertices[1]
    return v in tait.genuses
end

function check(tait::Tait, m::Match{:yang_baxter_star})
    has_subgraph(tait, m) || return false
    v = m.vertices[1]

    (is_open_vertex(tait, v) || is_genus(tait, v)) && return false
    hes = trace_vertex(tait, v)
    if length(hes) == 3
        nbs = Set(dst(tait, he) for he in hes)
        if length(nbs) == 3 && all(x->!is_open_vertex(tait, dst(tait, x)), hes)
            return Set(hes) == Set(m.half_edges)
        end
    end
    return false
end

function check(tait::Tait, m::Match{:yang_baxter_triangle})
    has_subgraph(tait, m) || return false
    hes = m.half_edges

    length(hes) == 3 || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false

    f = face(tait, hes[1])
    f != 0 || return false
    f == face(tait, hes[2]) == face(tait, hes[3]) || return false
    return length(trace_face(tait, f)) == 3
end

function check(tait::Tait, m::Match{:charge_rm_v})
    has_subgraph(tait, m) || return false
    hes = m.half_edges
    v = m.vertices[1]

    is_open_vertex(tait, v) && return false
    src(tait, hes[1]) == v || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false
    for i = lastindex(hes):2
        σ_inv(tait, hes[i]) == hes[i-1] || return false
        p = quon_param(tait, hes[i])
        is_para_pi(p) || return false
    end
    p1 = quon_param(tait, hes[1])
    return is_para_pi(p1)
end

function check(tait::Tait, m::Match{:charge_rm_f})
    has_subgraph(tait, m) || return false
    hes = m.half_edges

    face(tait, hes[1]) != 0 || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false
    for i = 1:(lastindex(hes)-1)
        next(tait, hes[i]) == hes[i+1] || return false
        p = quon_param(tait, hes[i])
        is_orth_pi(p) || return false
    end
    p_end = quon_param(tait, hes[end])
    return is_orth_pi(p_end)
end

function check(tait::Tait, m::Match{:z_fusion})
    has_subgraph(tait, m) || return false
    hes = m.half_edges

    length(hes) == 2 || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false
    next(tait, hes[1]) == hes[2] || return false
    next(tait, hes[2]) == hes[1] || return false
    return face(tait, hes[1]) != 0
end

function check(tait::Tait, m::Match{:x_fusion})
    has_subgraph(tait, m) || return false
    v = m.vertices[1]
    hes = m.half_edges

    is_genus(tait, v) && return false
    length(hes) == 2 || return false
    any(is_open_half_edge(tait, he) for he in hes) && return false
    src(tait, hes[1]) == v || return false
    σ_inv(tait, hes[1]) == hes[2] || return false
    return σ_inv(tait, hes[2]) == hes[1]
end

function check(tait::Tait, m::Match{:perm_rz})
    has_subgraph(tait, m) || return false
    (v, g2) = m.vertices
    he1, he2 = m.half_edges

    g1 = dst(tait, he1)
    is_open_half_edge(tait, he1) && return false
    src(tait, he2) == g2 || return false
    (is_genus(tait, g1) && is_genus(tait, g2)) || return false
    is_genus(tait, v) && return false
    hes_v = trace_vertex(tait, v)
    fs_v = Set(face(tait, he) for he in hes_v)
    return face(tait, he2) in fs_v
end

function check(tait::Tait, m::Match{:identity})
    has_subgraph(tait, m) || return false
    he = m.half_edges[1]
    return is_para_zero(quon_param(tait, he)) || is_orth_zero(quon_param(tait, he))
end

function check(tait::Tait, m::Match{:genus_fusion})
    has_subgraph(tait, m) || return false
    g1, g2 = m.vertices
    he1, he2 = m.half_edges

    (is_genus(tait, g1) && is_genus(tait, g2)) || return false
    (src(tait, he1) == g1 && src(tait, he2) == g2) || return false
    face(tait, he1) == face(tait, he2) || return false
    return face(tait, he1) != 0
end

function check(tait::Tait, m::Match{:swap_genus})
    has_subgraph(tait, m) || return false
    vs = m.vertices
    hes = m.half_edges

    all(is_genus(tait, g) for g in vs[1:4]) || return false
    all(!is_genus(tait, v) for v in vs[5:9]) || return false
    v0 = vs[9]
    all(src(tait, he) == v0 for he in hes) || return false
    all(length(trace_face(tait, face(tait, he))) == 4 for he in hes) || return false
    all(check_swap_parameters(tait, v) for v in vs[5:9]) || return false
    g1 = vs[1]
    return !is_genus_connected_to_open_edge(tait, g1)
end


text:
Tables.istable(::Type{<:AbstractDataFrame}) = true
Tables.columnaccess(::Type{<:AbstractDataFrame}) = true
Tables.columns(df::AbstractDataFrame) = df
Tables.rowaccess(::Type{<:AbstractDataFrame}) = true
Tables.rows(df::AbstractDataFrame) = eachrow(df)
Tables.rowtable(df::AbstractDataFrame) = Tables.rowtable(Tables.columntable(df))
Tables.namedtupleiterator(df::AbstractDataFrame) =
    Tables.namedtupleiterator(Tables.columntable(df))

function Tables.columnindex(df::Union{AbstractDataFrame, DataFrameRow}, idx::Symbol)
    ind = index(df)
    if ind isa Index
        return get(ind.lookup, idx, 0)
    else
        parent_ind = ind.parent
        loc = get(parent_ind.lookup, idx, 0)
        return loc == 0 || loc > length(ind.remap) ? 0 : max(0, ind.remap[loc])
    end
end

Tables.columnindex(df::Union{AbstractDataFrame, DataFrameRow}, idx::AbstractString) =
    columnindex(df, Symbol(idx))

Tables.schema(df::AbstractDataFrame) = Tables.Schema(propertynames(df), eltype.(eachcol(df)))
Tables.materializer(df::AbstractDataFrame) = DataFrame

Tables.getcolumn(df::AbstractDataFrame, i::Int) = df[!, i]
Tables.getcolumn(df::AbstractDataFrame, nm::Symbol) = df[!, nm]

Tables.getcolumn(dfr::DataFrameRow, i::Int) = dfr[i]
Tables.getcolumn(dfr::DataFrameRow, nm::Symbol) = dfr[nm]

getvector(x::AbstractVector) = x
getvector(x) = [x[i] for i = 1:length(x)]
# note that copycols is ignored in this definition (Tables.CopiedColumns implies copies have already been made)
fromcolumns(x::Tables.CopiedColumns, names; copycols::Bool=true) =
    DataFrame(AbstractVector[getvector(Tables.getcolumn(x, nm)) for nm in names],
              Index(names),
              copycols=false)
fromcolumns(x, names; copycols::Bool=true) =
    DataFrame(AbstractVector[getvector(Tables.getcolumn(x, nm)) for nm in names],
              Index(names),
              copycols=copycols)

function DataFrame(x::T; copycols::Bool=true) where {T}
    if !Tables.istable(x)
        if x isa AbstractVector && all(col -> isa(col, AbstractVector), x)
            return DataFrame(Vector{AbstractVector}(x), copycols=copycols)
        elseif (x isa AbstractVector || x isa Tuple) &&
            all(v -> v isa Pair{Symbol, <:AbstractVector}, x)
            return DataFrame(AbstractVector[last(v) for v in x], [first(v) for v in x],
                             copycols=copycols)
        end
    end
    cols = Tables.columns(x)
    names = collect(Symbol, Tables.columnnames(cols))
    return fromcolumns(cols, names, copycols=copycols)
end

function Base.append!(df::DataFrame, table; cols::Symbol=:setequal,
                      promote::Bool=(cols in [:union, :subset]))
    if table isa Dict && cols == :orderequal
        throw(ArgumentError("passing `Dict` as `table` when `cols` is equal to " *
                            "`:orderequal` is not allowed as it is unordered"))
    end
    append!(df, DataFrame(table, copycols=false), cols=cols, promote=promote)
end

# This supports the Tables.RowTable type; needed to avoid ambiguities w/ another constructor
DataFrame(x::AbstractVector{NamedTuple{names, T}}; copycols::Bool=true) where {names, T} =
    fromcolumns(Tables.columns(Tables.IteratorWrapper(x)), collect(names), copycols=false)

Tables.istable(::Type{<:Union{DataFrameRows,DataFrameColumns}}) = true
Tables.columnaccess(::Type{<:Union{DataFrameRows,DataFrameColumns}}) = true
Tables.rowaccess(::Type{<:Union{DataFrameRows,DataFrameColumns}}) = true
Tables.columns(itr::Union{DataFrameRows,DataFrameColumns}) = Tables.columns(parent(itr))
Tables.rows(itr::Union{DataFrameRows,DataFrameColumns}) = Tables.rows(parent(itr))
Tables.schema(itr::Union{DataFrameRows,DataFrameColumns}) = Tables.schema(parent(itr))
Tables.rowtable(itr::Union{DataFrameRows,DataFrameColumns}) = Tables.rowtable(parent(itr))
Tables.namedtupleiterator(itr::Union{DataFrameRows,DataFrameColumns}) =
    Tables.namedtupleiterator(parent(itr))
Tables.materializer(itr::Union{DataFrameRows,DataFrameColumns}) =
    Tables.materializer(parent(itr))

Tables.getcolumn(itr::Union{DataFrameRows,DataFrameColumns}, i::Int) =
    Tables.getcolumn(parent(itr), i)
Tables.getcolumn(itr::Union{DataFrameRows,DataFrameColumns}, nm::Symbol) =
    Tables.getcolumn(parent(itr), nm)

IteratorInterfaceExtensions.getiterator(df::AbstractDataFrame) =
    Tables.datavaluerows(Tables.columntable(df))
IteratorInterfaceExtensions.isiterable(x::AbstractDataFrame) = true
TableTraits.isiterabletable(x::AbstractDataFrame) = true

ext: jl


text:
using BinaryBuilder, Pkg

function configure(version_offset, min_julia_version)
    # The version of this JLL is decoupled from the upstream version.
    # Whenever we package a new upstream release, we initially map its
    # version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
    # So for example version 2.6.3 would become 200.600.300.

    name = "PROJ"
    upstream_version = v"7.2.1"
    version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                            upstream_version.minor * 100 + version_offset.minor,
                            upstream_version.patch * 100 + version_offset.patch)

    # Collection of sources required to build PROJ
    sources = [
        ArchiveSource("https://download.osgeo.org/proj/proj-$upstream_version.tar.gz",
            "b384f42e5fb9c6d01fe5fa4d31da2e91329668863a684f97be5d4760dbbf0a14"),
    ]

    # Bash recipe for building across all platforms
    script = raw"""
    cd $WORKSPACE/srcdir/proj-*/

    # Get rid of target sqlite3, to avoid it's picked up by the build system
    rm "${bindir}/sqlite3${exeext}"

    if [[ ${target} == *mingw* ]]; then
        SQLITE3_LIBRARY=${libdir}/libsqlite3-0.dll
        CURL_LIBRARY=${libdir}/libcurl-4.dll
        TIFF_LIBRARY_RELEASE=${libdir}/libtiff-5.dll
    else
        SQLITE3_LIBRARY=${libdir}/libsqlite3.${dlext}
        CURL_LIBRARY=${libdir}/libcurl.${dlext}
        TIFF_LIBRARY_RELEASE=${libdir}/libtiff.${dlext}
    fi

    mkdir build
    cd build
    cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release \
        -DBUILD_SHARED_LIBS=ON \
        -DBUILD_TESTING=OFF \
        -DSQLITE3_INCLUDE_DIR=${includedir} \
        -DSQLITE3_LIBRARY=$SQLITE3_LIBRARY \
        -DCURL_INCLUDE_DIR=${includedir} \
        -DCURL_LIBRARY=$CURL_LIBRARY \
        -DTIFF_INCLUDE_DIR=${includedir} \
        -DTIFF_LIBRARY_RELEASE=$TIFF_LIBRARY_RELEASE \
        ..
    make -j${nproc}
    make install
    """

    platforms = expand_cxxstring_abis(supported_platforms())

    # The products that we will ensure are always built
    products = [
        LibraryProduct(["libproj", "libproj_$(upstream_version.major)_$(upstream_version.minor)"], :libproj),

        # Excecutables
        ExecutableProduct("cct", :cct),
        ExecutableProduct("cs2cs", :cs2cs),
        ExecutableProduct("geod", :geod),
        ExecutableProduct("gie", :gie),
        ExecutableProduct("proj", :proj),
        ExecutableProduct("projinfo", :projinfo),
        ExecutableProduct("projsync", :projsync),

        # complete contents of share/proj, must be kept up to date
        FileProduct(joinpath("share", "proj", "CH"), :ch),
        FileProduct(joinpath("share", "proj", "GL27"), :gl27),
        FileProduct(joinpath("share", "proj", "ITRF2000"), :itrf2000),
        FileProduct(joinpath("share", "proj", "ITRF2008"), :itrf2008),
        FileProduct(joinpath("share", "proj", "ITRF2014"), :itrf2014),
        FileProduct(joinpath("share", "proj", "nad.lst"), :nad_lst),
        FileProduct(joinpath("share", "proj", "nad27"), :nad27),
        FileProduct(joinpath("share", "proj", "nad83"), :nad83),
        FileProduct(joinpath("share", "proj", "other.extra"), :other_extra),
        FileProduct(joinpath("share", "proj", "proj.db"), :proj_db),
        FileProduct(joinpath("share", "proj", "proj.ini"), :proj_ini),
        FileProduct(joinpath("share", "proj", "world"), :world),
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = [
        # Host SQLite needed to build proj.db
        HostBuildDependency("SQLite_jll"),
        Dependency("SQLite_jll"),
        # TODO: v4.3.0 is available, use that next time
        Dependency("Libtiff_jll"; compat="4.1.0"),
        Dependency("Zlib_jll"),
    ]

    jll_stdlibs = Dict(
        v"1.3" => [
            Dependency("LibCURL_jll", v"7.71.1"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll", v"1.9.0"),
            Dependency("MbedTLS_jll", v"2.16.8"),
            Dependency("nghttp2_jll", v"1.40.0"),
        ],
        v"1.6" => [
            Dependency("LibCURL_jll"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll"),
            Dependency("MbedTLS_jll", v"2.24.0"),
            Dependency("nghttp2_jll"),
        ]
    )

    append!(dependencies, jll_stdlibs[min_julia_version])

    # Build the tarballs, and possibly a `build.jl` as well.
    return name, version, sources, script, platforms, products, dependencies
end

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
module RandomWalks

using Reexport
import Distributions
using Distributions: Distribution, Exponential, Pareto

include("utils.jl")

include("Points.jl")
@reexport using .Points

include("LatticeVars.jl")
@reexport using .LatticeVars

include("Lattices.jl")
@reexport using .Lattices

include("WalksBase.jl")
@reexport using .WalksBase

include("Actors.jl")
@reexport using .Actors

include("Walks.jl")
@reexport using .Walks

include("LatticeWalks.jl")
@reexport using .LatticeWalks

end # module


text:
using DeepUncertainty
using Test
using Flux
using Flux.CUDA
using Flux: cpu, gpu
using DistributionsAD

@testset "Layers" begin
    # MC layers
    include("./layers/mclayers.jl")
    # Batch ensembe layers
    include("./layers/batchensemble.jl")
    # Variational layers
    include("./layers/variational.jl")
    # Variational Batch ensembe layers
    include("./layers/bayesian_be.jl")
    # Spectral layers
    include("./layers/spectrallayers.jl")
end

@testset "Bayesian utils" begin
    include("bayesian.jl")
end

@testset "CUDA" begin
    if CUDA.functional()
        include("cuda/runtests.jl")
    else
        @warn "CUDA unavailable, not testing GPU support"
    end
end


text:
"""
Implements configuration B, i.e. BLG in Bernal stacking with self energy on both the ZZ and
the AC edges. There is a different phase difference between left and right sides of the
weak link (-ϕ/2 for x<0 ϕ/2 for x>0).
Neither smoothness nor leads are implemented. Also disorder and rotation are not considered.
"""
function rectangle_weaklink(p)
    (; Ls, Δ, Ws) = p
    (; model0, field!, modelinter) = modelS(p)
    lat_top, lat_bot = MBSinBLG.latBLG(p, 0) #0 angle rotation

    # LOCAL SELF-ENERGY MODEL
    xmin, xmax = extrema(r->r[1], sitepositions(Quantica.combine(lat_top, lat_bot)))
    ymin, ymax = extrema(r->r[2], sitepositions(Quantica.combine(lat_top, lat_bot)))
    self_region(r) = !(xmin+Ls <= r[1] <= xmax - Ls) || !(ymin+Ws <= r[2] <= ymax - Ws)
    diagphi(phi) = Diagonal(SA[cis(phi), cis(phi), cis(-phi), cis(-phi)])
    sCself! = @onsite!((o, r; phi) -> o +
        self_region(r) * Δ * diagphi(sign(r[1])* phi/4) * σyτy * diagphi(sign(r[1])*phi/4)')

    # HAMILTONIAN BUILD
    h_top = lat_top |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5)
    h_bot = lat_bot |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5)
    ph = Quantica.combine(h_top, h_bot; coupling = modelinter) |>
        parametric(sCself!)
            return ph
end

"""
    spectrumvsphase(philist, p; nev = 16, kw...)
computes a number nev of eigenvalues as a function of the phase difference across the weak
link.
    see: rectangle_weakling(::Params()), Params(), modelS()
"""
function spectrumvsphase(philist, p; nev = 16, kw...)
    ph = rectangle_weaklink(p; kw...)
    elist = zeros(Float64, nev, length(philist))
    for i in 1:length(philist)
        s = spectrum(ph(phi = philist[i]), method = ArpackPackage(sigma = 1e-7, nev = nev))
        elist[:,i] = s.energies
    end
    return philist, elist
end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# yearmonthday is the opposite of totaldays
# taking Rata Die Day # and returning proleptic Gregorian date
@test Dates.yearmonthday(-306) == (0, 2, 29)
@test Dates.yearmonth(-306) == (0, 2)
@test Dates.monthday(-306) == (2, 29)
@test Dates.yearmonthday(-305) == (0, 3, 1)
@test Dates.yearmonth(-305) == (0, 3)
@test Dates.monthday(-305) == (3, 1)
@test Dates.yearmonthday(-2) == (0, 12, 29)
@test Dates.yearmonth(-2) == (0, 12)
@test Dates.monthday(-2) == (12, 29)
@test Dates.yearmonthday(-1) == (0, 12, 30)
@test Dates.yearmonth(-1) == (0, 12)
@test Dates.monthday(-1) == (12, 30)
@test Dates.yearmonthday(0) == (0, 12, 31)
@test Dates.yearmonth(-0) == (0, 12)
@test Dates.monthday(-0) == (12, 31)
@test Dates.yearmonthday(1) == (1, 1, 1)
@test Dates.yearmonth(1) == (1, 1)
@test Dates.monthday(1) == (1, 1)
# year, month, and day return the indivial components
# of yearmonthday, avoiding additional calculations when possible
@test Dates.year(-1) == 0
@test Dates.month(-1) == 12
@test Dates.day(-1) == 30
@test Dates.year(0) == 0
@test Dates.month(0) == 12
@test Dates.day(0) == 31
@test Dates.year(1) == 1
@test Dates.month(1) == 1
@test Dates.day(1) == 1
@test Dates.yearmonthday(730120) == (2000, 1, 1)
@test Dates.year(730120) == 2000
@test Dates.month(730120) == 1
@test Dates.day(730120) == 1

# Test totaldays and yearmonthday from January 1st of "from" to December 31st of "to"
# test_dates(-10000, 10000) takes about 15 seconds
# test_dates(year(typemin(Date)), year(typemax(Date))) is full range
# and would take.......a really long time
let from=0, to=2100, y=0, m=0, d=0
    test_day = Dates.totaldays(from, 1, 1)
    for y in from:to
        for m = 1:12
            for d = 1:Dates.daysinmonth(y, m)
                days = Dates.totaldays(y, m, d)
                @test days == test_day
                @test (y, m, d) == Dates.yearmonthday(days)
                test_day += 1
            end
        end
    end
end

# Test year, month, day, hour, minute
let y=0, m=0, d=0, h=0, mi=0
    for m = 1:12
        for d = 1:Dates.daysinmonth(y, m)
            for h = 0:23
                for mi = 0:59
                    dt = Dates.DateTime(y, m, d, h, mi)
                    @test y == Dates.year(dt)
                    @test m == Dates.month(dt)
                    @test d == Dates.day(dt)
                    @test d == Dates.dayofmonth(dt)
                    @test h == Dates.hour(dt)
                    @test mi == Dates.minute(dt)
                    @test (m, d) == Dates.monthday(dt)
                end
            end
        end
    end
end

# Test second, millisecond
let y=0, m=0, d=0, h=0, mi=0, s=0, ms=0
    for y in [-2013, -1, 0, 1, 2013]
        for m in [1, 6, 12]
            for d in [1, 15, Dates.daysinmonth(y, m)]
                for h in [0, 12, 23]
                    for s = 0:59
                        for ms in [0, 1, 500, 999]
                            dt = Dates.DateTime(y, m, d, h, mi, s, ms)
                            @test y == Dates.year(dt)
                            @test m == Dates.month(dt)
                            @test d == Dates.day(dt)
                            @test h == Dates.hour(dt)
                            @test s == Dates.second(dt)
                            @test ms == Dates.millisecond(dt)
                        end
                    end
                end
            end
        end
    end
end

let from=0, to=2100, y=0, m=0, d=0
    for y in from:to
        for m = 1:12
            for d = 1:Dates.daysinmonth(y, m)
                dt = Dates.Date(y, m, d)
                @test y == Dates.year(dt)
                @test m == Dates.month(dt)
                @test d == Dates.day(dt)
            end
        end
    end
end

# test hour, minute, second
let h=0, mi=0, s=0, ms=0, us=0, ns=0
    for h = (0, 23), mi = (0, 59), s = (0, 59),
        ms in (0, 1, 500, 999), us in (0, 1, 500, 999), ns in (0, 1, 500, 999)
        t = Dates.Time(h, mi, s, ms, us, ns)
        @test h == Dates.hour(t)
        @test mi == Dates.minute(t)
        @test s == Dates.second(t)
        @test ms == Dates.millisecond(t)
        @test us == Dates.microsecond(t)
        @test ns == Dates.nanosecond(t)
    end
end

# week function
# Tests from https://en.wikipedia.org/wiki/ISO_week_date
@test Dates.week(Dates.Date(2005, 1, 1)) == 53
@test Dates.week(Dates.Date(2005, 1, 2)) == 53
@test Dates.week(Dates.Date(2005, 12, 31)) == 52
@test Dates.week(Dates.Date(2007, 1, 1)) == 1
@test Dates.week(Dates.Date(2007, 12, 30)) == 52
@test Dates.week(Dates.Date(2007, 12, 31)) == 1
@test Dates.week(Dates.Date(2008, 1, 1)) == 1
@test Dates.week(Dates.Date(2008, 12, 28)) == 52
@test Dates.week(Dates.Date(2008, 12, 29)) == 1
@test Dates.week(Dates.Date(2008, 12, 30)) == 1
@test Dates.week(Dates.Date(2008, 12, 31)) == 1
@test Dates.week(Dates.Date(2009, 1, 1)) == 1
@test Dates.week(Dates.Date(2009, 12, 31)) == 53
@test Dates.week(Dates.Date(2010, 1, 1)) == 53
@test Dates.week(Dates.Date(2010, 1, 2)) == 53
@test Dates.week(Dates.Date(2010, 1, 2)) == 53
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=1999
dt = Dates.DateTime(1999, 12, 27)
dt1 = Dates.Date(1999, 12, 27)
check = (52, 52, 52, 52, 52, 52, 52, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2000
dt = Dates.DateTime(2000, 12, 25)
dt1 = Dates.Date(2000, 12, 25)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Test from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2030
dt = Dates.DateTime(2030, 12, 23)
dt1 = Dates.Date(2030, 12, 23)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end
# Tests from http://www.epochconverter.com/date-and-time/weeknumbers-by-year.php?year=2004
dt = Dates.DateTime(2004, 12, 20)
dt1 = Dates.Date(2004, 12, 20)
check = (52, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 53, 1, 1, 1, 1, 1, 1, 1)
for i = 1:21
    @test Dates.week(dt) == check[i]
    @test Dates.week(dt1) == check[i]
    dt = dt + Dates.Day(1)
    dt1 = dt1 + Dates.Day(1)
end

# Vectorized accessors
a = Dates.Date(2014, 1, 1)
dr = [a, a, a, a, a, a, a, a, a, a]
@test Dates.year.(dr) == repmat([2014], 10)
@test Dates.month.(dr) == repmat([1], 10)
@test Dates.day.(dr) == repmat([1], 10)

a = Dates.DateTime(2014, 1, 1)
dr = [a, a, a, a, a, a, a, a, a, a]
@test Dates.year.(dr) == repmat([2014], 10)
@test Dates.month.(dr) == repmat([1], 10)
@test Dates.day.(dr) == repmat([1], 10)
@test Dates.hour.(dr) == repmat([0], 10)
@test Dates.minute.(dr) == repmat([0], 10)
@test Dates.second.(dr) == repmat([0], 10)
@test Dates.millisecond.(dr) == repmat([0], 10)

b = Dates.Time(1, 2, 3, 4, 5, 6)
tr = [b, b, b, b, b, b, b, b, b, b]
@test Dates.hour.(tr) == repmat([1], 10)
@test Dates.minute.(tr) == repmat([2], 10)
@test Dates.second.(tr) == repmat([3], 10)
@test Dates.millisecond.(tr) == repmat([4], 10)
@test Dates.microsecond.(tr) == repmat([5], 10)
@test Dates.nanosecond.(tr) == repmat([6], 10)

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
import ParticleInCell.F90

@testset "Landau Damping Julia" begin

    nx = 128
    ny = 16
    alpha = 0.1
    kx = 0.5
    ky = 0.0
    dimx = 2pi / kx
    dimy = 1.0

    mesh1 = TwoDGrid(dimx, nx, dimy, ny)
    mesh2 = TwoDGrid(dimx, nx, dimy, ny)

    ex = zeros(nx+1, ny+1)
    ey = zeros(nx+1, ny+1)
    bz = zeros(nx+1, ny+1)
    jx = zeros(nx+1, ny+1)
    jy = zeros(nx+1, ny+1)

    nbpart = 100 * nx * ny

    group1 = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)
    group2 = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)

    alpha = 0.1
    kx = 0.5

    sampler = LandauDamping(alpha, kx)

    ParticleInCell.sample!(group1, mesh1, sampler)
    ParticleInCell.sample!(group2, mesh2, sampler)

    fdtd1 = FDTD(mesh1)
    fdtd2 = FDTD(mesh2)

    time = 0
    dt = 0.01

    for i = 1:nx, j = 1:ny+1
        fdtd1.ex[i, j] = alpha / kx * sin(kx * mesh1.x[i])
        fdtd2.ex[i, j] = alpha / kx * sin(kx * mesh2.x[i])
    end

    kernel = CloudInCell()

    particles1 = group1.array
    particles2 = group2.array

    for istep = 1:1

        istep > 1 && faraday!(fdtd1, mesh1, 0.5dt)
        istep > 1 && faraday!(fdtd2, mesh2, 0.5dt)

        update_fields!(ex, ey, bz, mesh1, fdtd1)
        update_fields!(ex, ey, bz, mesh2, fdtd2)

        push_v!(group1, kernel, mesh1, ex, ey, bz,  dt)
        F90.push_v!(group2, kernel, mesh2, ex, ey, bz, dt)
        @test all(particles1 .== particles2)

        push_x!(group1, mesh1, 0.5dt)
        F90.push_x!(group2, mesh2, 0.5dt)
        @test all(particles1 .== particles2)

        compute_current!(jx, jy, mesh1, kernel, group1)
        F90.compute_current!(jx, jy, mesh2, kernel, group2)

        push_x!(group1, mesh1, 0.5dt)
        F90.push_x!(group2, mesh2, 0.5dt)
        @test all(particles1 .== particles2)

        faraday!(fdtd1, mesh1, 0.5dt)
        faraday!(fdtd2, mesh2, 0.5dt)

        ampere_maxwell!(fdtd1, mesh1, jx, jy, dt)
        ampere_maxwell!(fdtd2, mesh2, jx, jy, dt)

        time = time + dt

    end

    @test true

end


text:
include("../btg.jl")
include("../datasets/load_abalone.jl")
i = 3 #LOOCV at this point

#ind = 350:370
ind = 1050:1097
#posx = 1:3 #
posx = 1:7
posc = 1:1
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean
Fx = data[ind, posc]
y = float(target[ind])
#pind = 10:10 #prediction index
pind = i:i
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
#rangeθ = [100.0 200]
rangeθ = [10.0 2000]
rangeλ = [-1.0 1] #we will always used 1 range scale for lambda
btg1 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
#θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node
#λ1 = btg1.nodesWeightsλ.nodes[3]
##################################################################################################
PyPlot.close("all") #close existing windows

function perturb(tup::Tuple{Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 1}})
    x, Fx, y = tup
    x = x + randn(size(x, 1), size(x, 2)) * 0.1
    Fx = Fx + randn(size(Fx, 1), size(Fx, 2)) * 0.1
    y = y + randn(length(y)) * 0.1
    return (x, Fx, y)
end

(pdf, cdf, dpdf) = solve(btg1); #initialize training_buffer_dicts, solve once so can use fast techiques to extrapolate submatrix determinants, etc.
#pdf(x0, Fx0, .4);
(pdf, cdf, dpdf) = solve(btg1; validate = i) #do LOOCV this time around
a1 = y-> pdf(x0, Fx0, y)

"""
leave-one-out-training-data
"""
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end

#btg 2, trained on subset of trainingData
(trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, i)
btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ)
(pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2)
a2 = y -> pdf_minus_i(x_i, Fx_i, y)

#@assert equals(btg2.θλbuffer_dict, btg1.validation_θλ_buffer_dict) #do the computed buffers match up, only match if quadtype is Gaussian, regular grid, not monte carlo

if false #plot
figure(5)
(x1, y1) = plt_data(a1, .01, 1.5, 150)
(x2, y2) = plt_data(a2, .01, 1.5, 150)
PyPlot.plot(x1, y1, color = "blue")
PyPlot.plot(x2, y2, color = "green")
end

if false
    #function calls to initialize debug_log
    a1(.5);
    a2(.5);

    h1 = sort(btg1.debug_log, by = first, rev = true)[1:5]
    h2 = sort(btg2.debug_log, by = first, rev = true)[1:5]

    f1 = sort(btg1.debug_log, by = last, rev = true)[1:5]
    f2 = sort(btg2.debug_log, by = last, rev = true)[1:5]

end

i=3
(pdf, cdf, dpdf) = solve(btg1; validate = i) #do LOOCV this time around
a3 = y-> pdf(x0, Fx0, y) #doesnt matter what x0 and Fx0 are

#btg 2, trained on subset of trainingData
(trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, i)
btg3 = btg(trainingdata_minus_i, rangeθ, rangeλ)
(pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg3)
a4 = y -> pdf_minus_i(x_i, Fx_i, y)

if false
    #function calls to initialize debug_log
    a3(.5);
    a4(.5);
    (histx, histy) = plt_data(a3, .01, 1.2, 100);
    (histx2, histy2) = plt_data(a4, .01, 1.2, 100);

    figure(27)
    PyPlot.plot(histx, histy)
    Pyplot.plot(histx2, histy2)

    b1 = sort(btg1.debug_log, by = first, rev = true)[1:5]
    b2 = sort(btg3.debug_log, by = first, rev = true)[1:5]

    n1 = sort(btg1.debug_log, by = last, rev = true)[1:5]
    n2 = sort(btg3.debug_log, by = last, rev = true)[1:5]

end

###
###
### plotting stuff
###
###

m = 6; n = 8
plt, axs = PyPlot.subplots(m, n)
#figure(1)
for j = 1:m*n
    (pdf, cdf, dpdf) = solve(btg1, validate = j);
    #a = y0 -> dpdf(x0, Fx0, y0);
    b1 = y0 -> pdf(x0, Fx0, y0);
    #c1 = y0 -> cdf(x0, Fx0, y0);
    (x, y) = plt_data(b1, .01, 1.2, 100)

    (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, j)
    btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
    (pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2)
    b2 = y -> pdf_minus_i(x_i, Fx_i, y)
    (x1, y1) = plt_data(b2, .01, 1.2, 100)

    ind1 = Int64(ceil(j/n))
    ind2 = Int64(j - n*(floor((j-.1)/n)))

    axs[ind1, ind2].plot(x, y, color = "red", linewidth = 2.0, linestyle = ":")
    axs[ind1, ind2].plot(x1, y1, color = "blue", linewidth = 1.0, linestyle = "--")
    axs[ind1, ind2].axvline(x =  getLabel(btg1.trainingData)[j])
    #PyPlot.plot(x, y)
    #PyPlot.plot(x1, y1)
    #PyPlot.axvline(x =  getLabel(btg1.trainingData)[j])
end
for ax in axs
    ax.set(xlabel="x-label", ylabel="y-label")
end
# Hide x labels and tick labels for top plots and y ticks for right plots.
for ax in axs
    ax.label_outer()
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    WebMercatorfromLLA(datum=wgs84)

Convert from LLA to Web Mercator / Pseudo Mercator, following the convention of
proj, which uses a scaling factor of the semi-major axis of the ellipsoid
https://proj.org/operations/projections/webmerc.html.

!!! warning
    For web mapping applications this projection is ubiquitous due to its
    simplicity of implementation, but this simplicity gives rise to poor
    mathematical properties: it doesn't preserve angles away from the equator.
    Other projections should be preferred if possible, and especially when
    measurement is important. See
    https://earth-info.nga.mil/GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf
    for an extended discussion.
"""
struct WebMercatorfromLLA <: Transformation
    el::Ellipsoid
end

WebMercatorfromLLA(d::Datum=wgs84) = WebMercatorfromLLA(ellipsoid(d))

"""
    LLAfromWebMercator

Inverse of WebMercatorfromLLA — see the docs for that transformation.
"""
struct LLAfromWebMercator <: Transformation
    el::Ellipsoid
end

LLAfromWebMercator(d::Datum=wgs84) = LLAfromWebMercator(ellipsoid(d))

function (trans::WebMercatorfromLLA)(lla::LLA)
    xy = trans(LatLon(lla.lat, lla.lon))
    SA[xy[1], xy[2], lla.alt]
end

function (trans::WebMercatorfromLLA)(ll::LatLon)
    lat_lim = 85.06 # according to https://epsg.io/3857
    if abs(ll.lat) > lat_lim
        throw(ArgumentError("Exceeded Web Mercator latitude bounds"))
    end
    x, y = web_mercator_forward(ll.lat, ll.lon, trans.el.a)
    SA[x, y]
end

function (trans::LLAfromWebMercator)(point::AbstractVector)
    x = point[1]
    y = point[2]
    lat, lon = web_mercator_reverse(x, y, trans.el.a)
    if length(point) == 2
        LatLon(lat, lon)
    elseif length(point) == 3
        LLA(lat, lon, point[3])
    else
        throw(ArgumentError("Expected input vector of length 2 or 3"))
    end
end

Base.inv(trans::WebMercatorfromLLA) = LLAfromWebMercator(trans.el)
Base.inv(trans::LLAfromWebMercator) = WebMercatorfromLLA(trans.el)

# Web / Psuedo Mercator. Following the scaling convention of proj, the scaling
# will be set to the ellipsoid semi-major axis.
#   https://proj.org/operations/projections/webmerc.html
function web_mercator_forward(lat, lon, scaling)
    x = scaling * deg2rad(lon)
    y = scaling * log(tand((90 + lat)/2))
    (x,y)
end

function web_mercator_reverse(x, y, scaling)
    lon = rad2deg(x / scaling)
    lat = 2 * atand(exp(y / scaling)) - 90
    (lat,lon)
end


text:
using Base.Test
using MathProgBase

function linprogtest(solver=MathProgBase.defaultLPsolver; objtol = 1e-7, primaltol = 1e-6)
    println("Testing linprog and subfunctions with solver ", string(typeof(solver)))
    # min -x
    # s.t. 2x + y <= 1.5
    # x,y >= 0
    # solution is (0.75,0) with objval -0.75

    # test buildlp and solvelp
    m = buildlp([-1,0],[2 1],'<',1.5,solver)
    sol = solvelp(m)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval -0.75 objtol
    @test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

    # test linprog
    sol = linprog([-1,0],[2 1],'<',1.5,solver)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval -0.75 objtol
    @test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

    # test buildlp and solvelp
    m = buildlp([-1,0],[2 1],'<',1.5,solver)
    sol = solvelp(m)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval -0.75 objtol
    @test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

    # test linprog
    sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval -0.75 objtol
    @test_approx_eq_eps norm(sol.sol - [0.75,0.0]) 0 primaltol

    # test infeasible problem:
    # min x
    # s.t. 2x+y <= -1
    # x,y >= 0

    # test buildlp and solvelp
    m = buildlp([1,0],[2 1],'<',-1,solver)
    sol = solvelp(m)
    @test sol.status == :Infeasible

    r = sol.attrs[:infeasibilityray][1]
    @test_approx_eq r/abs(r) -1.0

    # test linprog
    sol = linprog([1,0],[2 1],'<',-1,solver)
    @test sol.status == :Infeasible

    r = sol.attrs[:infeasibilityray][1]
    @test_approx_eq r/abs(r) -1.0

    # test unbounded problem:
    # min -x-y
    # s.t. -x+2y <= 0
    # x,y >= 0

    # test buildlp and solvelp
    m = buildlp([-1,-1],[-1 2],'<',[0],solver)
    sol = solvelp(m)
    @test sol.status == :Unbounded

    # test linprog
    sol = linprog([-1,-1],[-1 2],'<',[0],solver)
    @test sol.status == :Unbounded

    # unbounded problem with unique ray:
    # min -x-y
    # s.t. x-y == 0
    # x,y >= 0

    # test buildlp and solvelp
    m = buildlp([-1,-1],[1 -1],'=',0,solver)
    sol = solvelp(m)
    @test sol.status == :Unbounded
    @test sol.attrs[:unboundedray][1] > 1e-7
    @test_approx_eq sol.attrs[:unboundedray][1] sol.attrs[:unboundedray][2]

    # test linprog
    sol = linprog([-1,-1],[1 -1],'=',0,solver)
    @test sol.status == :Unbounded
    @test sol.attrs[:unboundedray][1] > 1e-7
    @test_approx_eq sol.attrs[:unboundedray][1] sol.attrs[:unboundedray][2]

    println("Passed")
end


text:
using LuaCall
using Test

f(x, y) = (x^2 * sin(y)) / (1 - x)
lua"""
    function f(x, y)
        return (x^2 * math.sin(y)) / (1 - x)
    end
"""
@test luacall(:f, Float64, 2.0, 3.0) == f(2.0, 3.0)


text:
module FashionMNIST

using ..MNIST: gzopen, imageheader, rawimage, labelheader, rawlabel
using ..Data: download_and_verify

const dir = joinpath(@__DIR__, "../../deps/fashion-mnist")

function load()
  mkpath(dir)
  cd(dir) do
    for (file, hash) in [("train-images-idx3-ubyte", "3aede38d61863908ad78613f6a32ed271626dd12800ba2636569512369268a84"),
                         ("train-labels-idx1-ubyte", "a04f17134ac03560a47e3764e11b92fc97de4d1bfaf8ba1a3aa29af54cc90845"),
                         ("t10k-images-idx3-ubyte" , "346e55b948d973a97e58d2351dde16a484bd415d4595297633bb08f03db6a073"),
                         ("t10k-labels-idx1-ubyte" , "67da17c76eaffca5446c3361aaab5c3cd6d1c2608764d35dfb1850b086bf8dd5")]
      isfile(file) && continue
      @info "Downloading Fashion-MNIST dataset"
      download_and_verify("http://fashion-mnist.s3-website.eu-central-1.amazonaws.com/$file.gz", "$file.gz", hash)
      open(file, "w") do io
        write(io, gzopen(read, "$file.gz"))
      end
    end
  end
end

const TRAINIMAGES = joinpath(dir, "train-images-idx3-ubyte")
const TRAINLABELS = joinpath(dir, "train-labels-idx1-ubyte")
const TESTIMAGES = joinpath(dir, "t10k-images-idx3-ubyte")
const TESTLABELS = joinpath(dir, "t10k-labels-idx1-ubyte")

"""
    images()
    images(:test)

Load the Fashion-MNIST images.

Each image is a 28×28 array of `Gray` colour values
(see [Colors.jl](https://github.com/JuliaGraphics/Colors.jl)).

Return the 60,000 training images by default; pass `:test` to retrieve the
10,000 test images.
"""
function images(set = :train)
  load()
  io = IOBuffer(read(set == :train ? TRAINIMAGES : TESTIMAGES))
  _, N, nrows, ncols = imageheader(io)
  [rawimage(io) for _ in 1:N]
end

"""
    labels()
    labels(:test)

Load the labels corresponding to each of the images returned from [`images()`](@ref).
Each label is a number from 0-9.

Return the 60,000 training labels by default; pass `:test` to retrieve the
10,000 test labels.
"""
function labels(set = :train)
  load()
  io = IOBuffer(read(set == :train ? TRAINLABELS : TESTLABELS))
  _, N = labelheader(io)
  [rawlabel(io) for _ = 1:N]
end

end


text:
@inline _serror() = throw(ArgumentError("Scale factor K has to be strictly positive."))

for KIND in (:MarginLoss, :DistanceLoss, :SupervisedLoss)
    SCALEDKIND = Symbol(:Scaled, KIND)
    @eval begin
        @doc """
            $($(string(SCALEDKIND))){L,K} <: $($(string(KIND)))

        Can an be used to represent a `K` times scaled version of a
        given type of loss `L`, which must be a subtype of
        `$($(string(KIND)))`.
        For example: To create a typealias for a `1.5` times scaled
        version of some loss `My$($(string(KIND)))`, type:

        ```julia
        const MyScaled$($(string(KIND))) = LossFunctions.$($(string(SCALEDKIND))){My$($(string(KIND))),1.5}
        ```

        This new loss-type can then be instantiated in the same
        manner and with the same parameters as the original unscaled
        loss-type.

        In contrast, in order to only create a `K` times scaled
        instance of some specific loss you can use
        `scaled(my$($(lowercase(string(KIND)))), Val(K))`.
        See `?scaled` for more information.
        """
        struct ($SCALEDKIND){L<:$KIND,K} <: $KIND
            loss::L
            ($SCALEDKIND){L,K}(loss::L) where {L<:$KIND, K} = new{L,K}(loss)
        end

        @generated function (::Type{($SCALEDKIND){L,K}})(args...) where {L<:$KIND, K}
            typeof(K) <: Number || _serror()
            :(($($SCALEDKIND)){L,K}(L(args...)))
        end
        ($SCALEDKIND)(loss::T, ::Val{K}) where {T,K} = ($SCALEDKIND){T,K}(loss)
        ($SCALEDKIND)(loss::T, k::Number) where {T} = ($SCALEDKIND)(loss, Val(k))
        (*)(::Val{K}, loss::$KIND) where {K} = ($SCALEDKIND)(loss, Val(K))
        (*)(k::Number, loss::$KIND) = ($SCALEDKIND)(loss, Val(k))
        scaled(loss::T, ::Val{K}) where {T<:$KIND,K} = ($SCALEDKIND)(loss, Val(K))

        @generated ($SCALEDKIND)(s::$SCALEDKIND{T,K1}, ::Val{K2}) where {T,K1,K2} = :(($($SCALEDKIND))(s.loss, Val($(K1*K2))))

    end
end

"""
    Union{ScaledSupervisedLoss, ScaledMarginLoss, ScaledDistanceLoss}

Mainly intended for dispatch. Look at `?ScaledMarginLoss`,
or `ScaledDistanceLoss` for more information.

use `scaled` to create a scaled version of some loss.
"""
const ScaledLoss{T,K} = Union{ScaledSupervisedLoss{T,K}, ScaledMarginLoss{T,K}, ScaledDistanceLoss{T,K}}

"""
    scaled(loss::SupervisedLoss, K)

Returns a version of `loss` that is uniformly scaled by `K`.
This function dispatches on the type of `loss` in order to choose
the appropriate type of scaled loss that will be used as the decorator.
For example, if `typeof(loss) <: DistanceLoss` then the given `loss`
will be boxed into a `ScaledDistanceLoss`.

Note: If `typeof(K) <: Number`, then this method will poison the
type-inference of the calling scope. This is because `K` will be
promoted to a type parameter. For a typestable version use the
following signature: `scaled(loss, Val(K))`
"""
scaled(l::Loss, k::Number) = scaled(l, Val(k))

for fun in (:value, :deriv, :deriv2)
    @eval @fastmath ($fun)(l::ScaledLoss{T,K}, num::Number) where {T,K} = K * ($fun)(l.loss, num)
    @eval @fastmath ($fun)(l::ScaledLoss{T,K}, target::Number, output::Number) where {T,K} = K * ($fun)(l.loss, target, output)
end

for prop in [:isminimizable, :isdifferentiable,
             :istwicedifferentiable,
             :isconvex, :isstrictlyconvex,
             :isstronglyconvex, :isnemitski,
             :isunivfishercons, :isfishercons,
             :islipschitzcont, :islocallylipschitzcont,
             :isclipable, :ismarginbased, :isclasscalibrated,
             :isdistancebased, :issymmetric]
    @eval ($prop)(l::ScaledLoss) = ($prop)(l.loss)
end

for prop_param in (:isdifferentiable, :istwicedifferentiable)
    @eval ($prop_param)(l::ScaledLoss, at) = ($prop_param)(l.loss, at)
end


text:
using BytePairEncoding: isolate_gloss
@testset "glossary" begin
    @testset "isolate_gloss" begin
        gloss = "like"
        isolike(x) = isolate_gloss(x, gloss)
        @test isolike("") == [""]
        @test isolike("word") == ["word"]
        @test isolike("like") == ["like"]
        @test isolike("likeword") == ["like", "word"]
        @test isolike("wordlikeword") == ["word", "like", "word"]
        @test isolike("likelike") == ["like", "like"]
        @test isolike("wordlikewordlike") == ["word", "like", "word", "like"]
    end

    @testset "multiple gloss" begin
        glossaries = ["like", "Manuel", "USA"]
        @test isolate_gloss("wordlikeUSAwordManuelManuelwordUSA" , glossaries) == ["word", "like", "USA", "word", "Manuel", "Manuel", "word", "USA"]
    end

    @testset "regex gloss" begin
        glossaries = [r"<country>\w*</country>", r"<name>\w*</name>", r"\d+"]
        @test isolate_gloss("wordlike<country>USA</country>word10001word<name>Manuel</name>word<country>USA</country>", glossaries) == ["wordlike", "<country>USA</country>", "word", "10001", "word", "<name>Manuel</name>", "word", "<country>USA</country>"]
    end

    @testset "bpe with gloss" begin
        glossaries = ["like", "Manuel", "USA"]

        #mocking bpe
        bpefile = joinpath(dirname(@__FILE__), "data/bpe.out")
        bpe = Bpe(bpefile; sepsym="@@", endsym="/", glossaries = glossaries)
        function BytePairEncoding.uncache_bpe(bpe::GenericBPE, x::String)
            if BytePairEncoding.isgloss(bpe.glossaries, x)
                return [x*bpe.endsym]
            else
                l = length(x)
                return [x[1:div(l,2)]*bpe.sepsym, x[div(l,2)+1:end]*bpe.endsym]
            end
        end

        @test join(segment(bpe, "wordlikeword likeManuelword"), " ") == "wo@@ rd/ like/ wo@@ rd/ like/ Manuel/ wo@@ rd/"

    end
end


text:
# doubly periodic dipole
# https://dx.doi.org/10.1103/PhysRevLett.112.145301
# should be eigenstate of Hamiltonian evolution

using Test, Plots, Revise, VortexDistributions

# test dipole
xp = 20
yp = 3
vp = PointVortex(xp,yp,1)

xn = 1.1
yn = -3
vn = PointVortex(xn,yn,-1)

Thetad(.12,.11,xp,yp,xn,yn)

x = LinRange(-π,π,300)
y = x

#TODO add missing method for Thetad(x,y,dipole::Array{PointVortex})
# test double periodicity
testpsi = exp.(im*Thetad(x,y,xp,yp,xn,yn))
@test testpsi[end,:] ≈ testpsi[1,:]
@test testpsi[:,1] ≈ testpsi[:,end]
heatmap(x,y,angle.(testpsi),transpose=true)

# test domain
x = LinRange(-3π,3π,300)
y = LinRange(-4π,4π,300)

# test dipole
xp = 2.1*pi
yp = 2*pi

xn = -2*pi
yn = -2.4*pi

testphase3 = Thetad(x,y,xp,yp,xn,yn)
heatmap(x,y,testphase3,transpose=true)

psivort = Torus(exp.(im*testphase3),x,y)
vortices = findvortices(psivort)

vraw = rawData(vortices)

testpsi = exp.(im*testphase3)
@test testpsi[:,1] ≈ testpsi[:,end]
@test testpsi[1,:] ≈ testpsi[end,:]

#--- add new method for making periodic dipoles
using Test, Plots, Revise, VortexDistributions

x = LinRange(-50,50,300)
y = x

# test dipole
xp = 20
yp = 3
vp = PointVortex(xp,yp,1)

xn = 1.1
yn = -3
vn = PointVortex(xn,yn,-1)

dip = ScalarVortex([vp;vn])

function periodic_dipole!(psi::F,dip::Array{ScalarVortex{T},1}) where {T <: CoreShape, F<:Field}
    @assert length(dip) == 2
    @assert dip[1].vort.qv + dip[2].vort.qv == 0
    @assert hypot(dip[1].vort.xv-dip[2].vort.xv,dip[1].vort.yv-dip[2].vort.yv) >= 2*pi
    @unpack ψ,x,y = psi
    (dip[1].vort.qv > 0) ? (jp = 1;jn = 2) : (jp = 2;jn = 1)
    vp = rawData(dip[jp].vort)[1:2]
    vn = rawData(dip[jn].vort)[1:2]
    @. ψ *= abs(dip[jn](x,y')*dip[jp](x,y'))
    ψ .*= exp.(im*Thetad(x,y,vp...,vn...))
    @pack! psi = ψ
end

psivort = Torus(one.(x.*y'),x,y)

heatmap(abs2.(psivort.ψ))

periodic_dipole!(psivort,dip)

heatmap(x,y,abs2.(psivort.ψ))

heatmap(x,y,angle.(psivort.ψ))

#check that phase is periodic (note we don't enforce periodic density)

testpsi = angle.(psivort.ψ)
@test testpsi[:,1] ≈ testpsi[:,end]
@test testpsi[1,:] ≈ testpsi[end,:]


text:
#= modular_kpp.jl

Here we implement a 'modular' KPP model, with three interchangable components:

1. A model for mixing depth, h
2. A model for the local diffusivity, K
3. A model for the nonlocal flux term, A

Note below the following acronyms:

* LMD94: Large, McWilliams, and Doney (1994) "Oceanic vertical mixing..."
* RH18: Riechl and Hallberg (2018) "ePBL"
* LMD07: Siebsma, Soares and Teixiera (2007) "An eddy diffusivity mass flux..."

For mixing depth models we have

1. The diagnostic bulk Richarson number scheme proposed by LMD94
2. The diagnostic TKE budget-based scheme proposed by RH18

For K-profile models we have

1. The K-profile proposed by LMD94
2. The K-profile proposed by Holtslag (1998), and described in SST07

For nonlocal flux models we have

1. The countergradient flux model proposed by LMD94
2. The mass flux model proposed by SST07

=#

module ModularKPP

export
    LMDMixingDepth,
    LMDCounterGradientFlux,
    LMDDiffusivity,
    ROMSMixingDepth,
    HoltslagDiffusivity,
    DiagnosticPlumeModel

using
    OceanTurb,
    LinearAlgebra

using OceanTurb: maxsqrt, ∂z⁺ # upwards-biased difference

import OceanTurb.KPP: 𝒲_unstable, 𝒲_stable, u★, w★, d,
                      isunstable, isforced, unresolved_kinetic_energy,
                      ∂B∂z

abstract type AbstractModularKPPModel{K, H, N, TS, G, T} <: AbstractModel{TS, G, T} end

const nsol = 4
@solution U V T S

"""
    Model{S, G, T, U, B, F} <: AbstractModel{S, G, T}

Struct for KPP models.
"""
mutable struct Model{KP, NP, HP, SP, SO, BC, ST, TS, G, T, F} <: AbstractModularKPPModel{KP, NP, HP, TS, G, T}
           clock :: Clock{T}
            grid :: G
     timestepper :: TS
        solution :: SO
             bcs :: BC
     diffusivity :: KP
    nonlocalflux :: NP
     mixingdepth :: HP
        kprofile :: SP
       constants :: Constants{T}
           state :: ST
         forcing :: F
end

include("state.jl")
include("model_boundary_conditions.jl")
include("forcing.jl")
include("mixing_depth.jl")
include("diffusivity_models.jl")
include("shape_functions.jl")
include("nonlocal_flux.jl")

"""
    Model(; kwargs...)

Construct a KPP Model.
"""
function Model(;
            grid = UniformGrid(3, 1.0),
       constants = Constants(),
     diffusivity = LMDDiffusivity(),
    nonlocalflux = LMDCounterGradientFlux(),
     mixingdepth = LMDMixingDepth(),
        kprofile = StandardCubicPolynomial(),
         stepper = :BackwardEuler,
             bcs = ModelBoundaryConditions(eltype(grid)),
         forcing = Forcing()
    )

     K = Accessory{Function}(KU, KV, KT, KS)
     R = Accessory{Function}(RU, RV, RT, RS)
     A = Accessory{Function}(AU, AV, AT, AS)
    eq = Equation(K=K, R=R, A=A, update=update_state!)

       state = State(diffusivity, nonlocalflux, mixingdepth, grid)
    solution = Solution((CellField(grid) for i=1:nsol)...)
         lhs = OceanTurb.build_lhs(solution)

    timestepper = Timestepper(stepper, eq, solution, lhs)

    return Model(Clock(), grid, timestepper, solution, bcs,
                 diffusivity, nonlocalflux, mixingdepth, kprofile, constants, state,
                 forcing)
end

#
# Equation specification
#

RU(m, i) =   m.constants.f * m.solution.V[i] + m.forcing.U(m, i)
RV(m, i) = - m.constants.f * m.solution.U[i] + m.forcing.V(m, i)

const LMDModel = AbstractModularKPPModel{<:LMDDiffusivity}
const HoltslagModel = AbstractModularKPPModel{<:HoltslagDiffusivity}

# K_{U,V,T,S} is calculated at face points
KU(m::LMDModel, i) = K_KPP(m.state.h, 𝒲_LMD_U(m, i), d(m, i), m.kprofile) + m.diffusivity.KU₀
KT(m::LMDModel, i) = K_KPP(m.state.h, 𝒲_LMD_T(m, i), d(m, i), m.kprofile) + m.diffusivity.KT₀
KS(m::LMDModel, i) = K_KPP(m.state.h, 𝒲_LMD_S(m, i), d(m, i), m.kprofile) + m.diffusivity.KS₀

KU(m::HoltslagModel, i) = K_KPP(m.state.h, 𝒲_Holtslag(m, i), d(m, i), m.kprofile) + m.diffusivity.KU₀
KT(m::HoltslagModel, i) = K_KPP(m.state.h, 𝒲_Holtslag(m, i), d(m, i), m.kprofile) + m.diffusivity.KT₀
KS(m::HoltslagModel, i) = K_KPP(m.state.h, 𝒲_Holtslag(m, i), d(m, i), m.kprofile) + m.diffusivity.KS₀

const KV = KU

RT(m, i) = - ∂z_explicit_nonlocal_flux_T(m, i) + m.forcing.T(m, i)
RS(m, i) = - ∂z_explicit_nonlocal_flux_S(m, i) + m.forcing.S(m, i)

AU(m, i) = 0
AV(m, i) = 0

# Using explicit advection for now:
AT(m, i) = 0 #- mass_flux(m, i)
AS(m, i) = 0 #- mass_flux(m, i)

end # module


text:
"""
    LogWeight(β,α,s::Space)

represents a function on `-1..1` weighted by `log((1+x)^β*(1-x)^α)`.
For other domains, the weight is inferred by mapping to `-1..1`.
"""
struct LogWeight{S,DD,RR,T} <: WeightSpace{S,DD,RR}
    β::T
    α::T
    space::S
end

LogWeight(β,α,space) =
    LogWeight{typeof(space),domaintype(space),rangetype(space),promote_type(typeof(β),typeof(α))}(β,α,space)

spacescompatible(A::LogWeight,B::LogWeight)=A.β==B.β && A.α == B.α && spacescompatible(A.space,B.space)
canonicalspace(A::LogWeight)=A

logweight(β,α,x)=log((1+x)^β*(1-x)^α)
weight(sp::LogWeight,x)=logweight(sp.β,sp.α,tocanonical(sp,x))

setdomain(sp::LogWeight,d::Domain)=LogWeight(sp.β,sp.α,setdomain(sp.space,d))

function coefficients(f::AbstractVector,sp1::LogWeight,sp2::LogWeight)
    β,α=sp1.β,sp1.α
    c,d=sp2.β,sp2.α

    if isapprox(c,β) && isapprox(d,α)
        coefficients(f,sp1.space,sp2.space)
    else
        (Conversion(sp1,sp2)*f)
    end
end

for (OPrule,OP) in ((:maxspace_rule,:maxspace),(:union_rule,:union))
    @eval begin
        function $OPrule(A::LogWeight,B::LogWeight)
            if isapprox(A.β,B.β) && isapprox(A.α,B.α)
                LogWeight(A.β,A.α,$OP(A.space,B.space))
            else
                NoSpace()
            end
        end
        # there are no other comatible spaces
        # this is mostly to overcome typing bug
        # in spacestes
        $OPrule(::LogWeight,::JacobiWeight)=NoSpace()
        $OPrule(::JacobiWeight,::LogWeight)=NoSpace()
    end
end

# Same as JacobiWeight

# avoid redundency
function Multiplication(f::Fun{JacobiWeight{SS,DD,RR,TT},T}, S::LogWeight{LWS,DD,RR}) where {SS,LWS,DD<:IntervalOrSegment,RR,T,TT}
    M=Multiplication(Fun(space(f).space,f.coefficients),S)
    rsp=JacobiWeight(space(f).β,space(f).α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun{<:LogWeight},S::JacobiWeight{SS,DD}) where {SS,DD<:IntervalOrSegment}
    M=Multiplication(f,S.space)
    rsp=JacobiWeight(S.β,S.α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun,S::LogWeight)
    M=Multiplication(f,S.space)
    rsp=LogWeight(S.β,S.α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun{<:LogWeight},S::LogWeight)
    M=Multiplication(f,S.space)
    rsp=LogWeight(S.β,S.α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun{<:LogWeight},S::Space)
    M=Multiplication(Fun(space(f).space,f.coefficients),S)
    rsp=LogWeight(space(f).β,space(f).α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end


text:
using ABCDeGraphGenerator
using Random

@info "Usage: julia com_sampler.jl filename τ₂ c_min c_max n max_iter [seed]"
@info "Example: julia com_sampler.jl community_sizes.dat 2 50 1000 10000 1000 42"

filename = ARGS[1]
τ₂ = parse(Float64, ARGS[2])
c_min = parse(Int, ARGS[3])
c_max = parse(Int, ARGS[4])
n = parse(Int, ARGS[5])
max_iter = parse(Int, ARGS[6])
length(ARGS) == 7 && Random.seed!(parse(Int, ARGS[7]))

@info "Expected value of community size: $(ABCDeGraphGenerator.get_ev(τ₂, c_min, c_max))"

coms = ABCDeGraphGenerator.sample_communities(τ₂, c_min, c_max, n, max_iter)

open(io -> foreach(d -> println(io, d), coms), filename, "w")


text:
export AU

const AU = 1.49597870700e8 # km


text:
# Define some points and data in 2D
arrayPoints = permutedims([0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0], (2,1))
adjointPoints = [0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0]'
data = [0.0; 0.5; 0.5; 1.0]

@testset "Shepard" begin

    @testset "Constructors" for idw in (:Shepard, )
        @eval @test $idw(2) == $idw()
    end

    @testset "Evaluation" for points in (arrayPoints, adjointPoints)

        itp = interpolate(Shepard(), points, data)

        # Check that we get back the original data at the sample points and that we get close
        # when evaluating near the sampling points
        ev = evaluate(itp, points)
        @test ev ≈ data
        ev = evaluate(itp, points .+ 0.001*randn(size(points)))
        @test all(isapprox.(data, ev, atol = 1e-2))
    end
end


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

@doc raw"""NetworkPolicy describes what network traffic is allowed for a set of Pods

    IoK8sApiNetworkingV1NetworkPolicy(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - spec::IoK8sApiNetworkingV1NetworkPolicySpec : Specification of the desired behavior for this NetworkPolicy.
"""
mutable struct IoK8sApiNetworkingV1NetworkPolicy <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1NetworkPolicySpec } # spec name: spec

    function IoK8sApiNetworkingV1NetworkPolicy(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1NetworkPolicy, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiNetworkingV1NetworkPolicy, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiNetworkingV1NetworkPolicy, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiNetworkingV1NetworkPolicy, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        o
    end
end # type IoK8sApiNetworkingV1NetworkPolicy

const _property_map_IoK8sApiNetworkingV1NetworkPolicy = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"))
const _property_types_IoK8sApiNetworkingV1NetworkPolicy = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiNetworkingV1NetworkPolicySpec")
Base.propertynames(::Type{ IoK8sApiNetworkingV1NetworkPolicy }) = collect(keys(_property_map_IoK8sApiNetworkingV1NetworkPolicy))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1NetworkPolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1NetworkPolicy[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1NetworkPolicy }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1NetworkPolicy[property_name]

function check_required(o::IoK8sApiNetworkingV1NetworkPolicy)
    true
end

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


text:
using Test

using SnoopCompile

@testset "@snoopl" begin

    @snoopl "func_names.csv" "llvm_timings.yaml" begin
        @eval module M
            i(x) = x+5
            h(a::Array) = i(a[1]::Integer) + 2
            g(y::Integer) = h(Any[y])
        end;
        @eval M.g(3)
    end;

    times, info = SnoopCompile.read_snoopl("func_names.csv", "llvm_timings.yaml")

    @test length(times) == 3  # i(), h(), g()
    @test length(info) == 3  # i(), h(), g()

    rm("func_names.csv")
    rm("llvm_timings.yaml")
end


text:
## ScalarAffineFunction-in-Set
## ScalarQuadraticFunction-in-Set

function _build_constraint(
    model::Optimizer,
    f::F,
    s::MOI.GreaterThan{T},
) where {
    T <: Real,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
}
    return cpo_java_ge(model.inner, _parse(model, f), s.lower)
end

function _build_constraint(
    model::Optimizer,
    f::F,
    s::MOI.LessThan{T},
) where {
    T <: Real,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
}
    return cpo_java_le(model.inner, _parse(model, f), s.upper)
end

function _build_constraint(
    model::Optimizer,
    f::F,
    s::MOI.EqualTo{T},
) where {
    T <: Real,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
}
    return cpo_java_eq(model.inner, _parse(model, f), s.value)
end

function _build_constraint(
    model::Optimizer,
    f::F,
    s::MOI.EqualTo{T},
) where {
    T <: Integer,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
}
    return cpo_java_eq(model.inner, _parse(model, f), Int32(s.value))
end

function _build_constraint(
    model::Optimizer,
    f::F,
    s::MOI.Interval{T},
) where {
    T <: Real,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
}
    if s.lower == Inf
        return _build_constraint(model, f, MOI.LessThan(s.upper))
    elseif s.upper == Inf
        return _build_constraint(model, f, MOI.GreaterThan(s.lower))
    elseif s.lower == s.upper
        return _build_constraint(model, f, MOI.EqualTo(s.lower))
    end

    return cpo_java_range(model.inner, s.lower, _parse(model, f), s.upper)
end

# No vector of constraints, there is no more efficient way to do it.

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{F, S},
) where {
    T <: Real,
    F <: Union{MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}},
    S <:
    Union{MOI.GreaterThan{T}, MOI.LessThan{T}, MOI.EqualTo{T}, MOI.Interval{T}},
}
    cpo_java_remove(model.inner, _info(model, c).constraint)
    delete!(model.constraint_info, c)
    return
end

# TODO: function MOI.set(model::Optimizer, ::MOI.ConstraintSet, c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}, s::S) where {S}

## VectorOfVariables-in-SOS{I|II}
# Not available. Bridge it?

## VectorOfVariables-in-SecondOrderCone

function _build_constraint(
    model::Optimizer,
    f::MOI.VectorOfVariables,
    s::MOI.SecondOrderCone,
)
    if length(f.variables) != s.dimension
        error("Dimension of $(s) does not match number of terms in $(f)")
    end

    # First, check the lower bound on t.
    t_info = _info(model, f.variables[1])
    if !_has_lb(model, t_info.index) || _get_lb(model, t_info.index) < 0.0
        if _get_lb(model, t_info.index) < 0.0
            t_info.old_lb = _get_lb(model, t_info.index)
        end

        t_info.n_socs += 1
        _set_lb(model, t_info.index, 0.0)
    end

    # Then, add the quadratic constraint.
    expr = _parse(model, f, s)
    return cpo_java_gt(model, expr, 0)
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SecondOrderCone},
)
    # Remove the constraint.
    cpo_java_remove(model, _info(model, c).constraint)

    # Maybe restore the old bound on t.
    t_info = _info(model, f.variables[1])
    t_info.n_socs -= 1

    if t_info.n_socs == 0
        _set_lb(model, t_info.index, t_info.old_lb, typeof(t_info.old_lb))
        t_info.old_lb = nothing
    end

    return
end


text:
using Documenter, Literate, PDELib

function make_all()

    makedocs(
        sitename="PDELib.jl",
        modules = [PDELib],
        clean = false,
        doctest = false,
        authors = "J. Fuhrmann, Ch.Merdon. T. Streckenbach",
        repo="https://github.com/j-fu/PDELib.jl",
        pages=[
            "Home"=>"index.md",
            "Introductory Material"=>"intro.md",
        ],
    )

    if !isinteractive()
        deploydocs(repo = "github.com/WIAS-BERLIN/PDELib.jl.git",
                   devbranch = "main"
                   )
    end
end

make_all()


text:
"""
   X, Y, p, t = gensquaremesh(nx::Int64,  ny::Int64 = 50)

to generate a square mesh
"""
function gensquaremesh(nx::Int64,  ny::Int64 = 50)
    X = [x for x ∈  range(-1, stop = 1., length = nx), y = 1:ny]
    Y = [y for x = 1:nx, y ∈  range(-1, stop = 1., length = ny)]
    p = hcat(X[:], Y[:])
    t = delaunay(p)

    nx, ny = 10, 10
    Xs = [x for x ∈  range(-1, stop = 1., length = nx), y = 1:ny]
    Ys = [y for x = 1:nx, y ∈  range(-1, stop = 1., length = ny)]
    ps = hcat(Xs[:], Ys[:])
    ts = delaunay(ps)
    # save display
    @save "/home/oudet/.julia/dev/FemAdjoint/test/data/squaremesh.jld2" p t ps ts
    Z = sin.(2 * X * pi) .* cos.(3 * Y * pi)
    show(heatmap(Z, xfact=.1, yfact=.1, xoffset=-1.5, colormap = :inferno))
    return nothing
end

function genjac(meshfile::String = "diskmesh.jld2";
                addstr::String = "")
    meshfile = "/home/oudet/.julia/dev/FemAdjoint/test/data/" * meshfile
    dictjld = FileIO.load(meshfile)
    p, t = dictjld["p" * addstr], dictjld["t" * addstr]
    if ("jacK" * addstr) ∉ keys(dictjld) # regen jac data
        y = similar(FemAdjoint.assembKM_P12D(p, t)[1])
        fK(y, x) = FemAdjoint.assembK_P12D_inplace!(x, t, y)
        println("recompute sparsity pattern for K")
        @time sparsity_pattern = Symbolics.jacobian_sparsity(fK, y, p[:])
        jac = Float64.(sparse(sparsity_pattern))
        dictjld["jacK" * addstr] = jac
        FileIO.save(meshfile, dictjld)
        display(dictjld)
        println("end export sparsity pattern")
    end
    if ("jacM" * addstr) ∉ keys(dictjld) # regen jac data
        y = similar(FemAdjoint.assembKM_P12D(p, t)[2])
        fM(y, x) = FemAdjoint.assembM_P12D_inplace!(x, t, y)
        println("recompute sparsity pattern for M")
        @time sparsity_pattern = Symbolics.jacobian_sparsity(fM, y, p[:])
        jac = Float64.(sparse(sparsity_pattern))
        dictjld["jacM" * addstr] = jac
        FileIO.save(meshfile, dictjld)
        display(dictjld)
        println("end export sparsity pattern")
    end
    return nothing
end

function plotunicode(U;
                     colormap::Symbol = :viridis, title::String = "")
    np = length(U)
    snp = ceil(Int64, sqrt(np))
    u = reshape(U, snp, snp)
    show(UnicodePlots.heatmap(real.(u), xfact = 2 / (snp - 1),
                              yfact = 2 / (snp - 1),
                              xoffset = -1., yoffset = -1.,
                              colormap = colormap, title = title,
                              #, colormap = :inferno,
                              width = 30, height = 30))
    print("\n")
    return nothing
end


text:
# # Intersection

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/intersection.ipynb)

# In this example we demonstrate how to define a T-shape intersection with `AutomotiveDrivingModels`.
# You will also learn how to define your own custom action type and driver model type.

# ## Generate a T-Shape intersection

# In order to generate the road network, we first initialize a Roadway object.

using AutomotiveDrivingModels
using AutoViz
AutoViz.colortheme["background"] = colorant"white"; # hide
using Random

roadway = Roadway()
#md nothing # hide

# Define coordinates of the entry and exit points to the intersection

r = 5.0 # turn radius
w = DEFAULT_LANE_WIDTH

A = VecSE2(0.0,w,-π)
B = VecSE2(0.0,0.0,0.0)
C = VecSE2(r,-r,-π/2)
D = VecSE2(r+w,-r,π/2)
E = VecSE2(2r+w,0,0)
F = VecSE2(2r+w,w,-π)
#md nothing # hide

# The next step consists in appending all the lanes to the road network.
# We can define a helper function to add a new lane to the roadway.

function append_to_curve!(target::Curve, newstuff::Curve)
    s_end = target[end].s
    for c in newstuff
        push!(target, CurvePt(c.pos, c.s+s_end, c.k, c.kd))
    end
    return target
end
#md nothing # hide

# Example of a lane that consists of 3 road segments, a straight curve
# (from the left to the center), a turning part (right turn) and a final
# straight curve.

# Append right turn coming from the left
curve = gen_straight_curve(convert(VecE2, B+VecE2(-100,0)), convert(VecE2, B), 2)
append_to_curve!(curve, gen_bezier_curve(B, C, 0.6r, 0.6r, 51)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, C), convert(VecE2, C+VecE2(0,-50.0)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

# visualize the current lane constellation
snapshot = render([roadway])
#md write("partial_intersection.svg", snapshot) # hide
#md # ![partial intersection](partial_intersection.svg)

# Let's repeat the process and complete the T-shape intersection

# Append straight left
curve = gen_straight_curve(convert(VecE2, B+VecE2(-100,0)), convert(VecE2, B), 2)
append_to_curve!(curve, gen_straight_curve(convert(VecE2, B), convert(VecE2, E), 2)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, E), convert(VecE2, E+VecE2(50,0)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

# Append straight right
curve = gen_straight_curve(convert(VecE2, F+VecE2(50,0)), convert(VecE2, F), 2)
append_to_curve!(curve, gen_straight_curve(convert(VecE2, F), convert(VecE2, A), 2)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, A), convert(VecE2, A+VecE2(-100,0)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

# Append left turn coming from the right
curve = gen_straight_curve(convert(VecE2, F+VecE2(50,0)), convert(VecE2, F), 2)
append_to_curve!(curve, gen_bezier_curve(F, C, 0.9r, 0.9r, 51)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, C), convert(VecE2, C+VecE2(0,-50)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

# Append right turn coming from below
curve = gen_straight_curve(convert(VecE2, D+VecE2(0,-50)), convert(VecE2, D), 2)
append_to_curve!(curve, gen_bezier_curve(D, E, 0.6r, 0.6r, 51)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, E), convert(VecE2, E+VecE2(50,0)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

# Append left turn coming from below
curve = gen_straight_curve(convert(VecE2, D+VecE2(0,-50)), convert(VecE2, D), 2)
append_to_curve!(curve, gen_bezier_curve(D, A, 0.9r, 0.9r, 51)[2:end])
append_to_curve!(curve, gen_straight_curve(convert(VecE2, A), convert(VecE2, A+VecE2(-100,0)), 2))
lane = Lane(LaneTag(length(roadway.segments)+1,1), curve)
push!(roadway.segments, RoadSegment(lane.tag.segment, [lane]))

snapshot = render([roadway])
#md write("intersection.svg", snapshot) # hide
#md # ![intersection](intersection.svg)

# We can identify each lane thanks to the following user-defined functions.
# We define a `LaneOverlay` object that indicate the lane to highlight.
# One could implement any custom type to display other information on the lane.
# We then add a new method to the `add_renderable!` function that execute the specific
# action (coloring in blue). Look at `Autoviz.jl` for more detail on the function
# `add_renderable!`.

# The following animation iterates over the individual lanes of the intersection
# layout and highlights them:

struct LaneOverlay
    roadway::Roadway
    lane::Lane
    color::Colorant
end
function AutoViz.add_renderable!(rendermodel::RenderModel, overlay::LaneOverlay)
    add_renderable!(rendermodel, overlay.lane, overlay.roadway, color_asphalt=overlay.color)
end

using Reel

animation = roll(fps=1.0, duration=length(roadway.segments)) do t, dt
    i = Int(floor(t/dt)) + 1
    renderables = [
        roadway,
        LaneOverlay(roadway, roadway[LaneTag(i,1)], RGBA(0.0,0.0,1.0,0.5))
    ]
    render(renderables)
end;

#md write("highlighted_lanes.gif", animation) # hide
#md nothing # hide
#md # ![highlighted lanes](highlighted_lanes.gif)

# ## Navigation in the new road network

# Let's populate the intersection

vs0 = VehicleState(B + polar(50.0,-π), roadway, 8.0) # initial state of the vehicle
scene = Scene([Entity(vs0, VehicleDef(), 1)])

snapshot = render([roadway, scene])
#md write("intersection_populated.svg", snapshot) # hide
#md # ![populated intersection](intersection_populated.svg)

# We will use lateral and longitudinal acceleration to control a car in the intersection. The first step is to define a corresponding action type that will contain the acceleration inputs.

struct LaneSpecificAccelLatLon
    a_lat::Float64
    a_lon::Float64
end

# Next, add a method to the propagate function to update the state using our new action type.

function AutomotiveDrivingModels.propagate(veh::Entity, action::LaneSpecificAccelLatLon, roadway::Roadway, Δt::Float64)
    lane_tag_orig = veh.state.posF.roadind.tag
    state = propagate(veh, LatLonAccel(action.a_lat, action.a_lon), roadway, Δt)
    roadproj = proj(state.posG, roadway[lane_tag_orig], roadway, move_along_curves=false)
    retval = VehicleState(Frenet(roadproj, roadway), roadway, state.v)
    return retval
end

# **Driver Model:**

# We define a driver model,
# which can be seen as a distribution over actions.  # TODO
# Here we will define the simplest model, which is to repeat the same action.

struct InterDriver <: DriverModel{LaneSpecificAccelLatLon}
    a::LaneSpecificAccelLatLon
end

AutomotiveDrivingModels.observe!(model::InterDriver, scene::Scene, roadway::Roadway, egoid::Int64) = model
Base.rand(::AbstractRNG, model::InterDriver) = model.a

# **Simulate:**

# First associate a model to each driver in the scene using a dictionary.
# Here we only have one driver identified by its ID: 1.
# Then everything is ready to run the `simulate!` function.

using Reel

timestep = 0.1
nticks = 100

vs0 = VehicleState(B + polar(50.0,-π), roadway, 8.0)
scene = Scene([Entity(vs0, VehicleDef(), 1)])
models = Dict(1 => InterDriver(LaneSpecificAccelLatLon(0.0,0.0)))
scenes = simulate(scene, roadway, models, nticks, timestep)

animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
    i = Int(floor(t/dt)) + 1
    renderables = [roadway, scenes[i]]
    render(renderables)
end

#md write("animated_intersection.gif", animation) # hide
#md # ![animated intersection](animated_intersection.gif)


text:
function downconvert!(
    downconverted_signal::StructArray{Complex{T}, 1},
    signal::StructArray{Complex{ST}, 1},
    carrier_replica::StructArray{Complex{T}, 1},
    start_sample::Integer,
    num_samples::Integer
) where {T, ST}
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    c_re = carrier_replica.re; c_im = carrier_replica.im
    @avx for i = start_sample:num_samples + start_sample - 1
        ds_re[i] = s_re[i] * c_re[i] + s_im[i] * c_im[i]
        ds_im[i] = s_im[i] * c_re[i] - s_re[i] * c_im[i]
    end
    downconverted_signal
end

# This function creates the carrier and downconverts
# directly. It is a little bit faster than first creating the
# replica and then downconverting. This is not in use at the moment
# to be consistent (see next function).
function downconvert!(
    downconverted_signal::StructArray{Complex{T}},
    signal::StructArray{Complex{TS}},
    carrier_frequency,
    sampling_frequency,
    start_phase,
    start_sample,
    num_samples
) where {T, TS}
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    carrier_freq = upreferred(carrier_frequency / Hz)
    sampling_freq = upreferred(sampling_frequency / Hz)
    @avx for i = start_sample:start_sample + num_samples - 1
        c_im, c_re = sincos(T(2π) * ((i - start_sample) * T(carrier_freq) / T(sampling_freq) + T(start_phase)))
        ds_re[i] = s_re[i] * c_re + s_im[i] * c_im
        ds_im[i] = s_im[i] * c_re - s_re[i] * c_im
    end
    downconverted_signal
end

# Same as above but for the multiple antenna case. It is faster
# for less than 4 antennas, but slower otherwise. This is not in
# use due to this circumstances.
function downconvert!(
    downconverted_signal::StructArray{Complex{T}, 2},
    signal::StructArray{Complex{TS}, 2},
    carrier_frequency,
    sampling_frequency,
    start_phase,
    start_sample,
    num_samples
) where {N, T, TS}
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    carrier_freq = upreferred(carrier_frequency / Hz)
    sampling_freq = upreferred(sampling_frequency / Hz)
    @avx for i = start_sample:start_sample + num_samples - 1
        c_im, c_re = sincos(T(2π) * ((i - start_sample) * T(carrier_freq) / T(sampling_freq) + T(start_phase)))
        for j = 1:size(s_re, 2)
            ds_re[i,j] = s_re[i,j] * c_re + s_im[i,j] * c_im
            ds_im[i,j] = s_im[i,j] * c_re - s_re[i,j] * c_im
        end
    end
    downconverted_signal
end

function downconvert!(
    downconverted_signal::StructArray{Complex{T}, 2},
    signal::StructArray{Complex{ST}, 2},
    carrier_replica::StructArray{Complex{T}, 1},
    start_sample::Integer,
    num_samples::Integer
) where {T, ST}
    ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
    s_re = signal.re; s_im = signal.im
    c_re = carrier_replica.re; c_im = carrier_replica.im
    @avx for i = start_sample:num_samples + start_sample - 1, j = 1:size(s_re, 2)
        ds_re[i, j] = s_re[i, j] * c_re[i] + s_im[i, j] * c_im[i]
        ds_im[i, j] = s_im[i, j] * c_re[i] - s_re[i, j] * c_im[i]
    end
    downconverted_signal
end

@static if VERSION >= v"1.6"
    function downconvert!(
        downconverted_signal::StructArray{Complex{T}, 1},
        signal::AbstractArray{Complex{ST}, 1},
        carrier_replica::StructArray{Complex{T}, 1},
        start_sample::Integer,
        num_samples::Integer
    ) where {T, ST}
        signal_real = reinterpret(reshape, ST, signal)
        ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
        c_re = carrier_replica.re; c_im = carrier_replica.im
        @avx for i = start_sample:num_samples + start_sample - 1
            ds_re[i] = signal_real[1, i] * c_re[i] + signal_real[2, i] * c_im[i]
            ds_im[i] = signal_real[2, i] * c_re[i] - signal_real[1, i] * c_im[i]
        end
        downconverted_signal
    end

    function downconvert!(
        downconverted_signal::StructArray{Complex{T}, 2},
        signal::AbstractArray{Complex{ST}, 2},
        carrier_replica::StructArray{Complex{T}, 1},
        start_sample::Integer,
        num_samples::Integer
    ) where {T, ST}
        signal_real = reinterpret(reshape, ST, signal)
        ds_re = downconverted_signal.re; ds_im = downconverted_signal.im
        c_re = carrier_replica.re; c_im = carrier_replica.im
        @avx for i = start_sample:num_samples + start_sample - 1, j = 1:size(signal_real, 3)
            ds_re[i, j] = signal_real[1, i, j] * c_re[i] + signal_real[2, i, j] * c_im[i]
            ds_im[i, j] = signal_real[2, i, j] * c_re[i] - signal_real[1, i, j] * c_im[i]
        end
        downconverted_signal
    end
else
    function downconvert!(
        downconverted_signal,
        signal::AbstractMatrix,
        carrier_replica,
        start_sample::Integer,
        num_samples::Integer
    )
        sample_range = start_sample:start_sample + num_samples - 1
        downconverted_signal[sample_range,:] .= @view(signal[sample_range,:]) .* conj.(@view(carrier_replica[sample_range]))
        downconverted_signal
    end

    function downconvert!(
        downconverted_signal,
        signal::AbstractVector,
        carrier_replica,
        start_sample::Integer,
        num_samples::Integer
    )
        sample_range = start_sample:start_sample + num_samples - 1
        downconverted_signal[sample_range,:] .= @view(signal[sample_range]) .* conj.(@view(carrier_replica[sample_range]))
        downconverted_signal
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "openPMD_api"
version = v"0.13.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openPMD/openPMD-api/archive/refs/tags/0.13.4.tar.gz",
                  "46c013be5cda670f21969675ce839315d4f5ada0406a6546a91ec3441402cf5e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openPMD-api-0.13.4
mkdir build
cd build
if [[ "$target" == *-apple-* ]]; then
    # Set up a wrapper script for the assembler. GCC's assembly
    # output isn't accepted by the LLVM assembler.
    as=$(which as)
    mv "$as" "$as.old"
    export AS="$as.old"
    ln -s "$WORKSPACE/srcdir/scripts/as.llvm" "$as"
fi
mpiopts=
if [[ "$target" == x86_64-w64-mingw32 ]]; then
    mpiopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI -DMPI_C_LIBRARIES=msmpi64 -DMPI_CXX_LIBRARIES=msmpi64"
elif [[ "$target" == *-mingw* ]]; then
    mpiopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI"
fi
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release ${mpiopts} ..
make -j${nproc}
make -j${nproc} install
install_license ../COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# Apparently, macOS doesn't use different C++ string APIs
platforms = expand_cxxstring_abis(platforms; skip=Sys.isapple)

# The products that we will ensure are always built
products = [
    LibraryProduct("libopenPMD", :libopenPMD),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="ADIOS2_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    # We would need a parallel version of HDF5
    # Dependency(PackageSpec(name="HDF5_jll")),
    Dependency(PackageSpec(name="MPICH_jll")),
    Dependency(PackageSpec(name="MicrosoftMPI_jll")),
]

# Build the tarballs, and possibly a `build.jl` as well.
# We need C++14, which requires at least GCC 5.
# GCC 5 an incompatible signatures for `posix_memalign` on linux/musl, fixed on GCC 6
# GCC 5 has a bug regarding `std::to_string` on freebsd, fixed on GCC 6
# macos encounters an ICE in GCC 6; switching to GCC 7 instead
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7")


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((1, 2, 2, 2),),
    ]
insts["fast"] = [
    ((2, 2, 3, 2),),
    ((3, 3, 3, 3),),
    ((3, 3, 5, 4),),
    ((5, 2, 5, 3),),
    ((1, 30, 2, 30),),
    ((10, 1, 3, 1),),
    ]
insts["various"] = [
    ((3, 3, 3, 3),),
    ((3, 3, 5, 4),),
    ((5, 2, 5, 3),),
    ((2, 30, 3, 30),),
    ((2, 15, 6, 15),),
    ((2, 30, 2, 40),),
    ((8, 3, 3, 3),),
    ]
return (PolyEnvelopeJuMP, insts)


text:
import Porta
import Test

samples = 36
markerradius = 0.01
segments = 36
basecenter = [0.0; 0.0; 0.0]
basecolor = [0.3; 0.3; 0.3]
baseradius = 1.0
markercenter = [rand(samples) rand(samples) rand(samples)]
markercolor = [rand(samples) rand(samples) rand(samples)]
basemap = Porta.🌐(basecenter,
                   basecolor,
                   baseradius,
                   markercenter,
                   markercolor,
                   markerradius,
                   segments)
Test.@test size(basemap.basemanifold) == (segments, segments, 3)
Test.@test size(basemap.basecolor) == size(basemap.basemanifold)
Test.@test size(basemap.markermanifold) == (samples, segments, segments, 3)
Test.@test size(basemap.markercolor) == size(basemap.markermanifold)


text:
using JLD
using DSGE
using DifferentialEquations
using ForwardDiff

# Get transition equations
m = KrusellSmithCT()
T, R, C, inverse_basis, basis = DSGE.solve(m)

# Make initial guess
n_vars = n_states(m)
n_expectation_errors = n_shocks_expectational(m)
n_shocks = n_shocks_exogenous(m)
x = zeros(Float64, n_vars)

# Set up SDE
f(u,p,t) = T*u
g(u,p,t) = R
dt = 1 ./ 90 # daily frequency when one period is one quarter
tspan = (0.0,200.) # 50 years of data
u0 = reshape(basis * x, (size(inverse_basis, 2), 1))
W = WienerProcess(0.0, 0., 0.)
prob = SDEProblem(f, g, u0, tspan, W)

# Simulate many trials
trials = Dict()
ntrials = 100
for n = 1:ntrials
    sol = DifferentialEquations.solve(prob, EM(), dt = dt)
    trials[n] = sol.u
end

if ntrials > 1
    save("many_fine_data_red_basis.jld", "trials", trials)
else
    save("fine_data_red_basis.jld", "states", trials[1])
end


text:
@testset "syms" begin
    @syms x u v w

    @testset "allowed number of args" begin
        # Probably, there is an easier test for this possible.
        p_single_arg = @rule P(~u) => ~u + ~u
        @eqtest p_single_arg(P(u)) == (2u)

        p_multiple_args = @rule P(~u, ~v, ~w) => ~u + ~v + ~w
        @eqtest p_multiple_args(P(u, v, w)) == (u + v + w)

        d_single_arg = @rule Do(~x) => ~x + ~x
        @eqtest d_single_arg(Do(x)) == (2x)

        given_single_arg = @rule ¦(~u, ~v) => ~u + ~v
        @eqtest given_single_arg(u ¦ v) == (u + v)
    end
end


text:
using MetacommunityDynamics
using DynamicGrids
using Dispersal
using Distributions
using ColorSchemes
using Plots

consumermodel =
    Eating{Tuple{:C,:R}}(
        functionalresponse=LotkaVolterra(0.3), dt=0.1) +
    AdjacentBernoulliDispersal{:C}(DispersalKernel(radius=2), 0.6) +
    LinearMortality{:C}(0.07);

resourcemodel =
    MetacommunityDynamics.LogisticGrowth{:R}(λ=1.4, K=10., dt=0.1)+
    AdjacentBernoulliDispersal{:R}(DispersalKernel(radius=6), 0.1);

model = resourcemodel + consumermodel

latticesize = 100
initconsumer = zeros(Biomass, latticesize,latticesize)
initresource = zeros(Biomass, latticesize,latticesize)

initr = 58000
initc = 5000

map(i -> initresource[i...] = 10., rand(DiscreteUniform(1, latticesize), (initr,2)))
map(i -> initconsumer[i...] = 5., rand(DiscreteUniform(1, latticesize), (initc,2)))

arrayout = ArrayOutput((C=initconsumer, R=initresource ), tspan=1:1500)
@time sim!(arrayout, model)

## Ploting
Cs = zeros(Float64, length(arrayout))
Rs = zeros(Float64, length(arrayout))
for (i,fr) in enumerate(arrayout)
    Cs[i] = sum(fr[:C])
    Rs[i] = sum(fr[:R])
end
plot(1:length(Cs),Cs, label="consumer")
plot!(1:length(Rs),Rs, label="resource")
xlabel!("time")
ylabel!("biomass")


text:
function assert_allow_step(state)
    if state.broke_on_error
        printstyled(stderr, "Cannot step after breaking on error\n"; color=:red)
        return false
    end
    if state.level != 1
        printstyled(stderr, "Cannot step in a non leaf frame\n"; color=:red)
        return false
    end
    return true
end

function show_breakpoint(io::IO, bp::BreakpointRef)
    outbuf = IOContext(IOBuffer(), io)
    if bp.err === nothing
        print(outbuf, "Hit breakpoint:\n")
    else
        print(outbuf, "Breaking for error:\n")
        Base.display_error(outbuf, bp.err, [])
    end
    print(io, String(take!(outbuf.io)))
end

function execute_command(state::DebuggerState, ::Union{Val{:c},Val{:nc},Val{:n},Val{:se},Val{:s},Val{:si},Val{:sg},Val{:so}}, cmd::AbstractString)
    assert_allow_step(state) || return false
    cmd == "so" && (cmd = "finish")
    ret = debug_command(state.mode, state.frame, Symbol(cmd))
    if ret === nothing
        state.overall_result = get_return(root(state.frame))
        state.frame = nothing
        return false
    else
        state.frame, pc = ret
        if pc isa BreakpointRef
            if pc.stmtidx != 0 # This is the dummy breakpoint to stop just after entering a call
                if state.terminal !== nothing # fix this, it happens when a test hits this and hasnt set a terminal
                    show_breakpoint(Base.pipe_writer(state.terminal), pc)
                end
            end
            if pc.err !== nothing
                state.broke_on_error = true
            end
        end
        if pc_expr(state.frame) === nothing # happens when stopping for @bp
            JuliaInterpreter.maybe_next_call!(state.frame)
        end
        return true
    end
end

function execute_command(state::DebuggerState, ::Val{:bt}, cmd)
    num = 0
    frame = state.frame
    while frame !== nothing
        num += 1
        print_frame(Base.pipe_writer(state.terminal), num, frame)
        frame = caller(frame)
    end
    println()
    return false
end

function stacklength(frame)
    s = 0
    while frame !== nothing
        s += 1
        frame = caller(frame)
    end
    return s
end

execute_command(state::DebuggerState, ::Val{:st}, cmd) = true

function execute_command(state::DebuggerState, ::Union{Val{:f}, Val{:fr}}, cmd)
    subcmds = split(cmd, ' ')
    if length(subcmds) == 1
        new_level = 1
    else
        new_level = tryparse(Int, subcmds[2])
        if new_level == nothing
            printstyled(stderr, "Failed to parse $(repr(subcmds[2])) as an integer\n"; color=:red)
            return false
        end
    end

    if new_level > stacklength(state.frame) || new_level < 1
        printstyled(stderr, "Not a valid frame index\n"; color=:red)
        return false
    end

    if subcmds[1] == "fr"
        print_frame(Base.pipe_writer(state.terminal), new_level, state.frame)
        return false
    else
        state.level = new_level
        return true
    end
end

function execute_command(state::DebuggerState, ::Val{:w}, cmd::AbstractString)
    # TODO show some info messages?
    cmds = split(cmd)
    if length(cmds) == 1
        io = Base.pipe_writer(state.terminal)
        show_watch_list(io, state)
        return false
    elseif length(cmds) >= 2
        if cmds[2] == "rm"
            if length(cmds) == 2
                clear_watch_list!(state)
                return false
            elseif length(cmds) == 3
                i = parse(Int, cmds[3])
                clear_watch_list!(state, i)
                return false
            end
        end
        if cmds[2] == "add"
            if add_watch_entry!(state, join(cmds[3:end]))
            end
            return false
        end
    end
    # Error
    return execute_command(state, Val(:_), cmd)
end

function execute_command(state::DebuggerState, _, cmd)
    println("Unknown command `$cmd`. Executing `?` to obtain help.")
    execute_command(state, Val{Symbol("?")}(), "?")
end

function execute_command(state::DebuggerState, ::Val{:?}, cmd::AbstractString)
    display(
            @md_str """
    Basic Commands:\\
    - `st`: show the status\\
    - `n`: step to the next line\\
    - `s`: step into the next call\\
    - `so`: step out of the current call\\
    - `c`: continue execution until a breakpoint is hit\\
    - `bt`: show a simple backtrace\\
    - ``` `stuff ```: run `stuff` in the current function's context\\
    - `fr [v::Int]`: show all variables in the current frame, `v` defaults to `1`\\
    - `f [n::Int]`: go to the `n`-th frame\\
    - `w`\\
        - `w add expr`: add an expression to the watch list\\
        - `w`: show all watch expressions evaluated in the current function's context\\
        - `w rm [i::Int]`: remove all or the `i`:th watch expression\\
    - `o`: open the current line in an editor\\
    - `q`: quit the debugger, returning `nothing`\\
    - `C`: toggle compiled mode\\
    Advanced commands:\\
    - `nc`: step to the next call\\
    - `se`: step one expression step\\
    - `si`: same as `se` but step into a call if a call is the next expression\\
    - `sg`: step into a generated function\\
    """)
    return false
end

function execute_command(state::DebuggerState, ::Val{:o}, cmd::AbstractString)
    frame = active_frame(state)
    loc = JuliaInterpreter.whereis(frame)
    if loc === nothing
        printstyled(stderr, "Could not find source location\n"; color=:red)
        return false
    end
    file, line = loc
    if !isfile(file)
        printstyled(stderr, "Could not find file: $(repr(file))\n"; color=:red)
        return false
    end
    InteractiveUtils.edit(file, line)
    return false
end