text:
# coding=utf-8
using Distributions: sample
using .ObjMesh
abstract type AbstractWorldObj end

mutable struct WorldObj <: AbstractWorldObj
    visible::Bool
    color::Vec3
    domain_rand::Bool
    angle::Float32
    y_rot::Float32
    pos::Vector{Float32}
    scale::Float32
    min_coords::Vector{Float32}
    max_coords::Vector{Float32}
    kind::String
    mesh::ObjectMesh
    optional::Bool
    static::Bool
    safety_radius::Float32
    obj_corners::Array{Float32, 2}
    obj_norm::Array{Float32, 2}
end

function WorldObj(obj, domain_rand, safety_radius_mult)
    ##
    #Initializes the object and its properties
    ##
    # XXX this is relied on by things but it is not always set
    # (Static analysis complains)
    visible = true
    # same
    color = Vec3([0f0])
    # maybe have an abstract method is_visible, get_color()

    dict_vals = process_obj_dict(obj, safety_radius_mult)
    pos, scale, min_coords, max_coords, y_rot = dict_vals[6:end]
    dict_vals = dict_vals[1:5]

    angle = deg2rad(y_rot)

    # Find corners and normal vectors assoc world object
    obj_corners = generate_corners(pos, min_coords, max_coords, angle, scale)
    obj_norm = generate_norm(obj_corners)

    WorldObj(visible, color, domain_rand, angle, y_rot, pos,
             scale, min_coords, max_coords, dict_vals...,
             obj_corners, obj_norm)
end

function process_obj_dict(obj, safety_radius_mult)
    kind = obj["kind"]
    mesh = obj["mesh"]
    pos = obj["pos"]
    scale = obj["scale"]
    optional = obj["optional"]
    min_coords = obj["mesh"].min_coords
    max_coords = obj["mesh"].max_coords
    static = obj["static"]
    safety_radius = safety_radius_mult * calculate_safety_radius(mesh, scale)
    y_rot = obj["y_rot"]

    return (kind, mesh, optional, static, safety_radius, pos, scale,
            min_coords, max_coords, y_rot)
end

_obj_corners(wobj::WorldObj) = wobj.obj_corners
_obj_corners(obj::AbstractWorldObj) = _obj_corners(obj.wobj)

_obj_norm(wobj::WorldObj) = wobj.obj_norm
_obj_norm(obj::AbstractWorldObj) = _obj_norm(obj.wobj)

function _set_color!(wobj::WorldObj, color::Vec3)
    wobj.color = color
end

_set_color!(obj::AbstractWorldObj, color::Vec3) = _set_color!(obj.wobj, color)

function _set_visible!(wobj::WorldObj, val::Bool)
    wobj.visible = val
end

_set_visible!(obj::AbstractWorldObj, val::Bool) = _set_visible!(obj.wobj, val)
_visible(wobj::WorldObj) = wobj.visible
_visible(obj::AbstractWorldObj) = _visible(obj.wobj)

_optional(wobj::WorldObj) = wobj.optional
_optional(obj::AbstractWorldObj) = _optional(obj.wobj)

_pos(wobj::WorldObj) = wobj.pos
_pos(obj::AbstractWorldObj) = _pos(obj.wobj)

_max_coords(wobj::WorldObj) = wobj.max_coords
_max_coords(obj::AbstractWorldObj) = _max_coords(obj.wobj)

_scale(wobj::WorldObj) = wobj.scale
_scale(obj::AbstractWorldObj) = _scale(obj.wobj)

_y_rot(wobj::WorldObj) = wobj.y_rot
_y_rot(obj::AbstractWorldObj) = _y_rot(obj.wobj)

_mesh(wobj::WorldObj) = wobj.mesh
_mesh(obj::AbstractWorldObj) = _mesh(obj.wobj)

_color(wobj::WorldObj) = wobj.color
_color(obj::AbstractWorldObj) = _color(obj.wobj)

_safety_radius(wobj::WorldObj) = wobj.safety_radius
_safety_radius(obj::AbstractWorldObj) = _safety_radius(obj.wobj)

_static(wobj::WorldObj) = wobj.static
_static(obj::AbstractWorldObj) = _static(obj.wobj)

_kind(wobj::WorldObj) = wobj.kind
_kind(obj::AbstractWorldObj) = _kind(obj.wobj)

function render(obj::AbstractWorldObj, draw_bbox::Bool)
    ##
    #Renders the object to screen
    ##
    !_visible(obj) && return

    # Draw the bounding box
    if draw_bbox
    #    gl.glColor3f(1, 0, 0)
    #    gl.glBegin(gl.GL_LINE_LOOP)
    #    gl.glVertex3f(self.obj_corners.T[0, 0], 0.01, self.obj_corners.T[1, 0])
    #    gl.glVertex3f(self.obj_corners.T[0, 1], 0.01, self.obj_corners.T[1, 1])
    #    gl.glVertex3f(self.obj_corners.T[0, 2], 0.01, self.obj_corners.T[1, 2])
    #    gl.glVertex3f(self.obj_corners.T[0, 3], 0.01, self.obj_corners.T[1, 3])
    #    gl.glEnd()
    end

    #gl.glTranslatef(*self.pos)
    #gl.glScalef(self.scale, self.scale, self.scale)
    #gl.glRotatef(self.y_rot, 0, 1, 0)
    #gl.glColor3f(*self.color)

    transformation_mat = get_transformation_mat(_pos(obj), _scale(obj), _y_rot(obj), (0,1,0))
    obj_mesh = _mesh(obj)
    transformed_vlists = transform_mesh(obj_mesh, transformation_mat)
    #let mesh store a matrix. make vec3 out of it only when you render
    Δed_faces = triangulate_faces.(transformed_vlists, obj_mesh.texclists, obj_mesh.clists, obj_mesh.textures)
    Δed_faces = vcat(Δed_faces...)
end

function make_Δ(v1::Vec3, vert2::Vector{Float32}, vert3::Vector{Float32}, mat::Material)
    v2 = Vec3(vert2[1:1], vert2[2:2], vert2[3:3])
    v3 = Vec3(vert3[1:1], vert3[2:2], vert3[3:3])
    Triangle(v1, v2, v3, mat)
end

function make_Δ(v1::Vec3, vert2::Vector{Float32}, vert3::Vector{Float32},
                tc1::Vector{Float32}, tc2::Vector{Float32}, tc3::Vector{Float32},
                texture::NamedTuple)
    v2 = Vec3(vert2[1:1], vert2[2:2], vert2[3:3])
    v3 = Vec3(vert3[1:1], vert3[2:2], vert3[3:3])

    uv_coordinates = [tc1, tc2, tc3]

    mat = Material(;uv_coordinates=uv_coordinates, texture...)
    Triangle(v1, v2, v3, mat)
end

function make_Δ(v1::Vec3, vert2::Vector{Float32}, vert3::Vector{Float32},
                tc1::Vector{Float32}, tc2::Vector{Float32}, tc3::Vector{Float32},
                color::Vec3)
    v2 = Vec3(vert2[1:1], vert2[2:2], vert2[3:3])
    v3 = Vec3(vert3[1:1], vert3[2:2], vert3[3:3])

    uv_coordinates = [tc1, tc2, tc3]

    mat = Material(;uv_coordinates=uv_coordinates, color_diffuse=color)
    Triangle(v1, v2, v3, mat)
end

function triangulate_faces(list_verts::Matrix{Float32}, texc::Matrix{Float32},
                           color::Vec3, texture::NamedTuple)
    v1 = Vec3(list_verts[1, 1:1], list_verts[1, 2:2], list_verts[1, 3:3])
    num_verts = size(list_verts, 1)

    texture_ = (color_diffuse = color,
                texture_diffuse = texture.texture_diffuse)

    Δs = map(i->make_Δ(v1, list_verts[i, :], list_verts[i+1, :], texc[1, :],
                       texc[i, :], texc[i+1, :], texture_), 2:num_verts-1)
    return Δs
end

function triangulate_faces(list_verts::Matrix{Float32}, texc::Matrix{Float32},
                           color::Vec3, texture::Nothing=nothing)
    v1 = Vec3(list_verts[1, 1:1], list_verts[1, 2:2], list_verts[1, 3:3])
    num_verts = size(list_verts, 1)

    Δs = map(i->make_Δ(v1, list_verts[i, :], list_verts[i+1, :], texc[1, :],
                       texc[i, :], texc[i+1, :], color), 2:num_verts-1)
    return Δs
end

function triangulate_faces(list_verts::Matrix{Float32}, texc::Nothing,
                           color::Vec3, texture::Nothing=nothing)
    v1 = Vec3(list_verts[1, 1:1], list_verts[1, 2:2], list_verts[1, 3:3])
    num_verts = size(list_verts, 1)

    mat = Material(;color_diffuse=color)

    Δs = map(i->make_Δ(v1, list_verts[i, :], list_verts[i+1, :], mat), 2:num_verts-1)
    return Δs
end

function triangulate_faces(list_verts::Matrix{Float32}, texc::Nothing,
                           color::Vec3, texture::NamedTuple)
    v1 = Vec3(list_verts[1, 1:1], list_verts[1, 2:2], list_verts[1, 3:3])
    num_verts = size(list_verts, 1)

    texture_ = (color_diffuse = color,
                texture_diffuse = texture.texture_diffuse)
    mat = Material(;texture_...)

    Δs = map(i->make_Δ(v1, list_verts[i, :], list_verts[i+1, :], mat), 2:num_verts-1)
    return Δs
end

function triangulate_faces(list_verts::Array{Float32, 3}, list_texc::Array{Float32, 3},
                           list_colors::Vector{Vec3}, texture::Union{NamedTuple,Nothing})
    # list_verts: 3D array representing list of faces along 3rd dim.
    # Every face has 3 vertices. Vertices are along 1st dim.
    # X, Y and Z component of vertices, along 2nd dim.
    # Similarly for texture coordinates

    @assert size(list_verts, 1) == 3
    @assert size(list_verts, 3) ==  size(list_colors, 1)

    num_faces = size(list_verts, 3)
    vt = Vector{Triangle}()
    for i in 1:num_faces
        Δ = triangulate_faces(list_verts[:, :, i], list_texc[:, :, i],
                              list_colors[i], texture)
        vt = vcat(vt, Δ)
    end
    return vt
end

function triangulate_faces(list_verts::Array{Float32,3}, list_texc::Array{Float32,3},
                           list_colors::Vector{Vec3}, texture::Vec3)
   texture_ = (texture_diffuse = texture,)
   triangulate_faces(list_verts, list_texc, list_colors, texture_)
end

get_transformation_mat(pos::Vector, scale, θ, rot_axis) =
    rotate_mat(θ) * scale_mat(scale) * translation_mat(pos)

translation_mat(pos...) = translation_mat(collect(Float32.(pos)))

function translation_mat(pos::Vector{Float32})
    #@assert length(pos) == 3
    #mat = Matrix{Float32}(I, 3, 3)
    mat = [
        1f0 0f0 0f0 0f0;
        0f0 1f0 0f0 0f0;
        0f0 0f0 1f0 0f0;
        pos[1] pos[2] pos[3] 1f0]

    return mat
end

scale_mat(scale...) = scale_mat(collect(Float32.(scale)))

scale_mat(scale::Float32) = scale_mat(ones(Float32, 3) * scale)

function scale_mat(scale::Vector{Float32})
    # scale is a vector of length 3
    scale_ = vcat(scale, [1f0])
    #mat = Matrix{Float32}(I, 4, 4)
    mat = [
        1f0 0f0 0f0 0f0;
        0f0 1f0 0f0 0f0;
        0f0 0f0 1f0 0f0;
        0f0 0f0 0f0 1f0]

    return mat .* scale_
end

function rotate_mat(θ, axis=(0,1,0))
    # axis: one-hot vector, each element corresponds to x, y or z axis
    θ = deg2rad(θ)
    #mat = Matrix{Float32}(I, 4, 4)
    mat = nothing
    if axis[2] == 1
        mat = [
        cos(θ) 0f0 -sin(θ) 0f0;
        0f0    1f0 0f0     0f0;
        sin(θ) 0f0 cos(θ)  0f0;
        0f0    0f0 0f0     1f0]
    elseif axis[1] == 1
        mat = [
        1f0    0f0    0f0     0f0;
        0f0    cos(θ) sin(θ)  0f0;
        0f0   -sin(θ) cos(θ)  0f0;
        0f0    0f0    0f0     1f0]
    else
        mat = [
         cos(θ)   sin(θ)   0f0  0f0;
        -sin(θ)   cos(θ)   0f0  0f0;
         0f0      0f0      1f0  0f0;
         0f0      0f0      0f0  1f0]
    end

    return mat
end

function transform_mat(mat::Matrix{T}, transformation_mat::Matrix{T}) where T
    numVecs, dim = size(mat)
    mat_ = hcat(mat, ones(eltype(mat), numVecs))
    return (mat_ * transformation_mat)[:, 1:dim]
end

function transform_mesh(vlist::Array{Float32,3}, transformation_mat::Array{Float32,2})
    w, h, num_faces = size(vlist)
    trans_vlist = map(i -> transform_mat(vlist[:, :, i], transformation_mat), 1:num_faces)
    tvlist = deepcopy(trans_vlist[1])
    for i in 2:num_faces
        tvlist = hcat(tvlist, trans_vlist[i])
    end
    return reshape(tvlist, w, h, num_faces)
end

transform_mesh(obj_mesh::ObjectMesh, transformation_mat::Array{Float32,2}) =
map(vlist->transform_mesh(vlist, transformation_mat), obj_mesh.vlists)

function projection_mat(fovy::T, aspect::T, zNear::T, zFar::T) where T<:Float32
    f = cot(fovy/2f0)
    r = zFar / zNear

    proj_mat = [f/aspect   0f0    0f0             0f0;
	        0f0        f      0f0             0f0;
	        0f0        0f0   (r+1f0)/(1f0-r) -1f0;
	        0f0        0f0  2f0*zFar/(1f0-r)  0f0]
end

# Below are the functions that need to
# be reimplemented for any dynamic object
function check_collision(wobj::WorldObj, agent_corners, agent_norm)
    ##
    #See if the agent collided with this object
    #For static, return false (static collisions checked w
    #numpy in a batch operation)
    ##
    !wobj.static && throw(NotImplementedError)

    return false
end

function proximity(wobj::WorldObj, agent_pos::Vector{Float32}, agent_safety_rad::Float32)
    ##
    #See if the agent is too close to this object
    #For static, return 0
    ##
    !wobj.static && throw(NotImplementedError)

    return 0f0
end

proximity(abs_wobj::AbstractWorldObj, agent_pos::Vector{Float32},
          agent_safety_rad::Float32) = proximity(abs_wobj.wobj, agent_pos, agent_safety_rad)

function step!(wobj::WorldObj, delta_time::Float32)
    ##
    #Use a motion model to move the object in the world
    ##
    !wobj.static && throw(NotImplementedError)
end

step!(abs_wobj::AbstractWorldObj, delta_time::Float32) = step!(abs_wobj.wobj, delta_time)

mutable struct DuckiebotObj <: AbstractWorldObj
    wobj::WorldObj
    follow_dist::Float32
    velocity::Float32
    max_iterations::Int
    gain::Float32
    trim::Float32
    radius::Float32
    k::Float32
    limit::Float32
    wheel_dist::Float32
    robot_width::Float32
    robot_length::Float32
end

function DuckiebotObj(obj, domain_rand::Bool, safety_radius_mult::Float32,
                      wheel_dist::Float32, robot_width::Float32, robot_length::Float32,
                      gain::Float32=2f0, trim::Float32=0f0, radius::Float32=0.0318f0,
                      k::Float32=27f0, limit::Float32=1f0)
    wobj = WorldObj(obj, domain_rand, safety_radius_mult)
    follow_dist, velocity = nothing, nothing
    if domain_rand
        follow_dist = Float32(rand(Uniform(0.3f0, 0.4f0)))
        velocity = Float32(rand(Uniform(0.05f0, 0.15f0)))
    else
        follow_dist = 0.3f0
        velocity = 0.1f0
    end

    max_iterations = 1000

    DuckiebotObj(wobj, follow_dist, velocity, max_iterations, gain, trim,
                 radius, k, limit, wheel_dist, robot_width, robot_length)
end

function get_dir_vec(db_obj::DuckiebotObj, angle::Float32)
    x = cos(angle)
    z = -sin(angle)
    return [x, 0f0, z]
end

function get_right_vec(db_obj::DuckiebotObj, angle::Float32)
    x = sin(angle)
    z = cos(angle)
    return [x, 0f0, z]
end

function check_collision(db_obj::DuckiebotObj, agent_corners::Matrix{Float32}, agent_norm::Matrix{Float32})
    ##
    #See if the agent collided with this object
    ##
    return intersects_single_obj(
        agent_corners,
        permutedims(db_obj.wobj.obj_corners),
        agent_norm,
        db_obj.wobj.obj_norm
    )
end

function _update_pos(db_obj::DuckiebotObj, action::Vector{Float32}, deltaTime::Float32)
    vel, angle = action

    # assuming same motor constants k for both motors
    k_r = db_obj.k
    k_l = db_obj.k

    # adjusting k by gain and trim
    k_r_inv = (db_obj.gain + db_obj.trim) / k_r
    k_l_inv = (db_obj.gain - db_obj.trim) / k_l

    omega_r = (vel + 0.5f0 * angle * db_obj.wheel_dist) / db_obj.radius
    omega_l = (vel - 0.5f0 * angle * db_obj.wheel_dist) / db_obj.radius

    # conversion from motor rotation rate to duty cycle
    u_r = omega_r * k_r_inv
    u_l = omega_l * k_l_inv

    # limiting output to limit, which is 1.0 for the duckiebot
    u_r_limited = max(min(u_r, db_obj.limit), -db_obj.limit)
    u_l_limited = max(min(u_l, db_obj.limit), -db_obj.limit)

    # If the wheel velocities are the same, then there is no rotation
    if u_l_limited == u_r_limited
        db_obj.wobj.pos = db_obj.wobj.pos .+ deltaTime * u_l_limited * get_dir_vec(db_obj, db_obj.wobj.angle)
        return
    end

    # Compute the angular rotation velocity about the ICC (center of curvature)
    w = (u_r_limited - u_l_limited) / db_obj.wheel_dist

    # Compute the distance to the center of curvature
    r = (db_obj.wheel_dist * (u_l_limited + u_r_limited)) / (2 * (u_l_limited - u_r_limited))

    # Compute the rotation angle for this time step
    rotAngle = w * deltaTime

    # Rotate the robot's position around the center of rotation
    r_vec = get_right_vec(db_obj, db_obj.wobj.angle)
    px, py, pz = db_obj.wobj.pos
    cx = px + r * r_vec[1]
    cz = pz + r * r_vec[3]
    npx, npz = rotate_point(px, pz, cx, cz, rotAngle)

    # Update position
    db_obj.wobj.pos = [npx, py, npz]

    # Update the robot's direction angle
    db_obj.wobj.angle += rotAngle
    db_obj.wobj.y_rot += rad2deg(rotAngle)

    # Recompute the bounding boxes (BB) for the duckiebot
    db_obj.wobj.obj_corners = agent_boundbox(
        db_obj.wobj.pos,
        db_obj.robot_width,
        db_obj.robot_length,
        get_dir_vec(db_obj, db_obj.wobj.angle),
        get_right_vec(db_obj, db_obj.wobj.angle)
    )
end

mutable struct DuckieObj <: AbstractWorldObj
    wobj::WorldObj
    pedestrian_wait_time::Float32
    vel::Float32
    heading::Vector{Float32}
    start::Vector{Float32}
    center::Vector{Float32}
    pedestrian_active::Bool
    wiggle::Float32
    time::Float32
end

function DuckieObj(obj, domain_rand::Bool, safety_radius_mult::Float32, walk_distance::Float32)
    wobj = WorldObj(obj, domain_rand, safety_radius_mult)

    walk_distance = walk_distance + 0.25f0

    # Dynamic duckie stuff

    # Randomize velocity and wait time
    if wobj.domain_rand
        pedestrian_wait_time = Float32(rand(3:19))
        vel = abs(Float32(rand(Normal(0.02f0, 0.005f0))))
    else
        pedestrian_wait_time = 8f0
        vel = 0.02f0
    end

    # Movement parameters
    heading = heading_vec(wobj.angle)
    start = copy(wobj.pos)
    center = wobj.pos
    pedestrian_active = false

    # Walk wiggle parameter
    wiggle = sample([14, 15, 16], 1)[1]
    wiggle = Float32(π / wiggle)

    time = 0f0

    DuckieObj(wobj, pedestrian_wait_time, vel, heading, start,
              center, pedestrian_active, wiggle, time)
end

function check_collision(dobj::DuckieObj, agent_corners::Matrix{Float32}, agent_norm::Matrix{Float32})
    ##
    #See if the agent collided with this object
    ##
    return intersects_single_obj(
        agent_corners,
        permutedims(dobj.wobj.obj_corners),
        agent_norm,
        dobj.wobj.obj_norm
    )
end

function proximity(db_obj::Obj, agent_pos::Vector{Float32}, agent_safety_rad::Float32) where Obj <: Union{DuckieObj, DuckiebotObj}
    ##
    #See if the agent is too close to this object
    #based on a heuristic for the "overlap" between
    #their safety circles
    ##
    d = norm(agent_pos .- db_obj.wobj.pos)
    score = d - agent_safety_rad - db_obj.wobj.safety_radius

    return min(0f0, score)
end

function step!(dobj::DuckieObj, delta_time::Float32)
    ##
    #Use a motion model to move the object in the world
    ##

    dobj.time += delta_time

    # If not walking, no need to do anything
    if !dobj.pedestrian_active
        dobj.pedestrian_wait_time -= delta_time
        if dobj.pedestrian_wait_time ≤ 0
            dobj.pedestrian_active = true
        end
        return
    end

    # Update centers and bounding box
    vel_adjust = dobj.heading * dobj.vel
    dobj.center += vel_adjust
    dobj.obj_corners += vel_adjust[[1, end]]

    distance = norm(dobj.center .- dobj.start)

    distance > dobj.walk_distance && finish_walk(dobj)

    dobj.pos = dobj.center
    angle_delta = dobj.wiggle * sin(48dobj.time)
    dobj.wobj.y_rot = rad2deg(dobj.wobj.angle + angle_delta)
    dobj.wobj.obj_norm = generate_norm(dobj.obj_corners)
end

function finish_walk(dobj::DuckieObj)
    ##
    #After duckie crosses, update relevant attributes
    #(vel, rot, wait time until next walk)
    ##
    dobj.start = copy(dobj.center)
    dobj.angle += Float32(π)
    dobj.pedestrian_active = false

    if dobj.wobj.domain_rand
        # Assign a random velocity (in opp. direction) and a wait time
        # TODO: Fix this: This will go to 0 over time
        dobj.vel = -sign(dobj.vel) * abs(Float32(rand(Normal(0.02f0, 0.005f0))))
        dobj.pedestrian_wait_time = rand(3:19)
    else
        # Just give it the negative of its current velocity
        dobj.vel *= -1f0
        dobj.pedestrian_wait_time = 8
    end
end

mutable struct TrafficLightObj <: AbstractWorldObj
    wobj::WorldObj
    texs::Vector{Vec3}
    time::Float32
    freq::Float32
    pattern::Int
end

function TrafficLightObj(obj, domain_rand, safety_radius_mult)
    wobj = WorldObj(obj, domain_rand, safety_radius_mult)

    texs = [
        load_texture(get_file_path("src/textures", "trafficlight_card0", "jpg")),
        load_texture(get_file_path("src/textures", "trafficlight_card1", "jpg"))
    ]
    time = 0f0

    # Frequency and current pattern of the lights
    if wobj.domain_rand
        freq = Float32(rand(Normal(4, 7)))
        pattern = rand(1:2)
    else
        freq = 5f0
        pattern = 1
    end

    # Use the selected pattern
    wobj.mesh.textures[1] = texs[pattern]

    TrafficLightObj(wobj, texs, time, freq, pattern)
end

function step!(tl_obj::TrafficLightObj, delta_time::Float32)
    ##
    #Changes the light color periodically
    ##

    tl_obj.time += delta_time
    if round(tl_obj.time, digits=3) % tl_obj.freq == 0  # Swap patterns
        tl_obj.pattern ⊻= 3
        tl_obj.wobj.mesh.textures[1] = tl_obj.texs[tl_obj.pattern]
    end
end

function is_green(tl_obj::TrafficLightObj, direction='N')
    if direction == 'N' || direction == 'S'
        if tl_obj.wobj.y_rot == 45 || tl_obj.wobj.y_rot == 135
            return tl_obj.pattern == 1
        elseif tl_obj.wobj.y_rot == 225 || tl_obj.wobj.y_rot == 315
            return tl_obj.pattern == 2
        end
    elseif direction == 'E' || direction == 'W'
        if tl_obj.wobj.y_rot == 45 || tl_obj.wobj.y_rot == 135
            return tl_obj.pattern == 2
        elseif tl_obj.wobj.y_rot == 225 || tl_obj.wobj.y_rot == 315
            return tl_obj.pattern == 1
        end
    end
    return false
end


text: for i = 99:-1:1 print("\n$i bottles of beer on the wall\n$i bottles of beer\nTake one down, pass it around\n$(i-1) bottles of beer on the wall\n") end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Base.Test
using TransitionParser

@test 1 == 1


text:
module ExampleBase

import DifferentialEquations: solve, solve!
using DiffEqBase: init
using Parameters: @with_kw

using ...LyapunovExponents: LEProblem, ContinuousLEProblem, DiscreteLEProblem,
    LESolver, lyapunov_exponents
import ...LyapunovExponents: dimension, report

"""
A type to hold an example dynamical system and its known Lyapunov exponents.
"""
@with_kw struct LEExample{ProblemType}
    name
    phase_dynamics
    u0
    t_renorm
    param
    tangent_dynamics
    t_attr
    known_exponents
    atol
    rtol
    terminator_options
    integrator_options = []
end

LEExample{ProblemType}(
            name,
            phase_dynamics,
            u0,
            t_renorm,
            param,
            tangent_dynamics,
            t_attr,
            known_exponents,
            atol,
            rtol,
            terminator_options;
            kwargs...
        ) where {ProblemType} =
    LEExample{ProblemType}(;
            name = name,
            phase_dynamics = phase_dynamics,
            u0 = u0,
            t_renorm = t_renorm,
            param = param,
            tangent_dynamics = tangent_dynamics,
            t_attr = t_attr,
            known_exponents = known_exponents,
            atol = atol,
            rtol = rtol,
            terminator_options = terminator_options,
            kwargs...)

const ContinuousExample = LEExample{ContinuousLEProblem}
const DiscreteExample = LEExample{DiscreteLEProblem}

function LEProblem(example::LEExample{P}; kwargs...) where {P <: LEProblem}
    t_tran = example.t_attr * 0.1
    if P <: DiscreteLEProblem
        t_tran = ceil(typeof(example.t_attr), digits=t_tran)
    end
    P(example.phase_dynamics,
      example.u0,
      example.param;
      t_renorm = example.t_renorm,
      t_attr = example.t_attr,
      t_tran = t_tran,
      tangent_dynamics = example.tangent_dynamics,
      kwargs...)
end

dimension(example::LEExample) = length(example.u0)

function solve(example::LEExample;
               dim_lyap=dimension(example), kwargs...)
    solve(LEProblem(example; dim_lyap=dim_lyap);
          record = true,
          terminator = true,
          terminator_options = example.terminator_options,
          integrator_options = example.integrator_options,
          kwargs...)
end

mutable struct LEDemo
    example::LEExample
    prob::LEProblem
    solver::LESolver

    LEDemo(example, prob) = new(example, prob)
end

"""
    LEDemo(example::LEExample; <keyword arguments>)

Here is an example code for constructing an example dynamical system,
calculate its LEs and plot them:
```julia
using LyapunovExponents
using Plots
demo = solve!(LyapunovExponents.lorenz_63())
plot(demo)
```

Create a `LEDemo` holding an `example` and an appropriate `LEProblem`
created from the `example`.
"""
function LEDemo(example::LEExample; kwargs...)
    LEDemo(example, LEProblem(example; kwargs...))
end

"""
    solve!(demo::LEDemo; progress=-1, <keyword arguments>)

Initialize `demo.solver` from `demo.prob` and run
`solve!(demo.solver)` to calculate the Lyapunov exponents.
"""
function solve!(demo::LEDemo; progress = -1, kwargs...)
    demo.solver = init(demo.prob;
                       record = true,
                       terminator = true,
                       terminator_options = demo.example.terminator_options,
                       integrator_options = demo.example.integrator_options,
                       kwargs...)
    solve!(demo.solver,
           progress = progress)
    return demo
end

function Base.show(io::IO, demo::LEDemo)
    print(io, "Demo: ", demo.example.name)
    if isdefined(demo, :solver)
        print(io, ", ", demo.solver)
    else
        print(io, " [solver not initialized]")
    end
end

function report(io::IO, demo::LEDemo;
                convergence::Bool = true)
    printstyled(io, "Lyapunov Exponents Demo", color=:blue)
    println(io, ": ", demo.example.name)

    if isdefined(demo, :solver)
        report(io, demo.solver; convergence = false)
    else
        printstyled(io, "[solver not initialized]", color=:red)
        println(io)
    end

    if ! isempty(demo.example.known_exponents)
        known_exponents = demo.example.known_exponents

        printstyled(io, "Known LEs", color=:yellow)
        print(io, ": ")
        show(IOContext(io, :limit => true), known_exponents)
        println(io)

        LEs = lyapunov_exponents(demo.solver)
        dim = min(length(known_exponents), length(LEs))
        actual = LEs[1:dim]
        desired = known_exponents[1:dim]
        abserr = abs.(actual .- desired)
        relerr = abserr ./ max.(abs.(actual), abs.(desired))

        printstyled(io, "Abs. Error", color=:yellow)
        print(io, ": ")
        show(IOContext(io, :limit => true), abserr)
        println(io)

        printstyled(io, "Rel. Error", color=:yellow)
        print(io, ": ")
        show(IOContext(io, :limit => true), relerr)
        println(io)
   end

    if convergence
        report(io, demo.solver.sol.convergence)
    end
end

end


text:
# Pseudofrequencies
# =================

"Parametric abstract type to define pseudofrequencies types"
abstract type Pseudofrequencies end

"""
You can use `NoPseudofrequencies()` to avoid pseudocount corrections where a
`Pseudofrequencies` type is needed.
"""
struct NoPseudofrequencies <: Pseudofrequencies end

# BLOSUM based pseudofrequencies
# ==============================

"""
`BLOSUM_Pseudofrequencies` type. It takes to arguments/fields:

- `α` : Usually the number of sequences or sequence clusters in the MSA.
- `β` : The weight of the pseudofrequencies, a value close to 8.512 when `α` is the number of sequence clusters.
"""
struct BLOSUM_Pseudofrequencies <: Pseudofrequencies
    α::Float64
    β::Float64
end

"""
`_calculate_blosum_pseudofrequencies!{T}(Pab::ContingencyTable{T,2,UngappedAlphabet})`

This function uses the conditional probability matrix `BLOSUM62_Pij` to fill the temporal
array field of `Pab` with pseudo frequencies (`Gab`). This function needs the real
frequencies/probabilities `Pab` because they are used to estimate the pseudofrequencies.

`Gab = Σcd  Pcd ⋅ BLOSUM62( a | c ) ⋅ BLOSUM62( b | d )`
"""
function _calculate_blosum_pseudofrequencies!(Pab::ContingencyTable{T,2,UngappedAlphabet}) where T
    @assert gettotal(Pab) ≈ one(T) "The input should be a probability table (normalized)"
    pab  = getarray(gettable(Pab))
    gab  = Pab.temporal
    bl62 = getarray(gettable(BLOSUM62_Pij))
    total = zero(T)
    @inbounds for b in 1:20, a in 1:20
        gab[a,b] = zero(T)
        for d in 1:20
            bl62_db = bl62[d,b]
            for c in 1:20
                P = pab[c,d]
                if P != 0
                    # BLOSUM62_Pij[c,a] is p(a|c)
                    gab[a,b] += ( P * bl62[c,a] * bl62_db )
                end
            end
        end
        total += gab[a,b]
    end
    if total ≉ one(T)
        @inbounds for col in 1:20
            @simd for row in 1:20
                gab[row,col] /= total
            end
        end
    end
    Pab
end

"""
`apply_pseudofrequencies!{T}(Pab::ContingencyTable{T,2,UngappedAlphabet}, pseudofrequencies::BLOSUM_Pseudofrequencies)`

When a `BLOSUM_Pseudofrequencies(α,β)` is used, this function applies pseudofrequencies
`Gab` over `Pab`, as a weighted mean of both. It uses the conditional probability
matrix `BLOSUM62_Pij` and the real frequencies/probabilities `Pab` to estimate the
pseudofrequencies `Gab`. α is the weight of the real frequencies `Pab` and β the weight
of the pseudofrequencies.

`Gab = Σcd  Pcd ⋅ BLOSUM62( a | c ) ⋅ BLOSUM62( b | d )`
`Pab = (α ⋅ Pab + β ⋅ Gab )/(α + β)`
"""
function apply_pseudofrequencies!(Pab::ContingencyTable{T,2,UngappedAlphabet},
                                  pseudofrequencies::BLOSUM_Pseudofrequencies) where T
    α = T(pseudofrequencies.α)
    β = T(pseudofrequencies.β)
    if β == 0.0
        return(Pab)
    end
    _calculate_blosum_pseudofrequencies!(Pab)
    pab  = getarray(gettable(Pab))
    gab  = Pab.temporal
    frac = one(T) / ( α + β )
    @inbounds for col in 1:20
        @simd for row in 1:20
            pab[row,col] = ( α * pab[row,col] + β * gab[row,col] ) * frac
        end
    end
    update_marginals!(Pab)
    normalize!(Pab)
end

@inline apply_pseudofrequencies!(Pab::ContingencyTable, pseudofrequencies::NoPseudofrequencies) = Pab


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using ConsistencyResampling

DocMeta.setdocmeta!(
    ConsistencyResampling, :DocTestSetup, :(using ConsistencyResampling); recursive=true
)

makedocs(;
    modules=[ConsistencyResampling],
    authors="David Widmann",
    repo="https://github.com/devmotion/ConsistencyResampling.jl/blob/{commit}{path}#{line}",
    sitename="ConsistencyResampling.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/ConsistencyResampling.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md", "supported.md"],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/devmotion/ConsistencyResampling.jl",
    push_preview=true,
    devbranch="main",
)


text:
# this replaces a Colon by div(size(Array), prod(other dims))
function _to_newdims(x::Tuple, dims::Tuple{Vararg{Union{Colon,Integer}}})
    return __to_newdims(axes2length(x), dims)
end
_to_newdims(x::Tuple, dims::Tuple{Vararg{Integer}}) = dims

function __to_newdims(x::Int, dims::Tuple{Colon,Vararg{Any}})
    return (div(x, prod(tail(dims))), tail(dims)...)
end
function __to_newdims(x::Int, dims::Tuple{Integer,Vararg{Any}})
    f = first(dims)
    return (f, __to_newdims(x, f, tail(dims))...)
end
function __to_newdims(x::Int, val::Int, dims::Tuple{Integer,Vararg{Any}})
    f = first(dims)
    return (f, __to_newdims(x, val * f, tail(dims))...)
end
function __to_newdims(x::Int, val::Int, dims::Tuple{Colon,Vararg{Any}})
    return (newval, div(x, prod(tail(dims)) * val), tail(dims)...)
end
__to_newdims(x::Int, val::Int, dims::Tuple{Colon}) = (div(x, val),)

"""
    reshape_axes(a, dims)
    reshape_axes(a, dims...)

## Examples
```jldoctest
julia> axs = (Index{:a}(1:10), Index{:b}(1:10), Index(1:10))
```
"""
reshape_axes(x::AbstractArray, dims...) = reshape_axes(x,Tuple(dims))
reshape_axes(x::Tuple, dims...) = reshape_axes(x, Tuple(dims))

reshape_axes(x::AbstractArray, dims::Tuple) = reshape_axes(axes(x), dims)
reshape_axes(x::Tuple, dims::Tuple) = _reshape_axes(x, _to_newdims(x, dims))

function _reshape_axes(axs::Tuple{Any,Vararg}, dims::Tuple{Integer,Vararg})
    (reshape_axis(first(axs), first(dims)), _reshape_axes(tail(axs), tail(dims))...)
end
_reshape_axes(axs::Tuple{}, dims::Tuple{}) = ()

"""
    reshape_axis(x::AbstractUnitRange, len) -> AbstractUnitRange

Returns an axis of the same type as `x` and length `len`.
"""
reshape_axis(x, len::Integer) = set_length(x, len)

"""
    resize_axis!(x::AbstractUnitRange, len) -> x

Returns `x` resized to the same length as `len`.
"""
resize_axis!(x, len::Integer) = (set_length!(x, len); return x)


text:
#=
Implement a file syncing algorithm for two computers over a low-bandwidth network. What if we know the files in the two computers are mostly the same?
=#

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Test
using LogicCircuits
using DataFrames

@testset "MNIST Loader tests" begin

    train, _, test = mnist();

    @test train isa DataFrame
    @test test isa DataFrame

    @test num_features(train) == 784
    @test num_examples(train) == 60000
    @test num_features(test) == 784
    @test num_examples(test) == 10000
    @test isfpdata(train)
    @test !isbinarydata(train)
    @test isfpdata(test)
    @test !isbinarydata(test)

    @test feature_values(train, 1) isa Vector{Float32}
    @test feature_values(train, 784) isa Vector{Float32}
    @test convert(Matrix,train) isa Array{Float32,2};

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

    train = shuffle_examples(train)
    test = shuffle_examples(test)

    @test num_features(train) == 784
    @test num_examples(train) == 60000
    @test num_features(test) == 784
    @test num_examples(test) == 10000
    @test isfpdata(train)
    @test !isbinarydata(train)
    @test isfpdata(test)
    @test !isbinarydata(test)

    @test feature_values(train, 1) isa Vector{Float32}
    @test feature_values(train, 784) isa Vector{Float32}
    @test convert(Matrix,train) isa Array{Float32,2};

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

    train, _, test = threshold(train, nothing, test)

    @test num_features(train) == 784
    @test num_examples(train) == 60000
    @test num_features(test) == 784
    @test num_examples(test) == 10000

    @test !isfpdata(train)
    @test isbinarydata(train)
    @test !isfpdata(test)
    @test isbinarydata(test)

    @test feature_values(train, 1) isa BitVector
    @test feature_values(train, 784) isa BitVector
    @test feature_values(test, 1) isa BitVector
    @test feature_values(test, 784) isa BitVector

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

    train, _, test = mnist(true) # labeled MNIST

    @test train isa DataFrame
    @test test isa DataFrame

    @test num_features(train) == 784+1
    @test num_examples(train) == 60000
    @test num_features(test) == 784+1
    @test num_examples(test) == 10000
    @test !isfpdata(train) # because of the label being Int
    @test !isbinarydata(train)
    @test !isfpdata(test) # because of the label being Int
    @test !isbinarydata(test)

    @test feature_values(train, 1) isa Vector{Float32}
    @test feature_values(train, 784) isa Vector{Float32}

    ####################################
    train, valid, test = sampled_mnist()

    @test num_features(train) == 784
    @test num_examples(train) == 50000
    @test num_features(valid) == 784
    @test num_examples(valid) == 10000
    @test num_features(test) == 784
    @test num_examples(test) == 10000

    @test !isfpdata(train)
    @test isbinarydata(train)
    @test !isfpdata(valid)
    @test isbinarydata(valid)
    @test !isfpdata(test)
    @test isbinarydata(test)

    @test feature_values(train, 1) isa BitVector
    @test feature_values(train, 784) isa BitVector
    @test feature_values(valid, 1) isa BitVector
    @test feature_values(valid, 784) isa BitVector
    @test feature_values(test, 1) isa BitVector
    @test feature_values(test, 784) isa BitVector
end


text:
function getRealVs(V,evals)

    # get real vectors made from the eigenvectors

    # by paying attention to this order, things will always rotate CCW
    if abs(evals[1])>0  # if the eigenvalue with negative imaginary component comes first
        Vr = [V[:,1] + V[:,2], (V[:,1] - V[:,2])*im];
    else
        Vr = [V[:,2] + V[:,1], (V[:,2] - V[:,1])*im]
    end
    Vr = Vr / sqrt(2)
    return Vr
    #=
    # now get axes aligned so that plan is spread mostly along the horizontal axis()
    testProj = (Vr"*Ared[1:numAnalyzedTimes:end,:]')"; # just picks out the plan times
    rotV = princomp[testProj]
    crossProd = cross([rotV[:,1];0], [rotV[:,2];0])
    if crossProd[3] < 0, rotV[:,2] = -rotV[:,2]; end   # make sure the second vector is 90 degrees clockwise from the first
    Vr = Vr*rotV;

    # flip both axes if necessary so that the maximum move excursion is in the positive direction
    testProj = (Vr"*Ared')";  # all the times
    if max(abs(testProj[:,2])) > max(testProj[:,2])  # 2nd column is the putative "muscle potent" direction.
        Vr = -Vr
    end
    =#
end

max_stars_repo_name: alexmorley/jPCA

max_issues_repo_name: alexmorley/jPCA

max_forks_repo_name: alexmorley/jPCA


text:
# Phases in a 1D Earth

"""
    Phase(model, name, delta, depth, time, dtdd, inc, takeoff)
    Phase(..., pierce, distance, radius)

Construct a `Phase`, which represents a single event-station path and one single phase
arrival.

`Phase` objects have the following fields:

- `model`: Name of model used to calculate phase properties
- `name`: Name of phase
- `delta`: Epicentral distance in °
- `depth`: Depth of event below surface in km
- `time`: Travel time in s
- `dtdd`: Horizontal slowness (ray parameter) in s/°
- `inc`: Incidence angle at receiver, measured from downwards in °
- `takeoff`: Takeoff angle at source, measured from downwards in °

The following fields are filled if the ray path has been calculated via `path`:

- `distance`: `Vector` of distances in ° along the path
- `radius`: `Vector` of radii in km along the path
"""
struct Phase{T} <: AbstractPhase
    model::String
    name::String
    delta::T
    depth::T
    time::T
    dtdd::T
    inc::T
    takeoff::T
    pierce::Vector{T}
    distance::Vector{T}
    radius::Vector{T}
end
Phase(model, name, delta, depth, time, dtdd, inc, takeoff) =
    Phase{TauPyFloat}(model, name, delta, depth, time, dtdd, inc, takeoff, [], [], [])
function Phase(p::AbstractPhase, pierce::AbstractVector, distance::AbstractVector, radius::AbstractVector)
    length(pierce) == length(distance) == length(radius) ||
        throw(ArgumentError("pierce, distance and radius must be the same length"))
    Phase{TauPyFloat}(p.model, p.name, p.delta, p.depth, p.time, p.dtdd, p.inc, p.takeoff, pierce, distance, radius)
end

"""
    path(depth, distance, phase="ttall"; model="$DEFAULT_MODEL", cache=true) -> p::Vector{Phase}

Create a set of `Phase`s `p` which contain the computed ray path for a set of `phase`s
from an event `depth` km deep and at `distance`° away.
Optionally specify the `model` (one of: $(AVAILABLE_MODELS)).

If no arrivals are found for the geometry, an empty `Vector{Phase}` is returned.

Set `cache` to `false` to always recompute the ray path and ignore the module cache.
"""
function path(depth, distance, phase="ttall"; model=DEFAULT_MODEL, cache=true)
    phase = phase isa AbstractString ? [phase] : phase
    _call_taup(:get_ray_paths, model, depth, distance, phase; cache=cache)
end

"""
    travel_time(depth, distance, phase="ttall"; model="$DEFAULT_MODEL", cache=true) -> p::Vector{Phase}

Return a `Vector` of `Phase`s, given an event `depth` km deep and `distance`°
away.  Optionally specify a `phase` name; otherwise all arrivals are returned.
Optionally specify the model (one of: $(AVAILABLE_MODELS)).

If no arrivals are found for the geometry, an empty `Vector{Phase}` is returned.

Set `cache` to `false` to always recompute the travel time and ignore the module cache.
"""
function travel_time(depth, distance, phase="ttall"; model=DEFAULT_MODEL, cache=true)
    phase = phase isa AbstractString ? [phase] : phase
    _call_taup(:get_travel_times, model, depth, distance, phase; cache=cache)
end

"""Helper function which takes `obspy.taup.Arrivals` and returns `Phase`s.
Arguments should be the same as those to osbpy.taup calls, prefixed
by `model`."""
function _phases_from_arrivals(arr, model, depth, distance, phase)
    p = Vector{Phase{TauPyFloat}}()
    for a in arr
        name = a.name
        delta = a.distance
        depth = a.source_depth
        time = a.time
        dtdd = deg2rad(a.ray_param)
        inc = a.incident_angle
        takeoff = a.takeoff_angle
        path = a.path
        distance, radius = if path === nothing
            TauPyFloat[], TauPyFloat[]
        else
            pathlist = path.tolist()
            rad2deg.([pp[3] for pp in pathlist]), RADIUS[model] .- [pp[4] for pp in pathlist]
        end
        push!(p, Phase{TauPyFloat}(model, name, delta, depth, time, dtdd, inc, takeoff,
              [], distance, radius))
    end
    p
end

max_stars_repo_name: anowacki/TauPy


text:
using Test
using ProbabilisticCircuits
using Suppressor
using LightGraphs: nv, ne
using MetaGraphs: get_prop

@testset "Test loaded clt file" begin
   # TODO add more clt files
   t = zoo_clt("4.clt")
   @test t isa CLT
   @test_nowarn @suppress_out print_tree(t)
   @test nv(t) == 4
   @test ne(t) == 0
   for (v, p) in enumerate([.2, .1, .3, .4])
      cpt = get_prop(t, v, :cpt)
      @test cpt[1] + cpt[0] ≈ 1.0
      @test cpt[1] ≈ p
   end
end


text:
# ==============================================================================
# Xavier Gandibleux - November 2021
#   Implemented in Julia 1.6

# ==============================================================================
# Using ϵ-constraint method with vOptGeneric, compute the set of non-dominated
# points of the following problem:
#
#   Max  sum{j=1,...,n} C(k,j) x(j)                k=1,...,2
#   st   sum{j=1,...,n} A(i,j) x(j) <= B(i)        i=1,...,m
#                       x(j) = 0 or 1

# ==============================================================================
using vOptGeneric, JuMP, GLPK

# ==============================================================================
# Setting a Bi01IP instance and solving with vOptGeneric
function vSolveBi01IP(solverSelected, C, A, B)

  m, n = size(A)

  # ---- setting the model
  println("Building...")
  Bi01IP = vModel( solverSelected )
  @variable( Bi01IP, x[1:n], Bin )
  @addobjective( Bi01IP, Max, sum(C[1,j] * x[j] for j=1:n) )
  @addobjective( Bi01IP, Max, sum(C[2,j] * x[j] for j=1:n) )
  @constraint( Bi01IP, cte[i=1:m], sum(A[i,j] * x[j] for j=1:n) <= B[i])

  # ---- Invoking the solver (epsilon constraint method)
  println("Solving...")
  vSolve( Bi01IP, method=:epsilon, step=0.5, verbose=true )

  # ---- Querying the results
  println("Querying...")
  Y_N = getY_N( Bi01IP )

  return Y_N, x
end

# ==============================================================================

# Example on how to use it :

#solverSelected = GLPK.Optimizer
#Y_N, X_E = setBi01IP(solverSelected, C, A, B)

# ---- Displaying the results (X_E and Y_N)
#for n = 1:length(Y_N)
#    X = value.(X_E, n)
#    print(findall(elt -> elt ≈ 1, X))
#    println("| z = ",Y_N[n])
#end

#for n = 1:length(Y_N)
#    println(Y_N[n])
#end

max_stars_repo_name: xgandibleux/momhMKPstu

max_issues_repo_name: xgandibleux/momhMKPstu

max_forks_repo_name: xgandibleux/momhMKPstu


text:
## dialogs --

## modal window without methods

## Modal Message box
##
## static (no methods) dialog to display information
##
## Arguments:
##
## * `parent::Widget` used to get toolkit and to locate dialog
## * `text::String` main message
## * `informativeText::String` extra detail for dialog
## * `icon::MaybSymbol` needs implementing
## * `title::MaybeString` title for dialog window, not all OSes display this
##
## Returns
##
## no return value
function messagebox(parent::Widget, text::String; icon::Union(Nothing, String)=nothing,
                    title::Union(Nothing, String)=nothing,
                    informativeText::Union(Nothing, String)=nothing) # detail
    messagebox(parent.toolkit, parent, text, icon=icon, title=title, informativeText=informativeText)
    nothing
end

## Modal message box for confirmation
##
## static (no methods) dialog to get yes/no answer
##
## Arguments:
##
## * `parent::Widget` used to get toolkit and to locate dialog
## * `text::String` main message
## * `informativeText::String` extra detail for dialog
## * `icon::MaybSymbol` one of `nothing`, `:question`, `:info`, `:warning`, `:critical`
## * `title::MaybeString` title for dialog window, not all OSes display this
##
## Returns:
##
## * returns `:accept` or `:reject` symbol
##
function confirmbox(parent::Widget, text::String; icon::Union(Nothing, Symbol)=nothing,
                    title::Union(Nothing, String)=nothing,
                    informativeText::Union(Nothing, String)=nothing) # detail
    confirmbox(parent.toolkit, parent, text, icon=icon, title=title, informativeText=informativeText)
end

##
## Arguments:
##
## * `mode::Symbol` one of `:open`, `:multiple`, `:directory` or `:save`
## * `:message::String` message
## * `:title::String` title
## * `:filetypes::MaybeVector{Tuple}` a vector of tuples in form
##   `(label, extension`, as in `("jpeg", ".jpg") (In tcltk. For Qt,
##   this can be relaxed to `("Description", "*.jpg *.gif *.png")`,
##   say
##
## Returns:
##
## depends on `mode`:
## * `:open` a file name or `nothing`
## * `:multiple` a vector of file names or `nothing`
## * `directory` a directory name or `nothing`
## * `save` a filename or `nothing`
##
## TODO:
## setFile, setDirectory..., but hard to do with TclTk
function filedialog(parent::Widget;
                    mode::Symbol=:open, ## :open, :multiple, :directory, :save
                    message::String="",
                    title::String="",
                    filetypes::Union(Nothing, Vector{Tuple})=nothing)

    if findfirst([:open, :multiple, :directory, :save], mode) == 0
        error("mode is one of :open, :multiple, :directory, and :save")
    end

    filedialog(parent.toolkit, parent, mode=mode, message=message, title=title, filetypes=filetypes)
end

## a modeless dialog window

## Esc key in dialog closes window "reject"
## exec() for modal
type Dialog <: BinContainer
    o
    block
    model
    children
    attrs::Dict
    toolkit
    state::Union(Nothing, Symbol) ## (nothing, :accept, :reject)
    open::Function
    exec::Function
    accept::Function
    reject::Function
    result::Function
    done::Function
    close::Function
    function Dialog(widget, block, toolkit, model)
        self = new(widget, block,  model, {}, Dict(), toolkit, :reject)
        self.exec = () -> begin
            show_dialog(self.toolkit, self, true)
            setModal(self, true)
        end
        self.open = () -> begin
            show_dialog(self.toolkit, self, true)
        end
        self.done = (value::Symbol) -> begin
            self.state = value
            self.close()
        end
        self.accept = () -> self.done(:accept)
        self.reject = () -> self.done(:reject)
        self.result = () -> self.state
        self.close = () -> begin
            state = self.state
            show_dialog(self.toolkit, self, false)
            notify(self.model, "finished", self.state)
            self.state == :accept ? notify(self.model, "accepted") :
                                    self.state == :reject ? notify(self.model, "rejected") : nothing
        end
        self
    end
end

## return bin-like container

## Modaless dialog
##
## This constructor produces a non-modal dialog that is used like a bin-container
## Its design follows Qt:
## * one creates a dialog instance with specified buttons, say `dlg`
## * one can add a widget to the dialog through `push!`. This is usually a container for other widgets.
## * One displays the dialog by calling `dlg.exec()` or `dlg.open()`. If possible, `.exec()` will be modal.
## * To get a return value one connects to either the `finished (state)` signal or
##   the `accepted` or `rejected` signals. Otherwise, `dlg.state` holds the last state of the dialog.
##
## The dialog is not actually closed, just hidden so one can call `dlg.open()` or `dlg.exec()` again.
##
## Arguments
## * `parent::Widget` used to locate dialog
## * `buttons::Vector{Symbol}=[:ok]` vector of buttons from `:ok`, `:cancel`, `:close`, `:apply`, `:reset`, and `:help`.
## * `default::MaybeSymbol=[:ok]` if specified makes button the default
## * `title::String=""` title for dialog
## * `kwargs...` passed on as properties
##
## Signals
## * `finished (state)` called when button pressed
## * `accepted` called when ok button pressed
## * `rejected` called when cancel or close button pressed.
##
## TODO:
## implement default
function dialog(parent::Widget;
                buttons::Vector{Symbol} = [:ok],
                default::Union(Symbol, Nothing)=nothing,
                title::String="",
                kwargs...)

    model = EventModel()

    widget, block = dialog(parent.toolkit, parent, model, buttons=buttons; default=default, title=title) # use parent for placement too!

    dlg = Dialog(widget, block,  parent.toolkit, model)
    connect(model, "done", state -> dlg.done(state))

    show_dialog(dlg.toolkit, dlg, false)
    add_bindings(dlg.toolkit, dlg)
    for (k, v) in kwargs
        dlg[k] = v
    end

    dlg
end

setModal(dlg::Dialog, value::Bool) = setModal(dlg.toolkit, dlg, value)
setModaless(dlg::Dialog, value::Bool) = setModaless(dlg.toolkit, dlg, value)
destroy(dlg::Dialog) = destroy_dialog(dlg.toolkit, dlg)


text:
using MacroTools
import MacroTools: @q, combinearg, combinedef

export assess!,
    AbstractFitness,
    @fitness,
    FitnessFunction,
    setup!

abstract type AbstractFitness{G} <: Function end

setup!(f::AbstractFitness) = f

"""
    assess!(individual, f::AbstractFitness) -> FitnessValue

Calculate fitness of `individual`, returning a `FitnessValue`.  Stores the result of assessment in
the individual, so the actual calculation is done at most once.
"""
function assess!(individual::Individual{G}, f::AbstractFitness{>:G}) where {G}
    if individual.fitnessvalue === nothing
        individual.fitnessvalue = f(individual.genome)
    end

    individual.fitnessvalue
end

struct FitnessFunction{G, F} <: AbstractFitness{G}
    evaluate::F

    FitnessFunction{G}(f::F) where {G, F} = new{G, F}(f)
end

(f::FitnessFunction{G})(genome::G) where {G} = f.evaluate(genome)

macro fitness(fundef::Expr)
    def = splitdef(fundef)
    @assert (length(def[:args]) == 1) "Only single argument functions allowed!"
    argname, argtype, slurp, default = splitarg(def[:args][1])

    fname = def[:name]
    def[:name] = gensym(def[:name])
    fun = combinedef(def)
    @esc fname argtype fun

    return @q const $fname = FitnessFunction{$argtype}($fun)
end

## can't do this because of https://github.com/JuliaLang/julia/issues/14919
# """
#     (f::Fitness)(individual[, generation]) -> number

# Calculate fitness of `individual`, returning a number.
# """
# (f::F)(individual) where {F<:Fitness} = NaN
# (f::F)(individual, generation) where {F<:Fitness} = f(individual)


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7, 0.5)
equations = LinearScalarAdvectionEquation3D(advection_velocity)

initial_condition = initial_condition_convergence_test
boundary_conditions = BoundaryConditionDirichlet(initial_condition)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping.
function mapping(xi, eta, zeta)
  # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries
  # xi = 1.5 * xi_ + 1.5
  # eta = 1.5 * eta_ + 1.5
  # zeta = 1.5 * zeta_ + 1.5

  y = eta + 1/6 * (cos(1.5 * pi * (2 * xi - 3)/3) *
                   cos(0.5 * pi * (2 * eta - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  x = xi + 1/6 * (cos(0.5 * pi * (2 * xi - 3)/3) *
                  cos(2 * pi * (2 * y - 3)/3) *
                  cos(0.5 * pi * (2 * zeta - 3)/3))

  z = zeta + 1/6 * (cos(0.5 * pi * (2 * x - 3)/3) *
                    cos(pi * (2 * y - 3)/3) *
                    cos(0.5 * pi * (2 * zeta - 3)/3))

  return SVector(x, y, z)
end

cells_per_dimension = (4, 4, 4)
mesh = StructuredMesh(cells_per_dimension, mapping, periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.2)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_restart,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
"""
    PerformanceOpt

Abstract supertype for performance optimisation flags.
"""
abstract type PerformanceOpt end

"""
    NoOpt <: PerformanceOpt

    NoOpt()

Flag to run a simulation without performance optimisations besides basic high performance
programming. Still fast, but not intelligent about the work that it does: all cells are run
for all rules.

`NoOpt` is the default `opt` method.
"""
struct NoOpt <: PerformanceOpt end

"""
    Processor

Abstract supertype for selecting a hardware processor, such as ia CPU or GPU.
"""
abstract type Processor end

"""
    CPU <: Processor

Abstract supertype for CPU processors.
"""
abstract type CPU <: Processor end

"""
    SingleCPU <: CPU

    SingleCPU()

[`Processor`](@ref) flag that specifies to use a single thread on a single CPU.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=SingleCPU())
# or
output = sim!(output, rule; proc=SingleCPU())
```
"""
struct SingleCPU <: CPU end

"""
    ThreadedCPU <: CPU

    ThreadedCPU()

[`Processor`](@ref) flag that specifies to use a `Threads.nthreads()` CPUs.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=ThreadedCPU())
# or
output = sim!(output, rule; proc=ThreadedCPU())
```
"""
struct ThreadedCPU{L} <: CPU
    spinlock::L
end
ThreadedCPU() = ThreadedCPU(Base.Threads.SpinLock())
Base.Threads.lock(opt::ThreadedCPU) = lock(opt.spinlock)
Base.Threads.unlock(opt::ThreadedCPU) = unlock(opt.spinlock)

"""
    BoundaryCondition

Abstract supertype for flags that specify the boundary conditions used in the simulation,
used in [`inbounds`](@ref) and to update [`NeighborhoodRule`](@ref) grid padding.
These determine what happens when a neighborhood or jump extends outside of the grid.
"""
abstract type BoundaryCondition end

"""
    Wrap <: BoundaryCondition

    Wrap()

[`BoundaryCondition`](@ref) flag to wrap cordinates that boundary boundaries back to the
opposite side of the grid.

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Wrap())
# or
output = sim!(output, rule; boundary=Wrap())
```
"""
struct Wrap <: BoundaryCondition end

"""
    Remove <: BoundaryCondition

    Remove()

[`BoundaryCondition`](@ref) flag that specifies to assign `padval` to cells that overflow
grid boundaries. `padval` defaults to `zero(eltype(grid))` but can be assigned as a keyword
argument to an [`Output`](@ref).

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Remove())
# or
output = sim!(output, rule; boundary=Remove())
```
"""
struct Remove <: BoundaryCondition end


text:
@doc raw"""
# Mathematical Pendulum

"""
module Pendulum

    using GeometricEquations

    export pendulum_ode, pendulum_pode, pendulum_iode, pendulum_idae

    function pendulum_ode_f(t, x, f)
        f[1] = x[2]
        f[2] = sin(x[1])
        nothing
    end

    function pendulum_ode(x₀=[acos(0.4), 0.0])
        ODE(pendulum_ode_f, x₀)
    end

    function pendulum_pode_v(t, q, p, v)
        v[1] = p[1]
        nothing
    end

    function pendulum_pode_f(t, q, p, f)
        f[1] = sin(q[1])
        nothing
    end

    function pendulum_pode(q₀=[acos(0.4)], p₀=[0.0])
        PODE(pendulum_pode_v, pendulum_pode_f, q₀, p₀)
    end

    function pendulum_iode_α(t, q, v, p)
        p[1] = q[2]
        p[2] = 0
        nothing
    end

    function pendulum_iode_f(t, q, v, f)
        f[1] = sin(q[1])
        f[2] = v[1] - q[2]
        nothing
    end

    function pendulum_iode_g(t, q, λ, g)
        g[1] = 0
        g[2] = λ[1]
        nothing
    end

    function pendulum_iode_v(t, q, p, v)
        v[1] = q[2]
        v[2] = sin(q[1])
        nothing
    end

    function pendulum_iode(q₀=[acos(0.4), 0.0], p₀=[0.0, 0.0])
        IODE(pendulum_iode_α, pendulum_iode_f,
             pendulum_iode_g, q₀, p₀; v̄=pendulum_iode_v)
    end

    # function pendulum_pdae_v(t, q, p, v)
    #     v[1] = 0
    #     nothing
    # end
    #
    # function pendulum_pdae_f(t, q, p, f)
    #     f[1] = sin(q[1])
    #     nothing
    # end
    #
    # function pendulum_pdae_u(t, q, p, λ, u)
    #     u[1] = λ[1]
    #     nothing
    # end
    #
    # function pendulum_pdae_g(t, q, p, λ, g)
    #     g[1] = 0
    #     nothing
    # end
    #
    # # TODO
    # function pendulum_pdae_ϕ(t, q, p, ϕ)
    #     ϕ[1] = p[1] - q[2]
    #     nothing
    # end
    #
    # # TODO
    # function pendulum_pdae(q₀=[acos(0.4)], p₀=[0.0], λ₀=[0.0, 0.0])
    #     PDAE(pendulum_pdae_v, pendulum_pdae_f, pendulum_pdae_u, pendulum_pdae_g, pendulum_pdae_ϕ, q₀, p₀, λ₀)
    # end

    function pendulum_idae_u(t, q, p, λ, u)
        u[1] = λ[1]
        u[2] = λ[2]
        nothing
    end

    function pendulum_idae_g(t, q, p, λ, g)
        g[1] = 0
        g[2] = λ[1]
        nothing
    end

    function pendulum_idae_ϕ(t, q, p, ϕ)
        ϕ[1] = p[1] - q[2]
        ϕ[2] = p[2]
        nothing
    end

    function pendulum_idae(q₀=[acos(0.4), 0.0], p₀=[0.0, 0.0], λ₀=[0.0, 0.0])
        IDAE(pendulum_iode_f, pendulum_iode_α,
             pendulum_idae_u, pendulum_idae_g,
             pendulum_idae_ϕ, q₀, p₀, λ₀;
             v̄=pendulum_iode_v)
    end

end


text:
to_dot(g::GraphStore) = to_dot(g.graph, g.attributes)
function to_dot(g::GraphStore, ex::ExVertex)
  n = inclist(ExVertex, ExEdge{ExVertex}, is_directed=g.graph.is_directed)
  v = ExVertex(1, ex.label)
  v.attributes = ex.attributes
  add_vertex!(n, v)
  to_dot(n)
end
function to_dot(g::GraphStore, ex::ExEdge)
  n = inclist(ExVertex, ExEdge{ExVertex}, is_directed=g.graph.is_directed)
  vtxfrom = ExVertex(1, ex.source.label)
  vtxfrom.attributes = ex.source.attributes
  add_vertex!(n, vtxfrom)
  vtxto = ExVertex(2, ex.target.label)
  vtxto.attributes = ex.target.attributes
  add_vertex!(n, vtxto)
  add_edge!(n, ExEdge(1, vtxfrom, vtxto, ex.attributes))
  to_dot(n)
end

import GraphViz: Graph
draw(g::GraphStore) = Graph(to_dot(g))
draw(g::GraphStore, ex::Union{ExVertex,ExEdge}) = Graph(to_dot(g, ex))
draw(vertices::Vector{ExVertex}) = println("aa")
draw(edges::Vector{Vector{ExEdge{ExVertex}}}) = println("i")


text:
"""
    x = Unknown(T)

Constructs the Julia singleton `Unknown{T}()` which semantically represents that
`x` is equal to some instance of `T`, but which one is unknown. It represents
missing data rather than non-existance.

`Unknown` is used by the `Optional` type, which itself is a typealias for
a `Switch`.
"""
immutable Unknown{T}
end

(::Type{Unknown}){T}(::Type{T}) = Unknown{T}()

Base.show{T}(io::IO, ::Unknown{T}) = print(io, "Unknown(", T, ")")

"""
    Optional{T}(value)
    Optional{T}()

Construct an `Optional` container that contains either `value::T` or
`Unknown{T}()`.

This differs semantically from `Base.Nullable{T}` in that an `Unknown` state
means that some value `T` exists but is not known to the program, while a
null-valued `Nullable` represents the absence of a value. This subtle difference
allows us to implement 3-valued logic on `Optional{Bool}`, for example.
"""
typealias Optional{T} Switch{T,Unknown{T}}

(::Type{Optional})(x) = Optional{typeof(x)}(x)
(::Type{Optional{T}}){T}() = Optional{T}(Unknown{T}())

hasvalue{T}(x::Optional{T}) = x.which
get{T}(x::Optional{T}) = x.a.value

Base.show{T}(io::IO, x::Optional{T}) = hasvalue(x) ? print(io, "Optional(", get(x), ")") : print(io, "Optional{", T ,"}()")

# Three-valued logic
@inline Base.:!(::Unknown{Bool}) = Unknown{Bool}()

@inline Base.:&(::Unknown{Bool}, ::Unknown{Bool}) = Unknown{Bool}()
@inline Base.:&(x::Bool, ::Unknown{Bool}) = x ? Optional{Bool}() : Optional{Bool}(false)
@inline Base.:&(::Unknown{Bool}, y::Bool) = y ? Optional{Bool}() : Optional{Bool}(false)

@inline Base.:|(::Unknown{Bool}, ::Unknown{Bool}) = Unknown{Bool}()
@inline Base.:|(x::Bool, ::Unknown{Bool}) = x ? Optional{Bool}(true) : Optional{Bool}()
@inline Base.:|(::Unknown{Bool}, y::Bool) = y ? Optional{Bool}(true) : Optional{Bool}()

@inline Base.:$(::Unknown{Bool}, ::Unknown{Bool}) = Unknown{Bool}()
@inline Base.:$(x::Bool, ::Unknown{Bool}) = Unknown{Bool}()
@inline Base.:$(::Unknown{Bool}, y::Bool) = Unknown{Bool}()

# Basic arrithmetic operators
@inline Base.:-{T<:Number}(::Unknown{T}) = Unknown{promote_op(-,T)}()

@inline Base.:+{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(+,T1,T2)}()
@inline Base.:+{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(+,T1,T2)}()
@inline Base.:+{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(+,T1,T2)}()

@inline Base.:-{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(-,T1,T2)}()
@inline Base.:-{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(-,T1,T2)}()
@inline Base.:-{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(-,T1,T2)}()

@inline Base.:*{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(*,T1,T2)}()
@inline Base.:*{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(*,T1,T2)}()
@inline Base.:*{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(*,T1,T2)}()

@inline Base.:/{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(/,T1,T2)}()
@inline Base.:/{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(/,T1,T2)}()
@inline Base.:/{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(/,T1,T2)}()

@inline Base.://{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(//,T1,T2)}()
@inline Base.://{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(//,T1,T2)}()
@inline Base.://{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(//,T1,T2)}()

@inline Base.:\{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(\,T1,T2)}()
@inline Base.:\{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(\,T1,T2)}()
@inline Base.:\{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(\,T1,T2)}()

@inline Base.:%{T1<:Number, T2<:Number}(::Unknown{T1}, ::Unknown{T2}) = Unknown{promote_op(%,T1,T2)}()
@inline Base.:%{T1<:Number, T2<:Number}(::T1, ::Unknown{T2}) = Unknown{promote_op(%,T1,T2)}()
@inline Base.:%{T1<:Number, T2<:Number}(::Unknown{T1}, ::T2) = Unknown{promote_op(%,T1,T2)}()


text:
"""
    SqExponentialARD(lnℒ::AbstractVector)
    <: RadialKernel{WeightedSqEuclidean}

The squared exponential kernel (aka the Gaussian kernel) with
vector valued length scale
``\\exp(2 \\ln\\mathscr{L}_i) = \\mathscr{L}_i^2 > 0``.

Finding an optimal ``\\mathscr{L}`` is refered to as automatic
relevance detetermination.

For the matrix ``\\Sigma_{ii} = \\mathscr{L}_i^{-2}``

```math
k_{\\mathscr{L}}(x, y) = k_{\\mathscr{L}}(\\lVert x - y\\rVert_\\Sigma^2) =
\\exp\\Big\\lbrace -\\frac{1}{2}\\lVert x - y\\rVert_\\Sigma^2\\Big\\rbrace
```

"""
mutable struct SqExponentialARD{
    F<:AbstractFloat,
    V<:AbstractVector{F},
} <: RadialKernel{WeightedSqEuclidean}
    dist::WeightedSqEuclidean
    lnℒ::V
    diagΣ::V
end

function SqExponentialARD(lnℒ::AbstractVector)
    diagΣ = exp.(-2 .* lnℒ)
    SqExponentialARD(WeightedSqEuclidean(diagΣ), lnℒ, diagΣ)
end

@inline (k::SqExponentialARD)(τ::AbstractFloat) = exp(-τ / 2)

GaussianARD = SqExponentialARD

numparams(k::SqExponentialARD) = (length(k.lnℒ),)
paramtypes(k::SqExponentialARD{F}) where F = (F,)
params(k::SqExponentialARD) = (lnℒ = k.lnℒ,)
function setparams!(k::SqExponentialARD{F}, lnℒ::AbstractVector{F}) where F
    diagΣ = exp.(-2 .* lnℒ)
    copyto!(k.lnℒ, lnℒ)
    copyto!(k.diagΣ, diagΣ)
    k.dist = WeightedSqEuclidean(diagΣ)
end


text:
#!/usr/bin/env julia

using EchogramImages
using Test

a = rand(100,100)
img = imagesc(a)
m,n = size(img)
@test m == 480
@test n == 640


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase, Pkg

name = "Infernal"
version = v"1.1.4"

easel_version = v"0.48"
hmmer_version = v"3.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/EddyRivasLab/infernal/archive/refs/tags/infernal-$(version).tar.gz",
                  "311163c0a21a216e90862df81cd6c13811032076535e44fed4eabdcf45093fea"),
    ArchiveSource("https://github.com/EddyRivasLab/easel/archive/refs/tags/easel-$(easel_version.major).$(easel_version.minor).tar.gz",
                  "c5d055acbe88fa834e81424a15fc5fa54ac787e35f2ea72d4ffd9ea2c1aa29cf"),
    ArchiveSource("https://github.com/EddyRivasLab/hmmer/archive/refs/tags/hmmer-$(hmmer_version).tar.gz",
                  "fab109c67fb8077b32f7907bf07efbc071147be0670aee757c9a3ca7e2d485be")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/infernal-*/

mv ../easel-*/ easel
mv ../hmmer-*/ hmmer

# Update the config.sub from infernal.  Otherwise we get an error when running
# configure: "Invalid configuration `x86_64-linux-musl'."
update_configure_scripts

# generate configure script
autoreconf -vi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-pic --enable-threads --with-gsl

make -j${nproc}
make install

install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# Notes
# - Infernal requires SSE or VMX vector instructions,
#   VMX only on big-endian platforms (ppc64 not ppc64le)
# - ARM vector instruction support coming soon
# - build fails on windows
#   easel.c:39:20: fatal error: syslog.h: No such file or directory
platforms = supported_platforms(; exclude = p -> Sys.iswindows(p) || proc_family(p) != "intel")

# The products that we will ensure are always built
products = [
    ExecutableProduct("cmemit", :cmemit),
    ExecutableProduct("cmconvert", :cmconvert),
    ExecutableProduct("cmfetch", :cmfetch),
    ExecutableProduct("cmpress", :cmpress),
    ExecutableProduct("cmbuild", :cmbuild),
    ExecutableProduct("cmstat", :cmstat),
    ExecutableProduct("cmsearch", :cmsearch),
    ExecutableProduct("cmscan", :cmscan),
    ExecutableProduct("cmcalibrate", :cmcalibrate),
    ExecutableProduct("cmalign", :cmalign),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module FiniteElements
include("Objekte.jl")

# Write your package code here.
export Abc

struct Abc
    hallo::Float64
end

end


text:
if isdefined(Base, :Dates)
    @test JSON.json(Date(2016, 8, 3)) == "\"2016-08-03\""
end

@test JSON.json(:x) == "\"x\""
@test_throws ArgumentError JSON.json(Base)

immutable Type151{T}
    x::T
end

@test JSON.parse(JSON.json(Type151)) == string(Type151)

JSON.lower{T}(v::Type151{T}) = Dict(:type => T, :value => v.x)
@test JSON.parse(JSON.json(Type151(1.0))) == Dict(
    "type" => "Float64",
    "value" => 1.0)

fixednum = Fixed{Int16, 15}(0.1234)
@test JSON.parse(JSON.json(fixednum)) == Float64(fixednum)


text:
module XlibTests

using Xlib

# Example based on tutorial by Philipp K. Janert
# (http://www.linuxjournal.com/article/4879)):

function test1()

    dpy = XOpenDisplay(C_NULL)
    dpy == C_NULL && return 1
    scr = DefaultScreen(dpy)
    white = WhitePixel(dpy, scr)
    black = BlackPixel(dpy, scr)

    println("default screen: ", DefaultScreen(dpy))
    println("protocol version: ", ProtocolVersion(dpy))
    println("protocol revision: ", ProtocolRevision(dpy))
    println("default screen of display: ", ScreenOfDisplay(dpy, scr))
    println("size of screen: ", DisplayWidth(dpy, scr), " x ", DisplayHeight(dpy, scr))
    println("black pixel: ", black)
    println("white pixel: ", white)

    win = XCreateSimpleWindow(dpy,
                              DefaultRootWindow(dpy),
                              50, 50,   # origin
                              200, 200, # size
                              0, black, # border
                              white )   # background
    XMapWindow(dpy, win)

    evt = Ref(XEvent())
    eventMask = StructureNotifyMask
    XSelectInput(dpy, win, eventMask)
    while true
        XNextEvent(dpy, evt) # calls XFlush
        EventType(evt) == MapNotify && break
    end

    gc = XCreateGC(dpy, win, 0, C_NULL)
    XSetForeground(dpy, gc, black)
    XDrawLine(dpy, win, gc, 10, 10,190,190) # from-to
    XDrawLine(dpy, win, gc, 10,190,190, 10) # from-to

    #=
    sleep(0.1)
    XWarpPointer(dpy, DefaultRootWindow(dpy), None, 0, 0, 0, 0, 120, 10)
    for k in 1:30
        sleep(0.1)
        XWarpPointer(dpy, DefaultRootWindow(dpy), None, 0, 0, 0, 0, 3, -1)
    end
    =#

    eventMask = ButtonPressMask|ButtonReleaseMask
    XSelectInput(dpy,win,eventMask) # override previous settings
    while true
        XNextEvent(dpy, evt) # calls XFlush
        EventType(evt) == ButtonRelease && break
    end

    XDestroyWindow(dpy, win)
    XCloseDisplay(dpy)
    return 0
end

# Example based on https://en.wikipedia.org/wiki/Xlib
function test2()
    # Open connection to the server.
    dpy = XOpenDisplay(C_NULL)
    dpy == C_NULL && error("unable to open display")
    scr = DefaultScreen(dpy)

    # Create a window.
    win = XCreateSimpleWindow(dpy, RootWindow(dpy, scr), 10, 10, 200, 100, 1,
                              BlackPixel(dpy, scr), WhitePixel(dpy, scr))

    # Select the kind of events we are interested in.
    XSelectInput(dpy, win, ExposureMask | KeyPressMask)

    # Map (show) the window.
    XMapWindow(dpy, win)

    # Run event loop.
    evt = Ref(XEvent())
    while true
        XNextEvent(dpy, evt)

        # Draw or redraw the window.
        if EventType(evt) == Expose
            XFillRectangle(dpy, win, DefaultGC(dpy, scr), 20, 20, 10, 10)
            XDrawString(dpy, win, DefaultGC(dpy, scr), 50, 50, "Hello, World!")
        end

        # Exit whenever a key is pressed.
        if EventType(evt) == KeyPress
            break
        end
    end

    # Shutdown server connection
    XCloseDisplay(dpy)

    nothing
end

end


text:
using PyCall
gym = pyimport("gym")

# %% swingup functions

function swingup_reset!(env::PyObject)
    env.reset()
    env.env.state += [0.0, 0.0, pi, 0.0]
    env.env.state[3] = (env.env.state[3]+pi)%2pi-pi
    return env.state
end

function swingup_step(env::PyObject, action::Int64, reward_shape::Bool = false)
    new_s, r , done_cart , _ = env.step(action)
    new_s[3] = mod(new_s[3]+pi, 2pi) -pi
    env.env.state = new_s
    if reward_shape
        if done_cart
            thresh = 12*pi/180
            rs_max = 0.1
            r = rs_max * (pi - abs(new_s[3]))/(pi - thresh)
        end
    end
    done = (abs(new_s[1]) > 2.4) | (abs(new_s[4]) > 4pi )
    return new_s, r, done, []
end

# %% tests
#
# env = gym.make("CartPole-v1")
# env.seed(0)
# swingup_reset!(env)
# @show env.state
# new_state, r ,done, _ = swingup_step(env, 1,true)
# @show r
# new_state, r ,done, _ = swingup_step(env, 0, true)
# @show r
# # for j in 1:100
# #     new_state, r ,done, _ = swingup_step(env, Int(((-1)^j+1)/2))
# #     @show done
# #     @show new_state
# # end
#
# env.close()
# env = nothing

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using POMDPModels
using POMDPTesting
using Test

let
    R = [-1. -100 10; -1 10 -100]

    T = zeros(2,3,2)
    T[:,:,1] = [1. 0.5 0.5; 0 0.5 0.5]
    T[:,:,2] = [0. 0.5 0.5; 1 0.5 0.5]

    O = zeros(2,3,2)
    O[:,:,1] = [0.85 0.5 0.5; 0.15 0.5 0.5]
    O[:,:,2] = [0.15 0.5 0.5; 0.85 0.5 0.5]

    pomdp1 = TigerPOMDP()

    pomdp2 = TabularPOMDP(T, R, O, 0.95)

    policy = RandomPolicy(pomdp1, rng=MersenneTwister(2))
    sim = RolloutSimulator(rng=MersenneTwister(3), max_steps=100)

    simulate(sim, pomdp1, policy, updater(policy), initialstate_distribution(pomdp1))

    # test gen(:o, ...)
    o = initialobs(pomdp1, true, MersenneTwister(1))
    @test o == 1
    # test vec
    ov = convert_o(Array{Float64}, true, pomdp1)
    @test ov == [1.]
    o = convert_o(Bool, ov, pomdp1)
    @test o == true

    POMDPTesting.probability_check(pomdp1)
    POMDPTesting.probability_check(pomdp2)
end


text:
using Documenter, PlutoSliderServer, MeshArrays

import OceanStateEstimation
OceanStateEstimation.get_ecco_velocity_if_needed()
MeshArrays.GRID_LL360_download()
MeshArrays.GRID_LLC90_download()
MeshArrays.GRID_CS32_download()

makedocs(
    sitename = "MeshArrays",
    format   = Documenter.HTML(),
    modules  = [MeshArrays],
    pages = [
    "Home" => "index.md",
    "Get Started" => "start.md",
    "Tutorials" => "tutorials.md",
    "Main Features" => "main.md",
    "API documentation" => "API.md",
    "Miscellaneous" => "detail.md",
    "Video Examples" => "videos.md",
    ]
)

lst=("basics.jl","geography.jl","vectors.jl")
pth=("tutorials","tutorials","tutorials")

for ii in 1:length(lst)
    i=lst[ii]
    fil_in=joinpath(@__DIR__,"..", "examples",i)
    fil_out=joinpath(@__DIR__,"build", pth[ii] ,i[1:end-2]*"html")
    PlutoSliderServer.export_notebook(fil_in)
    mv(fil_in[1:end-2]*"html",fil_out)
    cp(fil_in,fil_out[1:end-4]*"jl")
end

deploydocs(
    repo = "github.com/JuliaClimate/MeshArrays.jl",
)


text:
function fib30()
#
#  test polynomial suggested by Goedecker
#
    n = 30;
    p = fib(30)
    z = [-.9644254204388659;
         -.9439446810207413-.1984294062681023*im; -.9439446810207413+.1984294062681023*im;
         -.8833406981383851-.3885717512211873*im; -.8833406981383851+.3885717512211873*im;
         -.7850939960122262-.5624744046156477*im; -.7850939960122262+.5624744046156477*im;
         -.6532262984053521-.7128393846266873*im; -.6532262984053521+.7128393846266873*im;
         -.4931368697016981-.8333156264808382*im; -.4931368697016981+.8333156264808382*im;
         -.3113834310599106-.9187504784857605*im; -.3113834310599106+.9187504784857605*im;
         -.1154175582352257-.9653888511274815*im; -.1154175582352257+.9653888511274815*im;
         .8671260421202705e-1-.9710099019085114*im; .8671260421202705e-1+.9710099019085114*im;
         .2866757308842761-.9349927196966469*im; .2866757308842761+.9349927196966469*im;
         .4761639390711853-.8583045796885356*im; .4761639390711853+.8583045796885356*im;
         .6471568601471165-.7434108410885516*im; .6471568601471165+.7434108410885516*im;
         .7920805342279742-.5941279591518842*im; .7920805342279742+.5941279591518842*im;
         .9037637390940364-.4155315939426088*im; .9037637390940364+.4155315939426088*im;
         .9752028356220179-.2142999451221434*im; .9752028356220179+.2142999451221434*im;
         1.999999999068677];

    z = [z ones(n)]
    p, PolyZeros(z)
end


text:
using Test
using Polyhedra
using JuMP
const MOIB = MOI.Bridges

@testset "Continuous Linear problems with VRepOptimizer" begin
    optimizer = VRepOptimizer{Float64}()
    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "VRep"
    end
    cache = MOIU.UniversalFallback(Polyhedra._MOIModel{Float64}())
    cached = MOIU.CachingOptimizer(cache, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)
    config = MOI.Test.Config(
        exclude=Any[
            MOI.ConstraintBasisStatus,
            MOI.VariableBasisStatus,
            MOI.ConstraintDual,
            MOI.DualObjectiveValue,
            MOI.ObjectiveBound,
        ],
    )
    MOI.Test.runtests(
        bridged,
        config,
        exclude = String[
            "test_attribute_RawStatusString",
            "test_attribute_SolveTimeSec",
            "test_attribute_SolverVersion",
            #   MathOptInterface.jl issue #1431
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
end
@testset "simplex chebyshev center with $T" for T in [Float64, Rational{BigInt}]
    h = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HyperPlane([1, 1], 1)
    lib = DefaultLibrary{T}(VRepOptimizer{T})
    poly = polyhedron(h, lib)
    center, radius = chebyshevcenter(poly)
    @test center ≈ [1/2, 1/2]
    @test radius ≈ 1/2
    @test dim(poly) == 1
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using LightGraphs: AbstractSimpleGraph
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, outneighbors

function adjlist(g::AbstractSimpleGraph)
    N = nv(g)
    el = Vector{Int}[outneighbors(g, i) for i = 1:N]
    return el
end

function adjlist(g::AbstractSimpleWeightedGraph)
    N = nv(g)
    el = Vector{Int}[outneighbors(g, i) for i = 1:N]
    return el
end


text:
@testset "Data" begin
  @testset "Basics" begin
    # underlying table
    d = georef((z=rand(10,10),))
    t = values(d)
    @test size(t) == (100,1)
    @test propertynames(t) == [:z]

    # data equality
    d1 = georef((z=[1,2,3],))
    d2 = georef((z=[1,2,3],))
    d3 = georef((z=[4,5,6],))
    @test d1 == d2
    @test d1 != d3
    @test d2 != d3

    # bidimensional views
    d = georef((z=rand(10,10),))
    v = view(d, 1:3, [:z])
    td = DataFrame(values(d))
    tv = DataFrame(values(v))
    @test nelms(v) == 3
    @test td[1:3,:z] == tv[:,:z]

    # tables API
    x, y, z = rand(100), rand(100), rand(100)
    t = DataFrame(x=x,y=y,z=z)
    d = georef((x=x,y=y,z=z))
    v = view(d, 1:3)
    @test Tables.istable(d) == true
    @test Tables.istable(v) == true
    @test Tables.rowaccess(d) == Tables.rowaccess(t)
    @test Tables.columnaccess(d) == Tables.columnaccess(t)
    @test Tables.rows(d) == Tables.rows(t)
    @test Tables.columns(d) == Tables.columns(t)
    rv = DataFrame(Tables.rows(v))
    rt = DataFrame(Tables.rows(t[1:3,:]))
    cv = DataFrame(Tables.columns(v))
    ct = DataFrame(Tables.columns(t[1:3,:]))
    @test rv == rt
    @test cv == ct
  end

  @testset "GeoTable" begin
    # basic checks
    D = readgeotable(joinpath(datadir,"data3D.tsv"))
    @test variables(D) == (Variable(:value, Float64),)
    @test nelms(D) == 100

    # show methods
    df = DataFrame(x=[1,2,3],y=[4,5,6],z=[1.,2.,3.])
    sdata = georef(df, (:x,:y))
    @test sprint(show, sdata) == "3 SpatialData{Int64,2}"
    @test sprint(show, MIME"text/plain"(), sdata) == "3 PointSet{Int64,2}\n  variables\n    └─z (Float64)"

    if visualtests
      df = DataFrame(x=[25.,50.,75.],y=[25.,75.,50.],z=[1.,0.,1.])
      sdata = georef(df, (:x,:y))
      @plottest plot(sdata) joinpath(datadir,"geodf.png") !isCI
    end
  end

  @testset "PointSet" begin
    # basic checks
    ps = readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t')
    @test variables(ps) == (Variable(:value, Float64),)
    @test nelms(ps) == 100

    # show methods
    ps = georef((z=[1,2,3],), [1. 0. 1.; 0. 1. 1.])
    @test sprint(show, ps) == "3 SpatialData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), ps) == "3 PointSet{Float64,2}\n  variables\n    └─z (Int64)"

    if visualtests
      sdata = georef((z=[1.,0.,1.],), [25. 50. 75.; 25. 75. 50.])
      @plottest plot(sdata) joinpath(datadir,"pset-data.png") !isCI
    end
  end

  @testset "RegularGrid" begin
    # basic checks
    g = georef((Z=rand(100,100),))
    @test variables(g) == (Variable(:Z, Float64),)
    @test nelms(g) == 10000

    # show methods
    g = georef((z=[1,2,3,4],), RegularGrid(2,2))
    @test sprint(show, g) == "4 SpatialData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), g) == "2×2 RegularGrid{Float64,2}\n  variables\n    └─z (Int64)"

    if visualtests
      sdata = georef((z=[1 2; 3 4],))
      @plottest plot(sdata) joinpath(datadir,"grid2D-data1.png") !isCI
      sdata = georef((z=[1 2; 3 4],), (-10.,-10.), (10.,10.))
      @plottest plot(sdata) joinpath(datadir,"grid2D-data2.png") !isCI
    end
  end

  @testset "StructuredGrid" begin
    X = readdlm(joinpath(datadir,"HurricaneX.dat"))
    Y = readdlm(joinpath(datadir,"HurricaneY.dat"))
    P = readdlm(joinpath(datadir,"HurricaneP.dat"))
    g = georef((P=P,), StructuredGrid(X, Y))

    # basic checks
    @test variables(g) == (Variable(:P, Float64),)
    @test nelms(g) == 221*366

    # show methods
    @test sprint(show, g) == "80886 SpatialData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), g) == "221×366 StructuredGrid{Float64,2}\n  variables\n    └─P (Float64)"

    if visualtests
      # TODO
    end
  end
end


text:
mutable struct TreeProblem
    pd::PopulationData
    bt::BinaryTree
    outgroupnode::Int
    nlevels::Int
    model::JuMP.Model
    assign::JuMP.JuMPArray{JuMP.Variable}
    countedge::JuMP.JuMPDict{JuMP.Variable}
    weight::JuMP.JuMPArray{JuMP.Variable}
    f3formula::JuMP.JuMPDict{JuMP.GenericAffExpr{Float64,JuMP.Variable}}
    f3err::JuMP.JuMPDict{JuMP.Variable}
end

function TreeProblem(
    pd::PopulationData,
    bt::BinaryTree;
    binaryencoding::Bool = false,
    nlevels::Int = 1,
    nmixtures::Int = pd.npop,
    solver = Gurobi.GurobiSolver())

    npop = pd.npop
    edges = bt.edges
    leaves = getleaves(bt)
    outgroupnode = leaves[1]
    levels = 1:nlevels

    tree = JuMP.Model(solver=solver)
    if binaryencoding
        @assert nlevels == 1
        JuMP.@variable(tree, assign[1:npop,getleaves(bt),levels] >= 0)
        binaryencodingconstraints(pd, bt, tree, assign)
    else
        JuMP.@variable(tree, assign[1:npop,getleaves(bt),levels], Bin)
    end
    JuMP.@variable(tree, weight[edges] >= 0)
    JuMP.@variable(tree, weightaux[a=1:npop,b=a:npop,edges,levels,levels] >= 0)
    JuMP.@variable(tree, countedge[a=1:npop,b=a:npop,edges,levels,levels] >= 0)
    JuMP.@expression(tree,
        f3formula[a=1:npop,b=a:npop],
        sum(weightaux[a,b,edg,l,m] for edg in edges, l in levels, m in levels)/nlevels^2)
    JuMP.@variable(tree, f3err[a=1:npop,b=a:npop])

    validtreeconstraints(pd, bt, tree, assign, outgroupnode, nlevels, nmixtures)
    countedgeconstraints(pd, bt, tree, assign, countedge, outgroupnode, nlevels)
    errorconstraints(pd, bt, tree, weight, weightaux, countedge, f3formula, f3err, nlevels)

    keytoind = Dict()
    ind = 0
    for key in keys(f3err)
        key[1] == pd.outgroup && continue
        key[2] == pd.outgroup && continue
        ind += 1
        keytoind[key] = ind
    end
    mat = zeros(ind,ind)
    for (a1,b1) in keys(keytoind), (a2,b2) in keys(keytoind)
        mat[keytoind[a1,b1],keytoind[a2,b2]] = pd.cov[a1,b1,a2,b2]
    end
    matinv = inv(mat)

    JuMP.@objective(tree, Min,
        sum(matinv[keytoind[a1,b1],keytoind[a2,b2]]*f3err[a1,b1]*f3err[a2,b2]
            for (a1,b1) in keys(keytoind), (a2,b2) in keys(keytoind)))

    TreeProblem(
        pd, bt, outgroupnode, nlevels,
        tree,
        assign, countedge,
        weight,
        f3formula, f3err)
end

function countedgeconstraints(
    pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    assign::JuMP.JuMPArray{JuMP.Variable},
    countedge::JuMP.JuMPDict{JuMP.Variable},
    outgroupnode::Int,
    nlevels::Int)

    npop = pd.npop
    levels = 1:nlevels

    for (u,v) in bt.edges
        if in((u,v), bt.pathedges[outgroupnode,1])
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] <= 1-sum(assign[a,n,l] for n in getsubtreeleaves(bt,v)))
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] <= 1-sum(assign[b,n,m] for n in getsubtreeleaves(bt,v)))
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] >= 1 - sum(assign[a,n,l] + assign[b,n,m] for n in getsubtreeleaves(bt,v)))
        else
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] <= sum(assign[a,n,l] for n in getsubtreeleaves(bt,v)))
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] <= sum(assign[b,n,m] for n in getsubtreeleaves(bt,v)))
            JuMP.@constraint(tree,
                [a=1:npop, b=a:npop, l=levels, m=levels],
                countedge[a,b,(u,v),l,m] >= sum(assign[a,n,l] + assign[b,n,m] for n in getsubtreeleaves(bt,v)) - 1)
        end
    end

end

function errorconstraints(pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    weight::JuMP.JuMPArray{JuMP.Variable},
    weightaux::JuMP.JuMPDict{JuMP.Variable},
    countedge::JuMP.JuMPDict{JuMP.Variable},
    f3formula::JuMP.JuMPDict{JuMP.GenericAffExpr{Float64,JuMP.Variable}},
    f3err::JuMP.JuMPDict{JuMP.Variable},
    nlevels::Int)

    npop = pd.npop
    bigm = maximum(pd.f3)*2
    levels = 1:nlevels

    # set weight bilinear terms
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,edg=bt.edges, l=levels, m=levels],
        weightaux[a,b,edg,l,m] <= bigm*countedge[a,b,edg,l,m])
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,edg=bt.edges, l=levels, m=levels],
        weightaux[a,b,edg,l,m] <= weight[edg])
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,edg=bt.edges, l=levels, m=levels],
        weightaux[a,b,edg,l,m] >= weight[edg] + bigm*countedge[a,b,edg,l,m] - bigm)

    # set error terms
    JuMP.@constraint(tree, [a=1:npop,b=a:npop],
        f3err[a,b] == pd.f3[a,b] - f3formula[a,b])

end


text:
# # Integration with SimpleSDMLayers

# The `SpatialBoundaries` package works really well with `SimpleSDMLayers`, so
# that you can (i) apply wombling and boundaries finding to a `SimpleSDMLayer`
# object, and (ii) convert the output of a `Womble` object to a *pair* of
# `SimpleSDMLayer` corresponding to the rate and direction of change.

using SimpleSDMLayers
using SpatialBoundaries
using StatsPlots

# Note that the warning about dependencies is a side-effect of loading some
# functionalities for `SimpleSDMLayers` as part of `SpatialBoundaries`, and can
# safely be ignored.

# In this example, we will look at temperature data over North America, and
# measure the rate of change, as well as the direction of it. These data are
# extracted from the WorldClim database.

temperature = SimpleSDMPredictor(
    WorldClim, BioClim, 1; left=-130.0, right=-65.0, bottom=5.0, top=65.0
)

# We can have a look at this layer, after setting a few defaults:

default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)
plot(temperature)

# There is an overload of the `wombling` method for SDM layers, so we can call
# it directly -- this method might result in a bit more memory usage than
# expected, as it requires to transform the `nothing` values into `NaN`s, which
# in turn might require to convert the inner elements of the layer grid.

W = wombling(temperature)

# By default, this returns a `LatticeWomble`. Let's look at the direction of
# change -- mapping this information is difficult, so we will focus on the

stephist(
    deg2rad.(sort(vec(W.θ)));
    proj=:polar,
    lab="",
    c=:teal,
    fill=(0, 0.2, :teal),
    nbins=100,
)

# The values for the direction of change are concentrated around 180° - note
# that the direction of change is reported as a wind direction, meaning that
# values tend to increase on a south-north axis.

# We can also map the rate of change. This is far easier to do with a proper SDM
# layer, so we will convert the wombling output:

Lr, Ld = SimpleSDMPredictor(W)

# Note that we do not use `convert` here, because this call returns *two* layers
# in a tuple -- this is a slight deviation from what we expect with
# `SimpleSDMLayers`, but it makes the code a little easier to write, and so is
# considered an acceptable trade-off.

plot(Lr)

# With the `Ld` layer storing the direction of change, we *can* do a map - the
# main issue with mapping directions is colors, which need to have the same
# begin and endpoint, while remaining accessible to people with color vision
# defficiencies. [Cyclic color schemes](https://www.fabiocrameri.ch/colourmaps/)
# work for this purpose:

plot(Ld, c=:brocO, clim=(0, 360))


text:
## Exercise 8-10
## The following functions are all intended to check whether a string contains any lowercase letters, but at least some of them are wrong. For each function, describe what the function actually does (assuming that the parameter is a string).

"""
if first letter is lowercase return true, else return false;
"""
function anylowercase1(s)
    for c in s
        if islowercase(c)
            return true
        else
            return false
        end
    end
end

"""
return string true if 'c' is lowercase, i.e. will always return true
"""
function anylowercase2(s)
    for c in s
        if islowercase('c')
            return "true"
        else
            return "false"
        end
    end
end

"""
returns true if last char is lowercase, else false
"""
function anylowercase3(s)
    for c in s
        flag = islowercase(c)
    end
    flag
end

"""
yeilds expected results
"""
function anylowercase4(s)
    flag = false
    for c in s
        flag = flag || islowercase(c)
    end
    flag
end

"""
will return true if every char is lowercase, else false
"""
function anylowercase5(s)
    for c in s
        if !islowercase(c)
            return false
        end
    end
    true
end


text:
module NMRPulses

export Pulse,
    BasicPulse,
    IdlePulse,
    Sequence,
    PiPulse,
    PiHalfPulse,
    ArbitraryPulse,
    PWPulse,
    prepare_sequence

include("pulses.jl")
include("tools.jl")

end