text:
using DemeNet: DemeSpec, Deme, Signer, Certificate, Intent, save
using PeaceCypher

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:DemeNet)
save(demespec)
DemeSpec(demespec.uuid)==demespec # All fields are equal. Need to implement equality for this test

# # Now some testing

deme = Deme(demespec)
maintainer = Signer(deme,"maintainer")

notary = deme.notary

signer = notary.Signer()
msg = "Hello World"
signature = notary.Signature(msg,signer)
notary.verify(msg,signature)

### Let's test Certificates and intent
maintainer = Signer(deme,"maintainer")
cert = Certificate("hello world",maintainer)
intent = Intent(cert,notary)


text:
export PosTag

struct PosTag
  name::String
end

ext: jl

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using BlockBandedMatrices, Test, LinearAlgebra

include("test_blockbanded.jl")
include("test_blockskyline.jl")
include("test_bandedblockbanded.jl")
include("test_broadcasting.jl")
include("test_linalg.jl")
include("test_misc.jl")
include("test_triblockbanded.jl")
include("test_adjtransblockbanded.jl")
include("test_blockskylineqr.jl")


text:
# =================
# Dependencies
# =================

using ReachabilityAnalysis, CarlemanLinearization
using Plots, LaTeXStrings, LinearAlgebra, SparseArrays
using Plots.PlotMeasures

using LazySets: center
using CarlemanLinearization: _error_bound_specabs_R

include("../utils.jl")

# =================
# Model definition
# =================

function quadratic(; a, b, N)
    dl = zeros(N-1)
    d = [a*i for i in 1:N]
    du = [b*i for i in 1:N-1]
    return Tridiagonal(dl, d, du)
end

function logistic(; r, K, N)
    return quadratic(a=r, b=-r/K, N=N)
end

# define model
r = -0.5
K = 0.8

# standard form
F1 = hcat(r)
F2 = hcat(-r/K)
a = r
b = -r/K
A = logistic(r=r, K=K, N=4);

x0 = interval(0.5)
R, Re_λ1 = _error_bound_specabs_R(x0, F1, F2; check=true)

using Revise
Tmax = 10.0

x0set = interval(0.47, 0.53)

NSTEPS = 1000

# analytic solution
u0 = mid(x0)
tsp = range(0.0, Tmax, length=NSTEPS+1)
f_analytic(t, u0) = u0 * a * exp(a * t) / (a + b*(1 - exp(a*t))*u0)
sol_a = f_analytic.(tsp, u0)

# solution with Carleman, discrete time
function sol_carlin_discrete_time(; N)
    A = logistic(r=r, K=K, N=N)
    y0 = center(kron_pow_stack(x0, N)) |> Singleton
    prob = @ivp(y' = A * y, y(0) ∈ y0, dim=N)

    return solve(prob, T=Tmax, alg=ORBIT(δ=0.5, approx_model=NoBloating()))
end

# solution with Carleman, continuous time
function sol_carlin_continuous_time(; N)
    F1 = hcat(r); F2 = hcat(-r/K)
    dirs = _template(n=1, N=N)
    alg = LGG09(δ=0.1, template=dirs, approx_model=Forward())
    return _solve_CARLIN(x0set, F1, F2; alg=alg, N=N, T=Tmax, bloat=false)
end

# ============
# Figure 1a
# ============

function figure_1a()
    sol_discr_N1 = sol_carlin_discrete_time(N=1)
    sol_discr_N2 = sol_carlin_discrete_time(N=2)
    sol_discr_N3 = sol_carlin_discrete_time(N=3)
    sol_discr_N4 = sol_carlin_discrete_time(N=4)
    sol_discr_N5 = sol_carlin_discrete_time(N=5)
    sol_discr_N6 = sol_carlin_discrete_time(N=6);

    fig = plot(legend=:topright, xlab=L"t", ylab=L"x(t)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick=([0.0, 2.5, 5.0, 7.5, 10.0], [L"0", L"2.5", L"5.0", L"7.5", L"10.0"]),
               ytick=([0.0, 0.1, 0.2, 0.3, 0.4, 0.5], [L"0", L"0.1", L"0.2", L"0.3", L"0.4", L"0.5"]),
               xlims=(0.0, 10.1),
               ylims=(-0.01, 0.5),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    plot!(fig, tsp, sol_a, c=:magenta, lab=L"\textrm{Analytic}", lw=3.0, alpha=1.)

    plot!(fig, sol_discr_N1, vars=(0, 1), lab=L"N=1", seriestype=:path, linestyle=:dash, c=:green, markersize=8, markershape=:circle, lw=2.0, alpha=1.)
    plot!(fig, sol_discr_N2, vars=(0, 1), lab=L"N=2", seriestype=:path, linestyle=:dash, c=:red, markersize=8, markershape=:utriangle, lw=2.0, alpha=1.)
    plot!(fig, sol_discr_N4, vars=(0, 1), lab=L"N=4", seriestype=:path, linestyle=:dash, c=:blue, markersize=8, markershape=:rect, lw=2.0, alpha=1.)

    plot!(fig, sol_discr_N6, vars=(0, 1), lab=L"N=6", seriestype=:path, linestyle=:dash, c=:orange, markersize=8, markershape=:star, lw=2.0, alpha=1.)

end

fig = figure_1a()
savefig(fig, joinpath(TARGET_FOLDER, "figure_1a.pdf"))

# ============
# Figure 1b
# ============

function figure_1b()
    sol_discr_N1 = sol_carlin_discrete_time(N=1)
    sol_discr_N2 = sol_carlin_discrete_time(N=2)
    sol_discr_N3 = sol_carlin_discrete_time(N=3)
    sol_discr_N4 = sol_carlin_discrete_time(N=4)
    sol_discr_N5 = sol_carlin_discrete_time(N=5)
    sol_discr_N6 = sol_carlin_discrete_time(N=6);

    fig = plot(legend=:topright, xlab=L"t", ylab=L"x(t)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick=([0.0, 2.5, 5.0, 7.5, 10.0], [L"0", L"2.5", L"5.0", L"7.5", L"10.0"]),
               ytick=([0.0, 0.1, 0.2, 0.3, 0.4, 0.5], [L"0", L"0.1", L"0.2", L"0.3", L"0.4", L"0.5"]),
               xlims=(0.0, 10.1),
               ylims=(-0.01, 0.6),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    tsp = range(0.0, Tmax, length=NSTEPS+1)

    sol_a_high = f_analytic.(tsp, x0set.hi)
    sol_a_low = f_analytic.(tsp, x0set.lo)

    # initial set
    plot!(Singleton([0.0]) × Interval(x0set), lw=8.0, alpha=1., c=:black, lab="", markershape=:none, seriestype=:shape, xlims=(-0.05, 10.1))

    plot!(fig, sol_carlin_continuous_time(N=1), vars=(0, 1), lab="", c=:green, lc=:green)
    plot!(fig, sol_carlin_continuous_time(N=2), vars=(0, 1), lab="", c=:red, lc=:red)
    plot!(fig, sol_carlin_continuous_time(N=6), vars=(0, 1), lab="", c=:orange, lc=:orange)

    plot!(fig, sol_discr_N1, vars=(0, 1), lab=L"N=1", seriestype=:path, linestyle=:dash, c=:green, markersize=8, markershape=:circle, lw=2.0, alpha=1.)
    plot!(fig, sol_discr_N2, vars=(0, 1), lab=L"N=2", seriestype=:path, linestyle=:dash, c=:red, markersize=8, markershape=:utriangle, lw=2.0, alpha=1.)
    plot!(fig, sol_discr_N6, vars=(0, 1), lab=L"N=6", seriestype=:path, linestyle=:dash, c=:orange, markersize=8, markershape=:star, lw=2.0, alpha=1.)

    plot!(fig, tsp, sol_a_high, c=:magenta, lab=L"\textrm{Analytic}", lw=3.0, alpha=1.)
    plot!(fig, tsp, sol_a_low, c=:magenta, lab="", lw=3.0, alpha=1.)

end

fig = figure_1b()
savefig(fig, joinpath(TARGET_FOLDER, "figure_1b.pdf"))


text:
export slide, ideal_start_end

using Dates

# From SeisNoise
# https://github.com/tclements/SeisNoise.jl/blob/2a816653f5119c3276421938b136c141b8fa51da/src/slicing.jl#L55
"""
    slide(A, cc_len, cc_step, fs, starttime, endtime)
Cut `A` into sliding windows of length `cc_len` points and offset `cc_step` points.
# Arguments
- `A::AbstractArray`: 1D time series.
- `cc_len::Real`: Cross-correlation window length in seconds.
- `cc_step::Real`: Step between cross-correlation windows in seconds.
- `starttime::Float64`: Time of first sample in `A` in unix time.
- `endtime::Float64`: Time of last sample in `A` in unix time.
# Returns
- `out::Array`: Array of sliding windows
- `starts::Array`: Array of start times of each window, in Unix time. E.g to convert
        Unix time to date time, use u2d(starts[1]) = 2018-08-12T00:00:00
"""
function slide(A::AbstractArray, cc_len::Real, cc_step::Real, fs::AbstractFloat,
               starttime::Float64,endtime::Float64)
    N = size(A, 1)
    window_samples = Int(cc_len * fs)
    starts = Array(range(starttime, stop=endtime, step=cc_step))
    ends = starts .+ cc_len .- 1. / fs
    ind = findlast(x -> x <= endtime, ends)
    starts = starts[1:ind]
    ends = ends[1:ind]

    # fill array with overlapping windows
    if cc_step == cc_len && N % cc_len == 0
        return Array(reshape(A,window_samples, N ÷ window_samples)), starts
    elseif cc_step == cc_len # disregard data from edge
        return Array(reshape(A[1 : N - N % window_samples], window_samples, N ÷ window_samples)), starts
    else # need overlap between windows
        out = Array{eltype(A),2}(undef, window_samples,length(starts))
        s = convert.(Int, round.((hcat(starts,ends) .- starttime) .* fs .+ 1.))
        @inbounds for ii in eachindex(starts)
            #out[:,ii] .= @view(A[s[ii,1]:s[ii,1] + window_samples - 1]) # To make sure dimension is consistent
            out[:,ii] .= @view(A[s[ii,1]:s[ii,2]])
        end

    end
  return out, starts
end

"""
    ideal_start_end(S::DateTime, E::DateTime, cc_len::Int, cc_step::Int)
Return ideal start and end times.
"""
function ideal_start_end(S::DateTime, E::DateTime, fs::Float64, cc_len::Real, cc_step::Real)
    starts = Array(S:Second(cc_step):E)
    ends = starts .+ Second(cc_len)
    ind = findlast(x -> x <= E, ends)
    starts = starts[1:ind]
    ends = ends[1:ind]
    return starts, ends
end


text:
using Bridge, StaticArrays, Bridge.Models
const R = ℝ
Random.seed!(2)
Random.seed!(3)
πH = 2000. # prior

t = 1.0
T = 5.00
n = 50001 # total imputed length
m = 100 # number of segments
M = div(n-1,m)
skippoints = 2
dt = (T-t)/(n-1)
tt = t:dt:T
si = 3.
# 10, 20, 8/3 srand(2)
θ² = 1.
γ = 1.
P = Bridge.Models.Pendulum(θ², γ)
Psm = Bridge.Models.Pendulum(θ², 0.0)
x0 = Models.x0(P)

W = sample(tt, Wiener())
X = SamplePath(tt, zeros(ℝ{2}, length(tt)))
Bridge.solve!(Euler(), X, x0, W, P)
W = sample(tt, Wiener())
Xsm = SamplePath(tt, zeros(ℝ{2}, length(tt)))
Bridge.solve!(Euler(), Xsm, x0, W, Psm)

Xtrue = copy(X)

# Observation scheme and subsample
_pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
SV = ℝ{2}
SM = typeof(one(Bridge.outer(zero(SV))))
if isdefined(:partial)
    if !partial
        L = I
        Σ = SDiagonal(1., 1.)
        lΣ = chol(Σ)'
        RV = SV
        RM = SM

        V = SamplePath(collect(_pairs(Xtrue))[1:M:end])
        map!(y -> L*y + lΣ*randn(RV), V.yy, V.yy)
    else
        L = @SMatrix [1.0 0.0]
        Σ = isdefined(:Σ_) ? Σ_ : 1.00
        @show Σ
        lΣ = chol(Σ)'
        RV = Float64
        RM = Float64 #typeof(one(Bridge.outer(zero(RV))))

        V_ = SamplePath(collect(_pairs(Xtrue))[1:M:end])
        V = SamplePath(V_.tt, map(y -> (L*y)[1] + lΣ*randn(RV), V_.yy))
    end
else
    println("Info:  No observation scheme defined. Set `partial in [true, false]`")
end


text:
#=
# Builds a tiled ordering of a 2D or 3D grid
# Tiles are tiledim in size (except edge tiles which could be smaller)
# Grid has griddim size
# Returns the global order in which objects shall be indexed.
# example: 2D, griddim=(4,4), tiledim=(3,3)
#     13  14  15  16
#      7   8   9  12
#      4   5   6  11
#      1   2   3  10
# Invert=true will give  [1, 2, 3, 5, 6, 7 ...]   (encode?)
# Invert=false will give [1, 2, 3, 10, 4, 5, ...] (decode?)
=#
export reorder_grid_tiled, get_tiled_order

function tiled_order_2d(griddim, tiledim, invert = true)
    # griddim and tiledim are tuples of grid dimensions
    gridx = griddim[1];
    gridy = griddim[2];
    tilex = tiledim[1];
    tiley = tiledim[2];

    nnodes = gridx*gridy;
    tiled = zeros(Int64, nnodes); # The ordering

    # Count full tiles and renaining edge nodes.
    fullx = Int(floor(gridx / tilex));
    partialx = gridx - fullx*tilex;
    fully = Int(floor(gridy / tiley));
    partialy = gridy - fully*tiley;

    # Fill in ordering
    for ni=1:nnodes
        # Find destination tile index
        tx = Int(floor(mod(ni-1, gridx)/tilex)) + 1;
        ty = Int(floor((ni-1)/(gridx*tiley))) + 1;
    end

    tind = 0;
    for j = 1:(fully+1)
        ytill = j<=fully ? tiley : partialy;
        for i = 1:(fullx+1)
            xtill = i<=fullx ? tilex : partialx;

            # Add that tile's nodes one at a time
            for tj=1:ytill
                for ti=1:xtill
                    gind = ti + (i-1)*tilex + gridx*(tj + (j-1)*tiley - 1);
                    tind = tind + 1;
                    if invert
                        tiled[gind] = tind;
                    else
                        tiled[tind] = gind;
                    end
                end
            end

        end
    end

    return tiled;
end

function tiled_order_3d(griddim, tiledim, invert = true, tileorder=[])
    # griddim and tiledim are tuples of grid dimensions
    gridx = griddim[1];
    gridy = griddim[2];
    gridz = griddim[3];
    tilex = tiledim[1];
    tiley = tiledim[2];
    tilez = tiledim[3];

    nnodes = gridx*gridy*gridz;
    tiled = zeros(Int64, nnodes); # The ordering

    # Count full tiles and remaining edge nodes.
    fullx = Int(floor(gridx / tilex));
    partialx = gridx - fullx*tilex;
    fully = Int(floor(gridy / tiley));
    partialy = gridy - fully*tiley;
    fullz = Int(floor(gridz / tilez));
    partialz = gridz - fullz*tilez;

    # Fill in ordering
    if length(tileorder) == 0
        Ntiles = (fullx+1)*(fully+1)*(fullz+1);
        tileorder = 1:Ntiles;
    end

    tind = 0;
    for i=1:length(tileorder)
        # Get the tile coordinates in (x,y,z)
        tileind = tileorder[i];
        tx = Int(mod(tileind-1, fullx+1)) + 1;
        ty = Int(floor(mod(tileind-1, (fullx+1)*(fully+1))/(fullx+1))) + 1;
        tz = Int(floor((tileind-1)/((fullx+1)*(fully+1)))) + 1;
        xtill = tx<=fullx ? tilex : partialx;
        ytill = ty<=fully ? tiley : partialy;
        ztill = tz<=fullz ? tilez : partialz;

        # Add that tile's nodes one at a time
        for tk=1:ztill
            for tj=1:ytill
                for ti=1:xtill
                    gind = ti + (tx-1)*tilex + gridx*((tj + (ty-1)*tiley - 1) + gridy*(tk + (tz-1)*tilez - 1));
                    tind = tind + 1;
                    if invert
                        tiled[gind] = tind;
                    else
                        tiled[tind] = gind;
                    end
                end
            end
        end
    end

    return tiled;
end

function reorder_grid_tiled(grid, griddim, tiledim)
    dim = size(grid.allnodes,1);
    nel = size(grid.loc2glb,2);
    nfc = size(grid.face2glb,2);

    tiled = get_tiled_order(dim, griddim, tiledim, true);

    newnodes = zeros(size(grid.allnodes));
    newbdry = copy(grid.bdry);
    newloc2glb = copy(grid.loc2glb);
    newglbvertex = copy(grid.glbvertex);
    newface2glb = copy(grid.face2glb);
    newfaceVertex2glb = copy(grid.faceVertex2glb);

    for mi=1:length(tiled)
        for d=1:dim
            newnodes[d, tiled[mi]] = grid.allnodes[d,mi];
        end
    end

    for bid=1:length(newbdry)
        for i=1:length(newbdry[bid])
            newbdry[bid][i] = tiled[grid.bdry[bid][i]];
        end
    end

    for ei=1:nel
        for ni=1:size(newloc2glb,1)
            newloc2glb[ni,ei] = tiled[grid.loc2glb[ni,ei]];
        end
        for ni=1:size(newglbvertex,1)
            newglbvertex[ni,ei] = tiled[grid.glbvertex[ni,ei]];
        end
    end

    for fi=1:nfc
        for ni=1:size(newface2glb,1)
            newface2glb[ni,fi] = tiled[grid.face2glb[ni,fi]];
        end
        for ni=1:size(newfaceVertex2glb,1)
            newfaceVertex2glb[ni,fi] = tiled[grid.faceVertex2glb[ni,fi]];
        end
    end

    return Femshop.Grid(newnodes, newbdry, grid.bdryface, grid.bdrynorm, grid.bids, newloc2glb, newglbvertex, newface2glb, newfaceVertex2glb);
end

function get_tiled_order(dim, griddim, tiledim, invert = true)
    if dim == 3
        return tiled_order_3d(griddim, tiledim, invert);
    elseif dim == 2
        return tiled_order_2d(griddim, tiledim, invert);
    else
        return 1:griddim[1]; # 1D doesn't make sense
    end
end


text:
type Softmax <: BanditAlgorithm
  temperature::Float64
  counts::Vector{Int64}
  values::Vector{Float64}
end

function Softmax(n_arms::Int64, temperature::Float64)
  Softmax(temperature, zeros(Int64, n_arms), zeros(n_arms))
end

function initialize(algo::Softmax, n_arms::Int64)
  algo.counts = zeros(Int64, n_arms)
  algo.values = zeros(n_arms)
end

function initialize(algo::Softmax, arms::Vector)
  initialize(algo, length(arms))
end

function select_arm(algo::Softmax)
  z = sum(exp(algo.values ./ algo.temperature))
  probs = exp(algo.values ./ algo.temperature) ./ z
  return rand(Categorical(probs))
end

function update(algo::Softmax, chosen_arm::Int64, reward::Real)
  algo.counts[chosen_arm] += 1
  n = algo.counts[chosen_arm]

  value = algo.values[chosen_arm]
  algo.values[chosen_arm] = ((n - 1) / n) * value + (1 / n) * reward
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted, broadcastable, materialize!

struct TensorStyle <: BroadcastStyle end
struct TensorAsScalarStyle <: BroadcastStyle end

BroadcastStyle(::Type{<: AbstractTensor}) = TensorStyle()
BroadcastStyle(x::TensorStyle, ::Broadcast.DefaultArrayStyle{0}) = x
BroadcastStyle(x::TensorStyle, ::Broadcast.Style{Tuple}) = x
BroadcastStyle(x::TensorStyle, ::BroadcastStyle) = TensorAsScalarStyle()
BroadcastStyle(x::TensorAsScalarStyle, ::BroadcastStyle) = TensorAsScalarStyle()

broadcastable(bc::Broadcasted{<: TensorStyle}) = copy(bc)

@generated function _promote_space_for_broadcast(x::Tuple)
    spaces = [Space(t) for t in x.parameters if t <: AbstractTensor]
    quote
        @_inline_meta
        promote_space($(spaces...))
    end
end

_broadcastable_for_tensorstyle(x::Any, ::Type{TT}) where {TT} = x
_broadcastable_for_tensorstyle(x::AbstractTensor, ::Type{TT}) where {TT} = Tuple(convert(TT, x))
broadcastable_for_tensorstyle(x::Tuple{Any}, ::Type{TT}) where {TT} = (_broadcastable_for_tensorstyle(x[1], TT),)
broadcastable_for_tensorstyle(x::Tuple{Any, Any}, ::Type{TT}) where {TT} = (_broadcastable_for_tensorstyle(x[1], TT), _broadcastable_for_tensorstyle(x[2], TT))
broadcastable_for_tensorstyle(x::Tuple, ::Type{TT}) where {TT} = (_broadcastable_for_tensorstyle(x[1], TT), broadcastable_for_tensorstyle(Base.tail(x), TT)...)
@inline function Base.copy(bc::Broadcasted{TensorStyle})
    S = _promote_space_for_broadcast(bc.args)
    TT = tensortype(S)
    TT(broadcast(bc.f, broadcastable_for_tensorstyle(bc.args, TT)...))
end

@inline _ref(x::AbstractTensor) = Ref(x)
@inline _ref(x::Any) = x
@inline function broadcasted(::TensorAsScalarStyle, f, args...)
    broadcasted(f, map(_ref, args)...)
end

# for broadcast!(op, ::Array, ::AbstractTensor)
function materialize!(dest::AbstractArray{<: AbstractTensor}, bc::Broadcasted{TensorStyle})
    materialize!(dest, (copy(bc),))
end
function materialize!(dest::AbstractArray{<: Real}, bc::Broadcasted{TensorStyle})
    copy!(dest, copy(bc))
end


text:
using PowerSystems
using Pkg

const DATA_FOLDER = joinpath(dirname(dirname(Base.find_package("PowerSystems"))),"data")

function build_powersystems()

    if !isdir(DATA_FOLDER)
       Pkg.build("PowerSystems")
    end

end

build_powersystems()


text:
module TurbulenceClosures

export
    AbstractEddyViscosityClosure,
    IsotropicDiffusivity,
    AnisotropicDiffusivity,
    AnisotropicBiharmonicDiffusivity,
    TwoDimensionalLeith,
    ConstantSmagorinsky,
    SmagorinskyLilly,
    AnisotropicMinimumDissipation,
    HorizontallyCurvilinearAnisotropicDiffusivity,
    ConvectiveAdjustmentVerticalDiffusivity,

    DiffusivityFields,
    calculate_diffusivities!,

    ∇_dot_qᶜ,
    ∂ⱼ_τ₁ⱼ,
    ∂ⱼ_τ₂ⱼ,
    ∂ⱼ_τ₃ⱼ,

    cell_diffusion_timescale

using CUDA
using KernelAbstractions

import Oceananigans.Utils: with_tracers

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.BuoyancyModels
using Oceananigans.Utils

using Oceananigans.Architectures: AbstractArchitecture, device

#####
##### Abstract types
#####

"""
    AbstractTurbulenceClosure

Abstract supertype for turbulence closures.
"""
abstract type AbstractTurbulenceClosure{TimeDiscretization} end

# Fallbacks
add_closure_specific_boundary_conditions(closure, boundary_conditions, args...) = boundary_conditions

viscous_flux_ux(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_uy(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_uz(i, j, k, grid, args...) = zero(eltype(grid))

viscous_flux_vx(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_vy(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_vz(i, j, k, grid, args...) = zero(eltype(grid))

viscous_flux_wx(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_wy(i, j, k, grid, args...) = zero(eltype(grid))
viscous_flux_wz(i, j, k, grid, args...) = zero(eltype(grid))

diffusive_flux_x(i, j, k, grid, args...) = zero(eltype(grid))
diffusive_flux_y(i, j, k, grid, args...) = zero(eltype(grid))
diffusive_flux_z(i, j, k, grid, args...) = zero(eltype(grid))

#####
##### Include module code
#####

include("implicit_explicit_time_discretization.jl")
include("turbulence_closure_utils.jl")
include("diffusion_operators.jl")
include("viscous_dissipation_operators.jl")
include("velocity_tracer_gradients.jl")
include("abstract_isotropic_diffusivity_closure.jl")
include("abstract_eddy_viscosity_closure.jl")
include("closure_tuples.jl")

# Implementations:
include("turbulence_closure_implementations/nothing_closure.jl")
include("turbulence_closure_implementations/isotropic_diffusivity.jl")
include("turbulence_closure_implementations/anisotropic_diffusivity.jl")
include("turbulence_closure_implementations/horizontally_curvilinear_anisotropic_diffusivity.jl")
include("turbulence_closure_implementations/horizontally_curvilinear_anisotropic_biharmonic_diffusivity.jl")
include("turbulence_closure_implementations/anisotropic_biharmonic_diffusivity.jl")
include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl")
include("turbulence_closure_implementations/smagorinsky_lilly.jl")
include("turbulence_closure_implementations/anisotropic_minimum_dissipation.jl")
include("turbulence_closure_implementations/tke_based_vertical_diffusivity.jl")
include("turbulence_closure_implementations/convective_adjustment_vertical_diffusivity.jl")

# Miscellaneous utilities
include("diffusivity_fields.jl")
include("turbulence_closure_diagnostics.jl")
include("vertically_implicit_diffusion_solver.jl")

#####
##### Some value judgements here
#####

"""
    ConstantSmagorinsky

An alias for `SmagorinskyLilly`.
"""
const ConstantSmagorinsky = SmagorinskyLilly

end # module


text:
# BioGenerics.IO
# ==============
#
# I/O interfaces for BioJulia packages.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

module IO

# IO Types
# ---------

"Abstract formatted input/output type."
abstract type AbstractFormattedIO end

function (::Type{T})(f, io::Base.IO) where {T <: AbstractFormattedIO}
    fmt = T(io)
    try
        f(fmt)
    finally
        close(fmt)
    end
end

"""
    stream(io::AbstractFormattedIO)

Return the underlying `IO` object; subtypes of `AbstractFormattedIO` must
implement this method.
"""
function stream end

# delegate method call
for f in (:eof, :flush, :close)
    @eval function Base.$(f)(io::AbstractFormattedIO)
        return $(f)(stream(io))
    end
end

"""
Abstract data reader type.

See `subtypes(AbstractReader)` for all available data readers.
"""
abstract type AbstractReader <: AbstractFormattedIO end

Base.IteratorSize(::Type{T}) where T <: AbstractReader = Base.SizeUnknown()

function Base.open(::Type{T}, filepath::AbstractString, args...; kwargs...) where T <: AbstractReader
    return T(open(filepath), args...; kwargs...)
end

function Base.open(f::Function, ::Type{T}, args...; kwargs...) where T <: AbstractFormattedIO
    io = open(T, args...; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

function Base.read(input::AbstractReader)
    return read!(input, eltype(input)())
end

"""
    tryread!(reader::AbstractReader, output)

Try to read the next element into `output` from `reader`.

If the result could not be read, then `nothing` will be returned instead.
"""
function tryread!(reader::AbstractReader, output)
    try
        read!(reader, output)
        return output
    catch ex
        if isa(ex, EOFError)
            return nothing
        end
        rethrow()
    end
end

function Base.iterate(reader::AbstractReader, nextone = eltype(reader)())
    if tryread!(reader, nextone) === nothing
        return nothing
    else
        return copy(nextone), nextone
    end
end

"""
Abstract data writer type.

See `subtypes(AbstractWriter)` for all available data writers.
"""
abstract type AbstractWriter <: AbstractFormattedIO end

function Base.open(::Type{T}, filepath::AbstractString, args...; kwargs_...) where T <: AbstractWriter
    # Special case to improve inference by avoiding the kwarg stuff
    if isempty(kwargs_)
        return T(open(filepath, "w"), args...)
    end
    kwargs::Vector{<:Pair{Symbol}} = collect(kwargs_)
    i = findfirst(kwarg -> kwarg[1] == :append, kwargs)
    if i !== nothing
        append = kwargs[i][2]
        if !isa(append, Bool)
            throw(ArgumentError("append must be boolean"))
        end
        deleteat!(kwargs, i)
    else
        append = false
    end
    return T(open(filepath, append ? "a" : "w"), args...; kwargs...)
end

#

end  # module BioGenerics.IO

max_forks_repo_name: BioJulia/BioGenerics


text:
@inline function prop_decomp(A::SymTen,B::SymTen) #Decompose A regarding B
    α = (A:B)/(B:B)
    ApB = α*B
    return α, ApB
end

@inline function inph_decomp(A::SymTen,B::SymTen)
    λ,e = eigvec(B)

    m1 = e[1] ⋅ A ⋅ e[1]
    m2 = e[2] ⋅ A ⋅ e[2]
    m3 = e[3] ⋅ A ⋅ e[3]

    l1 = λ[1]
    l2 = λ[2]
    l3 = λ[3]
    l12 = l1 - l2
    l13 = l1 - l3
    l23 = l2 - l3
    aux = inv(l12*l13*l23)

    α0 = aux*(-l13*l1*l3*m2 + l2*l2*(l3*m1 - l1*m3) + l2*(l1*l1*m3 - l3*l3*m1))
    α1 = aux*(l3*l3*(m1-m2) + l1*l1*(m2-m3) + l2*l2*(m3-m1))
    α2 = aux*(l3*(m2-m1) + l2*(m1-m3) + l1*(m3-m2))

    AinB = m1*symouter(e[1],e[1]) + m2*symouter(e[2],e[2]) + m3*symouter(e[3],e[3])

    return α0, α1, α2, AinB

end


text:
using Mocking

describe_services_1_patch = @patch function describe_services(a...; b...)
    return Dict{String, Any}("NextToken" => "abcdefg123", "Services" => Any[Dict{String, Any}("AttributeNames" => Any["volumeType", "maxIopsvolume", "instanceCapacity10xlarge", "locationType", "operation"], "ServiceCode" => "AmazonEC2")], "FormatVersion" => "aws_v1")
end

get_attribute_values_1_patch = @patch function get_attribute_values(a...; b...)
    return Dict{String, Any}("AttributeValues" => Any[Dict{String, Any}("Value" => "Throughput Optimized HDD"), Dict{String, Any}("Value" => "Provisioned IOPS")], "NextToken" => "GpgauEXAMPLEezucl5LV0w==:7GzYJ0nw0DBTJ2J66EoTIIynE6O1uXwQtTRqioJzQadBnDVgHPzI1en4BUQnPCLpzeBk9RQQAWaFieA4+DapFAGLgk+Z/9/cTw9GldnPOHN98+FdmJP7wKU3QQpQ8MQr5KOeBkIsAqvAQYdL0DkL7tHwPtE5iCEByAmg9gcC/yBU1vAOsf7R3VaNN4M5jMDv3woSWqASSIlBVB6tgW78YL22KhssoItM/jWW+aP6Jqtq4mldxp/ct6DWAl+xLFwHU/CbketimPPXyqHF3/UXDw==")
end

get_products_1_patch = @patch function get_products(a...; b...)
    return Dict{String, Any}("NextToken" => "57r3EXAMPLEujbzWfHF7Ciw==:ywSmZsD3mtpQmQLQ5XfOsIMkYybSj+vAT+kGmwMFq+K9DGmIoJkz7lunVeamiOPgthdWSO2a7YKojCO+zY4dJmuNl2QvbNhXs+AJ2Ufn7xGmJncNI2TsEuAsVCUfTAvAQNcwwamtk6XuZ4YdNnooV62FjkV3ZAn40d9+wAxV7+FImvhUHi/+f8afgZdGh2zPUlH8jlV9uUtj0oHp8+DhPUuHXh+WBII1E/aoKpPSm3c=", "PriceList" => Any["{\"product\":{\"productFamily\":\"Storage\",\"attributes\":{\"storageMedia\":\"SSD-backed\",\"maxThroughputvolume\":\"320 MB/sec\",\"volumeType\":\"Provisioned IOPS\",\"maxIopsvolume\":\"20000\",\"servicecode\":\"AmazonEC2\",\"usagetype\":\"CAN1-EBS:VolumeUsage.piops\",\"locationType\":\"AWS Region\",\"location\":\"Canada (Central)\",\"servicename\":\"Amazon Elastic Compute Cloud\",\"maxVolumeSize\":\"16 TiB\",\"operation\":\"\"},\"sku\":\"WQGC34PB2AWS8R4U\"},\"serviceCode\":\"AmazonEC2\",\"terms\":{\"OnDemand\":{\"WQGC34PB2AWS8R4U.JRTCKXETXF\":{\"priceDimensions\":{\"WQGC34PB2AWS8R4U.JRTCKXETXF.6YS6EN2CT7\":{\"unit\":\"GB-Mo\",\"endRange\":\"Inf\",\"description\":\"\$0.138 per GB-month of Provisioned IOPS SSD (io1)  provisioned storage - Canada (Central)\",\"appliesTo\":[],\"rateCode\":\"WQGC34PB2AWS8R4U.JRTCKXETXF.6YS6EN2CT7\",\"beginRange\":\"0\",\"pricePerUnit\":{\"USD\":\"0.1380000000\"}}},\"sku\":\"WQGC34PB2AWS8R4U\",\"effectiveDate\":\"2017-08-01T00:00:00Z\",\"offerTermCode\":\"JRTCKXETXF\",\"termAttributes\":{}}}},\"version\":\"20170901182201\",\"publicationDate\":\"2017-09-01T18:22:01Z\"}"], "FormatVersion" => "aws_v1")
end


text:
const nums = parse.(Int, readlines("day01.txt"))

# Find two numbers that sum to 2020 and multiply them.
function part1(nums, target)
    set = Set{Int}()
    for v in nums
        complement = target - v
        complement in set && return complement * v
        push!(set, v)
    end
    return -1
end

part1(nums, 2020)

# Find three numbers that sum to 2020 and multiply them.
function part2(nums)
    for i in 1:length(nums)
        target = 2020 - nums[i]
        prod = part1(vcat(nums[1:i-1], nums[i+1:end]), target)
        prod > 0 && return nums[i] * prod
    end
    return -1
end

part2(nums)

# P.S. These are pretty much TwoSum and ThreeSum problems on LeetCode.


text:
module TestIntegrate
    using Calculus2
    using Base.Test

    @test isapprox(
        integrate(x -> 1 / x, 1.0, 2.0),
        log(2) - log(1)
    )

    @test isapprox(
        integrate(x -> -sin(x), 0.0, pi),
        (cos(pi) - cos(0.0))
    )

    @test isapprox(
        integrate(x -> 1, 0, 1),
        1.0 - 0.0,
    )

    @test isapprox(
        integrate(x -> x, 0, 1),
        1//2 * 1.0^2 - 1//2 * 0.0^2
    )

    @test isapprox(
        integrate(x -> x * x, 0, 1),
        1//3 * 1.0^3 - 1//3 * 0.0^3
    )

    @test isapprox(
        integrate(sin, 0, pi),
        -cos(pi) - (-cos(0))
    )

    @test isapprox(
        integrate(cos, 0, pi),
        sin(pi) - sin(0)
    )

    @test isapprox(
        integrate(x -> sin(x)^2 + sin(x)^2, 0, pi),
        pi
    )

    @test isapprox(
        integrate(x -> 1 / x, 0.01, 1),
        4.60517
    )

    @test isapprox(
        integrate(x -> cos(x)^8, 0, 2 * pi),
        35 * pi / 64
    )

    @test isapprox(
        integrate(x -> sin(x) - sin(x^2) + sin(x^3), pi, 2 * pi),
        -1.830467
    )

    @test isapprox(
        integrate(sin, 0.0, pi, method = :quadrature),
        sqrt(2),
        atol = 10e-1
    )

    @test isapprox(
        integrate(sin, 0.0, pi, method = :simpsons),
        sqrt(2),
        atol = 10e-1
    )

    @test isapprox(
        integrate(sin, 0.0, pi, method = :monte_carlo),
        sqrt(2),
        atol = 10e-1
    )

    @test isapprox(
        ∫(cos, 0.0, 1.0),
        sin(1.0) - sin(0.0)
    )
end


text:
using IRDumps
using Test

@testset "IRDumps.jl" begin
    # Write your tests here.
end


text:
module TestThreading

include("preamble.jl")

@testset "splitrange" begin
    @test splitrange(16, 4) == [1:4, 5:8, 9:12, 13:16]
    @test splitrange(4, 16) == [1:1, 2:2, 3:3, 4:4]
    @test splitrange(16, 5) == [1:4, 5:7, 8:10, 11:13, 14:16]
end

@testset "tseed!" begin
    rngs = [MersenneTwister() for _ = 1:Threads.nthreads()]
    x = zeros(1000)
    y = zeros(1000)
    tseed!(rngs, 1)
    Threads.@threads for i = 1:length(x)
        x[i] = rand(rngs[Threads.threadid()])
    end
    tseed!(rngs, 1)
    Threads.@threads for i = 1:length(x)
        y[i] = rand(rngs[Threads.threadid()])
    end
    @test x == y
    states = map(rng -> rng.state, rngs)
    @test length(unique(states)) == Threads.nthreads()
end

@testset "nblasthreads/@with_blasthreads" begin
    BLAS.set_num_threads(4)
    @test nblasthreads() == 4
    BLAS.set_num_threads(1)
    @test nblasthreads() == 1
    @with_blasthreads 6 begin
        @test nblasthreads() == 6
    end
    @test nblasthreads() == 1
end

end


text:
using NearestNeighbors
using BenchmarkTools
using JLD
import Printf: @sprintf
import Random: MersenneTwister
import SpecialFunctions: gamma

include("generate_report.jl")

const EXTENSIVE_BENCHMARK = false

const SUITE = BenchmarkGroup()
SUITE["build tree"] = BenchmarkGroup()
SUITE["knn"] = BenchmarkGroup()
SUITE["inrange"] = BenchmarkGroup()

for n_points in (EXTENSIVE_BENCHMARK ? (10^3, 10^5) : 10^5)
    for dim in (EXTENSIVE_BENCHMARK ? (1, 3) : 3)
        data = rand(MersenneTwister(1), dim, n_points)
        for leafsize in (EXTENSIVE_BENCHMARK ? (1, 10) : 10)
            for reorder in (true, false)
                for (tree_type, SUITE_name) in ((KDTree, "kd tree"),
                                                (BallTree, "ball tree"))
                    tree = tree_type(data; leafsize = leafsize, reorder = reorder)
                    SUITE["build tree"]["$(tree_type) $dim × $n_points, ls = $leafsize"] = @benchmarkable $(tree_type)($data; leafsize = $leafsize, reorder = $reorder)
                    for input_size in (1, 1000)
                        input_data = rand(MersenneTwister(1), dim, input_size)
                        for k in (EXTENSIVE_BENCHMARK ? (1, 10) : 10)
                            SUITE["knn"]["$(tree_type) $dim × $n_points, ls = $leafsize, input_size = $input_size, k = $k"] = @benchmarkable knn($tree, $input_data, $k)
                        end
                        perc = 0.01
                        V = π^(dim / 2) / gamma(dim / 2 + 1) * (1 / 2)^dim
                        r = (V * perc * gamma(dim / 2 + 1))^(1/dim)
                        r_formatted = @sprintf("%3.2e", r)
                        SUITE["inrange"]["$(tree_type) $dim × $n_points, ls = $leafsize, input_size = $input_size, r = $r_formatted"] = @benchmarkable inrange($tree, $input_data, $r)
                    end
                end
            end
        end
    end
end

function run_benchmarks(name)
    paramspath = joinpath(dirname(@__FILE__), "params.jld")
    if !isfile(paramspath)
        println("Tuning benchmarks...")
        tune!(SUITE)
        JLD.save(paramspath, "SUITE", params(SUITE))
    end
    loadparams!(SUITE, JLD.load(paramspath, "SUITE"), :evals, :samples)
    results = run(SUITE, verbose = true, seconds = 2)
    JLD.save(joinpath(dirname(@__FILE__), name * ".jld"), "results", results)
end

function generate_report(v1, v2)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    v2_res = load(joinpath(dirname(@__FILE__), v2 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_compare.md"), "w") do f
        printreport(f, judge(minimum(v1_res), minimum(v2_res)); iscomparisonjob = true)
    end
end

function generate_report(v1)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_single.md"), "w") do f
        printreport(f, minimum(v1_res); iscomparisonjob = false)
    end
end

# run_benchmarks("primary")
# generate_report("primary") # generate a report with stats about a run
# run_benchmarks("secondary")
# generate_report("secondary", "primary") # generate report comparing two runs


text:
using LinearAlgebra, OrdinaryDiffEq
function foo(du, u, (A, tmp), t)
    tmp = DiffEqBase.get_tmp(tmp, u)
    mul!(tmp, A, u)
    @. du = u + tmp
    nothing
end
prob = ODEProblem(foo, ones(5, 5), (0., 1.0), (ones(5,5), DiffEqBase.dualcache(zeros(5,5))))
solve(prob, TRBDF2())


text:
module Terra

using NCDatasets
using DimensionalData
using Rasters

const NCD = NCDatasets
const DD = DimensionalData

include("Rasters.jl")

end # module


text:
JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didOpen")}}, params) = DidOpenTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/didOpen")},DidOpenTextDocumentParams}, server)
    server.isrunning = true
    uri = r.params.textDocument.uri
    if hasdocument(server, URI2(uri))
        doc = getdocument(server, URI2(uri))
        set_text!(doc, r.params.textDocument.text)
        doc._version = r.params.textDocument.version
        set_open_in_editor(doc, true)
        get_line_offsets(doc, true)
    else
        doc = Document(uri, r.params.textDocument.text, false, server)
        setdocument!(server, URI2(uri), doc)
        doc._version = r.params.textDocument.version
        doc._workspace_file = any(i->startswith(uri, filepath2uri(i)), server.workspaceFolders)
        doc._runlinter = !is_ignored(uri, server)
        set_open_in_editor(doc, true)

        try_to_load_parents(uri2filepath(uri), server)
    end
    parse_all(doc, server)
end

# JSONRPC.parse_params(::Type{Val{Symbol("julia/reloadText")}}, params) = DidOpenTextDocumentParams(params)
# function process(r::JSONRPC.Request{Val{Symbol("julia/reloadText")},DidOpenTextDocumentParams}, server)
#     server.isrunning = true
#     uri = r.params.textDocument.uri
#     if hasdocument(server, URI2(uri))
#         doc = getdocument(server, URI2(uri))
#         set_text!(doc, r.params.textDocument.text)
#         doc._version = r.params.textDocument.version
#     else
#         doc = Document(uri, r.params.textDocument.text, false, server)
#         setdocument!(server, URI2(uri), doc)
#         doc._version = r.params.textDocument.version
#         if any(i->startswith(uri, filepath2uri(i)), server.workspaceFolders)
#             doc._workspace_file = true
#         end
#         set_open_in_editor(doc, true)
#         if is_ignored(uri, server)
#             doc._runlinter = false
#         end
#     end
#     get_line_offsets(doc)
#     parse_all(doc, server)
# end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didClose")}}, params) = DidCloseTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/didClose")},DidCloseTextDocumentParams}, server)
    uri = r.params.textDocument.uri
    doc = getdocument(server, URI2(uri))
    empty!(doc.diagnostics)
    publish_diagnostics(doc, server)
    if is_workspace_file(doc)
        set_open_in_editor(doc, false)
    else
        if any(d.root == doc.root && d._open_in_editor for (uri, d::Document) in getdocuments_pair(server) if d != doc)
            # If any other open document shares doc's root we just mark it as closed...
            set_open_in_editor(doc, false)
        else
            # ...otherwise we delete all documents that share root with doc.
            for (u,d) in getdocuments_pair(server)
                if d.root == doc.root
                    deletedocument!(server, u)
                end
            end
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didSave")}}, params) = DidSaveTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/didSave")},DidSaveTextDocumentParams}, server)
    uri = r.params.textDocument.uri
    doc = getdocument(server, URI2(uri))
    if r.params.text isa String
        if get_text(doc) != r.params.text
            error("Mismatch between server and client text for $(doc._uri). _open_in_editor is $(doc._open_in_editor). _workspace_file is $(doc._workspace_file). _version is $(doc._version).")
        end
    end
    parse_all(doc, server)
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/willSave")}}, params) = WillSaveTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/willSave")},WillSaveTextDocumentParams}, server) end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/willSaveWaitUntil")}}, params) = WillSaveTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/willSaveWaitUntil")},WillSaveTextDocumentParams}, server)
    return TextEdit[]
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didChange")}}, params) = DidChangeTextDocumentParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/didChange")},DidChangeTextDocumentParams}, server::LanguageServerInstance)
    doc = getdocument(server, URI2(r.params.textDocument.uri))
    if r.params.textDocument.version < doc._version
        error("The client and server have different textDocument versions for $(doc._uri).")
    end
    doc._version = r.params.textDocument.version

    if length(r.params.contentChanges) == 1 && !endswith(doc._uri, ".jmd") && !ismissing(first(r.params.contentChanges).range)
        tdcce = first(r.params.contentChanges)
        new_cst = _partial_update(doc, tdcce)
        scopepass(getroot(doc), doc)
        StaticLint.check_all(getcst(doc), server.lint_options, server)
        empty!(doc.diagnostics)
        mark_errors(doc, doc.diagnostics)
        publish_diagnostics(doc, server)
    else
        for tdcce in r.params.contentChanges
            applytextdocumentchanges(doc, tdcce)
        end
        parse_all(doc, server)
    end
end

function _partial_update(doc::Document, tdcce::TextDocumentContentChangeEvent)
    cst = getcst(doc)
    insert_range = get_offset(doc, tdcce.range)
    updated_text = edit_string(get_text(doc), insert_range, tdcce.text)
    set_text!(doc, updated_text)
    doc._line_offsets = nothing

    i1, i2, loc1, loc2 = get_update_area(cst, insert_range)
    is = insert_size(tdcce.text, insert_range)
    if isempty(updated_text)
        empty!(cst.args)
        StaticLint.clear_meta(cst)
        cst.span = cst.fullspan = 0
    elseif 0 < i1 <= i2
        old_span = cst_len(cst, i1, i2)
        ps = ParseState(updated_text, loc1)
        args = EXPR[]
        if i1 == 1 && (ps.nt.kind == CSTParser.Tokens.WHITESPACE || ps.nt.kind == CSTParser.Tokens.COMMENT)
            CSTParser.next(ps)
            push!(args, CSTParser.mLITERAL(ps.nt.startbyte, ps.nt.startbyte, "", Tokens.NOTHING))
        else
            push!(args, CSTParser.parse(ps)[1])
        end
        while ps.nt.startbyte < old_span + loc1 + is && !ps.done
            push!(args, CSTParser.parse(ps)[1])
        end
        new_span = 0
        for i = 1:length(args)
            new_span += args[i].fullspan
        end
        # remove old blocks
        while old_span + is < new_span && i2 < length(cst.args)
            i2+=1
            old_span += cst.args[i2].fullspan
        end
        for i = i1:i2
            StaticLint.clear_meta(cst.args[i])
        end
        deleteat!(cst.args, i1:i2)

        #insert new blocks
        for a in args
            insert!(cst.args, i1, a)
            CSTParser.setparent!(cst.args[i1], cst)
            i1 += 1
        end
    else
        StaticLint.clear_meta(cst)
        cst = CSTParser.parse(updated_text, true)
    end
    CSTParser.update_span!(cst)
    doc.cst = cst
    if typof(doc.cst) === CSTParser.FileH
        doc.cst.val = doc.path
        set_doc(doc.cst, doc)
    end
end

insert_size(inserttext, insertrange) = sizeof(inserttext) - max(last(insertrange) - first(insertrange), 0) # OK, used to adjust EXPR spans

function cst_len(x, i1 = 1, i2 = length(x.args))
    n = 0
    @inbounds for i = i1:i2
        n += x.args[i].fullspan
    end
    n
end

function get_update_area(cst, insert_range)
    loc1 = loc2 = 0
    i1 = i2 = 0

    while i1 < length(cst.args)
        i1 += 1
        a = cst.args[i1]
        if loc1 <= first(insert_range) <= loc1 + a.fullspan
            loc2 = loc1
            i2 = i1
            if !(loc1 <= last(insert_range) <= loc1 + a.fullspan)
                while i2 < length(cst.args)
                    i2 += 1
                    a = cst.args[i2]
                    if loc2 <= last(insert_range) <= loc2 + a.fullspan
                        if i2 < length(cst.args) && last(insert_range) <= loc2 + a.fullspan
                            i2+=1
                        end
                        break
                    end
                    loc2 += a.fullspan
                end
            elseif i2 < length(cst.args) && last(insert_range) == loc1 + a.fullspan
                i2 += 1
            end
            break
        end
        loc1 += a.fullspan
    end
    return i1, i2, loc1, loc2
end

function applytextdocumentchanges(doc::Document, tdcce::TextDocumentContentChangeEvent)
    if ismissing(tdcce.range) && ismissing(tdcce.rangeLength)
        # No range given, replace all text
        set_text!(doc, tdcce.text)
    else
        editrange = get_offset(doc, tdcce.range)
        set_text!(doc, edit_string(get_text(doc), editrange, tdcce.text))
    end
    doc._line_offsets = nothing
end

function edit_string(text, editrange, edit)
    if first(editrange) == last(editrange) == 0
        text = string(edit, text)
    elseif first(editrange) == 0 && last(editrange) == sizeof(text) # OK, this branch is probably not needed
        text = edit
    elseif first(editrange) == 0
        text = string(edit, text[nextind(text, last(editrange)):end])
    elseif first(editrange) == last(editrange) == sizeof(text) # OK, this branch is probably not needed
        text = string(text, edit)
    elseif last(editrange) == sizeof(text) # OK, this branch is probably not needed
        text = string(text[1:first(editrange)], edit)
    else
        text = string(text[1:first(editrange)], edit, text[min(lastindex(text), nextind(text, last(editrange))):end])
    end
end

function parse_all(doc::Document, server::LanguageServerInstance)
    ps = CSTParser.ParseState(get_text(doc))
    StaticLint.clear_meta(getcst(doc))
    if endswith(doc._uri, ".jmd")
        doc.cst, ps = parse_jmd(ps, get_text(doc))
    else
        doc.cst, ps = CSTParser.parse(ps, true)
    end
    if typof(doc.cst) === CSTParser.FileH
        doc.cst.val = doc.path
        set_doc(doc.cst, doc)
    end

    scopepass(getroot(doc), doc)
    StaticLint.check_all(getcst(doc), server.lint_options, server)
    empty!(doc.diagnostics)
    mark_errors(doc, doc.diagnostics)

    publish_diagnostics(doc, server)
end

function mark_errors(doc, out = Diagnostic[])
    line_offsets = get_line_offsets(doc)
    errs = StaticLint.collect_hints(getcst(doc))
    n = length(errs)
    n == 0 && return out
    i = 1
    start = true
    offset = errs[i][1]
    r = Int[0, 0]
    pos = 0
    nlines = length(line_offsets)
    if offset > last(line_offsets)
        line = nlines
    else
        line = 1
        io = IOBuffer(get_text(doc))
        while line < nlines
            seek(io, line_offsets[line])
            char = 0
            while line_offsets[line] <= offset < line_offsets[line + 1]
                while offset > position(io)
                    c = read(io, Char)
                    if UInt32(c) >= 0x010000
                        char += 1
                    end
                    char += 1
                end
                if start
                    r[1] = line
                    r[2] = char
                    offset += errs[i][2].span
                else
                    if typof(errs[i][2]) === CSTParser.ErrorToken
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 1, "Julia", "Julia", "Parsing error", missing))
                    elseif typof(errs[i][2]) === CSTParser.IDENTIFIER && !StaticLint.haserror(errs[i][2])
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 2, "Julia", "Julia", "Missing reference: $(errs[i][2].val)", missing))
                    elseif StaticLint.haserror(errs[i][2]) && StaticLint.errorof(errs[i][2]) isa StaticLint.LintCodes
                        push!(out, Diagnostic(Range(r[1] - 1, r[2], line - 1, char), 3, "Julia", "Julia", get(StaticLint.LintCodeDescriptions, StaticLint.errorof(errs[i][2]), ""), missing))
                    end
                    i += 1
                    i>n && break
                    offset = errs[i][1]
                end
                start = !start
                offset = start ? errs[i][1] : errs[i][1] + errs[i][2].span
            end
            line += 1
        end
        close(io)
    end
    return out
end

function publish_diagnostics(doc::Document, server)
    if server.runlinter && server.symbol_store_ready && is_workspace_file(doc)
        publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, doc._version, doc.diagnostics)
    else
        publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, doc._version, Diagnostic[])
    end
    JSONRPCEndpoints.send_notification(server.jr_endpoint, "textDocument/publishDiagnostics", publishDiagnosticsParams)
end

function clear_diagnostics(uri::URI2, server)
    doc = getdocument(server, uri)
    empty!(doc.diagnostics)
    publishDiagnosticsParams = PublishDiagnosticsParams(doc._uri, doc._version, Diagnostic[])
    JSONRPCEndpoints.send_notification(server.jr_endpoint, "textDocument/publishDiagnostics", publishDiagnosticsParams)
end

function clear_diagnostics(server)
    for uri in getdocuments_key(server)
        clear_diagnostics(uri, server)
    end
end

function parse_jmd(ps, str)
    currentbyte = 1
    blocks = []
    while ps.nt.kind != Tokens.ENDMARKER
        CSTParser.next(ps)
        if ps.t.kind == Tokens.CMD || ps.t.kind == Tokens.TRIPLE_CMD
            push!(blocks, (ps.t.startbyte, CSTParser.INSTANCE(ps)))
        end
    end
    top = EXPR(CSTParser.FileH, EXPR[])
    if isempty(blocks)
        return top, ps
    end

    for (startbyte, b) in blocks
        if typof(b) === CSTParser.LITERAL && kindof(b) == CSTParser.Tokens.TRIPLE_CMD && (startswith(b.val, "julia") || startswith(b.val, "{julia"))
            blockstr = b.val
            ps = CSTParser.ParseState(blockstr)
            # skip first line
            while ps.nt.startpos[1] == 1
                CSTParser.next(ps)
            end
            prec_str_size = currentbyte:startbyte + ps.nt.startbyte + 3

            push!(top.args, CSTParser.mLITERAL(length(prec_str_size), length(prec_str_size), "", CSTParser.Tokens.STRING))

            args, ps = CSTParser.parse(ps, true)
            append!(top.args, args.args)
            CSTParser.update_span!(top)
            currentbyte = top.fullspan + 1
        elseif typof(b) === CSTParser.LITERAL && kindof(b) == CSTParser.Tokens.CMD && startswith(b.val, "j ")
            blockstr = b.val
            ps = CSTParser.ParseState(blockstr)
            CSTParser.next(ps)
            prec_str_size = currentbyte:startbyte + ps.nt.startbyte + 1
            push!(top.args, CSTParser.mLITERAL(length(prec_str_size), length(prec_str_size), "", CSTParser.Tokens.STRING))

            args, ps = CSTParser.parse(ps, true)
            append!(top.args, args.args)
            CSTParser.update_span!(top)
            currentbyte = top.fullspan + 1
        end
    end

    prec_str_size = currentbyte:sizeof(str) # OK
    push!(top.args, CSTParser.mLITERAL(length(prec_str_size), length(prec_str_size), "", CSTParser.Tokens.STRING))

    return top, ps
end

function search_for_parent(dir::String, file::String, drop = 3, parents = String[])
    drop<1 && return parents
    !isdir(dir) && return parents
    !hasreadperm(dir) && return parents
    for f in readdir(dir)
        filename = joinpath(dir, f)
        if isvalidjlfile(filename)
            # Could be sped up?
            content = try
                s = read(filename, String)
                # We throw an error in the case of an invalid
                # UTF-8 sequence so that the same code path
                # is used that handles file IO problems
                isvalid(s) || error()
                s
            catch err
                continue
            end
            occursin(file, content) && push!(parents, joinpath(dir, f))
        end
    end
    search_for_parent(splitdir(dir)[1], file, drop - 1, parents)
    return parents
end

function is_parentof(parent_path, child_path, server)
    !isvalidjlfile(parent_path) && return false
    previous_server_docs = collect(getdocuments_key(server)) # additions to this to be removed at end
    # load parent file
    puri = filepath2uri(parent_path)
    if !hasdocument(server, URI2(puri))
        content = try
            s = read(parent_path, String)
            # We throw an error in the case of an invalid
            # UTF-8 sequence so that the same code path
            # is used that handles file IO problems
            isvalid(s) || error()
            s
        catch err
            return false
        end
        pdoc = Document(puri, content, false, server)
        setdocument!(server, URI2(puri), pdoc)
        CSTParser.parse(get_text(pdoc))
        if typof(pdoc.cst) === CSTParser.FileH
            pdoc.cst.val = pdoc.path
            set_doc(pdoc.cst, pdoc)
        end
    else
        pdoc = getdocument(server, URI2(puri))
    end
    scopepass(getroot(pdoc), pdoc)
    # check whether child has been included automatically
    if any(uri2filepath(k._uri) == child_path for k in getdocuments_key(server) if !(k in previous_server_docs))
        cdoc = getdocument(server,URI2(filepath2uri(child_path)))
        parse_all(cdoc, server)
        scopepass(getroot(cdoc))
        return true
    else
        # clean up
        foreach(k-> !(k in previous_server_docs) && deletedocument!(server, k), getdocuments_key(server))
        return false
    end
end

function try_to_load_parents(child_path, server)
    for p in search_for_parent(splitdir(child_path)...)
        p == child_path && continue
        success = is_parentof(p, child_path, server)
        if success
            return try_to_load_parents(p, server)
        end
    end
end


text:
"""
    bbands(x::Array{Float64}; n::Int64=10, sigma::Float64=2.0)::Matrix{Float64}

Bollinger bands (moving average with standard deviation bands)

*Output*
- Column 1: lower band
- Column 2: middle band
- Column 3: upper band
"""
function bbands(x::Array{Float64}; n::Int64=10, sigma::Float64=2.0, ma::Function=sma, args...)::Matrix{Float64}
    @assert n<size(x,1) && n>0 "Argument n is out of bounds."
    out = zeros(size(x,1), 3)  # cols := lower bound, ma, upper bound
    out[:,2] = ma(x, n=n, args...)
    sd = runsd(x, n=n, cumulative=false)
    out[:,1] = out[:,2] - sigma*sd
    out[:,3] = out[:,2] + sigma*sd
    return out
end

"""
    tr(hlc::Matrix{Float64})::Array{Float64}

True range
"""
function tr(hlc::Matrix{Float64})::Array{Float64}
    @assert size(hlc,2) == 3 "HLC array must have 3 columns."
    n = size(hlc,1)
    out = zeros(n)
    out[1] = NaN
    @inbounds for i=2:n
        out[i] = max(hlc[i,1]-hlc[i,2], hlc[i,1]-hlc[i-1,3], hlc[i-1,3]-hlc[i,2])
    end
    return out[:,1]
end

"""
    atr(hlc::Matrix{Float64}; n::Int64=14)::Array{Float64}

Average true range (uses exponential moving average)
"""
function atr(hlc::Matrix{Float64}; n::Int64=14, ma::Function=ema)::Array{Float64}
    @assert n<size(hlc,1) && n>0 "Argument n out of bounds."
    return [NaN; ma(tr(hlc)[2:end], n=n)]
end

"""
    keltner(hlc::Matrix{Float64}; nema::Int64=20, natr::Int64=10, mult::Int64=2)::Matrix{Float64}

Keltner bands

*Output*
Column 1: lower band
Column 2: middle band
Column 3: upper band
"""
function keltner(hlc::Array{Float64,2}; nema::Int64=20, natr::Int64=10, mult::Int64=2)::Matrix{Float64}
    @assert size(hlc,2) == 3 "HLC array must have 3 columns."
    out = zeros(size(hlc,1), 3)
    out[:,2] = ema(hlc[:,3], n=nema)
    out[:,1] = out[:,2] - mult*atr(hlc, n=natr)
    out[:,3] = out[:,2] + mult*atr(hlc, n=natr)
    return out
end


text:
"""
This is the parent type for all basic errors during ray tracing.
The type field describes the type of error.
If the maximum number of iterations to find a surface intersection
was too large, n_iterations contains the number of iterations.
"""
struct RayError
    type::Symbol
    n_iterations::Integer
end

function ray_miss_error()
    RayError(:ray_miss, 0)
end

function total_internal_reflection_error()
    RayError(:total_internal_reflection, 0)
end

function intersection_max_iterations_error()
    RayError(:intersection_max_iterations, 0)
end

function vignetted_error()
    RayError(:vignetted, 0)
end

iserror(x::RayError) = true
iserror(x) = false

export iserror

"""
This is the type returned from a ray trace when an error occured
error_type is the actual RayError with any additional specific information
global_index is the index of the surface that the error occured at
data is the ray trace result up to this point
"""
struct RaytraceError{D} <: Exception
    error_type::RayError
    global_index::Integer # global surface indcex from start of lens; object surface is index 1
    component_id::Symbol
    local_index::Integer # surface index within the component
    surface_id::Union{Nothing,Symbol}
    data::D
end

function iserror(x::RaytraceError)
    true
end

error_type(err::RaytraceError) = err.error_type.type


text:
# ---
# title: 1509. Minimum Difference Between Largest and Smallest Value in Three Moves
# id: problem1509
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/minimum-difference-between-largest-and-smallest-value-in-three-moves/description/>
# hidden: true
# ---
#
# Given an array `nums`, you are allowed to choose one element of `nums` and
# change it by any value in one move.
#
# Return the minimum difference between the largest and smallest value of `nums`
# after perfoming at most 3 moves.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,3,2,4]
#     Output: 0
#     Explanation: Change the array [5,3,2,4] to [ **2** , **2** ,2, **2** ].
#     The difference between the maximum and minimum is 2-2 = 0.
#
# **Example 2:**
#
#
#
#     Input: nums = [1,5,0,10,14]
#     Output: 1
#     Explanation: Change the array [1,5,0,10,14] to [1, **1** ,0, **1** , **1** ].
#     The difference between the maximum and minimum is 1-0 = 1.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [6,6,0,1,1,4,6]
#     Output: 2
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,5,6,14,15]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `-10^9 <= nums[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

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


text:
@testset "Manopt Cyclic Proximal Point" begin
    using Dates
    n = 100
    N = PowerManifold(Circle(),n)
    f = artificial_S1_signal(n)
    F = x -> costL2TV(N,f,0.5,x)
    proxes = [ (λ,x) -> prox_distance(N,λ,f,x), (λ,x) -> prox_TV(N,0.5*λ,x) ]
    o = cyclic_proximal_point(
        N,
        F,
        proxes,
        f;
        λ = i -> π/(2*i),
        stopping_criterion = StopAfterIteration(100),
        debug = [DebugIterate()," ",DebugCost()," ",DebugProximalParameter(),"\n",10000],
        record = [RecordProximalParameter(), RecordIterate(f), RecordCost()],
        return_options=true
        )
    fR = get_solver_result(o)
    fR2 = cyclic_proximal_point(
        N,
        F,
        proxes,
        f;
        λ = i -> π/(2*i),
        stopping_criterion = StopAfterIteration(100),
    )
    @test fR == fR2
    rec = get_record(o)
    @test F(f) > F(fR)
    #
    o = CyclicProximalPointOptions(f, StopAfterIteration(1), i -> π/(2*i))
    p = ProximalProblem(N,F,proxes,[1,2])
    @test_throws ErrorException getProximalMap(p,1.,f,3)

    order1 = Manopt.update_cpp_order(3,0,1:3,RandomEvalOrder())
    @test all(isinteger.(order1))
    @test minimum(order1) == 1
    @test maximum(order1) == 3

    # i=0 introduce new fixed ranom order
    order2 = Manopt.update_cpp_order(3,0,1:3,FixedRandomEvalOrder())
    @test all(isinteger.(order2))
    @test minimum(order2) == 1
    @test maximum(order2) == 3
    order3 = Manopt.update_cpp_order(3,1,order2,FixedRandomEvalOrder())
    @test order3 == order2 # only update on i=0

    @test Manopt.update_cpp_order(3,0,1:3,LinearEvalOrder()) == 1:3
    @test Manopt.update_cpp_order(3,1,1:3,LinearEvalOrder()) == 1:3
end


text:
# Generate random demand and covariate data for the resource allocation model

# Function to generate random correlation matrices. Based on https://stats.stackexchange.com/questions/124538/how-to-generate-a-large-full-rank-random-correlation-matrix-with-some-strong-cor
function generateRandomCorrMat(dim::Int64)

	betaparam::Float64 = 2.0

	partCorr = zeros(Float64,dim,dim)
	corrMat = eye(dim)

	for k = 1:dim-1
		for i = k+1:dim
			partCorr[k,i] = ((rand(Beta(betaparam,betaparam),1))[1] - 0.5)*2.0
			p::Float64 = partCorr[k,i]
			for j = (k-1):-1:1
				p = p*sqrt((1-partCorr[j,i]^2)*(1-partCorr[j,k]^2)) + partCorr[j,i]*partCorr[j,k]
			end
			corrMat[k,i] = p
			corrMat[i,k] = p
		end
	end

	permut = randperm(dim)
	corrMat = corrMat[permut, permut]

	return corrMat
end

# generate base parameters for the demand model given number of covariates
# generates a common set of coefficients across all replicates
function generateBaseDemandModel(numCovariates::Int64)

	# generate parameters for the covariate distribution
	covariate_mean = zeros(Float64, maxNumCovariates)
	covariate_corrMat = generateRandomCorrMat(maxNumCovariates)
	covariate_covMat = covariate_corrMat + 1E-09*eye(maxNumCovariates)

	# generate coefficients for the dependence of the demands on the covariates
	alpha = 50.0 + 5.0*rand(Normal(), numCustomers)
	beta = repeat([10.0,5.0,2.0],outer=[1,numCustomers])' + rand(Uniform(-4.0,4.0), numCustomers, 3)

	if(numCovariates < 3)
		throw(ErrorException("Expected numCovariates >= 3!"))
	end

	coeff_true = zeros(Float64,numCovariates+1,numCustomers)
	coeff_true[1,:] = alpha
	coeff_true[2,:] = beta[:,1]
	coeff_true[3,:] = beta[:,2]
	coeff_true[4,:] = beta[:,3]

	return covariate_mean, covariate_covMat, coeff_true
end

# generate demand and covariate data using a sparse (non)linear model
function generateDemandData(numCovariates::Int64,numSamples::Int64,degree::Float64,covariate_mean::Array{Float64},covariate_covMat::Array{Float64,2},coeff_true::Array{Float64,2})

	covariate_mean_case = covariate_mean[1:numCovariates]
	covariate_covMat_case = covariate_covMat[1:numCovariates,1:numCovariates]
	covariate_covMat_chol = (cholfact(Hermitian(covariate_covMat_case)))[:L]

	#*===========================================
	# first, construct the the covariates
	# first covariate is simply the constant one (for the intercept)
	# the next set of covariates are iid samples from a multivariate folded normal distribution

	covariate_data = zeros(Float64, numSamples, numCovariates+1)

	covariate_data[:,1] = ones(Float64,numSamples)
	covariate_tmp = rand(Normal(),maxNumCovariates,maxNumDataSamples)
	for s = 1:numSamples
		covariate_data[s,2:numCovariates+1] = covariate_covMat_chol*covariate_tmp[1:numCovariates,s] + covariate_mean_case
	end
	covariate_data = abs.(covariate_data)

	#*===========================================
	# next, construct the random demands from the covariates using a sparse (non)linear model

	demand_errors = demand_errors_scaling*rand(Normal(), maxNumDataSamples, numCustomers)
	demand_data = ((covariate_data).^(degree))*coeff_true + demand_errors[1:numSamples,:]

	return demand_data, covariate_data
end

# generate a new realization of the covariates from a multivariate folded normal distribution
function generateCovariateReal(numCovariates::Int64,covariate_mean::Array{Float64},covariate_covMat::Array{Float64,2})

	covariate_mean_case = covariate_mean[1:numCovariates]
	covariate_covMat_case = covariate_covMat[1:numCovariates,1:numCovariates]
	covariate_covMat_chol = (cholfact(Hermitian(covariate_covMat_case)))[:L]

	covariate_tmp = rand(Normal(),maxNumCovariates)
	covariate_obs = ones(Float64, numCovariates+1)
	covariate_obs[2:numCovariates+1] = covariate_covMat_chol*(covariate_tmp[1:numCovariates]) + covariate_mean_case
	covariate_obs = abs.(covariate_obs)

	return covariate_obs
end

# generate samples from the true conditional distribution of the demands given a covariate realization
function generateTrueCondScenarios(numScenarios::Int64,covariate_obs::Array{Float64,1},degree::Float64,coeff_true::Array{Float64,2})

	demand_errors = demand_errors_scaling*rand(Normal(), numScenarios, numCustomers)
	demand_scen = repeat((((covariate_obs).^(degree))'*coeff_true)', outer = [1,numScenarios])' + demand_errors

	return demand_scen
end

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
export Puzzle

struct Puzzle
    board::Matrix{Int} # n x n matrix
    function Puzzle(board::Matrix{Int})
        # verify input is a valid Sudoku puzzle (9 x 9)
        @assert size(board) == (9, 9) "Invalid puzzle dimensions!"
        return new(board);
    end
end

function Base.getindex(p::Puzzle, i::Int, j::Int)::Int
    rows, cols = size(p.board)
    if (1 ≤ i ≤ rows) && (1 ≤ j ≤ cols)
        return p.board[i,j];
    end
    throw(BoundsError(p));
end

function Base.setindex!(p::Puzzle, n::Int, i::Int, j::Int)
    rows, cols = size(p.board)
    if (1 ≤ i ≤ rows) && (1 ≤ j ≤ cols)
        p.board[i,j] = n;
    else
        throw(BoundsError(p));
    end
end

function Base.show(io::IO, p::Puzzle)::Nothing
    println(io, "\033[2J") # clear repl
    println(io, "\033[$(displaysize(stdout)[1])A")
    rows, cols = size(p.board)
    for i in 1:rows
        for j in 1:cols
            print(io, " $(p.board[i,j]) ")
            if (j % 3 == 0) && (j != cols)
                print(io, "|")
            end
        end
        println(io)
        if (i % 3 == 0) && (i != rows)
            # border lines
            println(io, ("-"^9 * "+")^2 * "-"^9)
        end
    end
end

max_stars_repo_name: Joshua-Robison/Sudoku

max_issues_repo_name: Joshua-Robison/Sudoku

max_forks_repo_name: Joshua-Robison/Sudoku


text:
type MedianPAC<:FiniteHorizonAgent
    policy::Matrix{UInt64}
    ns::Array{UInt64, 5}
    nsactive::Array{UInt64, 5}
    nactive::Array{UInt64, 3}
    nsnew::Array{UInt64, 4}
    cbucket::Array{UInt64, 3}
    #n::Array{UInt64, 3}
    rewards::Array{Float64, 3}
    #Vopt::Matrix{Float64}
    #δ::Float64
    ϵa::Float64
    maxobs::Float64
    maxRet::Float64
    maxR::Float64
    keep_updating::Bool
    trash_samples::Bool
end

function MedianPAC(S, A, H, rewards, δ, ϵ, keep_updating=false, trash_samples=false)
    policy = zeros(UInt64, S, H)
    maxRet = sum(maximum(reshape(rewards, S*A, H), 1))
    # set to satisfy ϵ0 = ϵa / H in Theorem 5.4
    ϵa = 1/2*H*ϵ
    km = ceil(Int64, 50 / 9 * log(4 / δ * 2 * log2(2 * maxRet / ϵa * S * A * H)))
    nbuckets = km
    ns = zeros(UInt64, S, nbuckets, S, A, H)
    nsnew = zeros(UInt64, nbuckets, S, A, H)
    nsactive = zeros(UInt64, S, nbuckets, S, A, H)
    nactive = zeros(UInt64, S, A, H)
    cbucket = ones(UInt64, S, A, H)

    # k in the paper
    maxobs = km * 2^ceil(Int64,log2(4 * maxRet^2 / ϵa^2))

    m = MedianPAC(policy, ns, nsactive, nactive, nsnew, cbucket, rewards, ϵa, maxobs,
                  maxRet, maximum(rewards[:]), keep_updating, trash_samples)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end

nbuck(m::MedianPAC) = size(m.ns, 2)
function observe!(m::MedianPAC, s, a, t, r, sn)
    if !m.keep_updating && m.nactive[s,a,t] >= m.maxobs
        return
    end
    m.nsnew[m.cbucket[s, a, t], s, a, t] = sn
    m.cbucket[s,a ,t] += 1
    if m.cbucket[s, a, t] > nbuck(m)
        m.cbucket[s,a,t] = 1
        for i ∈ 1:nbuck(m)
            m.ns[m.nsnew[i, s, a, t], i, s, a, t] += 1
        end
        if sum(m.ns[:,1,s,a,t]) >= max(1, 2 * m.nactive[s, a, t])
            # new batch of twice the size, replace existing batch
            m.nsactive[:,:, s,a ,t] = m.ns[:,:,s, a, t]
            m.nactive[s, a, t] = sum(m.ns[:,1, s, a, t])
            if m.trash_samples
                m.ns[:,:,s,a,t] = 0
            end
        end
    end
end
maxV(m::MedianPAC) = m.maxRet
maxR(m::MedianPAC) = m.maxR

function update_policy!(m::MedianPAC)
    S = nS(m)
    H = horizon(m)
    A = nA(m)
    ϵb = m.ϵa * sqrt(m.maxobs)
    Vmax = maxV(m)
    Rmax = maxR(m)

    Q = zeros(A)
    V = zeros(S)
    Vnew = zeros(S)
    means = zeros(nbuck(m))
    for t=H:-1:1
        V, Vnew = Vnew, V
        curmaxV = min(maximum(V) + Rmax, Vmax)
        for s ∈ 1:S
            for a ∈ 1:A
                if m.nactive[s,a, t] == 0
                    Q[a] = curmaxV
                    break
                end
                slack = ϵb / √(m.nactive[s, a, t])
                for i ∈ 1:nbuck(m)
                    means[i] = V ⋅ m.nsactive[:,i, s, a, t] / m.nactive[s, a, t]
                end
                Q[a] = min(curmaxV, slack + m.rewards[s, a, t] + median(means))
            end
            bestA = indmax(Q)
            m.policy[s, t] = bestA
            Vnew[s] = Q[bestA]
        end
    end
end


text:
using FromFile
@from "Core.jl" import CONST_TYPE, Node, copyNode, Options

# Count the operators, constants, variables in an equation
function countNodes(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countNodes(tree.l)
    else
        return 1 + countNodes(tree.l) + countNodes(tree.r)
    end
end

# Count the max depth of a tree
function countDepth(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countDepth(tree.l)
    else
        return 1 + max(countDepth(tree.l), countDepth(tree.r))
    end
end

function stringOp(op::F, tree::Node, options::Options;
                  bracketed::Bool=false,
                  varMap::Union{Array{String, 1}, Nothing}=nothing)::String where {F}
    if op in [+, -, *, /, ^]
        l = stringTree(tree.l, options, bracketed=false, varMap=varMap)
        r = stringTree(tree.r, options, bracketed=false, varMap=varMap)
        if bracketed
            return "$l $(string(op)) $r"
        else
            return "($l $(string(op)) $r)"
        end
    else
        l = stringTree(tree.l, options, bracketed=true, varMap=varMap)
        r = stringTree(tree.r, options, bracketed=true, varMap=varMap)
        return "$(string(op))($l, $r)"
    end
end

"""
    stringTree(tree::Node, options::Options; kws...)

Convert an equation to a string.

# Arguments

- `varMap::Union{Array{String, 1}, Nothing}=nothing`: what variables
    to print for each feature.
"""
function stringTree(tree::Node, options::Options;
                    bracketed::Bool=false,
                    varMap::Union{Array{String, 1}, Nothing}=nothing)::String
    if tree.degree == 0
        if tree.constant
            return string(tree.val)
        else
            if varMap == nothing
                return "x$(tree.feature)"
            else
                return varMap[tree.feature]
            end
        end
    elseif tree.degree == 1
        return "$(options.unaops[tree.op])($(stringTree(tree.l, options, bracketed=true, varMap=varMap)))"
    else
        return stringOp(options.binops[tree.op], tree, options, bracketed=bracketed, varMap=varMap)
    end
end

# Print an equation
function printTree(tree::Node, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    println(stringTree(tree, options, varMap=varMap))
end

# Count the number of unary operators in the equation
function countUnaryOperators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 1 + countUnaryOperators(tree.l)
    else
        return 0 + countUnaryOperators(tree.l) + countUnaryOperators(tree.r)
    end
end

# Count the number of binary operators in the equation
function countBinaryOperators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 0 + countBinaryOperators(tree.l)
    else
        return 1 + countBinaryOperators(tree.l) + countBinaryOperators(tree.r)
    end
end

# Count the number of operators in the equation
function countOperators(tree::Node)::Int
    return countUnaryOperators(tree) + countBinaryOperators(tree)
end

# Count the number of constants in an equation
function countConstants(tree::Node)::Int
    if tree.degree == 0
        if tree.constant
            return 1
        else
            return 0
        end
    elseif tree.degree == 1
        return 0 + countConstants(tree.l)
    else
        return 0 + countConstants(tree.l) + countConstants(tree.r)
    end
end

# Get all the constants from a tree
function getConstants(tree::Node)::AbstractVector{CONST_TYPE}
    if tree.degree == 0
        if tree.constant
            return [tree.val]
        else
            return CONST_TYPE[]
        end
    elseif tree.degree == 1
        return getConstants(tree.l)
    else
        both = [getConstants(tree.l), getConstants(tree.r)]
        return [constant for subtree in both for constant in subtree]
    end
end

# Set all the constants inside a tree
function setConstants(tree::Node, constants::AbstractVector{T}) where {T<:Real}
    if tree.degree == 0
        if tree.constant
            tree.val = convert(CONST_TYPE, constants[1])
        end
    elseif tree.degree == 1
        setConstants(tree.l, constants)
    else
        numberLeft = countConstants(tree.l)
        setConstants(tree.l, constants)
        setConstants(tree.r, constants[numberLeft+1:end])
    end
end


text:
"""
    Model for testing steady state solver with linear growth variables.
"""
module E6

using ModelBaseEcon

model = Model()

@parameters model begin
    p_dlp = 0.0050000000000000
    p_dly = 0.0045000000000000
end

@variables model begin
    dlp; dly; dlyn; lp; ly; lyn
end

@shocks model begin
    dlp_shk; dly_shk
end

@autoexogenize model begin
    ly = dly_shk
    lp = dlp_shk
end

@equations model begin
    dly[t]=(1-0.2-0.2)*p_dly+0.2*dly[t-1]+0.2*dly[t+1]+dly_shk[t]
    dlp[t]=(1-0.5)*p_dlp+0.1*dlp[t-2]+0.1*dlp[t-1]+0.1*dlp[t+1]+0.1*dlp[t+2]+0.1*dlp[t+3]+dlp_shk[t]
    dlyn[t]=dly[t]+dlp[t]
    ly[t]=ly[t-1]+dly[t]
    lp[t]=lp[t-1]+dlp[t]
    lyn[t]=lyn[t-1]+dlyn[t]
 end

@initialize model

end


text:
import LinearAlgebra: Symmetric

"""
NumDiffTools.hessian(f, x::Real, [abstol, reltol, h, t, maxiters])

Returns `H(f)`, i.e., the Hessian, evaluted at `x`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function hessian(f, x::Real; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x);

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 2:2:2*maxiters

    extrapolate(h->(f(x + h) + f(x - h) - 2*fx)/h^2,             # Function to extrapolate
                A1, A2,                                          # Work Arrays
                hh, k,                                           # Sequence of step sizes (h)
                abstol, reltol                                   # tolerances
                )
end

"""
NumDiffTools.hessian(f, x::AbstractArray, [abstol, reltol, h, t, maxiters])

Returns `H(f)`, i.e., the Hessian, evaluted at `x`, assuming `f` is called as `f(x)`.
"""
function hessian(f, x::AbstractArray; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x);

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 2:2:2*maxiters
    result = Matrix{typeof(fx)}(undef, length(x), length(x))
    dx1 = zero(x)
    dx2 = zero(x)
    for i = 1:length(x)
        # Diagonal element
        dx1[i] = 1
        result[i,i] =  extrapolate(h->(f(x + h*dx1) + f(x - h*dx1) - 2*fx)/h^2,             # Function to extrapolate
                A1, A2,                                                                     # Work Arrays
                hh, k,                                                                      # Sequence of step sizes (h)
                abstol, reltol                                                              # tolerances
                )
        for j = (i+1):length(x)
            dx2[j] = 1
            result[i,j] =  extrapolate(h->(f(x + h*dx1 + h*dx2) - f(x + h*dx1 - h*dx2) - f(x - h*dx1 + h*dx2) + f(x - h*dx1 - h*dx2))/(4*h^2),
                                       A1, A2,                                                                     # Work Arrays
                                       hh, k,                                                                      # Sequence of step sizes (h)
                                       abstol, reltol                                                              # tolerances
                                       )
            result[j,i] = result[i,j]
            dx2[j] = 0
        end
        dx1[i] = 0
    end
    result
end

"""
NumDiffTools.hessian!(result::AbstractArray, f, x::AbstractArray, [abstol, reltol, h, t, maxiters])

Stores `H(f)`, i.e., the Hessian, evaluted at `x` in `result`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function hessian!(result::AbstractArray, f, x::AbstractArray; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x);

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 2:2:2*maxiters

    dx1 = zero(x)
    dx2 = zero(x)
    for i = 1:length(x)
        # Diagonal element
        dx1[i] = 1
        result[i,i] =  extrapolate(h->(f(x + h*dx1) + f(x - h*dx1) - 2*fx)/h^2,             # Function to extrapolate
                A1, A2,                                                                     # Work Arrays
                hh, k,                                                                      # Sequence of step sizes (h)
                abstol, reltol                                                              # tolerances
                )
        for j = (i+1):length(x)
            dx2[j] = 1
            result[i,j] =  extrapolate(h->(f(x + h*dx1 + h*dx2) - f(x + h*dx1 - h*dx2) - f(x - h*dx1 + h*dx2) + f(x - h*dx1 - h*dx2))/(4*h^2),
                                       A1, A2,                                                                     # Work Arrays
                                       hh, k,                                                                      # Sequence of step sizes (h)
                                       abstol, reltol                                                              # tolerances
                                       )
            result[j,i] = result[i,j]
            dx2[j] = 0
        end
        dx1[i] = 0
    end
end


text:
"""
    jupytertolatex(notebook, targetdir="./build_latex"; template=:book, fontpath=nothing)
Function to convert Jupyter notebooks. The arguments are the same as the ones in
`notebooktolatex`.
"""
function jupytertolatex(notebook, targetdir="./build_latex"; template=:book, fontpath=nothing)

    createproject(targetdir, template, fontpath)

    notebookname = basename(notebook)[1:end-6]
    notebookdir  = dirname(notebook)*"/"
    jsonnb = JSON.parse(read(notebook, String))
    texfile = read(targetdir*"/main.tex", String)
    lineinsert = 1
    for (i,line) in enumerate(split(texfile, "\n"))
        if startswith(line, "% INCLUDE NOTEBOOKS")
            lineinsert = i
            break
        end
    end

    if !occursin("\\include{./notebooks/"*notebookname*"}",read(targetdir*"/main.tex", String))

        insertlinebelow(targetdir*"/main.tex",
            "\\include{./notebooks/"*notebookname*"}", lineinsert)
    end

    notebook    = targetdir*"/notebooks/"*notebookname*".tex"
    figureindex = Dict(:i=>0)
    open(notebook, "w") do f
        write(f,"\\newpage\n")
        for cell in jsonnb["cells"]

            # Checks whether the cell has markdown
            if get(cell,"cell_type", nothing) == "markdown" || get(cell,"cell_type", nothing) == "raw"
                parsed = markdowntolatex(strip(join(cell["source"])), targetdir, notebookdir, template=template)
                write(f,parsed)

            # Checks whether the cell has code and whether the code is hidden
            elseif get(cell,"cell_type", nothing) == "code" && nestedget(cell,["metadata","jupyter", "source_hidden"],nothing) === nothing
                if get(cell,"outputs", nothing) != []
                    write(f,"\n\\bigskip")
                    # The textcl=true option is for using words with accent in the comments. Useful for non-english.
                    write(f,"\n\\begin{lstlisting}[language=JuliaLocal, style=julia, texcl=true]\n")
                    write(f, strip(join(cell["source"])))
                    write(f,"\n\\end{lstlisting}\n")
                end
            end

            ## Collecting outputs

            # Checks if the output is an empty array or if the outputs_hidden is true
            check_output_hidden = get(cell, "outputs",[]) == [] || nestedget(cell, ["metadata","jupyter","outputs_hidden"], nothing) !== nothing

            # Collect the output if the cell has code and the output is not hidden
            if get(cell, "cell_type", nothing) == "code" && !check_output_hidden
                for output in get(cell, "outputs", nothing)
                    if get(output, "output_type", nothing) == "stream"
                        write(f,"\n\\begingroup\n")
                        write(f,"\n\\fontsize{10pt}{12pt}\\selectfont\n")
                        write(f,"\n\\begin{verbatim}\n")
                        write(f, join(output["text"]))
                        write(f,"\n\\end{verbatim}\n")
                        write(f,"\n\\endgroup\n")
                    elseif get(output, "output_type", nothing) == "execute_result"
                        if nestedget(output,["data","text/latex"], nothing) !== nothing
                            write(f, "\n"*join(output["data"]["text/latex"]))
                        elseif nestedget(output,["data","text/plain"], nothing) !== nothing
                            write(f,"\n\\begingroup\n")
                            write(f,"\n\\fontsize{10pt}{12pt}\\selectfont\n")
                            write(f,"\n\\begin{verbatim}\n")
                            write(f, "\n"*join(output["data"]["text/plain"]))
                            write(f,"\n\\end{verbatim}\n")
                            write(f,"\n\\endgroup\n")
                        end
                        if nestedget(output,["data","image/png"], nothing) !== nothing
                            png = base64decode(output["data"]["image/png"])

                            figureindex[:i]+=1
                            figurename = notebookname*"_figure"*string(figureindex[:i])*".png"
                            write(targetdir*"/figures/"*figurename, png)
                            write(f,"\n\\begin{figure}[H]\n")
                            write(f,"\t\\centering\n")
                            write(f,"\t\\includegraphics[width=0.8\\textwidth]{./figures/"*figurename*"}\n")
                            write(f,"\t\\label{fig:"*figurename*"}\n")
                            write(f,"\n\\end{figure}\n")
                        end
                        if nestedget(output,["data","image/svg+xml"], nothing) !== nothing
                            svg = join(output["data"]["image/svg+xml"])
                            figureindex[:i]+=1
                            figurename = notebookname*"_figure"*string(figureindex[:i])
                            figuresvg = targetdir*"/figures/"*figurename*".svg"

                            # Save svg figure
                            write(figuresvg, svg)

                            # Convert svg to pdf
                            figurepdf = targetdir*"/figures/"*figurename*".pdf"
                            rsvg_convert() do cmd
                               run(`$cmd $figuresvg -f pdf -o $figurepdf`)
                            end

                            # Parse pdf image to Latex
                            write(f,"\n\\begin{figure}[H]\n")
                            write(f,"\t\\centering\n")
                            write(f,"\t\\includegraphics[width=0.8\\textwidth]{./figures/"*figurename*".pdf}\n")
                            write(f,"\t\\label{fig:"*figurename*"}\n")
                            write(f,"\n\\end{figure}\n")
                        end
                    end
                end
            end
        end
    end
end


text:
using Test
using SafeTestsets

include("generate_example_tests.jl")

# Note: comment outer @testset to stop after first @safetestset failure
@time @testset verbose = true "GRAPE.jl Package" begin

    print("\n* Example 1 (examples/simple_state_to_state.jl):")
    @time @safetestset "Example 1" begin
        include(joinpath("examples", "simple_state_to_state.jl"))
    end

    println("")

end


text:
export rule

@rule Uninformative(:out, Marginalisation) () = missing


text:
# strongly inspired by https://github.com/s-macke/VoxelSpace
module VoxelSpace
using MiniFB, ColorTypes, PNGFiles, Printf
using Base: unsafe_trunc

include("map_utils.jl")
include("render.jl")
include("input.jl")

const VIEW_DIST = 1024
const WIN_WIDTH = 1200
const WIN_HEIGHT = 900
const MAP_WIDTH = 1024
const MAP_HEIGHT = 1024
const ZSCALE = 1024f0/2
const HORIZON = 1024f0/2

function run(map="C1W")
    # Load map and create framebuffer
    datac, datah = read_map(map)
    fbuffer = Matrix{RGB24}(undef, WIN_WIDTH, WIN_HEIGHT)
    hbuffer = Vector{Int}(undef, WIN_WIDTH)

    # Create window
    window = mfb_open("VoxelSpace.jl", WIN_WIDTH, WIN_HEIGHT)

    # Set up input callbacks
    key_cb_fun_c = @cfunction(key_cb_fun, Cvoid, (Ptr{Cvoid}, mfb_key, mfb_key_mod, Bool))
    mouse_cb_fun_c = @cfunction(mouse_cb_fun, Cvoid, (Ptr{Cvoid}, Int32, Int32))
    mfb_set_keyboard_callback(window, key_cb_fun_c)
    mfb_set_mouse_move_callback(window, mouse_cb_fun_c)

    # Set up FPS counter
    t1 = time_ns()
    @printf("FPS: %3.0f",0.0)

    # Initial position and orientation
    px = 512f0
    py = 512f0
    pz = 80f0
    θ  = 0f0

    while mfb_wait_sync(window)
        # Compute movement direction and apply to pos
        dx,dy,dz,dθ = compute_movement!(θ)
        px += dx
        py += dy
        pz += dz
        θ  += dθ

        # Render viewport
        render!(fbuffer,hbuffer,datac,datah,(px,py,pz),θ)

        # Compute FPS
        t2 = time_ns()
        @printf("\b\b\b%3.0f", 1E9/(t2-t1))
        t1 = t2

        # Update framebuffer
        state = mfb_update(window, fbuffer)
        state == MiniFB.STATE_OK || break
    end

    mfb_close(window)
    println()

    return nothing
end

end # module

max_stars_repo_name: jebej/VoxelSpace

max_issues_repo_name: jebej/VoxelSpace

max_forks_repo_name: jebej/VoxelSpace