text:
# Globals
const FONTSCALE = MJCore.FONTSCALE_200  # can be 100, 150, 200
const maxgeom = 5000                    # preallocated geom array in mjvScene

include("sensor.jl")

include("robot_commands.jl")

# Input handlers
include("keyboard_handler.jl")
include("mouse_handler.jl")
include("joystick_handler.jl")

# Utility functions for working with the model's cameras
function list_named_cameras(m::jlModel)
    println("Total number of named cameras: ", m.ncam)
    for i = 1:min(m.ncam, MJCore.MAXUIMULTI-2)
        println("\t$i: $(mj_id2name(m, MJCore.mjOBJ_CAMERA, i-1))")
    end
end

function get_named_cam_id(m::jlModel, camera_name::AbstractString)
    # returns -1 if model does not contain the named camera
    Int(mj_name2id(m, MJCore.mjOBJ_CAMERA, camera_name)::Int32)
end

include("imgui_setup.jl")

# Set up simulator and GLFW window environments
function setup(mm::jlModel, dd::jlData, width=1200, height=900) # TODO named args for callbacks
    GLFW.WindowHint(GLFW.SAMPLES, 4)
    GLFW.WindowHint(GLFW.VISIBLE, 1)

    s = mjSim(mm, dd, "Simulate"; width=width, height=height)

    @info("Refresh Rate: $(s.refreshrate)")
    @info("Resolution: $(width)x$(height)")

    # Make the window's context current
    GLFW.MakeContextCurrent(s.window)
    GLFW.SwapInterval(1) # enable vsync

    # init abstract visualization
    mjv_defaultCamera(s.cam)
    mjv_defaultOption(s.vopt)
    #profilerinit()
    sensorinit(s)

    # make empty scene
    mjv_defaultScene(s.scn)
    mjv_makeScene(s.m, s.scn, maxgeom)

    # mujoco setup
    mjv_defaultPerturb(s.pert)
    mjr_defaultContext(s.con)
    mjr_makeContext(s.m, s.con, FONTSCALE)  # model specific setup

    alignscale(s)
    mjv_updateScene(s.m, s.d, s.vopt, s.pert, s.cam, MJCore.mjCAT_ALL, s.scn)

    # List the named cameras
    # list_named_cameras(s.m)

    # MWB !!! TODO: Generalize this. Name should come from a config file, or  if there are
    # only 3 cameras total (-1, 0, and 1), then we probably know that s.camid should be 1.
    # Get the id of the head mounted camera
    s.camid = get_named_cam_id(s.m, "puppercam")

    # Set up the CImGui environment
    imgui_setup(s.window)

    # create texture for image drawing
    s.camviewid = ImGui_ImplOpenGL3_CreateImageTexture(s.camviewport.width, s.camviewport.height, format=GL_RGB)

    # Set up GLFW callbacks
    GLFW.SetInputMode(s.window, GLFW_LOCK_KEY_MODS, true)
    GLFW.SetKeyCallback(s.window, (w,k,sc,a,m)->keyboard(s,w,k,sc,a,m))
    GLFW.SetCursorPosCallback(s.window, (w,x,y)->mouse_move(s,w,x,y))
    GLFW.SetMouseButtonCallback(s.window, (w,b,a,m)->mouse_button(s,w,b,a,m))
    GLFW.SetScrollCallback(s.window, (w,x,y)->scroll(s,w,x,y))
    ##GLFW.SetDropCallback(s.window, drop)

   return s
end


text:
function sierpinski(n, token::AbstractString="*")
	x = fill(token, 1, 1)
	for _ in 1:n
		h, w = size(x)
		s = fill(" ", h,(w + 1) ÷ 2)
		t = fill(" ", h,1)
		x = [[s x s] ; [x t x]]
	end
	return x
end

function printsierpinski(m::Matrix)
    for r in 1:size(m, 1)
        println(join(m[r, :]))
    end
end

sierpinski(4) |> printsierpinski

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Pkg
using SafeTestsets
const LONGER_TESTS = false

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

#Start Test Script

@time begin

if GROUP == "All" || GROUP == "InterfaceI" || GROUP == "Interface"
  #@time @safetestset "Linear Solver Tests" begin include("interface/linear_solver_test.jl") end
  @time @safetestset "Discrete Algorithm Tests" begin include("interface/discrete_algorithm_test.jl") end
  @time @safetestset "Tstops Tests" begin include("interface/ode_tstops_tests.jl") end
  @time @safetestset "Backwards Tests" begin include("interface/ode_backwards_test.jl") end
  @time @safetestset "Initdt Tests" begin include("interface/ode_initdt_tests.jl") end
  @time @safetestset "Linear Tests" begin include("interface/ode_twodimlinear_tests.jl") end
  @time @safetestset "Differentiation Trait Tests" begin include("interface/differentiation_traits_tests.jl") end
  @time @safetestset "Inf Tests" begin include("interface/inf_handling.jl") end
  @time @safetestset "Jacobian Tests" begin include("interface/jacobian_tests.jl") end
  @time @safetestset "saveat Tests" begin include("interface/ode_saveat_tests.jl") end
  @time @safetestset "save_idxs Tests" begin include("interface/ode_saveidxs_tests.jl") end
  @time @safetestset "Scalar Handling Tests" begin include("interface/scalar_handling_tests.jl") end
  @time @safetestset "Static Array Tests" begin include("interface/static_array_tests.jl") end
  @time @safetestset "u_modifed Tests" begin include("interface/umodified_test.jl") end
  @time @safetestset "Composite Algorithm Tests" begin include("interface/composite_algorithm_test.jl") end
  @time @safetestset "Complex Tests" begin include("interface/complex_tests.jl") end
  @time @safetestset "Ndim Complex Tests" begin include("interface/ode_ndim_complex_tests.jl") end
  @time @safetestset "Number Type Tests" begin include("interface/ode_numbertype_tests.jl") end
  @time @safetestset "Stiffness Detection Tests" begin include("interface/stiffness_detection_test.jl") end
  @time @safetestset "Composite Interpolation Tests" begin include("interface/composite_interpolation.jl") end
  @time @safetestset "Export tests" begin include("interface/export_tests.jl") end
  @time @safetestset "Type Handling Tests" begin include("interface/type_handling.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceII" || GROUP == "Interface")
  @time @safetestset "Linear Nonlinear Solver Tests" begin include("interface/linear_nonlinear_tests.jl") end
  @time @safetestset "Linear Solver Tests" begin include("interface/linear_solver_test.jl") end
  @time @safetestset "Sparse Diff Tests" begin include("interface/sparsediff_tests.jl") end
  @time @safetestset "Enum Tests" begin include("interface/enums.jl") end
  @time @safetestset "Mass Matrix Tests" begin include("interface/mass_matrix_tests.jl") end
  @time @safetestset "Second Order with First Order Solver Tests" begin include("interface/second_order_with_first_order_solvers.jl") end
  @time @safetestset "AD Tests" begin include("interface/ad_tests.jl") end
  @time @safetestset "DAE Initialize Integration" begin include("interface/dae_initialize_integration.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceIII" || GROUP == "Interface")
  @time @safetestset "Derivative Utilities Tests" begin include("interface/utility_tests.jl") end
  @time @safetestset "DEStats Tests" begin include("interface/destats_tests.jl") end
  @time @safetestset "No Index Tests" begin include("interface/noindex_tests.jl") end
  @time @safetestset "No Jac Tests" begin include("interface/nojac.jl") end
  @time @safetestset "Units Tests" begin include("interface/units_tests.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Integrators_I" || GROUP == "Integrators")
  @time @safetestset "Reinit Tests" begin include("integrators/reinit_test.jl") end
  @time @safetestset "Events Tests" begin include("integrators/ode_event_tests.jl") end
  @time @safetestset "Alg Events Tests" begin include("integrators/alg_events_tests.jl") end
  @time @safetestset "Autodiff Events Tests" begin include("integrators/autodiff_events.jl") end
  @time @safetestset "Discrete Callback Dual Tests" begin include("integrators/discrete_callback_dual_test.jl") end
  @time @safetestset "Iterator Tests" begin include("integrators/iterator_tests.jl") end
  @time @safetestset "Integrator Interface Tests" begin include("integrators/integrator_interface_tests.jl") end
  @time @safetestset "Error Check Tests" begin include("integrators/check_error.jl") end
  @time @safetestset "Event Detection Tests" begin include("integrators/event_detection_tests.jl") end
  @time @safetestset "Event Repetition Detection Tests" begin include("integrators/event_repeat_tests.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Integrators_II" || GROUP == "Integrators")
  @time @safetestset "Reverse Directioned Event Tests" begin include("integrators/rev_events_tests.jl") end
  @time @safetestset "Differentiation Direction Tests" begin include("integrators/diffdir_tests.jl") end
  @time @safetestset "Resize Tests" begin include("integrators/resize_tests.jl") end
  @time @safetestset "DAE Initialization Tests" begin include("integrators/dae_initialization_tests.jl") end
  @time @safetestset "Cache Tests" begin include("integrators/ode_cache_tests.jl") end
  @time @safetestset "Add Steps Tests" begin include("integrators/ode_add_steps_tests.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Regression_I" || GROUP == "Regression")
  @time @safetestset "Dense Tests" begin include("regression/ode_dense_tests.jl") end
  @time @safetestset "Special Interp Tests" begin include("regression/special_interps.jl") end
  @time @safetestset "Inplace Tests" begin include("regression/ode_inplace_tests.jl") end
  @time @safetestset "Adaptive Tests" begin include("regression/ode_adaptive_tests.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Regression_II" || GROUP == "Regression")
  @time @safetestset "PSOS Energy Conservation Tests" begin include("regression/psos_and_energy_conservation.jl") end
  @time @safetestset "Unrolled Tests" begin include("regression/ode_unrolled_comparison_tests.jl") end
  @time @safetestset "Time derivative Tests" begin include("regression/time_derivative_test.jl") end
  @time @safetestset "IIP vs OOP Tests" begin include("regression/iipvsoop_tests.jl") end
  @time @safetestset "Inference Tests" begin include("regression/inference.jl") end
end

if !is_APPVEYOR && GROUP == "AlgConvergence_I"
  @time @safetestset "Partitioned Methods Tests" begin include("algconvergence/partitioned_methods_tests.jl") end
  @time @safetestset "Convergence Tests" begin include("algconvergence/ode_convergence_tests.jl") end
  @time @safetestset "DAE Convergence Tests" begin include("algconvergence/dae_convergence_tests.jl") end
  @time @safetestset "Non-autonomous Convergence Tests" begin include("algconvergence/non-autonomous_convergence_tests.jl") end
  @time @safetestset "Adams Variable Coefficients Tests" begin include("algconvergence/adams_tests.jl") end
  @time @safetestset "Nordsieck Tests" begin include("algconvergence/nordsieck_tests.jl") end
end

if !is_APPVEYOR && GROUP == "AlgConvergence_II"
  @time @safetestset "SSPRK Tests" begin include("algconvergence/ode_ssprk_tests.jl") end
  @time @safetestset "Low Storage RK Tests" begin include("algconvergence/ode_low_storage_rk_tests.jl") end
  @time @safetestset "OwrenZen Tests" begin include("algconvergence/owrenzen_tests.jl") end
  @time @safetestset "Runge-Kutta-Chebyshev Tests" begin include("algconvergence/rkc_tests.jl") end
end

if !is_APPVEYOR && GROUP == "AlgConvergence_III"
  @time @safetestset "Linear Methods Tests" begin include("algconvergence/linear_method_tests.jl") end
  @time @safetestset "Split Methods Tests" begin include("algconvergence/split_methods_tests.jl") end
  @time @safetestset "Rosenbrock Tests" begin include("algconvergence/ode_rosenbrock_tests.jl") end
  @time @safetestset "FIRK Tests" begin include("algconvergence/ode_firk_tests.jl") end
  @time @safetestset "Linear-Nonlinear Methods Tests" begin include("algconvergence/linear_nonlinear_convergence_tests.jl") end
  @time @safetestset "Linear-Nonlinear Krylov Methods Tests" begin include("algconvergence/linear_nonlinear_krylov_tests.jl") end
  @time @safetestset "Feagin Tests" begin include("algconvergence/ode_feagin_tests.jl") end
  @time @safetestset "Extrapolation Tests" begin include("algconvergence/ode_extrapolation_tests.jl") end
  @time @safetestset "Symplectic Tests" begin include("algconvergence/symplectic_tests.jl") end
end

if !is_APPVEYOR && GROUP == "Downstream"
  activate_downstream_env()
  @time @safetestset "DelayDiffEq Tests" begin include("downstream/delaydiffeq.jl") end
end

if !is_APPVEYOR && GROUP == "ODEInterfaceRegression"
  activate_downstream_env()
  @time @safetestset "Init dt vs dorpri tests" begin include("odeinterface/init_dt_vs_dopri_tests.jl") end
  @time @safetestset "ODEInterface Regression Tests" begin include("odeinterface/odeinterface_regression.jl") end
end

if !is_APPVEYOR && GROUP == "Multithreading"
  @time @safetestset "Extrapolation Tests" begin include("multithreading/ode_extrapolation_tests.jl") end
end

if !is_APPVEYOR && GROUP == "GPU"
    activate_downstream_env()
    @time @safetestset "Simple GPU" begin
    import OrdinaryDiffEq
    include(joinpath(dirname(pathof(OrdinaryDiffEq.DiffEqBase)), "..", "test/gpu/simple_gpu.jl"))
  end
  @time @safetestset "Linear LSRK GPU" begin include("gpu/linear_lsrk.jl") end
  @time @safetestset "Reaction-Diffusion Stiff Solver GPU" begin include("gpu/reaction_diffusion_stiff.jl") end
end

end # @time


text:
using OrdinaryDiffEq, Zygote
using DiffEqSensitivity, Test, ForwardDiff

abstol = 1e-12
reltol = 1e-12
savingtimes = 0.5

function test_discrete_callback(cb, tstops, g, dg!, cboop=nothing, tprev=false)
  function fiip(du, u, p, t)
    du[1] = dx = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = dy = -p[3] * u[2] + p[4] * u[1] * u[2]
  end
  function foop(u, p, t)
    dx = p[1] * u[1] - p[2] * u[1] * u[2]
    dy = -p[3] * u[2] + p[4] * u[1] * u[2]
    [dx, dy]
  end

  p = [1.5, 1.0, 3.0, 1.0]
  u0 = [1.0; 1.0]

  prob = ODEProblem(fiip, u0, (0.0, 10.0), p)
  proboop = ODEProblem(foop, u0, (0.0, 10.0), p)

  sol1 = solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)
  sol2 = solve(prob, Tsit5(), u0=u0, p=p, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)

  if cb.save_positions == [1, 1]
    @test length(sol1.t) != length(sol2.t)
  else
    @test length(sol1.t) == length(sol2.t)
  end

  du01, dp1 = Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
    u0, p)

  du01b, dp1b = Zygote.gradient(
    (u0, p) -> g(solve(proboop, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint())),
    u0, p)

  du01c, dp1c = Zygote.gradient(
    (u0, p) -> g(solve(proboop, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=BacksolveAdjoint(checkpointing=false))),
    u0, p)

  if cboop === nothing
    du02, dp2 = Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=ReverseDiffAdjoint())), u0, p)
  else
    du02, dp2 = Zygote.gradient(
      (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cboop, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=ReverseDiffAdjoint())), u0, p)
  end

  du03, dp3 = Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=InterpolatingAdjoint(checkpointing=true))),
    u0, p)

  du03c, dp3c = Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=InterpolatingAdjoint(checkpointing=false))),
    u0, p)

  du04, dp4 = Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes, sensealg=QuadratureAdjoint())),
    u0, p)

  dstuff = ForwardDiff.gradient(
    (θ) -> g(solve(prob, Tsit5(), u0=θ[1:2], p=θ[3:6], callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)),
    [u0; p])

  @info dstuff

  # tests wrt discrete sensitivities
  if tprev
    # tprev depends on stepping behaviour of integrator. Thus sensitivities are necessarily (slightly) different.
    @test du02 ≈ dstuff[1:2] rtol = 1e-3
    @test dp2 ≈ dstuff[3:6] rtol = 1e-3
    @test du01 ≈ dstuff[1:2] rtol = 1e-3
    @test dp1 ≈ dstuff[3:6] rtol = 1e-3
    @test du01 ≈ du02 rtol = 1e-3
    @test dp1 ≈ dp2 rtol = 1e-3
  else
    @test du02 ≈ dstuff[1:2]
    @test dp2 ≈ dstuff[3:6]
    @test du01 ≈ dstuff[1:2]
    @test dp1 ≈ dstuff[3:6]
    @test du01 ≈ du02
    @test dp1 ≈ dp2
  end

  # tests wrt continuous sensitivities
  @test du01b ≈ du01
  @test dp1b ≈ dp1
  @test du01c ≈ du01
  @test dp1c ≈ dp1
  @test du01 ≈ du03 rtol = 1e-7
  @test du01 ≈ du03c rtol = 1e-7
  @test du03 ≈ du03c
  @test du01 ≈ du04
  @test dp1 ≈ dp3
  @test dp1 ≈ dp3c
  @test dp1 ≈ dp4 rtol = 1e-7

  cb2 = DiffEqSensitivity.track_callbacks(CallbackSet(cb), prob.tspan[1], prob.u0, prob.p, BacksolveAdjoint())
  sol_track = solve(prob, Tsit5(), u0=u0, p=p, callback=cb2, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)
  #cb_adj = DiffEqSensitivity.setup_reverse_callbacks(cb2,BacksolveAdjoint())

  adj_prob = ODEAdjointProblem(sol_track, BacksolveAdjoint(), dg!, sol_track.t, nothing,
    callback=cb2,
    abstol=abstol, reltol=reltol)
  adj_sol = solve(adj_prob, Tsit5(), abstol=abstol, reltol=reltol)
  @test du01 ≈ -adj_sol[1:2, end]
  @test dp1 ≈ adj_sol[3:6, end]

  # adj_prob = ODEAdjointProblem(sol_track,InterpolatingAdjoint(),dg!,sol_track.t,nothing,
  # 						 callback = cb2,
  # 						 abstol=abstol,reltol=reltol)
  # adj_sol = solve(adj_prob, Tsit5(), abstol=abstol,reltol=reltol)
  #
  # @test du01 ≈ -adj_sol[1:2,end]
  # @test dp1 ≈ adj_sol[3:6,end]
end

@testset "Discrete callbacks" begin
  @testset "ODEs" begin
    println("ODEs")
    @testset "simple loss function" begin
      g(sol) = sum(sol)
      function dg!(out, u, p, t, i)
        (out .= -1)
      end
      @testset "callbacks with no effect" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 0.0
        cb = DiscreteCallback(condition, affect!, save_positions=(false, false))
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callbacks with no effect except saving the state" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 0.0
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 2.0
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callback at multiple time points" begin
        affecttimes = [2.03, 4.0, 8.0]
        condition(u, t, integrator) = t ∈ affecttimes
        affect!(integrator) = integrator.u[1] += 2.0
        cb = DiscreteCallback(condition, affect!)
        test_discrete_callback(cb, affecttimes, g, dg!)
      end
      @testset "state-dependent += callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (integrator.u .+= integrator.p[2] / 8 * sin.(integrator.u))
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "other callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (integrator.u[1] = 2.0; @show "triggered!")
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "parameter changing callback at single time point" begin
        condition(u, t, integrator) = t == 5.1
        affect!(integrator) = (integrator.p .= 2 * integrator.p .- 0.5)
        affect(integrator) = (integrator.p = 2 * integrator.p .- 0.5)
        cb = DiscreteCallback(condition, affect!)
        cboop = DiscreteCallback(condition, affect)
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.1]
        test_discrete_callback(cb, tstops, g, dg!, cboop)
      end
      @testset "tprev dependent callback" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (@show integrator.tprev; integrator.u[1] += integrator.t - integrator.tprev)
        cb = DiscreteCallback(condition, affect!)
        tstops = [4.999, 5.0]
        test_discrete_callback(cb, tstops, g, dg!, nothing, true)
      end
    end
    @testset "MSE loss function" begin
      g(u) = sum((1.0 .- u) .^ 2) ./ 2
      dg!(out, u, p, t, i) = (out .= 1.0 .- u)
      @testset "callbacks with no effect" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 0.0
        cb = DiscreteCallback(condition, affect!, save_positions=(false, false))
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callbacks with no effect except saving the state" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 0.0
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = integrator.u[1] += 2.0
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "callback at multiple time points" begin
        affecttimes = [2.03, 4.0, 8.0]
        condition(u, t, integrator) = t ∈ affecttimes
        affect!(integrator) = integrator.u[1] += 2.0
        cb = DiscreteCallback(condition, affect!)
        test_discrete_callback(cb, affecttimes, g, dg!)
      end
      @testset "state-dependent += callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (integrator.u .+= integrator.p[2] / 8 * sin.(integrator.u))
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "other callback at single time point" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (integrator.u[1] = 2.0; @show "triggered!")
        cb = DiscreteCallback(condition, affect!)
        tstops = [5.0]
        test_discrete_callback(cb, tstops, g, dg!)
      end
      @testset "parameter changing callback at single time point" begin
        condition(u, t, integrator) = t == 5.1
        affect!(integrator) = (integrator.p .= 2 * integrator.p .- 0.5)
        affect(integrator) = (integrator.p = 2 * integrator.p .- 0.5)
        cb = DiscreteCallback(condition, affect!)
        cboop = DiscreteCallback(condition, affect)
        tstops = [5.1]
        test_discrete_callback(cb, tstops, g, dg!, cboop)
      end
      @testset "tprev dependent callback" begin
        condition(u, t, integrator) = t == 5
        affect!(integrator) = (@show integrator.tprev; integrator.u[1] += integrator.t - integrator.tprev)
        cb = DiscreteCallback(condition, affect!)
        tstops = [4.999, 5.0]
        test_discrete_callback(cb, tstops, g, dg!, nothing, true)
      end
    end
  end
end


text:
"""
    remake_prob_with_jac(prob)

Add a jacobian function to the ODE function, using ModelingToolkit.jl.
"""
function remake_prob_with_jac(prob::ODEProblem)
    IIP = isinplace(prob.f)
    try
        p = prob.p isa DiffEqBase.NullParameters ? [] : collect(prob.p)
        prob = remake(prob, p=p)
        sys = modelingtoolkitize(prob)
        jac = eval(ModelingToolkit.generate_jacobian(sys)[IIP ? 2 : 1])
        f = ODEFunction{IIP}(
            prob.f.f,
            jac=jac;
            analytic=prob.f.analytic,
        )
        return remake(prob, f=f)
    catch
        error("Could not generate a jacobian for the problem")
    end
end


text:
using DelimitedFiles

mutable struct AceData
    #=
    Structure to contain information found in ace files

    IZ: List of 16 nuclide ZA identifiers. Needed for  S(α, β).
    AW: List of 16 atomic weight ratios to go with IZ. Needed for  S(α, β).
    NXS: List of 16 integers defining features of the ACE data file.
    JXS: List of 32 integers locating different parts of the ACE data file.
    XSS: List of floats containing the majority of the nuclear data.
    =#
    IZ::Vector{Int64}
    AW::Vector{Float64}
    NXS::Vector{Int64}
    JXS::Vector{Int64}
    XSS::Vector{Float64}
end

function read_ace(filename::String)
    #=
    Creates an AceData struct based on the path to an ACE data file.

    :param filename: String for the path to the ACE data file.
    :return acedata: AceData struct containing all the information in the file.
    =#
    file_lines = readlines(filename)
    pos = 1

    if length(split(file_lines[1])) == 3
        pos = pos + 2
    end
    HZ, AW, TZ, HD = split(file_lines[pos], " ", keepempty = false)
    HK = file_lines[pos+1][1:73]
    HM = file_lines[pos+1][74:80]
    pos = pos + 2
    IZAW = split(join(file_lines[pos:pos+3]))
    IZ = parse.(Int, IZAW[1:2:end])
    AW = parse.(Float64, IZAW[2:2:end])
    pos = pos + 4
    NXS = parse.(Int, split(join(file_lines[pos:pos+1])))
    pos = pos + 2
    JXS = parse.(Int, split(join(file_lines[pos:pos+3])))
    pos = pos + 4
    XSS = parse.(Float64, split(join(file_lines[pos:end])))
    #XSS = vec(readdlm(filename, Float64, skipstart = pos-1))
    # Return data
    acedata = AceData(IZ, AW, NXS, JXS, XSS)
    return acedata
end

function E_grid(acedata::AceData)
    #=
    Returns the full common energy grid from an AceData struct.

    :param acedata: AceData struct containing the relevant information
    :return: commom energy grid in the ACE data.
    =#
    NES = acedata.NXS[3]
    return acedata.XSS[acedata.JXS[1] : acedata.JXS[1] + NES - 1]
end

function E_idx(acedata::AceData, idx::Int)
    #=
    Returns energy at a specified index in the grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    if idx > get_NES(acedata)
        throw(ArgumentError("Index out of grid range."))
    end
    return acedata.XSS[acedata.JXS[1] + idx - 1]
end

function σ_T_idx(acedata::AceData, idx::Int)
    #=
    Returns microscopic total cross section at a specified index in the grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    NES = get_NES(acedata)
    if idx > NES
        throw(ArgumentError("Index out of grid range."))
    end
    return acedata.XSS[NES + acedata.JXS[1] + idx - 1]
end

function σ_a_idx(acedata::AceData, idx::Int)
    #=
    Returns microscopic absorption cross section at a specified index in the
    grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    NES = get_NES(acedata)
    if idx > NES
        throw(ArgumentError("Index out of grid range."))
    end
    return acedata.XSS[2*NES + acedata.JXS[1] + idx - 1]
end

function σ_el_idx(acedata::AceData, idx::Int)
    #=
    Returns microscopic elastic scattering cross section at a specified index in
    the grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    NES = get_NES(acedata)
    if idx > NES
        throw(ArgumentError("Index out of grid range."))
    end
    return acedata.XSS[4*NES + acedata.JXS[1] + idx - 1]
end

function σ_f_idx(acedata::AceData, idx::Int)
    #=
    Returns microscopic total fission cross section at a specified index in the
    grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    F = acedata.JXS[21]
    if F == 0
        return 0
    else
        IE = convert(Int, acedata.XSS[F])
        if idx < IE
            return 0
        else
            NE = convert(Int, acedata.XSS[F + 1])
            if idx > NE + IE - 1
                throw(ArgumentError("Index out of grid range for FIS."))
            end
        end
        return acedata.XSS[F + acedata.JXS[1] + idx - IE]
    end
end

function heat_idx(acedata::AceData, idx::Int)
    #=
    Returns microscopic heating cross section at a specified index in the grid.

    :param acedata: AceData struct containing relevant information
    :param idx: index of interset
    =#
    NES = get_NES(acedata)
    if idx > NES
        throw(ArgumentError("Index out of grid range."))
    end
    return acedata.XSS[5*NES + acedata.JXS[1] + idx - 1]
end

function get_NES(acedata::AceData)
    #=
    Returns the length of the energy grid.

    :param acedata: AceData struct containing relevant information
    =#
    return acedata.NXS[3]
end

function find_idx(
    acedata::AceData,
    E::Float64,
    lo::Int = 1,
    hi::Int = get_NES(acedata))
    #=
    Returns the index of the enegy space directly below or equal to the given
    energy. Uses recursive binary search algorithm.

    :param acedata: AceData struct containing relevant information
    :param E: Energy of interest
    :param lo: Lower bound index of search
    param hi: Upper bound index of search
    =#
    # Base case to terminate search
    if hi - lo == 1
        return lo
    end
    # Split search down middle
    mid = (lo + hi) ÷ 2
    # Recursive cases if mid point is too high or low
    if E_idx(acedata, mid) > E
        return find_idx(acedata, E, lo, mid)
    else
        return find_idx(acedata, E, mid, hi)
    end
end

function nu_prompt(acedata::AceData, E::Float64)

end

function interp_frac_log(
    acedata::AceData,
    E::Number,
    idx::Int)
    #=
    Returns the logarithmic fraction between the adjacent energy points in the
    grid used for interpolation

    :param acedata: AceData struct containing relevant information
    :param E: Energy of interest
    :param idx: Energy grid index directly below the energy of interest
    =#
    E0 = E_idx(acedata, idx)
    E1 = E_idx(acedata, idx + 1)
    return log(E/E0) / log(E1/E0)
end

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
function sieve2(n :: Int)
    ni = (n - 1) ÷ 2
    isprime = trues(ni)
    for i in 1:ni
        if isprime[i]
            j = 2i * (i + 1)
            if j > ni
                m = findall(isprime)
                map!((i::Int) -> 2i + 1, m, m)
                return pushfirst!(m, 2)
            else
                p = 2i + 1
                while j <= ni
                  isprime[j] = false
                  j += p
                end
            end
        end
    end
end


text:
# Add method to resolve exponentiation ambiguity.
^(n::Node{<:Real}, p::Integer) = invoke(^, Tuple{Node{<:Real}, Real}, n, p)

import Base: float
@explicit_intercepts float Tuple{∇ArrayOrScalar}
∇(::typeof(float), ::Type{Arg{1}}, p, y, ȳ, x) = float(ȳ)


text:
legend_spec =
[
	(:size, AbstractString, nothing),
	(:shape, AbstractString, nothing),
	(:fill, AbstractString, nothing),
	(:stroke, AbstractString, nothing),
	(:orient, AbstractString, nothing),
	(:offset, Number, nothing),
	(:title, AbstractString, nothing),
	(:format, AbstractString, nothing),
	(:formatType, AbstractString, nothing),
	(:values, Vector, nothing),
	(:properties, Union{Dict{Any, Any}, VegaMarkPropertySet}, nothing)
]

primitivefactory(:VegaLegend, legend_spec)

max_stars_repo_name: jpthor/vega


text:
# SPDX-License-Identifier: MIT

using ArgParse
using HDF5
using SpeechFeatures
using TOML
using WAV

#######################################################################
# Default values

const FEATYPES = Dict(
    "MFCC" => (MFCC,
               ("removedc", "dithering", "srate", "frameduration",
                "framestep", "preemphasis", "windowfn", "windowpower",
                "nfilters", "lofreq", "hifreq", "addenergy", "nceps",
                "liftering")
              ),
    "FBANK" => (LogMelSpectrum,
                ("removedc", "dithering", "srate", "frameduration",
                 "framestep", "preemphasis", "windowfn", "windowpower",
                 "nfilters", "lofreq", "hifreq")
               ),
    "MAGSPEC" => (LogMagnitudeSpectrum,
                  ("removedc", "dithering", "srate", "frameduration",
                   "framestep", "preemphasis", "windowfn", "windowpower")
                 )
)

const WINDOWS = Dict(
    "hann" => SpeechFeatures.HannWindow,
    "hamming" => SpeechFeatures.HammingWindow,
    "rectangular" => SpeechFeatures.RectangularWindow
)

const DEFAULTS = Dict(
    "removedc" => true,
    "dithering" => 0,
    "srate" => 16000,
    "frameduration" => 0.025,
    "framestep" => 0.01,
    "preemphasis" => 0.97,
    "windowfn" => "hann",
    "windowpower" => 0.85,
    "nfilters" => 26,
    "lofreq" => 80,
    "hifreq" => 7600,
    "addenergy" => true,
    "nceps" => 12,
    "liftering" => 22
)

const DEFAULTS_DELTA = Dict(
    "order" => 2,
    "deltawin" => 2
)

const DEFAULTS_NORM = Dict(
    "mean_norm" => true,
)

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

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--compression-level", "-c"
            arg_type = Int
            default = 0
            help = "compression level of the HDF5 archive"
        "configfile"
            required = true
            help = "features configuration file (TOML format)"
        "scp"
            required = true
            help = "SCP file to extract the features"
        "out"
            required = true
            help = "output HDF5 archive"
    end
    parse_args(s)
end

loadpath(path) = wavread(path, format = "double")

loadpipe(pipe) = channels, srate = wavread(IOBuffer(read(cmd)))
function loadpipe(pipe)
    subcmds = [`$(split(subcmd))` for subcmd in split(pipe[1:end-1], "|")]
    cmd = pipeline(subcmds...)
    wavread(IOBuffer(read(cmd)))
end

load(path_or_pipe) =
    endswith(path_or_pipe, "|") ? loadpipe(path_or_pipe) : loadpath(path_or_pipe)

function load_scpentry(line)
    tokens = split(strip(line))
    uttid = tokens[1]
    path_or_pipe = join(tokens[2:end], " ")
    String(uttid), path_or_pipe
end

function buildextractor(conf)
    featype, options = FEATYPES[conf["featype"]]
    feaopts = Dict()
    for opt in options
        if opt == "windowfn"
            feaopts[Symbol(opt)] = WINDOWS[conf[opt]]
        else
            feaopts[Symbol(opt)] = conf[opt]
        end
    end
    feaextractor = featype(;feaopts...)

    if "DeltaCoeffs" in keys(conf)
        deltaopts = merge(DEFAULTS_DELTA, conf["DeltaCoeffs"])
        Δ_extractor = DeltaCoeffs(order = deltaopts["order"],
                                  deltawin = deltaopts["deltawin"])
        feaextractor = Δ_extractor ∘ feaextractor
    end

    if "Normalization" in keys(conf)
        normopts = merge(DEFAULTS_NORM, conf["Normalization"])
        if normopts["mean_norm"]
            feaextractor = MeanNorm() ∘ feaextractor
        end
    end

    feaextractor
end

function main(args)
    conf = merge(DEFAULTS, TOML.parsefile(args["configfile"]))
    feaextractor = buildextractor(conf)

    @debug "extracting features to $(args["out"])"

    h5open(args["out"], "w") do f

        for (key, val) in conf
            if typeof(val) <: Dict
                for (skey, sval) in val
                    attributes(f)["$key.$skey"] = sval
                end
            else
                attributes(f)[key] = val
            end
        end

        for line in readlines(args["scp"])
            @debug "processing $line"
            uttid, path_or_pipe = load_scpentry(line)
            channels, srate = load(path_or_pipe)

            @assert srate == conf["srate"]
            @assert size(channels, 2) == 1

            fea = channels[:,1] |> feaextractor
            fea = convert(Array{Float32}, fea)

            f[uttid, compress = args["compression-level"]] = fea
        end
    end
end

run(`hostname`)
println(ARGS)
args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
abstract type AbstractTextEncoder end

"""
    process(::AbstractTextEncoder)

Get processing function of given encoder.
"""
process(e::AbstractTextEncoder) = e.process
process(::Type{AbstractTextEncoder}) = nestedcall(getvalue)

"""
    tokenize(e::AbstractTextEncoder, x)

Use encoder's tokenizer to tokenize `x`.
"""
tokenize(e::AbstractTextEncoder, x)  = e.tokenizer(x)

"""
    process(e::AbstractTextEncoder, x)

Use encoder's processing function to process `x`.
"""
process(e::AbstractTextEncoder, x) = process(e)(x)

"""
    lookup(e::AbstractTextEncoder, x)

Lookup `x` in encoder's vocabulary.
"""
lookup(e::AbstractTextEncoder, x) = lookup(OneHot, e.vocab, x)

"""
    encode(e::AbstractTextEncoder, x)

Encode `x`.
"""
encode(e::AbstractTextEncoder, x) = lookup(e, process(e, tokenize(e, x)))

"""
    decode(e::AbstractTextEncoder, x)

Decode `x`.
"""
decode(e::AbstractTextEncoder, x) = lookup(eltype(e.vocab), e.vocab, x)

"""
    TextEncoder(tokenizer, vocab, process = nestedcall(getvalue))

A simple encoder implementation.
"""
struct TextEncoder{T<:AbstractTokenizer, V<:AbstractVocabulary, P} <: AbstractTextEncoder
    tokenizer::T
    vocab::V
    process::P
end

TextEncoder(tkr::AbstractTokenizer, vocab::AbstractVocabulary) = TextEncoder(tkr, vocab, process(AbstractTextEncoder))

TextEncoder(builder, tkr::AbstractTokenizer, vocab::AbstractVocabulary) = TextEncoder(builder, TextEncoder(tkr, vocab))

"""
    TextEncoder(builder, e::TextEncoder)

Given an encoder, return a new encoder that has the same tokenizer and vocabulary. `builder` is
 a function that take a encoder and return a new processing function.
"""
TextEncoder(builder, e::TextEncoder) = TextEncoder(e.tokenizer, e.vocab, builder(e))


text:
module Visualize_Text

using Modia3D

@assembly TextDemo begin
   font1 = Modia3D.Font(charSize=0.2)
   font2 = Modia3D.Font(fontFamily="Arial", bold=true, charSize=0.2, color="Red", transparency=0.0)
   text1 = Modia3D.TextShape("This is a text line parallel to the screen"; font=font1)
   text2 = Modia3D.TextShape("This is text in the XY-plane"; font=font2, offset=[0.2, 0.0, 0.0],
                             axisAlignment=Modia3D.XY_Plane, alignment=Modia3D.Left)
   text3 = Modia3D.TextShape("This is text in the XZ-plane"; font=font2, offset=[0.0, 0.0, 0.2],
                             axisAlignment=Modia3D.XZ_Plane, alignment=Modia3D.Left)
   text4 = Modia3D.TextShape("This is text in the YZ-plane"; font=font2, offset=[0.0, 0.2, 0.0],
                             axisAlignment=Modia3D.YZ_Plane, alignment=Modia3D.Left)

   world      = Modia3D.Object3D(visualizeFrame=false)
   textShape1 = Modia3D.Object3D(world, text1)
   textShape2 = Modia3D.Object3D(world, text2; r=[1.0, 1.0, -1.0])
   textShape3 = Modia3D.Object3D(world, text3; r=[1.0, 1.0, -1.0])
   textShape4 = Modia3D.Object3D(world, text4; r=[1.0, 1.0, -1.0])
end

Modia3D.visualizeAssembly!( TextDemo(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.5) ) )

println("... success of Visualize_Text.jl!")

end


text:
const FILENAMES = ("epsg",
                   "esri")

parse_projection(line::String) = match(r"<(\d+)>\s+(.*?)\s+<>", line)

println("# Contents of this file is generated on $(now()). Do not edit by hand!")
println("# \$ julia generate_projection_codes.jl > projection_codes.jl\n")

for filename in FILENAMES
    if !isfile(filename)
        url = "https://raw.githubusercontent.com/OSGeo/proj.4/master/nad/$filename"
        download(url, filename)
    end
    proj_file = open(filename)
    lines = readlines(proj_file)
    close(proj_file)

    println("$filename = Dict(")
    for i=1:length(lines)
        m = parse_projection(lines[i])
        if is(m, nothing)
            continue
        else
            proj_code, proj_string = m.captures
            @printf("  %s => \"%s\", %s", proj_code, proj_string, lines[i-1])
        end
    end
    println(")\n")
end


text:
module EverestCoreModeTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Core Mode (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCoreModeTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.EverestCoreModeTrigger)
    return Dict{String, Any}(
        "mode" => sort(Maple.core_modes)
    )
end

end


text:
# ------------------------------------------------------------------------------
# Computations involving triangles
# -------------------------------------------------------------------------------

using SymPy

"""
# Function

    circumcircle(points)

Find the circle that goes through a list of points.
"""
function circumcircle(points)
    ptnum = length(points)
    if ptnum <= 2
        throw(ArgumentError("At least 3 points are needed, but only $(ptnum) are given."))
    end

    x, y = @vars x y
    c = Point(x, y)
    dist = [squaredist(pt, c) for pt in points]
    equations = [Eq(dist[i], dist[i+1]) for i in 1:length(points)-1]
    sol = solve(equations, [x, y])
    center = Point(simplify(sol[x]), simplify(sol[y]))
    radius = distance(center, points[1])
    radius = simplify(radius)
    Circle(center, radius)
end


text: println("Vem pro Emacs Logo")


text:
import Flux
import LinearAlgebra

include("TensorOps.jl")

struct MultiHeadAttention{T}
    WQueries::AbstractArray{T, 3}
    WKeys::AbstractArray{T, 3}
    WValues::AbstractArray{T, 3}
    HeadDense::Flux.TrackedArray
    L::Integer
end

function MultiHeadAttention(dims::Integer, layers::Integer)
    OutDims = Int32(floor(dims / layers))
    MultiHeadAttention(
        cat([randn(Float32, dims, OutDims) for _ in 1:layers]..., dims=3) |> Flux.gpu |> Flux.param,
        cat([randn(Float32, dims, OutDims) for _ in 1:layers]..., dims=3) |> Flux.gpu |> Flux.param,
        cat([randn(Float32, dims, OutDims) for _ in 1:layers]..., dims=3) |> Flux.gpu |> Flux.param,
        randn(Float32, dims, dims) |> Flux.gpu |> Flux.param,
        layers
    )
end
Flux.@treelike MultiHeadAttention

function ScalarDotAttention(Query::AbstractArray{T, 3}, Key::AbstractArray{T, 3}, Value::AbstractArray{T, 3}; futuremask = false) where T
    Numerator = BatchMatMul(Query, permutedims(Key, [2, 1, 3]))
    Denominator = Float32(sqrt(size(Key, 1)))

    if futuremask
        Mask = ones(Float32, size(Numerator)[1:2])

        # lower triangular matrix with very large negative numbers
        Mask = LinearAlgebra.tril!(trues(size(Mask)), -1) * Float32(-1e9)
        Mask = Flux.param(Mask) |> Flux.gpu
        Numerator = Numerator .* Mask
    end

    x = TensorSoftmax(Numerator / Denominator, dims=2)
    return BatchMatMul(x, Value)
end

function _MHALayer(MHA::MultiHeadAttention, MHAData::AbstractArray{T, 1} where T, L::Integer; futuremask = false)
    (Query, Key, Value) = MHAData
    Nbatches = size(Query, 3)

    HQueries = TensorDot(Query, MHA.WQueries[:, :, L])
    HKeys = TensorDot(Key, MHA.WKeys[:, :, L])
    HValues = TensorDot(Value, MHA.WValues[:, :, L])

    return ScalarDotAttention(HQueries, HKeys, HValues, futuremask=futuremask)
end

function (m::MultiHeadAttention)(Query::AbstractArray{T, 3}, Key::AbstractArray{T, 3}, Value::AbstractArray{T, 3}; layers = 1, futuremask = false) where T
    TotalHeads = cat([_MHALayer(m, [Query, Key, Value], layer, futuremask=futuremask) for layer in 1:m.L]..., dims=2)
    return TensorDot(TotalHeads, m.HeadDense)
end


text:
module ExtraTrees

using Statistics, Distributed, SharedArrays

export Model, FlattenModel

include("./Tree.jl")
include("./Model.jl")
include("./FlattenModel.jl")

end # module


text:
module SaveVTK2

using StaticArrays
using WriteVTK
using Random
using Parameters

export write_vtp

@with_kw mutable struct SimData
    Points::Array{StaticArrays.SArray{Tuple{3},Float32,1,3},1} =
        Array{StaticArrays.SArray{Tuple{3},Float32,1,3},1}()
    Idp::Array{Int32,1} = Array{Int32,1}()
    Vel::Array{StaticArrays.SArray{Tuple{3},Float32,1,3},1} =
        Array{StaticArrays.SArray{Tuple{3},Float32,1,3},1}()
    Rhop::Array{Float32,1} = Array{Float32,1}()
end

"""
    write_vtp(filename,points,attribute)

Saves a polydata file (.vtp) given a filename, a 2xN point array and and 1xN attribute array.
"""
function write_vtp(filename::String, sim_arr::SimData)

    #points   = convert_to_array(sim_arr.Points)
    points = reinterpret(eltype(eltype(sim_arr.Points)), sim_arr.Points')

    polys = empty(MeshCell{WriteVTK.PolyData.Polys,UnitRange{Int64}}[])
    verts = empty(MeshCell{WriteVTK.PolyData.Verts,UnitRange{Int64}}[])

    # Note: the order of verts, lines, polys and strips is not important.
    # One doesn't even need to pass all of them.
    all_cells = (verts, polys)

    fname = filename
    #allocations here, not my library
    vtk = vtk_grid(fname, points, all_cells..., compress = true, append = false)

    #vtk[@Name(attribute)] = attribute;

    if !isempty(sim_arr.Idp)
        vtk["Idp"] = sim_arr.Idp
    end
    if !isempty(sim_arr.Vel)
        velocity = reinterpret(eltype(eltype(sim_arr.Vel)), sim_arr.Vel')
        vtk["Velocity"] = velocity
    end
    if !isempty(sim_arr.Rhop)
        vtk["Rhop"] = sim_arr.Rhop
    end

    vtk_save(vtk)
end

end


text:
# functions related to beta distributions

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    betacdf,
    betaccdf,
    betalogcdf,
    betalogccdf,
    betainvcdf,
    betainvccdf,
    betainvlogcdf,
    betainvlogccdf

# Julia implementations
betapdf(α::Real, β::Real, x::Real) = exp(betalogpdf(α, β, x))

betalogpdf(α::Real, β::Real, x::Real) = betalogpdf(promote(α, β, x)...)
function betalogpdf(α::T, β::T, x::T) where {T<:Real}
    # we ensure that `log(x)` and `log1p(-x)` do not error
    y = clamp(x, 0, 1)
    val = xlogy(α - 1, y) + xlog1py(β - 1, -y) - logbeta(α, β)
    return x < 0 || x > 1 ? oftype(val, -Inf) : val
end


text:
include("xor.jl")

@testset "XOR data generator" begin
    a = rand_bin(20)
    @test length(a) == 20
    @test all(a.<=1)
    @test all(a.>=0)
    @test all(typeof.(a) .== Float64)

    b = xor(a)
    @test all(typeof.(b) .== Float64)
    @test length(b) == 1

    X, y = xor_dataset(20, 100)
    @test length(X)==100
    @test length(y)==100
    @test length(X[1])==20
    @test all(typeof.(y[1]) .== Float64)

    @test all(xor.(X) .== y)
end

@testset "Log Loss" begin
    @test log_loss(0, 0) < 0.01
    @test log_loss(0, 1) > 0
    @test log_loss(1, 1) < 0.01
    @test log_loss(1, 0) > 0

    X, y = xor_dataset(20, 100)
    @test sum(log_loss.(y, y)) < 0.1
end


text:
using KernelAbstractions, Test, CUDA

@kernel function mul2(A)
    I = @index(Global)
    A[I] = 2 * A[I]
end

@kernel function add3(A, B, C)
    I = @index(Global)
    A[I] = B[I] + C[I]
end

@kernel function addi(A, C, i)
    I = @index(Global)
    A[I] = i + C[I]
end

function test_typed_kernel_dynamic()
    A = ones(1024, 1024)
    kernel = mul2(CPU())
    res = @ka_code_typed kernel(A, ndrange=size(A), workgroupsize=16)
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})

    if has_cuda_gpu()
        A = CUDA.ones(1024, 1024)
        kernel = mul2(CUDADevice())
        res = @ka_code_typed kernel(A, ndrange=size(A), workgroupsize=(32, 32))
        @test isa(res, Pair{Core.CodeInfo, DataType})
        @test isa(res[1].code, Array{Any,1})
    end
end

function test_typed_kernel_dynamic_no_info()
    A = ones(1024, 1024)
    B = similar(A)
    C = similar(A)
    kernel = add3(CPU())
    res = @ka_code_typed kernel(A, B, C, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})

    if has_cuda_gpu()
        A = CUDA.ones(1024, 1024)
        B = similar(A)
        C = similar(A)
        kernel = add3(CUDADevice())
        res = @ka_code_typed kernel(A, B, C, ndrange=size(A))
        @test isa(res, Pair{Core.CodeInfo, DataType})
        @test isa(res[1].code, Array{Any,1})
    end
end

function test_typed_kernel_static()
    A = ones(1024, 1024)
    kernel = mul2(CPU(), 16)
    res = @ka_code_typed kernel(A, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})

    if has_cuda_gpu()
        A = CUDA.ones(1024, 1024)
        kernel = mul2(CUDADevice(), (32, 32))
        res = @ka_code_typed kernel(A, ndrange=size(A))
        @test isa(res, Pair{Core.CodeInfo, DataType})
        @test isa(res[1].code, Array{Any,1})
    end
end

function test_typed_kernel_no_optimize()
    A = ones(1024, 1024)
    B = similar(A)
    C = similar(A)
    kernel = add3(CPU(), 16)
    res = @ka_code_typed optimize=false kernel(A, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo,Core.TypeofBottom})
    res_opt = @ka_code_typed kernel(A, ndrange=size(A))
    @test size(res[1].code) < size(res_opt[1].code)

    if has_cuda_gpu()
        A = CUDA.ones(1024, 1024)
        B = similar(A)
        C = similar(A)
        kernel = add3(CUDADevice(), (32, 32))
        res = @ka_code_typed optimize=false kernel(A, ndrange=size(A))
        @test isa(res, Pair{Core.CodeInfo,Core.TypeofBottom})
        res_opt = @ka_code_typed kernel(A, ndrange=size(A))
        @test size(res[1].code) < size(res_opt[1].code)
    end
end

function test_expr_kernel()
    A = ones(1024, 1024)
    C = similar(A)
    kernel = addi(CPU())
    res = @ka_code_typed kernel(A, C, 1+2, ndrange=size(A))
    @test isa(res, Pair{Core.CodeInfo, DataType})
    @test isa(res[1].code, Array{Any,1})

    if has_cuda_gpu()
        A = CUDA.ones(1024, 1024)
        C = similar(A)
        kernel = addi(CUDADevice(), (32, 32))
        res = @ka_code_typed kernel(A, C, 1+2, ndrange=size(A))
        @test isa(res, Pair{Core.CodeInfo, DataType})
        @test isa(res[1].code, Array{Any,1})
    end
end

test_typed_kernel_dynamic()
test_typed_kernel_dynamic_no_info()
test_typed_kernel_static()
test_typed_kernel_no_optimize()
test_expr_kernel()


text:
using BinDeps

ENV["JULIA_ROOT"] = abspath(JULIA_HOME, "../../")

cd(joinpath(Pkg.dir(), "Clang", "deps", "src") )
run(`make`)
if (!ispath("../usr"))
  run(`mkdir ../usr`)
end
if (!ispath("../usr/lib"))
  run(`mkdir ../usr/lib`)
end
run(`mv libwrapclang.$(BinDeps.shlib_ext) ../usr/lib`)