text:
@present ExtendedFreeExtCalc2D(FreeExtCalc2D) begin
  X::Space
  proj₁_⁰⁰₀::Hom(Form0(X)⊗Form0(X),Form0(X))
  proj₂_⁰⁰₀::Hom(Form0(X)⊗Form0(X),Form0(X))
  proj₁_⁰⁰₀⁺::Hom(Form0(X)⊕Form0(X),Form0(X))
  proj₂_⁰⁰₀⁺::Hom(Form0(X)⊕Form0(X),Form0(X))
  proj₁_⁰¹₀::Hom(Form0(X)⊗Form1(X),Form0(X))
  proj₂_⁰¹₁::Hom(Form0(X)⊗Form1(X),Form1(X))
  proj₁_⁰¹₀⁺::Hom(Form0(X)⊕Form1(X),Form0(X))
  proj₂_⁰¹₁⁺::Hom(Form0(X)⊕Form1(X),Form1(X))
  proj₁_⁰²₀::Hom(Form0(X)⊗Form2(X),Form0(X))
  proj₂_⁰²₂::Hom(Form0(X)⊗Form2(X),Form2(X))
  proj₁_⁰²₀⁺::Hom(Form0(X)⊕Form2(X),Form0(X))
  proj₂_⁰²₂⁺::Hom(Form0(X)⊕Form2(X),Form2(X))
  proj₁_⁰⁰̃₀::Hom(Form0(X)⊗DualForm0(X),Form0(X))
  proj₂_⁰⁰̃₀̃::Hom(Form0(X)⊗DualForm0(X),DualForm0(X))
  proj₁_⁰⁰̃₀⁺::Hom(Form0(X)⊕DualForm0(X),Form0(X))
  proj₂_⁰⁰̃₀̃⁺::Hom(Form0(X)⊕DualForm0(X),DualForm0(X))
  proj₁_⁰¹̃₀::Hom(Form0(X)⊗DualForm1(X),Form0(X))
  proj₂_⁰¹̃₁̃::Hom(Form0(X)⊗DualForm1(X),DualForm1(X))
  proj₁_⁰¹̃₀⁺::Hom(Form0(X)⊕DualForm1(X),Form0(X))
  proj₂_⁰¹̃₁̃⁺::Hom(Form0(X)⊕DualForm1(X),DualForm1(X))
  proj₁_⁰²̃₀::Hom(Form0(X)⊗DualForm2(X),Form0(X))
  proj₂_⁰²̃₂̃::Hom(Form0(X)⊗DualForm2(X),DualForm2(X))
  proj₁_⁰²̃₀⁺::Hom(Form0(X)⊕DualForm2(X),Form0(X))
  proj₂_⁰²̃₂̃⁺::Hom(Form0(X)⊕DualForm2(X),DualForm2(X))
  proj₁_¹⁰₁::Hom(Form1(X)⊗Form0(X),Form1(X))
  proj₂_¹⁰₀::Hom(Form1(X)⊗Form0(X),Form0(X))
  proj₁_¹⁰₁⁺::Hom(Form1(X)⊕Form0(X),Form1(X))
  proj₂_¹⁰₀⁺::Hom(Form1(X)⊕Form0(X),Form0(X))
  proj₁_¹¹₁::Hom(Form1(X)⊗Form1(X),Form1(X))
  proj₂_¹¹₁::Hom(Form1(X)⊗Form1(X),Form1(X))
  proj₁_¹¹₁⁺::Hom(Form1(X)⊕Form1(X),Form1(X))
  proj₂_¹¹₁⁺::Hom(Form1(X)⊕Form1(X),Form1(X))
  proj₁_¹²₁::Hom(Form1(X)⊗Form2(X),Form1(X))
  proj₂_¹²₂::Hom(Form1(X)⊗Form2(X),Form2(X))
  proj₁_¹²₁⁺::Hom(Form1(X)⊕Form2(X),Form1(X))
  proj₂_¹²₂⁺::Hom(Form1(X)⊕Form2(X),Form2(X))
  proj₁_¹⁰̃₁::Hom(Form1(X)⊗DualForm0(X),Form1(X))
  proj₂_¹⁰̃₀̃::Hom(Form1(X)⊗DualForm0(X),DualForm0(X))
  proj₁_¹⁰̃₁⁺::Hom(Form1(X)⊕DualForm0(X),Form1(X))
  proj₂_¹⁰̃₀̃⁺::Hom(Form1(X)⊕DualForm0(X),DualForm0(X))
  proj₁_¹¹̃₁::Hom(Form1(X)⊗DualForm1(X),Form1(X))
  proj₂_¹¹̃₁̃::Hom(Form1(X)⊗DualForm1(X),DualForm1(X))
  proj₁_¹¹̃₁⁺::Hom(Form1(X)⊕DualForm1(X),Form1(X))
  proj₂_¹¹̃₁̃⁺::Hom(Form1(X)⊕DualForm1(X),DualForm1(X))
  proj₁_¹²̃₁::Hom(Form1(X)⊗DualForm2(X),Form1(X))
  proj₂_¹²̃₂̃::Hom(Form1(X)⊗DualForm2(X),DualForm2(X))
  proj₁_¹²̃₁⁺::Hom(Form1(X)⊕DualForm2(X),Form1(X))
  proj₂_¹²̃₂̃⁺::Hom(Form1(X)⊕DualForm2(X),DualForm2(X))
  proj₁_²⁰₂::Hom(Form2(X)⊗Form0(X),Form2(X))
  proj₂_²⁰₀::Hom(Form2(X)⊗Form0(X),Form0(X))
  proj₁_²⁰₂⁺::Hom(Form2(X)⊕Form0(X),Form2(X))
  proj₂_²⁰₀⁺::Hom(Form2(X)⊕Form0(X),Form0(X))
  proj₁_²¹₂::Hom(Form2(X)⊗Form1(X),Form2(X))
  proj₂_²¹₁::Hom(Form2(X)⊗Form1(X),Form1(X))
  proj₁_²¹₂⁺::Hom(Form2(X)⊕Form1(X),Form2(X))
  proj₂_²¹₁⁺::Hom(Form2(X)⊕Form1(X),Form1(X))
  proj₁_²²₂::Hom(Form2(X)⊗Form2(X),Form2(X))
  proj₂_²²₂::Hom(Form2(X)⊗Form2(X),Form2(X))
  proj₁_²²₂⁺::Hom(Form2(X)⊕Form2(X),Form2(X))
  proj₂_²²₂⁺::Hom(Form2(X)⊕Form2(X),Form2(X))
  proj₁_²⁰̃₂::Hom(Form2(X)⊗DualForm0(X),Form2(X))
  proj₂_²⁰̃₀̃::Hom(Form2(X)⊗DualForm0(X),DualForm0(X))
  proj₁_²⁰̃₂⁺::Hom(Form2(X)⊕DualForm0(X),Form2(X))
  proj₂_²⁰̃₀̃⁺::Hom(Form2(X)⊕DualForm0(X),DualForm0(X))
  proj₁_²¹̃₂::Hom(Form2(X)⊗DualForm1(X),Form2(X))
  proj₂_²¹̃₁̃::Hom(Form2(X)⊗DualForm1(X),DualForm1(X))
  proj₁_²¹̃₂⁺::Hom(Form2(X)⊕DualForm1(X),Form2(X))
  proj₂_²¹̃₁̃⁺::Hom(Form2(X)⊕DualForm1(X),DualForm1(X))
  proj₁_²²̃₂::Hom(Form2(X)⊗DualForm2(X),Form2(X))
  proj₂_²²̃₂̃::Hom(Form2(X)⊗DualForm2(X),DualForm2(X))
  proj₁_²²̃₂⁺::Hom(Form2(X)⊕DualForm2(X),Form2(X))
  proj₂_²²̃₂̃⁺::Hom(Form2(X)⊕DualForm2(X),DualForm2(X))
  proj₁_⁰̃⁰₀̃::Hom(DualForm0(X)⊗Form0(X),DualForm0(X))
  proj₂_⁰̃⁰₀::Hom(DualForm0(X)⊗Form0(X),Form0(X))
  proj₁_⁰̃⁰₀̃⁺::Hom(DualForm0(X)⊕Form0(X),DualForm0(X))
  proj₂_⁰̃⁰₀⁺::Hom(DualForm0(X)⊕Form0(X),Form0(X))
  proj₁_⁰̃¹₀̃::Hom(DualForm0(X)⊗Form1(X),DualForm0(X))
  proj₂_⁰̃¹₁::Hom(DualForm0(X)⊗Form1(X),Form1(X))
  proj₁_⁰̃¹₀̃⁺::Hom(DualForm0(X)⊕Form1(X),DualForm0(X))
  proj₂_⁰̃¹₁⁺::Hom(DualForm0(X)⊕Form1(X),Form1(X))
  proj₁_⁰̃²₀̃::Hom(DualForm0(X)⊗Form2(X),DualForm0(X))
  proj₂_⁰̃²₂::Hom(DualForm0(X)⊗Form2(X),Form2(X))
  proj₁_⁰̃²₀̃⁺::Hom(DualForm0(X)⊕Form2(X),DualForm0(X))
  proj₂_⁰̃²₂⁺::Hom(DualForm0(X)⊕Form2(X),Form2(X))
  proj₁_⁰̃⁰̃₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X))
  proj₂_⁰̃⁰̃₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X))
  proj₁_⁰̃⁰̃₀̃⁺::Hom(DualForm0(X)⊕DualForm0(X),DualForm0(X))
  proj₂_⁰̃⁰̃₀̃⁺::Hom(DualForm0(X)⊕DualForm0(X),DualForm0(X))
  proj₁_⁰̃¹̃₀̃::Hom(DualForm0(X)⊗DualForm1(X),DualForm0(X))
  proj₂_⁰̃¹̃₁̃::Hom(DualForm0(X)⊗DualForm1(X),DualForm1(X))
  proj₁_⁰̃¹̃₀̃⁺::Hom(DualForm0(X)⊕DualForm1(X),DualForm0(X))
  proj₂_⁰̃¹̃₁̃⁺::Hom(DualForm0(X)⊕DualForm1(X),DualForm1(X))
  proj₁_⁰̃²̃₀̃::Hom(DualForm0(X)⊗DualForm2(X),DualForm0(X))
  proj₂_⁰̃²̃₂̃::Hom(DualForm0(X)⊗DualForm2(X),DualForm2(X))
  proj₁_⁰̃²̃₀̃⁺::Hom(DualForm0(X)⊕DualForm2(X),DualForm0(X))
  proj₂_⁰̃²̃₂̃⁺::Hom(DualForm0(X)⊕DualForm2(X),DualForm2(X))
  proj₁_¹̃⁰₁̃::Hom(DualForm1(X)⊗Form0(X),DualForm1(X))
  proj₂_¹̃⁰₀::Hom(DualForm1(X)⊗Form0(X),Form0(X))
  proj₁_¹̃⁰₁̃⁺::Hom(DualForm1(X)⊕Form0(X),DualForm1(X))
  proj₂_¹̃⁰₀⁺::Hom(DualForm1(X)⊕Form0(X),Form0(X))
  proj₁_¹̃¹₁̃::Hom(DualForm1(X)⊗Form1(X),DualForm1(X))
  proj₂_¹̃¹₁::Hom(DualForm1(X)⊗Form1(X),Form1(X))
  proj₁_¹̃¹₁̃⁺::Hom(DualForm1(X)⊕Form1(X),DualForm1(X))
  proj₂_¹̃¹₁⁺::Hom(DualForm1(X)⊕Form1(X),Form1(X))
  proj₁_¹̃²₁̃::Hom(DualForm1(X)⊗Form2(X),DualForm1(X))
  proj₂_¹̃²₂::Hom(DualForm1(X)⊗Form2(X),Form2(X))
  proj₁_¹̃²₁̃⁺::Hom(DualForm1(X)⊕Form2(X),DualForm1(X))
  proj₂_¹̃²₂⁺::Hom(DualForm1(X)⊕Form2(X),Form2(X))
  proj₁_¹̃⁰̃₁̃::Hom(DualForm1(X)⊗DualForm0(X),DualForm1(X))
  proj₂_¹̃⁰̃₀̃::Hom(DualForm1(X)⊗DualForm0(X),DualForm0(X))
  proj₁_¹̃⁰̃₁̃⁺::Hom(DualForm1(X)⊕DualForm0(X),DualForm1(X))
  proj₂_¹̃⁰̃₀̃⁺::Hom(DualForm1(X)⊕DualForm0(X),DualForm0(X))
  proj₁_¹̃¹̃₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X))
  proj₂_¹̃¹̃₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X))
  proj₁_¹̃¹̃₁̃⁺::Hom(DualForm1(X)⊕DualForm1(X),DualForm1(X))
  proj₂_¹̃¹̃₁̃⁺::Hom(DualForm1(X)⊕DualForm1(X),DualForm1(X))
  proj₁_¹̃²̃₁̃::Hom(DualForm1(X)⊗DualForm2(X),DualForm1(X))
  proj₂_¹̃²̃₂̃::Hom(DualForm1(X)⊗DualForm2(X),DualForm2(X))
  proj₁_¹̃²̃₁̃⁺::Hom(DualForm1(X)⊕DualForm2(X),DualForm1(X))
  proj₂_¹̃²̃₂̃⁺::Hom(DualForm1(X)⊕DualForm2(X),DualForm2(X))
  proj₁_²̃⁰₂̃::Hom(DualForm2(X)⊗Form0(X),DualForm2(X))
  proj₂_²̃⁰₀::Hom(DualForm2(X)⊗Form0(X),Form0(X))
  proj₁_²̃⁰₂̃⁺::Hom(DualForm2(X)⊕Form0(X),DualForm2(X))
  proj₂_²̃⁰₀⁺::Hom(DualForm2(X)⊕Form0(X),Form0(X))
  proj₁_²̃¹₂̃::Hom(DualForm2(X)⊗Form1(X),DualForm2(X))
  proj₂_²̃¹₁::Hom(DualForm2(X)⊗Form1(X),Form1(X))
  proj₁_²̃¹₂̃⁺::Hom(DualForm2(X)⊕Form1(X),DualForm2(X))
  proj₂_²̃¹₁⁺::Hom(DualForm2(X)⊕Form1(X),Form1(X))
  proj₁_²̃²₂̃::Hom(DualForm2(X)⊗Form2(X),DualForm2(X))
  proj₂_²̃²₂::Hom(DualForm2(X)⊗Form2(X),Form2(X))
  proj₁_²̃²₂̃⁺::Hom(DualForm2(X)⊕Form2(X),DualForm2(X))
  proj₂_²̃²₂⁺::Hom(DualForm2(X)⊕Form2(X),Form2(X))
  proj₁_²̃⁰̃₂̃::Hom(DualForm2(X)⊗DualForm0(X),DualForm2(X))
  proj₂_²̃⁰̃₀̃::Hom(DualForm2(X)⊗DualForm0(X),DualForm0(X))
  proj₁_²̃⁰̃₂̃⁺::Hom(DualForm2(X)⊕DualForm0(X),DualForm2(X))
  proj₂_²̃⁰̃₀̃⁺::Hom(DualForm2(X)⊕DualForm0(X),DualForm0(X))
  proj₁_²̃¹̃₂̃::Hom(DualForm2(X)⊗DualForm1(X),DualForm2(X))
  proj₂_²̃¹̃₁̃::Hom(DualForm2(X)⊗DualForm1(X),DualForm1(X))
  proj₁_²̃¹̃₂̃⁺::Hom(DualForm2(X)⊕DualForm1(X),DualForm2(X))
  proj₂_²̃¹̃₁̃⁺::Hom(DualForm2(X)⊕DualForm1(X),DualForm1(X))
  proj₁_²̃²̃₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X))
  proj₂_²̃²̃₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X))
  proj₁_²̃²̃₂̃⁺::Hom(DualForm2(X)⊕DualForm2(X),DualForm2(X))
  proj₂_²̃²̃₂̃⁺::Hom(DualForm2(X)⊕DualForm2(X),DualForm2(X))
  sum₀::Hom(Form0(X)⊗Form0(X),Form0(X))
  sum₁::Hom(Form1(X)⊗Form1(X),Form1(X))
  sum₂::Hom(Form2(X)⊗Form2(X),Form2(X))
  sum₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X))
  sum₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X))
  sum₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X))
end

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
#based off https://github.com/Roboneet/CartPole.jl/blob/master/src/Envs/classic_control/CartPole.jl
@with_kw struct CartPoleMDP <: MDP{Array{Float32}, Array{Float32}}
    gravity::Float32 = 98f-1
    masscart::Float32 = 1f0
    masspole::Float32 = 1f-1
    total_mass::Float32 = masspole + masscart
    length::Float32 = 5f-1 # actually half the pole's length
    polemass_length::Float32 = (masspole * length)
    force_mag::Float32 = 1f1
    τ::Float32 = 2f-2 # seconds between state updates
    maxT::Float32 = τ*100
    kinematics_integrator::AbstractString = "euler"

    # Angle at which to fail the episode
    θ_threshold_radians::Float32 = Float32(12 * 2 * π / 360)
    x_threshold::Float32 = 24f-1
    γ::Float32 = 0.99f0
end

function POMDPs.gen(env::CartPoleMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG)
    x, ẋ, θ, θ̇ , t = s

    failed =  !(all(vcat(-env.x_threshold .≤ x .≤ env.x_threshold,
            -env.θ_threshold_radians .≤ θ .≤ env.θ_threshold_radians)))
    timeout = t > env.maxT || t ≈ env.maxT

    if failed
        return (sp=fill(-100, size(s)), r=0)
    elseif timeout
        return (sp=fill(100, size(s)), r=0)
    end

     force = a .* env.force_mag  # action is +1 or -1
     cosθ = cos.(θ)
     sinθ = sin.(θ)
     temp = (force .+ env.polemass_length * θ̇  .^ 2 .* sinθ) / env.total_mass
     θacc = (env.gravity*sinθ .- cosθ.*temp) ./
            (env.length * (4f0/3 .- env.masspole * cosθ .^ 2 / env.total_mass))
     xacc  = temp .- env.polemass_length * θacc .* cosθ / env.total_mass
     if env.kinematics_integrator == "euler"
         x_ = x  .+ env.τ * ẋ
         ẋ_ = ẋ  .+ env.τ * xacc
         θ_ = θ  .+ env.τ * θ̇
         θ̇_ = θ̇  .+ env.τ * θacc
     else # semi-implicit euler
         ẋ_ = ẋ  .+ env.τ * xacc
         x_ = x  .+ env.τ * ẋ_
         θ̇_ = θ̇  .+ env.τ * θacc
         θ_ = θ  .+ env.τ * θ̇_
     end

     sp = vcat(x_, ẋ_, θ_, θ̇_, t+env.τ)
     (sp=sp, r=1f0)
end

function POMDPs.initialstate(mdp::CartPoleMDP)
    ImplicitDistribution((rng) -> Float32[(rand(Float32, 4) * 1f-1 .- 5f-2)..., 0f0])
end

isfailure(mdp::CartPoleMDP, s) = s[1] == -100
# Terminate if there is a failure or if we exceed the maximum time (in the case that we are tracking such things)
POMDPs.isterminal(mdp::CartPoleMDP, s) = isfailure(mdp, s) || s[1] == 100
POMDPs.discount(mdp::CartPoleMDP) = mdp.γ

ext: jl

max_forks_repo_name: ancorso/POMDPGym


text:
using OceanParameterizations
using WindMixing
using Flux
using OrdinaryDiffEq
using Plots
using ArgParse

"Returns a dictionary of command line arguments."
function parse_command_line_arguments()
    settings = ArgParseSettings()

    @add_arg_table! settings begin
        "--reconstruct_fluxes"
            help = ""
            default = false
            arg_type = Bool

        "--enforce_surface_fluxes"
            help = ""
            default = true
            arg_type = Bool

        "--train_test_same"
            help = ""
            default = false
            arg_type = Bool

        "--subsample_frequency"
            help = ""
            default = 1
            arg_type = Int
    end

    return parse_args(settings)
end

@info "Parsing command line arguments..."
args = parse_command_line_arguments()

reconstruct_fluxes = args["reconstruct_fluxes"]
println("Reconstruct fluxes? $(reconstruct_fluxes)")

enforce_surface_fluxes = args["enforce_surface_fluxes"]
println("Enforce surface fluxes? $(enforce_surface_fluxes)")

subsample_frequency = args["subsample_frequency"]
println("Subsample frequency for training... $(subsample_frequency)")

train_test_same = args["train_test_same"]
println("Train and test on the same file? $(train_test_same)")

file_labels = Dict(
    "free_convection" => "Free convection",
    "strong_wind" => "Strong wind",
    "strong_wind_no_coriolis" => "Strong wind, no rotation",
    "weak_wind_strong_cooling" => "Weak wind, strong cooling",
    "strong_wind_weak_cooling" => "Strong wind, weak cooling",
    "strong_wind_weak_heating" => "Strong wind, weak heating"
)

## Pick training and test simulations

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis",
            "weak_wind_strong_cooling", "strong_wind_weak_cooling", "strong_wind_weak_heating"]

for i=1:length(files)

    if train_test_same
        # Train on only file i
        train_files=files[i]
    else
        # Train on all except file i
        train_files = files[1:end .!= i]
    end

    𝒟train = WindMixing.data(train_files,
                        scale_type=ZeroMeanUnitVarianceScaling,
                        reconstruct_fluxes=reconstruct_fluxes,
                        subsample_frequency=subsample_frequency,
                        enforce_surface_fluxes=enforce_surface_fluxes)
    # Test on file i
    test_file = files[i]
    𝒟test = WindMixing.data(test_file,
                        override_scalings=𝒟train.scalings, # use the scalings from the training data
                        reconstruct_fluxes=reconstruct_fluxes,
                        subsample_frequency=subsample_frequency,
                        enforce_surface_fluxes=enforce_surface_fluxes)
    les = read_les_output(test_file)

    output_gif_directory="GP/subsample_$(subsample_frequency)/reconstruct_$(reconstruct_fluxes)/enforce_surface_fluxes_$(enforce_surface_fluxes)/train_test_same_$(train_test_same)/test_$(test_file)"
    directory = pwd() * "/$(output_gif_directory)/"
    mkpath(directory)
    file = directory*"_output.txt"
    touch(file)
    o = open(file, "w")

    write(o, "= = = = = = = = = = = = = = = = = = = = = = = = \n")
    write(o, "Test file: $(test_file) \n")
    write(o, "Output will be written to: $(output_gif_directory) \n")

    ## Gaussian Process Regression

    # A. Find the kernel that minimizes the prediction error on the training data
    # * Sweeps over length-scale hyperparameter value in logγ_range
    # * Sweeps over covariance functions
    # logγ_range=-1.0:0.5:1.0 # sweep over length-scale hyperparameter
    # uncomment the next three lines to try this but just for testing the GPR use the basic get_kernel stuff below
    # uw_kernel = best_kernel(𝒟train.uw, logγ_range=logγ_range)
    # vw_kernel = best_kernel(𝒟train.vw, logγ_range=logγ_range)
    # wT_kernel = best_kernel(𝒟train.wT, logγ_range=logγ_range)

    # OR set the kernel manually here (to save a bunch of time):
    # Result of the hyperparameter search - optimize_GP_kernels.jl
    if reconstruct_fluxes
        uw_kernel = get_kernel(2,0.4,0.0,euclidean_distance)
        vw_kernel = get_kernel(2,0.5,0.0,euclidean_distance)
        wT_kernel = get_kernel(2,1.3,0.0,euclidean_distance)
    else
        uw_kernel = get_kernel(2,0.4,0.0,euclidean_distance)
        vw_kernel = get_kernel(2,0.4,0.0,euclidean_distance)
        wT_kernel = get_kernel(2,1.2,0.0,euclidean_distance)
    end

    # Report the kernels and their properties
    write(o, "Kernel for u'w'..... $(uw_kernel) \n")
    write(o, "Kernel for v'w'..... $(vw_kernel) \n")
    write(o, "Kernel for w'T'..... $(wT_kernel) \n")

    # Trained GP models
    uw_GP_model = gp_model(𝒟train.uw, uw_kernel)
    vw_GP_model = gp_model(𝒟train.vw, vw_kernel)
    wT_GP_model = gp_model(𝒟train.wT, wT_kernel)

    # GP predictions on test data
    uw_GP = predict(𝒟test.uw, uw_GP_model)
    vw_GP = predict(𝒟test.vw, vw_GP_model)
    wT_GP = predict(𝒟test.wT, wT_GP_model)

    uw_GP = predict(𝒟test.uw, uw_GP_model; scaled = true)
    vw_GP = predict(𝒟test.vw, vw_GP_model; scaled = true)
    wT_GP = predict(𝒟test.wT, wT_GP_model; scaled = true)

    # Report GP prediction error on the fluxes
    write(o, "GP prediction error on u'w'..... $(mse(uw_GP)) \n")
    write(o, "GP prediction error on v'w'..... $(mse(vw_GP)) \n")
    write(o, "GP prediction error on w'T'..... $(mse(wT_GP)) \n")

    # Compare GP predictions to truth
    myanimate(xs, name) = animate_prediction(xs, name, 𝒟test, test_file;
                            legend_labels=["GP(u,v,T)","Truth"], directory=output_gif_directory)
    myanimate(uw_GP, "uw")
    myanimate(vw_GP, "vw")
    myanimate(wT_GP, "wT")

    t  = 𝒟test.t
    Nz = 32
    f  = les.f⁰
    H  = Float32(abs(𝒟test.uw.z[end] - 𝒟test.uw.z[1]))
    τ  = Float32(abs(t[:,1][end] - t[:,1][1]))
    u_scaling = 𝒟test.scalings["u"]
    v_scaling = 𝒟test.scalings["v"]
    T_scaling = 𝒟test.scalings["T"]
    uw_scaling = 𝒟test.scalings["uw"]
    vw_scaling = 𝒟test.scalings["vw"]
    wT_scaling = 𝒟test.scalings["wT"]

    get_μ_σ(name) = (𝒟test.scalings[name].μ, 𝒟test.scalings[name].σ)
    μ_u, σ_u = get_μ_σ("u")
    μ_v, σ_v = get_μ_σ("v")
    μ_T, σ_T = get_μ_σ("T")
    μ_uw, σ_uw = get_μ_σ("uw")
    μ_vw, σ_vw = get_μ_σ("vw")
    μ_wT, σ_wT = get_μ_σ("wT")
    D_cell = Float32.(Dᶜ(Nz, 1/Nz))

    top_bottom(x) = (Float32(x[1,1]), Float32(x[end,1]))
    uw_top, uw_bottom = top_bottom(𝒟test.uw.scaled)
    vw_top, vw_bottom = top_bottom(𝒟test.vw.scaled)
    wT_top, wT_bottom = top_bottom(𝒟test.wT.scaled)

    # enforce surface fluxes in the predictions
    function predict_and_enforce_fluxes(model, x, top, bottom)
        return [top; model(x); bottom]
    end

    A = - τ / H
    B = f * τ

    function NDE_nondimensional_flux(x, p, t)
        u = x[1:Nz]
        v = x[Nz+1:2*Nz]
        T = x[2*Nz+1:96]
        dx₁ = A .* σ_uw ./ σ_u .* D_cell * uw_GP_model(x) .+ B ./ σ_u .* (σ_v .* v .+ μ_v) #nondimensional gradient
        dx₂ = A .* σ_vw ./ σ_v .* D_cell * vw_GP_model(x) .- B ./ σ_v .* (σ_u .* u .+ μ_u)
        dx₃ = A .* σ_wT ./ σ_T .* D_cell * wT_GP_model(x)
        return [dx₁; dx₂; dx₃]
    end

    function time_window(t, uvT, trange)
        return (Float32.(t[trange]), Float32.(uvT[:,trange]))
    end

    timesteps = 1:1:length(t) #1:5:100
    uvT₀ = Float32.(𝒟test.uvT_scaled[:,1])

    t_train, uvT_train = time_window(t, 𝒟test.uvT_scaled, timesteps)
    t_train = Float32.(t_train ./ τ)

    prob = ODEProblem(NDE_nondimensional_flux, uvT₀, (t_train[1], t_train[end]), saveat=t_train)
    sol = solve(prob, Tsit5())

    split_array(uvT) = (uvT[1:Nz,:], uvT[Nz+1:2*Nz,:], uvT[2*Nz+1:end,:])
    u_pred, v_pred, T_pred = split_array(sol)

    u_pair = (u_pred, 𝒟test.u.scaled)
    v_pair = (v_pred, 𝒟test.v.scaled)
    T_pair = (T_pred, 𝒟test.T.scaled)

    myanimate(xs, name) = animate_prediction(xs, name, 𝒟test, test_file;
                                legend_labels=["GP(u,v,T)","Truth"], directory=output_gif_directory)
    myanimate(u_pair, "u")
    myanimate(v_pair, "v")
    myanimate(T_pair, "T")

    write(o, "GP prediction error on u........ $(mse(u_pair)) \n")
    write(o, "GP prediction error on v........ $(mse(v_pair)) \n")
    write(o, "GP prediction error on T........ $(mse(T_pair)) \n")

    # Close output file
    close(o)
end

# tpoint = 100
# split_vector(uvT) = (uvT[1:Nz], uvT[Nz+1:2*Nz], uvT[2*Nz+1:end])
# u_pred, v_pred, T_pred = split_vector(sol[:,tpoint])
# p1 = plot(u_pred, zC_coarse, label="GP_DE")
# plot!(𝒟test.u.scaled[:,tpoint], zC_coarse, label="truth")
# p2 = plot(v_pred, zC_coarse, label="GP_DE")
# plot!(𝒟test.v.scaled[:,tpoint], zC_coarse, label="truth")
# p3 = plot(T_pred, zC_coarse, label="GP_DE")
# plot!(𝒟test.T.scaled[:,tpoint], zC_coarse, label="truth")
# layout = @layout [a b c]
# p = plot(p1, p2, p3, layout=layout)
# png(p, "hello.png")


text:
#md# It is the code to draw this repo's logo
using WordCloud
Luxor = WordCloud.Render.Luxor

function cloudshape(height, args...; backgroundcolor=(0, 0, 0, 0))
    height = ceil(height)
    r = height / 2
    d = Luxor.Drawing((2height, height)..., :svg)
    Luxor.origin()
    Luxor.background(parsecolor(backgroundcolor))
    Luxor.setcolor(parsecolor((0.22, 0.596, 0.149)))
    Luxor.pie(0, 0, r, 0, 2π, :fill)
    Luxor.setcolor(parsecolor((0.584, 0.345, 0.698)))
    Luxor.pie(r, r, r, -π, 0, :fill)
    Luxor.setcolor(parsecolor((0.796, 0.235, 0.20)))
    Luxor.Luxor.pie(-r, r, r, -π, 0, :fill)
    Luxor.finish()
    d
end
wc = wordcloud(
    repeat(string.('a':'z'), 5),
    repeat([1], 26 * 5),
    mask=cloudshape(500),
    transparent=(0, 0, 0, 0),
    colors=1,
    angles=-90:0,
    font="JuliaMono Black",
    density=0.7,
    )
generate!(wc, 2000, retry=1)
recolor!(wc, style=:main)
println("results are saved to logo.svg")
paint(wc, "logo.svg", background=false)
wc
#eval# runexample(:logo)
#md# ![](logo.svg)


text:
"""
    Keypoints{N}(sz) <: Block

A block representing an array of size `sz` filled with keypoints of type
`SVector{N}`.
"""
struct Keypoints{N,M} <: Block
    sz::NTuple{M,Int}
end
Keypoints{N}(n::Int) where {N} = Keypoints{N,1}((n,))
Keypoints{N}(t::NTuple{M,Int}) where {N,M} = Keypoints{N,M}(t)

function checkblock(
    ::Keypoints{N,M},
    ::AbstractArray{<:Union{SVector{N,T},Nothing},M},
) where {M,N,T}
    return true
end

mockblock(block::Keypoints{N}) where {N} = rand(SVector{N,Float32}, block.sz)

# ## Visualization

function showblock!(io, ::ShowText, block::Keypoints{2}, obs)
    print(io, UnicodePlots.scatterplot(first.(obs), last.(obs), marker=:cross))
end

function showblock!(io, ::ShowText, block::Bounded{2, <:Keypoints{2}}, obs)
    h, w = block.size
    plot = UnicodePlots.scatterplot(
        first.(obs), last.(obs),
        xlim=(0, w), ylim=(0, h), marker=:cross)
    print(io, plot)
end


text:
module TimesTables

using QML
using Observables
using Pkg.Artifacts
using Pkg.TOML
using Dates
using Logging
import FileIO

# Sound stuff
# import PortAudio
# import LibSndFile

export TimesOp, Multiplication, Division, Addition, Subtraction, OpGenerator, DefaultOpGenerator
export compute, check, generate

const artifact_version = v"0.0.1"

assetroot() = joinpath(artifact"timestables-assets", "timestables-assets-$(artifact_version)")
qmlfile() = joinpath(assetroot(), "qml", "times_problems.qml")
soundsdir() = joinpath(assetroot(), "sounds")

"""
Base type for all operationd, e.g. multiplication and division
"""
abstract type TimesOp end

"""
Encapsulate a × b
"""
struct Multiplication <: TimesOp
  a::Int32
  b::Int32
end

"""
Encapsulate a ÷ b
"""
struct Division <: TimesOp
  a::Int
  b::Int

  function Division(a,b)
    if a % b != 0
      error("Division constructor: $a is not divisible by $b")
    end
    return new(a,b)
  end
end

"""
Encapsulate a + b
"""
struct Addition <: TimesOp
  a::Int
  b::Int
end

"""
Encapsulate a - b
"""
struct Subtraction <: TimesOp
  a::Int
  b::Int
end

"""
Compute the result of an op
"""
compute(op::Multiplication) = op.a * op.b
compute(op::Division) = op.a ÷ op.b
compute(op::Addition) = op.a + op.b
compute(op::Subtraction) = op.a - op.b

"""
Check the result of an op
"""
check(op::TimesOp, x::Number) = (x == compute(op))

"""
Base type for op generators
"""
abstract type OpGenerator end

"""
Default OpGenerator that generates any type of op, limiting divisions and multiplication operands
between `times_min` and `times_max` and the sum of operands to addition and subtraction to addition_max.
These defaults correspond to second grade of primary school in Belgium.
"""
struct DefaultOpGenerator <: OpGenerator
  times_min::Int
  times_max::Int
  subtraction_min::Int
  addition_max::Int
end

DefaultOpGenerator() = DefaultOpGenerator(1,10,0,100)

"""
Generate a random multiplication op
"""
generate(::Type{Multiplication}, g::DefaultOpGenerator) = Multiplication(rand(g.times_min:g.times_max), rand(g.times_min:g.times_max))

"""
Generate a random division op
"""
function generate(::Type{Division}, g::DefaultOpGenerator)
  a = rand(g.times_min:g.times_max)
  b = rand(g.times_min:g.times_max)
  return Division(a*b, b)
end

"""
Generate a random addition op
"""
function generate(::Type{Addition}, g::DefaultOpGenerator)
  a = rand(0:g.addition_max)
  b = rand(0:g.addition_max - a)
  return Addition(a, b)
end

"""
Generate a random Subtraction op
"""
function generate(::Type{Subtraction}, g::DefaultOpGenerator)
  a = rand(0:g.addition_max)
  b = rand(0:a-g.subtraction_min)
  return Subtraction(a, b)
end

function generate_type(g::DefaultOpGenerator)
  return rand([Multiplication, Division, Addition, Subtraction])
end

"""
Generate a random op
"""
generate(g::OpGenerator = DefaultOpGenerator()) = generate(generate_type(g), g)

Base.print(io::IO, o::Multiplication) = write(io, "$(o.a) × $(o.b)")
Base.print(io::IO, o::Division) = write(io, "$(o.a) : $(o.b)")
Base.print(io::IO, o::Addition) = write(io, "$(o.a) + $(o.b)")
Base.print(io::IO, o::Subtraction) = write(io, "$(o.a) - $(o.b)")

function getconfig()
  configfile = joinpath(first(DEPOT_PATH), "prefs", "TimesTables.toml")
  if !isfile(configfile)
    mkpath(dirname(configfile))
    defaultconfig = Dict(
      "timesmin" => 1,
      "timesmax" => 10,
      "subtractionmin" => 0,
      "additionmax" => 100,
      "maxnumcorrect" => 3,
      "uselogfile" => false,
      "lockedwindow" => false
    )
    open(configfile, write=true) do f
      TOML.print(f, defaultconfig)
    end
  end
  return TOML.parsefile(configfile)
end

const generator = Ref{OpGenerator}()
const currentop = Observable{TimesOp}()
const currentopstring = Observable("")
const answerstring = Observable("")
const numcorrect = Observable(0)
const maxnumcorrect = Observable(0)
const state = Observable("STARTUP")

on(currentop) do op
  currentopstring[] = string(op)
end

on(answerstring) do s
  if isempty(s)
    return
  end
  answer = parse(Int, s)
  if check(currentop[], answer)
    numcorrect[] += 1
    @info "Correct answer for $(currentop[]): $s"
    if numcorrect[] == maxnumcorrect[]
      state[] = "FINISHED";
      playsound(cheer[])
    else
      state[] = "CORRECT";
      currentop[] = generate(generator[])
      playsound(woohoo[])
    end
  else
    numcorrect[] = 0;
    @info "Wrong answer for $(currentop[]): $s"
    state[] = "ERROR";
    playsound(ohno[])
  end
end

const cheer = Ref{Any}()
const ohno = Ref{Any}()
const woohoo = Ref{Any}()

#const audio_out = Ref{PortAudio.PortAudioStream}()
playsound(s) = nothing #PortAudio.write(audio_out[], s)

function julia_main()
  config = getconfig()

  cheer[] = nothing # FileIO.load(joinpath(soundsdir(), "cheer.wav"))
  ohno[] = nothing # FileIO.load(joinpath(soundsdir(), "ohno.wav"))
  woohoo[] = nothing # FileIO.load(joinpath(soundsdir(), "woohoo.wav"))

  #audio_out[] = PortAudio.PortAudioStream(0, 2)

  if config["uselogfile"]
    logio = io = open(joinpath(first(DEPOT_PATH), "logs", "TimesTables.log"), "a+")
    global_logger(SimpleLogger(logio))
  end

  generator[] = DefaultOpGenerator(config["timesmin"], config["timesmax"], config["subtractionmin"], config["additionmax"])
  numcorrect[] = 0
  maxnumcorrect[] = config["maxnumcorrect"]
  currentop[] = generate(generator[])
  state[] = "STARTUP"
  load(qmlfile(),
    φ = Float64(Base.MathConstants.golden),
    lockedwindow = config["lockedwindow"],
    problem = JuliaPropertyMap(
      "question" => currentopstring,
      "answer" => answerstring,
      "state" => state,
      "numcorrect" => numcorrect,
      "maxnumcorrect" => maxnumcorrect,
    ))

  exec()

  if numcorrect[] != maxnumcorrect[]
    @warn "Premature exit on $(now())"
  else
    @info "Normal exit on $(now())"
  end

  if config["uselogfile"]
    close(logio)
  end

  return 0
end

end # module


text:
using LinearAlgebra, Distributions

# stepsize: the stepsize
# n_steps: number of steps for the leapfrog
# mm: metric matrix
# post: the posterior distribution. Should be a function!
# grad: the gradient. Should be a function!
# mom: the moment for the beggining of the leapfrog
# par: the paramater for the beggining of the leapfrog

function sympl_integrator(stepsize,n_steps,mm,post,grad,mom,par)
    new_mom = mom
    new_par = par
    direction = 1
    for l in 1:n_steps
        new_mom = new_mom + 1/2*stepsize*grad(new_par)
        new_par = new_par + stepsize*mm*new_mom
        if isnan(post(new_par)) || isinf(post(new_par))
            @warn "Invalid value of the distribution. Going back"
            direction = -1
        end
        new_mom = new_mom + 1/2*stepsize*grad(new_par)
    end
    return new_par,new_mom
end

function sympl_integrator_window(stepsize,n_steps,mm,mom_dist,post,grad,mom,par)
    unif = Uniform(0,1)
    new_mom = mom
    new_par = par
    prob = 0
    prob_cur = 0
    cur_par, cur_mom = new_par,new_mom
    direction = 1
    for l in 1:n_steps
        new_mom = new_mom + 1/2*stepsize*grad(new_par)
        new_par = new_par + stepsize*mm*new_mom
        if isnan(loglike(new_par)) || isinf(loglike(new_par))
            direction = -1
        end
        new_mom = new_mom + 1/2*stepsize*grad(new_par)
        prob = exp(log_pdf(mom_dist,new_mom) + post(new_par))
        prob_cur += prob
        if prob/prob_cur > rand(unif)
            cur_par,cur_mom = new_par,new_mom
        end
    end
    return new_par,new_mom, cur_par,cur_mom
end

function metropolis_step(par_new,mom_new,par_old,mom_old,post,dist_mom)
    metropolis = post(par_new) + logpdf(dist_mom,mom_new) - post(par_old) -logpdf(dist_mom,mom_old)
    metropolis = min(exp(metropolis),1)
    unif = Uniform(0,1)
    if metropolis > rand(unif)
        return par_new, mom_new, 1
    else
        return par_old, mom_old, 0
    end
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using JsonGrinder, JSON, Test, SparseArrays, Flux, Random, HierarchicalUtils
using JsonGrinder: ExtractScalar, ExtractCategorical, ExtractArray, ExtractDict, ExtractVector
using JsonGrinder: extractempty
using Mill
using Mill: catobs, nobs, MaybeHotMatrix
using LinearAlgebra

function misequals(a, b)
	for i in 1:length(a)
		ismissing(b[i]) && !ismissing(a[i]) && return(false)
		ismissing(a[i]) && !ismissing(b[i]) && return(false)
		!ismissing(a[i]) && a[i] != b[i] && return(false)
	end
	return(true)
end

@testset "Testing ExtractScalar" begin
	sc = ExtractScalar(Float64,2,3)
	@test all(sc("5").data .== [9])
	@test all(sc(5).data .== [9])
	@test all(sc(nothing).data .=== [missing])
	@test all(sc(missing).data .=== [missing])
	@test nobs(sc(missing)) == 1
	@test nobs(sc(nothing)) == 1
	@test sc(extractempty).data isa Matrix{Float64}
	@test nobs(sc(extractempty)) == 0
	@test nobs(sc(5)) == 1

	sc = ExtractScalar(Float32, 0.5, 4.0)
	@test sc(1).data isa Matrix{Float32}
	@test sc(extractempty).data isa Matrix{Float32}
end

@testset "ExtractCategorical" begin
	e = ExtractCategorical(["a","b"])
	@test e("a").data ≈ [1, 0, 0]
	@test e("b").data ≈ [0, 1, 0]
	@test e("z").data ≈ [0, 0, 1]
	@test isequal(e(nothing).data, [missing missing missing]')
	@test isequal(e(missing).data, [missing missing missing]')
	@test typeof(e("a").data) == MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}
	@test typeof(e(nothing).data) == MaybeHotMatrix{Missing,Array{Missing,1},Int64,Missing}
	@test typeof(e(missing).data) == MaybeHotMatrix{Missing,Array{Missing,1},Int64,Missing}
	@test e(extractempty).data isa MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}
	@test nobs(e(extractempty)) == 0

	@test e(["a", "b"]).data ≈ [1 0; 0 1; 0 0]
	@test isequal(e(["a", missing]).data, [true missing; false missing; false missing])
	@test isequal(e(["a", missing, "x"]).data, [true missing false; false missing false; false missing true])
	@test typeof(e(["a", "b"]).data) == MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}
	@test typeof(e(["a", "b", nothing]).data) == MaybeHotMatrix{Union{Missing, Int64},Array{Union{Missing, Int64},1},Int64,Union{Missing, Bool}}

	@test isnothing(ExtractCategorical([]))
	e2 = ExtractCategorical(JsonGrinder.Entry(Dict("a"=>1,"c"=>1), 2))
	@test e2("a").data ≈ [1, 0, 0]
	@test e2("c").data ≈ [0, 1, 0]
	@test e2("b").data ≈ [0, 0, 1]
	@test isequal(e2(nothing).data, [missing missing missing]')
	@test isequal(e2(missing).data, [missing missing missing]')

	@test catobs(e("a"), e("b")).data ≈ [1 0; 0 1; 0 0]
	@test reduce(catobs, [e("a").data, e("b").data]) ≈ [1 0; 0 1; 0 0]
	# this is not working now, due to type stability issues
	# @test catobs(e("a").data, e("b").data) ≈ [1 0; 0 1; 0 0]
	@test isequal(e(Dict(1=>2)).data |> collect, [missing missing missing]')

	@test nobs(e("a")) == 1
	@test nobs(e("b")) == 1
	@test nobs(e("z")) == 1
	@test nobs(e(nothing)) == 1
	@test nobs(e(missing)) == 1
	@test nobs(e(missing).data) == 1
	@test nobs(e([missing, nothing])) == 2
	@test nobs(e([missing, nothing, "a"])) == 3

	e3 = ExtractCategorical(JsonGrinder.Entry(Dict(1=>1,2=>1), 2))
	@test e3(1).data ≈ [1, 0, 0]
	@test e3(2).data ≈ [0, 1, 0]
	@test e3(4).data ≈ [0, 0, 1]
	@test e3(1.).data ≈ [1, 0, 0]
	@test e3(2.).data ≈ [0, 1, 0]
	@test e3(4.).data ≈ [0, 0, 1]
	@test isequal(e3([]).data, [missing missing missing]')

	e4 = ExtractCategorical(JsonGrinder.Entry(Dict(1.0=>1,2.0=>1), 2))
	@test e4(1).data ≈ [1, 0, 0]
	@test e4(2).data ≈ [0, 1, 0]
	@test e4(4).data ≈ [0, 0, 1]
	@test e4(1.).data ≈ [1, 0, 0]
	@test e4(2.).data ≈ [0, 1, 0]
	@test e4(4.).data ≈ [0, 0, 1]
	@test isequal(e4([]).data, [missing missing missing]')
end

@testset "Testing array conversion" begin
	Mill.emptyismissing!(false)
	sc = ExtractArray(ExtractCategorical(2:4))
	@test all(sc([2,3,4]).data.data .== Matrix(1.0I, 4, 3))
	@test nobs(sc(nothing).data) == 0
	@test sc(nothing).data.data isa MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}
	@test nobs(sc(nothing).data.data) == 0
	@test all(sc(nothing).bags.bags .== [0:-1])

	Mill.emptyismissing!(true)
	@test sc(nothing).data isa Missing
	@test all(sc(nothing).bags.bags .== [0:-1])
	Mill.emptyismissing!(false)

	@test nobs(sc(extractempty).data.data) == 0
	@test nobs(sc(extractempty).data) == 0
	@test isempty(sc(extractempty).bags.bags)
	@test sc(extractempty).data.data isa MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}

	sc = ExtractArray(ExtractScalar(Float32))
	@test all(sc([2,3,4]).data.data .== [2 3 4])
	@test nobs(sc(nothing).data) == 0
	@test all(sc(nothing).bags.bags .== [0:-1])

	@test nobs(sc(extractempty).data.data) == 0
	@test nobs(sc(extractempty).data) == 0
	@test isempty(sc(extractempty).bags.bags)
	@test sc(extractempty).data.data isa Matrix{Float32}

	@test sc(nothing).data.data isa Matrix{Float32}
end

@testset "Testing feature vector conversion" begin
	sc = ExtractVector(5)
	@test sc([1, 2, 2, 3, 4]).data isa Matrix
	@test all(sc([1, 2, 2, 3, 4]).data .== [1, 2, 2, 3, 4])
	@test sc([1, 2, 2, 3, 4]).data isa Array{Float32, 2}
	@test sc(extractempty).data isa Matrix{Float32}
	@test nobs(sc(extractempty).data) == 0

	sc = ExtractVector{Int64}(5)
	@test all(sc([1, 2, 2, 3, 4]).data .== [1, 2, 2, 3, 4])
	@test sc([1, 2, 2, 3, 4]).data isa Array{Int64, 2}
	@test sc(nothing).data isa Matrix
	@test all(sc(nothing).data .=== missing)
	@test sc(extractempty).data isa Matrix{Int64}
	@test nobs(sc(extractempty).data) == 0

	# feature vector longer than expected
	sc = ExtractVector{Float32}(5)
	@test all(sc([1, 2, 2, 3, 4, 5]).data .== [1, 2, 2, 3, 4])
	@test typeof(sc([1, 2, 3, 4, 5]).data) == Array{Float32,2}
	@test sc([5, 6]).data[1:2] ≈ [5, 6]
	@test typeof(sc([1, 2]).data) == Array{Union{Missing,Float32},2}
	@test sc([1, 2]).data isa Matrix
	@test all(sc([5, 6]).data[3:5] .=== missing)
	@test all(sc(Dict(1=>2)).data .=== missing)
end

@testset "Testing ExtractDict" begin
	dict = Dict("a" => ExtractScalar(Float64,2,3),
				"b" => ExtractScalar(Float64),
				"c" => ExtractArray(ExtractScalar(Float64,2,3)))
	br = ExtractDict(dict)
	a1 = br(Dict("a" => 5, "b" => 7, "c" => [1,2,3,4]))
	a2 = br(Dict("a" => 5, "b" => 7))
	a3 = br(Dict("a" => 5, "c" => [1,2,3,4]))

	@test catobs(a1,a1)[:a].data ≈ [9 9]
	@test catobs(a1,a1)[:b].data ≈ [7 7]
	@test catobs(a1,a1)[:c].data.data ≈  [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a1,a1)[:c].bags .== [1:4,5:8])

	@test catobs(a1,a2)[:a].data ≈ [9 9]
	@test catobs(a1,a2)[:b].data ≈ [7 7]
	@test catobs(a1,a2)[:c].data.data ≈ [-3 0 3 6]
	@test all(catobs(a1,a2)[:c].bags .== [1:4,0:-1])

	@test catobs(a2,a3)[:a].data ≈ [9 9]
	@test misequals(catobs(a2,a3)[:b].data, [7.0 missing])
	@test catobs(a2,a3)[:c].data.data ≈ [-3 0 3 6]
	@test all(catobs(a2,a3)[:c].bags .== [0:-1,1:4])

	@test catobs(a1,a3)[:a].data ≈ [9 9]
	@test misequals(catobs(a1,a3)[:b].data, [7.0 missing])
	@test catobs(a1,a3)[:c].data.data ≈ [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a1,a3)[:c].bags .== [1:4,5:8])

	@test a1[:a].data ≈ [9]
	@test a1[:b].data ≈ [7]
	@test a2[:a].data ≈ [9]
	@test a2[:b].data ≈ [7]
	@test a3[:a].data ≈ [9]
	@test misequals(a3[:b].data, [missing])

	@test a1[:c].data.data ≈ [-3 0 3 6]
	@test all(a1[:c].bags .== [1:4])

	@test a3[:c].data.data ≈ [-3 0 3 6]
	@test all(a3[:c].bags .== [1:4])
	@test catobs(a3,a3)[:c].data.data ≈ [-3 0 3 6 -3 0 3 6]
	@test all(catobs(a3,a3)[:c].bags .== [1:4,5:8])

	a4 = br(extractempty)
	@test nobs(a4[:a]) == 0
	@test a4[:a].data isa Matrix{Float64}
	@test nobs(a4[:b]) == 0
	@test a4[:b].data isa Matrix{Float64}
	@test nobs(a4[:c]) == 0
	@test nobs(a4[:c].data) == 0
	@test a4[:c].data.data isa Matrix{Float64}
end

@testset "Testing Nested Missing Arrays" begin
	dict = Dict("a" => ExtractArray(ExtractScalar(Float32,2,3)),
		"b" => ExtractArray(ExtractScalar(Float32,2,3)))
	br = ExtractDict(dict)
	a1 = br(Dict("a" => [1,2,3], "b" => [1,2,3,4]))
	a2 = br(Dict("b" => [2,3,4]))
	a3 = br(Dict("a" => [2,3,4]))
	a4 = br(Dict{String,Any}())

	@test all(catobs(a1,a2).data[1].data.data .== [-3.0  0.0  3.0  6.0  0.0  3.0  6.0])
	@test all(catobs(a1,a2).data[1].bags .== [1:4, 5:7])
	@test all(catobs(a1,a2).data[2].data.data .== [-3.0  0.0  3.0])
	@test all(catobs(a1,a2).data[2].bags .== [1:3, 0:-1])

	@test all(catobs(a2,a3).data[1].data.data .== [0.0  3.0  6.0])
	@test all(catobs(a2,a3).data[1].bags .== [1:3, 0:-1])
	@test all(catobs(a2,a3).data[2].data.data .== [0 3 6])
	@test all(catobs(a2,a3).data[2].bags .== [0:-1, 1:3])

	@test all(catobs(a1,a4).data[1].data.data .== [-3.0  0.0  3.0  6.0])
	@test all(catobs(a1,a4).data[1].bags .== [1:4, 0:-1])
	@test all(catobs(a1,a4).data[2].data.data .== [-3.0  0.0  3.0])
	@test all(catobs(a1,a4).data[2].bags .== [1:3, 0:-1])

	@test all(a4.data[2].data.data isa Array{Float32,2})
	a4 = br(extractempty)
	@test nobs(a4[:a]) == 0
	@test nobs(a4[:a].data) == 0
	@test a4[:a].data.data isa Matrix{Float32}
	@test nobs(a4[:b]) == 0
	@test nobs(a4[:b].data) == 0
	@test a4[:b].data.data isa Matrix{Float32}
end

@testset "ExtractOneHot" begin
	samples = ["{\"name\": \"a\", \"count\" : 1}",
		"{\"name\": \"b\", \"count\" : 2}",]
	vs = JSON.parse.(samples)

	e = ExtractOneHot(["a","b"], "name", "count")
	@test e(vs).data[:] ≈ [1, 2, 0]
	@test e(nothing).data[:] ≈ [0, 0, 0]
	@test e(missing).data[:] ≈ [0, 0, 0]
	@test e(vs).data isa SparseMatrixCSC{Float32,Int64}
	@test e(nothing).data isa SparseMatrixCSC{Float32,Int64}
	@test e(missing).data isa SparseMatrixCSC{Float32,Int64}
	@test e(extractempty).data isa SparseMatrixCSC{Float32,Int64}
	@test nobs(e(extractempty)) == 0
	@test nobs(e(extractempty).data) == 0

	e = ExtractOneHot(["a","b"], "name", nothing)
	@test e(vs).data[:] ≈ [1, 1, 0]
	@test e(nothing).data[:] ≈ [0, 0, 0]
	@test e(missing).data[:] ≈ [0, 0, 0]
	@test typeof(e(vs).data) == SparseMatrixCSC{Float32,Int64}
	@test typeof(e(nothing).data) == SparseMatrixCSC{Float32,Int64}
	@test typeof(e(missing).data) == SparseMatrixCSC{Float32,Int64}
	vs = JSON.parse.(["{\"name\": \"c\", \"count\" : 1}"])
	@test e(vs).data[:] ≈ [0, 0, 1]
	@test typeof(e(vs).data) == SparseMatrixCSC{Float32,Int64}
	@test e(extractempty).data isa SparseMatrixCSC{Float32,Int64}
	@test nobs(e(extractempty)) == 0
	@test nobs(e(extractempty).data) == 0
end

@testset "ExtractString" begin
	e = ExtractString()
	@test e("Hello").data.s == ["Hello"]
	@test e(Symbol("Hello")).data.s == ["Hello"]
	@test e(["Hello", "world"]).data.s == ["Hello", "world"]
	@test e("Hello") isa ArrayNode{NGramMatrix{String,Vector{String},Int64},Nothing}
	@test e("Hello").data isa NGramMatrix{String,Vector{String},Int64}
	@test all(e(missing).data.s .=== [missing])
	@test e(missing).data isa NGramMatrix{Missing,Vector{Missing},Missing}
end

@testset "Extractor of keys as field" begin
	JsonGrinder.updatemaxkeys!(1000)
	js = [Dict(randstring(5) => rand()) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))

	b = ext(js[1])
	k = only(keys(js[1]))
	@test b.data[:item].data[1] ≈ (js[1][k] - ext.item.c) * ext.item.s
	@test b.data[:key].data.s[1] == k

	b = ext(nothing)
	@test nobs(b) == 1
	@test nobs(b.data) == 0
	b = ext(Dict())
	@test nobs(b) == 1
	@test nobs(b.data) == 0

	b = ext(extractempty)
	@test nobs(b) == 0
	@test nobs(b.data) == 0
	@test nobs(b.data[:item]) == 0
	@test b.data[:item].data isa Matrix{Float32}
	@test nobs(b.data[:key]) == 0
	@test b.data[:key].data isa Mill.NGramMatrix{String,Array{String,1},Int64}

	js = [Dict(randstring(5) => Dict(:a => rand(), :b => randstring(1))) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))

	b = ext(js[1])
	k = only(keys(js[1]))
	i = ext.item(js[1][k])
	@test b.data[:item][:a].data == i[:a].data
	@test b.data[:item][:b].data ==i[:b].data
	@test b.data[:key].data.s[1] == k

	b = ext(nothing)
	@test nobs(b) == 1
	@test nobs(b.data) == 0
	b = ext(Dict())
	@test nobs(b) == 1
	@test nobs(b.data) == 0
end

@testset "equals and hash test" begin
	other1 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
		"d" => ExtractVector(4),
	)
	br1 = ExtractDict(other1)
	other11 = Dict(
		"a" => ExtractArray(br1),
		"b" => ExtractScalar(Float64,2,3),
	)
	br11 = ExtractDict(other11)

	other2 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
		"d" => ExtractVector(4),
	)
	br2 = ExtractDict(other2)
	other22 = Dict("a" => ExtractArray(br2), "b" => ExtractScalar(Float64,2,3))
	br22 = ExtractDict(other22)

	@test hash(br11) === hash(br22)
	@test hash(br11) !== hash(br1)
	@test hash(br11) !== hash(br2)
	@test hash(br1) === hash(br2)

	@test br11 == br22
	@test br11 != br1
	@test br11 != br2
	@test br1 == br2
end
# todo: add tests for settings of suggestextractor!!!
# todo: add tests for single keyed dict, that model generated from such jsons behaves as we want
@testset "Extractor skip empty lists" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}], "b": [], "c": [[], []], "d": 1, "e": {"a": 1}, "f": {"a": []}}""")
	j2 = JSON.parse("""{"a": [{"a":3},{"b":4}], "b": [], "c": [[]], "d": 1, "e": {"a": 2}, "f": {"a": []}}""")
	j3 = JSON.parse("""{"a": [{"a":1},{"b":3}], "b": [], "c": [[], [], [], []], "d": 2, "e": {"a": 3}, "f": {"a": []}}""")
	j4 = JSON.parse("""{"a": [{"a":2},{"b":4}], "b": [], "c": [], "d": 4, "e": {"a": 3}, "f": {"a": []}}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch)
	@test ext[:a] isa ExtractArray
	@test isnothing(ext[:b])
	@test isnothing(ext[:c])
	@test ext[:d] isa ExtractScalar
	@test ext[:e] isa ExtractDict
	@test isnothing(ext[:f])
end

@testset "Extractor of numbers as strings" begin
	j1 = JSON.parse("""{"a": "1", "b": "a", "c": "1.1", "d": 1.1, "e": "1.2", "f": 1}""")
	j2 = JSON.parse("""{"a": "2", "b": "b", "c": "2", "d": 2, "e": "1.3", "f": 2}""")
	j3 = JSON.parse("""{"a": "3", "b": "c", "c": "2.3", "d": 2.3, "e": "1.4", "f": 3}""")
	j4 = JSON.parse("""{"a": "4", "b": "c", "c": "5", "d": 5, "e": "1.4", "f": 3}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch)

	@test ext[:a] isa ExtractCategorical
	@test ext[:b] isa ExtractString
	@test ext[:c] isa ExtractCategorical
	@test ext[:d] isa ExtractScalar{Float32}
	@test ext[:e] isa ExtractCategorical
	@test ext[:f] isa ExtractScalar{Float32}

	ext_j1 = ext(j1)
	ext_j2 = ext(j2)
	ext_j3 = ext(j3)
	ext_j4 = ext(j4)

	@test eltype(ext_j1[:a].data) <: Bool
	@test eltype(ext_j1[:b].data) <: Int64
	@test eltype(ext_j1[:c].data) <: Bool
	@test eltype(ext_j1[:d].data) <: Float32
	@test eltype(ext_j1[:e].data) <: Bool
	@test eltype(ext_j1[:f].data) <: Float32

	@test eltype(ext_j2[:a].data) <: Bool
	@test eltype(ext_j2[:b].data) <: Int64
	@test eltype(ext_j2[:c].data) <: Bool
	@test eltype(ext_j2[:d].data) <: Float32
	@test eltype(ext_j2[:e].data) <: Bool
	@test eltype(ext_j2[:f].data) <: Float32

	@test eltype(ext_j3[:a].data) <: Bool
	@test eltype(ext_j3[:b].data) <: Int64
	@test eltype(ext_j3[:c].data) <: Bool
	@test eltype(ext_j3[:d].data) <: Float32
	@test eltype(ext_j3[:e].data) <: Bool
	@test eltype(ext_j3[:f].data) <: Float32

	@test eltype(ext_j4[:a].data) <: Bool
	@test eltype(ext_j4[:b].data) <: Int64
	@test eltype(ext_j4[:c].data) <: Bool
	@test eltype(ext_j4[:d].data) <: Float32
	@test eltype(ext_j4[:e].data) <: Bool
	@test eltype(ext_j4[:f].data) <: Float32

	@test ext_j1["U"].data ≈ [0]
	@test ext_j2["U"].data ≈ [3/13]
	@test ext_j3["U"].data ≈ [4/13]
	@test ext_j4["U"].data ≈ [1]

	m = reflectinmodel(sch, ext)
	# this is for stable version
	# @test buf_printtree(m) == """
	# ProductModel ↦ ArrayModel(Dense(42, 10))
	#   ├── a: ArrayModel(Dense(5, 10))
	#   ├── b: ArrayModel(Dense(2053, 10))
	#   ├── c: ArrayModel(Dense(5, 10))
	#   ├── d: ArrayModel(identity)
	#   ├── e: ArrayModel(Dense(4, 10))
	#   └── f: ArrayModel(identity)"""
	  # this is test for master
	@test buf_printtree(m) == """
	ProductModel … ↦ ArrayModel(Dense(42, 10))
	  ├── a: ArrayModel(Dense(5, 10))
	  ├── b: ArrayModel(Dense(2053, 10))
	  ├── c: ArrayModel(Dense(5, 10))
	  ├── d: ArrayModel(identity)
	  ├── e: ArrayModel(Dense(4, 10))
	  └── f: ArrayModel(identity)"""
end
# todo: add separate tests for reflectinmodel(sch, ext) to test behavior with various missing and non-missing stuff in schema
# e.g. empty string, missing keys, irregular schema and MultiRepresentation

@testset "Suggest feature vector extraction" begin
	j1 = JSON.parse("""{"a": "1", "b": [1, 2, 3], "c": [1, 2, 3]}""")
	j2 = JSON.parse("""{"a": "2", "b": [2, 2, 3], "c": [1, 2, 3, 4]}""")
	j3 = JSON.parse("""{"a": "3", "b": [3, 2, 3], "c": [1, 2, 3, 4, 5]}""")
	j4 = JSON.parse("""{"a": "4", "b": [2, 3, 4], "c": [1, 2, 3]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch)

	@test ext[:a] isa ExtractCategorical
	@test ext[:b] isa ExtractVector
	@test ext[:b].n == 3
	@test ext[:c] isa ExtractArray{ExtractScalar{Float32}}

	ext_j1 = ext(j1)
	ext_j2 = ext(j2)
	ext_j3 = ext(j3)
	ext_j4 = ext(j4)

	@test ext_j1["U"].data ≈ [1, 2, 3]
	@test ext_j2["U"].data ≈ [2, 2, 3]
	@test ext_j3["U"].data ≈ [3, 2, 3]
	@test ext_j4["U"].data ≈ [2, 3, 4]

	@test ext_j1["s"].data ≈ [0 .25 .5]
	@test ext_j2["s"].data ≈ [0 .25 .5 .75]
	@test ext_j3["s"].data ≈ [0 .25 .5 .75 1.]
	@test ext_j4["s"].data ≈ [0 .25 .5]
end

@testset "Suggest complex" begin
	JsonGrinder.updatemaxkeys!(1000)
	js = [Dict("a" => rand(), "b" => Dict(randstring(5) => rand()), "c"=>[rand(), rand()], "d"=>[rand() for i in 1:rand(1:10)]) for _ in 1:1000]
	sch = JsonGrinder.schema(js)
	ext = JsonGrinder.suggestextractor(sch, (;key_as_field = 500))

	@test ext[:a] isa ExtractScalar{Float32}
	@test ext[:b] isa JsonGrinder.ExtractKeyAsField
	@test ext[:b].key isa ExtractString
	@test ext[:b].item isa ExtractScalar{Float32}
	@test ext[:c] isa ExtractVector
	@test ext[:c].n == 2
	@test ext[:d] isa ExtractArray
end

@testset "Extractor typed missing bags" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [1, 2]}""")
	j3 = JSON.parse("""{"a": [1, 2, 3]}""")
	j4 = JSON.parse("""{"a": [1, 2, 3, 4]}""")
	j5 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	ext = suggestextractor(sch)

	ext_j1 = ext(j1)
	ext_j2 = ext(j2)
	ext_j3 = ext(j3)
	ext_j4 = ext(j4)

	@test ext_j1[:a].data.data isa Array{Float32,2}
	@test ext_j2[:a].data.data isa Array{Float32,2}
	@test ext_j3[:a].data.data isa Array{Float32,2}
	@test ext_j4[:a].data.data isa Array{Float32,2}
end

@testset "testing irregular extractor" begin
	j1 = JSON.parse("""{"a": 4}""")
	j2 = JSON.parse("""{"a": { "a": "hello", "b":[5,6]}}""")
	j3 = JSON.parse("""{"a": [1, 2, 3, 4]}""")
	j4 = JSON.parse("""{"a": 5}""")
	j5 = JSON.parse("""{"a": 3}""")

	sch = schema([j1,j2,j3,j4,j5])
	ext = suggestextractor(sch)
	a = ext(j1)
	@test a[:a][:e1].data[1] == 0.5
	@test ismissing(a[:a][:e2].data[:a].data.s[1])
	@test nobs(a[:a][:e3]) == 1
	@test nobs(a[:a][:e3].data) == 1
end

@testset "Mixed scalar extraction" begin
	j1 = JSON.parse("""{"a": "1"}""")
	j2 = JSON.parse("""{"a": 4}""")
	j3 = JSON.parse("""{"a": "3.1"}""")
	j4 = JSON.parse("""{"a": 2.5}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4])
	sch_hash = hash(sch)
	ext = suggestextractor(sch)
	@test sch_hash == hash(sch)	# testing that my entry parkour does not modify schema
	@test nchildren(ext[:a]) == 1

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	@test e1["k"].data ≈ [0]
	@test e2["k"].data ≈ [1]
	@test e3["k"].data ≈ [0.7]
	@test e4["k"].data ≈ [0.5]
end

@testset "Mixed scalar extraction with other types" begin
	j1 = JSON.parse("""{"a": "1"}""")
	j2 = JSON.parse("""{"a": 2.5}""")
	j3 = JSON.parse("""{"a": "3.1"}""")
	j4 = JSON.parse("""{"a": 5}""")
	j5 = JSON.parse("""{"a": "4.5"}""")
	j6 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")
	j7 = JSON.parse("""{"a": {"Sylvanas is the worst warchief ever": "yes"}}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5, j6, j7])
	ext = suggestextractor(sch)

	@test buf_printtree(sch) ==
    """
    [Dict] (updated = 7)
      └── a: [MultiEntry] (updated = 7)
               ├── 1: [Scalar - String], 3 unique values, updated = 3
               ├── 2: [Scalar - Float64,Int64], 2 unique values, updated = 2
               ├── 3: [List] (updated = 1)
               │        └── [Scalar - Int64], 5 unique values, updated = 5
               └── 4: [Dict] (updated = 1)
                        └── Sylvanas is the worst warchief ever: [Scalar - String], 1 unique values, updated = 1"""

	@test buf_printtree(ext) ==
    """
	Dict
	  └── a: MultiRepresentation
	           ├── e1: FeatureVector with 5 items
	           ├── e2: Dict
	           │         └── Sylvanas is the worst warchief ever: String
	           └── e3: Float32"""

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	e5 = ext(j5)
	@test e1["s"].data ≈ [0]
	@test e2["s"].data ≈ [0.375]
	@test e3["s"].data ≈ [0.525]
	@test e4["s"].data ≈ [1.0]
	@test e5["s"].data ≈ [0.875]
	@test buf_printtree(e1) ==
	"""
	ProductNode with 1 obs
	  └── a: ProductNode with 1 obs
	           ├── e1: ArrayNode(5×1 Array with Missing elements) with 1 obs
	           ├── e2: ProductNode with 1 obs
	           │         └── Sylvanas is the worst warchief ever: ArrayNode(2053×1 NGramMatrix with Missing elements) with 1 obs
	           └── e3: ArrayNode(1×1 Array with Float32 elements) with 1 obs"""
end

@testset "mixing numeric and non-numeric strings" begin
	j1 = JSON.parse("""{"a": "hello"}""")
	j2 = JSON.parse("""{"a": "4"}""")
	j3 = JSON.parse("""{"a": 5}""")
	j4 = JSON.parse("""{"a": 3.0}""")
	j5 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5])
	ext = suggestextractor(sch)

	@test buf_printtree(sch) ==
	"""
	[Dict] (updated = 5)
	  └── a: [MultiEntry] (updated = 5)
	           ├── 1: [Scalar - String], 2 unique values, updated = 2
	           ├── 2: [Scalar - Float64,Int64], 2 unique values, updated = 2
	           └── 3: [List] (updated = 1)
	                    └── [Scalar - Int64], 5 unique values, updated = 5"""

	@test buf_printtree(ext) ==
	"""
	Dict
	  └── a: MultiRepresentation
	           ├── e1: String
	           ├── e2: Float32
	           └── e3: FeatureVector with 5 items"""

	e1 = ext(j1)
	e2 = ext(j2)
	e3 = ext(j3)
	e4 = ext(j4)
	e5 = ext(j5)
	@test all(e1["k"].data .=== [missing])
	@test e2["k"].data ≈ [0.5]
	@test e3["k"].data ≈ [1.0]
	@test e4["k"].data ≈ [0.0]
	@test all(e5["k"].data .=== [missing])

	@test hash(ext) !== hash(suggestextractor(JsonGrinder.schema([j1, j2, j4, j5])))
end

@testset "empty string and substring" begin
	a, b, c = split("a b ", " ")
	d = "d"
	e = ""
	f = 5.2
	@test a isa SubString{String}
	@test b isa SubString{String}
	@test c isa SubString{String}
	@test d isa String
	@test e isa String
	@test f isa AbstractFloat
	@test c == ""
	ext = JsonGrinder.extractscalar(AbstractString)
	@test SparseMatrixCSC(ext(c).data) == SparseMatrixCSC(ext(e).data)
	@test all(ext(f).data.s .=== ext(nothing).data.s)
end

@testset "key as field" begin
	j1 = JSON.parse("""{"a": 1}""")
	j2 = JSON.parse("""{"b": 2.5}""")
	j3 = JSON.parse("""{"c": 3.1}""")
	j4 = JSON.parse("""{"d": 5}""")
	j5 = JSON.parse("""{"e": 4.5}""")
	j6 = JSON.parse("""{"f": 5}""")
	sch = JsonGrinder.schema([j1, j2, j3, j4, j5, j6])
	e = JsonGrinder.key_as_field(sch, NamedTuple(), path = "")
	@test e isa JsonGrinder.ExtractKeyAsField

	e2 = JsonGrinder.key_as_field(sch, NamedTuple(), path = "")
	@test hash(e) === hash(e2)
	@test e == e2
end

@testset "AuxiliaryExtractor HUtils" begin
	e2 = ExtractCategorical(["a","b"])
	e = AuxiliaryExtractor(e2, (ext, sample)->ext(String(sample)))

	@test e("b") == e(:b)
	@test e("b").data ≈ [0, 1, 0]
	@test e(:b).data ≈ [0, 1, 0]

    @test buf_printtree(e) ==
	"""
    Auxiliary extractor with
      └── Categorical d = 3"""
end

@testset "Skipping single dict key" begin
	j1 = JSON.parse("""{"a": []}""")
	j2 = JSON.parse("""{"a": [1, 2]}""")
	j3 = JSON.parse("""{"a": [1, 2, 3]}""")
	j4 = JSON.parse("""{"a": [1, 2, 3, 4]}""")
	j5 = JSON.parse("""{"a": [1, 2, 3, 4, 5]}""")

	sch = JsonGrinder.schema([j1, j2, j3, j4, j5])
	ext = suggestextractor(sch)

	@test buf_printtree(ext) ==
	"""
	Dict
	  └── a: Array of
	           └── Float32"""

	ext_j2 = ext(j2)
	@test buf_printtree(ext_j2) ==
    """
	ProductNode with 1 obs
	  └── a: BagNode with 1 obs
	           └── ArrayNode(1×2 Array with Float32 elements) with 2 obs"""
		   # todo: add more tests for integration with Mill to make sure it's propagated well
end


text:
module ManualMemory

macro splice(iterator, body)
  @assert iterator.head == :call
  @assert iterator.args[1] == :in
  Expr(:..., :(($(esc(body)) for $(esc(iterator.args[2])) in $(esc(iterator.args[3])))))
end

include("manuals.jl")
include("manual_vectors.jl")
include("manual_bit_vectors.jl")
include("manual_strings.jl")
include("manual_alloc.jl")

export Manual, ManualVector, ManualBitVector, ManualString, @a, @v

end


text:
using BenchmarkTools
using Plots
using Base.Threads

f(x, p) = p*x*(one(x)-x)

function calc_attractor!(out, f, p, num_attract, warmup)

    x = 0.25
    for i in 1:warmup
        x = f(x, p)
    end
    @inbounds out[1] = x

    @inbounds for i in 1:num_attract-1
        out[i+1] = f(out[i], p)
    end
end

x₀ = 0.25
out = Vector{typeof(x₀)}(undef, 150)
out[1] = x₀
@btime calc_attractor!(out, f, 2.9, 150, 400)

@code_llvm calc_attractor!(out, f, 2.9, 150, 400)
@code_warntype calc_attractor!(out, f, 2.9, 150, 400)

function bifurcation_data!(data, f, r, n_cols)

    Threads.@threads for i in 1:n_cols
        calc_attractor!(@view(data[:, i]), f, r[i], 150, 400)
    end
end

r = 2.9:0.001:4
n_cols = length(r)
n_rows = 150
pl_data = Array{Float64}(undef, n_rows, n_cols)
pl_data[1, :] .= 0.25
@btime bifurcation_data!(pl_data, f, r, n_cols)

p = plot(xlims=(2.9,4), ylims=(0, 1), title="Bifurcation Plot", label=false)
for i in 1:n_rows
    scatter!(r, pl_data[i, :], label=false, color=:black, alpha=0.1)
end
p


text:
#
# Advent of Code
# Day 10: Balance Bots
#

println("Advent of Code")
println("Day 10: Balance Bots")
println("")

DEBUG = false

# Data file
# Input data has been saved to this file
inputDataFileName = "day10-input.txt"

inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

#inputData ="""
#value 5 goes to bot 2
#bot 2 gives low to bot 1 and high to bot 0
#value 3 goes to bot 1
#bot 1 gives low to output 1 and high to bot 0
#bot 0 gives low to output 2 and high to output 0
#value 2 goes to bot 2
#"""
#TARGET = [2, 5]

TARGET = [17, 61]

inputDataLines = split(strip(inputData), '\n')

NONE = 0
BOT = 1
OUTPUT = 2

DESTINATION = Dict(["bot" => BOT, "output" => OUTPUT])

type Bot
  inventory::Vector{Int}
  lowType::Int
  lowDestination::Int
  highType::Int
  highDestination::Int
end

bots = Dict{Int,Bot}()
outputs = Dict{Int, Vector{Int}}()

function runBot(bot::Int)
  #
  # Do the things the bot should do
  #

  # only do stuff if we have two items

  if haskey(bots, bot)
    myBot = bots[bot]
    if length(myBot.inventory) == 2
      if DEBUG println("\nBot $bot is processing its instructions.") end

      (lowValue, highValue) = sort(myBot.inventory)

      if [lowValue, highValue] == TARGET
        println("\n\nBot $bot is responsible for comparing ", lowValue, " to ", highValue, ".")
        #exit()
      end

      # process the values
      process(lowValue, myBot.lowType, myBot.lowDestination)
      process(highValue, myBot.highType, myBot.highDestination)

      # get rid of the tokens
      shift!(myBot.inventory)
      shift!(myBot.inventory)
    end

    if DEBUG println("") end

  else
    println("Tried to process bot $bot but it doesn't exist")
  end
end

function process(value::Int, myType::Int, destination::Int)
  if myType == BOT
    if DEBUG println("... giving $value to bot $destination") end

    if haskey(bots, destination)
      push!(bots[destination].inventory, value)
    else
      bots[destination] = Bot(Vector{Int}([value]), NONE, 0, NONE, 0)
    end
    runBot(destination)

  elseif myType == OUTPUT
    if DEBUG println("... putting $value in output $destination") end
    if haskey(outputs, destination)
      push!(outputs[destination], value)
    else
      outputs[destination] = Vector{Int}([value])
    end
  else
    println("Unknown type $myType")
  end
end

for line in inputDataLines
  # parse the instruction
  # use regex

  m = match(r"(value|bot) (\d+) (?:goes|gives low) to (bot|output) (\d+)(?: and high to (bot|output) (\d+))?", line)

  if m != nothing
    if m.captures[1] == "value"
      value = parse(m.captures[2])
      bot = parse(m.captures[4])
      if haskey(bots, bot)
        push!(bots[bot].inventory, value)
      else
        bots[bot] = Bot(Vector{Int}([value]), NONE, 0, NONE, 0)
      end

      if DEBUG println("Bot $bot was given $value.") end

      #runBot(bot)

    elseif m.captures[1] == "bot"
      bot = parse(m.captures[2])
      lowType = DESTINATION[m.captures[3]]
      lowDestination = parse(m.captures[4])
      highType = DESTINATION[m.captures[5]]
      highDestination = parse(m.captures[6])

      if haskey(bots, bot)
        inventory = bots[bot].inventory
      else
        inventory = Vector{Int}()
      end

      bots[bot] = Bot(inventory, lowType, lowDestination, highType, highDestination)

      if DEBUG println("Bot $bot was given new instructions.") end

      #runBot(bot)
    else
      println("Unknown directive: ", m.captures[1])
    end
  else
    println("No regex match to: $line")
  end
end

for botID in keys(bots)
  runBot(botID)
end

result = outputs[0][1] * outputs[1][1] * outputs[2][1]

println("Multiplication of 0, 1, 2: $result")


text:
module Clingo
include("LibClingo/LibClingo.jl")
using .LibClingo
include("handle.jl")
include("control.jl")
end


text:
using JuliaBlackBoard

mean(x) = sum(x)/length(x)

OPEN = JuliaBlackBoard.OPEN(@__FILE__)

OPEN() do io

    # -- Question ---
    # use a NUM type with answer computed
    q = mt"""
# The mean

Find the mean of

```
2 3 4
```
"""

    question(io, NUM,q,mean((2,3,4)))

    # -- Question ---
    # using some latex to test tth
    q = mt"""
# The mean

Find the mean, $(x_1 + \cdots x_n)/n$, of

```
{{{:xs}}}
```
"""

    for xs in ((2,3,4),
               (3,4,5))
        question(io, NUM, q(xs = join(xs, ",")), mean(xs))
    end

        # -- Question ---
    # using some latex to test markdown -> LaTeX
    q = mdltx"""
# The mean

Find the mean, $(x_1 + \cdots x_n)/n$, of

```
{{{:xs}}}
```
"""

    for xs in ((2,3,4),
               (3,4,5))
        question(io, NUM, q(xs = join(xs, ",")), mean(xs))
    end

    # -- Question ---
    # using all latex to test LaTeX and tectonic installation

     q = ltx"""
\section{The mean}

Find the mean, $\bar{x}$, of

\begin{verbatim}
2 3 4
\end{verbatim}
"""

    question(io, NUM,  q, mean((2,3,4)))

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: identitystore
using AWS.Compat
using AWS.UUIDs

"""
    describe_group(group_id, identity_store_id)
    describe_group(group_id, identity_store_id, params::Dict{String,<:Any})

Retrieves the group metadata and attributes from GroupId in an identity store.

# Arguments
- `group_id`: The identifier for a group in the identity store.
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string which contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

"""
describe_group(GroupId, IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("DescribeGroup", Dict{String, Any}("GroupId"=>GroupId, "IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
describe_group(GroupId, IdentityStoreId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("DescribeGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GroupId"=>GroupId, "IdentityStoreId"=>IdentityStoreId), params)); aws_config=aws_config)

"""
    describe_user(identity_store_id, user_id)
    describe_user(identity_store_id, user_id, params::Dict{String,<:Any})

Retrieves the user metadata and attributes from UserId in an identity store.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string which contains number and lower case letters. This value is generated at the time
  that a new identity store is created.
- `user_id`: The identifier for a user in the identity store.

"""
describe_user(IdentityStoreId, UserId; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("DescribeUser", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId, "UserId"=>UserId); aws_config=aws_config)
describe_user(IdentityStoreId, UserId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("DescribeUser", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId, "UserId"=>UserId), params)); aws_config=aws_config)

"""
    list_groups(identity_store_id)
    list_groups(identity_store_id, params::Dict{String,<:Any})

Lists the attribute name and value of the group that you specified in the search. We only
support DisplayName as a valid filter attribute path currently, and filter is required.
This API returns minimum attributes, including GroupId and group DisplayName in the
response.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string which contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: A list of Filter objects, which is used in the ListUsers and ListGroups
  request.
- `"MaxResults"`: The maximum number of results to be returned per request, which is used
  in the ListUsers and ListGroups request to specify how many results to return in one page.
  The length limit is 50 characters.
- `"NextToken"`: The pagination token used for the ListUsers and ListGroups APIs. This
  value is generated by the identity store service and is returned in the API response if the
  total results are more than the size of one page, and when this token is used in the API
  request to search for the next page.
"""
list_groups(IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("ListGroups", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
list_groups(IdentityStoreId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("ListGroups", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId), params)); aws_config=aws_config)

"""
    list_users(identity_store_id)
    list_users(identity_store_id, params::Dict{String,<:Any})

Lists the attribute name and value of the user that you specified in the search. We only
support UserName as a valid filter attribute path currently, and filter is required. This
API returns minimum attributes, including UserId and UserName in the response.

# Arguments
- `identity_store_id`: The globally unique identifier for the identity store, such as
  d-1234567890. In this example, d- is a fixed prefix, and 1234567890 is a randomly generated
  string which contains number and lower case letters. This value is generated at the time
  that a new identity store is created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: A list of Filter objects, which is used in the ListUsers and ListGroups
  request.
- `"MaxResults"`: The maximum number of results to be returned per request, which is used
  in the ListUsers and ListGroups request to specify how many results to return in one page.
  The length limit is 50 characters.
- `"NextToken"`: The pagination token used for the ListUsers and ListGroups APIs. This
  value is generated by the identity store service and is returned in the API response if the
  total results are more than the size of one page, and when this token is used in the API
  request to search for the next page.
"""
list_users(IdentityStoreId; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("ListUsers", Dict{String, Any}("IdentityStoreId"=>IdentityStoreId); aws_config=aws_config)
list_users(IdentityStoreId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = identitystore("ListUsers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("IdentityStoreId"=>IdentityStoreId), params)); aws_config=aws_config)


text:
# Step 2

D = zeros(Float64, size(C))
for i in eachindex(D)
	if C[i] > 0
		if L[C[i]] == "F"
			D[i] = 8.0
		end
		if L[C[i]] == "W"
			D[i] = 0.75
		end
		if L[C[i]] == "N"
			D[i] = 3.0
		end
		if L[C[i]] == "U"
			D[i] = 0.001
		end
	end
end

V = zeros(Int64, size(D))
for re in 1:50
	#x, y = sample(1:size(V,1)), sample(1:size(V,2))
	x, y, = 1, 1
	xD = copy(D)
	xV = zeros(Int64, size(D))
	for i in 1:300000
		xV[x,y] = xV[x,y]+1
		xD[x,y] = xD[x,y]*0.5
		xf, xt = x-1, x+1
		yf, yt = y-1, y+1
		xand = max(xf,1):min(xt,size(V,1))
		yand = max(yf,1):min(yt,size(V,2))
		tD = xD[xand,yand].+rand()*0.01
		keep = sample(eachindex(tD), weights(vec(tD)))
		ni = ind2sub(tD, keep)
		x, y = xand[ni[1]], yand[ni[2]]
		#=if ((1000 < x < 1200) & (1650 < y < 1850))
			V = V .+ xV
			info("#$re\thit after $i")
			break
		end=#
	end
	V = V .+ xV
end

V = V ./maximum(V);

bg = readpng("./img/test.png")
w = bg.width
h = bg.height

begin
	Drawing(w, h, "sim.png")
	origin()
	placeimage(bg, -w/2, -h/2)
	for i in 1:size(D,1)
		for j in 1:size(D,2)
			xinit = (i-1)
			yinit = (j-1)
			op = V[i,j] == 0.0 ? 0.0 : 0.9
			setopacity(op)
			sethue(V[i,j], 0.5, 0.7)
			rect(yinit-w/2, xinit-h/2, 1, 1, :fill)
		end
	end
	# NOTE the position is wrong
	setopacity(1.0)
	sethue("red")
	setline(4)
	rect(165*gsize-w/2, 80*gsize-h/2, 20*gsize, 20*gsize, :stroke)
	finish()
	preview()
end

max_stars_repo_name: tpoisot/WaterFinder

max_issues_repo_name: tpoisot/WaterFinder

max_forks_repo_name: tpoisot/WaterFinder


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"Quantum machine learning"
module QML

include("Optimization.jl")
include("CircuitOptimization.jl")

using QuantumCircuits.QCircuits.Gates: ParameterT
using QuantumCircuits.QCircuits.Circuit: getparameters, getRandParameters
using QuantumCircuits.QCircuits.QBase: setparameters!, bindparameters!
using QuantumCircuits.QML.CircuitOptimization

export ParameterT, getparameters, getRandParameters, setparameters!,
       bindparameters!, findparam

end # module


text:
__precompile__()
module RCall
using Compat, DataFrames, DataArrays
import Compat.String

import DataStructures: OrderedDict

import Base: eltype, show, convert, isascii,
    length, size, getindex, setindex!, start, next, done, names

export RObject,
   Sxp, NilSxp, StrSxp, CharSxp, LglSxp, IntSxp, RealSxp, CplxSxp,
   ListSxp, VecSxp, EnvSxp, LangSxp, ClosSxp, S4Sxp,
   getattrib, setattrib!, getnames, setnames!, getclass, setclass!, attributes,
   globalEnv,
   isna, anyna, isnull,
   rcopy, rparse, rprint, reval, rcall, rlang,
   @rimport, @rlibrary, @rput, @rget, @var_str, @R_str

include("setup.jl")
include("types.jl")
include("constants.jl")
include("methods.jl")
include("convert-base.jl")
include("convert-data.jl")
include("convert-default.jl")
include("eventloop.jl")
include("eval.jl")
include("io.jl")
include("functions.jl")
include("callback.jl")
include("operators.jl")
include("library.jl")
include("IJulia.jl")
include("render.jl")
include("macros.jl")
include("repl.jl")
include("deprecated.jl")

end # module


text:
"""

Automatically generated from SageMath calculations

Fergus Baker - 9th Nov 2021
             - 10th Feb 2022: updated to include Jacobian method
             - 1st Apr 2022: fix sign error in mu
"""
module MorrisThorneCoords

using ..ComputedGeodesicEquations

@inline function geodesic_eq(u, v, b)
    @fastmath ComputedGeodesicEquations.@let_unpack u v begin
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        out1 = 0
        out2 = r * v_phi^2 * sin_theta + r * v_theta^2
        out3 =
            1 / 2 * ((b^2 + r^2) * v_phi^2 * cos_theta - 4 * r * v_r * v_theta) /
            (b^2 + r^2)
        out4 =
            -((b^2 + r^2) * v_phi * v_theta * cos_theta + 2 * r * v_phi * v_r * sin_theta) /
            ((b^2 + r^2) * sin_theta)

        (out1, out2, out3, out4)
    end
end

@inline function constrain(μ, u, v, b)
    @fastmath ComputedGeodesicEquations.@let_unpack u v begin
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        sqrt(
            b^2 * v_phi^2 * sin_theta +
            r^2 * v_phi^2 * sin_theta +
            b^2 * v_theta^2 +
            r^2 * v_theta^2 +
            μ^2 +
            v_r^2,
        )
    end
end

@inline function jacobian(u, b)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        comp1 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 0
        ]
        comp2 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 2*r 0
            0 0 0 2*r*sin_theta
        ]
        comp3 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 (b^2+r^2)*cos_theta
        ]
        comp4 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 0
        ]
        (comp1, comp2, comp3, comp4)
    end
end

@inline function metric(u, b)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        ComputedGeodesicEquations.@SMatrix [
            -1 0 0 0
            0 1 0 0
            0 0 b^2+r^2 0
            0 0 0 (b^2+r^2)*sin_theta
        ]
    end
end

@inline function inverse_metric(u, b)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        ComputedGeodesicEquations.@SMatrix [
            -1 0 0 0
            0 1 0 0
            0 0 1/(b^2+r^2) 0
            0 0 0 1/((b^2+r^2)*sin_theta)
        ]
    end
end

end # module

@with_kw struct MorrisThorne{T} <: AbstractMetricParams{T}
    @deftype T
    b = 1.0
end

@with_kw struct MorrisThorneJac{T} <: AbstractMetricParams{T}
    @deftype T
    b = 1.0
end

geodesic_eq(m::MorrisThorne{T}, u, v) where {T} = MorrisThorneCoords.geodesic_eq(u, v, m.b)
geodesic_eq(m::MorrisThorneJac{T}, u, v) where {T} = jac_geodesic_eq(m, u, v)

constrain(m::MorrisThorne{T}, u, v; μ::T = 0.0) where {T} =
    MorrisThorneCoords.constrain(μ, u, v, m.b)
constrain(m::MorrisThorneJac{T}, u, v; μ::T = 0.0) where {T} =
    MorrisThorneCoords.constrain(μ, u, v, m.b)

# specialisations
metric(m::MorrisThorne{T}, u) where {T} = MorrisThorneCoords.metric(u, m.b)
inverse_metric(m::MorrisThorne{T}, u) where {T} = MorrisThorneCoords.inverse_metric(u, m.b)
jacobian(m::MorrisThorne{T}, u) where {T} = MorrisThorneCoords.jacobian(u, m.b)
metric(m::MorrisThorneJac{T}, u) where {T} = MorrisThorneCoords.metric(u, m.b)
inverse_metric(m::MorrisThorneJac{T}, u) where {T} =
    MorrisThorneCoords.inverse_metric(u, m.b)
jacobian(m::MorrisThorneJac{T}, u) where {T} = MorrisThorneCoords.jacobian(u, m.b)

export MorrisThorneCoords, MorrisThorne, MorrisThorneJac


text:
using Gadfly

function fh(x, y)::Real
    return 1
end

function dc(p)
    return -1 + abs(4 - sqrt(sum(p .^ 2)))
end

function dca(p)
    return sqrt(sum(p .^ 2)) - 2
end

function dcb(p, x1, x2, y1, y2)
    return -min(min(min(-y1 + p[2], y2 - p[2], -x1 + p[1], x2 - p[1])))
end

function dcrc(p, xc, yc, r)
    return sqrt((p[1] - xc) .^ 2 + (p[2] - yc) .^ 2) - r
end

function dcrcc(p)
    return max(dcb(p, -100, 100, -100, 100), -dcrc(p, 20, 10, 15))
end

function eplot(x,y)
    Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))
end

x, y = generate(dcrcc, fh, [-100 -100; 100 100], 7.0, [-100 -100; -100 100; 100 -100; 100 100])

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
function compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k, cold_start_s)
	sum_error = 0.; ndcg_sum = 0.;
	for i in collect(cold_start_s)
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		score = zeros(len)
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			score[j] = dot(ui,vj)
		end
		error_this = 0
		n_comps_this = 0
		for j in 1:(len - 1)
			jval = vals_d2_bar[j]
			for k in (j + 1):len
				kval = vals_d2_bar[k]
				if score[j] >= score[k] && jval < kval
					error_this += 1
				end
				if score[j] <= score[k] && jval > kval
					error_this += 1
				end
				n_comps_this += 1
			end
		end
		sum_error += error_this / n_comps_this

		p1 = sortperm(score, rev = true)
		p1 = p1[1:ndcg_k]
		M1 = vals_d2_bar[p1]
		p2 = sortperm(vals_d2_bar, rev = true)
		p2 = p2[1:ndcg_k]
		M2 = vals_d2_bar[p2]
		dcg = 0.; dcg_max = 0.
		for k = 1:ndcg_k
			dcg += (2 ^ M1[k] - 1) / log2(k + 1)
			dcg_max += (2 ^ M2[k] - 1) / log2(k + 1)
		end
		ndcg_sum += dcg / dcg_max
	end
	return sum_error / length(cold_start_s), ndcg_sum / length(cold_start_s)
end

function compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)
	sum_error = 0.; ndcg_sum = 0.;
	for i = 1:d1
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		score = zeros(len)
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			score[j] = dot(ui,vj)
		end
		error_this = 0
		n_comps_this = 0
		for j in 1:(len - 1)
			jval = vals_d2_bar[j]
			for k in (j + 1):len
				kval = vals_d2_bar[k]
				if score[j] >= score[k] && jval < kval
					error_this += 1
				end
				if score[j] <= score[k] && jval > kval
					error_this += 1
				end
				n_comps_this += 1
			end
		end
		sum_error += error_this / n_comps_this

		p1 = sortperm(score, rev = true)
		p1 = p1[1:ndcg_k]
		M1 = vals_d2_bar[p1]
		p2 = sortperm(vals_d2_bar, rev = true)
		p2 = p2[1:ndcg_k]
		M2 = vals_d2_bar[p2]
		dcg = 0.; dcg_max = 0.
		for k = 1:ndcg_k
			dcg += (2 ^ M1[k] - 1) / log2(k + 1)
			dcg_max += (2 ^ M2[k] - 1) / log2(k + 1)
		end
		ndcg_sum += dcg / dcg_max
	end
	return sum_error / d1, ndcg_sum / d1
end

function compute_precision(U, V, X, Y, d1, d2, rows, vals, rows_t, vals_t, zero_items)
	K = [1, 5, 10] # K has to be increasing order
	precision = [0, 0, 0]
	for i = 1:d1
    #for i = shuffle(1:d1)[1:500]
		tmp = nzrange(Y, i)
		test = Set{Int64}()
		for j in tmp
			push!(test, rows_t[j])
		end
		test = Set(rows_t[tmp])
		if isempty(test)
			continue
		end
		tmp = nzrange(X, i)
		vals_d2_bar = vals[tmp]
		train = Set(rows[tmp])
		score = zeros(d2)
		ui = U[:, i]
		for j = 1:d2
            #if j in train || !(j in zero_items)
			if j in train
                score[j] = -10e10
				continue
			end
			vj = V[:, j]
			score[j] = dot(ui,vj)
		end
		p = sortperm(score, rev = true)
		for c = 1: K[length(K)]
			j = p[c]
			if score[j] == -10e10
				break
			end
			if j in test
				for k in length(K):-1:1
					if c <= K[k]
						precision[k] += 1
					else
						break
					end
				end
			end
		end
	end
	precision = precision./K/d1
	#precision = precision./K/500
    return precision[1], precision[2], precision[3]
end

train = "ml1m-50-oc/ml1m_oc_50_train_ratings.csv"
test = "ml1m-50-oc/ml1m_oc_50_test_ratings.csv"
X = readdlm(train, ',' , Int64);
x = vec(X[:,1]);
y = vec(X[:,2]);
v = vec(X[:,3]);
Y = readdlm(test, ',' , Int64);
xx = vec(Y[:,1]);
yy = vec(Y[:,2]);
vv = vec(Y[:,3]);
n = max(maximum(x), maximum(xx));
msize = max(maximum(y), maximum(yy));
println("Training dataset ", train, " and test dataset ", test, " are loaded. \n There are ", n, " users and ", msize, " items in the dataset.")
X = sparse(x, y, v, n, msize); # userid by movieid
Y = sparse(xx, yy, vv, n, msize);
# julia column major
# now moveid by userid
X = X';
Y = Y';
rows = rowvals(X);
vals = nonzeros(X);
cols = zeros(Int, size(vals)[1]);
d2, d1 = size(X);
cc = 0;
cold_start_s = Set{Int64}()
for i = 1:d1
	tmp = nzrange(X, i);
	nowlen = size(tmp)[1];
	if nowlen < 11
		push!(cold_start_s, i)
	end
	for j = 1:nowlen
		cc += 1
		cols[cc] = i
	end
end

rows_t = rowvals(Y);
vals_t = nonzeros(Y);
cols_t = zeros(Int, size(vals_t)[1]);
cc = 0;
for i = 1:d1
	tmp = nzrange(Y, i);
	nowlen = size(tmp)[1];
	for j = 1:nowlen
		cc += 1
		cols_t[cc] = i
	end
end

U_mf = readdlm("W_bpr.txt");
U_mf = U_mf[:,2:size(U_mf)[2]]
U_mf = U_mf';
ndcg_k = 10
r = 100

V_tmp = readdlm("H_bpr.txt");
V_mf = zeros(r, d2);
zero_items = Set{Int64}();
for i = 1:size(V_tmp)[1]
    item_id = convert(Int64, V_tmp[i][1]);
    push!(zero_items, item_id);
    V_mf[:,item_id] = V_tmp[i,:][2: (r + 1)]
end

p1, p5, p10 = compute_precision(U_mf, V_mf, X, Y, d1, d2, rows, vals, rows_t, vals_t, zero_items)
println("top 1 precision: ", p1, " top 5 precision: ", p5, " top 10 precision: ", p10)

#outfile = "outfile_sse_bpr.dat"
# writing to files is very similar:
#f = open(outfile, "a")
# both print and println can be used as usual but with f as their first arugment
#println(f, "[", p1, ",", p5, ",", p10, "],")
#close(f)


text:
# MonoidAlternative
# =================

a = "a" => [1]
b = "b" => [2,3,4]

@test a ⊕ b == ("ab" => [1,2,3,4])


text:
## LEARNING CURVES

"""
    curve = learning_curve(mach; resolution=30,
                                 resampling=Holdout(),
                                 repeats=1,
                                 measure=rms,
                                 weights=nothing,
                                 operation=predict,
                                 range=nothing,
                                 acceleration=default_resource(),
                                 acceleration_grid=CPU1(),
                                 rngs=nothing,
                                 rng_name=nothing)

Given a supervised machine `mach`, returns a named tuple of objects
suitable for generating a plot of performance estimates, as a function
of the single hyperparameter specified in `range`. The tuple `curve`
has the following keys: `:parameter_name`, `:parameter_scale`,
`:parameter_values`, `:measurements`.

To generate multiple curves for a `model` with a random number
generator (RNG) as a hyperparameter, specify the name of the (possibly
nested) RNG field, and a vector `rngs` of RNG's, one for each
curve. Alternatively, set `rngs` to the number of curves desired, in
which case RNG's are automatically generated. The individual curve
computations can be distributed across multiple processes using
`acceleration=CPUProcesses()`. See the second example below for a
demonstration.

```julia
X, y = @load_boston;
atom = @load RidgeRegressor pkg=MultivariateStats
ensemble = EnsembleModel(atom=atom, n=1000)
mach = machine(ensemble, X, y)
r_lambda = range(ensemble, :(atom.lambda), lower=10, upper=500, scale=:log10)
curve = learning_curve(mach; range=r_lambda, resampling=CV(), measure=mav)
using Plots
plot(curve.parameter_values,
     curve.measurements,
     xlab=curve.parameter_name,
     xscale=curve.parameter_scale,
     ylab = "CV estimate of RMS error")
```

If using a `Holdout()` `resampling` strategy (with no shuffling) and
if the specified hyperparameter is the number of iterations in some
iterative model (and that model has an appropriately overloaded
`MLJBase.update` method) then training is not restarted from scratch
for each increment of the parameter, ie the model is trained
progressively.

```julia
atom.lambda=200
r_n = range(ensemble, :n, lower=1, upper=250)
curves = learning_curve(mach; range=r_n, verbosity=0, rng_name=:rng, rngs=3)
plot!(curves.parameter_values,
     curves.measurements,
     xlab=curves.parameter_name,
     ylab="Holdout estimate of RMS error")

```
    learning_curve(model::Supervised, X, y; kwargs...)
    learning_curve(model::Supervised, X, y, w; kwargs...)

Plot a learning curve (or curves) directly, without first constructing
a machine.

"""
learning_curve(mach::Machine{<:Supervised}; kwargs...) =
    learning_curve(mach.model, mach.args...; kwargs...)

# for backwards compatibility
learning_curve!(mach::Machine{<:Supervised}; kwargs...) =
    learning_curve(mach; kwargs...)

function learning_curve(model::Supervised, args...;
                        resolution=30,
                        resampling=Holdout(),
                        weights=nothing,
                        measures=nothing,
                        measure=measures,
                        operation=predict,
                        ranges::Union{Nothing,ParamRange}=nothing,
                        range::Union{Nothing,ParamRange},
                        repeats=1,
                        acceleration=default_resource(),
                        acceleration_grid=CPU1(),
                        verbosity=1,
                        rngs=nothing,
                        rng_name=nothing,
                        check_measure=true)

    range !== nothing || error("No param range specified. Use range=... ")

    if rngs != nothing
        rng_name == nothing &&
            error("Having specified `rngs=...`, you must specify "*
                  "`rng_name=...` also. ")
        if rngs isa Integer
            rngs = MersenneTwister.(1:rngs)
        elseif rngs isa AbstractRNG
            rngs = [rngs, ]
        elseif !(rngs isa AbstractVector{<:AbstractRNG})
            error("`rng` must have type `Integer` , `AbstractRNG` or "*
                  "`AbstractVector{<:AbstractRNG}`. ")
        end
    end

    tuned_model = TunedModel(model=model,
                             range=range,
                             tuning=Grid(resolution=resolution,
                                         shuffle=false),
                             resampling=resampling,
                             operation=operation,
                             measure=measure,
                             train_best=false,
                             weights=weights,
                             repeats=repeats,
                             acceleration=acceleration_grid)

    tuned = machine(tuned_model, args...)

    results = _tuning_results(rngs, acceleration, tuned, rng_name, verbosity)

    parameter_name=results.parameter_names[1]
    parameter_scale=results.parameter_scales[1]
    parameter_values=[results.parameter_values[:, 1]...]
    measurements = results.measurements

    return (parameter_name=parameter_name,
            parameter_scale=parameter_scale,
            parameter_values=parameter_values,
            measurements=measurements)
end

_collate(plotting1, plotting2) =
    merge(plotting1,
          (measurements=hcat(plotting1.measurements,
                             plotting2.measurements),))

# fallback:
_tuning_results(rngs, acceleration, tuned, rngs_name, verbosity) =
    error("acceleration=$acceleration unsupported. ")

# single curve:
_tuning_results(rngs::Nothing, acceleration, tuned, rngs_name, verbosity) =
    _single_curve(tuned, verbosity)

function _single_curve(tuned, verbosity)
    fit!(tuned, verbosity=verbosity, force=true)
    tuned.report.plotting
end

# CPU1:
function _tuning_results(rngs::AbstractVector, acceleration::CPU1,
                         tuned, rng_name, verbosity)
    old_rng = recursive_getproperty(tuned.model.model, rng_name)

    ret = reduce(_collate,
                 [begin
                  recursive_setproperty!(tuned.model.model, rng_name, rng)
                  fit!(tuned, verbosity=verbosity, force=true)
                  tuned.report.plotting
                  end
                  for rng in rngs])

    recursive_setproperty!(tuned.model.model, rng_name, old_rng)

    return ret
end

# CPUProcesses:
function _tuning_results(rngs::AbstractVector, acceleration::CPUProcesses,
    tuned, rng_name, verbosity)

    old_rng = recursive_getproperty(tuned.model.model, rng_name)
    ret = @distributed (_collate) for rng in rngs
        recursive_setproperty!(tuned.model.model, rng_name, rng)
        fit!(tuned, verbosity=-1, force=true)
        tuned.report.plotting
    end
    recursive_setproperty!(tuned.model.model, rng_name, old_rng)

    return ret
end


text:
answer = 0

println(answer)

max_stars_repo_name: zvlex/euler-manager

max_issues_repo_name: zvlex/euler-manager

max_forks_repo_name: zvlex/euler-manager


text:
using CircGeometry, Test

between_buffer = 5.0
vf = 0.4
n_bodies = 50

# test circle outline
material = MaterialParameters(vf,n_bodies)
outline = OutlineCircle(1.0,Point(0.0,0.0))
ps = generate_porous_structure(outline,material,between_buffer;log=false)
vf_exp = compute_volume_fraction(ps,outline)

CircGeometry.copyArraysToCenters!(ps)

# test rectangle outline
outline = OutlineRectangle(Point(-2.0,0.0),Point(1.0,2.0))
@test CircGeometry.compute_outline_area(outline) == 6.0
ps = generate_porous_structure(outline,material,between_buffer;log=false)
vf_exp = compute_volume_fraction(ps,outline)

# test polygon outline
outline = OutlinePolygon([Point(-1.0,0.0),Point(1.0,0.0),
            Point(1.0,2.0),Point(0.0,1.0),Point(-1.0,2.0)])
@test CircGeometry.compute_outline_area(outline) == 3.0
ps = generate_porous_structure(outline,material,between_buffer;log=false)
vf_exp = compute_volume_fraction(ps,outline)

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# Keep ODE solution in a domain specified by a function. Inspired by:
# Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne, "Non-negative solutions
# of ODEs," Applied Mathematics and Computation Vol. 170, 2005, pp. 556-569.

# type definitions

abstract type AbstractDomainAffect{T,S,uType} end

struct PositiveDomainAffect{T,S,uType} <: AbstractDomainAffect{T,S,uType}
    abstol::T
    scalefactor::S
    u::uType
end

struct GeneralDomainAffect{autonomous,F,T,S,uType} <: AbstractDomainAffect{T,S,uType}
    g::F
    abstol::T
    scalefactor::S
    u::uType
    resid::uType

    function GeneralDomainAffect{autonomous}(g::F, abstol::T, scalefactor::S, u::uType,
                                             resid::uType) where {autonomous,F,T,S,uType}
        new{autonomous,F,T,S,uType}(g, abstol, scalefactor, u, resid)
    end
end

# definitions of callback functions

# Workaround since it is not possible to add methods to an abstract type:
# https://github.com/JuliaLang/julia/issues/14919
(f::PositiveDomainAffect)(integrator) = affect!(integrator, f)
(f::GeneralDomainAffect)(integrator) = affect!(integrator, f)

# general method defintions for domain callbacks

"""
    affect!(integrator, f::AbstractDomainAffect)

Apply domain callback `f` to `integrator`.
"""
function affect!(integrator, f::AbstractDomainAffect{T,S,uType}) where {T,S,uType}
    # modify u
    u_modified!(integrator, modify_u!(integrator, f))
    # define array of next time step, absolute tolerance, and scale factor
    if uType <: Nothing
        if typeof(integrator.u) <: Union{Number,SArray}
            u = integrator.u
        else
            u = similar(integrator.u)
        end
    else
        u = f.u
    end
    abstol = T <: Nothing ? integrator.opts.abstol : f.abstol
    scalefactor = S <: Nothing ? 1//2 : f.scalefactor

    # setup callback and save addtional arguments for checking next time step
    args = setup(f, integrator)

    # cache current time step
    dt = integrator.dt
    dt_modified = false
    p = integrator.p

    # update time step of integrator to proposed next time step
    integrator.dt = get_proposed_dt(integrator)

    # adjust time step to bounds and time stops
    fix_dt_at_bounds!(integrator)
    modify_dt_for_tstops!(integrator)
    t = integrator.t + integrator.dt

    while integrator.tdir * integrator.dt > 0
        # calculate estimated value of next step and its residuals
        if typeof(u) <: Union{Number,SArray}
            u = integrator(t)
        else
            integrator(u, t)
        end

        # check whether time step is accepted
        isaccepted(u, p, t, abstol, f, args...) && break

        # reduce time step
        dtcache = integrator.dt
        integrator.dt *= scalefactor
        dt_modified = true

        # adjust new time step to bounds and time stops
        fix_dt_at_bounds!(integrator)
        modify_dt_for_tstops!(integrator)
        t = integrator.t + integrator.dt

        # abort iteration when time step is not changed
        if dtcache == integrator.dt
            if integrator.opts.verbose
                @warn("Could not restrict values to domain. Iteration was canceled since ",
                     "time step dt = ", integrator.dt, " could not be reduced.")
            end
            break
        end
    end

    # update current and next time step
    if dt_modified # add safety factor since guess is based on extrapolation
        set_proposed_dt!(integrator, 9//10*integrator.dt)
    else
        set_proposed_dt!(integrator, integrator.dt)
    end
    integrator.dt = dt
end

"""
    modify_u!(integrator, f::AbstractDomainAffect)

Modify current state vector `u` of `integrator` if required, and return whether it actually
was modified.
"""
modify_u!(integrator, ::AbstractDomainAffect) = false

"""
    setup(f::AbstractDomainAffect, integrator)

Setup callback `f` and return an arbitrary tuple whose elements are used as additional
arguments in checking whether time step is accepted.
"""
setup(::AbstractDomainAffect, integrator) = ()

"""
    isaccepted(u, abstol, f::AbstractDomainAffect, args...)

Return whether `u` is an acceptable state vector at the next time point given absolute
tolerance `abstol`, callback `f`, and other optional arguments.
"""
isaccepted(u, p, t, tolerance, ::AbstractDomainAffect, args...) = true

# specific method definitions for positive domain callback

function modify_u!(integrator, f::PositiveDomainAffect)
    # set all negative values to zero
    _set_neg_zero!(integrator,integrator.u) # Returns true if modified
end

function _set_neg_zero!(integrator,u::AbstractArray)
    modified = false
    @inbounds for i in eachindex(integrator.u)
        if integrator.u[i] < 0
            integrator.u[i] = 0
            modified = true
        end
    end
    modified
end

function _set_neg_zero!(integrator,u::Number)
    modified = false
    if integrator.u < 0
        integrator.u = 0
        modified = true
    end
    modified
end

function _set_neg_zero!(integrator,u::SArray)
    modified = false
    @inbounds for i in eachindex(integrator.u)
        if u[i] < 0
            u = setindex(u,zero(first(u)),i)
            modified = true
        end
    end
    modified && (integrator.u = u)
    modified
end

# state vector is accepted if its entries are greater than -abstol
isaccepted(u, p, t, abstol::Number, ::PositiveDomainAffect) = all(x -> x + abstol > 0, u)
isaccepted(u, p, t, abstol, ::PositiveDomainAffect) = all(x + y > 0 for (x,y) in zip(u, abstol))

# specific method definitions for general domain callback

# create array of residuals
setup(f::GeneralDomainAffect, integrator) =
    typeof(f.resid) <: Nothing ? (similar(integrator.u),) : (f.resid,)

function isaccepted(u, p, t, abstol, f::GeneralDomainAffect{autonomous,F,T,S,uType},
                    resid) where {autonomous,F,T,S,uType}
    # calculate residuals
    if autonomous
        f.g(resid, u)
    else
        f.g(resid, u, p, t)
    end

    # accept time step if residuals are smaller than the tolerance
    if typeof(abstol) <: Number
        all(x-> x < abstol, resid)
    else
        # element-wise comparison
        all(x < y for (x,y) in zip(resid, abstol))
    end
end

# callback definitions

function GeneralDomain(g, u=nothing; nlsolve=NLSOLVEJL_SETUP(), save=true,
                       abstol=nothing, scalefactor=nothing, autonomous=DiffEqBase.numargs(g)==2,
                       nlopts=Dict(:ftol => 10*eps()))
    if typeof(u) <: Nothing
        affect! = GeneralDomainAffect{autonomous}(g, abstol, scalefactor, nothing, nothing)
    else
        affect! = GeneralDomainAffect{autonomous}(g, abstol, scalefactor, deepcopy(u),
                                                  deepcopy(u))
    end
    condition = (u,t,integrator) -> true
    CallbackSet(ManifoldProjection(g; nlsolve=nlsolve, save=false,
                                   autonomous=autonomous, nlopts=nlopts),
                DiscreteCallback(condition, affect!; save_positions=(false, save)))
end

function PositiveDomain(u=nothing; save=true, abstol=nothing, scalefactor=nothing)
    if typeof(u) <: Nothing
        affect! = PositiveDomainAffect(abstol, scalefactor, nothing)
    else
        affect! = PositiveDomainAffect(abstol, scalefactor, deepcopy(u))
    end
    condition = (u,t,integrator) -> true
    DiscreteCallback(condition, affect!; save_positions=(false, save))
end

export GeneralDomain, PositiveDomain