text:
@inline function _colorediteration!(J,sparsity,rows_index,cols_index,vfx,colorvec,color_i,ncols)
    @inbounds for i in 1:length(cols_index)
        if colorvec[cols_index[i]] == color_i
            J[rows_index[i],cols_index[i]] = vfx[rows_index[i]]
        end
    end
end

@inline function _colorediteration!(J,sparsity::SparseMatrixCSC,rows_index,cols_index,vfx,colorvec,color_i,ncols)
    @inbounds for col_index in 1:ncols
        if colorvec[col_index] == color_i
            @inbounds for row_index in view(sparsity.rowval,sparsity.colptr[col_index]:sparsity.colptr[col_index+1]-1)
                J[row_index,col_index]=vfx[row_index]
            end
        end
    end
end

# fast version for the case where J and sparsity have the same sparsity pattern
@inline function _colorediteration!(Jsparsity::SparseMatrixCSC,vfx,colorvec,color_i,ncols)
    @inbounds for col_index in 1:ncols
        if colorvec[col_index] == color_i
            @inbounds for spidx in nzrange(Jsparsity, col_index)
                row_index = Jsparsity.rowval[spidx]
                Jsparsity.nzval[spidx]=vfx[row_index]
            end
        end
    end
end

#override default setting of using findstructralnz
_use_findstructralnz(sparsity) = ArrayInterface.has_sparsestruct(sparsity)
_use_findstructralnz(::SparseMatrixCSC) = false

# test if J, sparsity are both SparseMatrixCSC and have the same sparsity pattern of stored values
_use_sparseCSC_common_sparsity(J, sparsity) = false
_use_sparseCSC_common_sparsity(J::SparseMatrixCSC, sparsity::SparseMatrixCSC) =
    ((J.colptr == sparsity.colptr) && (J.rowval == sparsity.rowval))

function __init__()
    @require BlockBandedMatrices="ffab5731-97b5-5995-9138-79e8c1846df0" begin
        @require BlockArrays="8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin
            _use_findstructralnz(::BlockBandedMatrices.BandedBlockBandedMatrix) = false
            _use_findstructralnz(::BlockBandedMatrices.BlockBandedMatrix) = false

            @inline function _colorediteration!(Jac::BlockBandedMatrices.BandedBlockBandedMatrix,
                                                sparsity::BlockBandedMatrices.BandedBlockBandedMatrix,
                                                rows_index,cols_index,vfx,colorvec,color_i,ncols)
                λ,μ = BlockBandedMatrices.subblockbandwidths(Jac)
                rs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) # column block sizes
                cs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1))
                b = BlockBandedMatrices.BlockArray(vfx,rs)
                c = BlockBandedMatrices.BlockArray(colorvec,cs)
                @inbounds for J=BlockArrays.blockaxes(Jac,2)
                    c_v = c.blocks[J.n[1]]
                    @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J)
                        V = view(Jac,K,J)
                        b_v = b.blocks[K.n[1]]
                        data = BlockBandedMatrices.bandeddata(V)
                        p = pointer(data)
                        st = stride(data,2)
                        m,n = size(V)
                        @inbounds for j=1:n
                            if c_v[j] == color_i
                                @inbounds for k=max(1,j-μ):min(m,j+λ)
                                    unsafe_store!(p, b_v[k], (j-1)*st + μ + k - j + 1)
                                end
                            end
                        end
                    end
                end
            end

            @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix,
                                                sparsity::BlockBandedMatrices.BlockBandedMatrix,
                                                rows_index,cols_index,vfx,colorvec,color_i,ncols)
                rs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) # column block sizes
                cs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1))
                b = BlockBandedMatrices.BlockArray(vfx,rs)
                c = BlockBandedMatrices.BlockArray(colorvec,cs)
                @inbounds for J=BlockArrays.blockaxes(Jac,2)
                    c_v = c.blocks[J.n[1]]
                    blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J)
                    _,n = length.(getindex.(axes(Jac), (blockcolrange[1], J)))
                    @inbounds for j = 1:n
                        if c_v[j] == color_i
                            @inbounds for K = blockcolrange
                                V = view(Jac,K,J)
                                b_v = b.blocks[K.n[1]]
                                m = size(V,1)
                                @inbounds for k = 1:m
                                    V[k,j] = b_v[k]
                                end
                            end
                        end
                    end
                end
            end

        end

        @require BandedMatrices = "aae01518-5342-5314-be14-df237901396f" begin

            _use_findstructralnz(::BandedMatrices.BandedMatrix) = false

            @inline function _colorediteration!(Jac::BandedMatrices.BandedMatrix,
                                                sparsity::BandedMatrices.BandedMatrix,
                                                rows_index,cols_index,vfx,colorvec,color_i,ncols)
                nrows = size(Jac,1)
                l,u = BandedMatrices.bandwidths(Jac)
                #data = BandedMatrices.bandeddata(Jac)
                @inbounds for col_index in max(1,1-l):min(ncols,ncols+u)
                    if colorvec[col_index] == color_i
                        @inbounds for row_index in max(1,col_index-u):min(nrows,col_index+l)
                            #data[u+row_index-col_index+1,col_index] = vfx[row_index]
                            Jac[row_index,col_index]=vfx[row_index]
                        end
                    end
                end
            end
        end
    end
end


text:
  module vis

  using Reexport

  using Luxor
import Dates
using ProgressMeter

macro luxoremacs(body, width=600, height=600)
    quote
        path, io = Base.Filesystem.mktemp()
        close(io)

        path * ".png"
        lfname = path * ".png"
        # local lfname = Luxor._add_ext($(esc(fname)), :png)
        Drawing($width, $height, lfname)
        origin()
        background("white")
        sethue("black")
        $(esc(body))
        finish()
        # instead of preview, I'm printing out something
        println("$lfname")
        println("#<Image: $lfname>")
        # preview()
    end
end

function visualize(x, y, w, h, R)
    # shift all
    x = x .- R.xmax / 2
    y = y .- R.ymax / 2

    # apply scale
    scale = 500 / max(R.xmax, R.ymax)
    x = x .* scale
    y = y .* scale
    w = w .* scale
    h = h .* scale

    @svg begin
        setdash("dash")
        # the point is the center of box
        box(Luxor.Point(0,0), R.xmax * scale, R.ymax * scale, :stroke)
        setdash("solid")
        @showprogress 0.1 "drawing .." for i in 1:length(x)
            # this x y is center
            box(Luxor.Point(x[i],y[i]), w[i], h[i], :stroke)
        end
    end
end

  end

ext: jl

max_stars_repo_name: bhdl/bhdl

max_issues_repo_name: bhdl/bhdl

max_forks_repo_name: bhdl/bhdl


text:
export
  # NDTensors.jl
  insertblock!!,
  setindex,
  setindex!!,

  # blocksparse/blockdims.jl
  BlockDims,
  blockdim,
  blockdims,
  nblocks,
  blockindex,

  # blocksparse/blocksparse.jl
  # Types
  Block,
  BlockOffset,
  BlockOffsets,
  BlockSparse,
  # Methods
  blockoffsets,
  blockview,
  eachnzblock,
  findblock,
  isblocknz,
  nnzblocks,
  nnz,
  nzblock,
  nzblocks,

  # blocksparse/blocksparsetensor.jl
  # Types
  BlockSparseTensor,
  # Methods
  blockview,
  insertblock!,
  randomBlockSparseTensor,

  # dense.jl
  # Types
  Dense,
  DenseTensor,
  # Symbols
  ⊗,
  # Methods
  randomTensor,
  randomDenseTensor,
  array,
  contract,
  matrix,
  outer,
  permutedims!!,
  read,
  vector,
  write,

  # diag.jl
  # Types
  Diag,
  DiagTensor,

  # empty.jl
  EmptyStorage,
  EmptyTensor,
  EmptyBlockSparseTensor,

  # tensorstorage.jl
  data,
  TensorStorage,
  randn!,
  scale!,
  norm,

  # tensor.jl
  Tensor,
  tensor,
  inds,
  ind,
  store


text:
using MackeyGlass

X, T = MGGenerator();


text:
# Initial model in Julia

""" Parameters """

# Inertia of arm
M1  = 95.1089;

# Viscous friction
Fv1 = 203.5034;

# Coulomb friction
Fc1 = 20.3935;

# Offset
OF1 = -3.1648;

# Direct dynamic model (DDM)

# Inverse dynamic model (IDM)
τ_idm =

max_stars_repo_name: wmkouw/NSI-EMPS

max_issues_repo_name: wmkouw/NSI-EMPS

max_forks_repo_name: wmkouw/NSI-EMPS


text:
using ModelingToolkit
using Test

# Basic electric components
@parameters t
function Pin(name)
    @variables v(t) i(t)
    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
end

function Ground(name)
    p = Pin(name)
    eqs = [p.v ~ 0]
    ODESystem(eqs, t, [], [], systems=[p], name=name)
end

function OnePort(name)
    p = Pin(:p)
    n = Pin(:n)
    @variables v(t) i(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    v, i, ODESystem(eqs, t, [v, i], [], systems=[p, n], name=name)
end

function ConstantVoltage(name; V = 1.0)
    val = V
    v, i, op = OnePort(name)
    @parameters V
    push!(op.eqs, v ~ V)
    V = ModelingToolkit.value(V)
    push!(op.ps, V)
    op.default_p[V] = val
    op
end

function Resistor(name; R = 1.0)
    val = R
    v, i, op = OnePort(name)
    @parameters R
    push!(op.eqs, v ~ R * i)
    R = ModelingToolkit.value(R)
    push!(op.ps, R)
    op.default_p[R] = val
    op
end

function Capacitor(name; C = 1.0)
    val = C
    v, i, op = OnePort(name)
    @parameters C
    D = Differential(t)
    push!(op.eqs, D(v) ~ i / C)
    C = ModelingToolkit.value(C)
    push!(op.ps, C)
    op.default_p[C] = val
    op
end

function Inductor(name; L = 1.0)
    val = L
    v, i, op = OnePort(name)
    @parameters L
    D = Differential(t)
    push!(op.eqs, D(i) ~ v / L)
    L = ModelingToolkit.value(L)
    push!(op.ps, L)
    op.default_p[L] = val
    op
end

R1 = Resistor(:R1)
C1 = Capacitor(:C1)
I1 = Inductor(:I1)
@test isequal(R1.p.i, (@variables R1₊p₊i(t))[1])
@parameters R
@test ModelingToolkit.default_p(flatten(R1)) == Dict(R => 1.0)
@variables n₊v(t) p₊v(t) n₊i(t) p₊i(t)
@test ModelingToolkit.default_u0(flatten(R1)) == Dict(
                                                      n₊v => 1.0,
                                                      p₊v => 1.0,
                                                      n₊i => 1.0,
                                                      p₊i => 1.0,
                                                     )


text:
using Manopt, Manifolds, ManifoldsBase, LinearAlgebra, Test

@testset "Conjugate Gradient coefficient rules" begin
    F(x) = norm(x)^2
    ∇F(x) = 2 * x
    M = Euclidean(2)
    P = GradientProblem(M, F, ∇F)
    x0 = [0.0, 1.0]
    sC = StopAfterIteration(1)
    s = ConstantStepsize(1.0)
    retr = ExponentialRetraction()
    vtm = ParallelTransport()

    ∇1 = [1.0, 1.0]
    δ1 = [0.0, 2.0]
    ∇2 = [1.0, 1.5]
    δ2 = [0.5, 2.0]
    diff = ∇2 - ∇1

    dU = SteepestDirectionUpdateRule()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    @test O.coefficient(P, O, 1) == 0

    dU = ConjugateDescentCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == dot(∇2, ∇2) / dot(-δ2, ∇1)

    dU = DaiYuanCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == dot(∇2, ∇2) / dot(δ2, ∇2 - ∇1)

    dU = FletcherReevesCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P, O, 1) == 1.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == dot(∇2, ∇2) / dot(∇1, ∇1)

    dU = HagerZhangCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    denom = dot(δ1, diff)
    ndiffsq = dot(diff, diff)
    @test O.coefficient(P, O, 2) ==
          dot(diff, ∇2) / denom - 2 * ndiffsq * dot(δ1, ∇2) / denom^2

    dU = HeestenesStiefelCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == dot(diff, ∇2) / dot(δ1, diff)

    dU = LiuStoreyCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == -dot(∇2, diff) / dot(δ1, ∇1)

    dU = PolakRibiereCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P, O, 1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P, O, 2) == dot(∇2, diff) / dot(∇1, ∇1)
end
@testset "Conjugate Gradient runs – Low Rank matrix approx" begin
    A = Diagonal([2.0, 1.1, 1.0])
    M = Sphere(size(A, 1) - 1)
    F(x) = x' * A * x
    euclidean_∇F(x) = 2 * A * x
    ∇F(x) = project(M, x, euclidean_∇F(x))

    x0 = [2.0, 0.0, 2.0] / sqrt(8.0)
    x_opt = conjugate_gradient_descent(
        M,
        F,
        ∇F,
        x0;
        stepsize=ArmijoLinesearch(),
        coefficient=FletcherReevesCoefficient(),
        stopping_criterion=StopAfterIteration(15),
    )
    @test isapprox(F(x_opt), minimum(eigvals(A)); atol=2.0 * 10^-4)
    @test isapprox(x_opt, eigvecs(A)[:, size(A, 1)]; atol=3.0 * 10^-2)
    x_opt2 = conjugate_gradient_descent(
        M,
        F,
        ∇F,
        x0;
        stepsize=ArmijoLinesearch(),
        coefficient=FletcherReevesCoefficient(),
        stopping_criterion=StopAfterIteration(15),
        return_options=true,
    )
    @test get_solver_result(x_opt2) == x_opt
end


text:
using Test, aoc2019

@testset "fuel_required" begin
    @test fuel_required(12) == 2
    @test fuel_required(14) == 2
    @test fuel_required(1969) == 654
    @test fuel_required(100756) == 33583
    @test fuel_required(2) == 0
end

@testset "fuel_required_including_fuel" begin
    @test fuel_required_including_fuel(14) == 2
    @test fuel_required_including_fuel(1969) == 966
    @test fuel_required_including_fuel(100756) == 50346
end

@testset "d01_answers" begin
    @test aoc2019.d01_inputpath isa String
    input = aoc2019.d01_inputread(aoc2019.d01_inputpath)
    @test length(input) == 100
    @test aoc2019.d01_answer1(input) == 3334297
    @test aoc2019.d01_answer2(input) == 4998565
end


text:
module CookieRequest

import ..Layer, ..request
using ..URIs
using ..Cookies
using ..Pairs: getkv, setkv
import ..@debug, ..DEBUG_LEVEL

const default_cookiejar = Dict{String, Set{Cookie}}()

"""
    request(CookieLayer, method, ::URI, headers, body) -> HTTP.Response

Add locally stored Cookies to the request headers.
Store new Cookies found in the response headers.
"""

abstract type CookieLayer{Next <: Layer} <: Layer end
export CookieLayer

function request(::Type{CookieLayer{Next}},
                 method::String, url::URI, headers, body;
                 cookiejar::Dict{String, Set{Cookie}}=default_cookiejar,
                 kw...) where Next

    hostcookies = get!(cookiejar, url.host, Set{Cookie}())

    cookies = getcookies(hostcookies, url)
    if !isempty(cookies)
        setkv(headers, "Cookie", string(getkv(headers, "Cookie", ""), cookies))
    end

    res = request(Next, method, url, headers, body; kw...)

    setcookies(hostcookies, url.host, res.headers)

    return res
end

function getcookies(cookies, url)

    tosend = Vector{Cookie}()
    expired = Vector{Cookie}()

    # Check if cookies should be added to outgoing request based on host...
    for cookie in cookies
        if Cookies.shouldsend(cookie, url.scheme == "https",
                              url.host, url.path)
            t = cookie.expires
            if t != Dates.DateTime() && t < Dates.now(Dates.UTC)
                @debug 1 "Deleting expired Cookie: $cookie.name"
                push!(expired, cookie)
            else
                @debug 1 "Sending Cookie: $cookie.name to $url.host"
                push!(tosend, cookie)
            end
        end
    end
    setdiff!(cookies, expired)
    return tosend
end

function setcookies(cookies, host, headers)
    for (k,v) in filter(x->x[1]=="Set-Cookie", headers)
        @debug 1 "Set-Cookie: $v (from $host)"
        push!(cookies, Cookies.readsetcookie(host, v))
    end
end

end # module CookieRequest


text:
using Test, LinearAlgebra
using DiffEqOperators

@testset "utility functions" begin
    @test DiffEqOperators.unit_indices(2) == (CartesianIndex(1,0), CartesianIndex(0,1))
    @test DiffEqOperators.add_dims(zeros(2,2), ndims(zeros(2,2)) + 2) == [6. 6.; 0. 0.; 0. 0.]
    @test DiffEqOperators.perpindex(collect(1:5), 3) == [1, 2, 4, 5]
    @test DiffEqOperators.perpsize(zeros(2,2,3,2), 3) == (2, 2, 2)
end


text:
using RegularizedRegression: MatrixDataSet, indepvars, depvar

@testset "MatrixDataSet" begin

@testset "Default names for indep and dep var(s), unless set" begin
    N = 10
    P = 3
    X = rand(P, N)
    y = X[1, :] .* 2.0

    m1 = MatrixDataSet(X, y)
    @test indepvars(m1) == [:x1, :x2, :x3]
    @test depvar(m1)    == :y

    m2 = MatrixDataSet(X, y; depvar = :othername)
    @test indepvars(m2) == [:x1, :x2, :x3]
    @test depvar(m2)    == :othername

    m3 = MatrixDataSet(X, y; indepvars = [:myname1, :x2, :n3])
    @test indepvars(m3) == [:myname1, :x2, :n3]
    @test depvar(m3)    == :y
end

end


text:
# This file is part of Kpax3. License is MIT.

include("functions.jl")
include("histograms.jl")
include("images.jl")
include("diagnostics.jl")


text:
using CSV, DataFrames, JLD, FileIO

function unir(archivos::Array, columnas::Array, nombre_archivo::String = "new_PuMAC.csv", faltantes::Bool = true)

    # Asegurando que los valores introducidos sean los esperados
    if(size(archivos,1) < 1)
        error("Pon archivos csv para extraer, debe de ser un arreglo de tipo String")
    end
    #=
    for i in archivos
        if(contains(i, ".csv"))
            error("Algun archivo no es csv o esta mal escrito, asegurarse que sea de la forma [\"Archivo_1.csv\", \"Archivo_2.csv\", ...]")
        elseif(contains(i, ".jld"))
            error("Algun archivo no es jld o esta mal escrito, asegurarse que sea de la forma [\"Archivo_1.jld\", \"Archivo_2.jld\", ...]")
        end
    end
    =#
    if(size(columnas,1) < 1)
        error("Pon las columnas que quieres extraer, debe de ser un arreglo de tipo String")
    end
       #=
    if(!contains(nombre_archivo, ".jld"))
        error("El String que contiene al archivo necesita tener el tipo de archivo especificado o tipo de archivo no soportado")
    end
    =#
    # Se declara un arreglo con la misma cantidad de espacio que la cantidad de archivos pasados por el usuario
    arreglo = Array{DataFrame}(undef,size(archivos,1))
    tamanio = 0

    # Se extraen todos los archivos csv como DataFrames y se guardan en el arreglo,
    # tambien se obtiene la dimension del DataFrame mas grande
    for i in 1:size(archivos,1)
        arreglo[i] = DataFrame(CSV.File(archivos[i]))
        if tamanio < size(arreglo[i],1)
            tamanio = size(arreglo[i],1)
        end
    end

    # Se declara un Dataframe con la misma cardinalidad del DataFrame mas grande
    df_nuevo = DataFrame(ПУМАК = 1:tamanio)

    for i in 1:size(archivos,1), j in 1:size(names(arreglo[i]),1), k in columnas
        # Se obtienen solo las columnas deseadas por el usuario
        if (k == names(arreglo[i])[j])
            # Se asigna a x con los calores de la columna deseada por el usuario
            x = arreglo[i][:,k]
            # Se convierte a x en tipo Any
            x = convert(Array{Any,1}, x)

            # Se normalizan filas de DataFrames con cardinalidad menor que la mas grande,
            # si faltantes es true entonces los valores faltantes se normalizaran con el valor de missing, de otro modo se normalizara con ceros
            if(nrow(arreglo[i]) < tamanio)
                if (faltantes)
                    append!(x, Array{Missing}(undef, tamanio-nrow(arreglo[i])))
                else
                    append!(x, zeros(tamanio-nrow(arreglo[i])))
                end
            end

            # Se agrega la columna x al DataFrame
            df_nuevo = hcat(df_nuevo, DataFrame(Symbol(k) => x), makeunique=true)
        end
    end

    # Se retira la columna auxiliar "ПУМАК" del DataFrame
    select!(df_nuevo, Not("ПУМАК"))

    # Se exporta al tipo de archivo deseado
    if(contains(nombre_archivo, ".csv"))
        CSV.write(nombre_archivo , df_nuevo)
        print("el archivo se guardo en: ")
        println( joinpath( pwd(), nombre_archivo))

    elseif(contains(nombre_archivo, ".jld"))
        save(nombre_archivo, "df", df_nuevo)
        print("el archivo se guardo en: ")
        println(joinpath(pwd(), nombre_archivo))

    end

    return df_nuevo
end

function unir(carpeta::String, columnas::Array, nombre_archivo::String = "new_PuMAC.csv", faltantes::Bool = true)
    # Se extraen los archivos de tipo csv de la carpeta
    cd(carpeta)
    archivos = filter(x->endswith(x, ".csv"), readdir())

    unir(archivos, columnas, nombre_archivo, faltantes)
end


text:
module PkgOptions

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) PkgOptions

using Base: PkgId
using Pkg: TOML
using PkgOptionsLoader: optionspath
using UUIDs

aspkgid(pkg::PkgId) = pkg
aspkgid((name, uuid)::Tuple{Symbol, Any}) = PkgId(UUID(uuid), String(name))
aspkgid(pkg::Symbol) = Base.identify_package(String(pkg))

function _optionspath(pkg)
    pkgid = aspkgid(pkg)
    if pkgid === nothing
        error("Package $pkg cannot be identified.  Please check that package ",
              pkg, " is installed and importable.")
    end
    return optionspath(pkgid)
end

function load(pkg)
    tomlpath = _optionspath(pkg)
    if tomlpath === nothing || !isfile(tomlpath)
        return Dict{String,Any}()
    end
    return TOML.parsefile(tomlpath)
end

function set(pkg, options::AbstractDict)
    tomlpath = _optionspath(pkg)
    dirpath = dirname(tomlpath)
    isdir(dirpath) || mkpath(dirpath)
    open(tomlpath, "w") do io
        TOML.print(io, options)
    end
    return nothing
end

function set(pkg, option::Pair, options::Pair...)
    stroptions = (string(key) => value for (key, value) in (option, options...))
    newoptions = foldl(push!, stroptions; init=load(pkg))
    set(pkg, newoptions)
end

function set(pkg; options...)
    if isempty(options)
        error("Require options: `set(pkg; OPTION_1=VALUE_1, ...)`")
    end
    set(pkg, options...)
end

rm(pkg) = Base.rm(_optionspath(pkg), force=true)

"""
    PkgOptions.set(pkg::Symbol, options::Pair...)
    PkgOptions.set(pkg::Symbol; options...)

Set and/or add `options` to package options for package `pkg`.

    PkgOptions.set(pkg::Symbol, options::AbstractDict)

Set package options of package `pkg` to `options`, removing the ones not
specified by `options`.
"""
function set end

"""
    PkgOptions.rm(pkg::Symbol)

Remove package options (use the defaults) of package `pkg`.
"""
function rm end

"""
    PkgOptions.load(pkg::Symbol)
"""
function load end

end # module


text:
# ---
# title: 225. Implement Stack using Queues
# id: problem225
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Stack, Design
# link: <https://leetcode.com/problems/implement-stack-using-queues/description/>
# hidden: true
# ---
#
# Implement a last in first out (LIFO) stack using only two queues. The
# implemented stack should support all the functions of a normal queue (`push`,
# `top`, `pop`, and `empty`).
#
# Implement the `MyStack` class:
#
#   * `void push(int x)` Pushes element x to the top of the stack.
#   * `int pop()` Removes the element on the top of the stack and returns it.
#   * `int top()` Returns the element on the top of the stack.
#   * `boolean empty()` Returns `true` if the stack is empty, `false` otherwise.
#
# **Notes:**
#
#   * You must use **only** standard operations of a queue, which means only `push to back`, `peek/pop from front`, `size`, and `is empty` operations are valid.
#   * Depending on your language, the queue may not be supported natively. You may simulate a queue using a list or deque (double-ended queue), as long as you use only a queue's standard operations.
#
# **Follow-up:** Can you implement the stack such that each operation is
# **[amortized](https://en.wikipedia.org/wiki/Amortized_analysis)** `O(1)` time
# complexity? In other words, performing `n` operations will take overall `O(n)`
# time even if one of those operations may take longer.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MyStack", "push", "push", "top", "pop", "empty"]
#     [[], [1], [2], [], [], []]
#     **Output**
#     [null, null, null, 2, 2, false]
#
#     **Explanation**
#     MyStack myStack = new MyStack();
#     myStack.push(1);
#     myStack.push(2);
#     myStack.top(); // return 2
#     myStack.pop(); // return 2
#     myStack.empty(); // return False
#
#
#
#
# **Constraints:**
#
#   * `1 <= x <= 9`
#   * At most `100` calls will be made to `push`, `pop`, `top`, and `empty`.
#   * All the calls to `pop` and `top` are valid.
#
#
## @lc code=start
using LeetCode

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


text:
@testset "ShowMakie" begin
    @testset "ImageClassificationSingle" begin
        task = ImageClassificationSingle((16, 16), [1, 2])
        FastAI.test_task_show(task, ShowMakie())
    end
    @testset "ImageClassificationMulti" begin
        task = ImageClassificationMulti((16, 16), [1, 2])
        FastAI.test_task_show(task, ShowMakie())
    end
    @testset "ImageSegmentation" begin
        task = ImageSegmentation((16, 16), [1, 2])
        FastAI.test_task_show(task, ShowMakie())
    end
    @testset "ImageKeypointRegression" begin
        task = ImageKeypointRegression((16, 16), 10)
        FastAI.test_task_show(task, ShowMakie())
    end
end


text:
# --- Coupling Model Creation --- #

"""
    Coupling(aero::Wagner, stru::Section)

Coupling model for coupling an unsteady aerodynamic model based on Wagner's function (see
[`Wagner`](@ref)) and a two-degree of freedom typical section model (see [`Section`](@ref)).
This coupling introduces the freestream velocity ``U_\\infty``, air density
``\\rho_\\infty`` and air speed of sound ``c`` as additional parameters.
"""
function Coupling(aero::Wagner, ::Section)

    # model constants
    C1 = aero.C1
    C2 = aero.C2

    # coupling function
    g = (dx, x, p, t) -> wagner_section_inputs(dx, x, p, t; C1, C2)

    # number of states, inputs, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(6)
    ny = Val(5)
    np = Val(14)

    # number of parameters introduced by the coupling (use Val(N) to use inferrable dimensions)
    npc = Val(3)

    # jacobian definitions
    ratejac = Nonlinear() # TODO: define rate jacobian function
    statejac = Nonlinear() # TODO: define state jacobian function
    paramjac = Nonlinear() # TODO: define parameter jacobian function
    tgrad = Zeros()

    # convenience function for setting coupling parameters
    setparam = wagner_section_set_parameters!

    # convenience function for separating coupling parameters
    sepparam = wagner_section_separate_parameters

    # return resulting coupling
    return Coupling{false}(g, nx, ny, np, npc;
        ratejac = ratejac,
        statejac = statejac,
        paramjac = paramjac,
        tgrad = tgrad,
        setparam = setparam,
        sepparam = sepparam)
end

# --- Internal Methods --- #

# coupling function
function wagner_section_inputs(dx, x, p, t; C1, C2)
    # extract rate variables
    dλ1, dλ2, dh, dθ, dhdot, dθdot = dx
    # extract state variables
    λ1, λ2, h, θ, hdot, θdot = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ, c = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)
    # calculate loads
    N, A, M = wagner_loads(a, b, ρ, c, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)
    # lift is approximately normal force
    L = N
    # return inputs
    return SVector(u, v, ω, L, M)
end

# convenience function for defining the coupling function parameters
function wagner_section_set_parameters!(p; U, rho, c)
    p[1] = U
    p[2] = rho
    p[3] = c
    return p
end

# convenience function for separating the coupling function parameters
wagner_section_separate_parameters(p) = (U = p[1], rho = p[2], c = p[3])


text:
isdefined(:FourierOptics) || include("../src/FourierOptics.jl")

module TestFourierOptics

using Base.Test

using FourierOptics
using FourierOptics.CoordinateTransforms
using FourierOptics.Regions

Base.isapprox(x::NTuple{2,Real}, y::NTuple{2,Real}; kwds...) =
    (isapprox(x[1], y[1]; kwds...) && isapprox(x[2], y[2]; kwds...))

Base.isapprox(x::NTuple{4,Real}, y::NTuple{4,Real}; kwds...) =
    (isapprox(x[1], y[1]; kwds...) && isapprox(x[2], y[2]; kwds...) &&
     isapprox(x[3], y[3]; kwds...) && isapprox(x[4], y[4]; kwds...))

@testset "Coordinate Transforms" begin
    stp, x0, y0 = 1.2, -2.7, 3.1
    t = (0.7, -1.8)
    p = (0.9, 11.0)
    α = 0.17
    I = CoordinateTransform()
    A = CoordinateTransform(stp, x0, y0)
    B = inv(A)
    @test origin(A) ≈ (x0, y0)
    @test step(A) ≈ stp
    @test A ≈ (x0,y0) + stp*I
    @test A∘B ≈ I  atol=4*eps(Float64)
    @test B∘A ≈ I  atol=4*eps(Float64)
    @test A/A ≈ I  atol=4*eps(Float64)
    @test A\A ≈ I  atol=4*eps(Float64)
    B = t + A
    @test B(p) ≈ t .+ A(p)
    @test (t - A)(p) ≈ t .- A(p)
    @test (A + t)(p) ≈ A(p .+ t)
    @test (A - t)(p) ≈ A(p .- t)
    @test (α*A)(p) ≈ α.*A(p)
    @test (A*α)(p) ≈ A(α.*p)
end

@testset "Regions" begin
    stp, xc, yc = 1.2, -2.7, 3.1
    n1, n2 = 27, 52
    c = 0.7, -1.1
    grdxmin = xc - (n1 - 1)*stp/2
    grdxmax = xc + (n1 - 1)*stp/2
    grdymin = yc - (n2 - 1)*stp/2
    grdymax = yc + (n2 - 1)*stp/2
    boxxmin = xc - n1*stp/2
    boxxmax = xc + n1*stp/2
    boxymin = yc - n2*stp/2
    boxymax = yc + n2*stp/2
    X = linspace(grdxmin, grdxmax, n1)
    Y = linspace(grdymin, grdymax, n2)
    R1 = Region(X, Y)
    @test center(R1) ≈ (xc, yc)
    @test step(R1) ≈ stp
    @test size(R1) ≈ (n1, n2)
    @test extrema(R1) ≈ (grdxmin, grdxmax, grdymin, grdymax)
    @test boundingbox(R1) ≈ (boxxmin, boxxmax, boxymin, boxymax)
    R2 = recenter(R1)
    @test center(R2) ≈ (0,0)
    R2 = recenter(R1, c)
    @test center(R2) ≈ c
    @test center(R2 - c) ≈ (0,0)
end

end


text:
module IncrementalInference

@info "Multithreaded  convolutions possible, Threads.nthreads()=$(Threads.nthreads()).  See `addFactor!(.;threadmodel=MultiThreaded)`."

using Distributed
using Requires
using Reexport

@reexport using Distributions
@reexport using KernelDensityEstimate
@reexport using ApproxManifoldProducts
@reexport using Graphs
@reexport using LinearAlgebra

using
  Dates,
  DistributedFactorGraphs,
  Statistics,
  Random,
  NLsolve,
  StatsBase,
  JLD2,
  FileIO,
  ProgressMeter,
  DocStringExtensions,
  FunctionalStateMachine,
  Optim, # might be deprecated in favor for only NLsolve dependency
  JSON2

using Logging

const KDE = KernelDensityEstimate
const AMP = ApproxManifoldProducts
const DFG = DistributedFactorGraphs
const FSM = FunctionalStateMachine

import Base: convert
# import HDF5: root
import Distributions: sample
import Random: rand, rand!
import KernelDensityEstimate: getBW
import ApproxManifoldProducts: kde!, manikde!
import DistributedFactorGraphs: addVariable!, addFactor!, ls, lsf, isInitialized, hasOrphans, compare, compareAllSpecial

# missing exports
import DistributedFactorGraphs: PackedFunctionNodeData, FunctionNodeData

# TODO temporary for initial version of on-manifold products
KDE.setForceEvalDirect!(true)

export
  KDE,
  AMP,
  DFG,
  FSM,

  # DFG SpecialDefinitions
  AbstractDFG,
  hasVariable,
  getSolverParams,

  *,
  notifyCSMCondition,
  CSMHistory,
  getTreeCliqsSolverHistories,
  assignTreeHistory!,

  updateCliqSolvableDims!,
  fetchCliqSolvableDims,
  BayesTreeNodeData,

  # state machine methods
  StateMachine,
  exitStateMachine,
  getCliqSolveHistory,
  filterHistAllToArray,
  cliqHistFilterTransitions,
  printCliqSummary,
  printCliqHistorySummary,
  printGraphSummary,
  printSummary,
  getGraphFromHistory,
  getCliqSubgraphFromHistory,
  sandboxStateMachineStep,
  sandboxCliqResolveStep,
  # draw and animate state machine
  getStateLabel,
  histStateMachineTransitions,
  histGraphStateMachineTransitions,
  drawStateTransitionStep,
  drawStateMachineHistory,
  animateStateMachineHistoryByTime,
  animateStateMachineHistoryByTimeCompound,
  animateCliqStateMachines,
  csmAnimate,
  makeCsmMovie,
  lockUpStatus!,
  unlockUpStatus!,
  lockDwnStatus!,
  unlockDwnStatus!,
  getSiblingsDelayOrder,
  areSiblingsRemaingNeedDownOnly,

  # general types for softtyping of variable nodes
  InferenceVariable,
  ContinuousScalar,
  ContinuousMultivariate,
  SamplableBelief,
  Prior,
  PackedPrior,
  MsgPrior,
  PackedMsgPrior,
  PartialPrior,
  PackedPartialPrior,
  LinearConditional,
  PackedLinearConditional,
  MixturePrior,
  PackedMixturePrior,
  MixtureLinearConditional,
  PackedMixtureLinearConditional,

  ls2,
  lsRear,
  # from DFG
  ls,
  lsf,
  getVariableIds,
  calcVariablePPE,
  sortVarNested,
  hasOrphans,
  drawCopyFG,
  isVariable,
  isFactor,
  # from dfg
  getfnctype,
  getFactorType,
  getSofttype,
  getVariableType,
  lsfPriors,
  isPrior,
  lsTypes,
  lsfTypes,

  # using either dictionary or cloudgraphs
  # VariableNodeData,
  # PackedVariableNodeData,
  FactorMetadata,
  getpackedtype,
  encodePackedType,
  FunctionNodeData,
  PackedFunctionNodeData, # moved to DFG
  encodePackedType,
  decodePackedType,
  normalfromstring,
  categoricalfromstring,
  extractdistribution,

  FactorGraph,
  SolverParams,
  addNode!,
  addVariable!,
  deleteVariable!,
  addFactor!,
  deleteFactor!,
  addMsgFactors!,
  deleteMsgFactors!,
  factorCanInitFromOtherVars,
  doautoinit!,
  manualinit!,
  initVariable!,
  asyncTreeInferUp!,
  initInferTreeUp!,
  solveCliqWithStateMachine!,
  resetData!,
  resetTreeCliquesForUpSolve!,
  resetFactorGraphNewTree!,
  setVariableInitialized!,
  setVariableInferDim!,
  resetVariable!,
  resetCliqSolve!,
  getFactor,
  getFactorDim,
  getVariableDim,
  getVariableInferredDim,
  getVariableInferredDimFraction,
  getVariablePotentialDims,
  getVariableSolvableDim,
  getFactorSolvableDim,
  getFactorInferFraction,
  getCliqSiblingsPriorityInitOrder,
  isCliqFullDim,
  hasFactor,
  getVariable,
  # getVert, # deprecated use DFG.getVariable getFactor instead
  getData,
  setData!,
  getManifolds,
  getVarNode,
  getVal,
  getBW,
  setVal!,
  getNumPts,
  getBWVal,
  setBW!,
  setValKDE!,
  buildCliqSubgraphUp,
  buildCliqSubgraphDown,
  setCliqUpInitMsgs!,
  cliqInitSolveUpByStateMachine!,

  # state machine functions
  finishCliqSolveCheck_StateMachine,
  doCliqInferAttempt_StateMachine,
  determineCliqNeedDownMsg_StateMachine,
  blockUntilChildrenStatus_StateMachine,
  blockUntilSiblingsStatus_StateMachine,
  doesCliqNeeddownmsg_StateMachine,
  slowCliqIfChildrenNotUpsolved_StateMachine,
  whileCliqNotSolved_StateMachine,
  buildCliqSubgraph_StateMachine,
  isCliqUpSolved_StateMachine,
  determineAllChildrenNeedDownMsg_StateMachine,
  testCliqCanRecycled_StateMachine,
  buildCliqSubgraphForDown_StateMachine,

  #
  isPartial,
  isInitialized,
  isTreeSolved,
  isUpInferenceComplete,
  isCliqInitialized,
  isCliqReadyInferenceUp,
  isCliqUpSolved,
  areCliqVariablesAllInitialized,
  areCliqChildrenNeedDownMsg,
  areCliqChildrenAllUpSolved,
  ensureAllInitialized!,
  doCliqAutoInitUpPart1!,
  doCliqAutoInitUpPart2!,
  doCliqInitDown!,
  cycleInitByVarOrder!,
  doCliqUpSolve!,
  getCliqInitUpMsgs,
  getCliqStatus,
  setCliqStatus!,
  getSolveCondition,
  prepCliqInitMsgsUp,
  prepCliqInitMsgsDown!,
  updateFullVert!,
  getOutNeighbors,
  BayesTree,
  EasyMessage,
  NBPMessage,
  FullExploreTreeType,
  ExploreTreeType,
  FactorGraph,
  initfg,
  buildSubgraphFromLabels,
  subgraphShortestPath,
  transferUpdateSubGraph!,
  getEliminationOrder,
  buildBayesNet!,
  emptyBayesTree,
  buildTree!,
  buildTreeFromOrdering!,
  resetBuildTreeFromOrder!,
  prepBatchTree!,
  wipeBuildNewTree!,
  buildCliquePotentials,
  hasCliq,
  getCliq,
  whichCliq,
  getCliqDepth,
  getTreeAllFrontalSyms,
  getTreeCliqUpMsgsAll,
  stackCliqUpMsgsByVariable,
  getCliqChildMsgsUp,
  getCliqParentMsgDown,
  getCliqDownMsgsAfterDownSolve,
  isReadyCliqInferenceUp,
  childCliqs,
  getChildren,
  parentCliq,
  getParent,
  getCliqSiblings,
  getKDE,
  getVertKDE,
  initializeNode!,
  CliqStateMachineContainer,
  batchSolve!,
  solveTree!,
  solveCliq!,
  fifoFreeze!,
  getCurrentWorkspaceFactors,
  getCurrentWorkspaceVariables,

  # temp const types TODO
  TempBeliefMsg,
  TempUpMsgPlotting,

  #functors need
  getSample,
  freshSamples!,

  #Visualization
  writeGraphPdf, # deprecated, but first move code to drawGraph before deleting
  drawGraph,
  drawGraphCliq,
  drawCliqSubgraphUpMocking,
  drawTree,
  printgraphmax,
  # allnums,
  # isnestednum,
  # sortnestedperm,

  # Bayes (Junction) Tree
  evalPotential,
  evalFactor2,
  approxConv,
  approxConvBinary,

  # more debugging tools
  localProduct,
  treeProductUp,
  approxCliqMarginalUp!,
  dontMarginalizeVariablesAll!,
  unfreezeVariablesAll!,
  resetVariableAllInitializations!,
  isMarginalized,
  isMultihypo,
  getMultihypoDistribution,
  getHypothesesVectors,
  isCliqMarginalizedFromVars,
  isCliqParentNeedDownMsg,
  setCliqAsMarginalized!,
  updateTreeCliquesAsMarginalizedFromVars!,

  # weiged sampling
  AliasingScalarSampler,
  rand!,
  rand,
  fastnorm,

  # new wrapper (experimental)
  CommonConvWrapper,

  # solve inference
  inferOverTree!,
  inferOverTreeR!,
  inferOverTreeIterative!,

  #development interface
  getTreeCliqSolveOrderUp,
  getCliqOrderUpSolve,
  getCliqVarInitOrderUp,
  getCliqInitVarOrderDown,
  getCliqStatusUp,
  blockCliqUntilChildrenHaveUpStatus,
  blockCliqSiblingsParentNeedDown,
  getCliqNumAssocFactorsPerVar,
  upMsgPassingRecursive,
  downMsgPassingRecursive,

  upMsgPassingIterative!,
  downMsgPassingIterative!,

  # Inference types
  InferenceType,
  PackedInferenceType,
  Singleton,
  Pairwise,
  # introduced for approximate convolution operations
  setThreadModel!,
  SingleThreaded,
  MultiThreaded,

  # functor abstracts
  FunctorInferenceType,
  FunctorPairwise,
  FunctorPairwiseMinimize,
  FunctorSingleton,
  # FunctorPartialSingleton,
  FunctorPairwiseNH,
  FunctorSingletonNH,

  # Solving utils
  findRelatedFromPotential,
  shuffleXAltD!,
  numericRoot,
  numericRootGenericRandomized,
  numericRootGenericRandomizedFnc,
  numericRootGenericRandomizedFnc!,
  solveFactorMeasurements,

  # user functions
  proposalbeliefs,
  predictbelief,
  getCliq,
  getCliqMat,
  getCliqAssocMat,
  getCliqMsgMat,
  getCliqFrontalVarIds,
  getFrontals,                     # duplicate
  getCliqSeparatorVarIds,
  getCliqAllVarIds,
  getCliqVarIdsAll,
  getCliqVars,
  getCliqAllVarSyms,
  getCliqVarIdsPriors,
  getCliqVarSingletons,
  getCliqAllFactIds,
  getCliqFactorIdsAll,
  areCliqVariablesAllMarginalized,
  setTreeCliquesMarginalized!,

  # generic marginal used during elimitation game
  GenericMarginal,
  PackedGenericMarginal,

  uppA,
  convert,
  extractdistribution,

  # factor graph operating system utils (fgos)
  convert2packedfunctionnode,
  encodefg,
  decodefg,
  savejld,
  loadjld,
  landmarks,
  setCliqDrawColor,
  fetchCliqTaskHistoryAll!,
  fetchAssignTaskHistoryAll!,

  # Temp placeholder for evaluating string types to real types
  _evalType,
  saveDFG,
  loadDFG,
  rebuildFactorMetadata!,

  setUpMsg!,
  upMsg,
  getDwnMsgs,
  getUpMsgs,
  setDwnMsg!,
  dwnMsg,
  getDwnMsgs,
  getCliqMsgsUp,
  getCliqMsgsDown,
  getCliqVarSolveOrderUp,

  getSym,
  doCliqInferenceUp!,
  getFactorsAmongVariablesOnly,
  setfreeze!,

  #internal dev functions for recycling cliques on tree
  attemptTreeSimilarClique,
  getCliqSiblingsPartialNeeds,

  # some utils
  compare,
  compareAllSpecial,
  getIdx,
  showFactor,
  showVariable,
  getMeasurements,
  findFactorsBetweenFrom,
  addDownVariableFactors!,
  getDimension,

  # For 1D example,

  # TODO rename to L2 distance
  Ranged,
  PackedRanged,

  # development
  # dev exports
  addGraphsVert!,
  makeAddEdge!,
  shuffleXAltD,
  reshapeVec2Mat, # TODO deprecate

  # OBSOLETE TODO REMOVE #TODO TODO
  subGraphFromVerts,
  getMaxVertId,
  dlapi,
  localapi,
  showcurrentdlapi,
  setdatalayerAPI!,
  DataLayerAPI

# TODO should be deprecated
const NothingUnion{T} = Union{Nothing, T}

# non-free, but not currently use!
include("ccolamd.jl")

# regular
include("FactorGraphTypes.jl")

# const InMemDFGType = DFG.LightDFG{SolverParams} #swap out default in v0.8.0/v0.9.0?
const InMemDFGType = DFG.GraphsDFG{SolverParams} # JT TODO move to somewhere more fitting?

include("BeliefTypes.jl")
include("AliasScalarSampling.jl")
include("DefaultNodeTypes.jl")
include("FactorGraph01.jl")
include("SubGraphFunctions.jl")
include("SerializingDistributions.jl")
include("DispatchPackedConversions.jl")
include("FGOSUtils.jl")
include("CompareUtils.jl")

include("JunctionTreeTypes.jl")
include("JunctionTree.jl")
include("TreeBasedInitialization.jl")
include("GraphConstraintTypes.jl")
include("SolverUtilities.jl")
include("ExplicitDiscreteMarginalizations.jl")
include("InferDimensionUtils.jl")
include("ApproxConv.jl")
include("SolveTree01.jl")
include("CliqStateMachine.jl")
include("CliqStateMachineUtils.jl")

# special variables and factors, see RoME.jl for more examples
include("Variables/Sphere1D.jl")
include("Factors/Sphere1D.jl")

include("AdditionalUtils.jl")
include("SolverAPI.jl")

include("Deprecated.jl")

exportimg(pl) = error("Please do `using Gadfly` before IncrementalInference is used to allow image export.")
function __init__()
    @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" begin
      @info "Defining spyCliqMat(..) for visualizing association matrix of a clique in the Bayes (Juntion) tree"

      exportimg(pl) = Gadfly.PNG(pl)

      export spyCliqMat

      """
          $SIGNATURES

      Draw the clique association matrix, with keyword arguments for more or less console print outs.

      Notes
      * Columns are variables, rows are factors.
      * Drawn from up message passing perspective.
      * Blue color implies no factor association.
      * Frontal, separator, and upmessages are all drawn at different intensity of red.
      * Downward messages not shown, as they would just be singletons of the full separator set.
      """
      function spyCliqMat(cliq::Graphs.ExVertex; showmsg=true, suppressprint::Bool=false)
        mat = deepcopy(getCliqMat(cliq, showmsg=showmsg))
        # TODO -- add improved visualization here, iter vs skip
        mat = map(Float64, mat)*2.0.-1.0
        numlcl = size(getCliqAssocMat(cliq),1)
        mat[(numlcl+1):end,:] *= 0.9
        mat[(numlcl+1):end,:] .-= 0.1
        numfrtl1 = floor(Int,length(getData(cliq).frontalIDs) + 1)
        mat[:,numfrtl1:end] *= 0.9
        mat[:,numfrtl1:end] .-= 0.1
        if !suppressprint
          @show getData(cliq).itervarIDs
          @show getData(cliq).directvarIDs
          @show getData(cliq).msgskipIDs
          @show getData(cliq).directFrtlMsgIDs
          @show getData(cliq).directPriorMsgIDs
        end
        sp = Gadfly.spy(mat)
        push!(sp.guides, Gadfly.Guide.title("$(cliq.attributes["label"]) || $(cliq.attributes["data"].frontalIDs) :$(cliq.attributes["data"].conditIDs)"))
        push!(sp.guides, Gadfly.Guide.xlabel("fmcmcs $(cliq.attributes["data"].itervarIDs)"))
        push!(sp.guides, Gadfly.Guide.ylabel("lcl=$(numlcl) || msg=$(size(getCliqMsgMat(cliq),1))" ))
        return sp
      end
      function spyCliqMat(bt::BayesTree, lbl::Symbol; showmsg=true, suppressprint::Bool=false)
        spyCliqMat(whichCliq(bt,lbl), showmsg=showmsg, suppressprint=suppressprint)
      end
    end
end

# Old code that might be used again
# function getType(typestring::AS) where {AS <: AbstractString}
#  # eval(Meta.parse(typestring))()
#  # getfield(Main, Symbol(typestring))
#  getfield(@__MODULE__, Symbol(typestring))
# end

export setSerializationNamespace!, getSerializationModule, getSerializationModules

end


text:
using NeighbourLists
using JuLIP
using BenchmarkTools
using PyCall
using ASE

include("profile_pairlist.jl")
include("profile_nbody.jl")
include("profile_forces.jl")


text:
## %PREDICTMD_GENERATED_BY%

import PredictMDExtra
PredictMDExtra.import_all()

import PredictMD
PredictMD.import_all()

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "boston_housing", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin model output code

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

linear_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "linear_regression.jld2",
    )
random_forest_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "random_forest_regression.jld2",
    )
knet_mlp_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "knet_mlp_regression.jld2",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
linear_regression = nothing
Test.@test isnothing(linear_regression)

random_forest_regression = nothing
Test.@test isnothing(random_forest_regression)

knet_mlp_regression = nothing
Test.@test isnothing(knet_mlp_regression)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

linear_regression =
    PredictMD.load_model(linear_regression_filename)
random_forest_regression =
    PredictMD.load_model(random_forest_regression_filename)
knet_mlp_regression =
    PredictMD.load_model(knet_mlp_regression_filename)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

PredictMD.predict(linear_regression,training_features_df,)
PredictMD.predict(random_forest_regression,training_features_df,)
PredictMD.predict(knet_mlp_regression,training_features_df,)

PredictMD.predict(linear_regression,testing_features_df,)
PredictMD.predict(random_forest_regression,testing_features_df,)
PredictMD.predict(knet_mlp_regression,testing_features_df,)

### End model output code

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "boston_housing", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
immutable Laplace <: ContinuousUnivariateDistribution
    μ::Float64
    β::Float64

    function Laplace(μ::Real, β::Real)
        β > zero(β) || error("Laplace's scale must be positive")
        @compat new(Float64(μ), Float64(β))
    end

    @compat Laplace(μ::Real) = new(Float64(μ), 1.0)
    Laplace() = new(0.0, 1.0)
end

typealias Biexponential Laplace

@distr_support Laplace -Inf Inf

#### Parameters

location(d::Laplace) = d.μ
scale(d::Laplace) = d.β
params(d::Laplace) = (d.μ, d.β)

#### Statistics

mean(d::Laplace) = d.μ
median(d::Laplace) = d.μ
mode(d::Laplace) = d.μ

var(d::Laplace) = 2.0 * d.β^2
std(d::Laplace) = sqrt2 * d.β
skewness(d::Laplace) = 0.0
kurtosis(d::Laplace) = 3.0

entropy(d::Laplace) = log(2.0 * d.β) + 1.0

#### Evaluations

zval(d::Laplace, x::Float64) = (x - d.μ) / d.β
xval(d::Laplace, z::Float64) = d.μ + z * d.β

pdf(d::Laplace, x::Float64) = 0.5 * exp(-abs(zval(d, x))) / scale(d)
logpdf(d::Laplace, x::Float64) = - (abs(zval(d, x)) + log(2.0 * scale(d)))

cdf(d::Laplace, x::Float64) = (z = zval(d, x); z < 0.0 ? 0.5 * exp(z) : 1.0 - 0.5 * exp(-z))
ccdf(d::Laplace, x::Float64) = (z = zval(d, x); z > 0.0 ? 0.5 * exp(-z) : 1.0 - 0.5 * exp(z))
logcdf(d::Laplace, x::Float64) = (z = zval(d, x); z < 0.0 ? loghalf + z : loghalf + log2mexp(-z))
logccdf(d::Laplace, x::Float64) = (z = zval(d, x); z > 0.0 ? loghalf - z : loghalf + log2mexp(z))

quantile(d::Laplace, p::Float64) = p < 0.5 ? xval(d, log(2.0 * p)) : xval(d, -log(2.0 * (1.0 - p)))
cquantile(d::Laplace, p::Float64) = p > 0.5 ? xval(d, log(2.0 * (1.0 - p))) : xval(d, -log(2.0 * p))
invlogcdf(d::Laplace, lp::Float64) = lp < loghalf ? xval(d, logtwo + lp) : xval(d, -(logtwo + log1mexp(lp)))
invlogccdf(d::Laplace, lp::Float64) = lp > loghalf ? xval(d, logtwo + log1mexp(lp)) : xval(d, -(logtwo + lp))

function gradlogpdf(d::Laplace, x::Float64)
    μ, β = params(d)
    x == μ && error("Gradient is undefined at the location point")
    g = 1.0 / β
    x > μ ? -g : g
end

function mgf(d::Laplace, t::Real)
    st = d.β * t
    exp(t * d.μ) / ((1.0 - st) * (1.0 + st))
end
function cf(d::Laplace, t::Real)
    st = d.β * t
    cis(t * d.μ) / (1+st*st)
end

#### Sampling

function rand(d::Laplace)
    z = randexp()
    rand(Bool) ? d.μ + d.β * z : d.μ - d.β * z
end

#### Fitting

function fit_mle(::Type{Laplace}, x::Array)
    a = median(x)
    Laplace(a, mad(x, a))
end


text:
using Oceananigans.Grids: scalar_summary

prettysummary(x, args...) = summary(x)

function prettysummary(f::Function, showmethods=true)
    ft = typeof(f)
    mt = ft.name.mt
    name = mt.name
    n = length(methods(f))
    m = n==1 ? "method" : "methods"
    sname = string(name)
    isself = isdefined(ft.name.module, name) && ft == typeof(getfield(ft.name.module, name))
    ns = (isself || '#' in sname) ? sname : string("(::", ft, ")")
    if showmethods
        return string(ns, " (", "generic function", " with $n $m)")
    else
        return string(ns)
    end
end

prettysummary(x::AbstractFloat, args...) = scalar_summary(x)
prettysummary(x::Int, args...) = string(x)

# This is very important
function prettysummary(nt::NamedTuple, args...)
    n = nfields(nt)

    if n == 0
        return "NamedTuple()"
    else
        str = "("
        for i = 1:n
            f = nt[i]
            str = string(str, fieldname(typeof(nt), i), "=", prettysummary(getfield(nt, i)))
            if n == 1
                str = string(str, ",")
            elseif i < n
                str = string(str, ", ")
            end
        end
    end

    return string(str, ")")
end

function prettykeys(t)
    names = collect(keys(t))
    length(names) == 0 && return "()"
    length(names) == 1 && return string(first(names))
    return string("(", (string(n, ", ") for n in names[1:end-1])..., last(names), ")")
end


text:
function force_bond!(sys::System,bdinter::NoBond)
end

function force_bond!(sys::System,bdinter::Harmonic)
    force_bond_atomic!.(Ref(sys.atoms.forces),sys.atoms.bonds,Ref(sys.atoms.coords),bdinter,sys.pbc.Lx,sys.pbc.Ly,sys.pbc.Lz)
    #@show sum(sys.atoms.forces)
end

function force_bond_atomic!(forces::Array{Array{Float64,1}},
                            bond::Array{Int64,1},
                            coords::Array{Array{Float64,1}},
                            bdinter::Harmonic,lx::Float64,ly::Float64,lz::Float64)

    delta_x = coords[bond[2]][1]- coords[bond[3]][1]
    delta_y = coords[bond[2]][2]- coords[bond[3]][2]
    delta_z = coords[bond[2]][3]- coords[bond[3]][3]

    # consider pbc
    delta_x -= lx * RInt( delta_x / lx)
    delta_y -= ly * RInt( delta_y / ly)
    delta_z -= lz * RInt( delta_z / lz)

    r=(delta_x^2+delta_y^2+delta_z^2)^0.5

    κ=bdinter.k[bond[1]] #κ
    l0=bdinter.l0[bond[1]] #r_0

    force_=κ*(l0/r-1.0)
    force_divr = force_/r
    #@show r,bond
    forces[bond[2]][1] += delta_x * force_divr
    forces[bond[2]][2] += delta_y * force_divr
    forces[bond[2]][3] += delta_z * force_divr

    forces[bond[3]][1] += -delta_x * force_divr
    forces[bond[3]][2] += -delta_y * force_divr
    forces[bond[3]][3] += -delta_z * force_divr
    #@show bond[2],bond[3],forces[bond[2]][1],forces[bond[3]][1],force_,force_divr,r,delta_x,delta_y,delta_z,force_divr
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

using Plots

# +
plotly()

a=0.5

p = plot3d(xlabel="x",ylabel="y",zlabel="z",xlim=(-1.5,1.5),ylim=(-1.5,1.5))
# Draw slice along xy-plane
us = 0:0.05π:2π
ts = -1:0.05:1

for t in -1:0.05:1
    x=@. (t^2+a)*cos(us)
    y=@. (t^2+a)*sin(us)
    z=t*ones(length(x))
    plot3d!(p, x,y,z,label="",color=:blue)
end

for u in us
    x = @. (ts^2+a)*cos(u)
    y = @. (ts^2+a)*sin(u)
    z = ts
    plot3d!(p, x,y,z,label="",color=:blue)
end

display(p)

# +
gr()

as = hcat(reverse(0.0:0.01:0.5), 0.0:0.01:0.5)

anim=@animate for a in as
    IJulia.clear_output(true)
    p = plot3d(xlabel="x",ylabel="y",zlabel="z",xlim=(-1.5,1.5),ylim=(-1.5,1.5))
    # Draw slice along xy-plane
    us = 0:0.05π:2π
    ts = -1:0.05:1

    for t in -1:0.05:1
        x=@. (t^2+a)*cos(us)
        y=@. (t^2+a)*sin(us)
        z=t*ones(length(x))
        plot3d!(p, x,y,z,label="",color=:blue)
    end

    for u in us
        x = @. (ts^2+a)*cos(u)
        y = @. (ts^2+a)*sin(u)
        z = ts
        plot3d!(p, x,y,z,label="",color=:blue)
    end
    p |> display
end
# -

gif(anim)

# # Torus

# +
plotly()
a=3.

p = plot3d()
N = 45
u = range(0,2π,length=N)
for t in range(0,2π,length=N)
    x = @. cos(t) * (a + cos(u))
    y = @. sin(t) * (a + cos(u))
    z = @. sin(u)
    plot3d!(p, x, y, z, color=:blue, legend=false)
end

t=range(0,2π,length=N)
for u in range(0,2π,length=N)
    x = @. cos(t) * (a + cos(u))
    y = @. sin(t) * (a + cos(u))
    z = sin(u) * ones(length(x))
    plot3d!(p, x, y, z, color=:blue, legend=false)
end

p


text:
function rotate!(cur_player::AbstractPlayer, cur_piece::AbstractPiece, cur_grid::AbstractGrid, cur_direction::Integer)
  can_move(cur_player) || return false

  ( cur_piece.name == 'o' ) && return false

  initial_rotation = cur_piece.rotation

  cur_piece.rotation += cur_direction

  while cur_piece.rotation <= 0
    cur_piece.rotation += 4
  end

  cur_piece.rotation %= 4

  cur_blocks = map(
    cur_block -> calc_block_coords(cur_block, cur_direction),
    cur_piece.blocks
  )

  is_legal_move = check(cur_blocks, cur_grid)

  cur_x_offset = 0
  cur_y_offset = 0

  if !is_legal_move

    cur_offsets =
      get_rotate_offsets(
        initial_rotation,
        cur_piece.rotation,
        iseven(cur_piece.width)
      )

    for cur_offset in cur_offsets
      is_legal_move = check(
        cur_blocks,
        cur_grid,
        cur_offset
      )

      is_legal_move || continue

      cur_x_offset, cur_y_offset = cur_offset
      break
    end

  end

  is_legal_move || return false

  cur_piece.i += cur_y_offset
  cur_piece.j += cur_x_offset

  for cur_block in cur_piece.blocks
    cur_block.x, cur_block.y =
      ( cur_block.y * +sign(cur_direction) ),
      ( cur_block.x * -sign(cur_direction) )
  end

  reset_shadow!(cur_player)

  move!(cur_player)

  true
end

function rotate!(cur_player::AbstractPlayer, cur_piece::Nullable{AbstractPiece}, cur_grid::AbstractGrid, cur_direction::Integer)
  false
end


text:
const dust_regular = (5, 20, 100, 400)
const dust_golden = (50, 100, 400, 1600)

const craft_regular = (40, 100, 400, 1600)

function hssim(
        nreps::Int,
        decklist::DecklistArray,
        sets::Vector{Symbol},
        parameters::HearthstoneSimParameters = HearthstoneSimParameters()
    )

    sets_in_decklist = mapreduce(i -> collect(keys(decklist[i])), union, 1:4)

    @assert issubset(sets, sets_in_decklist)

    collection = Collection(sets, parameters.rarity_by_set)
    fresh_collection = deepcopy(collection)

    account_state::AccountState = AccountState()
    account_state.bonus_packs[:EXPERT1] = 7
    for s in parameters.sets_with_bonus
        account_state.bonus_packs[s] += 3
    end
    account_state.dust += parameters.bonus_dust
    if parameters.welcome_bundle
        account_state.dust += 400
        account_state.bonus_packs[:EXPERT1] += 10
    end
    fresh_account = deepcopy(account_state)

    ndeck = length(decklist)

    output = SimOutput(sets, nreps, ndeck)

    pack_table = Dict(
        map(x -> (x, packtablegen(ncard_db[x], parameters.rarity_probs)), sets)
    )

    pack_contents = map(i -> (0,0), 1:5)

    deckiter = collect(1:length(decklist))
    setiter = deepcopy(sets)

    for rep = 1:nreps

        if parameters.startfresh
            account_state = deepcopy(fresh_account)
            collection = deepcopy(fresh_collection)
        end

        for deckidx in deckiter
            for set in setiter
                sim_state::SimState = SimState(rep = rep, deck = deckidx, set = set, npacks = 0)

                subtarget = SingleSetDeckArray(decklist, set, deckidx)

                while missing_cards(collection, subtarget)

                    if account_state.dust >= remaining_dust(collection, subtarget)
                        spend_dust!(collection, subtarget, account_state, sim_state)
                        break
                    end

                    if account_state.bonus_packs[set] == 0
                        sim_state.npacks += 1
                    else
                        account_state.bonus_packs[set] -= 1
                    end

                    open_pack!(pack_contents,pack_table, collection, decklist, parameters,
                    account_state, sim_state, output)

                end

                output.total_packs_for_deck[set][deckidx, rep] = sim_state.npacks
                output.total_packs[set][rep] += sim_state.npacks

            end
        end

    end

    output
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
# ---
# title: 1388. Pizza With 3n Slices
# id: problem1388
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/pizza-with-3n-slices/description/>
# hidden: true
# ---
#
# There is a pizza with 3n slices of varying size, you and your friends will
# take slices of pizza as follows:
#
#   * You will pick **any** pizza slice.
#   * Your friend Alice will pick next slice in anti clockwise direction of your pick.
#   * Your friend Bob will pick next slice in clockwise direction of your pick.
#   * Repeat until there are no more slices of pizzas.
#
# Sizes of Pizza slices is represented by circular array `slices` in clockwise
# direction.
#
# Return the maximum possible sum of slice sizes which you can have.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/18/sample_3_1723.png)
#
#
#
#     Input: slices = [1,2,3,4,5,6]
#     Output: 10
#     Explanation: Pick pizza slice of size 4, Alice and Bob will pick slices with size 3 and 5 respectively. Then Pick slices with size 6, finally Alice and Bob will pick slice of size 2 and 1 respectively. Total = 4 + 6.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/02/18/sample_4_1723.png)**
#
#
#
#     Input: slices = [8,9,8,6,1,1]
#     Output: 16
#     Output: Pick pizza slice of size 8 in each turn. If you pick slice with size 9 your partners will pick slices of size 8.
#
#
# **Example 3:**
#
#
#
#     Input: slices = [4,1,2,5,8,3,1,9,7]
#     Output: 21
#
#
# **Example 4:**
#
#
#
#     Input: slices = [3,1,2]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= slices.length <= 500`
#   * `slices.length % 3 == 0`
#   * `1 <= slices[i] <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
using PyPlot
using Statistics

Statistics.mean(x::Real,y::Real) = (x+y)/2

x = [0.5,1.0,2.0]
xint = [mean(x[i],x[i+1]) for i in 1:length(x)-1]

fig,(ax1,ax2) = subplots(2,1,figsize=(5,2.2),sharex=true)

for i in 1:length(x)
    ax1.scatter(x[i],0,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax1.plot([0,xint[1]],[0,0],"C0")
ax1.plot([xint[1],xint[2]],[0,0],"C1")
ax1.plot([xint[2],3],[0,0],"C2")
ax1.set_yticks([])
ax1.set_frame_on(false)
ax1.set_title("Round nearest",loc="left")

ax1.plot([1.2,1.2],[-0.2,0.2],"k--",lw=0.5)
ax1.text(1.22,-0.2,L"x is round to $x_2$")

ax2.plot([-0.5,0,x[1]],[0,1,0],"C5")
ax2.plot([0,x[1],x[2]],[0,1,0],"C0")
ax2.plot([x[1],x[2],x[3]],[0,1,0],"C1")
ax2.plot([x[2],x[3],4],[0,1,0],"C2")
ax2.plot([x[3],4,5],[0,1,0],"C4")

for i in 1:length(x)
    ax2.scatter(x[i],1,20,color="C$(i-1)",edgecolor="k",zorder=10)
end

ax2.plot([1.2,1.2],[0,1],"k--",lw=0.5)
ax2.text(1.18,-0.3,L"x = ")
ax2.text(1.31,-0.15,L"$x_2$ at 80%")
ax2.text(1.31,-0.45,L"$x_3$ at 20% chance.")

ax2.set_xlim(x[1]-0.2,x[3]+0.2)
ax2.set_yticks([0,1])
ax2.set_yticklabels([L"0\%",L"100\%"])
ax2.set_xticks(x)
ax2.set_xticklabels([L"$x_1$",L"$x_2$",L"$x_3$"])
ax2.set_frame_on(false)
ax2.set_title("Stochastic rounding",loc="left")

ax1.set_title("a",loc="right",fontweight="bold")
ax2.set_title("b",loc="right",fontweight="bold")

tight_layout()
savefig("figs/schematic.pdf")
savefig("figs/schematic.png",dpi=200)
close(fig)


text:
import DataFrames.DataFrame

"""
Sets the data source using a DataFrame
`data_values(df::DataFrame)`

Adds all the columns of `df` to the data field.
"""
function data_values(df::DataFrame)
  vs = Any[]
  for n in names(df)
      push!(vs, (n, collect(df[n])) )
  end
  data_values(;vs...)
end


text:
using Oceananigans.Architectures: device_event
using Oceananigans.BoundaryConditions: OBC, CBC
using Oceananigans.Grids: parent_index_range, index_range_offset, default_indices, all_indices

using Adapt
using KernelAbstractions: @kernel, @index
using Base: @propagate_inbounds

import Oceananigans.BoundaryConditions: fill_halo_regions!
import Statistics: norm, mean, mean!

#####
##### The bees knees
#####

struct Field{LX, LY, LZ, O, G, I, D, T, B, S, F} <: AbstractField{LX, LY, LZ, G, T, 3}
    grid :: G
    data :: D
    boundary_conditions :: B
    indices :: I
    operand :: O
    status :: S
    boundary_buffers :: F

    # Inner constructor that does not validate _anything_!
    function Field{LX, LY, LZ}(grid::G, data::D, bcs::B, indices::I, op::O, status::S, buffers::F) where {LX, LY, LZ, G, D, B, O, S, I, F}
        T = eltype(data)
        return new{LX, LY, LZ, O, G, I, D, T, B, S, F}(grid, data, bcs, indices, op, status, buffers)
    end
end

#####
##### Constructor utilities
#####

function validate_field_data(loc, data, grid, indices)
    Fx, Fy, Fz = total_size(loc, grid, indices)

    if size(data) != (Fx, Fy, Fz)
        LX, LY, LZ = loc
        e = "Cannot construct field at ($LX, $LY, $LZ) with size(data)=$(size(data)). " *
            "`data` must have size ($Fx, $Fy, $Fz)."
        throw(ArgumentError(e))
    end

    return nothing
end

validate_boundary_condition_location(bc, ::Center, side) = nothing                        # anything goes for centers
validate_boundary_condition_location(::Union{OBC, Nothing, CBC}, ::Face, side) = nothing  # only open, connected or nothing on faces
validate_boundary_condition_location(::Nothing, ::Nothing, side) = nothing                # its nothing or nothing
validate_boundary_condition_location(bc, loc, side) = # everything else is wrong!
    throw(ArgumentError("Cannot specify $side boundary condition $bc on a field at $(loc)!"))

validate_boundary_conditions(loc, grid, ::Missing) = nothing
validate_boundary_conditions(loc, grid, ::Nothing) = nothing

function validate_boundary_conditions(loc, grid, bcs)
    sides = (:east, :west, :north, :south, :bottom, :top)
    directions = (1, 1, 2, 2, 3, 3)

    for (side, dir) in zip(sides, directions)
        topo = topology(grid, dir)()
        ℓ = loc[dir]()
        bc = getproperty(bcs, side)

        # Check that boundary condition jives with the grid topology
        validate_boundary_condition_topology(bc, topo, side)

        # Check that boundary condition is valid given field location
        topo isa Bounded && validate_boundary_condition_location(bc, ℓ, side)

        # Check that boundary condition arrays, if used, are on the right architecture
        validate_boundary_condition_architecture(bc, architecture(grid), side)
    end

    return nothing
end

function validate_index(idx, loc, topo, N, H)
    isinteger(idx) && return validate_index(Int(idx), loc, topo, N, H)
    return throw(ArgumentError("$idx are not supported window indices for Field!"))
end

validate_index(::Colon, loc, topo, N, H) = Colon()
validate_index(idx::UnitRange, ::Type{Nothing}, topo, N, H) = UnitRange(1, 1)

function validate_index(idx::UnitRange, loc, topo, N, H)
    all_idx = all_indices(loc, topo, N, H)
    (first(idx) ∈ all_idx && last(idx) ∈ all_idx) || throw(ArgumentError("The indices $idx must slice $I"))
    return idx
end

validate_index(idx::Int, args...) = validate_index(UnitRange(idx, idx), args...)

validate_indices(indices, loc, grid::AbstractGrid) =
    validate_index.(indices, loc, topology(grid), size(loc, grid), halo_size(grid))

#####
##### Some basic constructors
#####

# Common outer constructor for all field flavors that performs input validation
function Field(loc::Tuple, grid::AbstractGrid, data, bcs, indices, op=nothing, status=nothing)
    @apply_regionally validate_field_data(loc, data, grid, indices)
    @apply_regionally validate_boundary_conditions(loc, grid, bcs)
    indices = validate_indices(indices, loc, grid)
    buffers = FieldBoundaryBuffers(grid, data, bcs)
    LX, LY, LZ = loc
    return Field{LX, LY, LZ}(grid, data, bcs, indices, op, status, buffers)
end

"""
    Field{LX, LY, LZ}(grid::AbstractGrid,
                      T::DataType=eltype(grid); kw...) where {LX, LY, LZ}

Construct a `Field` on `grid` with data type `T` at the location `(LX, LY, LZ)`.
Each of `(LX, LY, LZ)` is either `Center` or `Face` and determines the field's
location in `(x, y, z)` respectively.

Keyword arguments
=================

- `data :: OffsetArray`: An offset array with the fields data. If nothing is providet the
  field is filled with zeros.
- `boundary_conditions`: If nothing is provided, then field is created using the default
  boundary conditions via [`FieldBoundaryConditions`](@ref).

Example
=======

```jldoctest
julia> using Oceananigans

julia> ω = Field{Face, Face, Center}(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)))
1×1×1 Field{Face, Face, Center} on RectilinearGrid on CPU
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 7×7×7 OffsetArray(::Array{Float64, 3}, -2:4, -2:4, -2:4) with eltype Float64 with indices -2:4×-2:4×-2:4
    └── max=0.0, min=0.0, mean=0.0
```
"""
function Field{LX, LY, LZ}(grid::AbstractGrid,
                           T::DataType=eltype(grid);
                           kw...) where {LX, LY, LZ}

    return Field((LX, LY, LZ), grid, T; kw...)
end

function Field(loc::Tuple,
               grid::AbstractGrid,
               T::DataType = eltype(grid);
               indices = default_indices(3),
               data = new_data(T, grid, loc, indices),
               boundary_conditions = FieldBoundaryConditions(grid, loc, indices))

    return Field(loc, grid, data, boundary_conditions, indices, nothing, nothing)
end

Field(z::ZeroField; kw...) = z
Field(f::Field; indices=f.indices) = view(f, indices...) # hmm...

"""
    CenterField(grid; kw...)

Returns `Field{Center, Center, Center}` on `arch`itecture and `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
CenterField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Center, Center}(grid, T; kw...)

"""
    XFaceField(grid; kw...)

Returns `Field{Face, Center, Center}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
XFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Face, Center, Center}(grid, T; kw...)

"""
    YFaceField(grid; kw...)

Returns `Field{Center, Face, Center}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
YFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Face, Center}(grid, T; kw...)

"""
    ZFaceField(grid; kw...)

Returns `Field{Center, Center, Face}` on `grid`.
Additional keyword arguments are passed to the `Field` constructor.
"""
ZFaceField(grid::AbstractGrid, T::DataType=eltype(grid); kw...) = Field{Center, Center, Face}(grid, T; kw...)

#####
##### Field utils
#####

# Canonical `similar` for Field (doesn't transfer boundary conditions)
function Base.similar(f::Field, grid=f.grid)
    loc = location(f)
    return Field(loc,
                 grid,
                 new_data(eltype(parent(f)), grid, loc, f.indices),
                 FieldBoundaryConditions(grid, loc, f.indices),
                 f.indices,
                 f.operand,
                 deepcopy(f.status))
end

"""
    offset_windowed_data(data, loc, grid, indices)

Return an `OffsetArray` of a `view` of `parent(data)` with `indices`.

If `indices === (:, :, :)`, return an `OffsetArray` of `parent(data)`.
"""
function offset_windowed_data(data, loc, grid, indices)
    halo = halo_size(grid)
    topo = topology(grid)

    if indices isa typeof(default_indices(3))
        windowed_parent = parent(data)
    else
        parent_indices = parent_index_range.(indices, loc, topo, halo)
        windowed_parent = view(parent(data), parent_indices...)
    end

    sz = size(grid)
    return offset_data(windowed_parent, loc, topo, sz, halo, indices)
end

"""
    view(f::Field, indices...)

Returns a `Field` with `indices`, whose `data` is
a view into `f`, offset to preserve index meaning.

Example
=======

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(2, 3, 4), x=(0, 1), y=(0, 1), z=(0, 1));

julia> c = CenterField(grid)
2×3×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×9×10 OffsetArray(::Array{Float64, 3}, -2:5, -2:6, -2:7) with eltype Float64 with indices -2:5×-2:6×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> c .= rand(size(c)...);

julia> v = view(c, :, 2:3, 1:2)
2×2×2 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 2×3×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 8×2×2 OffsetArray(view(::Array{Float64, 3}, :, 5:6, 4:5), -2:5, 2:3, 1:2) with eltype Float64 with indices -2:5×2:3×1:2
    └── max=0.854147, min=0.0109059, mean=0.520099

julia> size(v)
(2, 2, 2)

julia> v[2, 2, 2] == c[2, 2, 2]
true
```
"""
function Base.view(f::Field, i, j, k)
    grid = f.grid
    loc = location(f)

    # Validate indices (convert Int to UnitRange, error for invalid indices)
    window_indices = validate_indices((i, j, k), loc, f.grid)

    # Choice: OffsetArray of view of OffsetArray, or OffsetArray of view?
    #     -> the first retains a reference to the original f.data (an OffsetArray)
    #     -> the second loses it, so we'd have to "re-offset" the underlying data to access.
    #     -> we choose the second here, opting to "reduce indirection" at the cost of "index recomputation".
    #
    # OffsetArray around a view of parent with appropriate indices:
    windowed_data = offset_windowed_data(f.data, loc, grid, window_indices)

    return Field(loc,
                 grid,
                 windowed_data,
                 f.boundary_conditions, # keep original boundary conditions
                 window_indices,
                 f.operand,
                 f.status)
end

const WindowedData = OffsetArray{<:Any, <:Any, <:SubArray}
const WindowedField = Field{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:WindowedData}

# Conveniences
Base.view(f::Field, I::Vararg{Colon}) = f
Base.view(f::Field, i) = view(f, i, :, :)
Base.view(f::Field, i, j) = view(f, i, j, :)

boundary_conditions(not_field) = nothing

function boundary_conditions(f::Field)
    if f.indices === default_indices(3) # default boundary conditions
        return f.boundary_conditions
    else # filter boundary conditions in windowed directions
        return FieldBoundaryConditions(f.indices, f.boundary_conditions)
    end
end

immersed_boundary_condition(f::Field) = f.boundary_conditions.immersed

data(field::Field) = field.data

indices(obj, i=default_indices(3)) = i
indices(f::Field, i=default_indices(3)) = f.indices
indices(a::SubArray, i=default_indices(ndims(a))) = a.indices
indices(a::OffsetArray, i=default_indices(ndims(a))) = indices(parent(a), i)

"""Return indices that create a `view` over the interior of a Field."""
interior_view_indices(field_indices, interior_indices) = Colon()
interior_view_indices(::Colon,       interior_indices) = interior_indices

function interior(a::OffsetArray,
                  loc::Tuple,
                  topo::Tuple,
                  size::NTuple{N, Int},
                  halo_size::NTuple{N, Int},
                  ind::Tuple=default_indices(3)) where N

    i_interior = interior_parent_indices.(loc, topo, size, halo_size)
    i_view = interior_view_indices.(ind, i_interior)
    return view(parent(a), i_view...)
end

"""
    interior(f::Field)

Returns a view of `f` that excludes halo points."
"""
interior(f::Field) = interior(f.data, location(f), f.grid, f.indices)
interior(a::OffsetArray, loc, grid, indices) = interior(a, loc, topology(grid), size(grid), halo_size(grid), indices)
interior(f::Field, I...) = view(interior(f), I...)

# Don't use axes(f) to checkbounds; use axes(f.data)
Base.checkbounds(f::Field, I...) = Base.checkbounds(f.data, I...)

function Base.axes(f::Field)
    if f.indices === (:, : ,:)
        return Base.OneTo.(size(f))
    else
        return Tuple(f.indices[i] isa Colon ? Base.OneTo(size(f, i)) : f.indices[i] for i = 1:3)
    end
end

@propagate_inbounds Base.getindex(f::Field, inds...) = getindex(f.data, inds...)
@propagate_inbounds Base.getindex(f::Field, i::Int)  = parent(f)[i]
@propagate_inbounds Base.setindex!(f::Field, val, i, j, k) = setindex!(f.data, val, i, j, k)
@propagate_inbounds Base.lastindex(f::Field) = lastindex(f.data)
@propagate_inbounds Base.lastindex(f::Field, dim) = lastindex(f.data, dim)

Base.fill!(f::Field, val) = fill!(parent(f), val)
Base.parent(f::Field) = parent(f.data)
Adapt.adapt_structure(to, f::Field) = Adapt.adapt(to, f.data)

length_indices(N, i::Colon) = N
length_indices(N, i::UnitRange) = length(i)

total_size(f::Field) = length_indices.(total_size(location(f), f.grid), f.indices)
Base.size(f::Field)  = length_indices.(      size(location(f), f.grid), f.indices)

#####
##### Interface for field computations
#####

"""
    compute!(field)

Computes `field.data` from `field.operand`.
"""
compute!(field, time=nothing) = field # fallback

"""
    @compute(exprs...)

Call `compute!` on fields after defining them.
"""
macro compute(def)
    expr = Expr(:block)
    field = def.args[1]
    push!(expr.args, :($(esc(def))))
    push!(expr.args, :(compute!($(esc(field)))))
    return expr
end

# Computation "status" for avoiding unnecessary recomputation
mutable struct FieldStatus{T}
    time :: T
end

FieldStatus() = FieldStatus(0.0)
Adapt.adapt_structure(to, status::FieldStatus) = (; time = status.time)

"""
    compute_at!(field, time)

Computes `field.data` at `time`. Falls back to compute!(field).
"""
compute_at!(field, time) = compute!(field)

"""
    compute_at!(field, time)

Computes `field.data` if `time != field.status.time`.
"""
function compute_at!(field::Field, time)
    if isnothing(field.status) # then always compute:
        compute!(field, time)

    # Otherwise, compute only on initialization or if field.status.time is not current,
    elseif time == zero(time) || time != field.status.time
        compute!(field, time)
        field.status.time = time
    end

    return field
end

# This edge case occurs if `fetch_output` is called with `model::Nothing`.
# We do the safe thing here and always compute.
compute_at!(field::Field, ::Nothing) = compute!(field, nothing)

#####
##### Fields that are reduced along one or more dimensions
#####

const XReducedField = Field{Nothing}
const YReducedField = Field{<:Any, Nothing}
const ZReducedField = Field{<:Any, <:Any, Nothing}

const YZReducedField = Field{<:Any, Nothing, Nothing}
const XZReducedField = Field{Nothing, <:Any, Nothing}
const XYReducedField = Field{Nothing, Nothing, <:Any}

const XYZReducedField = Field{Nothing, Nothing, Nothing}

const ReducedField = Union{XReducedField, YReducedField, ZReducedField,
                           YZReducedField, XZReducedField, XYReducedField,
                           XYZReducedField}

reduced_dimensions(field::Field)           = ()
reduced_dimensions(field::XReducedField)   = tuple(1)
reduced_dimensions(field::YReducedField)   = tuple(2)
reduced_dimensions(field::ZReducedField)   = tuple(3)
reduced_dimensions(field::YZReducedField)  = (2, 3)
reduced_dimensions(field::XZReducedField)  = (1, 3)
reduced_dimensions(field::XYReducedField)  = (1, 2)
reduced_dimensions(field::XYZReducedField) = (1, 2, 3)

@propagate_inbounds Base.getindex(r::XReducedField, i, j, k) = getindex(r.data, 1, j, k)
@propagate_inbounds Base.getindex(r::YReducedField, i, j, k) = getindex(r.data, i, 1, k)
@propagate_inbounds Base.getindex(r::ZReducedField, i, j, k) = getindex(r.data, i, j, 1)

@propagate_inbounds Base.setindex!(r::XReducedField, v, i, j, k) = setindex!(r.data, v, 1, j, k)
@propagate_inbounds Base.setindex!(r::YReducedField, v, i, j, k) = setindex!(r.data, v, i, 1, k)
@propagate_inbounds Base.setindex!(r::ZReducedField, v, i, j, k) = setindex!(r.data, v, i, j, 1)

@propagate_inbounds Base.getindex(r::YZReducedField, i, j, k) = getindex(r.data, i, 1, 1)
@propagate_inbounds Base.getindex(r::XZReducedField, i, j, k) = getindex(r.data, 1, j, 1)
@propagate_inbounds Base.getindex(r::XYReducedField, i, j, k) = getindex(r.data, 1, 1, k)

@propagate_inbounds Base.setindex!(r::YZReducedField, v, i, j, k) = setindex!(r.data, v, i, 1, 1)
@propagate_inbounds Base.setindex!(r::XZReducedField, v, i, j, k) = setindex!(r.data, v, 1, j, 1)
@propagate_inbounds Base.setindex!(r::XYReducedField, v, i, j, k) = setindex!(r.data, v, 1, 1, k)

@propagate_inbounds Base.getindex(r::XYZReducedField, i, j, k) = getindex(r.data, 1, 1, 1)
@propagate_inbounds Base.setindex!(r::XYZReducedField, v, i, j, k) = setindex!(r.data, v, 1, 1, 1)

# Preserve location when adapting fields reduced on one or more dimensions
function Adapt.adapt_structure(to, reduced_field::ReducedField)
    LX, LY, LZ = location(reduced_field)
    return Field{LX, LY, LZ}(nothing,
                             adapt(to, reduced_field.data),
                             nothing,
                             nothing,
                             nothing,
                             nothing,
                             nothing)
end

#####
##### Field reductions
#####

# TODO: needs test
Statistics.dot(a::Field, b::Field) = mapreduce((x, y) -> x * y, +, interior(a), interior(b))

# TODO: in-place allocations with function mappings need to be fixed in Julia Base...
const SumReduction     = typeof(Base.sum!)
const MeanReduction    = typeof(Statistics.mean!)
const ProdReduction    = typeof(Base.prod!)
const MaximumReduction = typeof(Base.maximum!)
const MinimumReduction = typeof(Base.minimum!)
const AllReduction     = typeof(Base.all!)
const AnyReduction     = typeof(Base.any!)

initialize_reduced_field!(::SumReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), Base.add_sum, true, interior(c))
initialize_reduced_field!(::ProdReduction, f, r::ReducedField, c) = Base.initarray!(interior(r), Base.mul_prod, true, interior(c))
initialize_reduced_field!(::AllReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), &, true, interior(c))
initialize_reduced_field!(::AnyReduction,  f, r::ReducedField, c) = Base.initarray!(interior(r), |, true, interior(c))

initialize_reduced_field!(::MaximumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), interior(c))
initialize_reduced_field!(::MinimumReduction, f, r::ReducedField, c) = Base.mapfirst!(f, interior(r), interior(c))

filltype(f, c) = eltype(c)
filltype(::Union{AllReduction, AnyReduction}, grid) = Bool

function reduced_location(loc; dims)
    if dims isa Colon
        return (Nothing, Nothing, Nothing)
    else
        return Tuple(i ∈ dims ? Nothing : loc[i] for i in 1:3)
    end
end

reduced_indices(indices; dims) = Tuple(i ∈ dims ? Colon() : indices[i] for i in 1:3)

function reduced_dimension(loc)
    dims = ()
    for i in 1:3
        loc[i] == Nothing ? dims = (dims..., i) : dims
    end
    return dims
end

## Allow support for ConditionalOperation

get_neutral_mask(::Union{AllReduction, AnyReduction})  = true
get_neutral_mask(::Union{SumReduction, MeanReduction}) =   0
get_neutral_mask(::MinimumReduction) =   Inf
get_neutral_mask(::MaximumReduction) = - Inf
get_neutral_mask(::ProdReduction)    =   1

# If func = identity and condition = nothing, nothing happens
"""
    condition_operand(f::Function, op::AbstractField, condition, mask)

Wrap `f(op)` in `ConditionedOperand` with `condition` and `mask`. `f` defaults to `identity`.

If `f isa identity` and `isnothing(condition)` then `op` is returned without wrapping.

Otherwise return `ConditionedOperand`, even when `isnothing(condition)` but `!(f isa identity)`.
"""
@inline condition_operand(op::AbstractField, condition, mask) = condition_operand(identity, op, condition, mask)
@inline condition_operand(::typeof(identity), operand::AbstractField, ::Nothing, mask) = operand

@inline conditional_length(c::AbstractField)        = length(c)
@inline conditional_length(c::AbstractField, dims)  = mapreduce(i -> size(c, i), *, unique(dims); init=1)

# Allocating and in-place reductions
for reduction in (:sum, :maximum, :minimum, :all, :any, :prod)

    reduction! = Symbol(reduction, '!')

    @eval begin

        # In-place
        function Base.$(reduction!)(f::Function,
                                    r::ReducedField,
                                    a::AbstractField;
                                    condition = nothing,
                                    mask = get_neutral_mask(Base.$(reduction!)),
                                    kwargs...)

            return Base.$(reduction!)(identity,
                                      interior(r),
                                      condition_operand(f, a, condition, mask);
                                      kwargs...)
        end

        function Base.$(reduction!)(r::ReducedField,
                                    a::AbstractField;
                                    condition = nothing,
                                    mask = get_neutral_mask(Base.$(reduction!)),
                                    kwargs...)

            return Base.$(reduction!)(identity,
                                      interior(r),
                                      condition_operand(a, condition, mask);
                                      kwargs...)
        end

        # Allocating
        function Base.$(reduction)(f::Function,
                                   c::AbstractField;
                                   condition = nothing,
                                   mask = get_neutral_mask(Base.$(reduction!)),
                                   dims = :)

            T = filltype(Base.$(reduction!), c)
            loc = reduced_location(location(c); dims)
            r = Field(loc, c.grid, T; indices=indices(c))
            conditioned_c = condition_operand(f, c, condition, mask)
            initialize_reduced_field!(Base.$(reduction!), identity, r, conditioned_c)
            Base.$(reduction!)(identity, r, conditioned_c, init=false)

            if dims isa Colon
                return CUDA.@allowscalar first(r)
            else
                return r
            end
        end

        Base.$(reduction)(c::AbstractField; kwargs...) = Base.$(reduction)(identity, c; kwargs...)
    end
end

function Statistics._mean(f, c::AbstractField, ::Colon; condition = nothing, mask = 0)
    operator = condition_operand(f, c, condition, mask)
    return sum(operator) / conditional_length(operator)
end

function Statistics._mean(f, c::AbstractField, dims; condition = nothing, mask = 0)
    operand = condition_operand(f, c, condition, mask)
    r = sum(operand; dims)
    n = conditional_length(operand, dims)
    r ./= n
    return r
end

Statistics.mean(f::Function, c::AbstractField; condition = nothing, dims=:) = Statistics._mean(f, c, dims; condition)
Statistics.mean(c::AbstractField; condition = nothing, dims=:) = Statistics._mean(identity, c, dims; condition)

function Statistics.mean!(f::Function, r::ReducedField, a::AbstractField; condition = nothing, mask = 0)
    sum!(f, r, a; condition, mask, init=true)
    dims = reduced_dimension(location(r))
    n = conditional_length(condition_operand(f, a, condition, mask), dims)
    r ./= n
    return r
end

Statistics.mean!(r::ReducedField, a::AbstractArray; kwargs...) = Statistics.mean!(identity, r, a; kwargs...)

function Statistics.norm(a::AbstractField; condition = nothing)
    r = zeros(a.grid, 1)
    Base.mapreducedim!(x -> x * x, +, r, condition_operand(a, condition, 0))
    return CUDA.@allowscalar sqrt(r[1])
end

function Base.isapprox(a::AbstractField, b::AbstractField; kw...)
    conditioned_a = condition_operand(a, nothing, one(eltype(a)))
    conditioned_b = condition_operand(b, nothing, one(eltype(b)))
    # TODO: Make this non-allocating?
    return all(isapprox.(conditioned_a, conditioned_b; kw...))
end

#####
##### fill_halo_regions!
#####

function fill_halo_regions!(field::Field, args...; kwargs...)
    reduced_dims = reduced_dimensions(field)

    # To correctly fill the halo regions of fields with non-default indices, we'd have to
    # offset indices in the fill halo regions kernels.
    # For now we punt and don't support filling halo regions on windowed fields.
    # Note that `FieldBoundaryConditions` _can_ filter boundary conditions in
    # windowed directions:
    #
    #   filtered_bcs = FieldBoundaryConditions(field.indices, field.boundary_conditions)
    #
    # which will be useful for implementing halo filling for windowed fields in the future.
    if field.indices isa typeof(default_indices(3))
        fill_halo_regions!(field.data,
                           field.boundary_conditions,
                           instantiated_location(field),
                           field.grid,
                           args...;
                           reduced_dimensions = reduced_dims,
                           kwargs...)
    end

    return nothing
end


text:
using Revise
using Test
using LinearAlgebra
using MFGnet
using Printf

d = 2
m = 8
nex = 4
nTh = 3

Rs = [Float64;Float32]

for k=1:2
    @testset "derivative check $(Rs[k])" begin
    R = Rs[k]
    L = SingleLayer()
    Φ = PotentialNN(L)
    s = randn(R,d+1,nex)
    s[d+1,:] .= R(0)
    K = randn(R,m,d+1)
    b = randn(R,m)
    ΘN = (K,b)
    w = randn(R,m)
    A = randn(R,d+1,d+1)
    c = randn(R,d+1)
	b = randn(R,1)
    Θ = (w,ΘN,A,c,b)

    Φ0 = Φ(s,Θ)
    @test eltype(Φ0) == R
    ds = randn(R,size(s))
    dΦ,d2Φ = getGradAndHessian(Φ,s,Θ)
    @test eltype(dΦ) == R
    @test eltype(d2Φ)==R
    dΦds = sum(ds.*dΦ,dims=1)

    curv = zeros(1,nex)
    for j=1:nex
        curv[1,j] = ds[:,j]'*d2Φ[:,:,j]*ds[:,j]
    end
    trH = MFGnet.getTraceHess(Φ,s,Θ)
    for j=1:nex
        @test abs(trH[j] - tr(d2Φ[1:end-1,1:end-1,j]))/abs(tr(d2Φ[1:end-1,1:end-1,j])) < sqrt(eps(R))
    end

    for j=1:15
        h = R(2.0^(-j))
        Φt = Φ(s + h .* ds,Θ)

        E0 = norm(Φ0-Φt)
        E1 = norm(Φ0+h*dΦds - Φt)
        E2 = norm(Φ0+h*dΦds+0.5*h^2*curv - Φt)

        @printf("h=%1.3e\t\tE0=%1.3e\tE1=%1.3e\tE1=%1.3e\n",h,E0,E1,E2)
    end
    @test true
    end
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using ProgressMeter
using MATLAB
using Printf
include("PSS_function.jl")

# H=range(.11,stop=0.145,length=107)
Energy=0.112012
Energy=.13
# @showprogress 1 "Computing..." for Energy in H

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
# mat"""
# set(gca,'XTick',[])
# set(gca,'YTick',[])
# set(gca,'xticklabel',[])
# set(gca,'yticklabel',[])
# set(gca,'XColor','none')
# set(gca,'YColor','none')
# set(gca,'XColor','none')
# set(gca,'YColor','none')
# """
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/"
h=replace(@sprintf("%.15f",Energy),"."=>"_")
file_name=location*h*".fig"
file_name=h*".fig"
COLOR=["#393b79" ,"#5254a3","#6b6ecf","#9c9ede" ,"#637939","#8ca252" ,"#b5cf6b" ,"#cedb9c" ,"#8c6d31","#bd9e39" ,"#e7ba52","#e7cb94","#843c39","#ad494a" ,"#d6616b","#e7969c" ,"#7b4173" ,"#a55194","#ce6dbd" ,"#de9ed6"]
println(file_name)
mat"axis([ -2.5, 2.5,-1.5, 1.5])"
#defining colors for PSS
N_iter_Q=10;#50
Q_start=-.3
Q_end=.3
N_iter_P=10
P_start=-.01
P_end=.01
t_end=1e3
max_hit=1000
mat"axis([ $P_start, $P_end, $Q_start, $Q_end])"
#create an empty list to store dH
H_differences=zeros(0)
#julia's version of linspace
ArrP=range(P_start,stop=P_end,length=N_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=N_iter_Q)

@showprogress 1 "Computing..." for k=1:N_iter_P
    for j=1:N_iter_Q
        Q1,P1,dH=PSS_function(ArrQ[j], ArrP[k], Energy, t_end, max_hit)
            # Q2,P2,dH=PSS_function(-ArrQ[k], ArrP[j], Energy, t_end, max_hit)
            # Q3,P3,dH=PSS_function(ArrQ[k], -ArrP[j], Energy, t_end, max_hit)
            # Q4,P4,dH=PSS_function(-ArrQ[k], -ArrP[j], Energy, t_end, max_hit)
            # Q=vcat(Q1,Q2)
            # P=vcat(P1,P2)
            # Q_2=vcat(Q3,Q4)
            # P_2=vcat(P3,P4)
            if dH<1e-10
                current_color=COLOR[mod(k,length(COLOR))+1]
                mat"plot($P1,$Q1,'.','MarkerSize',.1,'color',$current_color); hold on;"
            end

    end
end

#
mat"savefig($file_name)"
# mat"close"
# end


text:
struct Porkchop
    startorbit::Orbit
    endorbit::Orbit
    starttimes::Vector{Float64}
    flighttimes::Vector{Float64}
    # transfers::Matrix{Transfer}
    departure_Δv::Matrix{Float64}
    arrival_Δv::Matrix{Float64}
    Δv::Matrix{Float64}
end

function Porkchop(startorb::Orbit, endorb::Orbit, starttimes::AbstractVector{<:Real}, flighttimes::AbstractVector{<:Real};
                  tfer_fun = Transfer, match_departure_plane=false, match_arrival_plane=false, kwargs...)
    departure_Δv = fill(NaN, length(starttimes), length(flighttimes))
    arrival_Δv = fill(NaN, length(starttimes), length(flighttimes))
    Δv = fill(NaN, length(starttimes), length(flighttimes))
    tferpath=path_to_body(startorb.primary, endorb.primary)
    cmnparent=closest_common_parent(startorb, endorb)
    for (i,st) in enumerate(starttimes)
        for (j,ft) in enumerate(flighttimes)
            tfer = tfer_fun(startorb, endorb, st, st+ft; transferpath=tferpath, commonparent=cmnparent, kwargs...)
            if match_departure_plane
                tfer = match_departure_inclination(tfer)
            end
            if match_arrival_plane
                tfer = match_arrival_inclination(tfer)
            end
            departure_Δv[i,j] = norm(tfer.burns[1].Δv̄)
            arrival_Δv[i,j] = norm(tfer.burns[end].Δv̄)
            Δv[i,j] = tfer.Δv
        end
    end
    # remove NaN and replace with highes value for Δv
    # this is done to avoid JSON serialization problems, and probably should be moved elsewhere
    departure_Δv[isnan.(departure_Δv)] .= maximum(x->isnan(x) ? -Inf : x, departure_Δv)
    arrival_Δv[isnan.(arrival_Δv)] .= maximum(x->isnan(x) ? -Inf : x, arrival_Δv)
    Δv[isnan.(Δv)] .= maximum(x->isnan(x) ? -Inf : x, Δv)
    return Porkchop(startorb, endorb, starttimes, flighttimes, departure_Δv, arrival_Δv, Δv)
end

Porkchop(startorb::Orbit, endorb::Orbit, stime1::Real, stime2::Real, ftime1::Real, ftime2::Real; npts=100, kwargs...
    ) = Porkchop(startorb, endorb, collect(range(stime1, stop=stime2, length=npts)), collect(range(ftime1, stop=ftime2, length=npts)); tfer_fun = Transfer, kwargs...)

function Porkchop(startorb::Orbit, endorb::Orbit, stime1::Real=0., stime2=nothing, ftime1=nothing, ftime2=nothing; kwargs...)
    transferbody = closest_common_parent(startorb, endorb)
    if isnothing(stime2) || isnothing(ftime1) || isnothing(ftime2)
        sorb = startorb
        while sorb.primary!=transferbody
            sorb = sorb.primary.orbit
        end
        eorb = endorb
        while eorb.primary!=transferbody
            eorb = eorb.primary.orbit
        end
    end
    # complete a full rotation of one orbit relative to the other to get start time range
    if isnothing(stime2)
        if isapprox(1/sorb.period, 1/eorb.period)
            stime2 = stime1 + sorb.period
        else
            stime2 = stime1 + 1.25*abs(1/(1/sorb.period - 1/eorb.period))   # a full relative rotation + 25%
        end
    end
    # choose flight times based on the apses of sorb and eorb
    if isnothing(ftime1) || isnothing(ftime2)
        hohmannaxis = max(abs(sorb.a*(1-sorb.e) + eorb.a*(1+eorb.e)), abs(sorb.a*(1+sorb.e) + eorb.a*(1-eorb.e)))
        ftimemid = period(hohmannaxis/2, transferbody.μ)/2
        if isnothing(ftime1) && isnothing(ftime2)
            ftime1 = ftimemid/2
            ftime2 = ftimemid*2
        elseif isnothing(ftime1)
            ftime1 = min(ftimemid/2, ftime2/2)
        else
            ftime2 = max(ftimemid*2, ftime1*2)
        end
    end
    return Porkchop(startorb, endorb, stime1, stime2, ftime1, ftime2; kwargs...)
end

fastPorkchop(args...; kwargs...) = Porkchop(args...; tfer_fun=fastTransfer, kwargs...)

### other methods ###

# returns best start time and end time
function best_start_end_times(pc::Porkchop; field=:Δv)
    (sidx, fidx) = Tuple(findmin(getfield(pc, field))[2])
    return pc.starttimes[sidx], pc.starttimes[sidx]+pc.flighttimes[fidx]
end

function best_transfer(pc::Porkchop; field=:Δv, tfer_fun=Transfer, kwargs...)
    stime, etime = best_start_end_times(pc; field=field)
    return tfer_fun(pc.startorbit, pc.endorbit, stime, etime; kwargs...)
end

### descriptive display ###

Base.show(io::IO, pc::Porkchop) = println(io,
    "Porkchop for transfer from ",
    "orbit around $(pc.startorbit.primary.name) to ",
    "orbit around $(pc.endorbit.primary.name)\n",
    "  Best Δv:\t$(minimum(pc.Δv)) m/s",
)