text:
using Unitful
using Unitful.DefaultSymbols
using Unitful: keV, h, c
using FFTW
FFTW.set_num_threads(Sys.CPU_THREADS)
using PyPlot

# a struct to hold the simulation info
struct PropagationInfo
    x
    y
    fx
    fy
end

# utilities
"The rectangle function."
rect(x) = one(x) .* (abs.(x) .< 0.5)

"frequencies for which the fft is calculated"
fftfreq(n) = ifftshift(-cld(n-1, 2) : fld(n-1, 2)) / n

# physics
"Transfer function for the Fresnel free propagation."
H(fx, fy, λ, z) = exp(im * 2π / λ * z) .* exp(-im * π * λ * z * (fx^2 + fy^2))
#H(fx, fy, λ, z) = exp(im * 2π / λ * z) .* exp(im * 2π/λ / z * (fx^2 + fy^2)) / (im * λ * z)

# gratings
singrating(x, y, m, L) = 0.5(1 + m * cos.(2π * x / L))
rectgrating(x, y, m, L; th = 0) = 0.5(1 + m * sign.(cos.(2π * x / L) - th))
rectgrating(pp::PropagationInfo, m, L; th = 0) = rectgrating.(
    pp.x, pp.y, m, L, th = th)
square(x, y, m, L) = m * rect(x / 2L) * rect(y / 2L)
circle(x, y, m, L) = m * rect(sqrt(x^2 + y^2) / 2L)
circle(pp::PropagationInfo, m, L) = circle.(pp.x, pp.y, m, L)
"Turns an amplitude grating into a phase one."
turn2phase(g) =  exp.(im * π / 2 * g)

function fresnelgrid(;δx = 0.1e-6, δy = δx,
        Nx = 2000, Ny = Nx, efficient = true)
    if efficient
        # FFTs are faster on arrays of lengths equal to powers of 2
        Nx = nextprod([2], Nx)
        Ny = nextprod([2], Ny)
    end

    # x are columns, y are rows
    x = δx * (collect(1:Nx) .- Nx/2)'
    y = δy * (collect(1:Ny) .- Ny/2)

    fx = (fftfreq(length(x)) ./ δx)'
    fy = fftfreq(length(y)) ./ δy

    return PropagationInfo(x, y, fx, fy)
end

tonormalspace(U) = ifft(U, (1,2))
tofourierspace(u) = fft(u, (1,2))
flatwave(pp) = ones(ComplexF64, length(pp.y), length(pp.x))

"Propagate a wavefront in the Fourier space"
propagate(U, fx, fy, λ, z) = U .* H.(fx, fy, λ, z)
propagate(U, pp::PropagationInfo, λ, z) = propagate(U, pp.fx, pp.fy, λ, z)

propagatenormal(u, pp::PropagationInfo, λ, z) = tonormalspace(
    propagate(tofourierspace(u), pp, λ, z))

pp = fresnelgrid(δx = 0.1μm, Nx = 4000, Ny = 1)
period = 10μm
g0 = rectgrating(pp, 1, period)
u1 = flatwave(pp) .* g0
λ = h * c / 46keV
u2(z) = propagatenormal(u1, pp, λ, z)

talbot_distance = 2period^2 / λ

z = reshape(collect(range(0mm, stop = 0.2talbot_distance, length = 2000)), 1, 1, :)
talbot_carpet = abs.(u2(z)).^2

plothunit = mm
plotvunit = μm

imshow(talbot_carpet[1, :, :],
    interpolation = "bicubic",
    origin="lower",
    extent=[
        ustrip(plothunit, z[1]), # left
        ustrip(plothunit, z[end]), # right
        ustrip(plotvunit, pp.x[1]), # bottom
        ustrip(plotvunit, pp.x[end]) # top
        ],
    aspect="auto")

xlabel("z position ($plothunit)")
ylabel("x position ($plotvunit)")

ylim(-20, 20)

savefig("$(dirname(dirname(@__FILE__)))/docs/carpets/julia.png")

max_stars_repo_name: rawlik/rosetta-carpet

max_issues_repo_name: rawlik/rosetta-carpet

max_forks_repo_name: rawlik/rosetta-carpet


text:
# Reset
reset_instr(obj::T) where {T <: AWG_SCPI} = write(obj, "*RST")

# Clear
clear_instr(obj::T) where {T <: AWG_SCPI} = write(obj, "*CLS")

# IDN
get_idn(obj::T) where {T <: AWG_SCPI} = query(obj, "*IDN?")

# Query action completed
query_complete(obj::T) where {T <: AWG_SCPI} = query(obj, "*OPC?")

# Wait for action to complete
wait_complete(obj::T) where {T <: AWG_SCPI} = query(obj, "*WAI")

"""
# Beeper
`set_beep(obj::T;st="off") where {T <: AWG_SCPI}`
- st = {"on", "off"}
"""
set_beep(obj::T;st="off") where {T <: AWG_SCPI} = write(obj, "SYST:BEEP:STAT $st")

# Read error
query_error(obj::T) where {T <: AWG_SCPI} = query(obj, "SYST:ERR?")

ext: jl


text:
function tensordict!(result, a, b)
    for (k,v) in a
        for (l,c) in b
            result[tensor(k,l)] = v*c
        end
    end
    return result
end

tensor_result{N,M,T,V}(::StateDict{N,T}, ::StateDict{M,V}) = StateDict{N+M, promote_type(T,V)}()
tensor_result{N,M,T,V}(::OpDict{N,T}, ::OpDict{M,V}) = OpDict{N+M, promote_type(T,V)}()
tensordict(a, b) = tensordict!(tensor_result(a, b), a, b)

function add_merge!(result, other)
    for (k,v) in other
        add_to_dict!(result, k, v)
    end
    return result
end

function sub_merge!(result, other)
    for (k,v) in other
        add_to_dict!(result, k, -v)
    end
    return result
end

merge_result{K,T,V}(::Dict{K,T}, ::Dict{K,V}) = Dict{K, promote_type(T,V)}()
add_merge(a, b) = add_merge!(merge(merge_result(a,b), a), b)
sub_merge(a, b) = sub_merge!(merge(merge_result(a,b), a), b)

function dscale!(result, d, c)
    for (k,v) in d
        result[k] = c*v
    end
    return result
end
dscale!(d, c) = dscale!(d, d, c)

scale_result{K,V,T}(::Dict{K,V}, ::T) = Dict{K, promote_type(T,V)}()
dscale(d, c) = dscale!(scale_result(d,c), d, c)

function add_to_dict!(dict, label, c)
    if c != 0
        dict[label] = get(dict, label, 0) + c
    end
    return dict
end


text:
using Test, WorldOceanAtlasTools
using NCDatasets
using OceanGrids, Unitful

# Alias for short name
WOA = WorldOceanAtlasTools

# For CI, make sure the download does not hang
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# include("test_urls.jl") # Only include locally, as it URLs randomly fail
include("test_functions.jl")
include("test_citations.jl")


text:
using Test
using LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators

using ClimaCore.Geometry

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())
# set up function space
function hvspace_2D(
    xlim = (-π, π),
    zlim = (0, 4π),
    helem = 10,
    velem = 40,
    npoly = 4,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_names = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.IntervalDomain(
        Geometry.XPoint{FT}(xlim[1]),
        Geometry.XPoint{FT}(xlim[2]),
        periodic = true,
    )
    horzmesh = Meshes.IntervalMesh(horzdomain; nelems = helem)
    horztopology = Topologies.IntervalTopology(horzmesh)

    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

# set up rhs!
hv_center_space, hv_face_space = hvspace_2D((-500, 500), (0, 1000))
#hv_center_space, hv_face_space = hvspace_2D((-500,500),(0,30000), 5, 30)

const MSLP = 1e5 # mean sea level pressure
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const γ = 1.4 # heat capacity ratio
const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
const C_v = R_d / (γ - 1) # heat capacity at constant volume
const R_m = R_d # moist R, assumed to be dry

function pressure(ρθ)
    if ρθ >= 0
        return MSLP * (R_d * ρθ / MSLP)^γ
    else
        return NaN
    end
end

Φ(z) = grav * z

# Reference: https://journals.ametsoc.org/view/journals/mwre/140/4/mwr-d-10-05073.1.xml, Section 5a
function init_dry_rising_bubble_2d(x, z)
    x_c = 0.0
    z_c = 350.0
    r_c = 250.0
    θ_b = 300.0
    θ_c = 0.5
    cp_d = C_p
    cv_d = C_v
    p_0 = MSLP
    g = grav

    # auxiliary quantities
    r = sqrt((x - x_c)^2 + (z - z_c)^2)
    θ_p = r < r_c ? 0.5 * θ_c * (1.0 + cospi(r / r_c)) : 0.0 # potential temperature perturbation

    θ = θ_b + θ_p # potential temperature
    π_exn = 1.0 - g * z / cp_d / θ # exner function
    T = π_exn * θ # temperature
    p = p_0 * π_exn^(cp_d / R_d) # pressure
    ρ = p / R_d / T # density
    ρθ = ρ * θ # potential temperature density

    return (ρ = ρ, ρθ = ρθ)
end

# initial conditions
coords = Fields.coordinate_field(hv_center_space)
face_coords = Fields.coordinate_field(hv_face_space)

Yc = map(coord -> init_dry_rising_bubble_2d(coord.x, coord.z), coords)
uₕ = map(_ -> Geometry.Covariant1Vector(0.0), coords)
w = map(_ -> Geometry.Covariant3Vector(0.0), face_coords)
Y = Fields.FieldVector(Yc = Yc, uₕ = uₕ, w = w)

function rhs_invariant!(dY, Y, _, t)

    cρ = Y.Yc.ρ # scalar on centers
    fw = Y.w # Covariant3Vector on faces
    cuₕ = Y.uₕ # Covariant1Vector on centers
    cρθ = Y.Yc.ρθ

    dρ = dY.Yc.ρ
    dw = dY.w
    duₕ = dY.uₕ
    dρθ = dY.Yc.ρθ

    # 0) update w at the bottom
    # fw = -g^31 cuₕ/ g^33

    hdiv = Operators.Divergence()
    hwdiv = Operators.Divergence()
    hgrad = Operators.Gradient()
    hwgrad = Operators.Gradient()
    hcurl = Operators.Curl()
    hwcurl = Operators.Curl()

    dρ .= 0 .* cρ

    ### HYPERVISCOSITY
    # 1) compute hyperviscosity coefficients

    χθ = @. dρθ = hwdiv(hgrad(cρθ / cρ)) # we store χθ in dρθ
    χuₕ = @. duₕ = hwgrad(hdiv(cuₕ))
    Spaces.weighted_dss!(dρθ)
    Spaces.weighted_dss!(duₕ)

    κ₄ = 100.0 # m^4/s
    @. dρθ = -κ₄ * hwdiv(cρ * hgrad(χθ))
    @. duₕ = -κ₄ * hwgrad(hdiv(χuₕ))

    # 1) Mass conservation
    If2c = Operators.InterpolateF2C()
    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    cw = If2c.(fw)
    fuₕ = Ic2f.(cuₕ)
    cuw = Geometry.Covariant13Vector.(cuₕ) .+ Geometry.Covariant13Vector.(cw)

    dw .= fw .* 0

    # 1.a) horizontal divergence
    dρ .-= hdiv.(cρ .* (cuw))

    # 1.b) vertical divergence
    vdivf2c = Operators.DivergenceF2C(
        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
    )
    # we want the total u³ at the boundary to be zero: we can either constrain
    # both to be zero, or allow one to be non-zero and set the other to be its
    # negation

    # explicit part
    dρ .-= vdivf2c.(Ic2f.(cρ .* cuₕ))
    # implicit part
    dρ .-= vdivf2c.(Ic2f.(cρ) .* fw)

    # 2) Momentum equation

    # curl term
    hcurl = Operators.Curl()
    # effectively a homogeneous Dirichlet condition on u₁ at the boundary
    vcurlc2f = Operators.CurlC2F(
        bottom = Operators.SetCurl(Geometry.Contravariant2Vector(0.0)),
        top = Operators.SetCurl(Geometry.Contravariant2Vector(0.0)),
    )
    #cω³ = hcurl.(cu) # zero because we're only in 1D: we can leave this off for the bubble
    fω¹² = hcurl.(fw) # Contravariant2Vector / Contravariant12Vector
    fω¹² .+= vcurlc2f.(cuₕ) # Contravariant2Vector / Contravariant12Vector

    # cross product
    # convert to contravariant
    # these will need to be modified with topography
    fu¹² =
        Geometry.Contravariant1Vector.(Geometry.Covariant13Vector.(Ic2f.(cuₕ))) # Contravariant12Vector in 3D
    fu³ = Geometry.Contravariant3Vector.(Geometry.Covariant13Vector.(fw))
    @. dw -= fω¹² × fu¹² # Covariant3Vector on faces
    @. duₕ -= If2c(fω¹² × fu³)

    # Needed for 3D:
    #cu¹² = Contravariant12Vector.(cu)
    #@. du += cω³ × cu¹²

    cp = @. pressure(cρθ)
    @. duₕ -= hgrad(cp) / cρ
    vgradc2f = Operators.GradientC2F(
        bottom = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
        top = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
    )
    @. dw -= vgradc2f(cp) / Ic2f(cρ)

    cE = @. (norm(cuw)^2) / 2 + Φ(coords.z)
    @. duₕ -= hgrad(cE)
    @. dw -= vgradc2f(cE)

    # 3) potential temperature

    @. dρθ -= hdiv(cuw * Yc.ρθ)
    @. dρθ -= vdivf2c(fw * Ic2f(cρθ))
    @. dρθ -= vdivf2c(Ic2f(cuₕ * cρθ))

    fcc = Operators.FluxCorrectionC2C(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )

    @. dρ += fcc(fw, cρ)
    @. dρθ += fcc(fw, cρθ)
    # dYc.ρuₕ += fcc(w, Yc.ρuₕ)

    Spaces.weighted_dss!(dY.Yc)
    Spaces.weighted_dss!(dY.uₕ)
    Spaces.weighted_dss!(dY.w)

    return dY
end

dYdt = similar(Y);
rhs_invariant!(dYdt, Y, nothing, 0.0);

# run!
using OrdinaryDiffEq
Δt = 0.025
prob = ODEProblem(rhs_invariant!, Y, (0.0, 1.0))
sol_invariant = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 5.0,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "bubble_invariant"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

# post-processing
anim = Plots.@animate for u in sol_invariant.u
    Plots.plot(u.Yc.ρθ ./ u.Yc.ρ, clim = (300.0, 300.8))
end
Plots.mp4(anim, joinpath(path, "theta.mp4"), fps = 20)

anim = Plots.@animate for u in sol_invariant.u
    Plots.plot(u.Yc.ρ)
end
Plots.mp4(anim, joinpath(path, "rho.mp4"), fps = 20)

If2c = Operators.InterpolateF2C()
anim = Plots.@animate for u in sol_invariant.u
    Plots.plot(Geometry.WVector.(Geometry.Covariant13Vector.(If2c.(u.w))))
end
Plots.mp4(anim, joinpath(path, "vel_w.mp4"), fps = 20)

anim = Plots.@animate for u in sol_invariant.u
    Plots.plot(Geometry.UVector.(Geometry.Covariant13Vector.(u.uₕ)))
end
Plots.mp4(anim, joinpath(path, "vel_u.mp4"), fps = 20)


text:
delays = 0:50
runs = 10
runs2 = 10
T_steps = 50

println("\nTesting MCDTS complete tree, Lorenz Prediction:")
@time begin
@testset "MCDTS single rollout on prediction of Lorenz" begin

    max_depth = 15
    x1 = data[1:end-T_steps,1]
    x2 = data[end-T_steps+1:end,1]
    y1 = data[1:end-T_steps,2]
    y2 = data[end-T_steps+1:end,2]

    # Prediction range-function, zeroth predictor first comp-MSE
    delays = 0:5
    runs = 10

    Random.seed!(1234)
    Tw_in = 1 #prediction horizon insample
    Tw_out = 5 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType, 0, 1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test length(τ_mcdts) == 2
    @test τ_mcdts[2] == 2
    @test 0.045 < L_mcdts < 0.046

    error_wheight_insample = 0.5
    error_wheight_oosample = 0.5
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType, 0, 1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())
    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test length(τ_mcdts) == 4
    @test τ_mcdts == [0, 1, 3, 5]
    @test 0.043 < L_mcdts < 0.045

    # Prediction range-function, linear predictor first comp-MSE
    delays = 0:5
    runs = 5
    Random.seed!(1234)
    Tw_in = 2 #prediction horizon insample
    Tw_out = 2 # prediction horizon out-of-sample
    KNN = 1 # nearest neighbors for pred method
    error_wheight_insample = 0.5
    error_wheight_oosample = 0.5

    PredMeth = MCDTS.local_model("linear", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType, 0, 1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts2 = best_node.τs
    L_mcdts2 = L(best_node)
    @test τ_mcdts2 == [0,5,1,2,3,4]
    @test 0.0002 < L_mcdts2 < 0.0003

    # Prediction range-function, zeroth predictor mean-KL
    delays = 0:10
    runs = 10

    Random.seed!(1234)
    KNN = 3 # nearest neighbors for pred method

    PredMeth = MCDTS.local_model("zeroth")
    PredLoss = MCDTS.PredictionLoss(4)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType), MCDTS.Range_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts3 = best_node.τs
    L_mcdts3 = L(best_node)
    @test τ_mcdts3 == [0,4]
    @test 0.069 < L_mcdts3 < 0.07

    # Prediction range-function, linear predictor first-comp-KL
    delays = 0:5
    runs = 5
    Random.seed!(1234)
    Tw_in = 2 #prediction horizon insample
    Tw_out = 2 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method

    PredMeth = MCDTS.local_model("linear", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(3)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType), MCDTS.Range_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts4 = best_node.τs
    L_mcdts4 = L(best_node)
    @test τ_mcdts4 == [0, 1, 3]
    @test 0.034 < L_mcdts4 < 0.035

    # multivariate prediction range-function, zeroth predictor first-comp-MSE
    delays = 0:5
    runs = 5
    data_sample = Dataset(hcat(x1,y1))

    Random.seed!(1234)
    Tw_in = 1 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType, 0, 1, [error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())

    tree = mcdts_embedding(data_sample, optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts5 = best_node.τs
    ts_mcdts5 = best_node.ts
    L_mcdts5 = L(best_node)
    @test τ_mcdts5 == [0, 2]
    @test ts_mcdts5 == [1, 2]
    @test 0.027 < L_mcdts5 < 0.028

    # multivariate prediction range-function, zeroth predictor first-comp-MSE, less fiducials

    Random.seed!(1234)
    Tw_in = 1 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0
    samplesize = 0.5

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,samplesize,[error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())

    tree = mcdts_embedding(data_sample, optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts5 = best_node.τs
    ts_mcdts5 = best_node.ts
    L_mcdts5 = L(best_node)
    @test τ_mcdts5 == [0, 4, 2]
    @test ts_mcdts5 == [1, 1, 1]
    @test 0.616 < L_mcdts5 < 0.617

    # Prediction Continuity-function, zeroth predictor first comp-MSE
    delays = 0:100
    runs = 10
    Random.seed!(1234)
    Tw_in = 1 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Continuity_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test τ_mcdts == [0, 9]
    @test 0.05 < L_mcdts < 0.06

    # Prediction Continuity-function, zeroth predictor first comp-MSE, more neighbors
    Random.seed!(1234)
    Tw_in = 1 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 8 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Continuity_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node2 = MCDTS.best_embedding(tree)
    τ_mcdts2 = best_node2.τs
    L_mcdts2 = L(best_node2)
    @test τ_mcdts2 == τ_mcdts
    @test 0.065 < L_mcdts2 < 0.066

    # Prediction Range-function, zeroth predictor first comp-MSE, Tw = 60
    delays = 0:5
    runs = 10
    Random.seed!(1234)
    Tw_in = 60 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 3 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(1)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Range_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test τ_mcdts == [0, 2]
    @test 1.026 < L_mcdts < 1.027

    # Prediction Continuity-function, zeroth predictor first all-comp-MSE, Tw = 50, more neighbors
    delays = 0:80
    runs = 10
    Random.seed!(1234)
    KNN = 6 # nearest neighbors for pred method

    PredMeth = MCDTS.local_model("zeroth", KNN)
    PredLoss = MCDTS.PredictionLoss(2)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType), MCDTS.Continuity_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test τ_mcdts == [0, 26, 7, 21]
    @test 0.09 < L_mcdts < 0.1

    # Prediction Continuity-function, linear predictor first all-comp-MSE, Tw = 50, more neighbors
    delays = 0:80
    runs = 10
    Random.seed!(1234)
    Tw_in = 5 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 6 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0

    PredMeth = MCDTS.local_model("linear", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(2)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,1,[error_wheight_insample; error_wheight_oosample]), MCDTS.Continuity_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test τ_mcdts == [0, 9, 5]
    @test 0.033 < L_mcdts < 0.034

    # Prediction Continuity-function, linear predictor first all-comp-MSE, Tw = 50, more neighbors, less fiducials
    delays = 0:80
    runs = 10
    Random.seed!(1234)
    Tw_in = 5 #prediction horizon insample
    Tw_out = 1 # prediction horizon out-of-sample
    KNN = 6 # nearest neighbors for pred method
    error_wheight_insample = 1
    error_wheight_oosample = 0
    samplesize = 0.5

    PredMeth = MCDTS.local_model("linear", KNN, Tw_out, Tw_in)
    PredLoss = MCDTS.PredictionLoss(2)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,samplesize,[error_wheight_insample; error_wheight_oosample]), MCDTS.Continuity_function())

    tree = mcdts_embedding(Dataset(x1), optmodel, w1, delays, runs; max_depth = max_depth)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    L_mcdts = L(best_node)
    @test τ_mcdts == [0, 9, 5]
    @test 0.164 < L_mcdts < 0.165

end
end

true


text:
"""
    Loop_subdivsion(mesh) -> refinement

Construct a refinement of `mesh` by Loop subdivision scheme.
Every face is subdived into four small faces and use weights to smooth the surface.

Only defined for 2D meshes.
"""
function Loop_subdivision(mesh::Mesh{U,3}) where U

    D1 = 3

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    connectMat = connectivity(edges,faces)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE

    verts = Array{vertextype(mesh)}(nv)

    # copy over the existing vertices
    for V in 1 : numvertices(mesh)
        # create a mask for new vertex point
        RingverticesIndices = []
        # loop over all the faces
        for (F, Face) in enumerate(cells(faces))
            for iv = 1 : 3
                # if V is a vertex of the face add all the vertices to the mask
                if Face[iv] == V
                    append!(RingverticesIndices, Face)
                end
            end
        end
       # delete duplicate vertices
        RingverticesIndices = unique(RingverticesIndices)
       # calculate valence for vertex V
        valence = length(RingverticesIndices) - 1
        # calculate weights
        if valence != 3
            weight = 3.0 / (8.0 * valence)
        else
            weight = 3.0 / 16.0
        end
        newvert = [0.0, 0.0, 0.0]
        for iv = 1 : length(RingverticesIndices)
             if RingverticesIndices[iv] == V
                     newvert += mesh.vertices[RingverticesIndices[iv]] * (1.0 - valence * weight)
             else
                     newvert += mesh.vertices[RingverticesIndices[iv]] * weight
             end
        end
        verts[V] = newvert
    end

    # add a vertex for each edge
    for (E, Edge) in enumerate(cells(edges))
        # create a mask for new edge point
        EndVerticesIndices = Edge
        FacesIndex = []
        RelatedVerticesIndices = []
        # Loop over all faces
        for (F, Face) in enumerate(cells(faces))
            if connectMat[F,E] != 0
                # Find out two faces in the mask
                append!(FacesIndex, F)
                ## Loop over vertices in the faces find out related vertices
                for iv = 1 : length(Face)
                    if Face[iv] != Edge[1] && Face[iv] != Edge[2]
                            append!(RelatedVerticesIndices, Face[iv])
                    end
                end
            end
        end
        # verts[NV+E] = cartesian(center(chart(edges, Edge)))
        # compute new edge vertices
        if length(RelatedVerticesIndices) != 2 || length(EndVerticesIndices) != 2
            print("The number of vertices in the mask for a edge point is not right!")
        end
        newvert = [0.0, 0.0, 0.0]
        for iv = 1 : length(EndVerticesIndices)
            newvert = newvert + mesh.vertices[EndVerticesIndices[iv]] * 0.375
        end
        for iv = 1: length(RelatedVerticesIndices)
            newvert = newvert + mesh.vertices[RelatedVerticesIndices[iv]] * 0.125
        end
        verts[NV+E] = newvert
    end

    # Build a matrix that given a coarse face gives
    # access to the coarse edges making up its boundary
    D = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add four faces in each coarse face
    nf = 4NF
    fcs = zeros(celltype(mesh), nf)

    for F in 1 : numcells(faces)

        # retrieve the indices of the fine vertices in the centroids
        # of the edges in the boundary of coarse face F
        es = zeros(Int,3)
        for i in nzrange(D, F)
            E = rows[i]
            j = abs(vals[i]) # local index of edge E in face F
            @assert 1 <= j <= 3
            e = NV + E
            es[j] = e
        end

        # build the four triangles
        r1 = mesh.faces[F][1]
        r2 = mesh.faces[F][2]
        r3 = mesh.faces[F][3]

        fcs[4(F-1)+1] = index(r1,es[3],es[2])
        fcs[4(F-1)+2] = index(r2,es[1],es[3])
        fcs[4(F-1)+3] = index(r3,es[2],es[1])
        fcs[4(F-1)+4] = index(es[1],es[2],es[3])
    end

    Mesh(verts, fcs)
end


text:
abstract type Pars end

"""
    struct summarising settings for updates on mcmc algorithm.

Besides initial values for (c,γ), the values that are most likely to require adjustments are `covθprop`, `δpos`, `δmom`

- `ρlowerbound`: lower bound for ρ, this is a a number in [0,1). Wiener innovations are updated according to
    W_new = ρinit * W_old + sqrt(1-ρinit^2) * W_ind
where W_ind is an independnet Wiener process. At each iteration, ρ is sampled uniformly on [ρlowerbound, 1]

- `η`: Stepsize cooling function for the adaptation to ensure diminishing adaptation.
Set at n -> min(0.1, 1.0/√n)

- `adaptskip`: to define guided proposals, we need to plug-in a momentum vector at time T. This is unobserved, and by default initialised as the zero vector. Now every `adaptskip` iterations, `mT` is updated to the value obtained in that iteration. In addition, every `adaptskip` iterations, adaptive tuning of mcmc
parameters is performed.

- `σobs`: If covariance matrix of extrinsic noise is not specified, then it is set to σobs^2*I

- `dt`: mesh width using for discretisation (the mesh always gets transformed by a time-change s  -> s*(2-s))

- `cinit`, : initial value for `c`. Note that we assume Hamltonain kernel x -> c exp(-|x|^2/(2a^2))

- `γinit`: initial value for `gamma` which quantifies the amount of intrinsic noise.

- `aprior`: prior on `a`

- `cprior`: prior on `c`

- `γprior`: prior on `γ`

- `κ`: in case of landmarksmatching, κ quantifies diffusiveness of the prior on the initial momenta

- `δpos`: stepsize for updating intial positions. Note that the value of δpos is adapted during mcmc-iterations.

- `δmom`: stepsize for updating intial momenta. Note that the value of δmom is adapted during mcmc-iterations.

- `skip_saveITER`: skip every skip_saveITER iteration in saving paths
"""
@with_kw struct Pars_ms <: Pars
    model:: Symbol = :ms
    ρlowerbound::Float64 = 0.99
    η::Any = n -> min(0.1, 1.0/√n)      #min(0.2, 10/n)
    adaptskip::Int64 = 20
    σobs:: Float64 = 0.01
    dt:: Float64 = 0.01
    aprior:: UnivariateDistribution = Exponential(1.0)
    γprior:: UnivariateDistribution = Exponential(1.0)
    cinit::  Float64 = 1.0
    γinit::  Float64 = 0.1
    κ:: Float64 = 100.0
    δpos:: Vector{Float64} = [0.001, 0.01]
    δmom:: Vector{Float64} = [0.001, 0.01, 0.1, 0.5]
    δsgd_mom:: Float64 = 0.01
    δa:: Float64 = 0.1
    δγ:: Float64 = 0.0
    skip_saveITER:: Int64 = 10
end

"""
    struct summarising settings for updates on mcmc algorithm.

Besides initial values for (c,γ), the values that are most likely to require adjustments are `covθprop`, `δpos`, `δmom`, `stdev`, `db`

- `ρlowerbound`: lower bound for ρ, this is a a number in [0,1). Wiener innovations are updated according to
        W_new = ρinit * W_old + sqrt(1-ρinit^2) * W_ind
    where W_ind is an independnet Wiener process. At each iteration, ρ is sampled uniformly on [ρlowerbound, 1]

- `η`: Stepsize cooling function for the adaptation to ensure diminishing adaptation.
Set at n -> min(0.1, 1.0/√n)

- `adaptskip`: to define guided proposals, we need to plug-in a momentum vector at time T. This is unobserved, and by default initialised as the zero vector. Now every `adaptskip` iterations, `mT` is updated to the value obtained in that iteration. In addition, every `adaptskip` iterations, adaptive tuning of mcmc
parameters is performed.

- `σobs`: If covariance matrix of extrinsic noise is not specified, then it is set to σobs^2*I

- `dt`: mesh width using for discretisation (the mesh always gets transformed by a time-change s  -> s*(2-s))

- `cinit`, : initial value for `c`. Note that we assume Hamltonain kernel x -> c exp(-|x|^2/(2a^2))

- `γinit`: initial value for `gamma` which quantifies the amount of intrinsic noise.

- `aprior`: prior on `a`

- `cprior`: prior on `c`

- `γprior`: prior on `γ`

- `stdev`: For the AHS-model we take noise fields centred at points that are both horizontally and vertically separaeted by a distance  that is an integer multiple of 2`stdev`

- `db`: assign noisefields within the box [-db[1], db[1]] x [-db[2], db[2]] (and similarly for other dimensions than 2)

- `κ`: in case of landmarksmatching, κ quantifies diffusiveness of the prior on the initial momenta

- `δpos`: stepsize for updating intial positions. Note that the value of δpos is adapted during mcmc-iterations.

- `δmom`: stepsize for updating intial momenta. Note that the value of δmom is adapted during mcmc-iterations.

- `skip_saveITER`: skip every skip_saveITER iteration in saving paths
"""
@with_kw struct Pars_ahs <: Pars
    model:: Symbol = :ahs
    ρlowerbound::Float64 = 0.99
    η::Any = n -> min(0.1, 1.0/√n)
    adaptskip::Int64 = 20
    σobs:: Float64 = 0.01
    dt:: Float64 = 0.01
    aprior:: UnivariateDistribution = Exponential(1.0)
    γprior:: UnivariateDistribution = Exponential(1.0)
    cinit::  Float64 = 1.0
    γinit::  Float64 = 0.1
    stdev:: Float64 = 0.75
    db::Array{Float64,1} = [2.0, 2.0]
    κ:: Float64 = 100.0
    δpos:: Vector{Float64} = [0.001, 0.01]
    δmom:: Vector{Float64} = [0.001, 0.01, 0.1, 0.5]
    δsgd_mom:: Float64 = 0.1
    δa:: Float64 = 0.1
    δγ:: Float64 = 0.0
    skip_saveITER:: Int64 = 10
end


text:
# Algorithm 5.3
# Runge Kutta Fehlberg

function rkf(f, a, b, α, max_step, min_step, tolerance)
	flag = 1
	result = Tuple{Float64, Float64, Float64}[]
	t = a
	k = Array{Float64}(undef, 6)
	h = max_step
	ω = α
	push!(result, (t, ω, h))
	while flag == 1
		k[1] = h*f(t, ω)
		k[2] = h*f(t + (1/4)*h, ω + (1/4)*k[1])
		k[3] = h*f(t + (3/8)*h, ω + (3/32)*k[1] + (9/32)*k[2])
		k[4] = h*f(t + (12/13)*h, ω + (1932/2197)*k[1] - (7200/2197)*k[2] + (7296/2197)*k[3])
		k[5] = h*f(t + h, ω + (439/216)*k[1] - 8*k[2] + (3680/513)*k[3] - (845/4104)*k[4])
		k[6] = h*f(t + (1/2)*h, ω - (8/27)*k[1] + 2*k[2] - (3544/2565)*k[3] + (1859/4104)*k[4] - (11/40)*k[5])
		R = (1/h)*abs((1/360)*k[1] - (128/4275)*k[3] - (2197/75240)*k[4] + (1/50)*k[5] + (2/55)k[6])
		if R <= tolerance
			t = t + h
			ω = ω + (25/216)*k[1] + (1408/25655)*k[3] + (21977/4104)*k[4] - (1/5)*k[5]
			push!(result, (t, ω, h))
		end
		δ = 0.84 * (tolerance/R)^(1/4)
		if δ <= 0.1
			h = 0.1*h
		elseif δ >= 4
			h = 4h
		else
			h = δ*h
		end
		if h > max_step
			h = max_step
		end
		if t >= b
			flag = 0
		elseif (t + h) > b
			h = b - t
		elseif h < min_step
			flag = 0
		end
	end
	return result
end

f(t, y) = y - t^2 + 1
println(rkf(f, 0, 2, 0.5, 0.25, 0, 10^(-5)))


text:
#####
##### `map`
#####

function rrule(::typeof(map), f, xs...)
    y = map(f, xs...)
    ∂xs = ntuple(length(xs)) do i
        Rule() do ȳ
            map(ȳ, xs...) do ȳi, xis...
                _, ∂xis = _checked_rrule(f, xis...)
                extern(∂xis[i](ȳi))
            end
        end
    end
    return y, (DNERule(), ∂xs...)
end

#####
##### `mapreduce`, `mapfoldl`, `mapfoldr`
#####

for mf in (:mapreduce, :mapfoldl, :mapfoldr)
    sig = :(rrule(::typeof($mf), f, op, x::AbstractArray{<:Real}))
    call = :($mf(f, op, x))
    if mf === :mapreduce
        insert!(sig.args, 2, Expr(:parameters, Expr(:kw, :dims, :(:))))
        insert!(call.args, 2, Expr(:parameters, Expr(:kw, :dims, :dims)))
    end
    body = quote
        y = $call
        ∂x = Rule() do ȳ
            broadcast(x, ȳ) do xi, ȳi
                _, ∂xi = _checked_rrule(f, xi)
                extern(∂xi(ȳi))
            end
        end
        return y, (DNERule(), DNERule(), ∂x)
    end
    eval(Expr(:function, sig, body))
end

#####
##### `sum`
#####

frule(::typeof(sum), x) = (sum(x), Rule(sum))

rrule(::typeof(sum), x) = (sum(x), Rule(cast))

function rrule(::typeof(sum), f, x::AbstractArray{<:Real}; dims=:)
    y, (_, _, ∂x) = rrule(mapreduce, f, Base.add_sum, x; dims=dims)
    return y, (DNERule(), ∂x)
end

function rrule(::typeof(sum), x::AbstractArray{<:Real}; dims=:)
    y, (_, ∂x) = rrule(sum, identity, x; dims=dims)
    return y, ∂x
end

function rrule(::typeof(sum), ::typeof(abs2), x::AbstractArray{<:Real}; dims=:)
    y = sum(abs2, x; dims=dims)
    ∂x = Rule(ȳ -> 2ȳ .* x)
    return y, (DNERule(), ∂x)
end

#####
##### `mean`
#####

_denom(x, dims::Colon) = length(x)
_denom(x, dims::Integer) = size(x, dims)
_denom(x, dims) = mapreduce(i->size(x, i), Base.mul_prod, unique(dims), init=1)

# TODO: We have `mean(f, x; dims)` as of 1.3.0-DEV.36

function rrule(::typeof(mean), x::AbstractArray{<:Real}; dims=:)
    _, dx = rrule(sum, x; dims=dims)
    n = _denom(x, dims)
    return mean(x; dims=dims), Rule(ȳ -> dx(ȳ) / n)
end

function rrule(::typeof(mean), f, x::AbstractArray{<:Real})
    _, (_, dx) = rrule(sum, f, x)
    n = _denom(x, :)
    return mean(f, x), (DNERule(), Rule(ȳ -> dx(ȳ) / n))
end


text:
using Base.Filesystem
using Dates

const TIMEOUT = 2*60*60 # seconds

const HANDLE = Ptr{Cvoid}
const LPVOID = Ptr{Cvoid}
const DWORD = UInt32
const BOOL = Cint
const UINT = Cuint
const JOBOBJECTINFOCLASS = Cint

function CloseHandle(handle)
    Base.windowserror(:CloseHandle,
        0 == ccall(:CloseHandle, stdcall, Cint, (HANDLE,), handle))
    nothing
end
CloseHandle(handle::File) = CloseHandle(handle.handle)

struct WindowsJobObject
    handle::HANDLE
end
Base.cconvert(::Type{HANDLE}, job::WindowsJobObject) = job.handle

struct SECURITY_ATTRIBUTES
    nLength::DWORD
    lpSecurityDescriptor::LPVOID
    bInheritHandle::Cint
end

function CreateJobObject(; sec_attrs = C_NULL, name = C_NULL)
    job = ccall((:CreateJobObjectA, "kernel32"), stdcall, HANDLE, (Ptr{SECURITY_ATTRIBUTES}, Cwstring), sec_attrs, name)
    Base.windowserror(:CreateJobObject, job == C_NULL)
    WindowsJobObject(job)
end

function TerminateJobObject(job, code=1)
    Base.windowserror(:TerminateJobObject,
        ccall((:TerminateJobObject, "kernel32"), BOOL,
            (HANDLE, UINT), job, code) == 0)
    nothing
end

const JobObjectBasicProcessIdList = 3

function QueryJobObjectBasicProcessIdList(job)
    nprocesses = 1024

    while true
        data = Vector{UInt8}(undef, 2 * sizeof(DWORD) + nprocesses * sizeof(UInt))
        outsize = Ref{DWORD}(sizeof(data))

        Base.windowserror(:QueryInformationJobObject, 0 ==
            ccall((:QueryInformationJobObject, "kernel32"), stdcall, BOOL,
                (HANDLE, JOBOBJECTINFOCLASS, Ptr{Cvoid}, DWORD, Ref{DWORD}),
                job, JobObjectBasicProcessIdList, data, sizeof(data), outsize))

        header_end = 2sizeof(DWORD)
        NumberOfAssignedProcesses, NumberOfProcessIdsInList = reinterpret(DWORD, data[1:header_end])
        if NumberOfProcessIdsInList < NumberOfAssignedProcesses
            # Give a bit of margin for new processes to be spawned
            nprocesses = NumberOfAssignedProcesses + 100
            continue
        end

        return reinterpret(UInt, data[header_end+1:end])[1:NumberOfProcessIdsInList]
    end
end

function AssignProcessToJobObject(job::WindowsJobObject, process)
    Base.windowserror(:AssignProcessToJobObject,
        0 == ccall((:AssignProcessToJobObject, "kernel32"), stdcall, BOOL, (HANDLE, HANDLE), job, process))
end

const MINIDUMP_TYPE = Cint
const MiniDumpNormal = 0x00000000
const MiniDumpWithDataSegs = 0x00000001
const MiniDumpWithFullMemory = 0x00000002
const MiniDumpWithHandleData = 0x00000004
const MiniDumpFilterMemory = 0x00000008
const MiniDumpScanMemory = 0x00000010
const MiniDumpWithUnloadedModules = 0x00000020
const MiniDumpWithIndirectlyReferencedMemory = 0x00000040
const MiniDumpFilterModulePaths = 0x00000080
const MiniDumpWithProcessThreadData = 0x00000100
const MiniDumpWithPrivateReadWriteMemory = 0x00000200
const MiniDumpWithoutOptionalData = 0x00000400
const MiniDumpWithFullMemoryInfo = 0x00000800
const MiniDumpWithThreadInfo = 0x00001000
const MiniDumpWithCodeSegs = 0x00002000
const MiniDumpWithoutAuxiliaryState = 0x00004000
const MiniDumpWithFullAuxiliaryState = 0x00008000
const MiniDumpWithPrivateWriteCopyMemory = 0x00010000
const MiniDumpIgnoreInaccessibleMemory = 0x00020000
const MiniDumpWithTokenInformation = 0x00040000
const MiniDumpWithModuleHeaders = 0x00080000
const MiniDumpFilterTriage = 0x00100000
const MiniDumpValidTypeFlags = 0x001fffff

const PROCESS_QUERY_INFORMATION = 0x0400
const PROCESS_VM_READ           = 0x0010
const PROCESS_DUP_HANDLE        = 0x0040
const PROCESS_TERMINATE         = 0x0001
const PROCESS_SET_QUOTA         = 0x0100

function OpenProcess(id::Integer, rights = PROCESS_QUERY_INFORMATION | PROCESS_VM_READ | PROCESS_DUP_HANDLE)
    proc = ccall((:OpenProcess, "kernel32"), stdcall, HANDLE, (DWORD, BOOL, DWORD),
        rights,
        false, id)
    Base.windowserror(:OpenProcess, proc == C_NULL)
    proc
end

function CreateMinidump(filename::String, proc_id::UInt; kind = MiniDumpWithProcessThreadData)
    proc_handle = OpenProcess(proc_id)
    file_handle = Filesystem.open(filename, JL_O_CREAT | JL_O_WRONLY, 0666)

    Base.windowserror(:MiniDumpWriteDump, 0 ==
        ccall((:MiniDumpWriteDump, "dbghelp"), stdcall, BOOL,
            (HANDLE, DWORD, HANDLE, MINIDUMP_TYPE,
                Ptr{Cvoid},
                Ptr{Cvoid},
                Ptr{Cvoid}),
            proc_handle, proc_id, file_handle.handle, kind,
            C_NULL, C_NULL, C_NULL))

    CloseHandle(file_handle)
    CloseHandle(proc_handle)
end

if length(ARGS) < 3
    println(stderr, "Usage: autodump.jl [buildnumber] [shortcommit] [command...]")
    exit(2)
end

run_id = popfirst!(ARGS)
shortcommit = popfirst!(ARGS)
proc = run(`$ARGS`, (stdin, stdout, stderr); wait=false)
job = CreateJobObject()
AssignProcessToJobObject(job, OpenProcess(getpid(proc), PROCESS_SET_QUOTA | PROCESS_TERMINATE))

@async begin
    sleep(TIMEOUT)
    try
        if isopen(proc)
            println(stderr, "\n\nProcess timed out. Creating core dump for each running process!")
            datestr = Dates.format(now(), dateformat"yyyy-mm-dd_HH_MM_SS")
            proc_ids = QueryJobObjectBasicProcessIdList(job)
            println(stderr, "Process IDs: $(Int.(proc_ids))")
            mkpath("dumps")
            foreach(proc_ids) do proc_id
                filename = "dumps\\dump-run_$(run_id)-gitsha_$(shortcommit)-pid_$(proc_id)-$(datestr).dmp"
                try
                    print(stderr, "Dumping $(proc_id) to $(filename)...")
                    CreateMinidump(filename, proc_id)
                    println(stderr, "Done!")
                catch e
                    bt = catch_backtrace()
                    println(stderr, "Failed to create minidump for process $(proc_id)")
                    Base.display_error(stderr, e, bt)
                end
            end
        end
    finally
        TerminateJobObject(job)
        exit(1)
    end
end

wait(proc)
exit(proc.exitcode)

max_stars_repo_name: ianshmean/julia-buildbot

max_issues_repo_name: ianshmean/julia-buildbot

max_forks_repo_name: ianshmean/julia-buildbot


text:
# Sample postprocessing script for plotting a tensor from one snapshot.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, LaTeXStrings, Printf

file = "/home/hongyang/Vlasiator/test/equatorial/timevarying/bulk.0001582.vlsv"

axisunit = EARTH
colorscale = Linear
cmap = matplotlib.cm.turbo

vars = ["proton/vg_ptensor_diagonal", "proton/vg_ptensor_offdiagonal"]

#################

meta = load(file)

pArgs1 = Vlasiator.set_args(meta, vars[1], axisunit; normal=:none)
pArgs2 = Vlasiator.set_args(meta, vars[2], axisunit; normal=:none)

x, y = Vlasiator.get_axis(pArgs1)
# [nPa]
pxx = Vlasiator.prep2d(meta, vars[1], :x)' .* 1e9
pyy = Vlasiator.prep2d(meta, vars[1], :y)' .* 1e9
pzz = Vlasiator.prep2d(meta, vars[1], :z)' .* 1e9
pyz = Vlasiator.prep2d(meta, vars[2], :1)' .* 1e9
pxz = Vlasiator.prep2d(meta, vars[2], :2)' .* 1e9
pxy = Vlasiator.prep2d(meta, vars[2], :3)' .* 1e9

P = (pxx, pyy, pzz, pyz, pxz, pxy)
P_str = ("Pxx", "Pyy", "Pzz", "Pyz", "Pxz", "Pxy")

vmin = minimum(minimum.(P))
vmax = maximum(maximum.(P))

cnorm1, cticks1 = Vlasiator.set_colorbar(colorscale, vmin, vmax)

fig, axs = subplots(3,2,
   figsize=(6, 12), sharex=true, sharey=true, constrained_layout=true)

c1 = axs[1].pcolormesh(x, y, pxx; norm=cnorm1, cmap)

for i in eachindex(axs)[2:end]
   axs[i].pcolormesh(x, y, P[i]; norm=cnorm1, cmap)
end

for (i, ax) in enumerate(axs)
   ax.set_aspect("equal")
   ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   ax.set_title(P_str[i])
end

for ax in axs[3,:]
   ax.set_xlabel(L"x [$R_E$]")
end

for ax in axs[:,1]
   ax.set_ylabel(L"y [$R_E$]")
end

# One colorbar for all subplots
cb1 = colorbar(c1; ax=axs, ticks=cticks1, fraction=0.046, pad=0.04)
cb1.ax.set_ylabel("[nPa]")
cb1.outline.set_linewidth(1.0)

str_title = @sprintf "Pressure Tensor at t = %4.1fs" meta.time

fig.suptitle(str_title, fontsize="xx-large")

savefig("ptensor.png", bbox_inches="tight")


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module PercStep

using Random
using ExtractMacro
using ..Interface
using ..Common

using ..DeltaE.ArraySets

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphPercStep

import ..Interface: energy, delta_energy, update_cache!, neighbors

function gen_ξ(N::Integer, P::Integer)
    ξ = BitArray(undef, P, N)
    ξv = Vector{BitVector}(undef, P)
    for a = 1:P
        ξv[a] = bitrand(N)
        ξ[a,:] = ξv[a]
    end
    # the ξv representation is for energy
    # the ξ representation is for delta_energy and update_cache!
    return ξ, ξv
end

mutable struct Stabilities
    p::ArraySet
    m::ArraySet
    Δs::Vector{Int}
    ξsi::BitVector
    last_move::Int
    function Stabilities(P::Integer)
        p = ArraySet(P)
        m = ArraySet(P)
        Δs = zeros(P)
        ξsi = BitArray(undef, P)
        return new(p, m, Δs, ξsi, 0)
    end
end

mutable struct GraphPercStep <: SimpleGraph{Int}
    N::Int
    P::Int
    ξ::BitMatrix
    ξv::Vector{BitVector}
    stab::Stabilities
    tmps::BitVector
    function GraphPercStep(ξ::BitMatrix, ξv::Vector{BitVector})
        P, N = size(ξ)
        #TODO: check
        isodd(N) || throw(ArgumentError("N must be odd, given: $N"))
        stab = Stabilities(P)
        tmps = BitVector(undef, N)
        return new(N, P, ξ, ξv, stab, tmps)
    end
end

"""
    GraphPercStep(N::Integer, P::Integer) <: SimpleGraph{Int}

A `SimpleGraph` implementing a single-layer binary perceptron with `N` binary (\$±1\$) synapses,
trained on `P` random i.i.d. \$±1\$ patterns.

The energy of the model is computed as the number of misclassified patterns.

See also [`GraphPercLinear`](@ref RRRMC.GraphPercLinear), [`GraphPercXEntr`](@ref RRRMC.GraphPercXEntr).
"""
GraphPercStep(N::Integer, P::Integer) = GraphPercStep(gen_ξ(N, P)...)

function Base.empty!(stab::Stabilities)
    @extract stab : p m Δs ξsi
    empty!(p)
    empty!(m)
    fill!(Δs, 0)
    stab.last_move = 0
    return stab
end

function energy(X::GraphPercStep, C::Config)
    @assert X.N == C.N
    @extract C : s
    @extract X : N P ξv stab tmps
    @extract stab : p m Δs

    E = 0
    empty!(stab)
    tmps = BitArray(undef, N)

    for a = 1:P
        map!(⊻, tmps, s, ξv[a])
        Δ = N - 2 * sum(tmps)
        Δs[a] = Δ
        if Δ == 1
            push!(p, a)
        elseif Δ < 0
            Δ == -1 && push!(m, a)
            E += 1
        end
    end

    return E
end

function update_cache!(X::GraphPercStep, C::Config, move::Int)
    @assert X.N == C.N
    @extract C : s
    @extract X : N P ξ stab
    @extract stab : p m Δs ξsi last_move
    @assert N == length(s)

    si = s[move]
    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    # @assert ξsi == ξ[:,move] ⊻ si
    @inbounds for a = 1:P
        oldΔ = Δs[a]
        newΔ = oldΔ + (2 - 4 * ξsi[a])

        if oldΔ > 1 && newΔ == 1
            push!(p, a)
        elseif oldΔ == 1
            delete!(p, a)
            newΔ == -1 && push!(m, a)
        elseif oldΔ == -1
            delete!(m, a)
            newΔ == 1 && push!(p, a)
        elseif oldΔ < -1 && newΔ == -1
            push!(m, a)
        end
        Δs[a] = newΔ
    end
    si && flipbits!(ξsi)

    return
end

function delta_energy_naive(X::GraphPercStep, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N P ξ stab
    @assert 1 ≤ move ≤ N

    spinflip!(C, move)
    E1 = energy(X, C)
    spinflip!(C, move)
    E0 = energy(X, C)
    return E1 - E0
end

function delta_energy(X::GraphPercStep, C::Config, move::Int)
    @assert C.N == X.N
    @extract C : s
    @extract X : N P ξ stab
    @extract stab : p m ξsi last_move
    @assert 1 ≤ move ≤ N

    si = s[move]
    last_move ≠ move && unsafe_copyto!(ξsi, 1, ξ, (move-1)*P + 1, P)
    stab.last_move = move
    si && flipbits!(ξsi)
    # @assert ξsi == ξ[:,move] ⊻ si
    ΔE = 0
    @inbounds for a in p
        ΔE += ~ξsi[a]
    end
    @inbounds for a in m
        ΔE -= ξsi[a]
    end
    si && flipbits!(ξsi)
    # td = delta_energy_naive(X, s, stab, move)
    # ΔE == td || @show ΔE, td
    # @assert ΔE == td
    return ΔE
end

function check_delta(X::GraphPercStep, C::Config, move::Int)
    delta = delta_energy(X, C, move)
    e0 = energy(X, s)
    spinflip!(C, move)
    e1 = energy(X, s)
    spinflip!(C, move)

    (e1-e0) == delta || (@show e1,e0,delta,e1-e0; error())
end

neighbors(X::GraphPercStep, i::Int) = AllButOne(X.N, i)

end


text:
function solve()
    a, b, c = [parse(Int, x) for x in split(readline())]
    if c == 0
        if a > b
            "Takahashi"
        else
            "Aoki"
        end
    else
        if a ≥ b
            "Takahashi"
        else
            "Aoki"
        end
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function invariantmeasure(pts, method::TriangulationEstimators.SimplexExact; tol = 1e-8)
    tog = transopergenerator(pts, method)
    to = tog(tol = tol)
    invariantmeasure(to)
end

function invariantmeasure(pts, method::TriangulationEstimators.SimplexPoint; tol = 1e-8, randomsampling::Bool = false, n::Int = 100)
    tog = transopergenerator(pts, method)
    to = tog(tol = tol, randomsampling = randomsampling, n = n)
    invariantmeasure(to)
end


text:
module TestShapeDiagrams
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra.ShapeDiagrams

A, B, C = Ob(FreeCategory, :A, :B, :C)

# Spans.
f, g = Hom(:f, C, A), Hom(:g, C, B)
span = Span(f,g)
@test apex(span) == C
@test left(span) == f
@test right(span) == g
f = Hom(:f, A, B)
@test_throws Exception Span(f,g)

# Cones.
f, g = Hom(:f, C, A), Hom(:g, C, B)
cone = Cone(C,[f,g])
@test apex(cone) == C
@test leg(cone,1) == f
@test leg(cone,2) == g
f = Hom(:f, A, B)
@test_throws Exception Cone(C,[f,g])

# Cospans.
f, g = Hom(:f, A, C), Hom(:g, B, C)
cospan = Cospan(f,g)
@test base(cospan) == C
@test left(cospan) == f
@test right(cospan) == g

# Cocones.
f, g = Hom(:f, A, C), Hom(:g, B, C)
cocone = Cocone(C,[f,g])
@test base(cocone) == C
@test leg(cocone,1) == f
@test leg(cocone,2) == g
f = Hom(:f, A, B)
@test_throws Exception Cocone(C,[f,g])

# Diagrams.

f, g, h = Hom(:f, A, C), Hom(:g, B, C), Hom(:h, A, B)
diag = Diagram([A,B,C],[(1,3,f),(2,3,g),(1,2,h)])
@test_throws Exception Diagram([A,B,C],[(1,2,f),(2,3,g),(1,2,h)])

end


text:
function generate_rand_cache(N, M; T=Float64, range=-100.0:100.0)
    cache = DS.PointCache{T}()
    for _ in 1:M
        DS.CachePush(cache, rand(range,N), rand(range))
    end
    return cache
end

@testset "Cache" begin
    T = Float64
    c = DS.PointCache{T}()
    @testset "Constructor" begin
        @test typeof(c.costs) == Dict{Vector{T},T}
        @test length(c.costs) == 0
        @test typeof(c.order) == Vector{Vector{T}}
        @test length(c.order) == 0
    end

    @testset "CachePush" begin
        @test_throws MethodError DS.CachePush(c, [1,2,3], 1)
        @test_throws MethodError DS.CachePush(c, [1.0,2.0,3.0], 1)
        @test_throws MethodError DS.CachePush(c, [1,2,3], 1.0)

        DS.CachePush(c, [1.0,2.0,3.0], 4.0)
        @test c.costs[[1.0,2.0,3.0]] == 4.0

        DS.CachePush(c, [1.0,2.0,3.0], 6.0)
        @test c.costs[[1.0,2.0,3.0]] == 6.0

        DS.CachePush(c, [4.21, 45.0, 1234321.0], 2112.0)
        @test c.costs[[1.0,2.0,3.0]] == 6.0
        @test c.costs[[4.21, 45.0, 1234321.0]] == 2112.0
    end

    @testset "CacheOrderPush" begin
        @test c.order == []
        DS.CacheOrderPush(c, nothing)
        @test c.order == []
        DS.CacheOrderPush(c, [1.0, 2.0, 3.0])
        @test c.order == [[1.0, 2.0, 3.0]]
        DS.CacheOrderPush(c, [2.0, 2.0, 2.0])
        @test c.order == [[1.0, 2.0, 3.0], [2.0, 2.0, 2.0]]
        DS.CacheOrderPush(c, [2.0, 2.0, 2.0])
        @test c.order == [[1.0, 2.0, 3.0], [2.0, 2.0, 2.0]]
        DS.CacheOrderPush(c, [1.0, 2.0, 3.0])
        @test c.order == [[1.0, 2.0, 3.0], [2.0, 2.0, 2.0],[1.0, 2.0, 3.0]]
    end

    @testset "CacheQuery" begin
        @test_throws MethodError DS.CacheQuery(c, [1,2,3])

        @test DS.CacheQuery(c, [1.0,2.0,3.0])
        @test !DS.CacheQuery(c, [1.0,1.0,1.0])
    end

    @testset "CacheGet" begin
        @test_throws MethodError DS.CacheGet(c, [1,2,3])

        @test DS.CacheGet(c, [4.21, 45.0, 1234321.0]) == 2112.0
        @test DS.CacheGet(c, [1.0,2.0,3.0]) ==  6.0
        @test_throws KeyError DS.CacheGet(c, [1.0,2.0,4.0])
     end

    @testset "CacheRandomSample" begin
        @test DS.CacheRandomSample(c, 0) == []
        @test DS.CacheRandomSample(c, 1)[1] in c.order
        @test length(DS.CacheRandomSample(c, 1)) == 1
        @test all(p in c.order for p in DS.CacheRandomSample(c, 10))
        @test length(DS.CacheRandomSample(c, 2)) == 2
    end

    @testset "CacheInitialPoint" begin
        @test_throws BoundsError DS.CacheInitialPoint(DS.PointCache{T}())
        @test DS.CacheInitialPoint(c) == ([1.0,2.0,3.0], 6.0)
    end

    @testset "CacheGetRange" begin
        @test DS.CacheGetRange(c, [[4.21, 45.0, 1234321.0], [1.0,2.0,3.0]]) == [2112.0, 6.0]
        @test DS.CacheGetRange(c, []) == []
        @test_throws KeyError DS.CacheGetRange(c, [[4.21, 45.0, 1234321.0], [1.0,2.0,4.0]])

    end

    @testset "CacheFilter" begin
        @test DS.CacheFilter(c, [[4.21, 45.0, 1234321.0], [1.0,2.0,3.0], [7.0,8.0,9.0]]) ==
              ([[4.21, 45.0, 1234321.0], [1.0,2.0,3.0]], [[7.0,8.0,9.0]])
        @test DS.CacheFilter(c, [[4.21, 45.0, 1234321.0], [1.0,2.0,3.0]]) ==
              ([[4.21, 45.0, 1234321.0], [1.0,2.0,3.0]], [])
        @test DS.CacheFilter(c, [[4.22, 68.0, -1114754.0], [8.0,1.0,2.0]]) ==
              ([], [[4.22, 68.0, -1114754.0], [8.0,1.0,2.0]])
        @test DS.CacheFilter(c, []) == ([], [])

    end
end


text:
const LYCEUM_REGISTRY = RegistrySpec(
    uuid = "96ca1a0e-015d-46b7-a12d-07d1320642ca",
    name = "LyceumRegistry",
    url = "https://github.com/Lyceum/LyceumRegistry.git",
)

const GENERAL_REGISTRY = RegistrySpec(
    uuid = "23338594-aafe-5451-b93e-139f81909106",
    name = "General",
    url = "https://github.com/JuliaRegistries/General.git",
)

const LYCEUM_PKGTEMPLATE = (
    user = "Lyceum",
    host = "github.com",
    license = "BSD3",
    authors = ["Colin Summers"],
    ssh = true,
    manifest = true,
    julia_version = VERSION,
    dev = false,
)

function namedgenerate(pkgname::String, preset::NamedTuple = LYCEUM_PKGTEMPLATE; kwargs...)
    push!(preset.authors, "The Contributors of $pkgname")
    t = Template(; preset..., kwargs...)
    generate(pkgname, t)
end

function lyceum_register(package_repo, commit, lyceumbot_pat)
    register(package_repo, commit, LYCEUM_REGISTRY.url, GitHub.authenticate(lyceumbot_pat), "lyceumbot")
end

function register(package_repo::AbstractString, commit::AbstractString, registry::AbstractString, auth::GitHub.Authorization, gituser::AbstractString; ignore_dirty::Bool=false, kwargs...)
    try
        Base.SHA1(commit)
    catch e
        throw(ArgumentError("\"$commit\" is not a valid SHA-1"))
    end

    package_repo = githttpsurl(package_repo) # always register https url

    cdtempdir() do
        gitconfig=Dict("user.name"=>gituser, "user.password"=>auth.token)
        git = DevTools.create_git_cmd(pairs(gitconfig)...)

        run(`$git clone $package_repo $(pwd())`)
        run(`$git checkout master`)
        success(`$git merge-base --is-ancestor $commit HEAD`) || error("$commit not found in master branch")
        run(`$git checkout $commit`)

        pkg = Pkg.Types.read_project(Pkg.Types.projectfile_path(pwd()))
        tree_hash = bytes2hex(Pkg.GitTools.tree_hash(pwd()))
        registry_deps = map(reg->reg.url, Pkg.Types.collect_registries())
        if !(registry in registry_deps)
            error("$registry not found in local registries")
        end

        rbrn = RegistryTools.register(package_repo, pkg, tree_hash; kwargs..., registry=registry, registry_deps=registry_deps, gitconfig=gitconfig, push=true, branch="master")
        haskey(rbrn.metadata, "error") && error(rbrn.metadata["error"])

        #name = rbrn.name
        #ver = rbrn.version
        #brn = rbrn.branch
        #meta = rbrn.metadata

        #repo = GitHub.repo(join(split(splitext(registry)[1], "/")[end-1:end], "/"), auth=auth)
        #params = Dict(
        #    "base"=>"master",
        #    "head"=>brn,
        #    "maintainer_can_modify"=>true,
        #    "title"=> "$(rbrn.metadata["kind"]): $name v$ver",
        #)
        #haskey(meta, "labels") && (params["labels"] = meta["labels"])
        #body = IOBuffer()
        #write(body, "Commit: $commit")
        #if get(rbrn.metadata, "warning", nothing) !== nothing
        #    write(body,
        #        """
        #        Warning: $(rbrn.metadata["warning"])
        #        This can be safely ignored. However, if you want to fix this you can do so. Call register() again after making the fix. This will update the Pull request.
        #        """
        #    )
        #end
        #params["body"] = String(take!(body))

        #pr, msg = create_or_find_pull_request(repo, auth, params)

        #return pr, rbrn
        return rbrn
    end
end

function create_or_find_pull_request(repo::GitHub.Repo, auth::GitHub.Authorization, params::Dict{<:AbstractString, Any})
    pr = nothing
    msg = ""
    try
        pr = GitHub.create_pull_request(repo; auth=auth, params=params)
        msg = "created"
        @debug("Pull request created")
    catch ex
        if is_pr_exists_exception(ex)
            @debug("Pull request already exists, not creating")
            msg = "updated"
        else
            rethrow(ex)
        end
    end

    if pr === nothing
        prs, _ = GitHub.pull_requests(repo; auth=auth, params=Dict(
            "state" => "open",
            "base" => params["base"],
            "head" => string(repo.owner.login, ":", params["head"]),
        ))
        if !isempty(prs)
            @assert length(prs) == 1 "PR lookup should only contain one result"
            @debug("PR found")
            pr = prs[1]
        end

        if pr === nothing
            error("PR already exists but unable to find it")
        else
            GitHub.update_pull_request(repo, pr; auth=auth, params=params)
        end
    end

    try # update labels
        if haskey(params, "labels")
            GitHub.edit_issue(repo, pr; auth = auth, params = Dict("labels"=>params["labels"]))
        end
    catch
        @debug "Failed to update labels, ignoring."
    end

    return pr, msg
end

function parse_github_exception(ex::ErrorException)
    msgs = map(strip, split(ex.msg, '\n'))
    d = Dict()
    for m in msgs
        a, b = split(m, ":"; limit=2)
        d[a] = strip(b)
    end
    return d
end

function is_pr_exists_exception(ex)
    d = parse_github_exception(ex)

    if d["Status Code"] == "422" &&
       match(r"A pull request already exists", d["Errors"]) !== nothing
        return true
    end

    return false
end

incrementversion!(pkg::Module, args...; kwargs...) =
    incrementversion!(pkgdir(pkg), args...; kwargs...)
function incrementversion!(pkgdir::String, args...; tag::Bool = false, kwargs...)
    isdirty(pkgdir) && error("$pkgdir is dirty")

    toml = parsetomls(pkgdir).project.dict
    newver = incrementversion(toml["version"], args...; kwargs...)
    toml["version"] = newver
    tomlpath = Pkg.Operations.projectfile_path(pkgdir)
    @info "Writing to $tomlpath"
    open(tomlpath, "w") do io
        Pkg.TOML.print(io, toml)
    end

    git = create_git_cmd(path = pkgdir)
    run(`$git add Project.toml`)
    message = "New version: v$(newver)"
    run(`$git commit -qm $message`)

    pkgdir
end
function incrementversion(
    version::Union{String,VersionNumber},
    which::Symbol;
    prerelease = :keep,
)
    v = version isa String ? VersionNumber(version) : version
    major = v.major
    minor = v.minor
    patch = v.patch
    pre = v.prerelease
    build = v.build

    build != () && @warn "Build not empty: $build"

    if which === :major
        major += 1
        minor = 0
        patch = 0
    elseif which === :minor
        minor += 1
        patch = 0
    elseif which === :patch
        patch += 1
    else
        error("which must be :major, :minor, or :patch. Got $which")
    end

    if prerelease === :keep
        nothing
    elseif prerelease === :dev
        pre = ("DEV",)
    elseif isnothing(prerelease)
        pre = ()
    else
        error("prerlease must be one of :keep, :dev, or nothing")
    end

    vnew = VersionNumber(major, minor, patch, pre, build)
    @info "Old version: $v. New version: $vnew"
    vnew
end

function isversionlatest(pkg::Module)
    project = parsetomls(pkgdir(pkg)).project.dict
    localver = VersionNumber(project["version"])
    latestver = latestversion(pkg)
    return isnothing(latestversion) || localver > latestver
end

latestversion(pkg::Module) = latestversion(UUID(parsetomls(pkgdir(pkg)).project.dict["uuid"]))
function latestversion(uuid::UUID)
    Pkg.Registry.update()
    registries = Pkg.Types.collect_registries()
    for reg in registries
        packages = Pkg.Types.read_registry(joinpath(reg.path, "Registry.toml"))["packages"]
        if haskey(packages, string(uuid))
            package_path = joinpath(reg.path, packages[string(uuid)]["path"])
            versions = collect(keys(Pkg.Operations.load_versions(package_path)))
            sort!(versions)
            return versions[end]
        end
    end
    nothing
end


text:
# ============================================================
# Approximation model in discrete time, i.e. without bloating
# ============================================================

"""
    NoBloating{EM, SO, IT} <: AbstractApproximationModel

No bloating, or discrete-time, approximation model.

### Fields

- `exp`     -- exponentiation method
- `setops`  -- set operations method
- `inv`     -- (optional, default: `false`) if `true`, assume that the state matrix
               is invertible and use its inverse in the `Φ` functions

### Algorithm

The transformations are:

- ``Φ ← \\exp(Aδ)``
- ``Ω_0 ← \\mathcal{X}_0``
- ``V(k) ← Φ₁(A, δ)U(k)``, ``k ≥ 0``.

The function ``Φ₁(A, δ)`` is defined in [`Φ₁`](@ref).
We allow ``U`` to be a sequence of time varying non-deterministic input sets.

See also Eqs.(14) in [[BFFPSV18]](@ref).
"""
struct NoBloating{EM, SO, IT} <: AbstractApproximationModel
    exp::EM
    setops::SO
    inv::IT
end

# convenience constructor using symbols
function NoBloating(; exp=BaseExp, setops=:lazy, inv=false)
    return NoBloating(_alias(exp), _alias(setops), Val(inv))
end

function Base.show(io::IO, alg::NoBloating)
    print(io, "`NoBloating` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
    print(io, "    - set operations method: $(alg.setops)\n")
    print(io, "    - invertibility assumption: $(alg.inv)")
end

Base.show(io::IO, m::MIME"text/plain", alg::NoBloating) = print(io, alg)

# homogeneous case
function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ, alg::NoBloating)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)

    Φ = _exp(A, δ, alg.exp)

    Ω0 = copy(X0) # alg.setops doesn't apply
    X = stateset(ivp)
    Sdiscr = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdiscr, Ω0)
end

# inhomogeneous case
function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ, alg::NoBloating)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    U = next_set(inputset(ivp), 1)

    Φ = _exp(A, δ, alg.exp)
    if isa(U, AbstractSingleton)
        Mu = _Φ₁_u(A, δ, alg.exp, alg.inv, element(U), Φ)
        V = Singleton(Mu)
    else
        M = _Φ₁(A, δ, alg.exp, alg.inv, Φ)
        V = _apply_setops(M * U, alg.setops)
    end

    Ω0 = _initial_state(X0)

    In = IdentityMultiple(one(eltype(A)), size(A, 1))
    X = stateset(ivp)
    Sdis = ConstrainedLinearControlDiscreteSystem(Φ, In, X, V)
    return InitialValueProblem(Sdis, Ω0)
end

# extension for no bloating case and singleton initial states
_initial_state(X0::CartesianProduct{N, <:Singleton{N}, <:Singleton{N}}) where {N} = convert(Singleton, X0)
_initial_state(X0) = X0


text:
using DEBTrait
using Distributions
using CSV, DataFrames
using Random
Random.seed!(123)

n_microbes          = 100
n_substrates        = 1
#
L_DNA               = collect(range(1e6, 8e6, length=n_microbes))
rRNA                = DEBTrait.genome_size_to_rRNA_copy_number(L_DNA)
min_gentime_distr   = Gamma(2.486, 0.97)
min_gentime         = rand(min_gentime_distr, n_microbes)
Isolate             = BaseGenome(L_DNA, rRNA, min_gentime, 1.0*ones(n_substrates,n_microbes))
D_S                 = DEBTrait.aqueous_diffusivity(59.04402*ones(1))
Monomer             = BaseMonomer(D_S[1]*ones(n_substrates), 6*ones(n_substrates), 3*ones(n_substrates), 0.3*ones(n_substrates))
Soil                = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC         = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
MetabolismC.y_EM    = ones(n_microbes)
MetabolismC.y_EX    = ones(n_microbes)
MetabolismC.α       = zeros(n_microbes)
AssimilationC       = PAssimilation(Isolate, Monomer, Soil, MetabolismC)

Acetate_dSize        = DataFrame()
Acetate_dSize.V_cell = DEBTrait.genome_size_to_cell_volume(L_DNA)
Acetate_dSize.min_gentime = min_gentime
Acetate_dSize.transporter_density = AssimilationC.ρ_p[1,:]
Acetate_dSize.K      = AssimilationC.K_D_0[1,:]
Acetate_dSize.V_max = AssimilationC.N_SB[1,:]*100

CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/biocrunch/Acetate_cell_size_variation.csv", Acetate_dSize)

#
n_microbes          = 1
n_substrates        = 1
L_DNA               = 3726411*ones(1)
rRNA                = DEBTrait.genome_size_to_rRNA_copy_number(L_DNA)
min_gentime         = rand(min_gentime_distr, n_microbes)
Isolate             = BaseGenome(L_DNA, rRNA, min_gentime, 1.0*ones(n_substrates,n_microbes))
D_S                 = DEBTrait.aqueous_diffusivity(59.04402*ones(1))
Monomer             = BaseMonomer(D_S, 6*ones(n_substrates), 3*ones(n_substrates), 0.3*ones(n_substrates))
Soil                = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
MetabolismC         = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
MetabolismC.y_EM    = ones(n_microbes)
MetabolismC.y_EX    = ones(n_microbes)
MetabolismC.α       = zeros(n_microbes)
AssimilationC       = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
#
Soil.pct_sand       = 70.0
Soil.pct_clay       = 30.0
Soil.N_cell         = 10.0*ones(n_microbes)

K_D_all = zeros(100)
s_sats = collect(range(0.01, 1.0, length=100))
for i in 1:100
    Soil.s_sat          = s_sats[i]
    AssimilationC       = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
    K_D_all[i]          = AssimilationC.K_D[1]
end

Acetate_Geobacter_dSat        = DataFrame()
Acetate_Geobacter_dSat.V_cell = fill(L_DNA[1], 100)
Acetate_Geobacter_dSat.min_gentime = fill(min_gentime[1], 100)
Acetate_Geobacter_dSat.transporter_density = fill(AssimilationC.ρ_p[1], 100)
N_SB = AssimilationC.N_SB[1,:]*100
Acetate_Geobacter_dSat.V_max = fill(N_SB[1], 100)
Acetate_Geobacter_dSat.s_sat = s_sats
Acetate_Geobacter_dSat.K = K_D_all

CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/biocrunch/Acetate_Geobacter_s_sat_variation.csv", Acetate_Geobacter_dSat)


text:
import Graphs: edge_index, source, target

type BayesianEdge
    index::Int
    source::BayesianNode
    target::BayesianNode

    BayesianEdge(s::BayesianNode, t::BayesianNode) = new(0,s,t)
end

function ==(e1::BayesianEdge, e2::BayesianEdge)
    e1.index == e2.index && e1.source == e2.source && e1.target == e2.target
end

BayesianEdge{T <: BayesianNode}(s::CBayesianNode, t::T) = throw("Source cannot be of type CBayesianNode")

edge_index(e::BayesianEdge) = e.index

source(e::BayesianEdge) = e.source
target(e::BayesianEdge) = e.target


text:
module Basis
include("abstractomega.jl")

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 8585ad90-d99f-11eb-0782-8f1b33d43b4d
begin
    import Pkg
    # activate a clean environment
    Pkg.activate(mktempdir())

    Pkg.add(name="Images")
    Pkg.add(name="FileIO")

	using Images
	using FileIO
end

# ╔═╡ 0c8998d9-d8c1-40c3-b147-f44003fef025
begin
	url = "https://i.imgur.com/VGPeJ6s.jpg"
	download(url, "philip.jpg")
	philip = load("philip.jpg");
end

# ╔═╡ 03a8f21d-c985-49b2-9a5e-9bd4355b2d72
decimate(img, k) = img[1:k:end, 1:k:end]

# ╔═╡ 5e7d5bf3-147e-44d1-9c25-8e5a830087e2
poor_phil = decimate(philip, 8)

# ╔═╡ 64857813-e2f0-45f3-992c-e3cb55a1bfaf
function clamp_at_boundary(M, i, j)
	return M[
		clamp(i, 1, size(M, 1)),
		clamp(j, 1, size(M, 2)),
	]
end

# ╔═╡ 568094ea-255d-4100-9c60-2641b78afddd
function convolve(M, kernel, M_index_func=clamp_at_boundary)
	size(M)
    height = size(kernel, 1)
    width = size(kernel, 2)

    half_height = height ÷ 2
    half_width = width ÷ 2
    new_image = similar(M)

	# (i, j) loop over the original image
    @inbounds for i in 1:size(M, 1)
        for j in 1:size(M, 2)
            # (k, l) loop over the neighbouring pixels
			new_image[i, j] = sum([
				kernel[k, l] * M_index_func(M, i - k, j - l)
				for k in -half_height:-half_height + height - 1
				for l in -half_width:-half_width + width - 1
			])
        end
    end

    return new_image
end

# ╔═╡ ece09f3d-9c33-486e-95b2-3a723ef27dc5
box_blur(n) = centered(ones(n, n) ./ (n^2))

# ╔═╡ f9bcaeea-6aaa-4ca2-9a33-f92dea365984
# convolve(poor_phil, box_blur(2))
size(poor_phil)

# ╔═╡ bae786fe-1bba-46e1-8728-778ad0bc349d
kernel = box_blur(5) # Kernel.gaussian((1, 0.1))

# ╔═╡ cb5a3e07-db6b-4751-9480-4f1c07ac7e95
sum(kernel)

# ╔═╡ cf71d105-e775-424f-9da8-9c552ede6fe7
convolve(poor_phil, kernel)

# ╔═╡ 2f573afa-0e25-453c-89ef-7e332252fe42

# ╔═╡ Cell order:
# ╠═8585ad90-d99f-11eb-0782-8f1b33d43b4d
# ╠═0c8998d9-d8c1-40c3-b147-f44003fef025
# ╠═5e7d5bf3-147e-44d1-9c25-8e5a830087e2
# ╟─03a8f21d-c985-49b2-9a5e-9bd4355b2d72
# ╟─64857813-e2f0-45f3-992c-e3cb55a1bfaf
# ╠═568094ea-255d-4100-9c60-2641b78afddd
# ╠═ece09f3d-9c33-486e-95b2-3a723ef27dc5
# ╠═f9bcaeea-6aaa-4ca2-9a33-f92dea365984
# ╠═bae786fe-1bba-46e1-8728-778ad0bc349d
# ╠═cb5a3e07-db6b-4751-9480-4f1c07ac7e95
# ╠═cf71d105-e775-424f-9da8-9c552ede6fe7
# ╠═2f573afa-0e25-453c-89ef-7e332252fe42

max_stars_repo_name: codewithkaranjeswani/my-pluto-notebooks

max_issues_repo_name: codewithkaranjeswani/my-pluto-notebooks

max_forks_repo_name: codewithkaranjeswani/my-pluto-notebooks


text:
using Flux
struct SirenDense
    weight
    bias
    ω0
end

Flux.@functor SirenDense

function SirenDense(din, dout; c=6, ω0=1, isfirst=false)
    w_std = isfirst ? (1 / din) : (sqrt(c / din) / ω0)
    dist = Uniform(-w_std, w_std)
    weight = Float32.(rand(dist, dout, din))
    bias = Float32.(rand(dist, dout))
    SirenDense(weight, bias, ω0)
end

function (a::SirenDense)(x)
  sin.(a.ω0 * (a.weight*x .+ a.bias))
end

struct ModulatedSiren
    sirens
    modulator
end

Flux.trainable(π::ModulatedSiren) = (Flux.trainable(π.sirens)..., Flux.trainable(π.modulator)...)

function (a::ModulatedSiren)(x, z)
    slayers = layers(a.sirens)
    mlayers = layers(a.modulator)
    @assert length(slayers) == 4
    @assert length(slayers) == length(mlayers) + 1

    h = z

    h = mlayers[1](h)
    x = slayers[1](x) .* h
    h = vcat(h, z)

    h = mlayers[2](h)
    x = slayers[2](x) .* h
    h = vcat(h, z)

    h = mlayers[3](h)
    x = slayers[3](x) .* h
    # h = vcat(h, z)
    #
    # h = mlayers[4](h)
    # x = slayers[4](x) .* h
    # h = vcat(h, z)
    #
    # h = mlayers[5](h)
    # x = slayers[5](x) .* h

    slayers[end](x)
end


text:
function kmeans(
    X::AbstractMatrix{<:Real},
    k::Integer;
    T::Float64 = 0.98,
    iseeds = randomcenters
)
    d, n = size(X)

    # assign initial cluster centers
    centers = iseeds(X, k)
    assignments = Vector{Int}(undef, n)

    dist = pairwise(Euclidean(), centers, X, dims = 2)
    D = eltype(dist)

    costs = Vector{D}(undef, n)
    counts = Vector{Int}(undef, k)
    sd = Array{D}(undef, d, k)
    iter = 0

    while true
        iter += 1
        if iter > 1000
            return ClusteringResult(centers, assignments, costs, counts, sd, iter)
        end

        dist = pairwise(Euclidean(), centers, X, dims = 2)
        prev_assignments = deepcopy(assignments)
        update_assignments!(dist, assignments, costs, counts)
        update_centers!(X, assignments, centers, counts)

        # check convergence threshold
        updates = assignments .== prev_assignments
        if (sum(updates)/length(updates)) > T
            return ClusteringResult(centers, assignments, costs, counts, sd, iter)
        end

        # compute std dev of each component in each cluster
        sd = Array{D}(undef, d, k)
        update_sd!(X, assignments, centers, counts, sd)
    end
end


text:
# enemydbcrud.jl (v0.2)

export createenemy
export insert!
export update!

export findenemies!
export deleteall!
export deletefirst!
export deletenext!

export printenemy

Enemy = Dict{String, Any};
export Enemy

Enemyarray = Array{Enemy, 1}
export Enemyarray

struct Enemy2

	uq::String
	name::String
	hp::Int32
	ac::Int32
	abl::Vector{Int32}
	att::Vector{String}
	loot::Vector{String}

	Enemy2() = new("null", "no name", 0, 0, [0, 0, 0, 0, 0, 0], "no attacks", "no loot")

end

struct Enemy3

	data::Dict{String, Any}
	Enemy3() = new(Dict())
	Enemy3(data::Dict{String, Any}) = new(data)

end

# create, read, update, destroy
#=
	createenemy()
		creates and enemy, either from console input or parameterlist
	insert!()
		inserts the enemy such that it appears in both
		lookup and enemytable
	findenemies!()
		returns a list of all enemies satisfying search criteria
		these are references to occurrences in the inputted array
	findall(), findfirst(), findnext()
		use functions implemented for Array{Dict{String, Any}, 1}
	update!()
		changes value of given enemy's key according to parameter
	deleteall(), deletefirst(), deletenext()
=#

function createenemy()
	d = Enemy();

	d["unique"] = cin(String, 1, "Give unique key: ");
	d["name"] = cin(String, 1, "Give enemy name: ");
	d["ac"] = cin(Int32, 1, "Give AC: ");
	d["hp"] = cin(Int32, 1, "Give HP: ");
	d["abl"] = cin(Int32, 6, "Give abilityscores, newline separates: \n");
	d["attacks"] = cin(String, -1, "Attacks: \n");
	d["loot"] = cin(String, -1, "Loot: \n");

	return d;

end

function createenemy(unique::String, name::String, ac::Int, hp::Int, abl::Array{Int, 1}, attacks::Array{String, 1}, loot::Array{String, 1})::Enemy

	d = Enemy();

	d["unique"] = unique;
	d["name"] = name;
	d["ac"] = ac;
	d["hp"] = hp;
	d["abl"] = abl;
	d["attacks"] = attacks;
	d["loot"] = loot;

	return d;

end

function insert!(lookup::Array{String, 1}, table::Enemyarray, e::Enemy)

	if haskey(e, "unique")
		push!(lookup, e["unique"]);
		push!(table, e);
	else
		println("Enemy had no unique key.");
	end

end

# Marked with ! because it returns a pointer to the array-element
# DEPRECATED
function findenemies!(f::Function, arr::Enemyarray)

	inds = findall(f, arr);

	enem::Enemyarray;

	for i in inds
		push!(enem, arr[i]);
	end

	return enem;

end

function update!(e::Enemy, key::String, value::Any)

	e[key] = value;

end

deletefirst!(f::Function, arr::Enemyarray) = deleteat!(arr, findfirst(f, arr));
deletenext!(f::Function, arr::Enemyarray, i::Int) = deleteat!(arr, findnext(f, arr, i));
deleteall!(f::Function, arr::Enemyarray) = deleteat!(arr, findall(f, arr));

function printenemy(e::Enemy)

	println(e["unique"]);
	println("Name: " * string(e["name"]));
	println("HP: " * string(e["hp"]) * ", AC: " * string(e["ac"]));
	println(Vector{Int}(e["abl"]));

	println("Attacks: ");
	for att in e["attacks"]
		println("\t" * att);
	end

	println("Loot: ");
	for loot in e["loot"]
		println("\t" * loot);
	end

end


text:
using Pkg

# Environment Stack knowledge from https://docs.julialang.org/en/v1/manual/code-loading/#Environment-stacks
DEVENV_LOC = joinpath(homedir(), ".julia", "config", "devenv")
if !(DEVENV_LOC in LOAD_PATH)
    @info "Adding $DEVENV_LOC to LOAD_PATH"
    push!(LOAD_PATH, DEVENV_LOC)
end

Pkg.activate(".")

try
    using Revise
catch e
    @warn "Error initializing Revise" exception=(e, catch_backtrace())
end

# utility
ls(x) = readdir(x)
ls() = readdir()
ty(x) = typeof(x)
fn(x) = fieldnames(x)
fnty = fn ∘ ty

unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))

function kill_julia_workers()
    julia_worker_procs = filter(x->occursin("--bind-to", x) && occursin("--worker", x), split(read(`pgrep -af julia`, String), "\n"))
    just_worker_pids = map(x->split(x, " ")[1], julia_worker_procs)
    worker_kill_command = `kill $just_worker_pids`
    (julia_worker_procs, worker_kill_command)
end


text:
# Machine Learning
# Data Preprocessing Array functions

using CSV
using Statistics
using Missings

"""
```
train_test_split(x::Array{T}; split_ratio::Float64=.7)
```
Split Data Into Train and Test Sets
"""
function train_test_split(x::Array{T}; split_ratio::Float64=.7)::Array{T} where {T<:Real}
    let len = length(len), split_no = Int64.(round(split_ratio * len))
    end
    train_set = data[1:split_no]
    test_set = data[split_no+1:len]
    return train_set,test_set
end

"""
```
standardization(x::Array{T})
```
Feature Scaling
standardization = x - mean(x) / standard deviation(x)
"""
function standardization(x::Array{T})::Array{T} where {T<:Real}
    let out = fill(NaN, size(x,1)), avg = mean(x), stdev = std(x)
    @inbounds for i = 1:size(x,1)
        out[i] = (x[i] - avg) / stdev
end
    return out
end
end

"""
```
normalization(x::Array{T})
```
Feature Scaling
normalization = x - min(x) / max(x) - min(x)
"""
function normalization(x::Array{T})::Array{T} where {T<:Real}
    let out = fill(NaN, size(x,1)), min = minimum(x), max = maximum(x)
    @inbounds for i = 1:size(x,1)
        out[i] = ((x[i] - min) / (max - min))
end
    return out
end
end

"""
```
data_fill_nan_inf(x::Array{T})
```
Convert NaN to Inf
"""
function data_fill_nan_inf(x::Array{T})::Array{T} where {T<:Real}
        @inbounds for i = 1:size(x,1)
        if (isnan(x[i]))
                x[i] = Inf
            end
end
return x
end

x = data_fill_nan_inf(x)

"""
```
fill_missing(x::Array{T}; method::String="mean")
```
Run function data_fill_nan_inf() this function thens replace Inf with the mean, median or previous values
"""
function fill_missing(x::Array{T}; method::String="mean")::Array{T} where {T<:Real}
    y = x[x .!= Inf] # Subset to remove InF, NaN
        @inbounds for i = 1:size(x,1)
        if method == "mean"
            avg = mean(y)
            if x[i] == Inf
                x[i] = avg
            end
        elseif method == "median"
            med = median(y)
            if x[i] == Inf
                x[i] = med
            end
        elseif method == "prev" && i >= 1
            prev = x[i-1]
            if x[i] == Inf
                x[i] = prev
            end
        elseif method == "prev" && i == 1
            prev = x[i+1]
            if x[i] == Inf
                x[i] = prev
            end
end
end
return x
end

# Data preprocessing - working with dataframes

# Read data
dataset = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Machine Learning/Machine-Learning/Machine-Learning/data_preprocessing/Data.csv", header=true,types=[String, Float64, Float64, String])

# Replace missing values with the mean
using Missings
dataset[:Age] = collect(Missings.replace(dataset[:Age], mean(skipmissing(dataset[:Age]))))
dataset[:Salary] = collect(Missings.replace(dataset[:Salary], mean(skipmissing(dataset[:Salary]))))

# One hot encoding categroial variables
for i in unique(dataset[:Country])
    dataset[Symbol(i)] = ifelse.(dataset[:Country] .== i, 1, 0)
end

dataset[:Purchased] = ifelse.(dataset[:Purchased] .== "Yes", 1, 0)

# Feature Scaling - note in order to not have multiple variables over powering each other in the model
dataset[:Age] = standardization(dataset[:Age])
dataset[:Salary] = standardization(dataset[:Salary])

# Split dataframe to train and test sets
split_ratio = .7
len = size(dataset,1)
split_no = Int64.(round(split_ratio * len))
train_set = dataset[1:split_no,:]
test_set = dataset[split_no+1:len,:]

max_stars_repo_name: JuliaTagBot/Machine-Learning

max_issues_repo_name: JuliaTagBot/Machine-Learning

max_forks_repo_name: JuliaTagBot/Machine-Learning