text:
module SingleRoomModule

import LinearAlgebra as LA
import MiniFB as MFB
import Random
import ..RayCastWorlds as RCW
import RayCaster as RC
import ReinforcementLearningBase as RLBase
import SimpleDraw as SD
import StaticArrays as SA

#####
##### game logic
#####

const NUM_OBJECTS = 2
const WALL = 1
const GOAL = 2
const NUM_ACTIONS = 4

mutable struct SingleRoomWorld{T, RNG, R}
    tile_map::BitArray{3}
    num_directions::Int
    player_position_wu::SA.SVector{2, T}
    player_direction_au::Int
    player_radius_wu::T
    position_increment_wu::T
    directions_wu::Vector{SA.SVector{2, T}}
    ray_stop_position_tu::Array{Int, 2}
    ray_hit_dimension::Array{Int, 1}
    ray_distance_wu::Array{T, 1}
    goal_position::CartesianIndex{2}
    rng::RNG
    reward::R
    goal_reward::R
    done::Bool
    semi_field_of_view_wu::T
    num_rays::Int
    ray_directions_wu::Vector{SA.SVector{2, T}}
end

function SingleRoomWorld(;
        T = Float32,
        height_tile_map_tu = 8,
        width_tile_map_tu = 16,
        num_directions = 128, # angles go from 0 to num_directions - 1 (0 corresponding to positive x-axes)
        player_radius_wu = convert(T, 1 / 8), # should be less than 0.5
        position_increment_wu = convert(T, 1 / 8),
        rng = Random.GLOBAL_RNG,
        R = Float32,
        semi_field_of_view_wu = convert(T, 2/3),
        num_rays = 512,
    )

    tile_map = falses(NUM_OBJECTS, height_tile_map_tu, width_tile_map_tu)

    tile_map[WALL, :, 1] .= true
    tile_map[WALL, :, width_tile_map_tu] .= true
    tile_map[WALL, 1, :] .= true
    tile_map[WALL, height_tile_map_tu, :] .= true

    goal_position = CartesianIndex(rand(rng, 2 : height_tile_map_tu - 1), rand(rng, 2 : width_tile_map_tu - 1))
    tile_map[GOAL, goal_position] = true

    directions_wu = Array{SA.SVector{2, T}}(undef, num_directions)
    for i in 1:num_directions
        theta_wu = (i - 1) * 2 * pi / num_directions
        directions_wu[i] = SA.SVector(convert(T, cos(theta_wu)), convert(T, sin(theta_wu)))
    end

    player_position_tu = RCW.sample_empty_position(rng, tile_map)
    player_position_wu = SA.SVector(convert(T, player_position_tu[1] - 0.5), convert(T, player_position_tu[2] - 0.5))

    player_direction_au = rand(rng, 0 : num_directions - 1)

    ray_directions_wu = Array{SA.SVector{2, T}}(undef, num_rays)
    ray_stop_position_tu = Array{Int}(undef, 2, num_rays)
    ray_hit_dimension = Array{Int}(undef, num_rays)
    ray_distance_wu = Array{T}(undef, num_rays)

    reward = zero(R)
    goal_reward = one(R)
    done = false

    world = SingleRoomWorld(tile_map,
                       num_directions,
                       player_position_wu,
                       player_direction_au,
                       player_radius_wu,
                       position_increment_wu,
                       directions_wu,
                       ray_stop_position_tu,
                       ray_hit_dimension,
                       ray_distance_wu,
                       goal_position,
                       rng,
                       reward,
                       goal_reward,
                       done,
                       semi_field_of_view_wu,
                       num_rays,
                       ray_directions_wu,
                      )

    RCW.reset!(world)

    return world
end

function RCW.reset!(world::SingleRoomWorld{T}) where {T}
    tile_map = world.tile_map
    rng = world.rng
    player_radius_wu = world.player_radius_wu
    goal_position = world.goal_position
    num_directions = world.num_directions
    _, height_tile_map_tu, width_tile_map_tu = size(tile_map)

    tile_map[GOAL, goal_position] = false

    new_goal_position = CartesianIndex(rand(rng, 2 : height_tile_map_tu - 1), rand(rng, 2 : width_tile_map_tu - 1))
    world.goal_position = new_goal_position
    tile_map[GOAL, new_goal_position] = true

    new_player_position_tu = RCW.sample_empty_position(rng, tile_map)
    new_player_position_wu = SA.SVector(convert(T, new_player_position_tu[1] - 0.5), convert(T, new_player_position_tu[2] - 0.5))
    world.player_position_wu = new_player_position_wu

    new_player_direction_au = rand(rng, 0 : num_directions - 1)
    world.player_direction_au = new_player_direction_au

    world.reward = zero(world.reward)
    world.done = false

    RCW.cast_rays!(world)

    return nothing
end

function RCW.act!(world::SingleRoomWorld, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action: $(action)"

    tile_map = world.tile_map
    player_direction_au = world.player_direction_au
    player_position_wu = world.player_position_wu
    player_radius_wu = world.player_radius_wu
    num_directions = world.num_directions
    num_rays = length(world.ray_directions_wu)
    goal_map = @view tile_map[GOAL, :, :]

    if action in Base.OneTo(2)
        directions_wu = world.directions_wu
        position_increment_wu = world.position_increment_wu
        player_direction_wu = directions_wu[player_direction_au + 1]
        wall_map = @view tile_map[WALL, :, :]

        if action == 1
            new_player_position_wu = RCW.move_forward(player_position_wu, player_direction_wu, position_increment_wu)
        else
            new_player_position_wu = RCW.move_backward(player_position_wu, player_direction_wu, position_increment_wu)
        end

        is_colliding_with_goal = RCW.is_player_colliding(goal_map, new_player_position_wu, player_radius_wu)
        is_colliding_with_wall = RCW.is_player_colliding(wall_map, new_player_position_wu, player_radius_wu)

        if is_colliding_with_goal || is_colliding_with_wall
            if is_colliding_with_goal
                world.reward = world.goal_reward
                world.done = true
            else
                world.reward = zero(world.reward)
                world.done = false
            end
        else
            world.player_position_wu = new_player_position_wu
            world.reward = zero(world.reward)
            world.done = false
        end
    else
        if action == 3
            new_player_direction_au = RCW.turn_left(player_direction_au, num_directions)
        else
            new_player_direction_au = RCW.turn_right(player_direction_au, num_directions)
        end

        world.player_direction_au = new_player_direction_au
        world.reward = zero(world.reward)
        world.done = false
    end

    return nothing
end

rotate_minus_90(vec::SA.SVector{2}) = typeof(vec)(vec[2], -vec[1])

function RCW.cast_rays!(world::SingleRoomWorld)
    tile_map = world.tile_map
    player_direction_au = world.player_direction_au
    player_position_wu = world.player_position_wu
    num_rays = length(world.ray_directions_wu)
    directions_wu = world.directions_wu
    num_directions = world.num_directions
    ray_stop_position_tu = world.ray_stop_position_tu
    ray_hit_dimension = world.ray_hit_dimension
    ray_distance_wu = world.ray_distance_wu
    ray_directions_wu = world.ray_directions_wu
    semi_field_of_view_wu = world.semi_field_of_view_wu

    _, height_tile_map_tu, width_tile_map_tu = size(tile_map)
    obstacle_map = @view any(tile_map, dims = 1)[1, :, :]
    field_of_view_start_au = player_direction_au - (num_rays - 1) ÷ 2
    field_of_view_end_au = field_of_view_start_au + num_rays - 1

    num_rays = length(ray_directions_wu)
    player_direction_wu = directions_wu[player_direction_au + 1]
    camera_direction_wu = rotate_minus_90(player_direction_wu)
    first_ray_direction_wu = player_direction_wu + semi_field_of_view_wu * camera_direction_wu
    last_ray_direction_wu = player_direction_wu - semi_field_of_view_wu * camera_direction_wu
    unnormalized_ray_directions_wu_range = range(first_ray_direction_wu, last_ray_direction_wu, length = num_rays)

    for i in 1:num_rays
        ray_direction_wu = LA.normalize(unnormalized_ray_directions_wu_range[i])
        ray_directions_wu[i] = ray_direction_wu
        i_hit_tu, j_hit_tu, hit_dimension, side_dist_wu = RC.cast_ray(obstacle_map, player_position_wu..., ray_direction_wu...)
        ray_stop_position_tu[1, i] = i_hit_tu
        ray_stop_position_tu[2, i] = j_hit_tu
        ray_hit_dimension[i] = hit_dimension
        ray_distance_wu[i] = side_dist_wu
    end

    return nothing
end

#####
##### drawing & playing
#####

const NUM_VIEWS = 2
const CAMERA_VIEW = 1
const TOP_VIEW = 2

struct SingleRoom{T, RNG, R, C} <: RCW.AbstractGame
    world::SingleRoomWorld{T, RNG, R}
    top_view::Array{C, 2}
    camera_view::Array{C, 2}
    tile_map_colors::NTuple{NUM_OBJECTS + 1, C}
    ray_color::C
    player_color::C
    floor_color::C
    ceiling_color::C
    wall_dim_1_color::C
    wall_dim_2_color::C
    goal_dim_1_color::C
    goal_dim_2_color::C
    camera_height_tile_wu::T
    height_camera_view_pu::Int
end

function SingleRoom(;
        T = Float32,
        height_tile_map_tu = 8,
        width_tile_map_tu = 16,
        num_directions = 128,
        player_radius_wu = convert(T, 1 / 8),
        position_increment_wu = convert(T, 1 / 8),
        rng = Random.GLOBAL_RNG,
        R = Float32,
        semi_field_of_view_wu = convert(T, 2/3),
        num_rays = 512,
        pu_per_tu = 32,
        camera_height_tile_wu = convert(T, 1),
        height_camera_view_pu = 256,
    )

    C = UInt32

    world = SingleRoomWorld(T = T,
                           height_tile_map_tu = height_tile_map_tu,
                           width_tile_map_tu = width_tile_map_tu,
                           num_directions = num_directions,
                           player_radius_wu = player_radius_wu,
                           position_increment_wu = position_increment_wu,
                           rng = rng,
                           R = R,
                           semi_field_of_view_wu = semi_field_of_view_wu,
                           num_rays = num_rays,
                          )

    tile_map_colors = (0x00FFFFFF, 0x00FF0000, 0x00000000)
    ray_color = 0x00808080
    player_color = 0x00c0c0c0
    floor_color = 0x00404040
    ceiling_color = 0x00FFFFFF
    wall_dim_1_color = 0x00808080
    wall_dim_2_color = 0x00c0c0c0
    goal_dim_1_color = 0x00800000
    goal_dim_2_color = 0x00c00000

    RCW.cast_rays!(world)

    camera_view = Array{C}(undef, height_camera_view_pu, num_rays)

    top_view = Array{C}(undef, height_tile_map_tu * pu_per_tu, width_tile_map_tu * pu_per_tu)

    env = SingleRoom(world,
                          top_view,
                          camera_view,
                          tile_map_colors,
                          ray_color,
                          player_color,
                          floor_color,
                          ceiling_color,
                          wall_dim_1_color,
                          wall_dim_2_color,
                          goal_dim_1_color,
                          goal_dim_2_color,
                          camera_height_tile_wu,
                          height_camera_view_pu,
                         )

    RCW.update_camera_view!(env)
    RCW.update_top_view!(env)

    return env
end

function RCW.reset!(env::SingleRoom)
    RCW.reset!(env.world)
    RCW.update_top_view!(env)
    RCW.update_camera_view!(env)
    return nothing
end

function RCW.act!(env::SingleRoom, action)
    world = env.world
    RCW.act!(world, action)
    RCW.cast_rays!(world)
    RCW.update_top_view!(env)
    RCW.update_camera_view!(env)
    return nothing
end

function draw_tile_map!(top_view, tile_map, colors)
    _, height_tile_map_tu, width_tile_map_tu = size(tile_map)
    height_top_view_pu, width_top_view_pu = size(top_view)

    pu_per_tu = height_top_view_pu ÷ height_tile_map_tu

    for j in 1:width_tile_map_tu
        for i in 1:height_tile_map_tu
            i_top_left = (i - 1) * pu_per_tu + 1
            j_top_left = (j - 1) * pu_per_tu + 1

            shape = SD.FilledRectangle(SD.Point(i_top_left, j_top_left), pu_per_tu, pu_per_tu)

            object = findfirst(@view tile_map[:, i, j])
            if isnothing(object)
                color = colors[end]
            else
                color = colors[object]
            end

            SD.draw!(top_view, shape, color)

            top_view[i_top_left, j_top_left : j_top_left + pu_per_tu - 1] .= 0x00cccccc
            top_view[i_top_left + pu_per_tu - 1, j_top_left : j_top_left + pu_per_tu - 1] .= 0x00cccccc
            top_view[i_top_left : i_top_left + pu_per_tu - 1, j_top_left] .= 0x00cccccc
            top_view[i_top_left : i_top_left + pu_per_tu - 1, j_top_left + pu_per_tu - 1] .= 0x00cccccc
        end
    end

    return nothing
end

function RCW.update_camera_view!(env::SingleRoom)
    world = env.world
    camera_view = env.camera_view
    floor_color = env.floor_color
    ceiling_color = env.ceiling_color
    wall_dim_1_color = env.wall_dim_1_color
    wall_dim_2_color = env.wall_dim_2_color
    goal_dim_1_color = env.goal_dim_1_color
    goal_dim_2_color = env.goal_dim_2_color
    camera_height_tile_wu = env.camera_height_tile_wu

    tile_map = world.tile_map
    player_direction_au = world.player_direction_au
    player_position_wu = world.player_position_wu
    ray_directions_wu = world.ray_directions_wu
    num_rays = length(ray_directions_wu)
    num_directions = world.num_directions
    ray_stop_position_tu = world.ray_stop_position_tu
    ray_hit_dimension = world.ray_hit_dimension
    ray_distance_wu = world.ray_distance_wu
    directions_wu = world.directions_wu
    semi_field_of_view_wu = world.semi_field_of_view_wu

    _, height_tile_map_tu, width_tile_map_tu = size(tile_map)
    height_camera_view_pu, width_camera_view_pu = size(camera_view)

    player_direction_wu = directions_wu[player_direction_au + 1]
    for i in 1:num_rays
        ray_direction_wu = ray_directions_wu[i]

        projected_distance_wu = ray_distance_wu[i] * sum(player_direction_wu .* ray_direction_wu)

        height_line = camera_height_tile_wu * num_rays / (2 * semi_field_of_view_wu * projected_distance_wu)
        if isfinite(height_line)
            height_line_pu = floor(Int, height_line)
        else
            height_line_pu = height_camera_view_pu
        end

        hit_dimension = ray_hit_dimension[i]
        ray_stop_position_i_tu = ray_stop_position_tu[1, i]
        ray_stop_position_j_tu = ray_stop_position_tu[2, i]

        if tile_map[WALL, ray_stop_position_i_tu, ray_stop_position_j_tu]
            if hit_dimension == 1
                color = wall_dim_1_color
            else hit_dimension == 2
                color = wall_dim_2_color
            end
        else
            if hit_dimension == 1
                color = goal_dim_1_color
            else hit_dimension == 2
                color = goal_dim_2_color
            end
        end

        k = width_camera_view_pu - i + 1

        if height_line_pu >= height_camera_view_pu - 1
            camera_view[:, k] .= color
        else
            padding_pu = (height_camera_view_pu - height_line_pu) ÷ 2
            camera_view[1:padding_pu, k] .= ceiling_color
            camera_view[padding_pu + 1 : end - padding_pu, k] .= color
            camera_view[end - padding_pu + 1 : end, k] .= floor_color
        end
    end

    return nothing
end

function RCW.update_top_view!(env::SingleRoom)
    world = env.world
    top_view = env.top_view
    tile_map_colors = env.tile_map_colors
    ray_color = env.ray_color
    player_color = env.player_color
    tile_map = world.tile_map
    player_direction_au = world.player_direction_au
    player_position_wu = world.player_position_wu
    ray_directions_wu = world.ray_directions_wu
    num_rays = length(ray_directions_wu)
    num_directions = world.num_directions
    player_radius_wu = world.player_radius_wu
    ray_stop_position_tu = world.ray_stop_position_tu
    ray_hit_dimension = world.ray_hit_dimension
    ray_distance_wu = world.ray_distance_wu
    directions_wu = world.directions_wu

    _, height_tile_map_tu, width_tile_map_tu = size(tile_map)
    height_top_view_pu, width_top_view_pu = size(top_view)

    pu_per_tu = height_top_view_pu ÷ height_tile_map_tu

    i_player_position_pu, j_player_position_pu = RCW.wu_to_pu.(player_position_wu, pu_per_tu)
    player_radius_pu = RCW.wu_to_pu(player_radius_wu, pu_per_tu)

    draw_tile_map!(top_view, tile_map, tile_map_colors)

    for i in 1:num_rays
        ray_direction_wu = ray_directions_wu[i]
        i_ray_stop_pu, j_ray_stop_pu = RCW.wu_to_pu.(player_position_wu + ray_distance_wu[i] * ray_direction_wu, pu_per_tu)
        SD.draw!(top_view, SD.Line(SD.Point(i_player_position_pu, j_player_position_pu), SD.Point(i_ray_stop_pu, j_ray_stop_pu)), ray_color)
    end

    SD.draw!(top_view, SD.Circle(SD.Point(i_player_position_pu - player_radius_pu, j_player_position_pu - player_radius_pu), 2 * player_radius_pu + 1), player_color)

    return nothing
end

RCW.get_action_keys(env::SingleRoom) = (MFB.KB_KEY_W, MFB.KB_KEY_S, MFB.KB_KEY_A, MFB.KB_KEY_D)
RCW.get_action_names(env::SingleRoom) = (:MOVE_FORWARD, :MOVE_BACKWARD, :TURN_LEFT, :TURN_RIGHT)

function RCW.play!(game::SingleRoom)
    world = game.world
    top_view = game.top_view
    camera_view = game.camera_view
    tile_map_colors = game.tile_map_colors
    ray_color = game.ray_color
    player_color = game.player_color
    floor_color = game.floor_color
    ceiling_color = game.ceiling_color
    wall_dim_1_color = game.wall_dim_1_color
    wall_dim_2_color = game.wall_dim_2_color

    height_top_view_pu, width_top_view_pu = size(top_view)
    height_camera_view_pu, width_camera_view_pu = size(camera_view)

    height_image = max(height_top_view_pu, height_camera_view_pu)
    width_image = max(width_top_view_pu, width_camera_view_pu)

    frame_buffer = zeros(UInt32, width_image, height_image)

    window = MFB.mfb_open(String(nameof(typeof(game))), width_image, height_image)

    action_keys = RCW.get_action_keys(game)

    current_view = CAMERA_VIEW
    if current_view == CAMERA_VIEW
        RCW.copy_image_to_frame_buffer!(frame_buffer, camera_view)
    elseif current_view == TOP_VIEW
        RCW.copy_image_to_frame_buffer!(frame_buffer, top_view)
    end

    steps_taken = 0

    function keyboard_callback(window, key, mod, is_pressed)::Cvoid
        if is_pressed
            println("*******************************")
            @show key

            if key == MFB.KB_KEY_Q
                MFB.mfb_close(window)
                return nothing
            elseif key == MFB.KB_KEY_R
                RCW.reset!(game)
                steps_taken = 0
            elseif key == MFB.KB_KEY_V
                current_view = mod1(current_view + 1, NUM_VIEWS)
                fill!(frame_buffer, 0x00000000)
            elseif key in action_keys
                action = findfirst(==(key), action_keys)
                RCW.act!(game, action)
                steps_taken += 1
            else
                @warn "No keybinding exists for $(key)"
            end

            if current_view == CAMERA_VIEW
                RCW.copy_image_to_frame_buffer!(frame_buffer, camera_view)
            elseif current_view == TOP_VIEW
                RCW.copy_image_to_frame_buffer!(frame_buffer, top_view)
            end

            @show steps_taken
            @show world.reward
            @show world.done
        end

        return nothing
    end

    MFB.mfb_set_keyboard_callback(window, keyboard_callback)

    while MFB.mfb_wait_sync(window)
        state = MFB.mfb_update(window, frame_buffer)

        if state != MFB.STATE_OK
            break;
        end
    end

    return nothing
end

#####
##### RLBase API
#####

RLBase.StateStyle(env::RCW.RLBaseEnv{E}) where {E <: SingleRoom} = RLBase.Observation{Any}()
RLBase.state_space(env::RCW.RLBaseEnv{E}, ::RLBase.Observation) where {E <: SingleRoom} = nothing
RLBase.state(env::RCW.RLBaseEnv{E}, ::RLBase.Observation) where {E <: SingleRoom} = env.env.camera_view

RLBase.reset!(env::RCW.RLBaseEnv{E}) where {E <: SingleRoom} = RCW.reset!(env.env)

RLBase.action_space(env::RCW.RLBaseEnv{E}) where {E <: SingleRoom} = Base.OneTo(NUM_ACTIONS)
(env::RCW.RLBaseEnv{E})(action) where {E <: SingleRoom} = RCW.act!(env.env, action)

RLBase.reward(env::RCW.RLBaseEnv{E}) where {E <: SingleRoom} = env.env.world.reward
RLBase.is_terminated(env::RCW.RLBaseEnv{E}) where {E <: SingleRoom} = env.env.world.done

end # module


text:
using Gen

##################################################
# mixture of uniform and gaussians 2D coordinate #
##################################################

struct UniformGaussianMixtureCoord <: Distribution{Vector{Float64}}
end

const mixture_of_gaussians = UniformGaussianMixtureCoord()

function (::UniformGaussianMixtureCoord)(uniform_prob, normal_means, normal_std)
    return random(mixture_of_gaussians, uniform_prob, normal_means, normal_std)
end

function Gen.random(
        ::UniformGaussianMixtureCoord,
        uniform_prob::Float64, normal_means::Vector{Tuple{Float64,Float64}}, normal_std::Float64)
    if bernoulli(uniform_prob)
        return uniform_coord()
    else
        i = uniform_discrete(1, length(normal_means))
        pt = normal_means[i]
        @assert length(pt) == 2
        x = normal(pt[1], normal_std)
        y = normal(pt[2], normal_std)
        return Float64[x, y]
    end
end

function Gen.logpdf(
        ::UniformGaussianMixtureCoord, pt::Vector{Float64},
        uniform_prob::Float64, normal_means::Vector{Tuple{Float64,Float64}}, normal_std::Float64)
    @assert length(pt) == 2
    ls = Vector{Float64}(undef, length(normal_means) + 1)
    for i in 1:length(normal_means)
        @assert length(normal_means[i]) == 2
        mu_x, mu_y = normal_means[i]
        ls[i] = log(1 - uniform_prob) - log(length(normal_means))
        ls[i] += logpdf(normal, pt[1], mu_x, normal_std) + logpdf(normal, pt[2], mu_y, normal_std)
    end
    ls[end] = log(uniform_prob)
    return logsumexp(ls)
end

#########################
# uniform 2D coordinate #
#########################

struct UniformCoord <: Distribution{Vector{Float64}}
end

const uniform_coord = UniformCoord()

function (::UniformCoord)()
    return random(uniform_coord)
end

function Gen.random(::UniformCoord)
    x = uniform(0, 1)
    y = uniform(0, 1)
    return Float64[x, y]
end

function Gen.logpdf(::UniformCoord, pt::Vector{Float64})
    @assert length(pt) == 2
    if (0.0 <= pt[1] <= 1.0) && (0.0 <= pt[2] <= 1.0)
        return 0.0
    else
        return -Inf
    end
end

############################################
# uniform 2d coordinate within a rectangle #
############################################

struct UniformCoordRect <: Distribution{Vector{Float64}}
end

const uniform_coord_rect = UniformCoordRect()

function (::UniformCoordRect)()
    return random(uniform_coord_rect)
end

function Gen.random(::UniformCoordRect, xmin, xmax, ymin, ymax)
    x = uniform(xmin, xmax)
    y = uniform(ymin, ymax)
    return Float64[x, y]
end

function Gen.logpdf(::UniformCoordRect, pt::Vector{Float64}, xmin, xmax, ymin, ymax)
    @assert length(pt) == 2
    if (xmin <= pt[1] <= xmax) && (ymin <= pt[2] <= ymax)
        return -log(xmax - xmin) - log(ymax - ymin)
    else
        return -Inf
    end
end

Gen.has_argument_grads(::UniformCoordRect) = (false, false, false, false)
Gen.has_output_grad(::UniformCoordRect) = false

export mixture_of_gaussians, uniform_coord, uniform_coord_rect

ext: jl


text:
using NLOptControl
using Base.Test

# NOTE: for some reason if a test fails, then later tests are effected.
# This behaviour can be very perplexing.
# If a test fails ignore all remaining errors they may be erroneous!
# Either fix the broken test or put it at the bottom of the stack.
const tol = 5e-2
const big_tol = 0.5 # can reduce if the number of points are increased
 const integrationConfigs = [:lgrExplicit,:lgrImplicit,:trapezoidal,:bkwEuler]
#const integrationConfigs = [:lgrExplicit]

include("ocp.jl")
include("mpc.jl")


text:
@testset "smo svdd" begin
    dummy_data, labels = generate_mvn_with_outliers(2, 1000, 123, true, true)
    pools = fill(:U, size(dummy_data, 2))

    init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(0.5), 0.5)

    model = SVDD.VanillaSVDD(dummy_data)
    SVDD.initialize!(model, init_strategy)
    solver = SMOSolver(1e-4, 10)
    status = SVDD.solve!(model, solver)
    @test status == :UserLimit

    solver = SMOSolver(1e-2, 1000)
    status = SVDD.fit!(model, solver)
    @test status == JuMP.MathOptInterface.OPTIMAL
    @test model.R > 0
end


text:
using MPIReco
using Test

@testset "test saving and loading of MPI images to and from MDF's" begin
  bSF = MPIFile(joinpath(datadir, "calibrations", "12.mdf"))
  b = MPIFile(joinpath(datadir, "measurements", "20211226_203916_MultiPatch", "1.mdf"))
  r = defaultRecoParams()
  r[:measPath] = filepath(b)
  r[:SFPath] = filepath(bSF)
  r[:frames] =  1
  r[:minFreq] = 80e3
  r[:SNRThresh] = 4
  r[:lambd] = 0.001
  r[:iterations] = 1

  # standard reconstruction
  c = reconstruction(r)
  names =  axisnames(c)
  values =  axisvalues(c)

  # save c as MDF
  saveRecoData(joinpath(imgdir, "reco.mdf"),c)

  # load MDF
  cmdf = loadRecoData(joinpath(imgdir, "reco.mdf"))
  @test axisnames(cmdf) == names
  @test axisvalues(cmdf) == values
end


text:
module CCBladeOpenMDAOExample

include("derivatives.jl")
include("openmdao.jl")

end # module


text:
"""
    submat(A::AbstractPDMat, inds)

Select a square submatrix of an `AbstractPDMat` specified by `inds`.
Returns a `AbstractPDMat` of the same type.
"""
submat(A::PDMat, inds) = PDMat(A[inds, inds])
submat(A::PSDMat, inds) = PSDMat(A[inds, inds])
submat(A::PDiagMat, inds) = PDiagMat(A.diag[inds])
submat(A::PDSparseMat, inds) = PDSparseMat(sparse(A[inds, inds]))
submat(A::WoodburyPDMat, inds) = WoodburyPDMat(A.A[inds, :], A.D, Diagonal(A.S.diag[inds]))
function submat(A::ScalMat, inds)
    # inds should not be able to "access" non-existent indices
    checkbounds(Bool, A, inds) || throw(BoundsError(A, inds))
    return ScalMat(length(inds), A.value)
end


text:
push!(LOAD_PATH,"./")
using Fractal
using Images

# fspecstr="0x5c82f6 0x720bd8 0x417b14 (0.303251,-0.0213448) zoom2.08799e-05 pow0.79505 affine1"
fspecstr="0x9f8f41 0x3f1e83 0xd3f1e8 (0.0,0.0) zoom3.0 pow0.79505 affine1"
camspec=parse(CameraSpec,fspecstr)
imagespec=ImageSpec(2^9,2^9,2)
for fracname in fractal_names()
	println("Generating $fracname...")
	filename="test_$(fracname).png"
	if isfile(filename)
		continue
	end
	fracspec=FractalSpec(create_fractal(fracname))
	fracgen=FractalGenerator(camspec,imagespec,fracspec)
	@time img=generate_fractal(fracgen)
	save(filename,img)
end

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal


text:
# dm[SingleItemIndex, SingleItemIndex)
function Base.getindex(d::DataMatrix,
                       i::SingleIndex,
                       j::SingleIndex)
    if d.na[i, j]
        return NA
    else
        return d.data[i, j]
    end
end

# dm[SingleItemIndex, MultiItemIndex]
function Base.getindex(x::DataMatrix,
                       i::SingleIndex,
                       col_inds::AbstractDataVector{Bool})
    getindex(x, i, find(col_inds))
end
function Base.getindex(x::DataMatrix,
                       i::SingleIndex,
                       col_inds::AbstractDataVector)
    getindex(x, i, dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       i::SingleIndex,
                       col_inds::MultiIndex)
    DataArray(x.data[i, col_inds], x.na[i, col_inds])
end

# dm[MultiItemIndex, SingleItemIndex]
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector{Bool},
                       j::SingleIndex)
    getindex(x, find(row_inds), j)
end
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector,
                       j::SingleIndex)
    getindex(x, dropna(row_inds), j)
end
# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::MultiIndex,
                       j::SingleIndex)
    DataArray(x.data[row_inds, j], x.na[row_inds, j])
end

# dm[MultiItemIndex, MultiItemIndex]
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector{Bool},
                       col_inds::AbstractDataVector{Bool})
    return getindex(x,
                    find(row_inds),
                    find(col_inds))
end
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector{Bool},
                       col_inds::AbstractDataVector)
    return getindex(x,
                    find(row_inds),
                    dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector{Bool},
                       col_inds::MultiIndex)
    return getindex(x,
                    find(row_inds),
                    col_inds)
end
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector,
                       col_inds::AbstractDataVector{Bool})
    return getindex(x,
                    dropna(row_inds),
                    find(col_inds))
end
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector,
                       col_inds::AbstractDataVector)
    return getindex(x,
                    dropna(row_inds),
                    dropna(col_inds))
end

# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::AbstractDataVector,
                       col_inds::MultiIndex)
    return getindex(x, dropna(row_inds), col_inds)
end

# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::MultiIndex,
                       col_inds::AbstractDataVector{Bool})
    return getindex(x,
                    row_inds,
                    find(col_inds))
end

# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::MultiIndex,
                       col_inds::AbstractDataVector)
    return getindex(x, row_inds, dropna(col_inds))
end

# TODO: Make inds::AbstractVector
function Base.getindex(x::DataMatrix,
                       row_inds::MultiIndex,
                       col_inds::MultiIndex)
    return DataArray(x.data[row_inds, col_inds],
                     x.na[row_inds, col_inds])
end

# dm[SingleItemIndex, SingleItemIndex] = NA
function Base.setindex!(dm::DataMatrix,
                        val::NAtype,
                        i::SingleIndex,
                        j::SingleIndex)
    dm.na[i, j] = true
    return dm
end

# dm[SingleItemIndex, SingleItemIndex] = Single Item
function Base.setindex!(dm::DataMatrix,
                        val::Any,
                        i::SingleIndex,
                        j::SingleIndex)
    dm.data[i, j] = val
    dm.na[i, j] = false
    return dm
end

# dm[MultiItemIndex, SingleItemIndex] = NA
function Base.setindex!(dm::DataMatrix,
                        val::NAtype,
                        row_inds::MultiIndex,
                        j::SingleIndex)
    dm.na[row_inds, j] = true
    return dm
end

# dm[MultiItemIndex, SingleItemIndex] = Multiple Items
function Base.setindex!{S, T}(dm::DataMatrix{S},
                              vals::Vector{T},
                              row_inds::MultiIndex,
                              j::SingleIndex)
    dm.data[row_inds, j] = vals
    dm.na[row_inds, j] = false
    return dm
end

# dm[MultiItemIndex, SingleItemIndex] = Single Item
function Base.setindex!(dm::DataMatrix,
                        val::Any,
                        row_inds::MultiIndex,
                        j::SingleIndex)
    dm.data[row_inds, j] = val
    dm.na[row_inds, j] = false
    return dm
end

# dm[SingleItemIndex, MultiItemIndex] = NA
function Base.setindex!(dm::DataMatrix,
                        val::NAtype,
                        i::SingleIndex,
                        col_inds::MultiIndex)
    dm.na[i, col_inds] = true
    return dm
end

# dm[SingleItemIndex, MultiItemIndex] = Multiple Items
function Base.setindex!{S, T}(dm::DataMatrix{S},
                              vals::Vector{T},
                              i::SingleIndex,
                              col_inds::MultiIndex)
    dm.data[i, col_inds] = vals
    dm.na[i, col_inds] = false
    return dm
end

# dm[SingleItemIndex, MultiItemIndex] = Single Item
function Base.setindex!(dm::DataMatrix,
                        val::Any,
                        i::SingleIndex,
                        col_inds::MultiIndex)
    dm.data[i, col_inds] = val
    dm.na[i, col_inds] = false
    return dm
end

# dm[MultiItemIndex, MultiItemIndex] = NA
function Base.setindex!(dm::DataMatrix,
                        val::NAtype,
                        row_inds::MultiIndex,
                        col_inds::MultiIndex)
    dm.na[row_inds, col_inds] = true
    return dm
end

# dm[MultiIndex, MultiIndex] = Multiple Items
function Base.setindex!{S, T}(dm::DataMatrix{S},
                              vals::Vector{T},
                              row_inds::MultiIndex,
                              col_inds::MultiIndex)
    dm.data[row_inds, col_inds] = vals
    dm.na[row_inds, col_inds] = false
    return dm
end

# dm[MultiItemIndex, MultiItemIndex] = Single Item
function Base.setindex!(dm::DataMatrix,
                        val::Any,
                        row_inds::MultiIndex,
                        col_inds::MultiIndex)
    dm.data[row_inds, col_inds] = val
    dm.na[row_inds, col_inds] = false
    return dm
end

# Extract the matrix diagonal
Base.diag{T}(dm::DataMatrix{T}) = DataArray(diag(dm.data), diag(dm.na))


text:
using ArcadeLearningEnvironment
using CartesianGeneticProgramming
using IICGP
using Random
using Test

const cfg_filename = string(@__DIR__, "/dualcgp_test.yaml")
const rom = "assault"
const roms = ["assault", "atlantis", "centipede", "tennis", "pong"]
const all_roms = setdiff(getROMList(), ["pacman", "surround"])
const seed = 0
Random.seed!(seed)

const max_frames = 10
const stickiness = 0.25 # main_cfg["stickiness"]
const grayscale = true # main_cfg["grayscale"]
const downscale = true # main_cfg["downscale"]
const lck = ReentrantLock()

global STTS = [] # Vector{Vector{Matrix{UInt8}}}(undef, max_frames)
global RGBS = []
global GRSS = []
global ACTS = [] # Vector{Int32}(undef, max_frames)
global REWS = [] # Vector{Int64}(undef, max_frames)
global STKS = [] # Vector{Float64}(undef, max_frames)
global RAMS = [] # Vector{Float64}(undef, max_frames)

function play_atari(
    encoder::CGPInd,
    reducer::Reducer,
    controller::CGPInd,
    lck::ReentrantLock,
    is_ref::Bool,
    rom::String,
    rom_state_ref::Ptr{Nothing};
    seed=seed,
    max_frames=max_frames,
    grayscale=grayscale,
    downscale=downscale,
    stickiness=stickiness,
    verbose=false
)
    # Random.seed!(seed)
    mt = MersenneTwister(seed)
    g = Game(rom, seed, lck=lck, state_ref=rom_state_ref)
    IICGP.reset!(reducer) # clean buffer
    # reward = 0.0
    frames = 0
    prev_action = 0
    # Buffers
    stts = []
    rgbs = []
    grss = []
    acts = []
    rews = []
    stks = []
    rams = []
    while ~game_over(g.ale)
        s = get_state(g, grayscale, downscale)
        rgb = get_rgb(g)
        gs = get_grayscale(g)
        ram = get_ram(g)
        stickydraw = rand(mt)
        if stickydraw > stickiness || frames == 0
            output = IICGP.process(encoder, reducer, controller, s)
            action = g.actions[argmax(output)]
            # action = g.actions[rand(mt, 1:length(g.actions))]
        else
            action = prev_action
        end
        frames += 1
        # action = a_seq[frames]
        r = act(g.ale, action)
        # reward += r
        if frames > max_frames
            break
        end
        # Fill buffers
        push!(stts, s)
        push!(rgbs, rgb)
        push!(grss, gs)
        push!(acts, action)
        push!(rews, r)
        push!(stks, stickydraw)
        push!(rams, ram)
    end
    close!(g)
    # Test
    global STTS
    global RGBS
    global GRSS
    global ACTS
    global REWS
    global STKS
    global RAMS
    if is_ref
        STTS = stts
        RGBS = rgbs
        GRSS = grss
        ACTS = acts
        REWS = rews
        STKS = stks
        RAMS = rams
    else
        @test stts == STTS
        @test rgbs == RGBS
        @test grss == GRSS
        @test acts == ACTS
        @test rews == REWS
        @test stks == STKS
        @test rams == RAMS
    end
end

@testset "Determinism on all Atari games for dualCGP" begin
    for r in all_roms
        main_cfg, enc_cfg, con_cfg, red, _ = IICGP.dualcgp_config(cfg_filename, r)
        state_ref = get_state_ref(r, seed)

        # Reference
        enc = IPCGPInd(enc_cfg)
        con = CGPInd(con_cfg)
        play_atari(enc, red, con, lck, true, r, state_ref)

        # Test if same individuals produce the same trajectories
        for _ in 1:2
            enc = IPCGPInd(enc_cfg, enc.chromosome)
            con = CGPInd(con_cfg, con.chromosome)
            play_atari(enc, red, con, lck, false, r, state_ref)
        end
    end
end

function act_for(e::Ptr{Nothing}, as::Vector{Int64}; verbose::Bool=false)
    gs = []
    w = getScreenWidth(e)
    h = getScreenHeight(e)
    for i in 1:length(as)
        act(e, as[i])
        gs_i = reshape(getScreenGrayscale(e), (w,h))
        push!(gs, gs_i)
        if game_over(e)
            if verbose
                println("\nReached frame number : ", i)
            end
            break
        end
    end
    gs
end

@testset "Determinism on ALE only" begin
    n_a = 50000
    for romfile in roms
        ale = ALE_new()
        loadROM(ale, romfile)
        actions = getMinimalActionSet(ale)
        w = getScreenWidth(ale)
        h = getScreenHeight(ale)
        ast = cloneSystemState(ale)
        a_seq = actions[rand(1:length(actions), n_a)]
        ags = act_for(ale, a_seq)
        # ags = reshape(getScreenGrayscale(ale), (w,h))
        # display(IICGP.implot(ags))

        for _ in 1:2
            ble = ALE_new()
            loadROM(ble, romfile)
            restoreSystemState(ble, ast)
            bst = cloneSystemState(ble)
            bgs = act_for(ble, a_seq)
            # bgs = reshape(getScreenGrayscale(ble), (w,h))
            # display(IICGP.implot(bgs))
            # println()
            # println("Same gs : ", bgs == ags)
            @test bgs == ags
            ALE_del(ble)
        end
        ALE_del(ale)
    end
end


text:
using SurrogateCube
# Define some Parameters
Ncomp = 3 #We want 3 independent components
Nvar  = 8 # we want 10 measured variables
noise = RedNoise(0.2,1.0,1.0,1.0) #Each component has a noise with sigma=0.5
means = [ConstantBaseline(0.0),ConstantBaseline(0.0),ConstantBaseline(0.0)] # All Baselines are flat
dists = [CubeEvent(0.2),EmptyEvent(),EmptyEvent()] # Only one of the components is disturbed with a size of 20% in each dimension
varNoise = WhiteNoise(0.1) # each variable generated has some independent white noise
k=10.0  # How strongly do we want to perturb
#Set dimensions
Nlon = 10
Nlat = 10
Ntime =100
ds = genDataCube(means,noise,dists,varNoise,Ncomp,Nvar,Ntime,Nlat,Nlon,k);
#Save the datacube to a file
save_cube(ds)
#Plot something
using PyPlot
pygui(true)
f=figure()
#plot(int(d[1,5,5,:][:]))
for i=1:Nvar
plot(ds.variables[:,5,5,i][:])
end

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
using DataFrames

export
    LogisticRegressionModel,
    fit!,
    predict

mutable struct LogisticRegressionModel{T, S}
    y::Vector{T}
    xs::AbstractArray{S}
    argv::Vector{Float64}
    n::Integer
end

function LogisticRegressionModel(df::DataFrame, label::Symbol, features::Vector{Symbol})
    n = nrow(df)
    y = df[!, label]
    xs = Matrix(df[!, features])
    β = zeros(length(features))
    return LogisticRegressionModel(y, xs, β, n)
end

function LogisticRegressionModel(df::DataFrame, label::Symbol, feature::Symbol)
    n = nrow(df)
    y = df[!, label]
    x = df[!, feature]
    β = zeros(1)
    return LogisticRegressionModel(y, x, β, n)
end

sigmoid(z::Real) = 1.0 / (1.0 + exp(-z))

z(model::LogisticRegressionModel, xs::Vector{<:Real}) = model.argv' * xs

predict(model::LogisticRegressionModel, xs::Vector{<:Real}) = sigmoid(z(model, xs))

predict(model::LogisticRegressionModel, xs::Real) = sigmoid(z(model, [xs]))

function log_likelyhood(model::LogisticRegressionModel)
    xs = model.xs
    y = model.y

    ll = 0
    for (i, yⁱ) in enumerate(y)
        xsⁱ = xs[i, :]
        ll += yⁱ * log(sigmoid(z(model, xsⁱ)))
        ll += (1 - yⁱ) * log(1 - sigmoid(z(model, xsⁱ)))
    end

    return ll ./ model.n
end

function ∇log_likelyhood(model::LogisticRegressionModel)
    # ∇ll = ∑(yⁱ - σ(zⁱ))xⱼⁱ
    xs = model.xs
    y = model.y

    ∇ll = zeros(length(model.argv))
    for (i, yⁱ) in enumerate(y)
        xsⁱ = xs[i, :]
        ∇ll .+= (yⁱ - sigmoid(z(model, xsⁱ))) .* xsⁱ
    end

    return ∇ll ./ model.n
end

function gradient_ascend(
    model::LogisticRegressionModel,
    η::Real,
    atol::Real,
    show::Bool
)
    β = model.argv
    while (ll = log_likelyhood(model)) < atol
        show && println("Log Likelyhood: $ll")
        β .+= η .* ∇log_likelyhood(model)
    end

    return β
end

function fit!(
    model::LogisticRegressionModel;
    method=gradient_ascend,
    η::Real=1e-2,
    atol::Real=-1e-4,
    show=false
)
    β = method(model, η, atol, show)
    model.argv .= β

    return model
end


text:
struct MyDijkstraState{T<:Integer,W<:Real} <: Graphs.AbstractPathState
    parents::Vector{T}
    dists::Vector{W}
end

function my_dijkstra_shortest_paths(
    g::AbstractGraph{T},
    s::Integer;
    edge_indices::Dict,
    edge_weights::AbstractVector{W},
) where {T<:Integer, W<:Real}
    dists = fill(typemax(W), nv(g))
    parents = zeros(T, nv(g))
    Q = PriorityQueue{T,W}()

    dists[s] = zero(W)
    enqueue!(Q, s, zero(W))

    while !isempty(Q)
        u = dequeue!(Q)
        d_u = dists[u]
        for v in outneighbors(g, u)
            w_uv = edge_weights[edge_indices[u, v]]
            dist_through_u = d_u + w_uv
            if dist_through_u < dists[v]
                dists[v] = dist_through_u
                parents[v] = u
                enqueue!(Q, v, dist_through_u)
            end
        end
    end

    return MyDijkstraState{T,W}(parents, dists)
end


text:
function single_ion_adsorption_water_flood_single_shock(core_props, fluids_ls, fluids_hs, rel_perms_hs,
    rel_perms_ls, core_flood, eq_const)
    # ====================================================================
    # This function is not done yet. The trick is to draw a tangent from
    # the high sal 1-sor_hs to the low sal curve to find the shock front
    # saturation. The low sal shock speed is calculated at this saturation
    # by fw/(sw+eq_const)
    # The implementation needs some if statements to make sure that the
    # shock front saturation calculated by the tangent line from -eq_const
    # to the low sal fw is indeed lower than 1-sor_hs
    # ====================================================================
    # construct the fractional flow curves
    fw_ls, dfw_ls = fractional_flow_function(rel_perms_ls, fluids_ls)
    fw_hs, dfw_hs = fractional_flow_function(rel_perms_hs, fluids_hs)
    # low sal shock (tangent line from (0,0))
    sw_shock_ls = tangent_line_saturation(rel_perms_ls, fluids_ls, (-eq_const, 0.0))
    # println("low sal sw_shock = $sw_shock_ls")
    t_D_BT_ls = 1/dfw_ls(sw_shock_ls) # breakthrough (BT) time [#PV]
    # println("low sal breakthrough time = $t_D_BT_ls")
    # High sal shock (cross point between the ls tangent and the hs fw)
    sw_shock_hs = cross_point_saturation(fw_hs, rel_perms_hs, (-eq_const, 0.0),
            (sw_shock_ls, fw_ls(sw_shock_ls)))
    # println("high sal sw_shock = $sw_shock_hs")
    sw_init = core_flood.initial_water_saturation
    t_D_BT_hs = (sw_shock_hs-sw_init)/(fw_hs(sw_shock_hs)-fw_hs(sw_init)) # breakthrough (BT) time [#PV]
    # println("high sal breakthrough time = $t_D_BT_hs")
    if t_D_BT_hs<t_D_BT_ls
        # @info "Using the single_ion_adsorption_water_flood function."
        return single_ion_adsorption_water_flood(core_props, fluids_ls, fluids_hs, rel_perms_hs,
                rel_perms_ls, core_flood, eq_const)
    end
    sor_ls = rel_perms_ls.sor
    sw_shock = cross_point_saturation(fw_ls, rel_perms_ls, (sw_shock_hs, fw_hs(sw_shock_hs)),
        (sw_init, fw_hs(sw_init)), sw_left = sw_shock_ls, sw_right = 1-sor_ls)
    # sw_shock = cross_point_saturation(fw_ls, rel_perms_ls, (-eq_const, 0.0),
    #         (sw_init, fw_hs(sw_init)), sw_left = sw_shock_ls, sw_right = 1-sor_ls) # just a test
    # sw_shock = cross_point_saturation(fw_ls, rel_perms_ls, (sw_shock_ls, fw_ls(sw_shock_ls)),
    #     (sw_init, fw_hs(sw_init)), sw_left = sw_shock_ls, sw_right = 1-sor_ls)
    if sw_shock<sw_shock_ls
        sw_shock = cross_point_saturation(fw_ls, rel_perms_ls, (sw_shock_hs, fw_hs(sw_shock_hs)),
            (sw_init, fw_hs(sw_init)), sw_left = sw_shock_ls, sw_right = 1-sor_ls)
    end
    sor_hs = rel_perms_hs.sor
    # t_D_BT_hs = (sw_shock-sw_init)/(fw_ls(sw_shock)-fw_hs(sw_init))
    t_D_BT_hs = (sw_shock+eq_const)/(fw_ls(sw_shock)-0.0)
    # println("sw_shock = $sw_shock")
    t_D_BT_ls = 1/dfw_ls(sw_shock)
    # println("low sal breakthrough time = $t_D_BT_ls")

    # tracer
    # I'm not sure about this, but the numerical solution suggests that
    # the tracer front follows the shock front for this special case
    # sw_tracer = tangent_line_saturation(rel_perms_ls, fluids_ls, (0.0, 0.0))
    # sw_tracer = max(sw_tracer, sw_shock)
    t_D_tracer = t_D_BT_hs # (sw_shock-sw_init)/(fw_ls(sw_shock)-fw_hs(sw_init))
    xt_tracer_shock = 1/t_D_tracer

    # construct the recovery factor curve versus the # of PV
    # R = zeros(1)
    # pv_R = zeros(1)
    # at breakthrough of the hs brine
    sw_init = core_flood.initial_water_saturation
    # push!(R, 0.0) # recovery at BT
    # push!(pv_R, t_D_BT_ls) # BT time
    # at breakthrough of the ls brine
    # push!(R, (1-fw_hs(sw_shock))*(t_D_BT_ls-t_D_BT_hs)/(1-sw_init)) # recovery at BT
    # push!(pv_R, t_D_BT_ls) # BT time

    # after breakthrough
    pv_inj = max(core_flood.injected_pore_volume, 2.0) # at least inject 2 pv
    # f_sw = sw -> (pv_inj-1/dfw_ls(sw)) # find the outlet saturation at the end of injection
    # sw_max = fzero(f_sw, sw_shock)
    # println("swmax $sw_max")
    # sw_tmp = linspace(sw_shock, sw_max, 100)
    # t_D_tmp = 1./dfw_ls.(sw_tmp)

    # s_av_tmp = sw_tmp-(fw_ls.(sw_tmp)-1).*t_D_tmp
    # R_tmp = (s_av_tmp-sw_init)/(1-sw_init)
    # append!(R, R_tmp)
    # append!(pv_R, t_D_tmp)

    # saturation profile
    sw_inj = core_flood.injected_water_saturation
    phi = core_props.porosity
    ut = core_flood.injection_velocity
    L = core_props.length
    pv_to_t = phi*L/ut
    s1 = collect(linspace(min(sw_inj, 1-sor_ls-eps()), sw_shock, 300))
    xt_s1 = dfw_ls.(s1)
    xt_shock_ls = 1/t_D_BT_ls
    xt_shock_hs = 1/t_D_BT_hs
    xt_prf=[xt_s1; xt_shock_ls+eps(); xt_shock_hs; xt_shock_hs+eps(); 1/0.3]
    sw_prf=[s1; sw_shock; sw_shock; sw_init; sw_init]
    xt_tracer = [0.0, xt_tracer_shock, xt_tracer_shock+eps(), xt_prf[end]]
    c_tracer = [1.0, 1.0, 0.0, 0.0]

    # extract data for calculating the pressure drop
    k = core_props.permeability
    krw_ls, kro_ls, dkrw_ls, dkro_ls = rel_perm_functions(rel_perms_ls)
    muo_ls, muw_ls = fluids_ls.oil_viscosity, fluids_ls.water_viscosity
    krw_hs, kro_hs, dkrw_hs, dkro_hs = rel_perm_functions(rel_perms_hs)
    muo_hs, muw_hs = fluids_hs.oil_viscosity, fluids_hs.water_viscosity
    λ_ls = s -> krw_ls(s)/muw_ls+kro_ls(s)/muo_ls
    λ_hs = s -> krw_hs(s)/muw_hs+kro_hs(s)/muo_hs
    # λ_t = [λ_ls.(s1); λ_hs(sw_shock); λ_hs(sw_shock); λ_hs(sw_init); λ_hs(sw_init)]
    λ_t = [λ_ls.(s1); λ_ls(sw_shock); λ_ls(sw_shock); λ_hs(sw_init); λ_hs(sw_init)]

    λ_w_ls = s -> krw_ls(s)/muw_ls
    λ_w_hs = s -> krw_hs(s)/muw_hs
    # λ_w = [λ_w_ls.(s1); λ_w_hs(sw_shock); λ_w_hs(sw_shock); λ_w_hs(sw_init); λ_w_hs(sw_init)]
    λ_w = [λ_w_ls.(s1); λ_w_ls(sw_shock); λ_w_ls(sw_shock); λ_w_hs(sw_init); λ_w_hs(sw_init)]

    xt = copy(xt_prf)
    sw = copy(sw_prf)
    i=1
    while(true)
        if i==length(xt)
            break
        elseif xt[i]>=xt[i+1]
            deleteat!(xt, i)
            deleteat!(sw, i)
            deleteat!(λ_t, i)
            deleteat!(λ_w, i)
        else
            i+=1
        end
    end

    x = collect(linspace(0,L,200))
    sw_int = Spline1D(xt, sw, k=1, bc="nearest")
    λ_int = Spline1D(xt, λ_t, k=1, bc="nearest")
    λ_w_int = Spline1D(xt, λ_w, k=1, bc="nearest")
    t_inj=pv_inj*pv_to_t
    t = collect([linspace(0,t_D_BT_hs*pv_to_t, 10);
                linspace((t_D_BT_hs+eps())*pv_to_t, t_D_BT_ls*pv_to_t, 40);
                linspace((t_D_BT_ls+eps())*pv_to_t,t_inj, 100)]) # [s] time
    p_inj = zeros(length(t))
    R_int = zeros(length(t))
    WC_int = zeros(length(t))
    p_inj[1]=L*ut./(k*λ_hs(sw_init))
    WC_int[1]=λ_w_hs(sw_init)/λ_hs(sw_init)
    for i in 2:length(t)
        xt_real = x/t[i]/(ut/phi)
        p_inj[i] = trapz(x, ut./(k*λ_int(xt_real)))
        R_int[i] = (trapz(x, sw_int(xt_real))/L-sw_init)/(1-sw_init)
        WC_int[i] = λ_w_int(xt_real[end])/λ_int(xt_real[end])
    end

    R_int[R_int.<0.0] .= 0.0

    return FracFlowResults([fw_hs, fw_ls], [Line([-eq_const, 0.0], [sw_shock_ls, fw_ls(sw_shock_ls)]),
                            Line([sw_init, fw_hs(sw_init)], [sw_shock_hs, fw_hs(sw_shock_hs)]),
                            Line([sw_shock, fw_ls(sw_shock)], [sw_init, fw_hs(sw_init)])],
                            [t/pv_to_t R_int], [t R_int], [xt sw], [xt_tracer c_tracer],
                            [t/pv_to_t p_inj], [t p_inj], [t/pv_to_t WC_int], [t WC_int])
    # return pv_R, R, xt_prf, sw_prf # for the time being to test the code
end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
"""
    pyrange([[start], [stop]], [step])

Construct a Python `range`. Unspecified arguments default to `None`.
"""
pyrange(x, y, z) = pybuiltins.range(x, y, z)
pyrange(x, y) = pybuiltins.range(x, y)
pyrange(y) = pybuiltins.range(y)
export pyrange

pyrange_fromrange(x::AbstractRange) = pyrange(first(x), last(x) + sign(step(x)), step(x))

pyisrange(x) = pytypecheck(x, pybuiltins.range)

function pyconvert_rule_range(::Type{R}, x::Py, ::Type{StepRange{T0,S0}}=Utils._type_lb(R), ::Type{StepRange{T1,S1}}=Utils._type_ub(R)) where {R<:StepRange,T0,S0,T1,S1}
    a = @pyconvert(Utils._typeintersect(Integer, T1), x.start)
    b = @pyconvert(Utils._typeintersect(Integer, S1), x.step)
    c = @pyconvert(Utils._typeintersect(Integer, T1), x.stop)
    a′, c′ = promote(a, c - oftype(c, sign(b)))
    T2 = Utils._promote_type_bounded(T0, typeof(a′), typeof(c′), T1)
    S2 = Utils._promote_type_bounded(S0, typeof(c′), S1)
    pyconvert_return(StepRange{T2, S2}(a′, b, c′))
end

function pyconvert_rule_range(::Type{R}, x::Py, ::Type{UnitRange{T0}}=Utils._type_lb(R), ::Type{UnitRange{T1}}=Utils._type_ub(R)) where {R<:UnitRange,T0,T1}
    b = @pyconvert(Utils._typeintersect(Integer, S1), x.step)
    b == 1 && return pyconvert_unconverted()
    a = @pyconvert(Utils._typeintersect(Integer, T1), x.start)
    c = @pyconvert(Utils._typeintersect(Integer, T1), x.stop)
    a′, c′ = promote(a, c - oftype(c, 1))
    T2 = Utils._promote_type_bounded(T0, typeof(a′), typeof(c′), T1)
    pyconvert_return(UnitRange{T2}(a′, c′))
end

pyconvert_rule_fast(::Type{T}, x::Py) where {T<:StepRange} = pyisrange(x) ? pyconvert_rule_range(T, x) : pyconvert_unconverted()
pyconvert_rule_fast(::Type{T}, x::Py) where {T<:UnitRange} = pyisrange(x) ? pyconvert_rule_range(T, x) : pyconvert_unconverted()


text:
@enum PINS begin
   EMPTY
   PIN
   INVALID
end

start_pos = [
   INVALID,INVALID,PIN,PIN,  PIN,INVALID,INVALID,
   INVALID,INVALID,PIN,PIN,  PIN,INVALID,INVALID,
   PIN,    PIN,    PIN,PIN,  PIN,PIN,    PIN,
   PIN,    PIN,    PIN,EMPTY,PIN,PIN,    PIN,
   PIN,    PIN,    PIN,PIN,  PIN,PIN,    PIN,
   INVALID,INVALID,PIN,PIN,  PIN,INVALID,INVALID,
   INVALID,INVALID,PIN,PIN,  PIN,INVALID,INVALID,
]

# inversion of starting position
solution_pos = map(x -> if x == PIN EMPTY elseif x == EMPTY PIN else INVALID end, start_pos)

mutable struct position_t
    moves::Vector{Tuple{Int64,Int64}}
    pos::Vector{PINS}
end

//u d l r

function valid_left(pos::Vector{PINS}, pin::Int64)
    endv = pin - 2;
    midv = pin - 1;
    return  ((endv > 1) && (((pin - 1) ÷ 7) == ((endv - 1) ÷ 7)) && (pos[pin] == PIN) && (pos[midv] == PIN) && (pos[endv] == EMPTY))
end

function valid_right(pos::Vector{PINS}, pin::Int64)
    endv = pin + 2;
    midv = pin + 1;
    return  ((endv < 50) && (((pin - 1) ÷ 7) == ((endv - 1) ÷ 7)) && (pos[pin] == PIN) && (pos[midv] == PIN) && (pos[endv] == EMPTY))
end

function valid_up(pos::Vector{PINS}, pin::Int64)
    endv = pin - 14
    midv = pin - 7
    return  ((endv >= 1) && (pos[pin] == PIN) && (pos[midv] == PIN) && (pos[endv] == EMPTY))
end

function valid_dn(pos::Vector{PINS}, pin::Int64)
    endv = pin + 14;
    midv = pin + 7;
    return  ((endv < 50) && (pos[pin] == PIN) && (pos[midv] == PIN) && (pos[endv] == EMPTY))
end

//fast d u l r
//slow l r d u
function get_moves(pos::Vector{PINS})
    moves = []

    for i = 1:49
        if valid_dn(pos, i)
            push!(moves, (i, i + 14))
        end
        if valid_up(pos, i)
            push!(moves, (i, i - 14))
        end
        if valid_left(pos, i)
            push!(moves, (i, i - 2))
        end
        if valid_right(pos, i)
            push!(moves, (i, i + 2))
        end
    end

    return moves
end

function solution_p(a::Vector{PINS}, b::Vector{PINS})
    return (a == b)
end

function make_move(pos::Vector{PINS}, move::Tuple{Int64,Int64})
    new_pos = copy(pos)

    new_pos[move[1]] = EMPTY;
    new_pos[move[2]] = PIN;
    new_pos[(move[1] + move[2]) ÷ 2] = EMPTY;

    return new_pos
end

function print_pos(pos)
    for i = 1:49
        if pos[i] == INVALID
            print(" ")
        elseif pos[i] == PIN
         print("X")
      elseif pos[i] == EMPTY
        print("O")
        end
        if (i % 7) == 0
            println("")
        end
    end
end

function find_solution(spos)
    solution = false
    sol_path = []
    push!(sol_path, spos)

    while !isempty(sol_path) && !solution
        spos = pop!(sol_path)
        solution = solution_p(spos.pos, solution_pos)
        if !solution
            moves = get_moves(spos.pos)
            for move in moves
                npos = position_t([], [])
                npos.pos = make_move(spos.pos, move)
                npos.moves = copy(spos.moves)
                push!(npos.moves, move)
                push!(sol_path, npos)
            end
        else
            println("Length = ", length(spos.moves))
        end
    end

    return (solution, spos)
end

function main()
    start = position_t([], start_pos)
    found, sol = find_solution(start)
    if found
        println("Found solution")
        print_pos(sol.pos)
      # Print solution path
        for move in sol.moves
            println("From ", move[1], " To ", move[2])
        end
    else
        println("There is no solution !")
    end
end

#get_moves(start_pos)

#print_pos(make_move(start_pos, (11, 25)))
#print_pos(make_move(start_pos, (23, 25)))
#print_pos(make_move(start_pos, (27, 25)))
#print_pos(make_move(start_pos, (39, 25)))
#solution_p(start_pos, solution_pos)
#solution_p(solution_pos, solution_pos)

main()


text:
function DynamicOutlierFilter(arg0::jint, arg1::jdouble)
    return DynamicOutlierFilter((jint, jdouble), arg0, arg1)
end

function get_parameters_drivers(obj::OutlierFilter)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_sigma(obj::DynamicOutlierFilter)
    return jcall(obj, "getSigma", Vector{jdouble}, ())
end

function modify(obj::DynamicOutlierFilter, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end

function set_sigma(obj::DynamicOutlierFilter, arg0::Vector{jdouble})
    return jcall(obj, "setSigma", void, (Vector{jdouble},), arg0)
end


text:
## Part 1
function _load_input(name::AbstractString)
    path = joinpath("inputs", name)
    return map(x -> parse(Int64, x), readlines(path))
end

function num_depth_increases(depths::AbstractVector{<:Real})
    length(depths) > 1 || throw(ArgumentError("Need at least two depths"))
    count = 0
    previous = @inbounds first(depths)
    for x in @view(depths[2:end])
        if x > previous
            count += 1
        end
        previous = x
    end
    return count
end

example = [199, 200, 208, 210, 200, 207, 240, 269, 260, 263]
@assert num_depth_increases(example) == 7
depths = _load_input("day_01.txt")
num_depth_increases(depths)

## Part 2
using Dates
using TimeDag

function num_depth_increases(depths::AbstractVector{<:Real}, window::Int64)
    t_start = DateTime(2000)
    block = Block(
        TimeDag.unchecked, t_start:Day(1):(t_start + Day(1) * (length(depths) - 1)), depths
    )

    _eval(x) = evaluate(x, t_start, last(block.times) + Day(1))

    # Form a node whose last value will be the answer.
    depths_node = block_node(block)
    window_sum = sum(depths_node, window)
    previous = lag(window_sum, 1)
    # There's a bug in TimeDag where you can't do sums of Bools correctly right now:
    #   https://github.com/invenia/TimeDag.jl/issues/22
    result_node = sum(TimeDag.apply((x, y) -> x > y ? 1 : 0, window_sum, previous))
    return last(_eval(result_node).values)
end

@assert num_depth_increases(example, 3) == 5
@time num_depth_increases(depths, 3)


text:
using SnpArrays
x = SnpArray("test")
snpmatrix = convert(Array{Float64,2}, x)

for i in 1:size(snpmatrix, 2)
	snpmatrix[:, i] = (snpmatrix[:, i] .- mean(snpmatrix[:, i])) / std(snpmatrix[:, i])
end

using SnpArrays
srand(2018)
x = SnpArray(rand(0:2, 5, 3))
view(x, :, [true, false, true])
#how to do the following?
view(x, :, [true, false, true]) * rand(2)

using IHT
# x = MendelIHT("test_control.txt")
x = MendelIHT("gwas 1 Control.txt")

using MendelIterHardThreshold
# x = IterHardThreshold("test_control.txt")
x = IterHardThreshold("gwas 1 Control.txt")

function plus()
	s = 0
	@simd for i = 1:10^6
		@inbounds s += i
	end
	return s
end

function plus_new()
	s = 0
	for i = 1:10^6
		@inbounds s += i
	end
	return s
end

#benchmark std computation
using SnpArrays, BenchmarkTools
x = SnpData("gwas 1 data")
std_vec = zeros(x.snps)
function hi()
    for i in 1:x.snps
        storage = view(x.snpmatrix, :, i)
        std_vec[i] .= 1.0 / std( convert(Vector{Float64}, storage) )
    end
	return std_vec
end

storage = zeros(x.people)
function hii()
    for i in 1:x.snps
    	copy!(storage, x.snpmatrix[:, i]) #in place version of convert
        std_vec[i] .= 1.0 ./ std(storage)
    end
	return std_vec
end

@benchmark hi()
@benchmark hii()

function hiii()
    mean_vec = 2.0x.maf #multiply by 2 because mean of each snp = 2.0*maf
    std_vec = 1.0 ./ sqrt.(2.0 .* x.maf .* (1.0 .- x.maf))
    return mean_vec, std_vec
end

using SnpArrays, BenchmarkTools
x = SnpArray("gwas 1 data")
z = bitrand(10000)
y = SnpArray(2200, sum(z))
function hii()
	x = SnpArray("gwas 1 data")
	z = bitrand(10000)
	y = SnpArray(2200, sum(z))
	y .= view(x, :, z)
end

using SnpArrays, PLINK
x = SnpData("test")
y = BEDFile("test.bed")
x.snpmatrix
y[:, 1:2] #obmit 3rd column, which is the intercept
convert(Matrix{Float64}, x.snpmatrix)


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, Pkg

# See https://github.com/JuliaLang/Pkg.jl/issues/2942
# Once this Pkg issue is resolved, this must be removed
uuid = Base.UUID("a83860b7-747b-57cf-bf1f-3e79990d037f")
delete!(Pkg.Types.get_last_stdlibs(v"1.6.3"), uuid)

name = "z3"
version = v"4.8.14"
julia_versions = [v"1.6.3", v"1.7.0", v"1.8.0", v"1.9.0"]

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Z3Prover/z3.git", "df8f9d7dcb8b9f9b3de1072017b7c2b7f63f0af8"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.15.sdk.tar.xz",
                  "2408d07df7f324d3beea818585a6d990ba99587c218a3969f924dfcc4de93b62"),
]

macfix = raw"""
# See https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/1185
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    # work around macOS SDK issue
    #     /workspace/srcdir/z3/src/ast/ast.h:: 189error:: 47:'get<unsigned int, int, ast *,
    #         symbol, zstring *, rational *, double, unsigned int>' is unavailable:
    #         introduced in macOS 10.14
    export MACOSX_DEPLOYMENT_TARGET=10.15
    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd
fi
"""

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

mkdir z3-build && cd z3-build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_FIND_ROOT_PATH="${prefix}" \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DZ3_BUILD_JULIA_BINDINGS=True \
    -DJulia_PREFIX="${prefix}" \
    ..
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/z3/LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = vcat(libjulia_platforms.(julia_versions)...)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libz3", :libz3),
    LibraryProduct("libz3jl", :libz3jl),
    ExecutableProduct("z3", :z3)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("libjulia_jll"),
    Dependency("libcxxwrap_julia_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms,
               products, dependencies; preferred_gcc_version=v"9",
               julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using InvertibleNetworks, LinearAlgebra, Test, Random

# Data
nx = 16
ny = 16
n_in = 3
batchsize = 4
X = randn(Float32, nx, ny, n_in, batchsize)

# Network
nscales = 2
steps_per_scale = 4
nc = 2
hf = 32
α0 = .2f0
H = NetworkHyperbolic(nx, ny, n_in, batchsize, nscales, steps_per_scale; α=α0, hidden_factor=hf, ncenter=nc)

###################################################################################################
# Test invertibility

X_ = H.inverse(H.forward(X)[1])
@show norm(X - X_), norm(X)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1e-2)

X_ = H.forward(H.inverse(X))[1]
@show norm(X - X_), norm(X)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1e-2)

####################################################################################################
# Training

# Loss
function loss(H, X, ind)
    Y, logdet = H.forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    ΔX, X = H.backward(ΔY, Y)
    return f, ΔX, H.HL[ind].W.grad, H.AL.s.grad
end

# Data
X = randn(Float32, nx, ny, n_in, batchsize)
X0 = randn(Float32, nx, ny, n_in, batchsize)
dX = X - X0

# Gradient test w.r.t. input X0
Y = H.forward(X)
f0, ΔX = loss(H, X0, 1)[1:2]
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test hyperbolic network\n")
for j=1:maxiter
    f = loss(H, X0 + h*dX, 1)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

@show rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
@show rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)

# Gradient test w.r.t. weights of hyperbolic network
H = NetworkHyperbolic(nx, ny, n_in, batchsize, nscales, steps_per_scale; α=α0, hidden_factor=hf, ncenter=nc)
H0 = NetworkHyperbolic(nx, ny, n_in, batchsize, nscales, steps_per_scale; α=α0, hidden_factor=hf, ncenter=nc)
H.forward(X)
H0.forward(X)   # evaluate to initialize actnorm layer

s0 = H0.AL.s.data
ds = H.AL.s.data - H0.AL.s.data

maxiter = 6
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

for i=1:length(H0.HL)
    h1 = 0.1f0
    print("\nGradient test invertible layer for layer $(i)\n")
    W0s = H0.HL[i].W.data
    dWg = H.HL[i].W.data - W0s
    f01, ΔX1, ΔWg, _ = loss(H0, X, i)
    for j=1:maxiter
        H0.HL[i].W.data = W0s + h1*dWg
        f = loss(H0, X, i)[1]
        err3[j] = abs(f - f01)
        err4[j] = abs(f - f01 - h1*dot(dWg, ΔWg))
        print(err3[j], "; ", err4[j], "\n")
        h1 = h1/2f0
    end

    @show local rate_1 = sum(err3[1:end-1]./err3[2:end])/(maxiter - 1)
    @show local rate_2 = sum(err4[1:end-1]./err4[2:end])/(maxiter - 1)
    H0.HL[i].W.data = W0s
    @test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=2f1)
    @test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=2f1)
end

print("\nGradient test invertible layer for affine layer\n")
f0, ΔX, ΔW, Δs = loss(H0, X, 1)
for j=1:maxiter
    H0.AL.s.data = s0 + h*ds
    f = loss(H0, X, 1)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(ds, Δs))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@show rate_1 = sum(err3[1:end-1]./err3[2:end])/(maxiter - 1)
@show rate_2 = sum(err4[1:end-1]./err4[2:end])/(maxiter - 1)

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=2f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=2f1)


text:
# This is the union of every color defined in X11 and in SVG, prefering the SVG
# definition when they clash.

const color_names = Dict(
    "aliceblue"            => (240, 248, 255),
    "antiquewhite"         => (250, 235, 215),
    "antiquewhite1"        => (255, 239, 219),
    "antiquewhite2"        => (238, 223, 204),
    "antiquewhite3"        => (205, 192, 176),
    "antiquewhite4"        => (139, 131, 120),
    "aqua"                 => (  0, 255, 255),
    "aquamarine"           => (127, 255, 212),
    "aquamarine1"          => (127, 255, 212),
    "aquamarine2"          => (118, 238, 198),
    "aquamarine3"          => (102, 205, 170),
    "aquamarine4"          => ( 69, 139, 116),
    "azure"                => (240, 255, 255),
    "azure1"               => (240, 255, 255),
    "azure2"               => (224, 238, 238),
    "azure3"               => (193, 205, 205),
    "azure4"               => (131, 139, 139),
    "beige"                => (245, 245, 220),
    "bisque"               => (255, 228, 196),
    "bisque1"              => (255, 228, 196),
    "bisque2"              => (238, 213, 183),
    "bisque3"              => (205, 183, 158),
    "bisque4"              => (139, 125, 107),
    "black"                => (  0,   0,   0),
    "blanchedalmond"       => (255, 235, 205),
    "blue"                 => (  0,   0, 255),
    "blue1"                => (  0,   0, 255),
    "blue2"                => (  0,   0, 238),
    "blue3"                => (  0,   0, 205),
    "blue4"                => (  0,   0, 139),
    "blueviolet"           => (138,  43, 226),
    "brown"                => (165,  42,  42),
    "brown1"               => (255,  64,  64),
    "brown2"               => (238,  59,  59),
    "brown3"               => (205,  51,  51),
    "brown4"               => (139,  35,  35),
    "burlywood"            => (222, 184, 135),
    "burlywood1"           => (255, 211, 155),
    "burlywood2"           => (238, 197, 145),
    "burlywood3"           => (205, 170, 125),
    "burlywood4"           => (139, 115,  85),
    "cadetblue"            => ( 95, 158, 160),
    "cadetblue1"           => (152, 245, 255),
    "cadetblue2"           => (142, 229, 238),
    "cadetblue3"           => (122, 197, 205),
    "cadetblue4"           => ( 83, 134, 139),
    "chartreuse"           => (127, 255,   0),
    "chartreuse1"          => (127, 255,   0),
    "chartreuse2"          => (118, 238,   0),
    "chartreuse3"          => (102, 205,   0),
    "chartreuse4"          => ( 69, 139,   0),
    "chocolate"            => (210, 105,  30),
    "chocolate1"           => (255, 127,  36),
    "chocolate2"           => (238, 118,  33),
    "chocolate3"           => (205, 102,  29),
    "chocolate4"           => (139,  69,  19),
    "coral"                => (255, 127,  80),
    "coral1"               => (255, 114,  86),
    "coral2"               => (238, 106,  80),
    "coral3"               => (205,  91,  69),
    "coral4"               => (139,  62,  47),
    "cornflowerblue"       => (100, 149, 237),
    "cornsilk"             => (255, 248, 220),
    "cornsilk1"            => (255, 248, 220),
    "cornsilk2"            => (238, 232, 205),
    "cornsilk3"            => (205, 200, 177),
    "cornsilk4"            => (139, 136, 120),
    "crimson"              => (220,  20,  60),
    "cyan"                 => (  0, 255, 255),
    "cyan1"                => (  0, 255, 255),
    "cyan2"                => (  0, 238, 238),
    "cyan3"                => (  0, 205, 205),
    "cyan4"                => (  0, 139, 139),
    "darkblue"             => (  0,   0, 139),
    "darkcyan"             => (  0, 139, 139),
    "darkgoldenrod"        => (184, 134,  11),
    "darkgoldenrod1"       => (255, 185,  15),
    "darkgoldenrod2"       => (238, 173,  14),
    "darkgoldenrod3"       => (205, 149,  12),
    "darkgoldenrod4"       => (139, 101,   8),
    "darkgray"             => (169, 169, 169),
    "darkgreen"            => (  0, 100,   0),
    "darkgrey"             => (169, 169, 169),
    "darkkhaki"            => (189, 183, 107),
    "darkmagenta"          => (139,   0, 139),
    "darkolivegreen"       => ( 85, 107,  47),
    "darkolivegreen1"      => (202, 255, 112),
    "darkolivegreen2"      => (188, 238, 104),
    "darkolivegreen3"      => (162, 205,  90),
    "darkolivegreen4"      => (110, 139,  61),
    "darkorange"           => (255, 140,   0),
    "darkorange1"          => (255, 127,   0),
    "darkorange2"          => (238, 118,   0),
    "darkorange3"          => (205, 102,   0),
    "darkorange4"          => (139,  69,   0),
    "darkorchid"           => (153,  50, 204),
    "darkorchid1"          => (191,  62, 255),
    "darkorchid2"          => (178,  58, 238),
    "darkorchid3"          => (154,  50, 205),
    "darkorchid4"          => (104,  34, 139),
    "darkred"              => (139,   0,   0),
    "darksalmon"           => (233, 150, 122),
    "darkseagreen"         => (143, 188, 143),
    "darkseagreen1"        => (193, 255, 193),
    "darkseagreen2"        => (180, 238, 180),
    "darkseagreen3"        => (155, 205, 155),
    "darkseagreen4"        => (105, 139, 105),
    "darkslateblue"        => ( 72,  61, 139),
    "darkslategray"        => ( 47,  79,  79),
    "darkslategray1"       => (151, 255, 255),
    "darkslategray2"       => (141, 238, 238),
    "darkslategray3"       => (121, 205, 205),
    "darkslategray4"       => ( 82, 139, 139),
    "darkslategrey"        => ( 47,  79,  79),
    "darkturquoise"        => (  0, 206, 209),
    "darkviolet"           => (148,   0, 211),
    "deeppink"             => (255,  20, 147),
    "deeppink1"            => (255,  20, 147),
    "deeppink2"            => (238,  18, 137),
    "deeppink3"            => (205,  16, 118),
    "deeppink4"            => (139,  10,  80),
    "deepskyblue"          => (  0, 191, 255),
    "deepskyblue1"         => (  0, 191, 255),
    "deepskyblue2"         => (  0, 178, 238),
    "deepskyblue3"         => (  0, 154, 205),
    "deepskyblue4"         => (  0, 104, 139),
    "dimgray"              => (105, 105, 105),
    "dimgrey"              => (105, 105, 105),
    "dodgerblue"           => ( 30, 144, 255),
    "dodgerblue1"          => ( 30, 144, 255),
    "dodgerblue2"          => ( 28, 134, 238),
    "dodgerblue3"          => ( 24, 116, 205),
    "dodgerblue4"          => ( 16,  78, 139),
    "firebrick"            => (178,  34,  34),
    "firebrick1"           => (255,  48,  48),
    "firebrick2"           => (238,  44,  44),
    "firebrick3"           => (205,  38,  38),
    "firebrick4"           => (139,  26,  26),
    "floralwhite"          => (255, 250, 240),
    "fuchsia"              => (255,   0, 255),
    "forestgreen"          => ( 34, 139,  34),
    "gainsboro"            => (220, 220, 220),
    "ghostwhite"           => (248, 248, 255),
    "gold"                 => (255, 215,   0),
    "gold1"                => (255, 215,   0),
    "gold2"                => (238, 201,   0),
    "gold3"                => (205, 173,   0),
    "gold4"                => (139, 117,   0),
    "goldenrod"            => (218, 165,  32),
    "goldenrod1"           => (255, 193,  37),
    "goldenrod2"           => (238, 180,  34),
    "goldenrod3"           => (205, 155,  29),
    "goldenrod4"           => (139, 105,  20),
    "gray"                 => (128, 128, 128),
    "gray0"                => (  0,   0,   0),
    "gray1"                => (  3,   3,   3),
    "gray2"                => (  5,   5,   5),
    "gray3"                => (  8,   8,   8),
    "gray4"                => ( 10,  10,  10),
    "gray5"                => ( 13,  13,  13),
    "gray6"                => ( 15,  15,  15),
    "gray7"                => ( 18,  18,  18),
    "gray8"                => ( 20,  20,  20),
    "gray9"                => ( 23,  23,  23),
    "gray10"               => ( 26,  26,  26),
    "gray11"               => ( 28,  28,  28),
    "gray12"               => ( 31,  31,  31),
    "gray13"               => ( 33,  33,  33),
    "gray14"               => ( 36,  36,  36),
    "gray15"               => ( 38,  38,  38),
    "gray16"               => ( 41,  41,  41),
    "gray17"               => ( 43,  43,  43),
    "gray18"               => ( 46,  46,  46),
    "gray19"               => ( 48,  48,  48),
    "gray20"               => ( 51,  51,  51),
    "gray21"               => ( 54,  54,  54),
    "gray22"               => ( 56,  56,  56),
    "gray23"               => ( 59,  59,  59),
    "gray24"               => ( 61,  61,  61),
    "gray25"               => ( 64,  64,  64),
    "gray26"               => ( 66,  66,  66),
    "gray27"               => ( 69,  69,  69),
    "gray28"               => ( 71,  71,  71),
    "gray29"               => ( 74,  74,  74),
    "gray30"               => ( 77,  77,  77),
    "gray31"               => ( 79,  79,  79),
    "gray32"               => ( 82,  82,  82),
    "gray33"               => ( 84,  84,  84),
    "gray34"               => ( 87,  87,  87),
    "gray35"               => ( 89,  89,  89),
    "gray36"               => ( 92,  92,  92),
    "gray37"               => ( 94,  94,  94),
    "gray38"               => ( 97,  97,  97),
    "gray39"               => ( 99,  99,  99),
    "gray40"               => (102, 102, 102),
    "gray41"               => (105, 105, 105),
    "gray42"               => (107, 107, 107),
    "gray43"               => (110, 110, 110),
    "gray44"               => (112, 112, 112),
    "gray45"               => (115, 115, 115),
    "gray46"               => (117, 117, 117),
    "gray47"               => (120, 120, 120),
    "gray48"               => (122, 122, 122),
    "gray49"               => (125, 125, 125),
    "gray50"               => (127, 127, 127),
    "gray51"               => (130, 130, 130),
    "gray52"               => (133, 133, 133),
    "gray53"               => (135, 135, 135),
    "gray54"               => (138, 138, 138),
    "gray55"               => (140, 140, 140),
    "gray56"               => (143, 143, 143),
    "gray57"               => (145, 145, 145),
    "gray58"               => (148, 148, 148),
    "gray59"               => (150, 150, 150),
    "gray60"               => (153, 153, 153),
    "gray61"               => (156, 156, 156),
    "gray62"               => (158, 158, 158),
    "gray63"               => (161, 161, 161),
    "gray64"               => (163, 163, 163),
    "gray65"               => (166, 166, 166),
    "gray66"               => (168, 168, 168),
    "gray67"               => (171, 171, 171),
    "gray68"               => (173, 173, 173),
    "gray69"               => (176, 176, 176),
    "gray70"               => (179, 179, 179),
    "gray71"               => (181, 181, 181),
    "gray72"               => (184, 184, 184),
    "gray73"               => (186, 186, 186),
    "gray74"               => (189, 189, 189),
    "gray75"               => (191, 191, 191),
    "gray76"               => (194, 194, 194),
    "gray77"               => (196, 196, 196),
    "gray78"               => (199, 199, 199),
    "gray79"               => (201, 201, 201),
    "gray80"               => (204, 204, 204),
    "gray81"               => (207, 207, 207),
    "gray82"               => (209, 209, 209),
    "gray83"               => (212, 212, 212),
    "gray84"               => (214, 214, 214),
    "gray85"               => (217, 217, 217),
    "gray86"               => (219, 219, 219),
    "gray87"               => (222, 222, 222),
    "gray88"               => (224, 224, 224),
    "gray89"               => (227, 227, 227),
    "gray90"               => (229, 229, 229),
    "gray91"               => (232, 232, 232),
    "gray92"               => (235, 235, 235),
    "gray93"               => (237, 237, 237),
    "gray94"               => (240, 240, 240),
    "gray95"               => (242, 242, 242),
    "gray96"               => (245, 245, 245),
    "gray97"               => (247, 247, 247),
    "gray98"               => (250, 250, 250),
    "gray99"               => (252, 252, 252),
    "gray100"              => (255, 255, 255),
    "green"                => (  0, 128,   0),
    "green1"               => (  0, 255,   0),
    "green2"               => (  0, 238,   0),
    "green3"               => (  0, 205,   0),
    "green4"               => (  0, 139,   0),
    "greenyellow"          => (173, 255,  47),
    "grey"                 => (128, 128, 128),
    "grey0"                => (  0,   0,   0),
    "grey1"                => (  3,   3,   3),
    "grey2"                => (  5,   5,   5),
    "grey3"                => (  8,   8,   8),
    "grey4"                => ( 10,  10,  10),
    "grey5"                => ( 13,  13,  13),
    "grey6"                => ( 15,  15,  15),
    "grey7"                => ( 18,  18,  18),
    "grey8"                => ( 20,  20,  20),
    "grey9"                => ( 23,  23,  23),
    "grey10"               => ( 26,  26,  26),
    "grey11"               => ( 28,  28,  28),
    "grey12"               => ( 31,  31,  31),
    "grey13"               => ( 33,  33,  33),
    "grey14"               => ( 36,  36,  36),
    "grey15"               => ( 38,  38,  38),
    "grey16"               => ( 41,  41,  41),
    "grey17"               => ( 43,  43,  43),
    "grey18"               => ( 46,  46,  46),
    "grey19"               => ( 48,  48,  48),
    "grey20"               => ( 51,  51,  51),
    "grey21"               => ( 54,  54,  54),
    "grey22"               => ( 56,  56,  56),
    "grey23"               => ( 59,  59,  59),
    "grey24"               => ( 61,  61,  61),
    "grey25"               => ( 64,  64,  64),
    "grey26"               => ( 66,  66,  66),
    "grey27"               => ( 69,  69,  69),
    "grey28"               => ( 71,  71,  71),
    "grey29"               => ( 74,  74,  74),
    "grey30"               => ( 77,  77,  77),
    "grey31"               => ( 79,  79,  79),
    "grey32"               => ( 82,  82,  82),
    "grey33"               => ( 84,  84,  84),
    "grey34"               => ( 87,  87,  87),
    "grey35"               => ( 89,  89,  89),
    "grey36"               => ( 92,  92,  92),
    "grey37"               => ( 94,  94,  94),
    "grey38"               => ( 97,  97,  97),
    "grey39"               => ( 99,  99,  99),
    "grey40"               => (102, 102, 102),
    "grey41"               => (105, 105, 105),
    "grey42"               => (107, 107, 107),
    "grey43"               => (110, 110, 110),
    "grey44"               => (112, 112, 112),
    "grey45"               => (115, 115, 115),
    "grey46"               => (117, 117, 117),
    "grey47"               => (120, 120, 120),
    "grey48"               => (122, 122, 122),
    "grey49"               => (125, 125, 125),
    "grey50"               => (127, 127, 127),
    "grey51"               => (130, 130, 130),
    "grey52"               => (133, 133, 133),
    "grey53"               => (135, 135, 135),
    "grey54"               => (138, 138, 138),
    "grey55"               => (140, 140, 140),
    "grey56"               => (143, 143, 143),
    "grey57"               => (145, 145, 145),
    "grey58"               => (148, 148, 148),
    "grey59"               => (150, 150, 150),
    "grey60"               => (153, 153, 153),
    "grey61"               => (156, 156, 156),
    "grey62"               => (158, 158, 158),
    "grey63"               => (161, 161, 161),
    "grey64"               => (163, 163, 163),
    "grey65"               => (166, 166, 166),
    "grey66"               => (168, 168, 168),
    "grey67"               => (171, 171, 171),
    "grey68"               => (173, 173, 173),
    "grey69"               => (176, 176, 176),
    "grey70"               => (179, 179, 179),
    "grey71"               => (181, 181, 181),
    "grey72"               => (184, 184, 184),
    "grey73"               => (186, 186, 186),
    "grey74"               => (189, 189, 189),
    "grey75"               => (191, 191, 191),
    "grey76"               => (194, 194, 194),
    "grey77"               => (196, 196, 196),
    "grey78"               => (199, 199, 199),
    "grey79"               => (201, 201, 201),
    "grey80"               => (204, 204, 204),
    "grey81"               => (207, 207, 207),
    "grey82"               => (209, 209, 209),
    "grey83"               => (212, 212, 212),
    "grey84"               => (214, 214, 214),
    "grey85"               => (217, 217, 217),
    "grey86"               => (219, 219, 219),
    "grey87"               => (222, 222, 222),
    "grey88"               => (224, 224, 224),
    "grey89"               => (227, 227, 227),
    "grey90"               => (229, 229, 229),
    "grey91"               => (232, 232, 232),
    "grey92"               => (235, 235, 235),
    "grey93"               => (237, 237, 237),
    "grey94"               => (240, 240, 240),
    "grey95"               => (242, 242, 242),
    "grey96"               => (245, 245, 245),
    "grey97"               => (247, 247, 247),
    "grey98"               => (250, 250, 250),
    "grey99"               => (252, 252, 252),
    "grey100"              => (255, 255, 255),
    "honeydew"             => (240, 255, 240),
    "honeydew1"            => (240, 255, 240),
    "honeydew2"            => (224, 238, 224),
    "honeydew3"            => (193, 205, 193),
    "honeydew4"            => (131, 139, 131),
    "hotpink"              => (255, 105, 180),
    "hotpink1"             => (255, 110, 180),
    "hotpink2"             => (238, 106, 167),
    "hotpink3"             => (205,  96, 144),
    "hotpink4"             => (139,  58,  98),
    "indianred"            => (205,  92,  92),
    "indianred1"           => (255, 106, 106),
    "indianred2"           => (238,  99,  99),
    "indianred3"           => (205,  85,  85),
    "indianred4"           => (139,  58,  58),
    "indigo"               => ( 75,   0, 130),
    "ivory"                => (255, 255, 240),
    "ivory1"               => (255, 255, 240),
    "ivory2"               => (238, 238, 224),
    "ivory3"               => (205, 205, 193),
    "ivory4"               => (139, 139, 131),
    "khaki"                => (240, 230, 140),
    "khaki1"               => (255, 246, 143),
    "khaki2"               => (238, 230, 133),
    "khaki3"               => (205, 198, 115),
    "khaki4"               => (139, 134,  78),
    "lavender"             => (230, 230, 250),
    "lavenderblush"        => (255, 240, 245),
    "lavenderblush1"       => (255, 240, 245),
    "lavenderblush2"       => (238, 224, 229),
    "lavenderblush3"       => (205, 193, 197),
    "lavenderblush4"       => (139, 131, 134),
    "lawngreen"            => (124, 252,   0),
    "lemonchiffon"         => (255, 250, 205),
    "lemonchiffon1"        => (255, 250, 205),
    "lemonchiffon2"        => (238, 233, 191),
    "lemonchiffon3"        => (205, 201, 165),
    "lemonchiffon4"        => (139, 137, 112),
    "lightblue"            => (173, 216, 230),
    "lightblue1"           => (191, 239, 255),
    "lightblue2"           => (178, 223, 238),
    "lightblue3"           => (154, 192, 205),
    "lightblue4"           => (104, 131, 139),
    "lightcoral"           => (240, 128, 128),
    "lightcyan"            => (224, 255, 255),
    "lightcyan1"           => (224, 255, 255),
    "lightcyan2"           => (209, 238, 238),
    "lightcyan3"           => (180, 205, 205),
    "lightcyan4"           => (122, 139, 139),
    "lightgoldenrod"       => (238, 221, 130),
    "lightgoldenrod1"      => (255, 236, 139),
    "lightgoldenrod2"      => (238, 220, 130),
    "lightgoldenrod3"      => (205, 190, 112),
    "lightgoldenrod4"      => (139, 129,  76),
    "lightgoldenrodyellow" => (250, 250, 210),
    "lightgray"            => (211, 211, 211),
    "lightgreen"           => (144, 238, 144),
    "lightgrey"            => (211, 211, 211),
    "lightpink"            => (255, 182, 193),
    "lightpink1"           => (255, 174, 185),
    "lightpink2"           => (238, 162, 173),
    "lightpink3"           => (205, 140, 149),
    "lightpink4"           => (139,  95, 101),
    "lightsalmon"          => (255, 160, 122),
    "lightsalmon1"         => (255, 160, 122),
    "lightsalmon2"         => (238, 149, 114),
    "lightsalmon3"         => (205, 129,  98),
    "lightsalmon4"         => (139,  87,  66),
    "lightseagreen"        => ( 32, 178, 170),
    "lightskyblue"         => (135, 206, 250),
    "lightskyblue1"        => (176, 226, 255),
    "lightskyblue2"        => (164, 211, 238),
    "lightskyblue3"        => (141, 182, 205),
    "lightskyblue4"        => ( 96, 123, 139),
    "lightslateblue"       => (132, 112, 255),
    "lightslategray"       => (119, 136, 153),
    "lightslategrey"       => (119, 136, 153),
    "lightsteelblue"       => (176, 196, 222),
    "lightsteelblue1"      => (202, 225, 255),
    "lightsteelblue2"      => (188, 210, 238),
    "lightsteelblue3"      => (162, 181, 205),
    "lightsteelblue4"      => (110, 123, 139),
    "lightyellow"          => (255, 255, 224),
    "lightyellow1"         => (255, 255, 224),
    "lightyellow2"         => (238, 238, 209),
    "lightyellow3"         => (205, 205, 180),
    "lightyellow4"         => (139, 139, 122),
    "lime"                 => (  0, 255,   0),
    "limegreen"            => ( 50, 205,  50),
    "linen"                => (250, 240, 230),
    "magenta"              => (255,   0, 255),
    "magenta1"             => (255,   0, 255),
    "magenta2"             => (238,   0, 238),
    "magenta3"             => (205,   0, 205),
    "magenta4"             => (139,   0, 139),
    "maroon"               => (128,   0,   0),
    "maroon1"              => (255,  52, 179),
    "maroon2"              => (238,  48, 167),
    "maroon3"              => (205,  41, 144),
    "maroon4"              => (139,  28,  98),
    "mediumaquamarine"     => (102, 205, 170),
    "mediumblue"           => (  0,   0, 205),
    "mediumorchid"         => (186,  85, 211),
    "mediumorchid1"        => (224, 102, 255),
    "mediumorchid2"        => (209,  95, 238),
    "mediumorchid3"        => (180,  82, 205),
    "mediumorchid4"        => (122,  55, 139),
    "mediumpurple"         => (147, 112, 219),
    "mediumpurple1"        => (171, 130, 255),
    "mediumpurple2"        => (159, 121, 238),
    "mediumpurple3"        => (137, 104, 205),
    "mediumpurple4"        => ( 93,  71, 139),
    "mediumseagreen"       => ( 60, 179, 113),
    "mediumslateblue"      => (123, 104, 238),
    "mediumspringgreen"    => (  0, 250, 154),
    "mediumturquoise"      => ( 72, 209, 204),
    "mediumvioletred"      => (199,  21, 133),
    "midnightblue"         => ( 25,  25, 112),
    "mintcream"            => (245, 255, 250),
    "mistyrose"            => (255, 228, 225),
    "mistyrose1"           => (255, 228, 225),
    "mistyrose2"           => (238, 213, 210),
    "mistyrose3"           => (205, 183, 181),
    "mistyrose4"           => (139, 125, 123),
    "moccasin"             => (255, 228, 181),
    "navajowhite"          => (255, 222, 173),
    "navajowhite1"         => (255, 222, 173),
    "navajowhite2"         => (238, 207, 161),
    "navajowhite3"         => (205, 179, 139),
    "navajowhite4"         => (139, 121,  94),
    "navy"                 => (  0,   0, 128),
    "navyblue"             => (  0,   0, 128),
    "oldlace"              => (253, 245, 230),
    "olive"                => (128, 128,   0),
    "olivedrab"            => (107, 142,  35),
    "olivedrab1"           => (192, 255,  62),
    "olivedrab2"           => (179, 238,  58),
    "olivedrab3"           => (154, 205,  50),
    "olivedrab4"           => (105, 139,  34),
    "orange"               => (255, 165,   0),
    "orange1"              => (255, 165,   0),
    "orange2"              => (238, 154,   0),
    "orange3"              => (205, 133,   0),
    "orange4"              => (139,  90,   0),
    "orangered"            => (255,  69,   0),
    "orangered1"           => (255,  69,   0),
    "orangered2"           => (238,  64,   0),
    "orangered3"           => (205,  55,   0),
    "orangered4"           => (139,  37,   0),
    "orchid"               => (218, 112, 214),
    "orchid1"              => (255, 131, 250),
    "orchid2"              => (238, 122, 233),
    "orchid3"              => (205, 105, 201),
    "orchid4"              => (139,  71, 137),
    "palegoldenrod"        => (238, 232, 170),
    "palegreen"            => (152, 251, 152),
    "palegreen1"           => (154, 255, 154),
    "palegreen2"           => (144, 238, 144),
    "palegreen3"           => (124, 205, 124),
    "palegreen4"           => ( 84, 139,  84),
    "paleturquoise"        => (175, 238, 238),
    "paleturquoise1"       => (187, 255, 255),
    "paleturquoise2"       => (174, 238, 238),
    "paleturquoise3"       => (150, 205, 205),
    "paleturquoise4"       => (102, 139, 139),
    "palevioletred"        => (219, 112, 147),
    "palevioletred1"       => (255, 130, 171),
    "palevioletred2"       => (238, 121, 159),
    "palevioletred3"       => (205, 104, 137),
    "palevioletred4"       => (139,  71,  93),
    "papayawhip"           => (255, 239, 213),
    "peachpuff"            => (255, 218, 185),
    "peachpuff1"           => (255, 218, 185),
    "peachpuff2"           => (238, 203, 173),
    "peachpuff3"           => (205, 175, 149),
    "peachpuff4"           => (139, 119, 101),
    "peru"                 => (205, 133,  63),
    "pink"                 => (255, 192, 203),
    "pink1"                => (255, 181, 197),
    "pink2"                => (238, 169, 184),
    "pink3"                => (205, 145, 158),
    "pink4"                => (139,  99, 108),
    "plum"                 => (221, 160, 221),
    "plum1"                => (255, 187, 255),
    "plum2"                => (238, 174, 238),
    "plum3"                => (205, 150, 205),
    "plum4"                => (139, 102, 139),
    "powderblue"           => (176, 224, 230),
    "purple"               => (128,   0, 128),
    "purple1"              => (155,  48, 255),
    "purple2"              => (145,  44, 238),
    "purple3"              => (125,  38, 205),
    "purple4"              => ( 85,  26, 139),
    "rebeccapurple"        => (102,  51, 153),
    "red"                  => (255,   0,   0),
    "red1"                 => (255,   0,   0),
    "red2"                 => (238,   0,   0),
    "red3"                 => (205,   0,   0),
    "red4"                 => (139,   0,   0),
    "rosybrown"            => (188, 143, 143),
    "rosybrown1"           => (255, 193, 193),
    "rosybrown2"           => (238, 180, 180),
    "rosybrown3"           => (205, 155, 155),
    "rosybrown4"           => (139, 105, 105),
    "royalblue"            => ( 65, 105, 225),
    "royalblue1"           => ( 72, 118, 255),
    "royalblue2"           => ( 67, 110, 238),
    "royalblue3"           => ( 58,  95, 205),
    "royalblue4"           => ( 39,  64, 139),
    "saddlebrown"          => (139,  69,  19),
    "salmon"               => (250, 128, 114),
    "salmon1"              => (255, 140, 105),
    "salmon2"              => (238, 130,  98),
    "salmon3"              => (205, 112,  84),
    "salmon4"              => (139,  76,  57),
    "sandybrown"           => (244, 164,  96),
    "seagreen"             => ( 46, 139,  87),
    "seagreen1"            => ( 84, 255, 159),
    "seagreen2"            => ( 78, 238, 148),
    "seagreen3"            => ( 67, 205, 128),
    "seagreen4"            => ( 46, 139,  87),
    "seashell"             => (255, 245, 238),
    "seashell1"            => (255, 245, 238),
    "seashell2"            => (238, 229, 222),
    "seashell3"            => (205, 197, 191),
    "seashell4"            => (139, 134, 130),
    "sienna"               => (160,  82,  45),
    "sienna1"              => (255, 130,  71),
    "sienna2"              => (238, 121,  66),
    "sienna3"              => (205, 104,  57),
    "sienna4"              => (139,  71,  38),
    "silver"               => (192, 192, 192),
    "skyblue"              => (135, 206, 235),
    "skyblue1"             => (135, 206, 255),
    "skyblue2"             => (126, 192, 238),
    "skyblue3"             => (108, 166, 205),
    "skyblue4"             => ( 74, 112, 139),
    "slateblue"            => (106,  90, 205),
    "slateblue1"           => (131, 111, 255),
    "slateblue2"           => (122, 103, 238),
    "slateblue3"           => (105,  89, 205),
    "slateblue4"           => ( 71,  60, 139),
    "slategray"            => (112, 128, 144),
    "slategray1"           => (198, 226, 255),
    "slategray2"           => (185, 211, 238),
    "slategray3"           => (159, 182, 205),
    "slategray4"           => (108, 123, 139),
    "slategrey"            => (112, 128, 144),
    "snow"                 => (255, 250, 250),
    "snow1"                => (255, 250, 250),
    "snow2"                => (238, 233, 233),
    "snow3"                => (205, 201, 201),
    "snow4"                => (139, 137, 137),
    "springgreen"          => (  0, 255, 127),
    "springgreen1"         => (  0, 255, 127),
    "springgreen2"         => (  0, 238, 118),
    "springgreen3"         => (  0, 205, 102),
    "springgreen4"         => (  0, 139,  69),
    "steelblue"            => ( 70, 130, 180),
    "steelblue1"           => ( 99, 184, 255),
    "steelblue2"           => ( 92, 172, 238),
    "steelblue3"           => ( 79, 148, 205),
    "steelblue4"           => ( 54, 100, 139),
    "tan"                  => (210, 180, 140),
    "tan1"                 => (255, 165,  79),
    "tan2"                 => (238, 154,  73),
    "tan3"                 => (205, 133,  63),
    "tan4"                 => (139,  90,  43),
    "teal"                 => (  0, 128, 128),
    "thistle"              => (216, 191, 216),
    "thistle1"             => (255, 225, 255),
    "thistle2"             => (238, 210, 238),
    "thistle3"             => (205, 181, 205),
    "thistle4"             => (139, 123, 139),
    "tomato"               => (255,  99,  71),
    "tomato1"              => (255,  99,  71),
    "tomato2"              => (238,  92,  66),
    "tomato3"              => (205,  79,  57),
    "tomato4"              => (139,  54,  38),
    "turquoise"            => ( 64, 224, 208),
    "turquoise1"           => (  0, 245, 255),
    "turquoise2"           => (  0, 229, 238),
    "turquoise3"           => (  0, 197, 205),
    "turquoise4"           => (  0, 134, 139),
    "violet"               => (238, 130, 238),
    "violetred"            => (208,  32, 144),
    "violetred1"           => (255,  62, 150),
    "violetred2"           => (238,  58, 140),
    "violetred3"           => (205,  50, 120),
    "violetred4"           => (139,  34,  82),
    "wheat"                => (245, 222, 179),
    "wheat1"               => (255, 231, 186),
    "wheat2"               => (238, 216, 174),
    "wheat3"               => (205, 186, 150),
    "wheat4"               => (139, 126, 102),
    "white"                => (255, 255, 255),
    "whitesmoke"           => (245, 245, 245),
    "yellow"               => (255, 255,   0),
    "yellow1"              => (255, 255,   0),
    "yellow2"              => (238, 238,   0),
    "yellow3"              => (205, 205,   0),
    "yellow4"              => (139, 139,   0),
    "yellowgreen"          => (154, 205,  50)
)