text:
using DashBootstrapComponents

badges = html_span([
    dbc_badge("Primary", pill = true, color = "primary", className = "me-1"),
    dbc_badge("Secondary", pill = true, color = "secondary", className = "me-1"),
    dbc_badge("Success", pill = true, color = "success", className = "me-1"),
    dbc_badge("Warning", pill = true, color = "warning", className = "me-1"),
    dbc_badge("Danger", pill = true, color = "danger", className = "me-1"),
    dbc_badge("Info", pill = true, color = "info", className = "me-1"),
    dbc_badge(
        "Light",
        pill = true,
        color = "light",
        text_color = "dark",
        className = "me-1",
    ),
    dbc_badge("Dark", pill = true, color = "dark"),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module Burghers
export flux,minf,maxf
function flux(x::Float64)
    0.5*x^2
end
@inline function minf(a::Float64,b::Float64)
    if b<=0.0
        return 0.5*b^2
        #return flux(b)
    elseif a>=0.0
        return 0.5*a^2
        #return flux(a)
    else
        return 0.0
    end
end
@inline function maxf(a::Float64,b::Float64)
    #max(flux(a),flux(b))
    0.5*max(a^2,b^2)
end
end

ext: jl

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
# Compact unwind data structure

@struct immutable unwind_info_section_header
    version::Uint32
    commonEncodingsArraySectionOffset::Uint32
    commonEncodingsArrayCount::Uint32
    personalityArraySectionOffset::Uint32
    personalityArrayCount::Uint32
    indexSectionOffset::Uint32
    indexCount::Uint32
end

function show(io::IO,h::unwind_info_section_header)
    println(io,"Unwind Info Section Header:")

    printentry(io,"version",h.version)
    printentry(io,"CEAoffset","0x",hex(h.commonEncodingsArraySectionOffset,2*sizeof(Uint32)))
    printentry(io,"CEACount",h.commonEncodingsArrayCount)
    printentry(io,"persOffset","0x",hex(h.personalityArraySectionOffset,2*sizeof(Uint32)))
    printentry(io,"persCount",h.personalityArrayCount)
    printentry(io,"indexOffset","0x",hex(h.indexSectionOffset,2*sizeof(Uint32)))
    printentry(io,"indexCount",h.indexCount)
end

@struct immutable unwind_info_section_header_index_entry
    functionOffset::Uint32
    secondLevelPagesSectionOffset::Uint32
    lsdaIndexArraySectionOffset::Uint32
end

function show(io::IO,h::unwind_info_section_header_index_entry)
    println(io,"Unwind Info Section Header:")

    printentry(io,"functionOffset","0x",hex(h.functionOffset,2*sizeof(Uint32)))
    printentry(io,"secondLevelOff","0x",hex(h.secondLevelPagesSectionOffset,2*sizeof(Uint32)))
    printentry(io,"lsdaIndexOffset","0x",hex(h.lsdaIndexArraySectionOffset,2*sizeof(Uint32)))
end

@struct immutable unwind_info_section_header_lsda_index_entry
    functionOffset::Uint32
    lsdaOffset::Uint32
end

@struct immutable unwind_info_regular_second_level_page_header
    entryPageOffset::Uint16
    entryCount::Uint16
end

@struct immutable unwind_info_compressed_second_level_page_header
    entryPageOff::Uint16
    entryCount::Uint16
    encodingsPageOffset::Uint16
    encodingsCount::Uint16
end

function show(io::IO,h::unwind_info_compressed_second_level_page_header)
    println(io,"Compressed Second Level Page:")

    printentry(io,"entryPageOff","0x",hex(h.entryPageOff,2*sizeof(Uint16)))
    printentry(io,"entryCount",h.entryCount)
    printentry(io,"encPageOff","0x",hex(h.encodingsPageOffset,2*sizeof(Uint16)))
    printentry(io,"encCount",h.encodingsCount)
end

function readSecondLevelHeader(io::IO)
    kind = read(io,Uint32)
    if kind == UNWIND_SECOND_LEVEL_REGULAR
        return unpack(io,unwind_info_regular_second_level_page_header)
    elseif kind == UNWIND_SECOND_LEVEL_COMPRESSED
        return unpack(io,unwind_info_compressed_second_level_page_header)
    else
        error("Unrecognized second level page header")
    end
end

function findEntry(io::IO, off)
    header = unpack(io,unwind_info_section_header)
    if header.indexCount == 0
        error("Unwind section is empty")
    end
    seek(io,header.indexSectionOffset)
    idxs = [unpack(io,unwind_info_section_header_index_entry) for _=1:header.indexCount]
    local idx
    for i = 1:header.indexCount
        if i == header.indexCount || idxs[i].functionOffset <= off <= idxs[i+1].functionOffset
            idx = i
            break
        end
    end
    if idxs[idx].functionOffset > off
        error("How could this happen?")
    end
    seek(io,idxs[idx].secondLevelPagesSectionOffset)
    slh = readSecondLevelHeader(io)
    if isa(slh, unwind_info_compressed_second_level_page_header)
        seek(io,idxs[idx].secondLevelPagesSectionOffset+slh.entryPageOff)
        entries = read!(io,Array(Uint32,slh.entryCount))
        for i = 1:slh.entryCount
            if i == slh.entryCount || (entries[i] & 0x00FFFFFF) <= off <= (entries[i+1] & 0x00FFFFFF)
                return entries[i]
            end
        end
        error("Entry not found")
    else
        error("Unimplemented")
    end
end


text:
######################################################################
# ukf.jl
# All about that unscented Kalman filter
######################################################################
#include("kalman.jl")

type UKF <: GaussianFilter
	mu::Vector{Float64}			# belief
	Sigma::Matrix{Float64}		# belief
	length::Float64
	lambda::Float64
	n::Int

	w0::Float64
	wi::Float64
	initialized::Bool
	initializer::Initializer

	#function UKF(m::SearchDomain, alpha::Float64, k::Float64)
	function UKF(m::SearchDomain)
		mu = m.length/2 * ones(2)
		S = 1e5 * eye(2)

		n = 2
		lambda = 2.

		# create w
		w0 = lambda / (n + lambda)
		wi = 0.5 / (n + lambda)

		initer = NaiveInitializer(m)
		return new(mu, S, m.length, lambda, n, w0, wi,false,initer)
	end
end

function update!(ukf::UKF, x::Vehicle, o::Float64)

	# check if initialized..
	if !ukf.initialized
		check_initialization(ukf, x, o)
		return nothing
	end

	# for ease of writing
	w0 = ukf.w0
	wi = ukf.wi

	# Prediction step
	# nothing changes, because jammer doesn't go anywhere

	# Update step
	# Compute X_bar
	mu = ukf.mu
	sqrtSigma = sqrtm((ukf.n+ukf.lambda) * ukf.Sigma)
	sqrtSigma1 = sqrtSigma[:,1]
	sqrtSigma2 = sqrtSigma[:,2]
	X0 = mu
	X1 = mu + sqrtSigma1
	X2 = mu + sqrtSigma2
	X3 = mu - sqrtSigma1
	X4 = mu - sqrtSigma2

	# TODO: check if the imaginary parts are too big
	Z0 = true_bearing(x, real(X0))
	Z1 = true_bearing(x, real(X1))
	Z2 = true_bearing(x, real(X2))
	Z3 = true_bearing(x, real(X3))
	Z4 = true_bearing(x, real(X4))

	# Compute z_hat and differences
	# TODO: I'm not sure this is valid for circular domain
	#zhat = w0*Z0 + wi*(Z1 + Z2 + Z3 + Z4)
	zhat = angle_mean((Z0,Z1,Z2,Z3,Z4), (w0,wi,wi,wi,wi))

	Z0_diff = fit_180(Z0 - zhat)
	Z1_diff = fit_180(Z1 - zhat)
	Z2_diff = fit_180(Z2 - zhat)
	Z3_diff = fit_180(Z3 - zhat)
	Z4_diff = fit_180(Z4 - zhat)

	# Compute St
	# TODO: I'm not sure this is valid for circular domain
	Qt = x.sensor.noise_sigma * x.sensor.noise_sigma
	St = w0*Z0_diff^2 + wi*(Z1_diff^2 + Z2_diff^2 + Z3_diff^2 + Z4_diff^2)
	St += Qt

	# Compute Sigma_xz
	# TODO: Not valid for circular domain?
	Sigma_xz = wi * (sqrtSigma1) * Z1_diff
	Sigma_xz += wi * (sqrtSigma2) * Z2_diff
	Sigma_xz -= wi * (sqrtSigma1) * Z3_diff
	Sigma_xz -= wi * (sqrtSigma2) * Z4_diff

	# Compute gain, mean, and standard deviation
	Kt = Sigma_xz * inv(St)

	ukf.mu = mu + Kt*fit_180(o - zhat)
	#ukf.Sigma = Sigma_bar - Kt*St*Kt'
	# below is same as above, but simpler
	ukf.Sigma = ukf.Sigma - Sigma_xz * inv(St) * Sigma_xz'

	#println("ukf.mu = ", round(ukf.mu,2))
	#println("ukf.Sigma = ", round(ukf.Sigma,2))
	return nothing
end


text:
using Plots
using LaTeXStrings

pyplot()

x = -π:0.01:π

p0=plot(x,sin.(x),title="naive plot")

p1=plot(
    x,
    sin.(x),
    xticks=(
        [-pi,-pi/2,0,pi/2,pi],
        ["-\$\\pi\$","-\$\\pi\$/2","0","\$\\pi\$/2","-\$\\pi\$"]
    ),
    title="sine-curve",
    label="\$\\sin\$ function",
    xlabel="\$x\$",
    ylabel="\$y\$",
    color="green",
)

p2=plot(
    x,
    cos.(x),
    xticks=(
        [-pi,-pi/2,0,pi/2,pi],
        [L"-\pi",L"-\frac{\pi}{2}","0",L"\frac{\pi}{2}",L"\pi"],
    ),
    title="cosine-curve",
    label=L"y=\cos(x)",
    xlabel=L"x",
    ylabel=L"y",
    color=:purple,
)

plot(p0, p1, p2, layout=(3,1))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# in this attempt, I shall play a set of games with the policy and move the policy using
# the data
using Game2048
using Flux, CUDA
CUDA.allowscalar(false)

policy = Chain(
    Conv((2,2), 1 => 256, relu),
    Conv((2,2), 256 => 128, relu),
    Conv((2,2), 128 => 64, relu),
    flatten,
    Dense(64, 4),
    softmax
) |> gpu

board = initboard()

makeinput(board)  = reshape(board, 4, 4, 1, 1)

cb = cu(Float32.(board) |> makeinput)

using BenchmarkTools
@benchmark policy(cb)
# BenchmarkTools.Trial:
#   memory estimate:  60.33 KiB
#   allocs estimate:  1421
#   --------------
#   minimum time:     418.799 μs (0.00% GC)
#   median time:      470.701 μs (0.00% GC)
#   mean time:        520.493 μs (2.57% GC)
#   maximum time:     36.679 ms (43.01% GC)
#   --------------
#   samples:          9507
#   evals/sample:     1

policy_cpu = Chain(
    Conv((2,2), 1 => 256, relu),
    Conv((2,2), 256 => 128, relu),
    Conv((2,2), 128 => 64, relu),
    flatten,
    Dense(64, 4),
    softmax
)

board = initboard()
fb = Float32.(board) |> makeinput
@benchmark policy_cpu(fb)
# BenchmarkTools.Trial:
#   memory estimate:  153.36 KiB
#   allocs estimate:  211
#   --------------
#   minimum time:     232.399 μs (0.00% GC)
#   median time:      295.501 μs (0.00% GC)
#   mean time:        310.653 μs (2.14% GC)
#   maximum time:     4.700 ms (91.46% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1


text:
module Logging

using Match

export debug, info, warn, err, critical, log,
       Logger,
       LogLevel, DEBUG, INFO, WARNING, ERROR, CRITICAL

include("enum.jl")

@enum LogLevel DEBUG INFO WARNING ERROR CRITICAL

type Logger
    name::String
    level::LogLevel
    output::IO
end

function configure(logger=_root; args...)
    for (tag, val) in args
        @match tag begin
            :io       => logger.output = val::IO
            :filename => logger.output = open(val, "w")
            :level    => logger.level  = val::LogLevel
            unk       => Base.error("Logging: unknown configure argument \"$unk\"")
        end
    end

    logger
end

Logger(name::String;args...) = configure(Logger(name, WARNING, STDERR); args...)

const _root = Logger("root")

for (fn,lvl,clr) in ((:debug,    DEBUG,    :cyan),
                     (:info,     INFO,     :blue),
                     (:warn,     WARNING,  :magenta),
                     (:err,      ERROR,    :red),
                     (:critical, CRITICAL, :red))

    @eval function $fn(msg::String, logger = _root)
        if $lvl >= logger.level
            Base.print_with_color($(Expr(:quote, clr)), logger.output, string($lvl), ":", logger.name, ":", msg, "\n")
        end
    end
end

end # module


text:
using PyPlot

plt.style.use("seaborn-paper")

PyPlot.rc("font", family="serif")
PyPlot.rc("text", usetex=true)
PyPlot.matplotlib.rcParams["axes.titlesize"] = 10
PyPlot.matplotlib.rcParams["axes.labelsize"] = 10
PyPlot.matplotlib.rcParams["xtick.labelsize"] = 9
PyPlot.matplotlib.rcParams["ytick.labelsize"] = 9

const sw = 3.40457
const dw = 7.05826

x = -3pi:pi/1000:+3pi
y = sin.(x)./x
z = tanh.(x)./x

fig, axs = plt.subplots(2,1, figsize=(sw,dw), sharex=true, sharey=true, constrained_layout=true)
axs[1].plot(x, y, "-C0")
axs[2].plot(x, z, "--C1")
axs[2].set_xlabel("\$ x \$")
plt.savefig("./julia_example.pdf", bbox_inches="tight")
plt.close()

max_stars_repo_name: cmarcotte/pyplot-paper-shortcuts

max_issues_repo_name: cmarcotte/pyplot-paper-shortcuts

max_forks_repo_name: cmarcotte/pyplot-paper-shortcuts


text:
# ----------
# Problem Gernerator for setting Up the problem
# ----------

#Declare nothing function
nothingfunction(args...) = nothing;

"""
function Problem(dev::Device=CPU();
  # Numerical parameters
                nx = 64,
                ny = nx,
                nz = nx,
                Lx = 2π,
                Ly = Lx,
                Lz = Lx,
   # Drag and/or hyper-viscosity for velocity/B-field
                 ν = 0,
                nν = 1,
                 μ = 0,
                 η = 0,
                nμ = 0,
   # Declare if turn on magnetic field/VP method/Dye Module
    	     B_field = false,
         VP_method = false,
        Dye_Module = false
  # Timestepper and equation options
           stepper = "RK4",
             calcF = nothingfunction,
  # Float type and dealiasing
                 T = Float32,
  aliased_fraction = 1/3)

Construct a 3D `Problem` on device `dev`.
Keyword arguments
=================
  - `dev`: (required) `CPU()` or `GPU()`; computer architecture used to time-step `problem`.
  - `nx`: Number of grid points in ``x``-domain.
  - `ny`: Number of grid points in ``y``-domain.
  - `nz`: Number of grid points in ``z``-domain.
  - `Lx`: Extent of the ``x``-domain.
  - `Ly`: Extent of the ``y``-domain.
  - `Lz`: Extent of the ``z``-domain.
  - `ν` : Kinematic viscosity coefficient.
  - `η` : Viscosity coefficient for magnetic field.
  - `nν`: (Hyper)-viscosity order, `nν```≥ 1``, not available right now
  - `B_field` :  Declaration of B-field
  - `VP_method`: Declaration of Volume penalization method
  - `Dye_Module`: Declaration of Dye, Passive tracer of the flow;
  - `stepper`: Time-stepping method.
  - `calcF`: Function that calculates the Fourier transform of the forcing, ``F̂``.
  - `aliased_fraction`: the fraction of high-wavenumbers that are zero-ed out by `dealias!()`.
  - `T`: `Float32` or `Float64`; floating point type used for `problem` data.
  - `usr_vars` : user defined variables, cloud either be emtpy or strcut
  - `usr_params` : user defined parameters, cloud either be emtpy or strcut
  - `usr_func` : user defined functions, cloud either be emtpy or strcut
"""

function Problem(dev::Device;
  # Numerical parameters
                nx = 64,
                ny = nx,
                nz = nx,
                Lx = 2π,
                Ly = Lx,
                Lz = Lx,
                dt = 0.0,
   # Drag and/or hyper-viscosity for velocity/B-field
                 ν = 0.0,
                nν = 0,
                 η = 0.0,
                nη = 0,
   # Declare if turn on magnetic field, VP method, Dye module
    	   B_field = false,
	     VP_method = false,
      Dye_Module = false,
  # Timestepper and equation options
           stepper = "RK4",
             calcF = nothingfunction,
  # Float type and dealiasing
                 T = Float32,
  aliased_fraction = 1/3,
  # User defined params/vars
          usr_vars = [],
        usr_params = [],
          usr_func = [])

  if (VP_method)
    error("VP Method is not implemented in this version!")
  end

  # Declare the grid
  grid = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz; T=T)

  # Declare vars
  vars = SetVars(dev, grid, usr_vars; B = B_field, VP =VP_method);

  # Delare params
  params = SetParams(dev, grid, calcF, usr_params;
             B = B_field, VP = VP_method, ν = ν, η = η, nν = nν);

  # Declare Fiuld Equations that will be iterating
  equation = Equation_with_forcing(dev, grid; B = B_field, VP = VP_method);

  # Return the Problem
  return MHDFLowsProblem(equation, stepper, dt, grid, vars, params, dev;
          BFlag = B_field, VPFlag = VP_method, DyeFlag = Dye_Module, usr_func = usr_func)

end

function Equation_with_forcing(dev, grid::AbstractGrid; B = false, VP= false)
  T = eltype(grid);
  Nₗ = ifelse(B,6,3)
  L = zeros(dev, T, (grid.nkr, grid.nl, grid.nm, Nₗ));

  if (B)
  	calcN! = ifelse(VP,MHDcalcN_VP!,MHDcalcN!);
  else
  	calcN! = ifelse(VP, HDcalcN_VP!, HDcalcN!);
  end

  return FourierFlows.Equation(L,calcN!, grid);
end

function MHDcalcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  MHDSolver.MHDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

function HDcalcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  HDSolver.HDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

function MHDcalcN_VP!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  MHDSolver_VP.MHDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

function HDcalcN_VP!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)

  HDSolver_VP.HDcalcN_advection!(N, sol, t, clock, vars, params, grid)

  addforcing!(N, sol, t, clock, vars, params, grid)

  return nothing
end

function addforcing!(N, sol, t, clock, vars, params, grid)
  params.calcF!(N, sol, t, clock, vars, params, grid)
  return nothing
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
#####
##### Checkpointer tests
#####

function test_model_equality(test_model, true_model)
    CUDA.@allowscalar begin
        @test all(test_model.velocities.u.data     .≈ true_model.velocities.u.data)
        @test all(test_model.velocities.v.data     .≈ true_model.velocities.v.data)
        @test all(test_model.velocities.w.data     .≈ true_model.velocities.w.data)
        @test all(test_model.tracers.T.data        .≈ true_model.tracers.T.data)
        @test all(test_model.tracers.S.data        .≈ true_model.tracers.S.data)
        @test all(test_model.timestepper.Gⁿ.u.data .≈ true_model.timestepper.Gⁿ.u.data)
        @test all(test_model.timestepper.Gⁿ.v.data .≈ true_model.timestepper.Gⁿ.v.data)
        @test all(test_model.timestepper.Gⁿ.w.data .≈ true_model.timestepper.Gⁿ.w.data)
        @test all(test_model.timestepper.Gⁿ.T.data .≈ true_model.timestepper.Gⁿ.T.data)
        @test all(test_model.timestepper.Gⁿ.S.data .≈ true_model.timestepper.Gⁿ.S.data)
        @test all(test_model.timestepper.G⁻.u.data .≈ true_model.timestepper.G⁻.u.data)
        @test all(test_model.timestepper.G⁻.v.data .≈ true_model.timestepper.G⁻.v.data)
        @test all(test_model.timestepper.G⁻.w.data .≈ true_model.timestepper.G⁻.w.data)
        @test all(test_model.timestepper.G⁻.T.data .≈ true_model.timestepper.G⁻.T.data)
        @test all(test_model.timestepper.G⁻.S.data .≈ true_model.timestepper.G⁻.S.data)
    end
    return nothing
end

"""
Run two coarse rising thermal bubble simulations and make sure

1. When restarting from a checkpoint, the restarted moded matches the non-restarted
   model to machine precision.

2. When using set!(new_model) to a checkpoint, the new model matches the non-restarted
   simulation to machine precision.

3. run!(new_model, pickup) works as expected
"""
function test_thermal_bubble_checkpointer_output(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
    closure = IsotropicDiffusivity(ν=4e-2, κ=4e-2)
    true_model = IncompressibleModel(architecture=arch, grid=grid, closure=closure)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    CUDA.@allowscalar true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    checkpointed_model = deepcopy(true_model)

    true_simulation = Simulation(true_model, Δt=Δt, stop_iteration=9)
    run!(true_simulation) # for 9 iterations

    checkpointed_simulation = Simulation(checkpointed_model, Δt=Δt, stop_iteration=5)
    checkpointer = Checkpointer(checkpointed_model, schedule=IterationInterval(5), force=true)
    push!(checkpointed_simulation.output_writers, checkpointer)

    # Checkpoint should be saved as "checkpoint5.jld" after the 5th iteration.
    run!(checkpointed_simulation) # for 5 iterations

    # model_kwargs = Dict{Symbol, Any}(:boundary_conditions => SolutionBoundaryConditions(grid))
    restored_model = restore_from_checkpoint("checkpoint_iteration5.jld2")

    #restored_simulation = Simulation(restored_model, Δt=Δt, stop_iteration=9)
    #run!(restored_simulation)

    for n in 1:4
        update_state!(restored_model)
        time_step!(restored_model, Δt, euler=false) # time-step for 4 iterations
    end

    test_model_equality(restored_model, true_model)

    #####
    ##### Test `set!(model, checkpoint_file)`
    #####

    new_model = IncompressibleModel(architecture=arch, grid=grid, closure=closure)

    set!(new_model, "checkpoint_iteration5.jld2")

    @test new_model.clock.iteration == checkpointed_model.clock.iteration
    @test new_model.clock.time == checkpointed_model.clock.time
    test_model_equality(new_model, checkpointed_model)

    #####
    ##### Test `run!(sim, pickup=true)
    #####

    new_simulation = Simulation(new_model, Δt=Δt, stop_iteration=9)

    # Pickup from explicit checkpoint path
    run!(new_simulation, pickup="checkpoint_iteration0.jld2")
    test_model_equality(new_model, true_model)

    run!(new_simulation, pickup="checkpoint_iteration5.jld2")
    test_model_equality(new_model, true_model)

    # Pickup using existing checkpointer
    new_simulation.output_writers[:checkpointer] =
        Checkpointer(new_model, schedule=IterationInterval(5), force=true)

    run!(new_simulation, pickup=true)
    test_model_equality(new_model, true_model)

    run!(new_simulation, pickup=0)
    test_model_equality(new_model, true_model)

    run!(new_simulation, pickup=5)
    test_model_equality(new_model, true_model)

    rm("checkpoint_iteration0.jld2", force=true)
    rm("checkpoint_iteration5.jld2", force=true)

    return nothing
end

function test_checkpoint_output_with_function_bcs(arch)
    grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1))

    @inline some_flux(x, y, t) = 2x + exp(y)
    top_u_bc = top_T_bc = FluxBoundaryCondition(some_flux)
    u_bcs = UVelocityBoundaryConditions(grid, top=top_u_bc)
    T_bcs = TracerBoundaryConditions(grid, top=top_T_bc)

    model = IncompressibleModel(architecture=arch, grid=grid, boundary_conditions=(u=u_bcs, T=T_bcs))
    set!(model, u=π/2, v=ℯ, T=Base.MathConstants.γ, S=Base.MathConstants.φ)

    checkpointer = Checkpointer(model, schedule=IterationInterval(1))
    write_output!(checkpointer, model)
    model = nothing

    restored_model = restore_from_checkpoint("checkpoint_iteration0.jld2")
    @test  ismissing(restored_model.velocities.u.boundary_conditions)
    @test !ismissing(restored_model.velocities.v.boundary_conditions)
    @test !ismissing(restored_model.velocities.w.boundary_conditions)
    @test  ismissing(restored_model.tracers.T.boundary_conditions)
    @test !ismissing(restored_model.tracers.S.boundary_conditions)

    CUDA.@allowscalar begin
        @test all(interior(restored_model.velocities.u) .≈ π/2)
        @test all(interior(restored_model.velocities.v) .≈ ℯ)
        @test all(interior(restored_model.velocities.w) .== 0)
        @test all(interior(restored_model.tracers.T) .≈ Base.MathConstants.γ)
        @test all(interior(restored_model.tracers.S) .≈ Base.MathConstants.φ)
    end
    restored_model = nothing

    properly_restored_model = restore_from_checkpoint("checkpoint_iteration0.jld2",
                                                      boundary_conditions=(u=u_bcs, T=T_bcs))

    CUDA.@allowscalar begin
        @test all(interior(properly_restored_model.velocities.u) .≈ π/2)
        @test all(interior(properly_restored_model.velocities.v) .≈ ℯ)
        @test all(interior(properly_restored_model.velocities.w) .== 0)
        @test all(interior(properly_restored_model.tracers.T) .≈ Base.MathConstants.γ)
        @test all(interior(properly_restored_model.tracers.S) .≈ Base.MathConstants.φ)
    end

    @test !ismissing(properly_restored_model.velocities.u.boundary_conditions)
    @test !ismissing(properly_restored_model.velocities.v.boundary_conditions)
    @test !ismissing(properly_restored_model.velocities.w.boundary_conditions)
    @test !ismissing(properly_restored_model.tracers.T.boundary_conditions)
    @test !ismissing(properly_restored_model.tracers.S.boundary_conditions)

    u, v, w = properly_restored_model.velocities
    T, S = properly_restored_model.tracers

    @test u.boundary_conditions.x.left  isa PBC
    @test u.boundary_conditions.x.right isa PBC
    @test u.boundary_conditions.y.left  isa PBC
    @test u.boundary_conditions.y.right isa PBC
    @test u.boundary_conditions.z.left  isa ZFBC
    @test u.boundary_conditions.z.right isa FBC
    @test u.boundary_conditions.z.right.condition isa ContinuousBoundaryFunction
    @test u.boundary_conditions.z.right.condition.func(1, 2, 3) == some_flux(1, 2, 3)

    @test T.boundary_conditions.x.left  isa PBC
    @test T.boundary_conditions.x.right isa PBC
    @test T.boundary_conditions.y.left  isa PBC
    @test T.boundary_conditions.y.right isa PBC
    @test T.boundary_conditions.z.left  isa ZFBC
    @test T.boundary_conditions.z.right isa FBC
    @test T.boundary_conditions.z.right.condition isa ContinuousBoundaryFunction
    @test T.boundary_conditions.z.right.condition.func(1, 2, 3) == some_flux(1, 2, 3)

    # Test that the restored model can be time stepped
    time_step!(properly_restored_model, 1)
    @test properly_restored_model isa IncompressibleModel

    return nothing
end

function run_cross_architecture_checkpointer_tests(arch1, arch2)
    grid = RegularRectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch1, grid=grid)
    set!(model, u=π/2, v=ℯ, T=Base.MathConstants.γ, S=Base.MathConstants.φ)

    checkpointer = Checkpointer(model, schedule=IterationInterval(1))
    write_output!(checkpointer, model)
    model = nothing

    restored_model = restore_from_checkpoint("checkpoint_iteration0.jld2", architecture=arch2)

    @test restored_model.architecture == arch2

    ArrayType = array_type(restored_model.architecture)
    CUDA.@allowscalar begin
        @test restored_model.velocities.u.data.parent isa ArrayType
        @test restored_model.velocities.v.data.parent isa ArrayType
        @test restored_model.velocities.w.data.parent isa ArrayType
        @test restored_model.tracers.T.data.parent isa ArrayType
        @test restored_model.tracers.S.data.parent isa ArrayType

        @test all(interior(restored_model.velocities.u) .≈ π/2)
        @test all(interior(restored_model.velocities.v) .≈ ℯ)
        @test all(interior(restored_model.velocities.w) .== 0)
        @test all(interior(restored_model.tracers.T) .≈ Base.MathConstants.γ)
        @test all(interior(restored_model.tracers.S) .≈ Base.MathConstants.φ)
    end

    # Test that the restored model can be time stepped
    time_step!(restored_model, 1)
    @test restored_model isa IncompressibleModel

    return nothing
end

function run_checkpointer_cleanup_tests(arch)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, grid=grid)
    simulation = Simulation(model, Δt=0.2, stop_iteration=10)

    simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(3), cleanup=true)
    run!(simulation)

    [@test !isfile("checkpoint_iteration$i.jld2") for i in 1:10 if i != 9]
    @test isfile("checkpoint_iteration9.jld2")

    return nothing
end

for arch in archs
    @testset "Checkpointer [$(typeof(arch))]" begin
        @info "  Testing Checkpointer [$(typeof(arch))]..."

        test_thermal_bubble_checkpointer_output(arch)
        test_checkpoint_output_with_function_bcs(arch)

        if CUDA.has_cuda()
            run_cross_architecture_checkpointer_tests(CPU(), GPU())
            run_cross_architecture_checkpointer_tests(GPU(), CPU())
        end

        run_checkpointer_cleanup_tests(arch)
    end
end


text:
module RRStrace

    using NativeDebugger
    using DebuggerFramework
    using RR

    include("syscalls.jl")
    include("generated/alltheioctls.jl")
    using Cxx

    signal_names = [
        "HUP", "INT", "QUIT", "ILL", "TRAP", "ABRT", "BUS", "FPE",
        "KILL", "USR1", "SEGV", "USR2", "PIPE", "ALRM", "TERM", "STKFLT",
        "CHLD", "CONT", "STOP", "TSTP", "TTIN", "TTOU", "URG", "XCPU", "XFSZ",
        "VTALRM", "PROF", "WINCH", "POLL", "IO", "PWR", "SYS"
    ]

    function apply_property!(task, property, arguments)
        if isa(property, CstringArgument)
            data = NativeDebugger.fallible_load(task, RemotePtr{UInt8}(arguments[property.which]), 40)
            idx = findfirst(data, 0); more = false;
            if idx == 0
                more = true
            else
                data = data[1:idx-1]
            end
            arguments[property.which] = string(repr(String(data)),more ? "..." : "")
        elseif isa(property, CustomTransform)
            arguments[property.which] = property.F(arguments[property.which])
        elseif isa(property, FDArgument)
            arguments[property.which] = sprint((args...)->print_with_color(:yellow,args...), string(arguments[property.which]))
        elseif isa(property, SigmaskPtrArgument) && arguments[property.which] != 0
            mask = NativeDebugger.load(task, RemotePtr{UInt64}(arguments[property.which]))
            sigs = signal_names[filter(x->((1<<(x-1))&mask != 0),1:length(signal_names))]
            arguments[property.which] =
                string("[",join(map(s->string("SIG",s),sigs), " "),"]")
        elseif isa(property, BufferArgument)
            more = arguments[property.whichLength] > 40
            size = min(arguments[property.whichLength], 40)
            data = NativeDebugger.fallible_load(task, RemotePtr{UInt8}(arguments[property.whichData]), size)
            arguments[property.whichData] = string(repr(String(data)),more ? "..." : "")
        end
    end

    function print_syscall_arguments(io, task, ev)
        arch, num = icxx"$ev.arch();", icxx"$ev.Syscall().number;"
        syscall = syscall_table[arch][num]
        if syscall[4] == -1
            print_with_color(:red, io, "(<printing not implemented>)")
        else
            arguments = Any[icxx"$task->regs().arg($i);" for i in 1:syscall[4]]
            if syscall[1] == :ioctl
                category = Char((arguments[2]&0xff00)>>8)
                code = arguments[2]&0xff
                size = (arguments[2]&0xffff0000) >> 16
                idxs = searchsorted(ioctls, (nothing,category,code,nothing), by=x->(x[2],x[3]))
                # Try to find ioctls whose size match exactly. This is useful
                # because some ioctls are only distinguised by size
                idxs′ = filter(x->ioctls[x][4] == size, idxs)
                idxs = isempty(idxs′) ? idxs : idxs′
                if length(idxs) > 1
                    print(io, "({", join(map(x->ioctls[x][1], idxs), " or "), "}, ")
                else
                    print(io, "(", ioctls[idxs[1]][1], ", ")
                end
                ioctl = first(ioctls[idxs])
                print(io,"{",join(map(ioctl[5]) do field
                    name, T, offset = field
                    string(".",name,"=",NativeDebugger.load(task, RemotePtr{T}(arguments[3]+offset)))
                end,", "),"}")
                print(io,")")
                return
            end
            for property in syscall[5]
                apply_property!(task, property, arguments)
            end
            print(io, "(", join(arguments, ", "),")")
        end
    end

    const syscalls_with_hex_printing = [:mmap, :brk]
    function print_syscall_return(io, task, ev)
        arch, num = icxx"$ev.arch();", icxx"$ev.Syscall().number;"
        result = icxx"$task->regs().syscall_result_signed();"
        syscallname = syscall_table[arch][num][1]
        if result < 0
          print(io, " = ")
          print_with_color(:red, io, unsafe_string(icxx"rr::errno_name(-$result);"))
        else
          print(io, " = ", (syscallname in syscalls_with_hex_printing) ?
              string("0x",hex(result, 2sizeof(UInt64))) :
              result)
        end
    end

    const x86_arch = icxx"rr::SupportedArch::x86;"
    const x64_arch = icxx"rr::SupportedArch::x86_64;"

    syscall_table = Dict(
      x86_arch => Dict(x[2] => x for x in filter(x->x[2] != -1, syscalls)),
      x64_arch => Dict(x[3] => x for x in filter(x->x[3] != -1, syscalls))
    )
    function strace(timeline)
        result = icxx"rr::REPLAY_CONTINUE;"
        while icxx"$result == rr::REPLAY_CONTINUE;"
            session = RR.current_session(timeline)
            ev = icxx"rr::Event ev = $session->current_trace_frame().event(); ev;"
            tid = RR.tid(icxx"$session->current_trace_frame();")
            time = RR.time(icxx"$session->current_trace_frame();")
            result = icxx"$session->replay_step(rr::RUN_CONTINUE).status;"
            if icxx"""$ev.is_syscall_event() &&
                      !$session->current_step_key().in_execution() &&
                      $ev.Syscall().state == rr::EXITING_SYSCALL;"""
                print(STDOUT,"[",tid," ",time,"] ")
                arch, num = icxx"$ev.arch();", icxx"$ev.Syscall().number;"
                task = RR.current_task(session)
                print_with_color(:bold, STDOUT, string(syscall_table[arch][num][1]))
                print_syscall_arguments(STDOUT, task, ev)
                print_syscall_return(STDOUT, task, ev)
                println()
            end
        end
    end

end # module


text:
ccount(a) = (ndims(a) == 1 ? length(a) : size(a, ndims(a)))

indbatch(x, b, offset = 0) = (C = ccount(x); min(i + offset, C):min(i + offset + b -1, C) for i in 1:b:C)

macro pthreads(flag, ex)
    esc(:($flag ? Threads.@threads($ex) : $ex))
end

function copy_batch!(dest, src; desc = "copy_batch: ")
    @assert !isempty(src)
    dmax = ndims(src) # argmax(size(src))
    mem = prod(size(src)) * sizeof(eltype(src))
    batchmem = parse(Int, get(ENV, "HDF5_BATCH_MEM", "1024")) * 1024^2
    nbatch = ceil(Int,  mem / min(batchmem, Sys.free_memory() / 10))
    batchsize = ceil(Int, size(src, dmax) / nbatch)
    slices = collect(indbatch(1:size(src, dmax), batchsize))
    prog = Progress(length(slices), desc = desc)
    flag = isa(dest, AbstractArray) && get(ENV, "HDF5_USE_THREADS", "0") == "1"
    verbose = parse(Int, get(ENV, "HDF5_PROGRESS", "1"))
    @pthreads flag for slice in slices
        is = ntuple(d -> d == dmax ? slice : (:), ndims(src))
        dest[is...] = convert(Array, src[is...])
        verbose == 1 && next!(prog)
    end
    return dest
end

function d_zeros(parent, path, T, dims, a...)
    has(parent, path) && o_delete(path)
    @assert all(x -> x > 0, dims)
    dset = d_create(parent, path, datatype(T), dataspace(tuple(dims...)), a...)
    copy_batch!(dset, Zeros{T}(dims...), desc = "zeros.$path ")
    flush(dset)
    return dset
end

write_batch(parent, name, data, a...) = parent[name, a...] = data

function write_batch(parent, name, data::AbstractArray, a...)
    has(parent, name) && o_delete(parent, name)
    T, dims = eltype(data), size(data)
    if get(ENV, "HDF5_USE_THREADS", "0") == "1" &&
        !any(in(("chunk", "compress", "blosc")), a)
        dset = d_zeros(parent, name, T, dims, a...)
        arr = readmmap(dset)
        copy_batch!(arr, data, desc = "write.$name ")
        Mmap.sync!(arr)
    else
        dset = d_create(parent, name, datatype(T), dataspace(dims), a...)
        copy_batch!(dset, data, desc = "write.$name ")
    end
end

write_batch(parent, dict::AbstractDict, a...) = for (k, v) in dict write_batch(parent, k, v, a...) end

read_nonarray(fid, s) = eval(Meta.parse(String(read(fid["nonarray"][s]))))

function write_nonarray(fid, s, x)
    !exists(fid, "nonarray") && g_create(fid, "nonarray")
    has(fid["nonarray"], s) && o_delete(fid["nonarray"], s)
    fid["nonarray"][s] = isa(x, AbstractString) ? "\"" * x * "\"" : string(x)
end


text:
# Unit Test for `SolveQuadraticProgram()`

using LinearAlgebra;
using SparseArrays;
# using Convex
# using SCS
# using ECOS #<! Seems to be very inaccurate for this problem
using Gurobi; #<! Solves problems which SCS struggles with
using OSQP;
using Plots;
using Random;
using StableRNGs;
using MAT;

seedNumber = 1234;
Random.seed!(seedNumber);
# Random.seed!(StableRNG(seedNumber), seedNumber);

include("GenerateQuadraticProgram.jl");
include("SolveQuadraticProgram.jl");
include("LinearSystemSolvers.jl");
include("SolveQuadraticProgramJump.jl");

# @enum ProblemClass randomQp = 1 inequalityConstrainedQp equalityConstrainedQp optimalControl portfolioOptimization lassoOptimization huberFitting supportVectorMachine isotonicRegression
@enum DataSource dataSourceGenerated = 1 dataSourceLoaded

## Parameters

# Simulaion
numSimulations  = 10;
numElements     = 800;
numConstraints  = 0; #<! Set to 0 for OSQP Paper dimensions

dataSource      = dataSourceGenerated;
dataFileName    = "QpModel.mat";

# problemClass    = rand(instances(ProblemClass));
problemClass    = isotonicRegression;

# Solver
numIterations   = 5000;
ρ               = 0.1;
adptΡ           = true;
numItrPolish    = 0;

if (dataSource == dataSourceGenerated)
    mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, numElements; numConstraints = numConstraints);
else
    dMatFile = matread(dataFileName);
    mP = dMatFile["mP"];
    vQ = dropdims(dMatFile["vQ"]; dims = 2);
    mA = dMatFile["mA"];
    vL = dropdims(dMatFile["vL"]; dims = 2);
    vU = dropdims(dMatFile["vU"]; dims = 2);
end

# Some problems might change the actual data
numElements     = size(mP, 1);
numConstraints  = size(mA, 1);

vX = zeros(numElements);

fObjFun(vX) = 0.5 * dot(vX, mP, vX) + dot(vQ, vX);

hLinSolInit = QDLdlInit;
hLinSol     = QDLdl!;

runTimeJump = @elapsed vT = SolveQpJump(mP, vQ, mA, vL, vU; hOptFun = Gurobi.Optimizer);
runTimeOsqp = @elapsed begin
    osqPModel = OSQP.Model();
    OSQP.setup!(osqPModel; P = mP, q = vQ, A = mA, l = vL, u = vU, rho = ρ, eps_abs = 1e-6, eps_rel = 1e-6, scaling = 0);
    osqpRes = OSQP.solve!(osqPModel);
end
runTime     = @elapsed convFlag = SolveQuadraticProgram!(vX, mP, vQ, mA, vL, vU, hLinSolInit, hLinSol; numIterations = numIterations, ρ = ρ, adptΡ = adptΡ, numItrPolish = numItrPolish);

maxAbsDev = norm(vT - vX, Inf);

display(scatter([vX, vT], title = "Solver Solution\nMax Absolute Deviation: $(maxAbsDev)\nRun Time: $(runTime) [Sec], Solver: $(hLinSol)", label = ["Solver" "Reference"]));
display(scatter([mA * vT, vL, vU], title = "Constraints Map", label = ["Solver" "Lower Boundary" "Upper Boundary"]));

println("The max absolute error is: $(maxAbsDev)");
println("The run time is: $(runTime) [Sec]");
println("The output flag is: $(convFlag)");
println("The JuMP.jl run time is: $(runTimeJump) [Sec]");
println("The OSQP run time is: $(runTimeOsqp) [Sec]");
println("The OSQP max absolute error is: $(norm(vT - osqpRes.x, Inf))");

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
#-----------------------------------------------------------------------# AutoCov and Lag
# Lag
"""
    Lag{T}(b::Integer)

Store the last `b` values for a data stream of type `T`.  Values are stored as

``v(t), v(t-1), v(t-2), …, v(t-b+1)``

so that `value(o::Lag)[1]` is the most recent observation and `value(o::Lag)[end]` is the
`b`-th most recent observation.

# Example

    o = fit!(Lag{Int}(10), 1:12)
    o[1]
    o[end]
"""
mutable struct Lag{T} <: OnlineStat{T}
    value::Vector{T}
    b::Int
    n::Int
end
Lag(T::Type, b::Integer) = Lag(T[], b, 0)
function _fit!(o::Lag, y)
    o.n += 1
    pushfirst!(o.value, y)
    length(o.value) > o.b && pop!(o.value)
end
Base.length(o::Lag) = length(o.value)
Base.getindex(o::Lag, i) = o.value[i]
Base.lastindex(o::Lag) = length(o)

"""
    AutoCov(b, T = Float64; weight=EqualWeight())

Calculate the auto-covariance/correlation for lags 0 to `b` for a data stream of type `T`.

# Example

    y = cumsum(randn(100))
    o = AutoCov(5)
    fit!(o, y)
    autocov(o)
    autocor(o)
"""
struct AutoCov{T, W} <: OnlineStat{T}
    cross::Vector{Float64}
    m1::Vector{Float64}
    m2::Vector{Float64}
    lag::Lag{T}         # y_{t-1}, y_{t-2}, ...
    wlag::Lag{Float64}  # γ_{t-1}, γ_{t-2}, ...
    v::Variance{W}
end
function AutoCov(k::Integer, T = Float64; kw...)
    d = k + 1
    AutoCov(zeros(d), zeros(d), zeros(d), Lag(T, d), Lag(Float64, d), Variance(;kw...))
end
nobs(o::AutoCov) = nobs(o.v)

function _fit!(o::AutoCov{T}, y) where {T}
    γ = o.v.weight(o.v.n + 1)
    _fit!(o.v, y)
    _fit!(o.lag, y)     # y_t, y_{t-1}, ...
    _fit!(o.wlag, γ)    # γ_t, γ_{t-1}, ...
    # M1 ✓
    for k in reverse(2:length(o.m2))
        @inbounds o.m1[k] = o.m1[k - 1]
    end
    @inbounds o.m1[1] = smooth(o.m1[1], y, γ)
    # Cross ✓ and M2 ✓
    @inbounds for k in 1:length(o.m1)
        γk = k ≤ length(o.wlag) ? o.wlag[k] : 0.0
        lagk = k ≤ length(o.lag) ? o.lag[k] : zero(T)
        o.cross[k] = smooth(o.cross[k], y * lagk, γk)
        o.m2[k] = smooth(o.m2[k], y, γk)
    end
end

function value(o::AutoCov)
    μ = mean(o.v)
    n = nobs(o)
    cr = o.cross
    m1 = o.m1
    m2 = o.m2
    [(n - k + 1) / n * (cr[k] + μ * (μ - m1[k] - m2[k])) for k in 1:length(m1)]
end
autocov(o::AutoCov) = value(o)
autocor(o::AutoCov) = value(o) ./ value(o)[1]

#-----------------------------------------------------------------------# Bootstrap
"""
    Bootstrap(o::OnlineStat, nreps = 100, d = [0, 2])

Calculate an online statistical bootstrap of `nreps` replicates of `o`.  For each call to `fit!`,
any given replicate will be updated `rand(d)` times (default is double or nothing).

# Example

    o = Bootstrap(Variance())
    fit!(o, randn(1000))
    confint(o, .95)
"""
struct Bootstrap{T, O <: OnlineStat{T}, D} <: OnlineStat{T}
    stat::O
    replicates::Vector{O}
    rnd::D
end
function Bootstrap(o::OnlineStat{T}, nreps::Integer = 100, d = [0, 2]) where {T}
    Bootstrap{T, typeof(o), typeof(d)}(o, [copy(o) for i in 1:nreps], d)
end
Base.show(io::IO, b::Bootstrap) = print(io, "Bootstrap($(length(b.replicates))): $(b.stat)")
"""
    confint(b::Bootstrap, coverageprob = .95)

Return a confidence interval for a Bootstrap `b`.
"""
function confint(b::Bootstrap, coverageprob = 0.95)
    states = value.(b.replicates)
    α = 1 - coverageprob
    return (quantile(states, α / 2), quantile(states, 1 - α / 2))
end
function fit_replicates!(b::Bootstrap, yi)
    for r in b.replicates
        for _ in 1:rand(b.rnd)
            _fit!(r, yi)
        end
    end
end
function _fit!(b::Bootstrap, y)
    _fit!(b.stat, y)
    fit_replicates!(b, y)
    b
end

#-----------------------------------------------------------------------# CallFun
"""
    CallFun(o::OnlineStat, f::Function)

Call `f(o)` every time the OnlineStat `o` gets updated.

# Example

    o = CallFun(Mean(), println)
    fit!(o, [0,0,1,1])
"""
struct CallFun{T, O <: OnlineStat{T}, F <: Function} <: OnlineStat{T}
    stat::O
    f::F
end
CallFun(o::O, f::F) where {T, O<:OnlineStat{T}, F} = CallFun{T, O, F}(o, f)
nobs(o::CallFun) = nobs(o.stat)
value(o::CallFun) = value(o.stat)
Base.show(io::IO, o::CallFun) = print(io, "CallFun: $(o.stat) |> $(o.f)")
_fit!(o::CallFun, arg)  = (_fit!(o.stat, arg); o.f(o.stat))
_merge!(o::CallFun, o2::CallFun) = _merge!(o.stat, o2.stat)

#-----------------------------------------------------------------------# CovMatrix
"""
    CovMatrix(p=0; weight=EqualWeight())
    CovMatrix(::Type{T}, p=0; weight=EqualWeight())

Calculate a covariance/correlation matrix of `p` variables.  If the number of variables is
unknown, leave the default `p=0`.

# Example

    o = fit!(CovMatrix(), randn(100, 4))
    cor(o)
    cov(o)
    mean(o)
    var(o)
"""
mutable struct CovMatrix{T,W} <: OnlineStat{VectorOb} where T<:Number
    value::Matrix{T}
    A::Matrix{T}  # x'x/n
    b::Vector{T}  # 1'x/n
    weight::W
    n::Int
end
function CovMatrix(::Type{T}, p::Int=0; weight = EqualWeight()) where T<:Number
    CovMatrix(zeros(T,p,p), zeros(T,p,p), zeros(T,p), weight, 0)
end
CovMatrix(p::Int=0; weight = EqualWeight()) = CovMatrix(zeros(p,p), zeros(p,p), zeros(p), weight, 0)
function _fit!(o::CovMatrix{T}, x) where {T}
    γ = o.weight(o.n += 1)
    if isempty(o.A)
        p = length(x)
        o.b = zeros(T, p)
        o.A = zeros(T, p, p)
        o.value = zeros(T, p, p)
    end
    smooth!(o.b, x, γ)
    smooth_syr!(o.A, x, γ)
end
nvars(o::CovMatrix) = size(o.A, 1)
function value(o::CovMatrix; corrected::Bool = true)
    o.value[:] = Matrix(Hermitian((o.A - o.b * o.b')))
    corrected && rmul!(o.value, bessel(o))
    o.value
end
function _merge!(o::CovMatrix, o2::CovMatrix)
    γ = o2.n / (o.n += o2.n)
    smooth!(o.A, o2.A, γ)
    smooth!(o.b, o2.b, γ)
    o
end
Statistics.cov(o::CovMatrix; corrected::Bool = true) = value(o; corrected=corrected)
Statistics.mean(o::CovMatrix) = o.b
Statistics.var(o::CovMatrix; kw...) = diag(value(o; kw...))
function Statistics.cor(o::CovMatrix; kw...)
    value(o; kw...)
    v = 1.0 ./ sqrt.(diag(o.value))
    rmul!(o.value, Diagonal(v))
    lmul!(Diagonal(v), o.value)
    o.value
end

#-----------------------------------------------------------------------# Diff
"""
    Diff(T::Type = Float64)

Track the difference and the last value.

# Example

    o = Diff()
    fit!(o, [1.0, 2.0])
    last(o)
    diff(o)
"""
mutable struct Diff{T <: Number} <: OnlineStat{Number}
    diff::T
    lastval::T
    n::Int
end
Diff(T::Type = Float64) = Diff(zero(T), zero(T), 0)
function _fit!(o::Diff{T}, x) where {T<:AbstractFloat}
    v = convert(T, x)
    o.diff = v - last(o)
    o.lastval = v
    o.n += 1
end
function _fit!(o::Diff{T}, x) where {T<:Integer}
    v = round(T, x)
    o.diff = v - last(o)
    o.lastval = v
    o.n += 1
end
Base.last(o::Diff) = o.lastval
Base.diff(o::Diff) = o.diff

#-----------------------------------------------------------------------# StatLag
"""
    StatLag(stat, b)

Track a moving window (previous `b` copies) of `stat`.

# Example

    fit!(StatLag(Mean(), 10), 1:20)
"""
struct StatLag{T, O<:OnlineStat{T}} <: OnlineStat{T}
    lag::Lag{O}
    stat::O
end
function StatLag(stat::O, b::Integer) where {T, O<:OnlineStat{T}}
    StatLag{T,O}(Lag(O,b), stat)
end
nobs(o::StatLag) = nobs(o.stat)
function _fit!(o::StatLag, y)
    _fit!(o.stat, y)
    _fit!(o.lag, copy(o.stat))
end
function Base.show(io::IO, o::StatLag)
    print(io, name(o, false, true))
    OnlineStatsBase.print_stat_tree(io, o.lag.value)
end

Base.@deprecate_binding StatHistory StatLag

#-----------------------------------------------------------------------# HyperLogLog
# https://storage.googleapis.com/pub-tools-public-publication-data/pdf/40671.pdf

"""
    HyperLogLog(T = Number)
    HyperLogLog{P}(T = Number)

Approximate count of distinct elements of a data stream of type `T`, using `2 ^ P`
"registers".  `P` must be an integer between 4 and 16 (default).

Ref: https://storage.googleapis.com/pub-tools-public-publication-data/pdf/40671.pdf

# Example

    o = HyperLogLog()
    fit!(o, rand(1:100, 10^6))

    using Random
    o2 = HyperLogLog(String)
    fit!(o2, [randstring(20) for i in 1:1000])
"""
mutable struct HyperLogLog{p, T} <: OnlineStat{T}
    M::Vector{Int}
    n::Int
    function HyperLogLog{p}(T::Type=Number) where {p}
        4 ≤ p ≤ 16 || throw(ArgumentError("Number of registers must be in 4:16"))
        new{p,T}(zeros(Int, 2^p), 0)
    end
end
HyperLogLog(T::Type=Number) = HyperLogLog{16}(T)

@deprecate HyperLogLog(p::Number, T::Type) HyperLogLog{p}(T::Type)

function Base.show(io::IO, o::HyperLogLog{p,T}) where {p,T}
    print(io, "HyperLogLog{$p, $T}: n=$(nobs(o)) | value=", value(o))
end

function _fit!(o::HyperLogLog, v)
    o.n += 1
    x = hash(v) % UInt32
    i = (x & mask(o)) + UInt32(1)
    w = (x & ~mask(o))
    o.M[i] = max(o.M[i], UInt32(leading_zeros(w) + 1))
end

function value(o::HyperLogLog)
    E = α(o) * _m(o) * _m(o) * inv(sum(x -> inv(2 ^ x), o.M))
    if E ≤ 5 * _m(o) / 2
        V = sum(==(0), o.M)
        return V == 0 ? E : _m(o) * log(_m(o) / V)
    elseif E ≤ 2 ^ 32 / 30
        return E
    else
        return -2 ^ 32 * log(1 - E / 2 ^ 32)
    end
end

function _merge!(o::HyperLogLog, o2::HyperLogLog)
    length(o.M) == length(o2.M) ||
        error("Merge failed. HyperLogLog objects have different number of registers.")
    o.n += o2.n
    for j in eachindex(o.M)
        @inbounds o.M[j] = max(o.M[j], o2.M[j])
    end
    o
end

@generated _m(o::HyperLogLog{p}) where {p} = 2 ^ p

@generated mask(o::HyperLogLog{p}) where {p} = UInt32(2 ^ p - 1)

@generated α(o::HyperLogLog{4}) = 0.673
@generated α(o::HyperLogLog{5}) = 0.697
@generated α(o::HyperLogLog{6}) = 0.709
@generated α(o::HyperLogLog{p}) where {p} = 0.7213 / (1 + 1.079 / 2 ^ p)

#-----------------------------------------------------------------------# KMeans
"Cluster center and the number of observations"
mutable struct Cluster
    value::Vector{Float64}
    n::Int
end
Base.show(io::IO, o::Cluster) = print(io, "Cluster: nobs=$(o.n), value=$(o.value)")
Base.isless(a::Cluster, b::Cluster) = isless(a.n, b.n)

"""
    KMeans(p, k; rate=LearningRate(.6))

Approximate K-Means clustering of `k` clusters and `p` variables.

# Example

    x = [randn() + 5i for i in rand(Bool, 10^6), j in 1:2]

    o = fit!(KMeans(2, 2), x)

    sort!(o; rev=true)  # Order clusters by number of observations
"""
mutable struct KMeans{T, W} <: OnlineStat{VectorOb}
    value::T
    buffer::Vector{Float64}
    rate::W
    n::Int
end
function KMeans(p::Integer, k::Integer; rate=LearningRate())
    KMeans(Tuple(Cluster(zeros(p), 0) for i in 1:k), zeros(k), rate, 0)
end
nobs(o::KMeans) = sum(x -> x.n, o.value)
function Base.show(io::IO, o::KMeans)
    print(io, "KMeans")
    OnlineStatsBase.print_stat_tree(io, o.value)
end
function Base.sort!(o::KMeans; kw...)
    o.value = Tuple(sort!(collect(o.value); kw...))
    o
end
function _fit!(o::KMeans, x)
    o.n += 1
    if o.n ≤ length(o.value)
        cluster = o.value[o.n]
        cluster.value[:] = collect(x)
        cluster.n += 1
    else
        fill!(o.buffer, 0.0)
        for k in eachindex(o.buffer)
            cluster = o.value[k]
            for j in eachindex(x)
                o.buffer[k] = norm(x[j] - cluster.value[j])
            end
        end
        k_star = argmin(o.buffer)
        cluster = o.value[k_star]
        smooth!(cluster.value, x, o.rate(cluster.n += 1))
    end
end

#-----------------------------------------------------------------------# MovingTimeWindow
"""
    MovingTimeWindow{T<:TimeType, S}(window::DatePeriod)
    MovingTimeWindow(window::DatePeriod; valtype=Float64, timetype=Date)

Fit a moving window of data based on time stamps.  Each observation must be a `Tuple`,
`NamedTuple`, or `Pair` where the first item is `<: Dates.TimeType`.  Only observations
with time stamps in the range

``most_recent_datetime - window <= time_stamp <= most_recent_datetime``

are kept track of.

# Example

    using Dates
    dts = Date(2010):Day(1):Date(2011)
    y = rand(length(dts))

    o = MovingTimeWindow(Day(4); timetype=Date, valtype=Float64)
    fit!(o, zip(dts, y))
"""
mutable struct MovingTimeWindow{T<:TimeType, S, D<:DatePeriod} <: OnlineStat{TwoThings{T,S}}
    values::Vector{Pair{T,S}}
    window::D
    n::Int
end
function MovingTimeWindow{T,S}(window::DatePeriod) where {T<:TimeType, S}
    MovingTimeWindow(Pair{T,S}[], window, 0)
end
function MovingTimeWindow(window::DatePeriod; valtype=Float64, timetype=Date)
    MovingTimeWindow{timetype, valtype}(window)
end
value(o::MovingTimeWindow) = sort!(o.values)
function _fit!(o::MovingTimeWindow, y)
    o.n += 1
    push!(o.values, Pair(y...))
    cutoff = maximum(first, o.values) - o.window
    filter!(x -> x[1] >= cutoff, o.values)
end
function _merge!(a::MovingTimeWindow, b::MovingTimeWindow)
    a.n += (nobs(b) - length(value(b)))
    for y in value(b)
        fit!(a, y)
    end
end

#-----------------------------------------------------------------------# MovingWindow
"""
    MovingWindow(b, T)
    MovingWindow(T, b)

Track a moving window of `b` items of type `T`.

# Example

    o = MovingWindow(10, Int)
    fit!(o, 1:14)
"""
mutable struct MovingWindow{T} <: OnlineStat{T}
    value::Vector{T}
    b::Int
    first::Int
    n::Int
end
MovingWindow(b::Int, T::Type) = MovingWindow(T[], b, 1, 0)
MovingWindow(T::Type, b::Int) = MovingWindow(b, T)
function value(o::MovingWindow)
    perm = vcat(collect(o.first:o.b), collect(1:(o.first-1)))
    o.first = 1
    permute!(o.value, perm)
end
function _fit!(o::MovingWindow, y)
    o.n += 1
    if length(o.value) < o.b
        push!(o.value, y)
    else
        o.value[o.first] = y
        o.first = (o.first == o.b) ? 1 : o.first + 1
    end
end
function Base.getindex(o::MovingWindow, i::Int)
    i2 = i + o.first - 1
    i2 = i2 > o.b ? i2 - o.b : i2
    o.value[i2]
end
Base.lastindex(o::MovingWindow) = o.b

#-----------------------------------------------------------------------# OrderStats
"""
    OrderStats(b::Int, T::Type = Float64; weight=EqualWeight())

Average order statistics with batches of size `b`.

# Example

    o = fit!(OrderStats(100), randn(10^5))
    quantile(o, [.25, .5, .75])
"""
mutable struct OrderStats{T, W} <: OnlineStat{Number}
    value::Vector{T}
    buffer::Vector{T}
    weight::W
    n::Int
end
function OrderStats(p::Integer, T::Type = Float64; weight=EqualWeight())
    OrderStats(zeros(T, p), zeros(T, p), weight, 0)
end
function _fit!(o::OrderStats, y)
    p = length(o.value)
    buffer = o.buffer
    i = (o.n % p) + 1
    o.n += 1
    buffer[i] = y
    if i == p
        sort!(buffer)
        smooth!(o.value, buffer, o.weight(o.n / p))
    end
end
function _merge!(o::OrderStats, o2::OrderStats)
    length(o.value) == length(o2.value) ||
        error("Merge failed.  OrderStats track different batch sizes")
    o.n += o2.n
    smooth!(o.value, o2.value, o2.n / o.n)
end
Statistics.quantile(o::OrderStats, arg...) = quantile(value(o), arg...)

# # tree/nbc help:
# function pdf(o::OrderStats, x)
#     if x ≤ first(o.value)
#         return 0.0
#     elseif x > last(o.value)
#         return 0.0
#     else
#         i = searchsortedfirst(o.value, x)
#         b = nobs(o) / (length(o.value) + 1)
#         return b / (o.value[i] - o.value[i-1])
#     end
# end
# split_candidates(o::OrderStats) = midpoints(value(o))

#-----------------------------------------------------------------------# ProbMap
"""
    ProbMap(T::Type; weight=EqualWeight())
    ProbMap(A::AbstractDict{T, Float64}; weight=EqualWeight())

Track a dictionary that maps unique values to its probability.  Similar to
[`CountMap`](@ref), but uses a weighting mechanism.

# Example

    o = ProbMap(Int)
    fit!(o, rand(1:10, 1000))
    probs(o)
"""
mutable struct ProbMap{T, A<:AbstractDict{T,Float64}, W} <: OnlineStat{T}
    value::A
    weight::W
    n::Int
end
function ProbMap(T::Type; weight = EqualWeight())
    ProbMap{T,OrderedDict{T,Float64}, typeof(weight)}(OrderedDict{T, Float64}(), weight, 0)
end
function ProbMap(d::AbstractDict{T, Float64}; weight=EqualWeight()) where {T}
    ProbMap{T,OrderedDict{T,Float64},typeof(weight)}(d, weight, 0)
end
function _fit!(o::ProbMap, y)
    γ = o.weight(o.n += 1)
    get!(o.value, y, 0.0)   # initialize class probability at 0 if it isn't present
    for ky in keys(o.value)
        if ky == y
            o.value[ky] = smooth(o.value[ky], 1.0, γ)
        else
            o.value[ky] *= (1 - γ)
        end
    end
end
function _merge!(o::ProbMap, o2::ProbMap)
    o.n += o2.n
    merge!((a, b)->smooth(a, b, o2.n / o.n), o.value, o2.value)
    o
end
function probs(o::ProbMap, levels = keys(o.value))
    out = zeros(length(levels))
    for (i, ky) in enumerate(levels)
        out[i] = get(o.value, ky, 0.0)
    end
    sum(out) == 0.0 ? out : out ./ sum(out)
end

#-----------------------------------------------------------------------# P2Quantile
"""
    P2Quantile(τ = 0.5)

Calculate the approximate quantile via the P^2 algorithm.  It is more computationally
expensive than the algorithms used by [`Quantile`](@ref), but also more exact.

Ref: [https://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf](https://www.cse.wustl.edu/~jain/papers/ftp/psqr.pdf)

# Example

    fit!(P2Quantile(.5), rand(10^5))
"""
mutable struct P2Quantile <: OnlineStat{Number}
    q::Vector{Float64}  # marker heights
    n::Vector{Int}      # marker position
    nprime::Vector{Float64}
    τ::Float64
    nobs::Int
end
function P2Quantile(τ::Real = 0.5)
    0 < τ < 1 || error("specified quantile must be in (0, 1)")
    nprime = [1, 1 + 2τ, 1 + 4τ, 3 + 2τ, 5]
    P2Quantile(zeros(5), collect(1:5), nprime, τ, 0)
end
Base.show(io::IO, o::P2Quantile) = print(io, "P2Quantile ($(o.τ)): n=$(nobs(o)) | value=$(value(o))")
value(o::P2Quantile) = o.q[3]
nobs(o::P2Quantile) = o.nobs
# function _merge!(a::P2Quantile, b::P2Quantile)
#     a.τ == b.τ || error("Quantiles are not the same: $(a.τ) != $(b.τ)")
#     a.nobs += b.nobs
#     # q
#     a.q[1] = min(a.q[1], b.q[1])
#     a.q[5] = max(a.q[5], b.q[5])
# end
function _fit!(o::P2Quantile, y::Real)
    o.nobs += 1
    q = o.q
    n = o.n
    nprime = o.nprime
    @inbounds if o.nobs > 5
        # B1
        k = searchsortedfirst(q, y) - 1
        k = min(k, 4)
        k = max(k, 1)
        q[1] = min(q[1], y)
        q[5] = max(q[5], y)
        # B2
        for i in (k+1):5
            n[i] += 1
        end
        nprime[2] += o.τ / 2
        nprime[3] += o.τ
        nprime[4] += (1 + o.τ) / 2
        nprime[5] += 1
        # B3
        for i in 2:4
            _interpolate!(o.q, o.n, nprime[i] - n[i], i)
        end
    # A
    elseif o.nobs < 5
        @inbounds o.q[o.nobs] = y
    else
        @inbounds o.q[o.nobs] = y
        sort!(o.q)
    end
end

function _interpolate!(q, n, di, i)
    @inbounds q1, q2, q3 = q[i-1], q[i], q[i+1]
    @inbounds n1, n2, n3 = n[i-1], n[i], n[i+1]
    @inbounds if (di ≥ 1 && n3 - n2 > 1) || (di ≤ -1 && n1 - n2 < -1)
        d = Int(sign(di))
        df = sign(di)
        v1 = (n2 - n1 + d) * (q3 - q2) / (n3 - n2)
        v2 = (n3 - n2 - d) * (q2 - q1) / (n2 - n1)
        qi = q2 + df / (n3 - n1) * (v1 + v2)
        if q1 < qi < q3
            q[i] = qi
        else
            q[i] += df * (q[i + d] - q2) / (n[i + d] - n2)
        end
        n[i] += d
    end
end

# function parabolic_interpolate(q1, q2, q3, n1, n2, n3, d)
#     qi = q2 + d / (n3 - n1) *
#         ((n2 - n1 + d) * (q3 - q2) / (n3 - n2) + (n3 - n2 - d) * (q2 - q1) / (n2 - n1))
# end

#-----------------------------------------------------------------------# Quantile
"""
    Quantile(q = [.25, .5, .75]; alg=OMAS(), rate=LearningRate(.6))

Calculate quantiles via a stochastic approximation algorithm `OMAS`, `SGD`, `ADAGRAD`, or
`MSPI`.  For better (although slower) approximations, see [`P2Quantile`](@ref) and
[`Hist`](@ref).

# Example

    fit!(Quantile(), randn(10^5))
"""
mutable struct Quantile{T <: Algorithm, W} <: OnlineStat{Number}
    value::Vector{Float64}
    τ::Vector{Float64}
    rate::W
    n::Int
    alg::T
end
function Quantile(τ::AbstractVector = [.25, .5, .75]; alg=OMAS(), rate=LearningRate(.6))
    init!(alg, length(τ))
    Quantile(zeros(length(τ)), sort!(collect(τ)), rate, 0, alg)
end
function _fit!(o::Quantile, y)
    γ = o.rate(o.n += 1)
    len = length(o.value)
    if o.n > len
        qfit!(o, y, γ)
    else
        o.value[o.n] = y
        o.n == len && sort!(o.value)
    end
end
function _merge!(o::Quantile, o2::Quantile)
    o.τ == o2.τ || error("Merge failed. Quantile objects track different quantiles.")
    o.n += o2.n
    γ = nobs(o2) / nobs(o)
    merge!(o.alg, o2.alg, γ)
    smooth!(o.value, o2.value, γ)
end

function qfit!(o::Quantile{SGD}, y, γ)
    for j in eachindex(o.value)
        o.value[j] -= γ * Float64((o.value[j] > y) - o.τ[j])
    end
end
function qfit!(o::Quantile{ADAGRAD}, y, γ)
    for j in eachindex(o.value)
        g = Float64((o.value[j] > y) - o.τ[j])
        o.alg.h[j] = smooth(o.alg.h[j], g * g, 1 / nobs(o))
        o.value[j] -= γ * g / sqrt(o.alg.h[j] + ϵ)
    end
end
function qfit!(o::Quantile{MSPI}, y, γ)
    for i in eachindex(o.τ)
        w = inv(abs(y - o.value[i]) + ϵ)
        halfyw = .5 * y * w
        b = o.τ[i] - .5 + halfyw
        o.value[i] = (o.value[i] + γ * b) / (1 + .5 * γ * w)
    end
end
function qfit!(o::Quantile{OMAS}, y, γ)
    s, t = o.alg.a, o.alg.b
    @inbounds for j in eachindex(o.τ)
        w = inv(abs(y - o.value[j]) + ϵ)
        s[j] = smooth(s[j], w * y, γ)
        t[j] = smooth(t[j], w, γ)
        o.value[j] = (s[j] + (2.0 * o.τ[j] - 1.0)) / t[j]
    end
end

# # OMAP...why is this bad?
# q_init(u::OMAP, p) = u
# function qfit!(o::Quantile{<:OMAP}, y, γ)
#     for j in eachindex(o.τ)
#         w = abs(y - o.value[j]) + ϵ
#         θ = y + w * (2o.τ[j] - 1)
#         o.value[j] = smooth(o.value[j], θ, γ)
#     end
# end

#-----------------------------------------------------------------------# ReservoirSample
"""
    ReservoirSample(k::Int, T::Type = Float64)

Create a sample without replacement of size `k`.  After running through `n` observations,
the probability of an observation being in the sample is `1 / n`.

# Example

    fit!(ReservoirSample(100, Int), 1:1000)
"""
mutable struct ReservoirSample{T<:Number} <: OnlineStat{Number}
    value::Vector{T}
    n::Int
end
function ReservoirSample(k::Int, T::Type = Float64)
    ReservoirSample(zeros(T, k), 0)
end
function _fit!(o::ReservoirSample, y)
    o.n += 1
    if o.n <= length(o.value)
        o.value[o.n] = y
    else
        j = rand(1:o.n)
        if j < length(o.value)
            o.value[j] = y
        end
    end
end
function _merge!(o::T, o2::T) where {T<:ReservoirSample}
    length(o.value) == length(o2.value) || error("Can't merge different-sized samples.")
    p = o.n / (o.n += o2.n)
    for j in eachindex(o.value)
        if rand() > p
            o.value[j] = o2.value[j]
        end
    end
end

# #-----------------------------------------------------------------------# Summarizer
# mutable struct Summarizer{T} <: OnlineStat{T}
#     group::Group
# end
# nobs(o::Summarizer) = nobs(o.group)


text:
export Grid, Coordinate ,Coordinates , dimension , extent , spacing , ncells , range
export directions

struct Grid{N}
    extent::NTuple{N,Float}
    spacing::NTuple{N,Float}
	size::NTuple{N,Int}
end

function  Grid(; extent ,  spacing)
		 length(extent) == length(spacing) ||
         throw(error("spacing and extent must have equal length"))
         N = length(extent)
		 size = round.(Int,extent ./ spacing)
return   Grid(NTuple{N,Float}(extent),NTuple{N,Float}(spacing),NTuple{N,Int}(size))
end

function Base.show(io::IO, g::Grid)
        print(io, "dimension:  $(dimension(g))\n",
			      "size     :  ",size(g), '\n',
                  "extent   :  ",extent(g), '\n',
                  "spacing  :  ",spacing(g))
end

dimension(grid::Grid{ND}) where {ND} = ND

extent(g::Grid) = g.extent
extent(g::Grid,d::Int)  =  g.extent[d]
extent(g::Grid,d::Direction{T}) where T =  g.extent[T]

spacing(g::Grid) = g.spacing
spacing(g::Grid,d::Int)  =  g.spacing[d]
spacing(g::Grid,d::Direction{T}) where T = g.spacing[T]

Base.size(g::Grid) =  g.size
Base.size(g::Grid,d::Int)  =  g.size[d]
Base.size(g::Grid,d::Direction{T}) where T = g.size[T]
ncells(g::Grid)  =  prod(size(g))

Base.range(g::Grid,dir::Direction) =   range(0,extent(g,dir),length = size(g,dir))

Base.CartesianIndices(g::Grid) = CartesianIndices(size(g))

"""
	  Coordinates(g::Grid,gridtype::GridType{T}) where T
Returns the primal or dual coordinates of all points on a grid.

"""
function Coordinates(g::Grid,gridtype::GridType{T}) where T
	Coordinates(g , gridtype , CartesianIndices(g))
end

"""
	 Coordinates(g::Grid,gridtype::GridType{T},CIns::CartesianIndices ) where T
Returns the primal or dual coordinates of a set of cartesian indecies.

"""
function Coordinates(g::Grid,gridtype::GridType{T},CIns::CartesianIndices ) where T
  ( Coordinate(g , gridtype, i )  for i in CIns)
end

"""
	Coordinate(g::Grid ,::GridType{:Dual},Ind::CartesianIndex)

Returns the equivalent dual coordinate of of a `CartesianIndex`.

The function doesnot chack for being inside the grid.
"""
function Coordinate(g::Grid ,::GridType{:Dual},Ind::CartesianIndex)
	δ =  spacing(g)
	(Ind.I .- 1)  .* δ  .+  0.5 .* δ
end
"""
	Coordinate(g::Grid , ::GridType{:Primal}, Ind::CartesianIndex)

Returns the equivalent primal coordinate of of a `CartesianIndex`.

The function doesnot chack for being inside the grid.
"""
function Coordinate(g::Grid , ::GridType{:Primal}, Ind::CartesianIndex)
	δ =  spacing(g)
	(Ind.I .- 1) .* δ
end

_ind(gt::GridType,g::Grid{ND},coor::NTuple{ND,Number}) where ND =  ntuple(i-> _ind(gt,g,coor[i],Direction{i}()),ND)
_ind(::GridType{:Dual},g::Grid,coordinate::Number,dir::Direction) =  1 + round(Int,(coordinate - 0.5 .* spacing(g,dir))  ./ spacing(g,dir))
_ind(::GridType{:Primal},g::Grid,coordinate::Number, dir::Direction) =  1 + round(Int,coordinate / spacing(g,dir))

function Base.CartesianIndex(g::Grid{ND},gt::GridType,coordinate::NTuple{ND,Number}) where ND
 	CartesianIndex(_ind(gt,g,coordinate))
end

function Base.CartesianIndex(g::Grid{ND},gt::GridType,coordinate::Real,dir::Direction) where ND
 	CartesianIndex(_ind(gt,g,coordinate,dir))
end

# FIXME
function Base.CartesianIndices(g::Grid{ND},dir::Direction{D}, r::UnitRange{Int}) where {ND,D}
	selectdim(CartesianIndices(g),D,r)
end
#####

function _convert2ind(gt::GridType,r::Tuple{Direction{D},Real,Real},g::Grid{N}) where {N,D}
	a,b = r[2],r[3] ; ax = r[1]
	ia = _ind(gt,g,a,ax)
	ib = _ind(gt,g,b,ax)
	return(ia:ib)
end
function Base.CartesianIndices(gt::GridType,g::Grid{N},r) where {N}
	 inds = [_convert2ind(gt, r[i],g) for i in 1:N]
	CartesianIndices((inds...,))
end

directions(g::Grid{1}) = (x̂,)
directions(g::Grid{2}) = (x̂,ŷ)
directions(g::Grid{3}) = (x̂,ŷ,ẑ)


text:
"""Minimal working example.
"""

using LibSerialPort

function main(args)
    if length(args) != 2
        println("Usage: $(basename(@__FILE__)) port baudrate")
        println("Available ports:")
        list_ports()
        return
    end
	portname, baudrate = args

	baudrate = parse(Int, baudrate)

	LibSerialPort.open(portname, baudrate) do sp
		sleep(2)

		if bytesavailable(sp) > 0
        	println(String(read(sp)))
    	end

	    write(sp, "hello\n")
	    sleep(0.1)
	    println(readline(sp))
	end

end

main(ARGS)


text:
#=――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
    NOTE
    The article data (see /data/articles.json) and algorithm to parse it are
    taken from AvsAn by Eamon Nerbonne, licensed under the Apache License 2.0.

    For more details about the AvsAn project, see
    <https://github.com/EamonNerbonne/a-vs-an>.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――=#

module Articulate

export indefinite

using JSON

const RAW_DATA = JSON.parsefile(
    joinpath(dirname(@__DIR__), "data", "articles.json"))

@enum Article a an

mutable struct ArticleTable
    article::Article
    prefixes::Dict{Char, ArticleTable}
end

function process(rawdata, default)
    pd = ArticleTable(default, Dict())
    if haskey(rawdata, "article")
        pd.article = rawdata["article"] == "a" ? a : an
    end
    for (c, v) in rawdata
        if c ≠ "article"
            @assert length(c) == 1
            pd.prefixes[c[1]] = process(v, pd.article)
        end
    end
    pd
end

const PROCESSED_DATA = process(RAW_DATA, a)

"""
    indefinite(word)

Determine the correct indefinite article, from “a” or “an”, for the given noun.

```jldoctest
julia> using EnglishText

julia> indefinite("hour")
"an"

julia> indefinite("hand")
"a"
```
"""
function indefinite(word, table=PROCESSED_DATA)
    if isempty(word)
        table.article |> string
    else
        f, r = first(word), word[nextind(word, 1):end]
        if haskey(table.prefixes, f)
            indefinite(r, table.prefixes[f])
        else
            table.article |> string
        end
    end
end

end  # module Articulate

import .Articulate: indefinite
export indefinite


text:
"""
    Bessel3Bridge(t, v, σ)

Bessel(3) bridge from below or above to the point `v` at time `t`,
not crossing `v`, with dispersion σ.
"""
struct Bessel3Bridge <: ContinuousTimeProcess{Float64}
    t::Float64
    v::Float64
    σ::Float64
end

"""
    Bessel{N}(σ)

N-dimensional Bessel process with dispersion σ.
Sample with
```
u = 0.0
t = 0:0.1:1
σ = 1.0
sample(u, t, Bridge.Bessel{3}(σ))
```
"""
struct Bessel{N} <: ContinuousTimeProcess{Float64}
    σ::Float64
end

function sample(u, tt, P::Bessel{N}) where N
    tt2 = tt*P.σ^2
    w = sample(tt2, Wiener(), u).yy.^2
    for i in 2:N
        w .+= (sample(tt2, Wiener()).yy).^2
    end
    SamplePath(tt, sqrt.(w))
end

function bessel3(u, tt, t, v, si)
    @assert u != v
    tt2 = tt*si^2
    w2 = sample(tt2, WienerBridge(t*si^2, 0.), 0.).yy
    w3 = sample(tt2, WienerBridge(t*si^2, 0.), 0.).yy

    if u > v
        w1 = sample(tt2, WienerBridge(t*si^2, 0.), u - v).yy
        SamplePath(tt, v .+ sqrt.(w1.^2 + w2.^2 + w3.^2))
    else
        w1 = sample(tt2, WienerBridge(t*si^2, 0.), v - u).yy
        SamplePath(tt, v .- sqrt.(w1.^2 + w2.^2 + w3.^2))
    end

end

b(t, x, P::Bessel3Bridge) = outer(P.σ)*inv(x - P.v) + inv(P.t - t)*(P.v - x)
σ(t, x, P::Bessel3Bridge) = P.σ
a(t, x, P::Bessel3Bridge) = outer(P.σ)
Γ(t, x, P::Bessel3Bridge) = inv(outer(P.σ))
constdiff(P::Bessel3Bridge) = true

sample(u, tt, P::Bessel3Bridge) = bessel3(u, tt, P.t, P.v, P.σ)

"""
    BesselProp

Bessel type proposal
"""
mutable struct BesselProp <: ContinuousTimeProcess{Float64}
    Target
    t; v
    BesselProp(Target::ContinuousTimeProcess{Float64}, t, v) = new(Target, t, v)
end

r(t, x, P::BesselProp) = inv(x - P.v) + inv((P.t - t)*a(P.t, P.v, P.Target))*(P.v-x)
H(t, x, P::BesselProp) =  1/(x - P.v)^2 + inv((P.t - t)*a(P.t, P.v, P.Target))

b(t, x, P::BesselProp) = b(t, x, P.Target) + a(t, x, P.Target)*r(t, x, P)
σ(t, x, P::BesselProp) = σ(t, x, P.Target)
a(t, x, P::BesselProp) = a(t, x, P.Target)
Γ(t, x, P::BesselProp) = Γ(t, x, P.Target)
constdiff(P::BesselProp) = constdiff(P.Target)

btilde(t, x, P::BesselProp) = 0*x
atilde(t, x, P::BesselProp) = a(P.t,P.v,P.Target)

function lptilde(s,u, P::BesselProp)
        v = P.v
        t = P.t
        at = a(t,v,P.Target)
        #abs(v-u)/sqrt(2pi*at*(t-s)^3)*exp(-(v-u)^2/(2*at*(t-s)))

        1/2*((u - v)^2/(at * (s - t)) - log(2π*at) - 3log(t - s) + 2log(abs(v - u)))
end

"""
    aeuler(u, s:dtmax:t, P, tau=0.5)

Adaptive Euler-Maruyama scheme from https://arxiv.org/pdf/math/0601029.pdf
sampling a path from u at s to t with adaptive stepsize of 2.0^(-k)*dtmax
"""
function aeuler(u, r, P::ContinuousTimeProcess{Float64}, tau = 0.5, kmax = 10)

    s = first(r)
    t = last(r)
    dtmax = step(r)
    k = 1
    yy = zeros(0)
    tt = zeros(0)

    y = u
    dt = 2.0^(-k)*dtmax
    while s < t
        append!(yy, y)
        append!(tt, s)
        k = max(0,k - 1)
        B = b(s, y, P)
        while k < kmax && abs(B-b(s, y + B* 2.0^(-k)*dtmax, P)) > tau
            k = k + 1
        end
        dt = 2.0^(-k)*dtmax
        y = y + B*dt + σ(s, y, P)*sqrt(dt)*randn()
        s = s + dt
    end
    SamplePath(tt,yy)
end

llikelihood(Xcirc::SamplePath{T}, Po::BesselProp) where {T} = llikelihoodleft(Xcirc, Po)


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

@doc raw"""VolumeError captures an error encountered during a volume operation.

    IoK8sApiStorageV1beta1VolumeError(;
        message=nothing,
        time=nothing,
    )

    - message::String : String detailing the error encountered during Attach or Detach operation. This string may be logged, so it should not contain sensitive information.
    - time::IoK8sApimachineryPkgApisMetaV1Time : Time the error was encountered.
"""
mutable struct IoK8sApiStorageV1beta1VolumeError <: SwaggerModel
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    time::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: time

    function IoK8sApiStorageV1beta1VolumeError(;message=nothing, time=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1VolumeError, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiStorageV1beta1VolumeError, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type IoK8sApiStorageV1beta1VolumeError

const _property_map_IoK8sApiStorageV1beta1VolumeError = Dict{Symbol,Symbol}(Symbol("message")=>Symbol("message"), Symbol("time")=>Symbol("time"))
const _property_types_IoK8sApiStorageV1beta1VolumeError = Dict{Symbol,String}(Symbol("message")=>"String", Symbol("time")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1VolumeError }) = collect(keys(_property_map_IoK8sApiStorageV1beta1VolumeError))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1VolumeError }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1beta1VolumeError[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1VolumeError }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1VolumeError[property_name]

function check_required(o::IoK8sApiStorageV1beta1VolumeError)
    true
end

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


text:
struct PipelineTessellationStateCreateInfo
    mHandleRef::Ref{vk.VkPipelineTessellationStateCreateInfo}

    function PipelineTessellationStateCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineTessellationStateCreateFlags = vk.VkFlags(0),
        patchControlPoints::UInt32 = UInt32(0)
    )

        this = new(Ref(vk.VkPipelineTessellationStateCreateInfo(
            vk.VK_STRUCTURE_TYPE_PIPELINE_TESSELLATION_STATE_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkPipelineTessellationStateCreateFlags
            patchControlPoints #::UInt32
        )))
    end
end

function handleRef(this::PipelineTessellationStateCreateInfo)::Ref{vk.VkPipelineTessellationStateCreateInfo}
    this.mHandleRef
end


text:
# try out some ideas for sequences
#   starting with the direct sequence
using Polylogarithms
using DataFrames, CSV
using PyPlot
using Printf
include("utilities.jl")
using SpecialFunctions

# chose a "hard value"
s = -6.084227492516975 - 7.500776841841313im
z = 1.4913747102091242 - 3.6642351854073443im
Li = 26.960114195028023 - 13.176237391440013im

μ = log(z)
(-μ)^(s-1) * SpecialFunctions.gamma(1-s)

result1 = polylog(s, z; min_iterations=1000 )
result2 = Polylogarithms.polylog_series_2(s, z; min_iterations=0 )
result2a = Polylogarithms.polylog_series_2(s, z; max_iterations=1 )
result3 = Polylogarithms.polylog_series_2(s, z; min_iterations=100 )
error1 = Li - result1[1]
error2 = Li - result2[1]
error3 = Li - result3[1]
rel_error1 = abs(error1) / abs(Li)
rel_error2 = abs(error2) / abs(Li)
rel_error3 = abs(error3) / abs(Li)

# a = z.^k ./ k.^s
# k_star = Int(ceil( real(s)/real(log(z)) ))

max_n = 100
results = zeros(Complex{Float64}, max_n)
rel_errors = zeros(Float64, max_n)
for k = 1:1:max_n
    results[k] = Polylogarithms.polylog_series_2(s, z; min_iterations=k, max_iterations=k)[1]
    rel_errors[k] = abs( results[k] - Li) / abs(Li)
end

a = diff(results)

figure()
k = 2:1:max_n
semilogy(k, abs.(a), "-")
plot(k, abs.(real.(a)), "-.")
plot(k, abs.(imag.(a)),"--")
# plot(k_star, abs.(a[k_star]), "o")

# let k=0
#     tmp = z
#     converged = false
#     b = zeros(Complex{Float64}, size(a))
#     max_iterations = length(a)
#     min_iterations = k_star
#     accuracy= 1.0e-12
#     while k<max_iterations && ~converged
#         k = k+1
#         b[k+1] = b[k] + tmp
#         tmp *= z * ( k/(k+1.0) )^s
#         # println("   total = $total")
#         if k > min_iterations && abs(tmp)/abs(b[k+1]) < 0.5*accuracy
#             converged = true
#         end
#     end
#     figure()
#     plot(abs.(b[1:k]))
#     println(b)
# end


text:
@inline function vec2ang_ver2(x, y, z)
    norm = sqrt(x^2 + y^2 + z^2)
    theta = acos(z / norm)
    phi = atan(y, x)
    phi = ifelse(phi > π, π-phi, phi)
    return (theta, phi)
end

@inline function quaternion_rotator(omega, t, rotate_axis)
    #= Generate a quaternion that rotates by the angle omega*t around the rotate_axis axis. =#
    rot_ang = omega * t
    Quaternion(@SVector [cos(rot_ang/2.0), rotate_axis[1]*sin(rot_ang/2.0), rotate_axis[2]*sin(rot_ang/2.0), rotate_axis[3]*sin(rot_ang/2.0)])
end

@inline function vector_rotator(vec, rot_ang, rotate_axis)
    #= Generate a quaternion that rotates by the angle omega*t around the rotate_axis axis. =#
    q = Quaternion(@SVector [cos(rot_ang/2.0), rotate_axis[1]*sin(rot_ang/2.0), rotate_axis[2]*sin(rot_ang/2.0), rotate_axis[3]*sin(rot_ang/2.0)])
    vec_q = Quaternion(0.0, vec)
    rot_vec = q * vec_q / q
    return @SVector [rot_vec.q1, rot_vec.q2, rot_vec.q3]
end


text:
@testset "constructors" begin

    sphere(v) = sqrt(sum(v.^2)) < 10
    @test randomSphereFilling(sphere,1,[0,0,0])

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using FinancialMonteCarlo,Nabla;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=0.02;
T=1.0;
d=0.01;
D=90.0;

Nsim=10000;
Nstep=30;
sigma=0.2
mc=MonteCarloConfiguration(Nsim,Nstep);
toll=0.8

rfCurve=ZeroRate(S0,r,d);

FwdData=Forward(T)
EUData=EuropeanOption(T,K)
AMData=AmericanOption(T,K)
BarrierData=BarrierOptionDownOut(T,K,D)
AsianFloatingStrikeData=AsianFloatingStrikeOption(T)
AsianFixedStrikeData=AsianFixedStrikeOption(T,K)
Model=BlackScholesProcess(sigma);

#f__(x) = pricer(BlackScholesProcess(x[1]),ZeroRate(x[2],x[3],d),mc,EUData);
f__(x) = pricer(BlackScholesProcess(x[1]),ZeroRate(x[2],x[3],d),mc,FwdData);
#x=Float64[sigma,S0]
x=Float64[sigma,S0,r]
gradient(x -> f__(x), x)


text:
using JuliaPackageTutorial
using Test

@testset "JuliaPackageTutorial.jl" begin
    #2x+3y
    @test JuliaPackageTutorial.my_f(2,1) == 7
    @test JuliaPackageTutorial.my_f(3,3) == 15
end

@testset "Derivative Test" begin
    @test JuliaPackageTutorial.derivative_of_my_f(2,1) == 2
end

max_stars_repo_name: HazilMohamed/JuliaPackageTutorial

max_issues_repo_name: HazilMohamed/JuliaPackageTutorial

max_forks_repo_name: HazilMohamed/JuliaPackageTutorial


text:
using Anamnesis
using Test, Random, LinearAlgebra

Random.seed!(999)

const NUMBER_TYPES = [Int, UInt, Float64, Complex{Float64}]

const N_CALL_LOOPS = 32

@testset "Scribe" begin
    for dtype ∈ NUMBER_TYPES
        callcount = [0, 0]
        f = (x::AbstractVector{<:Number}) -> (callcount[1] += 1; x⋅x + 1)
        sf = Scribe(f)
        g = (x::Number; y::String) -> (callcount[2] += 1; x + length(y))
        sg = Scribe(g)

        for i ∈ 1:N_CALL_LOOPS
            x = rand(dtype, rand(64:256))
            ξ = rand(dtype)
            y = randstring(rand(0:31))
            @test sf(x) == f(x)
            @test sg(ξ, y=y) == g(ξ, y=y)
            @test sf(x) == f(x)
            @test sg(ξ, y=y) == g(ξ, y=y)
        end
        @test callcount[1] == 3N_CALL_LOOPS
        @test callcount[2] == 3N_CALL_LOOPS
    end
end

@testset "Macros1" begin
    callcount = [0, 0]
    f(x::AbstractVector{<:Number}) = (callcount[1] += 1; x⋅x - 1)
    g(x::Number; y::String) = (callcount[2] += 1; x - length(y))

    for i ∈ 1:N_CALL_LOOPS
        x = rand(Float64, rand(64:256))
        ξ = rand(Float64)
        y = randstring(rand(0:31))
        @mem f(x)
        @test @scribeof(f)(x) == f(x)
        @mem g(ξ, y=y)
        @test @scribeof(g)(ξ, y=y) == g(ξ, y=y)
    end
    @test callcount[1] == 2N_CALL_LOOPS
    @test callcount[2] == 2N_CALL_LOOPS
end

@testset "Macros2" begin
    callcount = [0, 0]
    @anamnesis f1(x::AbstractVector{<:Number}, y=2) = (callcount[1] += 1; x⋅x - y)
    @anamnesis g1(x::Number; y::String) = (callcount[2] += 1; string(x, y))

    for i ∈ 1:N_CALL_LOOPS
        x = rand(Float64, rand(64:256))
        ξ = rand(Float64)
        y = randstring(rand(0:31))
        @test f1(x) == @rawfunc(f1)(x)
        @test f1(x) == @rawfunc(f1)(x)
        @test g1(ξ, y=y) == @rawfunc(g1)(ξ, y=y)
        @test @inferred g1(ξ, y=y) == @rawfunc(g1)(ξ, y=y)
    end
    @test callcount[1] == 3N_CALL_LOOPS
    @test callcount[2] == 3N_CALL_LOOPS
end

@testset "Macros3" begin
    callcount = [0, 0, 0]
    f2(x::AbstractVector{<:Number}) = (callcount[1] += 1; x⋅x - 1)
    g2(x::Number; y::String) = (callcount[2] += 1; string(x, y))
    h2(x) = (callcount[3] += 1; x)
    @anamnesis begin
        x = [1,2]
        f2(x) + h2(3)
        z = [f2(x), h2(3)]
        string(g2(1, y="fire"), g2(2, y="walk"), h2("with"), h2("me"))
        A = [h2("with") g2(1,y="fire")
             g2(2,y="walk") h2("me")]
    end f2 g2
    @test callcount[1] == 1
    @test callcount[2] == 2
    @test callcount[3] == 6
end


text:
# Julia-specific generators
include("generators.jl")

# Julia-specific update implementations
include("update_rules/equality.jl")
include("update_rules/addition.jl")
include("update_rules/multiplication.jl")
include("update_rules/exponential.jl")
include("update_rules/gaussian_mean_variance.jl")
include("update_rules/gaussian_mean_precision.jl")
include("update_rules/gaussian_weighted_mean_precision.jl")
include("update_rules/gamma.jl")
include("update_rules/log_normal.jl")
include("update_rules/wishart.jl")
include("update_rules/bernoulli.jl")
include("update_rules/categorical.jl")
include("update_rules/transition.jl")
include("update_rules/beta.jl")
include("update_rules/dirichlet.jl")
include("update_rules/gaussian_mixture.jl")
include("update_rules/transition_mixture.jl")
include("update_rules/probit.jl")
include("update_rules/logit.jl")
include("update_rules/softmax.jl")
include("update_rules/nonlinear_unscented.jl")
include("update_rules/nonlinear_sampling.jl")
include("update_rules/nonlinear_extended.jl")
include("update_rules/dot_product.jl")
include("update_rules/poisson.jl")
include("update_rules/moment_constraint.jl")
include("update_rules/chance_constraint.jl")
include("update_rules/point_mass_constraint.jl")
include("update_rules/sample_list.jl")