text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

mutable struct P4estMPICache{uEltype}
  mpi_neighbor_ranks::Vector{Int}
  mpi_neighbor_interfaces::Vector{Vector{Int}}
  mpi_send_buffers::Vector{Vector{uEltype}}
  mpi_recv_buffers::Vector{Vector{uEltype}}
  mpi_send_requests::Vector{MPI.Request}
  mpi_recv_requests::Vector{MPI.Request}
  n_elements_by_rank::OffsetArray{Int, 1, Array{Int, 1}}
  n_elements_global::Int
  first_element_global_id::Int
end

function P4estMPICache(uEltype)
  # MPI communication "just works" for bitstypes only
  if !isbitstype(uEltype)
    throw(ArgumentError("P4estMPICache only supports bitstypes, $uEltype is not a bitstype."))
  end

  mpi_neighbor_ranks = Vector{Int}(undef, 0)
  mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, 0)
  mpi_send_buffers = Vector{Vector{uEltype}}(undef, 0)
  mpi_recv_buffers = Vector{Vector{uEltype}}(undef, 0)
  mpi_send_requests = Vector{MPI.Request}(undef, 0)
  mpi_recv_requests = Vector{MPI.Request}(undef, 0)
  n_elements_by_rank = OffsetArray(Vector{Int}(undef, 0), 0:-1)
  n_elements_global = 0
  first_element_global_id = 0

  P4estMPICache{uEltype}(mpi_neighbor_ranks, mpi_neighbor_interfaces,
                         mpi_send_buffers, mpi_recv_buffers,
                         mpi_send_requests, mpi_recv_requests,
                         n_elements_by_rank, n_elements_global,
                         first_element_global_id)
end

function start_mpi_send!(mpi_cache::P4estMPICache, mesh, equations, dg, cache)
  data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1)

  for d in 1:length(mpi_cache.mpi_neighbor_ranks)
    send_buffer = mpi_cache.mpi_send_buffers[d]

    for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last  = (index - 1) * data_size + data_size
      local_side = cache.mpi_interfaces.local_sides[interface]
      @views send_buffer[first:last] .= vec(cache.mpi_interfaces.u[local_side, .., interface])
    end
  end

  for (index, d) in enumerate(mpi_cache.mpi_neighbor_ranks)
    mpi_cache.mpi_send_requests[index] = MPI.Isend(
      mpi_cache.mpi_send_buffers[index], d, mpi_rank(), mpi_comm())
  end

  return nothing
end

function start_mpi_receive!(mpi_cache::P4estMPICache)
  for (index, d) in enumerate(mpi_cache.mpi_neighbor_ranks)
    mpi_cache.mpi_recv_requests[index] = MPI.Irecv!(
      mpi_cache.mpi_recv_buffers[index], d, d, mpi_comm())
  end

  return nothing
end

function finish_mpi_send!(mpi_cache::P4estMPICache)
  MPI.Waitall!(mpi_cache.mpi_send_requests)
end

function finish_mpi_receive!(mpi_cache::P4estMPICache, mesh, equations, dg, cache)
  data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1)

  # Start receiving and unpack received data until all communication is finished
  d, _ = MPI.Waitany!(mpi_cache.mpi_recv_requests)
  while d != 0
    recv_buffer = mpi_cache.mpi_recv_buffers[d]

    for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last  = (index - 1) * data_size + data_size

      if cache.mpi_interfaces.local_sides[interface] == 1 # local element on primary side
        @views vec(cache.mpi_interfaces.u[2, .., interface]) .= recv_buffer[first:last]
      else # local element at secondary side
        @views vec(cache.mpi_interfaces.u[1, .., interface]) .= recv_buffer[first:last]
      end
    end

    d, _ = MPI.Waitany!(mpi_cache.mpi_recv_requests)
  end

  return nothing
end

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::ParallelP4estMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype<:Real}
  # Make sure to balance the p4est before creating any containers
  # in case someone has tampered with the p4est after creating the mesh
  balance!(mesh)

  elements       = init_elements(mesh, equations, dg.basis, uEltype)
  interfaces     = init_interfaces(mesh, equations, dg.basis, elements)
  mpi_interfaces = init_mpi_interfaces(mesh, equations, dg.basis, elements)
  boundaries     = init_boundaries(mesh, equations, dg.basis, elements)
  mortars        = init_mortars(mesh, equations, dg.basis, elements)
  mpi_cache      = init_mpi_cache(mesh, elements, mpi_interfaces,
                                  nvariables(equations), nnodes(dg), uEltype)

  cache = (; elements, interfaces, mpi_interfaces, boundaries, mortars, mpi_cache)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (; cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
  cache = (; cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

function init_mpi_cache(mesh::ParallelP4estMesh, elements, mpi_interfaces, nvars, nnodes, uEltype)
  mpi_cache = P4estMPICache(uEltype)
  init_mpi_cache!(mpi_cache, mesh, elements, mpi_interfaces, nvars, nnodes, uEltype)

  return mpi_cache
end

function init_mpi_cache!(mpi_cache::P4estMPICache, mesh::ParallelP4estMesh, elements, mpi_interfaces,
                         nvars, n_nodes, uEltype)
  mpi_neighbor_ranks, mpi_neighbor_interfaces =
    init_mpi_neighbor_connectivity(mpi_interfaces, mesh)

  mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests =
    init_mpi_data_structures(mpi_neighbor_interfaces, ndims(mesh), nvars, n_nodes, uEltype, nothing)

  # Determine local and total number of elements
  n_elements_global = Int(mesh.p4est.global_num_quadrants)
  n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks())),
                            n_elements_global) |> diff # diff sufficient due to 0-based quad indices
  n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1))
  # Account for 1-based indexing in Julia
  first_element_global_id = Int(unsafe_load(mesh.p4est.global_first_quadrant, mpi_rank() + 1)) + 1
  @assert n_elements_global == sum(n_elements_by_rank) "error in total number of elements"

  # TODO reuse existing structures
  @pack! mpi_cache = mpi_neighbor_ranks, mpi_neighbor_interfaces,
                     mpi_send_buffers, mpi_recv_buffers,
                     mpi_send_requests, mpi_recv_requests,
                     n_elements_by_rank, n_elements_global,
                     first_element_global_id

end

function init_mpi_neighbor_connectivity(mpi_interfaces, mesh::ParallelP4estMesh)
  # Let p4est iterate over all interfaces and call init_neighbor_rank_connectivity_iter_face
  # to collect connectivity information
  iter_face_c = cfunction(init_neighbor_rank_connectivity_iter_face, Val(ndims(mesh)))
  user_data = InitNeighborRankConnectivityIterFaceUserData(mpi_interfaces, mesh)

  # Ghost layer is required to determine owner ranks of quadrants on neighboring processes
  ghost_layer = ghost_new_p4est(mesh.p4est)
  @assert ghost_is_valid_p4est(mesh.p4est, ghost_layer) == 1
  iterate_p4est(mesh.p4est, user_data; ghost_layer=ghost_layer, iter_face_c=iter_face_c)
  ghost_destroy_p4est(ghost_layer)

  # Build proper connectivity data structures from information gathered by iterating over p4est
  @unpack global_interface_ids, neighbor_ranks_interface = user_data
  mpi_neighbor_ranks = neighbor_ranks_interface |> sort |> unique
  p = sortperm(global_interface_ids)
  neighbor_ranks_interface .= neighbor_ranks_interface[p]
  interface_ids = collect(1:nmpiinterfaces(mpi_interfaces))[p]

  # For each neighbor rank, init connectivity data structures
  mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, length(mpi_neighbor_ranks))
  for (index, d) in enumerate(mpi_neighbor_ranks)
    mpi_neighbor_interfaces[index] = interface_ids[findall(==(d), neighbor_ranks_interface)]
  end

  # Check that all interfaces were counted exactly once
  @assert sum(length(v) for v in mpi_neighbor_interfaces) == nmpiinterfaces(mpi_interfaces)

  return mpi_neighbor_ranks, mpi_neighbor_interfaces
end

mutable struct InitNeighborRankConnectivityIterFaceUserData{MPIInterfaces, Mesh}
  interfaces::MPIInterfaces
  interface_id::Int
  global_interface_ids::Vector{Int}
  neighbor_ranks_interface::Vector{Int}
  mesh::Mesh
end

function InitNeighborRankConnectivityIterFaceUserData(mpi_interfaces, mesh)
  global_interface_ids = fill(-1, nmpiinterfaces(mpi_interfaces))
  neighbor_ranks_interface = fill(-1, nmpiinterfaces(mpi_interfaces))

  return InitNeighborRankConnectivityIterFaceUserData{typeof(mpi_interfaces), typeof(mesh)}(
    mpi_interfaces, 1, global_interface_ids, neighbor_ranks_interface, mesh)
end

function init_neighbor_rank_connectivity_iter_face(info, user_data)
  data = unsafe_pointer_to_objref(Ptr{InitNeighborRankConnectivityIterFaceUserData}(user_data))

  # Function barrier because the unpacked user_data above is not type-stable
  init_neighbor_rank_connectivity_iter_face_inner(info, data)
end

# 2D
cfunction(::typeof(init_neighbor_rank_connectivity_iter_face), ::Val{2}) = @cfunction(init_neighbor_rank_connectivity_iter_face, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))
# 3D
cfunction(::typeof(init_neighbor_rank_connectivity_iter_face), ::Val{3}) = @cfunction(init_neighbor_rank_connectivity_iter_face, Cvoid, (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid}))

# Function barrier for type stability
function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
  @unpack interfaces, interface_id, global_interface_ids, neighbor_ranks_interface, mesh = user_data

  # Get the global interface ids and neighbor rank if current face belongs to an MPI interface
  if info.sides.elem_count == 2 # MPI interfaces have two neighboring elements
    # Extract surface data
    sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0 # No hanging nodes for MPI interfaces
      if sides[1].is.full.is_ghost == 1
        remote_side = 1
        local_side = 2
      elseif sides[2].is.full.is_ghost == 1
        remote_side = 2
        local_side = 1
      else # both sides are on this rank -> skip since it's a regular interface
        return nothing
      end

      # Sanity check, current face should belong to current MPI interface
      local_tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) # one-based indexing
      local_quad_id = local_tree.quadrants_offset + sides[local_side].is.full.quadid
      @assert interfaces.local_element_ids[interface_id] == local_quad_id + 1 # one-based indexing

      # Get neighbor ID from ghost layer
      proc_offsets = unsafe_wrap(Array, info.ghost_layer.proc_offsets, mpi_nranks() + 1)
      ghost_id = sides[remote_side].is.full.quadid # indexes the ghost layer, 0-based
      neighbor_rank = findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r+1],
                                1:mpi_nranks()) - 1 # mpi ranks are 0-based
      neighbor_ranks_interface[interface_id] = neighbor_rank

      # Global interface id is the globally unique quadrant id of the quadrant on the primary
      # side (1) multiplied by the number of faces per quadrant plus face
      if local_side == 1
        offset = unsafe_load(mesh.p4est.global_first_quadrant, mpi_rank() + 1) # one-based indexing
        primary_quad_id = offset + local_quad_id
      else
        offset = unsafe_load(mesh.p4est.global_first_quadrant, neighbor_rank + 1) # one-based indexing
        primary_quad_id = offset + sides[1].is.full.quad.p.piggy3.local_num
      end
      global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides[1].face
      global_interface_ids[interface_id] = global_interface_id

      user_data.interface_id += 1
    end
  end

  return nothing
end

# TODO: ::Any is a temporary hack. This method is already defined in the TreeMesh
# code. It uses `mpi_neighbor_mortars` as an argument but not uEltype, thus the number
# of arguments would be the same without ::Any which would break precompilation.
# See https://github.com/trixi-framework/Trixi.jl/pull/977#discussion_r793694635 for further
# discussion.
function init_mpi_data_structures(mpi_neighbor_interfaces, n_dims, nvars, n_nodes, uEltype, ::Any)
  data_size = nvars * n_nodes^(n_dims - 1)
  mpi_send_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_buffers = Vector{Vector{uEltype}}(undef, length(mpi_neighbor_interfaces))
  for index in 1:length(mpi_neighbor_interfaces)
    mpi_send_buffers[index] = Vector{uEltype}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
    mpi_recv_buffers[index] = Vector{uEltype}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
  end

  mpi_send_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))

  return mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests
end

include("dg_2d_parallel.jl")
include("dg_3d_parallel.jl")

end # muladd


text:
@testset "Exp Atoms: $solver" for solver in solvers
    if can_solve_exp(solver)
        @testset "exp atom" begin
            y = Variable()
            p = Problem{BigFloat}(:minimize, exp(y), y>=0)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 1 atol=TOL
            @test evaluate(exp(y)) ≈ 1 atol=TOL

            y = Variable()
            p = Problem{BigFloat}(:minimize, exp(y), y>=1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ exp(1) atol=TOL
            @test evaluate(exp(y)) ≈ exp(1) atol=TOL

            y = Variable(5)
            p = Problem{BigFloat}(:minimize, sum(exp(y)), y>=0)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 5 atol=TOL
            @test evaluate(sum(exp(y))) ≈ 5 atol=TOL

            y = Variable(5)
            p = Problem{BigFloat}(:minimize, sum(exp(y)), y>=0)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 5 atol=TOL
        end

        @testset "log atom" begin
            y = Variable()
            p = Problem{BigFloat}(:maximize, log(y), y<=1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0 atol=TOL

            y = Variable()
            p = Problem{BigFloat}(:maximize, log(y), y<=2)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ log(2) atol=TOL

            y = Variable()
            p = Problem{BigFloat}(:maximize, log(y), [y<=2, exp(y)<=10])
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ log(2) atol=TOL
        end

        @testset "log sum exp atom" begin
            y = Variable(5)
            p = Problem{BigFloat}(:minimize, logsumexp(y), y>=1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ log(exp(1) * 5) atol=TOL
        end

        @testset "logistic loss atom" begin
            y = Variable(5)
            p = Problem{BigFloat}(:minimize, logisticloss(y), y>=1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ log(exp(1) + 1) * 5 atol=TOL
        end

        @testset "entropy atom" begin
            y = Variable(5, Positive())
            p = Problem{BigFloat}(:maximize, entropy(y), sum(y)<=1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ -(log(1 / 5)) atol=TOL
        end

        @testset "relative entropy atom" begin
            x = Variable(1)
            y = Variable(1)
            # x log (x/y)
            p = Problem{BigFloat}(:minimize, relative_entropy(x,y), y==1, x >= 2)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 2 * log(2) atol=TOL
        end

        @testset "log perspective atom" begin
            x = Variable(1)
            y = Variable(1)
            # y log (x/y)
            p = Problem{BigFloat}(:maximize, log_perspective(x,y), y==5, x <= 10)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 5 * log(2) atol=TOL
        end
    end
end

ext: jl


text:
module DPMeansClustering

	using Reexport

    @reexport using StatsBase
    @reexport using Distances

	include("dp.jl")

    export dp_centers, dp_means, cluster, get_initial_centroids, get_initial_metacentroids, get_possible_centroids, C2MC_dp_centers
end


text:
module Each

export each

"""
    each(Type, collection)

Return the collection unchanged.

    each(Type, element::Type)

Return the element wrapped in a tuple.

Useful for functions where the arguments can be a list or a single element which should be
iterated over, but the elements could be iterable themselves.
Extra care should be taken when the list and the elements are similar types: When working
with tuples of tuples of integers, use `each(Tuple{Vararg{Integer}},collection)` rather than
`each(Tuple,collection)` as the latter will always use the entire list as one element.

# Examples
```julia-repl
julia>  function f(a::Union{String,AbstractArray{String}})
            # This will fail for "scalar" Strings
            for s in a
                println(s)
            end
        end
julia>  f(["apples","oranges"])
apples
oranges
julia>  f("apples")
a
p
p
l
e
s
julia>  function f(a::Union{String,AbstractArray{String}})
            # This handles strings correctly
            for s in each(String,a)
                println(s)
            end
        end
julia>  f(["apples","oranges"])
apples
oranges
julia>  f("apples")
apples

julia>  function f(a::Tuple)
            for t in each(Tuple,a)
                println(t)
            end
        end
julia>  f( (1, 2) ) # argument is a tuple, wrap in tuple
(1, 2)
julia>  f( ( (1, 2), (3, 4) ) ) # same thing
((1, 2), (3, 4))
julia>  function f(a::Tuple)
            for t in each(Tuple{Vararg{Number}},a)
                println(t)
            end
        end
julia>  f( (1, 2) )
(1, 2)
julia>  f( ( (1,2), (3, 4) ) )
(1, 2)
(3, 4)
```
"""
function each(::Type{T}, collection) where T
    eltype(collection) <:T && return collection
    error(typeof(collection)," does not contain elements of type ",T)
end

each(::Type{T}, element::T) where T = (element,)

end # Module


text:
import RayCaster as RC
import Test

Test.@testset "RayCaster.jl" begin
    obstacle_tile_map = BitArray([
                                  1 1 1 1 1
                                  1 0 0 0 1
                                  1 0 0 0 1
                                  1 0 0 0 1
                                  1 1 1 1 1
                                 ])

    height_obstacle_tile_map = size(obstacle_tile_map, 1)
    width_obstacle_tile_map = size(obstacle_tile_map, 2)

    Test.@testset "Continuous world" begin
        T1 = Float32
        T2 = Float32

        x_start = convert(T1, 5 / 2)
        y_start = convert(T1, 5 / 2)

        Test.@testset "theta = 0" begin
            theta = convert(T2, 0)
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 5
            Test.@test j_stop == 3
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ 1.5
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = atan(2 / 3)" begin
            theta = convert(T2, atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 5
            Test.@test j_stop == 4
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = (pi / 2) - atan(2 / 3)" begin
            theta = convert(T2, (pi / 2) - atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 4
            Test.@test j_stop == 5
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = pi / 2" begin
            theta = convert(T2, pi / 2)
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 3
            Test.@test j_stop == 5
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ 1.5
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = (pi / 2) + atan(2 / 3)" begin
            theta = convert(T2, (pi / 2) + atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 2
            Test.@test j_stop == 5
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = pi - atan(2 / 3)" begin
            theta = convert(T2, pi - atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 1
            Test.@test j_stop == 4
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = pi" begin
            theta = convert(T2, pi)
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 1
            Test.@test j_stop == 3
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ 1.5
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = pi + atan(2 / 3)" begin
            theta = convert(T2, pi + atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 1
            Test.@test j_stop == 2
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = (3 * pi / 2) - atan(2 / 3)" begin
            theta = convert(T2, (3 * pi / 2) - atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 2
            Test.@test j_stop == 1
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = 3 * pi / 2" begin
            theta = convert(T2, 3 * pi / 2)
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 3
            Test.@test j_stop == 1
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ 1.5
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = (3 * pi / 2) + atan(2 / 3)" begin
            theta = convert(T2, (3 * pi / 2) + atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 4
            Test.@test j_stop == 1
            Test.@test hit_dimension == 2
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end

        Test.@testset "theta = 2 * pi - atan(2 / 3)" begin
            theta = convert(T2, 2 * pi - atan(2 / 3))
            cos_theta = cos(theta)
            sin_theta = sin(theta)
            i_stop, j_stop, hit_dimension, delta_x_world_units_to_exit_start_tile, delta_y_world_units_to_exit_start_tile = RC.cast_ray_continous_world(obstacle_tile_map, x_start, y_start, cos_theta, sin_theta)
            Test.@test i_stop == 5
            Test.@test j_stop == 2
            Test.@test hit_dimension == 1
            # Test.@test total_euclidean ≈ sqrt(1.5 ^ 2 + 1 ^ 2)
            Test.@test delta_x_world_units_to_exit_start_tile ≈ 0.5
            Test.@test delta_y_world_units_to_exit_start_tile ≈ 0.5
        end
    end

    Test.@testset "Discrete world" begin
        I = Int
        world_units_per_tile_unit = convert(I, 256)

        i_start_world_units = convert(I, height_obstacle_tile_map * world_units_per_tile_unit ÷ 2 + 1)
        j_start_world_units = convert(I, width_obstacle_tile_map * world_units_per_tile_unit ÷ 2 + 1)

        Test.@testset "delta_i = 1, delta_j = 0" begin
            delta_i_world_units = convert(I, 1)
            delta_j_world_units = convert(I, 0)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 5)
            Test.@test j_stop_tile_units == convert(I, 3)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = 2, delta_j = 1" begin
            delta_i_world_units = convert(I, 2)
            delta_j_world_units = convert(I, 1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 5)
            Test.@test j_stop_tile_units == convert(I, 4)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = 1, delta_j = 1" begin
            delta_i_world_units = convert(I, 1)
            delta_j_world_units = convert(I, 1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 5)
            Test.@test j_stop_tile_units == convert(I, 4)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = 1, delta_j = 2" begin
            delta_i_world_units = convert(I, 1)
            delta_j_world_units = convert(I, 2)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 4)
            Test.@test j_stop_tile_units == convert(I, 5)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = 0, delta_j = 1" begin
            delta_i_world_units = convert(I, 0)
            delta_j_world_units = convert(I, 1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 3)
            Test.@test j_stop_tile_units == convert(I, 5)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = -1, delta_j = 2" begin
            delta_i_world_units = convert(I, -1)
            delta_j_world_units = convert(I, 2)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 2)
            Test.@test j_stop_tile_units == convert(I, 5)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = -1, delta_j = 1" begin
            delta_i_world_units = convert(I, -1)
            delta_j_world_units = convert(I, 1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 2)
            Test.@test j_stop_tile_units == convert(I, 5)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = -2, delta_j = 1" begin
            delta_i_world_units = convert(I, -2)
            delta_j_world_units = convert(I, 1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 1)
            Test.@test j_stop_tile_units == convert(I, 4)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = -1, delta_j = 0" begin
            delta_i_world_units = convert(I, -1)
            delta_j_world_units = convert(I, 0)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 1)
            Test.@test j_stop_tile_units == convert(I, 3)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
        end

        Test.@testset "delta_i = -2, delta_j = -1" begin
            delta_i_world_units = convert(I, -2)
            delta_j_world_units = convert(I, -1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 1)
            Test.@test j_stop_tile_units == convert(I, 2)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = -1, delta_j = -1" begin
            delta_i_world_units = convert(I, -1)
            delta_j_world_units = convert(I, -1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 1)
            Test.@test j_stop_tile_units == convert(I, 2)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = -1, delta_j = -2" begin
            delta_i_world_units = convert(I, -1)
            delta_j_world_units = convert(I, -2)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 2)
            Test.@test j_stop_tile_units == convert(I, 1)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = 0, delta_j = -1" begin
            delta_i_world_units = convert(I, 0)
            delta_j_world_units = convert(I, -1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 3)
            Test.@test j_stop_tile_units == convert(I, 1)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = 1, delta_j = -2" begin
            delta_i_world_units = convert(I, 1)
            delta_j_world_units = convert(I, -2)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 4)
            Test.@test j_stop_tile_units == convert(I, 1)
            Test.@test hit_dimension == 2
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = 1, delta_j = -1" begin
            delta_i_world_units = convert(I, 1)
            delta_j_world_units = convert(I, -1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 5)
            Test.@test j_stop_tile_units == convert(I, 2)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end

        Test.@testset "delta_i = 2, delta_j = -1" begin
            delta_i_world_units = convert(I, 2)
            delta_j_world_units = convert(I, -1)
            i_stop_tile_units, j_stop_tile_units, hit_dimension, delta_i_world_units_to_exit_start_tile, delta_j_world_units_to_exit_start_tile = RC.cast_ray_discrete_world(obstacle_tile_map, i_start_world_units, j_start_world_units, delta_i_world_units, delta_j_world_units, world_units_per_tile_unit)
            Test.@test i_stop_tile_units == convert(I, 5)
            Test.@test j_stop_tile_units == convert(I, 2)
            Test.@test hit_dimension == 1
            Test.@test delta_i_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2
            Test.@test delta_j_world_units_to_exit_start_tile == world_units_per_tile_unit ÷ 2 + 1
        end
    end
end


text:
using GLFW
using VulkanCore
using VulkanCore.LibVulkan

@assert GLFW.VulkanSupported()

include("vkhelper.jl")

const WIDTH = 800
const HEIGHT = 600

## init GLFW window
GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)    # not to create an OpenGL context
GLFW.WindowHint(GLFW.RESIZABLE, 0)
window = GLFW.CreateWindow(WIDTH, HEIGHT, "Vulkan")

## init Vulkan
## creating instance
appInfoRef = VkApplicationInfo(
    "Application Name: Create Instance",
    v"1.0.0",
    "No Engine Name",
    v"1.0.0",
    VK_API_VERSION_1_2,
) |> Ref

extensions = GLFW.GetRequiredInstanceExtensions()
@test check_extensions(extensions)

layers = get(ENV, "JULIA_GITHUB_ACTIONS_CI", "OFF") == "ON" ? String[] : ["VK_LAYER_KHRONOS_validation"]
@test check_layers(layers)

createInfoRef = VkInstanceCreateInfo(appInfoRef, layers, extensions) |> Ref

instanceRef = Ref(VkInstance(C_NULL))
result = GC.@preserve appInfoRef layers extensions vkCreateInstance(createInfoRef, C_NULL, instanceRef)

@static if get(ENV, "JULIA_GITHUB_ACTIONS_CI", "OFF") == "ON"
    @test_broken result == VK_SUCCESS
else
    @test result == VK_SUCCESS
end

## cleaning up
if result == VK_SUCCESS
    vkDestroyInstance(instanceRef[], C_NULL)
end
GLFW.DestroyWindow(window)


text:
using ReTest
using FluxTraining
using Animations
using Colors
using ImageIO
using FluxTraining: protect, Events, Phases, runstep, runepoch, epoch!, step!, testlearner, SanityCheckException
using FluxTraining: getcallback, setcallbacks!, replacecallback!, removecallback!, addcallback!
using .Phases: Phase, AbstractTrainingPhase

using Flux: trainable
using Flux.Optimise: Descent
using Suppressor
using FluxTraining: CHECKS
using Zygote


text:
using Random
Random.seed!(Random.GLOBAL_RNG, 123)

using Turing

include("posterior_predictive.jl")
include("bayesplot.jl")
include("arviz.jl")

using Dates
using Feather
const pest_data = let data = Feather.read(joinpath(datadir(), "pest_data.feather"))
    data.date = convert.(Date, data.date);
    DataFrame((n => Array(data[:, n]) for n in names(data))...)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using RCall
macro sym_str(x) :(Symbol($(esc(x)))) end
@rimport binom as binom
@rimport exactci as exactci

# %%
using Distributions
using StatsPlots
using Roots
x ⪅ y = x < y || x ≈ y

pval_exact(dist, k) = sum(pdf(dist, j) for j in support(dist) if pdf(dist, j) ⪅ pdf(dist, k))
pval_exact(n, p, k) = pval_exact(Binomial(n, p), k)
pval_2x1sd(dist, k) = min(1, 2cdf(dist, k), 2ccdf(dist, k - 1))
pval_2x1sd(n, p, k) = pval_2x1sd(Binomial(n, p), k)
chisq_stat(n, p, k) = (k - n*p)^2/(n*p*(1 - p))
pval_chisq(n, p, k) = ccdf(Chisq(1), chisq_stat(n, p, k))
ci(pval, n, k, α=0.05) = find_zeros(p -> pval(n, p, k) - α, 0, 1)

# %%
n, k = 20, 3

# %%
binom.binom_confint(k, n)

# %%
@show exactci.binom_exact(k, n, tsmethod = "central")[sym"conf.int"][1:2]
@show exactci.binom_exact(k, n, tsmethod = "minlik")[sym"conf.int"][1:2]
@show exactci.binom_exact(k, n, tsmethod = "blaker")[sym"conf.int"][1:2];

# %%
@show ci(pval_exact, n, k)
@show ci(pval_2x1sd, n, k)
@show ci(pval_chisq, n, k);

# %%
n, k = 20, 3
plot(; xtick=0:0.1:1)
plot!(p -> pval_exact(n, p, k), 0, 1; label="exact two-sided (or minlik)")
plot!(p -> pval_2x1sd(n, p, k), 0, 1; label="exact doubled one-sided (or central)")
plot!(p -> pval_chisq(n, p, k), 0, 1; label="chi-squared (or Wilson)")

# %%

max_stars_repo_path: 0013/binomial confidence intervals.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/binomial confidence intervals.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/binomial confidence intervals.jl

max_forks_repo_name: genkuroki/public


text:
using RecipesBase

function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray = nothing) where PT
    hull = PT[]
    if length(ps) == 0
        return hull
    end
    prev = sign_sense == 1 ? first(ps) : last(ps)
    cur = prev
    for j in (sign_sense == 1 ? (2:length(ps)) : ((length(ps)-1):-1:1))
        while prev != cur && counterclockwise(cur - prev, ps[j] - prev) >= 0
            cur = prev
            pop!(hull)
            if !isempty(hull)
                prev = last(hull)
            end
        end
        if yray !== nothing && counterclockwise(ps[j] - cur, yray) >= 0
            break
        else
            push!(hull, cur)
            prev = cur
            cur = ps[j]
        end
    end
    push!(hull, cur)
    hull
end

@recipe function f(p::Polyhedron{2})
    removevredundancy!(p)
    if hasrays(p)
        warn("Rays not supported yet in the 2D plotting recipe")
    end
    ps = collect(points(p))
    sort!(ps, by = x -> x[1])
    counterclockwise(p1, p2) = dot(cross([p1; 0], [p2; 0]), [0, 0, 1])
    hull = [getsemihull(ps, 1, counterclockwise); getsemihull(ps, -1, counterclockwise)]
    seriestype --> :shape
    legend --> false
    [p[1] for p in hull], [p[2] for p in hull]
end


text:
# Fonction qu'on veut tester
# include("../NLPAtXmod.jl")
# importall NLPAtXmod

# On vérifie que le constructeur pour problème sans contrainte fonctionne
uncons_nlp_at_x = NLPAtX(zeros(10))

@test uncons_nlp_at_x.x == zeros(10)
@test uncons_nlp_at_x.fx == nothing
@test uncons_nlp_at_x.gx == nothing
@test uncons_nlp_at_x.Hx == nothing
@test uncons_nlp_at_x.mu == nothing
@test uncons_nlp_at_x.cx == nothing
@test uncons_nlp_at_x.Jx == nothing

@test uncons_nlp_at_x.lambda == zeros(0)
@test uncons_nlp_at_x.current_time == nothing
@test uncons_nlp_at_x.current_score == nothing
@test (!(uncons_nlp_at_x.evals == nothing))

# On vérifie que le constucteur pour problème avec contrainte fonctionne
cons_nlp_at_x = NLPAtX(zeros(10), zeros(10))

@test cons_nlp_at_x.x == zeros(10)
@test cons_nlp_at_x.fx == nothing
@test cons_nlp_at_x.gx == nothing
@test cons_nlp_at_x.Hx == nothing
@test cons_nlp_at_x.mu == nothing
@test cons_nlp_at_x.cx == nothing
@test cons_nlp_at_x.Jx == nothing
@test (false in (cons_nlp_at_x.lambda .== 0.0)) == false
@test cons_nlp_at_x.current_time == nothing
@test cons_nlp_at_x.current_score == nothing

# On vérifie que la fonction update! fonctionne
update!(uncons_nlp_at_x, x = ones(10), fx = 1.0, gx = ones(10))
update!(uncons_nlp_at_x, lambda = ones(10), current_time = 1.0)
update!(uncons_nlp_at_x, Hx = ones(10,10), mu = ones(10), cx = ones(10), Jx = ones(10,10))

@test (false in (uncons_nlp_at_x.x .== 1.0)) == false #assez bizarre comme test...
@test uncons_nlp_at_x.fx == 1.0
@test (false in (uncons_nlp_at_x.gx .== 1.0)) == false
@test (false in (uncons_nlp_at_x.Hx .== 1.0)) == false
@test uncons_nlp_at_x.mu == ones(10)
@test uncons_nlp_at_x.cx == ones(10)
@test (false in (uncons_nlp_at_x.Jx .== 1.0)) == false
@test (false in (uncons_nlp_at_x.lambda .== 1.0)) == false
@test uncons_nlp_at_x.current_time == 1.0
@test uncons_nlp_at_x.current_score == nothing

reinit!(uncons_nlp_at_x)
@test uncons_nlp_at_x.x == ones(10)
@test uncons_nlp_at_x.fx == nothing
reinit!(uncons_nlp_at_x, x = zeros(10))
@test uncons_nlp_at_x.x == zeros(10)
@test uncons_nlp_at_x.fx == nothing

nlp_64 = NLPAtX(ones(10))
nlp_64.x = ones(10)
nlp_64.fx = 1.0
nlp_64.gx = ones(10)

# nlp_32 = convert_nlp(Float32, nlp_64)
# @test typeof(nlp_32.x[1]) == Float32
# @test typeof(nlp_32.fx[1]) == Float32
# @test typeof(nlp_32.gx[1]) == Float32
# @test isnan(nlp_32.mu[1])
# @test isnan(nlp_32.current_time)
#
# @test typeof(nlp_64.x[1]) == Float64
# @test typeof(nlp_64.fx[1]) == Float64
# @test typeof(nlp_64.gx[1]) == Float64
# @test isnan(nlp_64.mu[1])
# @test isnan(nlp_64.current_time)

#Test the _size_check:
try
    NLPAtX(ones(5), gx = zeros(4))
    @test false
catch
    @test true
end
try
    NLPAtX(ones(5), mu = zeros(4))
    @test false
catch
    @test true
end
try
    NLPAtX(ones(5), Hx = zeros(4,4))
    @test false
catch
    @test true
end
try
    NLPAtX(ones(5), zeros(1), cx = zeros(0))
    @test false
catch
    @test true
end


text:
module LyceumBase

using Adapt
using AutoHashEquals

using Base: @propagate_inbounds
using Base.Threads: Atomic, atomic_add!, atomic_sub!

using Dates
using Distributions: Distributions, Sampleable, sample, sample!
using DocStringExtensions
using ElasticArrays
using Future: randjump
using IterTools
using LinearAlgebra
using MacroTools
using Pkg
using Printf
using Random
using Shapes

using SpecialArrays
using SpecialArrays: True, False

using StaticArrays
using StructArrays
using UnicodePlots: UnicodePlots
using UnsafeArrays
using UnPack

include("util.jl")

####
#### Interfaces
####

"""
    tconstruct(T::Type, n::Integer, args...; kwargs...)

Return a `AbstractVector` of `n` instances of `T`. By default, this returns
`[T(args...; kwargs...) for _=1:n], but this function can be extended for types that can share
data across instances for greater cache efficiency/performance.
"""
function tconstruct(T::Type, n::Integer, args...; kwargs...)
    n > 0 || throw(ArgumentError("n must be > 0"))
    ntuple(_ -> T(args...; kwargs...), n)
end
export tconstruct

export AbstractEnvironment
export statespace, getstate!, getstate, setstate!
export observationspace, getobservation!, getobservation
export actionspace, getaction!, getaction, setaction!
export rewardspace, getreward
export reset!, randreset!
export step!, isdone, timestep
include("abstractenvironment.jl")

####
#### Tools
####

include("setfield.jl")
using .SetfieldImpl: @set!

include("math.jl")

export tseed!, getrange, splitrange, nblasthreads, @with_blasthreads
include("threading.jl")

export SPoint3D, MPoint3D
include("geometry.jl")

export Trajectory, TrajectoryBuffer, ntrajectories, nsamples, rollout!, rollout
include("trajectory.jl")

export EnvironmentSampler
include("environmentsampler.jl")

include("plotting.jl")

end # module


text:
# Script to get a feeling for speed vs accuracy of quadrature schemes
using QuadGK,FastGaussQuadrature
using Plots

# Get the "truth" function
# bare FD distribution (for decoupled neutrinos the form of the bare FD doesn't change with time)
# The analytic value of the density integral at early times is 7/8 π²/15 * Tν^4
# FD function multiplied by a reasonable pert (use the massive nu ICs)
# For "truth" use high-resolution integration
# Choose qmin, qmax, n_q

# Compare tolerances

#One-liner to map finite interval to unit interval (and back)
to_ui(lq,lqmi,lqma) = -1 + (1- (-1)) / (lqma-lqmi) * (lq-lqmi)
from_ui(x,lqmi,lqma) = lqmi + (lqma- lqmi) / (1- (-1)) * (x- (-1))

#how to choose this seems a bit up in the air, but CAMB uses qmax of 30T
T=2e-4
logqmin,logqmax=log10(T/30),log10(T*30)#1e-6,1e-1
dxdq(q) = (1+to_ui(1+logqmin,logqmin,logqmax))/(q*log(10)) #chovar factor, dx/dq=dx/dlogq dlogq/dq

plot()
f0(q) = 4π *2(2π)^-3 *q^2 / (1 + exp(q/T))
ζ=1.202056903159594
answer = 3 *T^3 *ζ / (2π^2)
logqq=collect(logqmin:.01:logqmax)
plot(logqq,f0.(10.0 .^logqq))
plot(to_ui.(logqq,logqmin,logqmax),
     f0.(10.0 .^from_ui.(to_ui.(logqq,logqmin,logqmax),logqmin,logqmax)))
ft(x) = f0(10.0 ^ from_ui(x,logqmin,logqmax)) / dxdq(10.0 ^ from_ui(x,logqmin,logqmax))
plot!(to_ui.(logqq,logqmin,logqmax),ft.(to_ui.(logqq,logqmin,logqmax)))

#QuadGK
# The lazy thing I was doing before - just run quadgk at low tolerance
@time Iqgk = quadgk(q ->  f0(q) ,10.0 .^logqmin, 10.0 .^logqmax,rtol=1e-6)[1]
println("Relative error on quadGK: ", abs.(1-Iqgk/answer))
#this is already very slow - it takes .73 seconds at rtol=1e-2
#going to 1e-6 (which seems about what fast quadrature is getting), gives ~.26 s
#Should get average timings b/c seeing big difference from first call
#Comparison with FastGaussQuadrature should take into account saved quantities.

#FastGaussQuadrature.jl routines
quad_sum(x,w,f) = sum(f.(x).*w)

# Legendre quadrature
n=10
@time xGLe, wGLe = gausslegendre( n )
@time IGLe = quad_sum(xGLe,wGLe,ft)
println("Relative error on GLeg: ", abs.(1-IGLe/answer), "  with n= ", n," points.")
#just from changing n it looks like rel error bottoms out at 1e-11
#^I guess this is the error we induce by not integrating the tails - truly negligible for bg...
#with 10 pts using the T/30 30T range get 1e-3 error, will do this for testing

# CC quadrature
ncc = n+3
@time xGCC, wGCC = gausslobatto( ncc )
@time IGCC = quad_sum(xGCC,wGCC,ft)
println("Relative error on CC: ", abs.(1-IGCC/answer), "  with n= ", ncc," points.")
#30 pts is pretty good - relative tolerance of 1e-6 for Leg, 1e-4 for cheby zeros

#FIXME: try using asymptotic expansion instead of just truncating the integral (standard practice)

#FIXME: Need to take into account possible integration over splines instead of f
#If integrate over spline is faster but probably introduces interp error, especially at ends
#Also we won't spline the perts - but should check this though for bg

# More realistic test on perts
#Look at a the result of the un-normalized density integral at z=0 using dlnf0dlnq to test q-dep
#"truth" as a higher-order low-tolerance quadgk over larger q-range
dlnf0dlnq0(q) = -q / T /(1 + exp(-q/T))

#high-tolerance QuadGK
a=1#10^-4
@time Iqgk_truth = quadgk(q ->  √(q^2 + a^2 * 0.06^2)* q^2 * f0(q)*dlnf0dlnq0(q) ,10.0 .^logqmin, 10.0 .^logqmax,rtol=1e-12)[1]

#fast QuadGK at 1e-3
@time Iqgk_fast = quadgk(q ->  √(q^2 + a^2 * 0.06^2)* q^2 * f0(q)*dlnf0dlnq0(q) ,10.0 .^logqmin, 10.0 .^logqmax,rtol=1e-3)[1]
println("Relative error on realistic 1e-3 gk: ", abs.(1-Iqgk_fast/Iqgk_truth))

#vs FastGaussQuadrature
function ftp(x)
    return (√((10.0 ^ from_ui(x,logqmin,logqmax))^2 + a^2 * 0.06^2)
            * f0(10.0 ^ from_ui(x,logqmin,logqmax)) * (10.0 ^ from_ui(x,logqmin,logqmax))^2
            *dlnf0dlnq0(10.0 ^ from_ui(x,logqmin,logqmax)) )/ dxdq(10.0 ^ from_ui(x,logqmin,logqmax))
end

#Gauss-Legendre
m = 15
@time xGLe2, wGLe2 = gausslegendre( m )
@time IGLe2 = quad_sum(xGLe2,wGLe2,ftp)
println("Relative error on realistic GLeg: ", abs.(1-IGLe2/Iqgk_truth), "  with m= ", m," points.")

#CC
mcc = m#+1
@time xGCC2, wGCC2 = gausslobatto( mcc )
@time IGCC2 = quad_sum(xGCC2,wGCC2,ftp)
println("Relative error on CC: ", abs.(1-IGCC2/Iqgk_truth), "  with n= ", mcc," points.")

#comparing to the f0 test, where we got 1e-3 error for 10 and 13 pts for GL and CC respectively, here we have
#12 and 13 pts, respectively, and similarly, we have that it is faster
#Need to do a careful comparison at some point where we can average over a number of times,
#as it is hard to know when the quad points are reusing p/re-computed values...
#Either way so long as the computation of quadrature points (which occurs in the background once)
#plus the sum operations that happen frequently are faster than quadgk * number of times it is called,
#then we are fine - and I think this is an easily attainable threshold for a small number of quad pts (<15)

#Would be interesting to see if cosmology affects this accuracy, because could reduce cost even further by
#just passing the q-pts (computed externally) once instead of computing in bg

#hmm I forgot a factor of q before, but even at higher redshift the performance
#is about the same on 15 pts, only around 1e-3 rerr


text:
using Gnuplot
x = LinRange(-10,10,200)
fg(x,μ,σ) = exp.(.-(x.-μ).^2 ./(2σ^2))./(σ*√(2π))
@gp x fg(x, 0.25, 1.5) "w filledcu lc '#E69F00' dt 1 t '0.25,1.5'"
@gp :- "set style fill transparent solid 0.3 noborder"
@gp :- "set key title 'μ,σ' box 3" xlab = "x" ylab = "P(x)"
@gp :- x fg(x, 2, 1) "w filledcu lc '#56B4E9' dt 1 t '2,1'"
@gp :- x fg(x, -1, 2) "w filledcu lc '#009E73' dt 1 t '-1,2'"
save(term="pngcairo size 600,400", output="./code/plt_ex13.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using MAT
export
save,
load,
Diary,
scalar,
writepb,
psave,
pload,
activate,
logging

pybytes(b) = PyObject(
    ccall(PyCall.@pysym(PyCall.PyString_FromStringAndSize),
                                             PyPtr, (Ptr{UInt8}, Int),
                                             b, sizeof(b)))

"""
    psave(o::PyObject, file::String)

Saves a Python objection `o` to `file`.
See also [`pload`](@ref)
"""
function psave(o::PyObject, file::String)
    pickle = pyimport("pickle")
    f = open(file, "w")
    pickle.dump(o, f)
    close(f)
end

"""
    pload(file::String)

Loads a Python objection from `file`.
See also [`psave`](@ref)
"""
function pload(file::String)
    r = nothing
    pickle = pyimport("pickle")
    @pywith pybuiltin("open")(file,"rb") as f begin
        r = pickle.load(f)
    end
    return r
end

"""
    save(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)

Saves the values of `vars` in the session `sess`. The result is written into `file` as a dictionary. If `vars` is nothing, it saves all the trainable variables.
See also [`save`](@ref), [`load`](@ref)
"""
function save(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject},  Array{Any}}=nothing, args...; kwargs...)
    if vars==nothing
        vars = get_collection(TRAINABLE_VARIABLES)
    elseif isa(vars, PyObject)
        vars = [vars]
    elseif isa(vars, Array{Any})
        vars = Array{PyObject}(vars)
    end
    d = Dict{String, Any}()
    vals = run(sess, vars, args...;kwargs...)
    for i = 1:length(vars)
        name = replace(vars[i].name, ":"=>"colon")
        name = replace(name, "/"=>"backslash")
        d[name] = vals[i]
    end
    matwrite(file, d)
    d
end

function save(sess::PyObject, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)
    if vars==nothing
        vars = get_collection(TRAINABLE_VARIABLES)
    elseif isa(vars, PyObject)
        vars = [vars]
    end
    d = Dict{String, Any}()
    vals = run(sess, vars, args...;kwargs...)
    for i = 1:length(vars)
        name = replace(vars[i].name, ":"=>"colon")
        name = replace(name, "/"=>"backslash")
        d[name] = vals[i]
    end
    d
end

"""
    load(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)

Loads the values of variables to the session `sess` from the file `file`. If `vars` is nothing, it loads values to all the trainable variables.
See also [`save`](@ref), [`load`](@ref)
"""
function load(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)
    if vars==nothing
        vars = get_collection(TRAINABLE_VARIABLES)
    elseif isa(vars, PyObject)
        vars = [vars]
    end
    d = matread(file)
    ops = PyObject[]
    for i = 1:length(vars)
        name = replace(vars[i].name, ":"=>"colon")
        name = replace(name,  "/"=>"backslash")
        if !(name in keys(d))
            @warn "$(vars[i].name) not found in the file, skipped"
        else
            if occursin("bias", name) && isa(d[name], Number)
                d[name] = [d[name]]
            end
            push!(ops, assign(vars[i], d[name]))
        end
    end
    run(sess, ops, args...; kwargs...)
end

function load(sess::PyObject, d::Dict, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)
    if vars==nothing
        vars = get_collection(TRAINABLE_VARIABLES)
    elseif isa(vars, PyObject)
        vars = [vars]
    end
    ops = PyObject[]
    for i = 1:length(vars)
        name = replace(vars[i].name, ":"=>"colon")
        name = replace(name,  "/"=>"backslash")
        if !(name in keys(d))
            @warn "$(vars[i].name) not found in the file, skipped"
        else
            push!(ops, assign(vars[i], d[name]))
        end
    end
    run(sess, ops, args...; kwargs...)
end

mutable struct Diary
    writer::PyObject
    tdir::String
end

"""
    Diary(suffix::Union{String, Nothing}=nothing)

Creates a diary at a temporary directory path. It returns a writer and the corresponding directory path
"""
function Diary(suffix::Union{String, Nothing}=nothing)
    tdir = mktempdir()
    printstyled("tensorboard --logdir=\"$(tdir)\" --port 0\n", color=:blue)
    Diary(tf.summary.FileWriter(tdir, tf.get_default_graph(),filename_suffix=suffix), tdir)
end

"""
    save(sw::Diary, dirp::String)

Saves [`Diary`](@ref) to `dirp`.
"""
function save(sw::Diary, dirp::String)
    cp(sw.tdir, dirp, force=true)
end

"""
    load(sw::Diary, dirp::String)

Loads [`Diary`](@ref) from `dirp`.
"""
function load(sw::Diary,dirp::String)
    sw.writer = tf.summary.FileWriter(dirp, tf.get_default_graph())
    sw.tdir = dirp
end

"""
    activate(sw::Diary, port::Int64=6006)

Running [`Diary`](@ref) at http://localhost:port.
"""
function activate(sw::Diary, port::Int64=0)
    printstyled("tensorboard --logdir=\"$(sw.tdir)\" --port $port\n", color=:blue)
    run(`tensorboard --logdir="$(sw.tdir)" --port $port '&'`)
end

"""
    scalar(o::PyObject, name::String)

Returns a scalar summary object.
"""
function scalar(o::PyObject, name::Union{String,Missing}=missing)
    if ismissing(name)
        name = tensorname(o)
    end
    tf.summary.scalar(name, o)
end

"""
    write(sw::Diary, step::Int64, cnt::Union{String, Array{String}})

Writes to [`Diary`](@ref).
"""
function Base.:write(sw::Diary, step::Int64, cnt::Union{String, Array{String}})
    if isa(cnt, String)
        sw.writer.add_summary(pybytes(cnt), step)
    else
        for c in cnt
            sw.writer.add_summary(pybytes(c), step)
        end
    end
end

function writepb(writer::PyObject, sess::PyObject)
    output_path = writer.get_logdir()
    tf.io.write_graph(sess.graph_def, output_path, "model.pb")
    return
end

"""
    logging(file::Union{Nothing,String}, o::PyObject...; summarize::Int64 = 3, sep::String = " ")

Logging `o` to `file`. This operator must be used with [`bind`](@ref).
"""
function logging(file::Union{Nothing,String}, o::PyObject...; summarize::Int64 = 3, sep::String = " ")
    if isnothing(file)
        tf.print(o..., summarize=summarize, sep=sep)
    else
        filepath = "file://$(abspath(file))"
        tf.print(o..., output_stream=filepath, summarize=summarize, sep=sep)
    end
end

logging(o::PyObject...; summarize::Int64 = 3, sep::String = " ") = logging(nothing, summarize=summarize, sep=sep)


text:
module ZonalFlow

using OrdinaryDiffEq
using DiffEqCallbacks
using RecursiveArrayTools
using FFTW
using LinearAlgebra

include("coefficients.jl")
include("equations.jl")
include("ics.jl")
include("tools.jl")
include("solvers.jl")
include("analysis.jl")

export nl,gql,gce2,ic_pert_eqm,ic_eqm
export energy,zonalenergy,inversefourier,meanvorticity,zonalvelocity,fourierenergy

end


text:
module markup

import Term: replace_double_brackets, OPEN_TAG_REGEX, replace_text
import Term: get_last_valid_str_idx, CLOSE_TAG_REGEX, GENERIC_CLOSER_REGEX

export extract_markup, MarkupTag, pairup_tags, clean_nested_tags

# -------------------------------- single tag -------------------------------- #
"""
    SingleTag

Represents a single tag `[style]` or `[/style]`
"""
mutable struct SingleTag
    markup::String
    start::Int
    stop::Int
    isclose::Bool
end

"""
    SingleTag(match::RegexMatch)

Construct a `SingleTag` out of a `RegexMatch`
"""
function SingleTag(match::RegexMatch)
    _start = match.offset
    _stop = match.offset + length(match.match) - 1

    return SingleTag(match.match[2:(end - 1)], _start, _stop, match.match[2] == '/')
end

# -------------------------------- markup tag -------------------------------- #
"""
    MarkupTag

Represents a complete markup tag.

It stores two `SingleTag`, the text inbetween and any other `MarkupTag`
that was detected in that text.
"""
mutable struct MarkupTag
    markup::String
    open::SingleTag
    close::SingleTag
    text::String
    inner_tags::Vector
end
function MarkupTag(text::AbstractString, open::SingleTag, close::SingleTag, inner::Vector)
    return MarkupTag(open.markup, open, close, text, inner)
end

function Base.show(io::IO, tag::MarkupTag)
    return print(io, "Markup tag ($(length(tag.inner_tags)) inner tags)")
end

# ---------------------------- extract markup tags --------------------------- #

"""
    has_markup(text::AbstractString)

Returns `true` if `text` includes a `MarkupTag`
"""
has_markup(text::AbstractString) = occursin(OPEN_TAG_REGEX, replace_double_brackets(text))

"""
    extract_markup(input_text::AbstractString; firstonly=false)

Extracts `MarkupTag`s from a piece of text.
"""
function extract_markup(input_text::AbstractString; firstonly = false)
    text = input_text  # copy so that we can edit it
    tags = []

    while has_markup(input_text)
        # for open_match in eachmatch(OPEN_TAG_REGEX, text)
        open_match = match(OPEN_TAG_REGEX, input_text)

        # get tag declaration
        tag_open = SingleTag(open_match)

        # get tag closing
        close_regex = r"\[\/+" * tag_open.markup * r"\]"
        if !occursin(close_regex, text[(tag_open.stop):end])
            # check if a generig closer [/] is present
            if occursin(GENERIC_CLOSER_REGEX, text[(tag_open.stop):end])
                tag_close = SingleTag(collect(eachmatch(GENERIC_CLOSER_REGEX, text))[end])
            else
                # otherwise inject a closer tag at the end
                text = text * "[/$(tag_open.markup)]"
                input_text = input_text * "[/$(tag_open.markup)]"
                tag_close = SingleTag(match(close_regex, text, tag_open.start))
            end
        else
            tag_close = SingleTag(match(close_regex, text, tag_open.start))
        end

        # crate tag and keep
        _start = tag_open.stop + 1
        _stop = tag_close.start - 1

        _start = !isvalid(text, _start) ? get_next_valid_str_idx(text, _start) : _start
        _stop = !isvalid(text, _stop) ? get_last_valid_str_idx(text, _stop) : _stop
        contained = text[_start:_stop]

        tag_open.stop = _start - 1
        tag_close.start = _stop + 1

        # get nested tags
        inner_tags = extract_markup(input_text[(tag_open.stop + 1):(tag_close.start - 1)])

        # remove tag's inside text
        input_text = replace_text(input_text, tag_open.start - 1, tag_close.stop + 1)

        # @info "TAG" input_text  text contained length(contained)
        push!(tags, MarkupTag(contained, tag_open, tag_close, inner_tags))
        if firstonly
            return tags[1]
        end
    end

    return tags
end

# ----------------------------- helper functions ----------------------------- #

"""
    clean_nested_tags(text::AbstractString)::AbstractString

Given a text with nested string like:
`[red]aaaa [green]bbbb[/green] cccc [blue] ddddd [/blue]eeee[/red]`

it adds extra tags to ensure that text within inner tags is handled properly, giving:
`[red]aaaa [green]bbbb[/green][red] cccc [/red][blue] ddddd [/blue][red]eeee[/red]`
"""
function clean_nested_tags(text::AbstractString)
    if !has_markup(text)
        return text
    end

    # @info "cleaning nested tags" text
    tags = extract_markup(text)
    for tag in tags
        if length(tag.inner_tags) > 0
            text = clean_nested_tags(tag, text)
        end
    end
    # text = join_lines(clean_nested_tags.(, text))
    # @info "after cleaning" text
    return text
end

"""
    clean_nested_tags(tag, text::AbstractString)

recursively applies to inner tags
"""
function clean_nested_tags(tag::MarkupTag, text::AbstractString)
    # @info "     inner cleaning nested tags" text
    for inner in tag.inner_tags
        text = clean_nested_tags(inner, text)

        _tag = "[$(inner.open.markup)]$(inner.text)[$(inner.close.markup)]"
        replacement = "[$(tag.close.markup)]" * _tag * "[$(tag.open.markup)]"
        text = replace(text, _tag => replacement)
    end

    return text
end

"""
    pairup_tags(text::Vector{AbstractString})

Given a vector of string with markup tags not properly closed/opened across lines,
it fixes things up.
"""
function pairup_tags(text::Vector{AbstractString})::Vector{AbstractString}
    # sweep lines to add tags starts/end
    for i in length(text):-1:1
        # get all open tags
        for opener in reverse(collect(eachmatch(OPEN_TAG_REGEX, text[i])))
            tag_open = SingleTag(opener)

            # get line where it closes
            j = i
            close_regex = r"\[\/+" * tag_open.markup * r"\]"
            for closing_line in text[i:end]
                occursin(close_regex, closing_line) ? break : j += 1
            end

            # if they are in the same line, continue
            if j == i
                continue
            end

            # if no close tag was found, add it at the end
            if j > length(text)
                text[end] = text[end] * "[/$(tag_open.markup)]"
                j = length(text)
            end

            # add correct close/open markups
            for idx in i:(j - 1)
                text[idx] = text[idx] * "[/$(tag_open.markup)]"
                text[1 + idx] = "[$(tag_open.markup)]" * text[1 + idx]
            end
        end
    end

    # do another sweep to close up orphaned tags
    for (i, line) in enumerate(text)
        for match in eachmatch(CLOSE_TAG_REGEX, line)
            markup = match.match[3:(end - 1)]

            open_regex = r"\[" * markup * r"\]"
            if !occursin(open_regex, line[1:(match.offset)])
                line = "[$markup]" * line
                # @info "injected" markup i line
            end
        end
        text[i] = line
    end
    return text
end
end


text:
function runExperimentAndWriteResults(kappaSquared::Array{Float64},h::Array{Float64,1},src::Array{Int64,1},n::Array{Int64,1},T_exact::Array{Float64},WU::Float64 )
Omega = zeros(2*length(n));
Omega[2:2:end] = (n.-1).*h;
Mesh = getRegularMesh(Omega,n.-1);
##########################################################

# if length(n)==2
	# matshow(sqrt(kappaSquared)); colorbar();
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],0:1:4);
# end

pMem = getEikonalTempMemory(n);
pEik = getEikonalParam(Mesh,kappaSquared,src,false);
pEik.T1 = zeros(Float64,tuple(n...));

tictimes = 1000000;
times = div(tictimes,prod(n)) .+ 1;
println("Timing is based on averaged ",times," runs.");
tt = time_ns()
for k=1:times
	solveFastMarchingUpwindGrad(pEik,pMem);
end
t1 = (time_ns() - tt)/(10^9)
t1 = t1/times;

if Mesh.dim==2
	selfMultiplyWithAnalyticSolution2D(Mesh.n.+1,Mesh.h,src,pEik.T1);
else
	selfMultiplyWithAnalyticSolution3D(Mesh.n.+1,Mesh.h,src,pEik.T1)
end

T1 = copy(pEik.T1)

pEik.HO = true;
ERR1 = T_exact .- pEik.T1;

tt = time_ns()
for k=1:times
	solveFastMarchingUpwindGrad(pEik,pMem);
end
t2 = (time_ns() - tt)/(10^9)
t2 = t2/times;

if Mesh.dim==2
	selfMultiplyWithAnalyticSolution2D(Mesh.n.+1,Mesh.h,src,pEik.T1);
else
	selfMultiplyWithAnalyticSolution3D(Mesh.n.+1,Mesh.h,src,pEik.T1)
end

sqrtN = sqrt(prod(n));
T2 = copy(pEik.T1)

ERR2 = T_exact - pEik.T1;

@printf("[%3.4f,%3.4f]\t[%d,%d]\t[%3.2e,%3.2e]\t%3.2fs(%3.2f)\t",h[1],h[2],n[1],n[2],norm(ERR1,Inf),norm(ERR1,2)/sqrtN,t1,t1/WU);

@printf("[%3.2e,%3.2e]\t%3.2fs(%3.2f)\n",norm(ERR2,Inf),norm(ERR2,2)/sqrtN,t2,t2/WU);

# if length(n)==2
	# figure()
	# contour(T_exact[end:-1:1,:],20,linestyles = "-",colors = "black",linewidths=1.5)
	# contour(T1[end:-1:1,:],20,linestyles = ":",colors = "red", linewidths=1.5)
	# contour(T2[end:-1:1,:],20,linestyles = "--",colors = "blue", linewidths=1.5)
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],4:-1:0);

	# matshow(abs(T1 - T_exact)); colorbar();
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],0:1:4);

	# matshow(abs(T2 - T_exact)); colorbar();
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],8):n[2],0:1:8);
	# yticks(0:div(n[1],4):n[1],0:1:4);
# end

# if length(n)==2
	# figure()
	# contour(T_exact[end:-1:1,:],100,linestyles = "-",colors = "black",linewidths=2.0)
	# contour(T1[end:-1:1,:],100,linestyles = ":",colors = "red", linewidths=2.0)
	# contour(T2[end:-1:1,:],100,linestyles = "--",colors = "blue", linewidths=2.0)
	# xlabel("y");
	# ylabel("x");
	# xticks(0:div(n[2],16):n[2],0:0.5:8);
	# yticks(0:div(n[1],40):n[1],4:-0.1:0);

# end

return;
end


text:
function constantdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000, dev=CPU())
   τ = 1/kappa  # time-scale for diffusive decay
  dt = 1e-9 * τ # dynamics are resolved

  prob = Diffusion.Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper, dev=dev)
  g = prob.grid

  # a gaussian initial condition c(x, t=0)
  c0ampl, σ = 0.01, 0.2
  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
  c0 = c0func.(g.x)

  # analytic solution for for 1D heat equation with constant κ
  tfinal = nsteps*dt
  σt = sqrt(2*kappa*tfinal + σ^2)
  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))

  Diffusion.set_c!(prob, c0)
  tcomp = @elapsed stepforward!(prob, nsteps)
  Diffusion.updatevars!(prob)

  prob, c0, cfinal, nsteps, tcomp
end

function constantdiffusionproblem2D(stepper; nx=128, ny=64, Lx=2π, Ly=4π, kappa=1e-2, nsteps=1000, dev=CPU())
    τ = 1/kappa  # time-scale for diffusive decay
    dt = 1e-9 * τ # dynamics are resolved

prob = Diffusion2D.Problem(nx=nx, ny=ny, Lx=Lx, Ly=Ly, kappa=kappa, dt=dt, stepper=stepper, dev=dev)
g = prob.grid

#Gaussian Initial Condition
c0, σ0 = 0.1, 0.2
c1, σ1 = 0.1, 0.2

c0_x_func(x) = @. c0*exp(-x^2/(2*(σ0)^2))
c0_y_func(y) = @. c1*exp(-y^2/(2*(σ1)^2))

c0_val = c0_x_func.(g.x) .* c0_y_func.(g.y)

#Analytic Solution
tfinal = nsteps*dt

σ0t = sqrt(2*kappa*tfinal + (σ0)^2)
σ1t = sqrt(2*kappa*tfinal + (σ1)^2)

cfinal_x = @. c0 * (σ0/σ0t) * exp(-g.x^2/(2*(σ0t)^2))
cfinal_y = @. c1 * (σ1/σ1t) * exp(-g.y^2/(2*(σ1t)^2))

cfinal_val = cfinal_x .* cfinal_y

#Numerical Solution
Diffusion2D.set_c!(prob, c0_val)
tcomp = @elapsed stepforward!(prob, nsteps)
Diffusion2D.updatevars!(prob)

prob, c0_val, cfinal_val, nsteps, tcomp
end

function varyingdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000, dev=CPU())
   τ = 1/kappa  # time-scale for diffusive decay
  dt = 1e-9 * τ # dynamics are resolved

  kappa = kappa*ones(nx) # this is actually a constant diffusion but defining it
                         # as an array makes stepforward! call function calcN!
                         # instead of just the linear coefficients L*sol

  prob = Diffusion.Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper, dev=dev)
  g = prob.grid

  # a gaussian initial condition c(x, t=0)
  c0ampl, σ = 0.01, 0.2
  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
  c0 = c0func.(g.x)

  # analytic solution for for 1D heat equation with constant κ
  tfinal = nsteps*dt
  σt = sqrt(2*kappa[1]*tfinal + σ^2)
  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))

  Diffusion.set_c!(prob, c0)
  tcomp = @elapsed stepforward!(prob, nsteps)
  Diffusion.updatevars!(prob)

  prob, c0, cfinal, nsteps, tcomp
end

function constantdiffusiontest(stepper, dev::Device=CPU(); kwargs...)
  prob, c0, c1, nsteps, tcomp = constantdiffusionproblem(stepper; kwargs...)
  normmsg = "$stepper: relative error ="
  @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
  isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)
end

function constantdiffusiontest2D(stepper, dev::Device=CPU(); kwargs...)
    prob, ci, cf, nsteps, tcomp = constantdiffusionproblem2D(stepper; kwargs...)
    normmsg = "$stepper: relative error ="
    @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(cf-prob.vars.c)/norm(cf), tcomp)
    isapprox(cf, prob.vars.c, rtol=nsteps*rtol_timesteppers)
end

function varyingdiffusiontest(stepper, dev::Device=CPU(); kwargs...)
  prob, c0, c1, nsteps, tcomp = varyingdiffusionproblem(stepper; kwargs...)
  normmsg = "$stepper: relative error ="
  @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
  isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)
end


text:
function Nn_series_coeff!(t::Transit_Struct{T}) where {T <: Real}
# Use series expansion to compute N_n.
# Computing leading coefficient (n=0).  Need to compute for
# n_max-3 to n_max.  Also, need just for k^2 < 1
# Nn_coeff is 2-dimensional array:
# Nn_coeff[2,jmax], where first dimension is for n_max-1 (1) to n_max (2),
# while third dimension is for the series coefficients.
# Computing leading coefficient (n=0):
coeff = zero(T)
# Loop over n_max-1 to n_max:
for j=1:2
  # m is n_max-1 to n_max (we need both since downward recursion requires two):
  m = t.n_max+j-2
  mhalf = 0.5*m
  coeff = 0.5*sqrt(pi)*exp(logabsgamma(mhalf+1.0)[1]-logabsgamma(mhalf+2.5)[1])
  # Add leading term to N_n:
  t.Nn_coeff[j,1] = coeff
  # Now, compute higher order terms until desired precision is reached:
  @inbounds for i=1:t.jmax-1
    coeff *= (4*i^2-1)/(2*i*(3+m+2*i))
    t.Nn_coeff[j,i+1] = coeff
  end
end
end


text:
# Command: Take Screenshot
"""
    screenshot(session::Session)::String
    screenshot(element::Element)::String

Take Screenshot.
"""
function screenshot end

function screenshot(session::Session)::String
    @unpack addr, id = session
    response =
        HTTP.get("$addr/session/$id/screenshot", [("Content-Type" => "application/json")])
    @assert response.status == 200
    JSON3.read(response.body).value
end
function screenshot(element::Element)::String
    @unpack addr, id = element.session
    response = HTTP.get(
        "$addr/session/$id/element/$(element.id)/screenshot",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
# Utilities for extracting information from relocation entries
r_sym(rel::Union(ELF32.Rel,ELF32.Rela)) = rel.r_info >> 8
r_type(rel::Union(ELF32.Rel,ELF32.Rela)) = rel.r_info % UInt8
r_sym(rel::Union(ELF64.Rel,ELF64.Rela)) = rel.r_info >> 32
r_type(rel::Union(ELF64.Rel,ELF64.Rela)) = rel.r_info % Uint32
addend(rel::Union(ELF32.Rela,ELF64.Rela)) = rel.r_addend
addend(rel::Union(ELF32.Rel,ELF64.Rel)) = 0

r_sym(x) = r_sym(deref(x))
r_type(x) = r_type(deref(x))
addend(x) = addend(deref(x))

# The relocation to apply
immutable RelocToApply
    value::UInt64
    size::Uint8
end

function compute_value(h, symbols, reloc)
    # Here one could to external symbol lookups, etc, but I'm not interested
    deref(symbols[r_sym(reloc) + 1]).st_value
end

# Apply relocations in `buffer`. `h` should be the buffer being relocated
function relocate!(buffer, h; LOI = nothing)
    sects = Sections(h)
    rel_sects = filter(sects) do sec
        sht = deref(sec).sh_type
        sht == SHT_RELA || sht == SHT_REL
    end
    for sec in rel_sects
        symbols = Symbols(link_sec(sec))
        relocated_sec = info_sec(sec)
        for reloc in Relocations(sec)
            if r_type(reloc) == 0
                continue
            end
            Value = compute_value(h, symbols, reloc)
            if LOI === nothing
                SectionLoadAddress = deref(relocated_sec).sh_addr
            else
                SectionLoadAddress = getSectionLoadAddress(LOI,relocated_sec)
            end
            rta = compute_relocation(h, reloc, SectionLoadAddress, Value)
            seek(buffer, deref(relocated_sec).sh_offset + deref(reloc).r_offset)
            write(buffer, rta.size == 8 ? rta.value :
                rta.size == 4 ? convert(UInt32,rta.value) :
                rta.size == 2 ? convert(UInt16,rta.value) :
                rta.size == 1 ? convert(UInt8,rta.value) :
                error("Unsupported Relocation Size"))
        end
    end
end

function compute_relocation(h, reloc, SectionLoadAddress, Value)
    if h.file.header.e_machine == EM_X86_64
        return compute_X86_64(reloc, SectionLoadAddress, Value)
    else
        error("Unknown Machine!")
    end
end

function compute_X86_64(reloc, SectionLoadAddress, Value)
    # Any error generated by the conversions below indicates a bad
    # value (the conversions are checked)
    kind = r_type(reloc)
    Offset = deref(reloc).r_offset
    Value += addend(reloc)
    if kind == R_X86_64_64
        return RelocToApply(Value, 8)
    elseif kind == R_X86_64_32
        return RelocToApply(convert(Uint32,Value), 4)
    elseif kind == R_X86_64_32S
        return RelocToApply(convert(Int32,reinterpret(Int,Value)), 4)
    else
        PCOffset = Value - (SectionLoadAddress + Offset)
        if kind == R_X86_64_PC64
            return RelocToApply(PCOffset, 8)
        elseif kind == R_X86_64_PC32
            return RelocToApply(convert(Int32,PCOffset), 4)
        else
            error("Unknown relocation ($(R_X86_64[kind])")
        end
    end
end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
@testset "create_surfaces_neuron.jl" begin

end


text:
# # Structural Model
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/ISS.ipynb)
#
#md # !!! note "Overview"
#md #     System type: linear time-invariant continuous system\
#md #     State dimension: 270\
#md #     Application domain: Aerospace Engineering
#
# ## Model description
#
# This is a model of component 1r (Russian service module) of the International Space Station [^ASG00].
# It has 270 states, 3 inputs and 3 outputs  The model consists of a continuous
# linear time-invariant system
#
# ```math
#   \begin{array}{lcl}
#   \dot{x}(t) &=& Ax(t) + Bu(t) \\
#   y(t) & = & C x(t)
#   \end{array}
# ```
# It was proposed as a benchmark in ARCH 2016 [^TLT16]. The matrix dimensiones are
# ``A ∈ \mathbb{R}^{270\times 270}``, ``B ∈ \mathbb{R}^{270\times 3}`` and
# ``C ∈ \mathbb{R}^{3\times 270}``.

# The matrices A, B, and C are available in MATLAB format
# ![](slicot.org/objects/software/shared/bench-data/iss.zip) (`iss.mat`). For convenience
# such `.mat` file has been converted to the [JLD2 format](https://github.com/JuliaIO/JLD2.jl)
# and stored in the file `iss.jld2`.

using ReachabilityAnalysis, JLD2
using ReachabilityAnalysis: add_dimension

LazySets.set_ztol(Float64, 1e-15);

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
ISS_path = joinpath(examples_dir, "ISS", "ISS.jld2")

@load ISS_path C;
const C3 = C[3, :]; # variable y₃
const C3_ext = vcat(C3, fill(0.0, 3));

# ## Reachability settings
#
# Initially all the variables are in the range ``[-0.0001, 0.0001]``, and the inputs
# are bounded: ``u_1(t) ∈ [0, 0.1]``, ``u_2(t) ∈ [0.8, 1]``, and ``u_3`` is in ``[0.9, 1]``.
# The time bound is 20.

# There are two versions of this benchmark, ISSF01 (time-varying inputs) and
# ISSC01 (constant inputs).

# - **ISSF01** (time-varying inputs): In this setting, the inputs can change
#   arbitrarily over time: ``\forall t: u(t) \in \mathcal{U}``.
#
# - **ISSC01** (constant inputs): The inputs are uncertain only in their initial
#   value, and constant over time: ``u(0)\in \mathcal{U}``, ``\dot u(t)= 0``.

# ## Specifications

# The verification goal is to check the ranges reachable by the output ``y_3(t)``,
# which is a linear combination of the state variables. In addition to the safety specification,
# for each version there is an UNSAT instance that serves as sanity checks to ensure
# that the model and the tool work as intended. But there is a caveat: In principle,
# verifying an UNSAT instance only makes sense formally if a witness is provided
# (counter-example, under-approximation, etc.).

# Since this library currently doesn't have the capability of verifying
# UNSAT instances, we run the tool with the same accuracy settings on an SAT-UNSAT pair of instances.
# The SAT instance demonstrates that the over-approximation is not too coarse,
# and the UNSAT instance demonstrates that the over-approximation is indeed conservative,
# at least in the narrow sense of the specification.

# - ISS01: Bounded time, safe property: For all ``t ∈ [0, 20]``, ``y_3(t) ∈ [−0.0007, 0.0007]``.
#          This property is used with the uncertain input case (ISSF01) and assumed to be satisfied.
#
# - ISS02: Bounded time, safe property: For all ``t ∈ [0, 20]``, ``y_3(t) ∈ [−0.0005, 0.0005]``.
#          This property is used with the constant input case (ISSC01) and assumed to be satisfied.
#
# - ISU01: Bounded time, unsafe property: For all ``t ∈ [0, 20]``, ``y_3 (t) ∈ [−0.0005, 0.0005]``.
#          This property is used with the uncertain input case (ISSF01) and assumed to be unsatisfied.
#
# - ISU02: Bounded time, unsafe property: For all ``t ∈ [0, 20]``, ``y_3 (t) ∈ [−0.00017, 0.00017]``.
#          This property is used with the constant input case (ISSC01) and assumed to be unsatisfied.

# ## Results

# The specification involves only the output ``y(t) := C_3 x(t)``, where ``C_3`` denotes
# the third row of the output matrix ``C ∈ \mathbb{R}^{3\times 270}``. Hence,
# it is sufficient to compute the flowpipe associated to ``y(t)`` directly,
# without the need of actually computing the full 270-dimensional flowpipe associated to
# all state variables. The flowpipe associated to a linear combination of state
# variables can be computed efficiently using the support-function based algorithm
# [[LGG09]](@ref). The idea is to define a template polyhedron with only two
# supporting directions, namely ``C_3`` and ``-C_3``.

# The step sizes chosen are ``6×10^{-4}`` for ISSF01 and ``1×10^{-2}`` for ISSC01.

# ### ISSF01

# In this scenario the inputs can change arbitrarily over time and are only bound
# to the given input range ``U``.

function ISSF01()
    @load ISS_path A B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
    X0 = BallInf(zeros(size(A, 1)), 0.0001)
    return @ivp(x' = A*x + B*u, x(0) ∈ X0, u ∈ U, x ∈ Universe(270))
end

#-

dirs = CustomDirections([C3, -C3]);
prob_ISSF01 = ISSF01();
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false));

# The solution `sol_ISSF01` is now a 270-dimensional set that (only) contains
# template reach-sets for the linear combination `C_3 x(t)`,

dim(sol_ISSF01)

# For visualization, it is necessary to specify that we want to plot "time" vs. $y_3(t)$.
# We can transform the flowpipe on the output ``y_3(t)`` by "flattening" the flowpipe
# along directions `dirs`.

πsol_ISSF01 = flatten(sol_ISSF01);

# Now `πsol_ISSF01` is a one-dimensional flowpipe,

dim(πsol_ISSF01)

#md # !!! tip "Performance tip"
#md #      Note that projecting the solution along direction ``C_3`` corresponds to computing
#md #      the min and max bounds for each reach-set `X`, that is, `Interval(-ρ(-C3, X), ρ(C3, X))`.
#md #      However, the method `flatten(sol_ISSF01, rows=(1, 2))` is more efficient since it uses
#md #      the matrix of support function evaluations obtained by `LGG09` along directions
#md #      ``C_3`` and ``-C_3``.

using Plots, Plots.PlotMeasures, LaTeXStrings

fig = Plots.plot();
Plots.plot!(fig, πsol_ISSF01[1:10:end], vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8,
     xlab=L"t",
     ylab=L"y_{3}",
     xtick=[0, 5, 10, 15, 20.], ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075],
     xlims=(0., 20.), ylims=(-0.00075, 0.00075),
     bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm);
fig

# ### ISSC01

# In this scenario the inputs are uncertain only in their initial
# value and constant over time.

function ISSC01()
    @load ISS_path A B

    A_ext = add_dimension(A, 3)
    A_ext[1:270, 271:273] = B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
    X0 = BallInf(zeros(size(A, 1)), 0.0001)
    X0 = X0 * U
    return @ivp(x' = A_ext*x, x(0) ∈ X0)
end

#-

dirs = CustomDirections([C3_ext, -C3_ext]);
prob_ISSC01 = ISSC01();
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false));

#-

# We can flatten the flowpipe to the output ``y_3(t)`` as before:

πsol_ISSC01 = flatten(sol_ISSC01);

fig = Plots.plot();
Plots.plot!(fig, πsol_ISSC01, vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8, lw=1.0,
    xlab=L"t",
    ylab=L"y_{3}",
    xtick=[0, 5, 10, 15, 20.], ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002],
    xlims=(0., 20.), ylims=(-0.0002, 0.0002),
    bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)
fig

# ## References

# [^ASG00]: Antoulas, Athanasios C., Danny C. Sorensen, and Serkan Gugercin. A survey of model reduction methods for large-scale systems. 2000.

# [^TLT16]: Tran, Hoang-Dung, Luan Viet Nguyen, and Taylor T. Johnson. *Large-scale linear systems from order-reduction (benchmark proposal).* 3rd Applied Verification for Continuous and Hybrid Systems Workshop (ARCH), Vienna, Austria. 2016.


text:
using Iris
using Test

@testset "Boundaries" begin
    @testset "1D" begin
        shape = Interval(0,1)
        @test_nowarn Boundary(shape,DirichletBC)
        @test_nowarn Boundary(shape,PML)
        @test_nowarn Boundary(shape,DirichletBC,PML{1}(.1),noBL{2}())
        bcs = (DirichletBC{1}(),DirichletBC{2}())
        bls = (PML{1}(.1),PML{1}(.1),PML{1}(.1))
        @test_throws ArgumentError Boundary(shape,bcs,bls)
        bls = (noBL{1}(.1),PML{2}(.2))
        bcs = (DirichletBC{1}(),DirichletBC{2}(),DirichletBC{3}())
        @test_throws ArgumentError Boundary(shape,bcs,bls)
        bcs = (DirichletBC{1}(),DirichletBC{2}())
        @test_nowarn bnd = Boundary(shape,bcs,bls)
        bnd = Boundary(shape,bcs,bls)
        @test ndims(bnd)==1
        @test Iris.Common.Boundaries.nsides(bnd)==2
    end
    @testset "2D" begin
        @testset "Square" begin
            shape = Square(1,0,0)
            @test_nowarn Boundary(shape,DirichletBC)
            @test_nowarn Boundary(shape,PML)
            @test_nowarn Boundary(shape,DirichletBC,PML{1}(.1),noBL{2}())
            bcs = (DirichletBC{1}(),DirichletBC{2}())
            bls = (PML{1}(),PML{1}(.1),PML{1}(.1),cPML{2}(.3),noBL{4}())
            @test_throws ArgumentError Boundary(shape,bcs,bls)
            bls = (noBL{1}(.1),PML{2}(.2))
            bcs = (DirichletBC{1}(),DirichletBC{2}(),DirichletBC{3}(),MatchedBC{3}(),NeumannBC{1}())
            @test_throws ArgumentError Boundary(shape,bcs,bls)
            bls = (noBL{3}(.1),PML{4}(.2),cPML{2}(),noBL{1}())
            bcs = (DirichletBC{2}(),NeumannBC{1}(),DirichletBC{4}(),DirichletBC{3}())
            @test_nowarn bnd = Boundary(shape,bcs,bls)
            bnd = Boundary(shape,bcs,bls)
            @test ndims(bnd)==2
            @test Iris.Common.Boundaries.nsides(bnd)==4
        end
        @testset "Circle" begin
            shape = Circle(1,0,0)
            @test_nowarn Boundary(shape)
            @test_nowarn Boundary(shape,DirichletBC)
            @test_nowarn Boundary(shape,PML)
            @test_throws ArgumentError Boundary(shape,DirichletBC,PML{1}(.1),noBL{2}())
            bcs = (DirichletBC{1}(),)
            bls = (PML{1}(),)
            bnd = Boundary(shape,bcs,bls)
            @test_nowarn Boundary(shape,bcs,bls)
            bls = (noBL{1}(.1),)
            bcs = (DirichletBC{1}(),)
            @test_nowarn Boundary(shape,bcs,bls)
            bls = (noBL{3}(.1),PML{4}(.2),cPML{2}(),noBL{1}())
            bcs = (DirichletBC{2}(),NeumannBC{1}(),DirichletBC{4}(),DirichletBC{3}())
            @test_throws ArgumentError bnd = Boundary(shape,bcs,bls)
            @test ndims(bnd)==2
            @test Iris.Common.Boundaries.nsides(bnd)==1
        end
    end
end


text:
module DynaUtils

include("DynaPlanning.jl")

using LightGraphs, LightGraphsFlows, SimpleWeightedGraphs
using GraphPlot, Plots
using Colors
using Compose, Cairo
using SparseArrays
using DataStructures
using Random
using Logging

import .DynaPlanning

"""
price_of_anarchy(r, c)

Compute price of anarchy ρ(G,r,c) for instances with linear cost functions. r and c are vectors of traffic rates and cost functions, respectively.

"""
function price_of_anarchy(r, c)

    C_f = sum([c[i](r[i]) * r[i] for i = 1:length(r)])

    # Optimal flow for the instance (G,r/2,c)
    r_star = 0.5 * r
    C_f_star = sum([c[i](r_star[i]) * r_star[i] for i = 1:length(r_star)])

    # Augmentation to (G,r,c)
    δ = 1
    C_f_star_augm = δ * (C_f - C_f_star)

    return C_f / C_f_star_augm

end

"""
    find_flow(G)

Compute the edge flow for graph `G`.

"""
function find_flow(SInstance)

    A = adjacency_matrix(SInstance.network)

    I, J, vals = findnz(A)

    # init f_e
    f_e = sparse(I, J, zeros(length(vals)))

    for a in SInstance.agents
        thisAgent = a
        thisPath = thisAgent.path
        edge_k = first(thisPath)
        f_e[edge_k.src, edge_k.dst] += 1.0
    end

    # normalize
    f_e = sparse(I, J, findnz(A)[3] / sum(findnz(A)[3]))

    return f_e
end

"""
    find_supply_demand(G)

Compute the supply and the demand for `SInstance`.

"""
function find_supply_demand(SInstance)

    # init d
    d = zeros(nv(SInstance.network), 1)

    for a in SInstance.agents
        thisAgent = a
        thisPath = thisAgent.path

        @debug thisPath

        edge_src = first(thisPath)
        d[edge_src.src] += 1.0

        edge_sink = last(thisPath)
        d[edge_sink.dst] -= 1.0
    end

    # normalize
    total_demand = sum([x] for x in d if x > zero(x))

    d *= 1 / total_demand

    @debug "Supply/demand:" d

    return d
end

"""
    recursive_DFS!(genericgraph, visited, sink)

Compute -recursively- the list of paths between `visited` (source) and `sink` (destination) pairs for `genericgraph`.

"""
function recursive_DFS!(pathList, genericgraph, visited, sink)

    nodes = LightGraphs.neighbors(genericgraph, visited[length(visited)])

    for i in nodes
        if i == sink
            push!(visited, i)
            push!(pathList, visited)
            pop!(visited)
        elseif !in(i, visited)
            push!(visited, i)
            recursive_DFS!(pathList, genericgraph, visited, sink)
            pop!(visited)
        end
    end

end

"""
    find_all_paths(G, src, dst)

Compute -recursively- all paths between `src` (source) and `dst` (destination) pairs for graph `G`.

"""
function find_all_paths(G, src, dst)

    A = adjacency_matrix(G)

    all_paths = Array{Array{Any},1}()

    DynaPlanning.recursive_DFS!(all_paths, SimpleDiGraph(A), [src], dst)

    return all_paths

end

function plot_graph(G, locs_x, locs_y)

    A = adjacency_matrix(G)

    nodelabel = collect(1:nv(SimpleDiGraph(A)))

    draw(
        PNG("graph.png", 10cm, 10cm),
        gplot(SimpleDiGraph(A), locs_x, locs_y, nodelabel = nodelabel),
    )

end

end

max_stars_repo_name: mgoar/dynaplanning

max_issues_repo_name: mgoar/dynaplanning

max_forks_repo_name: mgoar/dynaplanning


text:
export A_mul_B!

function A_mul_B!( alpha::Float64,
                   A::SparseMatrixCSC{Float64,Int},
                   x::Array{Float64},
                   beta::Float64,
                   y::Array{Float64},
                   nthreads::Int64=0 )
# Real:  y = beta*y  +  alpha * A*x

   if nthreads == 0
      Base.A_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

	n,m = size(A)
    nvec = size(x,2)

   if size(x,1) != m || size(y,1) != n
      throw(DimensionMismatch("length(x) != m || length(y) != n"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:a_mul_b_rr_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}),
                &nthreads, &nvec, &m, &n,    &alpha,   &beta,              A.nzval,      A.rowval,   A.colptr,   x,   y);

end  # function A_mul_B!

#------------------------------------------------------------------------------

function A_mul_B!( alpha::Complex128,
                   A::SparseMatrixCSC{Float64,Int},
                   x::Array{Complex128},
                   beta::Complex128,
                   y::Array{Complex128},
                   nthreads::Int64=0 )
# Real, Complex A:  y = beta*y  +  alpha * A*x

   if nthreads == 0
      Base.A_mul_B!( alpha, complex(A), x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

	n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != m || size(y,1) != n
      throw(DimensionMismatch("length(x) != m || length(y) != n"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:a_mul_b_rc_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                   &nthreads, &nvec, &m, &n,     &alpha,   &beta,              A.nzval,      A.rowval,   A.colptr,   convert(Ptr{Complex128}, pointer(x)),  convert(Ptr{Complex128}, pointer(y)));

end  # function A_mul_B!

#------------------------------------------------------------------------------

function A_mul_B!( alpha::Complex128,
                   A::SparseMatrixCSC{Complex128,Int},
                   x::Array{Complex128},
                   beta::Complex128,
                   y::Array{Complex128},
                   nthreads::Int64=0 )
# Complex A:  y = beta*y  +  alpha * A*x

   if nthreads == 0
      Base.A_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

	n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != m || size(y,1) != n
      throw(DimensionMismatch("length(x) != m || length(y) != n"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:a_mul_b_cc_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}, Ptr{Complex128}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                   &nthreads, &nvec, &m, &n,     &alpha,   &beta,              convert(Ptr{Complex128}, pointer(A.nzval)),      A.rowval,   A.colptr,   convert(Ptr{Complex128}, pointer(x)),  convert(Ptr{Complex128}, pointer(y)));

end  # function A_mul_B!


text:
__precompile__(true)

module RTFOOL

using Combinatorics, DataStructures

export StateSpace, Context, degeneracies, timestep
export Monotones

include("space.jl")
include("boltzmann.jl")
include("context.jl")
include("monotones.jl")

end