text:
module ConstructiveSolidGeometry

export Coord
export Ray
export Surface
export Plane
export Cone
export Sphere
export InfCylinder
export Box
export Region
export Cell
export Geometry
export +,-,*,^,|,~
export reflect
export generate_random_ray
export raytrace
export find_intersection
export halfspace
export is_in_cell
export find_cell_id
export plot_geometry_2D
export plot_cell_2D
export dot
export magnitude
export unitize
export cross

using Plots

"""
    type Coord

An {x,y,z} coordinate type. Used throughout the ConstructiveSolidGeometry.jl package for speed.

# Constructors
* `Coord(x::Float64, y::Float64, z::Float64)`
"""
type Coord
    x::Float64
    y::Float64
    z::Float64
end

"""
    type Ray

A ray is defined by its origin and a unitized direction vector

# Constructors
* `Ray(origin::Coord, direction::Coord)`
"""
type Ray
    origin::Coord
    direction::Coord
end

"""
    abstract Surface

An abstract class that all surfaces (`Sphere`, `Plane`, `InfCylinder`) inherit from. Implementation of new shapes should inherit from `Surface`.
"""
abstract Surface

"""
    type Plane <: Surface

Defined by a point on the surface of the plane, its unit normal vector, and an optional boundary condition.

# Constructors
* `Plane(point::Coord, normal::Coord)`
* `Plane(point::Coord, normal::Coord, boundary::String)`

# Arguments
* `point::Coord`: Any point on the surface of the plane
* `normal::Coord`: A unit normal vector of the plane. Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are "transmission" (default), "vacuum", and "reflective".
"""
type Plane <: Surface
    point::Coord
    normal::Coord
    reflective::Bool
	vacuum::Bool
	Plane(point::Coord, normal::Coord, ref::Bool, vac::Bool) = new(point, normal, ref, vac)
	function Plane(point::Coord, normal::Coord, boundary::String)
		if boundary == "reflective"
			new(point, normal, true, false)
		elseif boundary == "vacuum"
			new(point, normal, false, true)
		else
			new(point, normal, false, false)
		end
	end
	Plane(point::Coord, normal::Coord) = new(point, normal, false, false)
end

"""
    type Cone <: Surface

Defined by the tip of the cone, its direction axis vector, the angle between the central axis and the cone surface, and an optional boundary condition.

# Constructors
* `Cone(tip::Coord, axis::Coord, theta::Float64)`
* `Cone(tip::Coord, axis::Coord, theta::Float64, boundary::String)`

# Arguments
* `tip::Coord`: The vertex (tip) of the cone
* `axis::Coord`: A unit vector representing the central axis of the cone. As the cone equation actually defines two cones eminating in a mirrored fashion from the tip, this direction vector also indicates which cone is the true cone. I.e., following the direction of the axis vector when starting from the tip should lead inside the cone you actually want.  Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `theta::Float64`: The angle (in radians) between the central axis (must be between 0 and pi/2)
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type Cone <: Surface
    tip::Coord
    axis::Coord
	theta::Float64
    reflective::Bool
	vacuum::Bool
	Cone(tip::Coord, axis::Coord, theta::Float64, ref::Bool, vac::Bool) = new(tip, axis, theta, ref, vac)
	function Cone(tip::Coord, axis::Coord, theta::Float64, boundary::String)
		if boundary == "reflective"
			new(tip, axis, theta, true, false)
		elseif boundary == "vacuum"
			new(tip, axis, theta, false, true)
		else
			new(tip, axis, theta, false, false)
		end
	end
	Cone(tip::Coord, axis::Coord, theta::Float64) = new(tip, axis, theta, false, false)
end

"""
    type Sphere <: Surface

Defined by the center of the sphere, its radius, and an optional boundary condition.

# Constructors
* `Sphere(center::Coord, radius::Float64)`
* `Sphere(center::Coord, radius::Float64, boundary::String)`

# Arguments
* `center::Coord`: The center of the sphere
* `radius::Float64`: The radius of the sphere
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type Sphere <: Surface
    center::Coord
    radius::Float64
    reflective::Bool
	vacuum::Bool
	Sphere(c::Coord, r::Float64, ref::Bool, vac::Bool) = new(c, r, ref, vac)
	function Sphere(c::Coord, r::Float64, boundary::String)
		if boundary == "reflective"
			new(c, r, true, false)
		elseif boundary == "vacuum"
			new(c, r, false, true)
		else
			new(c, r, false, false)
		end
	end
	Sphere(c::Coord, r::Float64) = new(c, r, false, false)
end

"""
    type InfCylinder <: Surface

An arbitrary direction infinite cylinder defined by any point on its central axis, its radius, the unit normal direction of the cylinder, and an optional boundary condition. A finite cylinder can be generated by defining the intersection of an infinite cylinder and two planes.

# Constructors
* `InfCylinder(center::Coord, normal::Coord, radius::Float64)`
* `InfCylinder(center::Coord, normal::Coord, radius::Float64, boundary::String)`

# Arguments
* `center::Coord`: The center of the infinite cylinder
* `normal::Coord`: A unit normal direction vector of the cylinder (i.e., a vector along its central axis), Recommended to use `unitize(c::Coord)` if normalizing is needed.
* `radius::Float64`: The radius of the infinite cylinder
* `boundary::String`: Optional boundary condition, defined as a `String`. Options are \"transmission\" (default) or \"vacuum\".
"""
type InfCylinder <: Surface
    center::Coord
    normal::Coord
    radius::Float64
    reflective::Bool
	vacuum::Bool
	InfCylinder(c::Coord, n::Coord, r::Float64, ref::Bool, vac::Bool) = new(c, n, r, ref, vac)
	function InfCylinder(c::Coord, n::Coord, r::Float64, boundary::String)
		if boundary == "reflective"
			new(c, n, r, true, false)
		elseif boundary == "vacuum"
			new(c, n, r, false, true)
		else
			new(c, n, r, false, false)
		end
	end
	InfCylinder(c::Coord, n::Coord, r::Float64) = new(c, n, r, false, false)
end

"""
    type Box

An axis aligned box is defined by the minimum `Coord` and maximum `Coord` of the box. Note that a Box is only used by ConstructiveSolidGeometry.jl for bounding box purposes, and is not a valid surface to define CSG cells with. Instead, you must define all six planes of a box independently.

# Constructors
* `Box(min::Coord, max::Coord)`
"""
type Box
    lower_left::Coord
    upper_right::Coord
end

"""
    type Region

The volume that is defined by a surface and one of its halfspaces

# Constructors
* `Region(surface::Surface, halfspace::Int64)`

# Arguments
* `surface::Surface`: A `Sphere`, `Plane`, or `InfCylinder`
* `halfspace::Int64`: Either +1 or -1
"""
type Region
    surface::Surface
    halfspace::Int64
end

"""
    type Cell

Defined by an array of regions and the logical combination of those regions that define the cell

# Constructors
* `Cell(regions::Array{Region}, definition::Expr)`

# Arguments
* `regions::Array{Region}`: An array of regions that are used to define the cell
* `definition::Expr`: A logical expression that defines the volume of the cell. The intersection operator is ^, the union operator is |, and the complement operator is ~. Regions are defined by their integer indices in the regions array.
"""
type Cell
    regions::Array{Region}
    definition::Expr
end

"""
    type Geometry

The top level object that holds all the cells in the problem. This object contains all data regarding the geometry within a system.

# Constructors
* `Geometry(cells::Array{Cell}, bounding_box::Box)`

# Arguments
* `cells::Array{Cell}`: All cells inside the geometry. The cells must combine to fill the entire space of the bounding box. No two cells should overlap.
* `bounding_box::Box`: The bounding box around the problem.
"""
type Geometry
    cells::Array{Cell}
    bounding_box::Box
end

_p = Coord(0,0,0)
typeassert(_p, Coord)

import Base: +, -, *, ^, |, ~
+(a::Coord, b::Coord)     = Coord(a.x+b.x, a.y+b.y, a.z+b.z)
-(a::Coord, b::Coord)     = Coord(a.x-b.x, a.y-b.y, a.z-b.z)
*(a::Float64, b::Coord)   = Coord(a*b.x, a*b.y, a*b.z)
*(b::Coord, a::Float64,)  = Coord(a*b.x, a*b.y, a*b.z)
*(a::Int, b::Coord)       = Coord(a*b.x, a*b.y, a*b.z)
*(b::Coord, a::Int)       = Coord(a*b.x, a*b.y, a*b.z)
dot(a::Coord, b::Coord)   = (a.x*b.x + a.y*b.y + a.z*b.z)
"""
    magnitude(a::Coord)

A utility function to determine the magnitude of a `Coord` object. Typical use case is to subtract two Coord objects and check the resulting Coord object's magnitude to determine the distance between the two Coords.
"""
magnitude(a::Coord)       = sqrt(dot(a,a))
"""
    unitize(a::Coord)

A utility function to unitize a `Coord`
"""
unitize(a::Coord)         = (1. / magnitude(a) * a)
cross(a::Coord, b::Coord) = Coord(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x)

# Ray - 3D Plane Intersection
# Returns: hit, distance
# hit: a boolean indicating if an intersection occurred
# distance: the distance to intersection
# Edge case policies:
#    1. Ray is inside the plane: Returns (true, NaN)
#    2. Ray is parallel to the plan, but not inside: Returns (false, NaN)
#    3. Ray never hits plane: Returns (false, NaN)
"""
    function raytrace(ray::Ray, surface::Surface)

Determines if a `Ray` and a `Surface` intersect, and the distance to that intersection.

# Returns
* `Bool`: Indicates if the ray intersects the surface or not
* `Float64`: The distance between the ray's origin and the point of intersection
"""
function raytrace(ray::Ray, plane::Plane)
    dist::Float64 = dot( plane.point - ray.origin, plane.normal) / dot( ray.direction, plane.normal)
    # Check if parallel
    if dist < 0 || dist == Inf
        return false, NaN
    end
    return true, dist
end

# Ray - 3D Cone Intersection
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, cone::Cone)

	cos_theta_squared::Float64 = (cos(cone.theta))^2
	CO::Coord = ray.origin - cone.tip

	a::Float64 = dot(ray.direction, cone.axis)^2 - cos_theta_squared
	b::Float64 = 2.0 * (dot(ray.direction, cone.axis) * dot(CO, cone.axis) - dot(ray.direction, CO) * cos_theta_squared)
	c::Float64 = dot(CO, cone.axis)^2 - dot(CO, CO) * cos_theta_squared

	determinant::Float64 = b^2 - 4.0*a*c

	if determinant < 0
		return false, NaN
	elseif determinant == 0
		return true, -b/(2*a)
	end

	# Now we need to verify we are not intersecting with the shadow cone
	one_over_two_a::Float64 = 1.0 / (2.0 * a)
	t1::Float64 = (-b - sqrt(determinant)) * one_over_two_a
	t2::Float64 = (-b + sqrt(determinant)) * one_over_two_a

	p1::Coord = ray.origin + ray.direction * t1
	p2::Coord = ray.origin + ray.direction * t2

	valid1::Bool = false
	valid2::Bool = false

	# Check against shadow cone & ensure intersection is in front of ray
	if t1 >= 0 && dot((p1-cone.tip),cone.axis) > 0
		valid1 = true
	end
	if t2 >= 0 && dot((p2-cone.tip),cone.axis) > 0
		valid2 = true
	end

	# If both points hit real cone, return closer one
	if valid1 && valid2
		if t1 < t2
			return true, t1
		else
			return true, t2
		end
	end

	if valid1 && !valid2
		return true, t1
	end

	if !valid1 && valid2
		return true, t2
	end

	return false, NaN
end

# Ray - 3D Sphere Intersection
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, sphere::Sphere)
    d::Coord = ray.origin - sphere.center
    t::Float64 = -dot(ray.direction, d)
    discriminant::Float64 = t^2
    discriminant -= magnitude(d)^2
    discriminant += sphere.radius^2

    # If the discriminant is less than zero, they don't hit
    if discriminant < 0
        return false, Inf
    end
    sqrt_val::Float64 = sqrt(discriminant)
    pos::Float64 = t - sqrt_val
    neg::Float64 = t + sqrt_val

    if pos < 0  && neg < 0
        return false, NaN
    end
    if pos < 0 && neg > 0
        return true, neg
    end
    if pos < neg && pos > 0
        return true, pos
    end

    return true, neg
end

# Ray - 3D Infinite Cylinder Intersection (works for cylinder direction)
# Returns hit, distance
# hit: a boolean indicating if an intersection occurred (false if parallel or negative)
# dist: distance to closest intersection point
function raytrace(ray::Ray, infcylinder::InfCylinder)
    A = infcylinder.center
    # Generate point new point in cylinder for math
    B = infcylinder.center + infcylinder.normal
    O = ray.origin
    r = infcylinder.radius
    AB = B - A
    AO = O - A
    AOxAB = cross(AO, AB)
    VxAB = cross(ray.direction, AB)
    ab2::Float64 = dot(AB, AB)
    a::Float64 = dot(VxAB, VxAB)
    b::Float64 = 2.0 * dot(VxAB, AOxAB)
    c::Float64 = dot(AOxAB, AOxAB) - (r*r * ab2)

    # Check Determininant
    det::Float64 = b^2 - 4.0 * a * c
    if det < 0
        return false, Inf
    end

    pos::Float64 = (-b + sqrt(det)) / (2.0 * a)
    neg::Float64 = (-b - sqrt(det)) / (2.0 * a)

    if pos < 0
        if neg < 0
            return false, NaN
        end
        return true, neg
    end
    if neg < 0
        return true, pos
    end
    if pos < neg
        return true, pos
    else
        return true, neg
    end
end

"""
    reflect(ray::Ray, plane::Plane)

Reflects a ray off a plane.

# Return
* `Ray`: A new ray with the same origin as input, but with the new reflected direction
"""
function reflect(ray::Ray, plane::Plane)
    a = dot(ray.direction, plane.normal)
    b = plane.normal * (2.0 * a)
    c = ray.direction - b
    reflected_ray::Ray = Ray(ray.origin, c)
    return reflected_ray
end

"""
    generate_random_ray(box::Box)

Returns a randomly sampled ray from within an axis aligned bounding box.
"""
function generate_random_ray(box::Box)
    ray = Ray(Coord(0.0, 0.0, 0.0), Coord(0.0, 0.0, 0.0))

    # Sample Origin
    width = Coord(box.upper_right.x - box.lower_left.x, box.upper_right.y - box.lower_left.y, box.upper_right.z - box.lower_left.z)
    ray.origin.x = box.lower_left.x + rand(Float64)*width.x
    ray.origin.y = box.lower_left.y + rand(Float64)*width.y
    ray.origin.z = box.lower_left.z + rand(Float64)*width.z

    # Sample Direction From Sphere
    theta::Float64 = rand(Float64) * 2.0 * pi
    z::Float64 = -1.0 + 2.0 * rand(Float64)
    zo::Float64 = sqrt(1.0 - z*z)
    ray.direction.x = zo * cos(theta);
    ray.direction.y = zo * sin(theta);
    ray.direction.z = z;

    # Normalize Direction
    ray.direction = unitize(ray.direction)

    return ray
end

"""
    find_intersection(ray::Ray, regions::Array{Region})

Performs ray tracing on an array of regions.

# Return
* `Ray`: A new Ray that has been moved just accross the point of intersection.
* `Int64`: The surface id that was hit.
* `String`: The boundary condition of the surface that was hit.
"""
function find_intersection(ray::Ray, regions::Array{Region})
    BUMP::Float64 = 1.0e-9
    min::Float64 = 1e30
    id::Int64 = -1
    for i = 1:length(regions)
        hit, dist = raytrace(ray, regions[i].surface)
        if hit == true
            if dist < min
                min = dist
                id = i
            end
        end
    end

    new_ray::Ray = Ray(ray.origin + ray.direction * (min + BUMP), ray.direction)

    if regions[id].surface.reflective == true
        new_ray = reflect(new_ray, regions[id].surface)
        new_ray.origin = new_ray.origin + new_ray.direction * (2.0 * BUMP)
		return new_ray, id, "reflective"
    end

	if regions[id].surface.vacuum == true
		return new_ray, id, "vacuum"
    end

    return new_ray, id, "transmission"

end

"""
    find_intersection(ray::Ray, geometry::Geometry)

Performs ray tracing on a Geometry

# Return
* `Ray`: A new Ray that has been moved just accross the point of intersection.
* `Int64`: The surface id that was hit.
* `String`: The boundary condition of the surface that was hit.
"""
function find_intersection(ray::Ray, geometry::Geometry)
	cell_id = find_cell_id(ray.origin, geometry)
	regions::Array{Region} = geometry.cells[cell_id].regions
	return find_intersection(ray, regions)
end

# The halfpsace functions are private methods, and all take a coordinate and a surface to determine which halfspace the coordinate is in.

# Plane halfspace determination
function halfspace(c::Coord, plane::Plane)
    d::Float64 = -dot(plane.normal, plane.point)
    half::Float64 = dot(plane.normal, c) + d
    if half <= 0
        return -1
    else
        return 1
    end
end

# Cone halfspace determination
function halfspace(c::Coord, cone::Cone)
	p_minus_c::Coord = c - cone.tip
    half::Float64 = dot(p_minus_c,cone.axis)^2 - dot(p_minus_c,p_minus_c) * (cos(cone.theta))^2
    if half <= 0
        return -1
    else
        return 1
    end
end

# Sphere halfspace determination
function halfspace(c::Coord, sphere::Sphere)
    half::Float64 = (c.x - sphere.center.x)^2 + (c.y - sphere.center.y)^2 + (c.z - sphere.center.z)^2 - sphere.radius^2
    if half <= 0
        return -1
    else
        return 1
    end
end

# Infinite cylinder halfspace
function halfspace(c::Coord, cyl::InfCylinder)
    tmp::Coord = cross((c-cyl.center), cyl.normal)
    half::Float64 = dot(tmp, tmp) - cyl.radius^2
    if half <= 0
        return -1
    else
        return 1
    end
end

function ^(a::Region, b::Region)
    if halfspace(_p, a.surface) == a.halfspace
        if halfspace(_p, b.surface) == b.halfspace
            return true
        end
    end
    return false
end

function ^(a::Region, b::Bool)
    if halfspace(_p, a.surface) == a.halfspace
        if b == true
            return true
        end
    end
    return false
end

function ^(b::Bool, a::Region)
    if halfspace(_p, a.surface) == a.halfspace
        if b == true
            return true
        end
    end
    return false
end

function |(a::Region, b::Region)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if halfspace(_p, b.surface) == b.halfspace
        return true
    end
    return false
end

function |(a::Region, b::Bool)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if b == true
        return true
    end
    return false
end

function |(b::Bool, a::Region)
    if halfspace(_p, a.surface) == a.halfspace
        return true
    end
    if b == true
        return true
    end
    return false
end

function ~(a::Region)
    b::Region = Region(a.surface, a.halfspace)
    if a.halfspace == -1
        b.halfspace = 1
    else
        b.halfspace = -1
    end
    return b
end

"""
    is_in_cell(p::Coord, cell::Cell)

Determines if a point (such as a Ray origin) is inside a given cell
"""
function is_in_cell(p::Coord, cell::Cell)
    result = navigate_tree(p, cell.regions, cell.definition)
    return result
end

function navigate_tree(p::Coord, r::Array{Region}, ex::Expr)
    global _p = Coord(p.x, p.y, p.z)

	# Check if Complement
	if ex.args[1] == :~
		if typeof(ex.args[2]) == typeof(1)
			return ~ r[ex.args[2]]
		else
			return ~ navigate_tree(p, r, ex.args[2])
		end
	end

	if typeof(ex.args[2]) == typeof(1)
		# Case 1 - Both operands are leaves
		if typeof(ex.args[3]) == typeof(1)
			if ex.args[1] == :^
            	return r[ex.args[2]] ^ r[ex.args[3]]
			end
			if ex.args[1] == :|
            	return r[ex.args[2]] | r[ex.args[3]]
			end
		end
		# Case 2 - Left operand is leaf, right is not
		if typeof(ex.args[3]) != typeof(1)
			if ex.args[1] == :^
            	return r[ex.args[2]] ^ navigate_tree(p, r, ex.args[3])
			end
			if ex.args[1] == :|
            	return r[ex.args[2]] | navigate_tree(p, r, ex.args[3])
			end
		end
	end

	if typeof(ex.args[2]) != typeof(1)
		# Case 3 - left operand is not leaf, but right is
		if typeof(ex.args[3]) == typeof(1)
			if ex.args[1] == :^
            	return navigate_tree(p, r, ex.args[2]) ^ r[ex.args[3]]
			end
			if ex.args[1] == :|
            	return navigate_tree(p, r, ex.args[2]) | r[ex.args[3]]
			end
		end
		# Case 4 - Neither operand is a leaf
		if typeof(ex.args[3]) != typeof(1)
			if ex.args[1] == :^
            	return navigate_tree(p, r, ex.args[2]) ^ navigate_tree(p, r, ex.args[3])
			end
			if ex.args[1] == :|
            	return navigate_tree(p, r, ex.args[2]) | navigate_tree(p, r, ex.args[3])
			end
		end
	end
end

"""
    find_cell_id(p::Coord, geometry::Geometry)

Finds the cell id that a point resides within
"""
function find_cell_id(p::Coord, geometry::Geometry)
    for i = 1:length(geometry.cells)
        if is_in_cell(p, geometry.cells[i]) == true
            return i
        end
    end
    return -1
end

"""
    plot_geometry_2D(geometry::Geometry, view::Box, dim::Int64)

Plots a 2D x-y slice of a geometry.

# Arguments
* `geometry::Geometry`: the geometry we want to plot
* `view::Box`: The view box is an axis aligned box that defines where the picture will be taken, with both min and max z dimensions indicating the single z elevation the slice is taken at.
* `dim::Int64`: The dimension is the number of pixels along the x and y axis to use, which determines the resolution of the picture.
"""
function plot_geometry_2D(geometry::Geometry, view::Box, dim::Int64)
    delta_x = (view.upper_right.x - view.lower_left.x) / (dim)
    delta_y = (view.upper_right.y - view.lower_left.y) / (dim)

    x_coords = collect(view.lower_left.x + delta_x/2.0:delta_x:view.upper_right.x - delta_x/2.0)
    y_coords = collect(view.lower_left.y + delta_y/2.0:delta_y:view.upper_right.y - delta_y/2.0)

    pixels = Array{Int64, 2}(dim, dim)

    for i=1:dim
        for j=1:dim
            pixels[i,j] = find_cell_id(Coord(x_coords[i], y_coords[j], view.lower_left.z), geometry)
        end
    end
    pixels = rotl90(pixels)
    colors = Array{RGBA}(0)
    for i=1:length(geometry.cells)
        if (i-1)%4 == 0
            push!(colors, RGBA(1.0, 0.0, 0.0, 1.0))
        elseif (i-1)%4 == 1
            push!(colors, RGBA(0, 1.0, 0, 1.0))
        elseif (i-1)%4 == 2
            push!(colors, RGBA(1.0, 0.0, 1.0, 1.0))
        elseif (i-1)%4 == 3
            push!(colors, RGBA(0, 0, 1.0, 1.0))
        end
        #push!(colors, RGBA(rand(),rand(),rand(),1.0) )
    end
    gradient = ColorGradient(colors)
    heatmap(x_coords,y_coords,pixels,aspect_ratio=1, color=gradient, leg=false)
end

"""
    plot_cell_2D(geometry::Geometry, view::Box, dim::Int64, cell_id::Int64)

Plots a 2D x-y slice of a geometry, highlighting a specific cell in black.

# Arguments
* `geometry::Geometry`: the geometry we want to plot
* `view::Box`: The view box is an axis aligned box that defines where the picture will be taken, with both min and max z dimensions indicating the single z elevation the slice is taken at.
* `dim::Int64`: The dimension is the number of pixels along the x and y axis to use, which determines the resolution of the picture.
* `cell_id::Int64`: The index of the cell we wish to view
"""
function plot_cell_2D(geometry::Geometry, view::Box, dim::Int64, cell_id::Int64)
    delta_x = (view.upper_right.x - view.lower_left.x) / (dim)
    delta_y = (view.upper_right.y - view.lower_left.y) / (dim)

    x_coords = collect(view.lower_left.x + delta_x/2.0:delta_x:view.upper_right.x - delta_x/2.0)
    y_coords = collect(view.lower_left.y + delta_y/2.0:delta_y:view.upper_right.y - delta_y/2.0)

    pixels = Array{Int64, 2}(dim, dim)

    for i=1:dim
        for j=1:dim
            pixels[i,j] = find_cell_id(Coord(x_coords[i], y_coords[j], view.lower_left.z), geometry)
            if pixels[i, j] == cell_id
                pixels[i, j] = 0
            else
                pixels[i, j] = 1
            end
        end
    end
    pixels = rotl90(pixels)
    colors = Array{RGBA}(0)
    push!(colors, RGBA(0.0, 0.0, 0.0, 1.0))
    push!(colors, RGBA(1.0, 1.0, 1.0, 1.0))
    gradient = ColorGradient(colors)
    heatmap(x_coords,y_coords,pixels,aspect_ratio=1, color=gradient, leg=false)
end

end # module


text:
promote_graph(graph::AbstractMatrix, nf::AbstractMatrix) = graph
promote_graph(graph::AbstractMatrix, nf::Fill{T,S,Axes}) where {T,S,Axes} = graph

ext: jl


text:
module BenchFull

using BenchmarkTools

module SeqFull
using Base.Cartesian: @nexprs

function f(height = 5)
    @nexprs 32 i -> g(i, height)
end

function g(i, height)
    height <= 0 && return
    i < 32 && return
    @nexprs 32 i -> g(i, height - 1)
end

end  # module SeqFull

module BaseFull
using Base.Cartesian: @nexprs

function f(height = 5)
    @sync begin
        @nexprs 32 i -> Threads.@spawn g(i, height)
    end
end

function g(i, height)
    height <= 0 && return
    i < 32 && return
    @sync begin
        @nexprs 32 i -> Threads.@spawn g(i, height - 1)
    end
end

end  # module BaseFull

module TapirFull
using Base.Cartesian: @nexprs
using Base.Experimental: Tapir

function f(height = 5)
    Tapir.@sync begin
        @nexprs 32 i -> Tapir.@spawn g(i, height)
    end
end

function g(i, height)
    height <= 0 && return
    i < 32 && return
    Tapir.@sync begin
        @nexprs 32 i -> Tapir.@spawn g(i, height - 1)
    end
end

end  # module TapirFull

module WSFull
using Base.Cartesian: @nexprs
using Base.Experimental: Tapir
using TapirSchedulers

function f(height = 5)
    Tapir.@sync WorkStealingTaskGroup() begin
        @nexprs 32 i -> Tapir.@spawn g(i, height)
    end
end

function g(i, height)
    height <= 0 && return
    i < 32 && return
    Tapir.@sync WorkStealingTaskGroup() begin
        @nexprs 32 i -> Tapir.@spawn g(i, height - 1)
    end
end

end  # module WSFull

function setup(heights = [5])
    suite = BenchmarkGroup()
    for height in heights
        s0 = suite["height=$height"] = BenchmarkGroup()
        s0["seq"] = @benchmarkable SeqFull.f($height)
        s0["base"] = @benchmarkable BaseFull.f($height)
        s0["tapir"] = @benchmarkable TapirFull.f($height)
        s0["ws"] = @benchmarkable WSFull.f($height)
    end
    return suite
end

function clear() end

end  # module BenchFull


text:
module Neighborhoods

using ..PCG.Types
using ..PCG.Topologies

export Neighborhood

# TODO: make a template argument for all functions?
#       since distance can not be more complex value, than scalar
const Distance = Float64

struct NeighborhoodRange
    min_distance::Distance
    max_distance::Distance
end

function NeighborhoodRange()
    return new(one(Distance), one(Distance))
end

function NeighborhoodRange(distance::Distance)
    return new(distance, distance)
end

struct NeighborhoodTemplate{I<:TopologyIndex}
    range::NeighborhoodRange
    template::Vector{I}
end

struct Neighborhood{I<:TopologyIndex}
    index_type::Type{I}
    templates::Vector{NeighborhoodTemplate{I}}
    distance::Function

    function Neighborhood(topology::Topology, distance::Function)
        I = index_type(topology)
        return new{I}(I, NeighborhoodTemplate{I}[], distance)
    end
end

const Distances{I<:TopologyIndex} = Dict{I, Distance}

function get_min_distance(distances::Distances)
    min_index = nothing
    min_distance = Inf

    for (index, distance) in distances
        if distance < min_distance
            min_index = index
            min_distance = distance
        end
    end

    return (min_index, min_distance)
end

# TODO: rewrite to faster algorithm
function neighbors_before_distance(index::I, distance::Function, max_distance::Distance) where {I<:TopologyIndex}
    processed = Distances{I}()

    queue = Distances{I}()

    # theoretically, "zero" distance can be not equal to «number 0»
    queue[index] = distance(index, index)

    while !isempty(queue)

        min_index, min_distance = get_min_distance(queue)

        pop!(queue, min_index)

        processed[min_index] = min_distance

        for neighbor in neighborsof(min_index)

            if haskey(processed, neighbor) || haskey(queue, neighbor)
                continue
            end

            distance_to = distance(index, neighbor)

            if max_distance < distance_to
                continue
            end

            queue[neighbor] = distance_to
        end

    end

    return processed
end

function construct_template(neighborhood::Neighborhood{I}, range::NeighborhoodRange)::NeighborhoodTemplate{I} where {I<:TopologyIndex}

    neighbors = neighbors_before_distance(zero(neighborhood.index_type),
                                          neighborhood.distance,
                                          range.max_distance)

    area = Vector{neighborhood.index_type}()
    sizehint!(area, length(neighbors))

    for (index, distance_to) in neighbors
        if distance_to < range.min_distance
            continue
        end

        if range.max_distance < distance_to
            continue
        end

        push!(area, index)
    end

    return NeighborhoodTemplate{neighborhood.index_type}(range, area)
end

function (neighborhood::Neighborhood{I})(min_distance, max_distance)::NeighborhoodTemplate{I} where {I<:TopologyIndex}
    range = NeighborhoodRange(convert(Distance, min_distance),
                              convert(Distance, max_distance))

    for template in neighborhood.templates
        if template.range == range
            return template
        end
    end

    template = construct_template(neighborhood, range)

    push!(neighborhood.templates, template)

    return template
end

function (neighborhood::Neighborhood)(min_distance)
    return neighborhood(min_distance, min_distance)
end

function (neighborhood::Neighborhood)()
    return neighborhood(one(Distance), one(Distance))
end

function (neighborhood::NeighborhoodTemplate)(element::E) where E

    if !isenabled(element)
        return reserve_area!(element, 0)
    end

    template = neighborhood.template

    elements = reserve_area!(element, length(template))

    # TODO: optimize?
    fill!(elements, disable(element))

    for (i, delta) in enumerate(template)
        coordinates = element.topology_index + delta

        # TODO: do smth with that
        if is_valid(element.universe.topology, coordinates)
            elements[i] = construct_current_element(E, element.universe, coordinates)
        end
    end

    return elements
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


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

# %% [markdown]
# https://docs.julialang.org/en/v1/manual/running-external-programs/

# %%
pr(x) = print(read(x, String))

# %%
pr(`ls -l`)

# %%
ENV["HOGE"] = "hoge"

# %%
pr(pipeline(`env`, `grep '^H'`))

# %%

max_stars_repo_path: 0007/run commands.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/run commands.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/run commands.jl

max_forks_repo_name: genkuroki/public


text:
using Clang

# LIBXPRESS_HEADERS are those headers to be wrapped.
const LIBXPRESS_INCLUDE = joinpath(ENV["XPRESSDIR"], "include") |> expanduser |> normpath
const LIBXPRESS_HEADERS = [
                           joinpath(LIBXPRESS_INCLUDE, header) for header in [
                                                                               # "bindrv.h",
                                                                               # "mmhttp.h",
                                                                               # "mmquad.h",
                                                                               # "mmsystem.h",
                                                                               # "mmxprs.h",
                                                                               # "mosel_dso.h",
                                                                               # "mosel_mc.h",
                                                                               # "mosel_ni.h",
                                                                               # "mosel_nifct.h",
                                                                               # "mosel_rt.h",
                                                                               # "xprb.h",
                                                                               # "xprb_cpp.h",
                                                                               # "xprd.h",
                                                                               # "xprm_mc.h",
                                                                               # "xprm_ni.h",
                                                                               # "xprm_rt.h",
                                                                               # "xprnls.h",
                                                                               "xprs.h",
                                                                               # "xprs_mse_defaulthandler.h",
                                                                               # "XPRSdefaultMipSolEnumHandler.java",
                                                                               # "xslp.h",
                                                                            ]
                          ]

wc = init(; headers = LIBXPRESS_HEADERS,
            output_file = joinpath(@__DIR__, "../src/lib.jl"),
            common_file = joinpath(@__DIR__, "../src/common.jl"),
            clang_includes = vcat(LIBXPRESS_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(LIBXPRESS_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libxprs",
            clang_diagnostics = true,
            )

run(wc)


text:
export DeviceInfo, defaultInputDevice, defaultOutputDevice, deviceCount

function DeviceInfo(devid::Integer)
    devinfo = ccall((:Pa_GetDeviceInfo, portaudio), Ptr{DeviceInfo}, (Int32,), devid)
    if devinfo == C_NULL
        error("Could not get device info")
    end
    return unsafe_load(devinfo)
end

function defaultInputDevice()
    devid = ccall((:Pa_GetDefaultInputDevice, portaudio), Int32, ())
    if devid < 0
        error("No default input device")
    end
    return devid
end

function defaultOutputDevice()
    devid = ccall((:Pa_GetDefaultOutputDevice, portaudio), Int32, ())
    if devid < 0
        error("No default output device")
    end
    return devid
end

function deviceCount()
    devcount = ccall((:Pa_GetDeviceCount, portaudio), Int32, ())
    check_pa_error(devcount)
    return devcount
end


text:
import GL
import GLFW
importall BSP
importall ImmutableArrays
importall Input
importall FileSystem
importall Player

const near = 4
const far = 16384
projMatrix = eye(Matrix4x4{Float32})
fov = 60.0

function translationMatrix(pos::GL.Vec3)
	x = pos.e1
	y = pos.e2
	z = pos.e3
	Matrix4x4{Float32}(
	Vector4{Float32}(1, 0, 0, 0),
	Vector4{Float32}(0, 1, 0, 0),
	Vector4{Float32}(0, 0, 1, 0),
	Vector4{Float32}(x, y, z, 1))
end

modelMatrix = translationMatrix(GL.Vec3(0, 0, 0))

function updateProjMatrix(width::Cint, height::Cint)
	GL.Viewport(width, height)

	aspect = width / height
	ymax = near * tan(fov * pi / 360)
	ymin = -ymax
	xmin = ymin * aspect
	xmax = ymax * aspect

	x = (2*near) / (xmax-xmin)
	y = (2*near) / (ymax-ymin)
	a = (xmax+xmin) / (xmax-xmin)
	b = (ymax+ymin) / (ymax-ymin)
	c = -(far+near) / (far-near)
	d = -(2*far*near) / (far-near)

	global projMatrix = Matrix4x4{Float32}(
	Vector4{Float32}(x, 0, 0, 0),
	Vector4{Float32}(0, y, 0, 0),
	Vector4{Float32}(a, 0, c,-1),
	Vector4{Float32}(0, 0, d, 0))

	return
end

GLFW.Init()
GLFW.OpenWindowHint(GLFW.OPENGL_VERSION_MAJOR, 3)
GLFW.OpenWindowHint(GLFW.OPENGL_VERSION_MINOR, 2)
GLFW.OpenWindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.OpenWindowHint(GLFW.OPENGL_FORWARD_COMPAT, 1)
GLFW.OpenWindow(0, 0, 8, 8, 8, 0, 24, 0, GLFW.WINDOW)
GLFW.SetWindowTitle("Quake 2.jl")
GLFW.SetWindowSizeCallback(updateProjMatrix)
GLFW.SwapInterval(0)

println("Red bits:     ", GLFW.GetWindowParam(GLFW.RED_BITS))
println("Green bits:   ", GLFW.GetWindowParam(GLFW.GREEN_BITS))
println("Blue bits:    ", GLFW.GetWindowParam(GLFW.BLUE_BITS))
println("Alpha bits:   ", GLFW.GetWindowParam(GLFW.ALPHA_BITS))
println("Depth bits:   ", GLFW.GetWindowParam(GLFW.DEPTH_BITS))
println("Stencil bits: ", GLFW.GetWindowParam(GLFW.STENCIL_BITS))

GL.Enable(GL.CULL_FACE)
GL.Enable(GL.DEPTH_TEST)

println("Vendor:   ", GL.GetString(GL.VENDOR))
println("Renderer: ", GL.GetString(GL.RENDERER))
println("Version:  ", GL.GetString(GL.VERSION))
println("GLSL:     ", GL.GetString(GL.SHADING_LANGUAGE_VERSION))

bspFile = qopen(string("maps/", ARGS[1], ".bsp"))
bsp = read(bspFile, Bsp)
close(bspFile)

maxLights = bsp.max_lights + 1

const vertex_shader_src = string("
#version 150

uniform mat4 ModelMatrix;
uniform mat4 ViewMatrix;
uniform mat4 ProjMatrix;

uniform vec4 TexU;
uniform vec4 TexV;
uniform uint TexW;
uniform uint TexH;

uniform vec3 CameraPosition;

in vec3 VertexPosition;

out vec3 FragPosition;
out vec2 TexCoords;
out vec3 ViewDir;

void main()
{
	FragPosition = VertexPosition;
	ViewDir = normalize(VertexPosition - CameraPosition);
	vec4 pos = vec4(VertexPosition, 1.0);
	TexCoords = vec2(dot(TexU, pos) / TexW, dot(TexV, pos) / TexH);
	gl_Position = ProjMatrix * ViewMatrix * ModelMatrix * pos;
}
")

const fragment_shader_src = string("
#version 150

struct light_t
{
	vec3 Position;
	vec3 Color;
	float Power;
};

uniform int DrawMode;

uniform vec3 FaceNormal;

uniform bool DiffuseLighting;
uniform bool SpecularLighting;

uniform vec3 AmbientLight;
uniform int NumLights;
uniform light_t Light[", maxLights, "];

uniform bool DiffuseMapping;
uniform sampler2D DiffuseMap;
uniform sampler2D NormalMap;
uniform uint SurfFlags;

uniform vec4 TexU;
uniform vec4 TexV;

in vec3 FragPosition;
in vec2 TexCoords;
in vec3 ViewDir;

out vec4 FragColor;

void main()
{
	vec3 tang = normalize(TexU.xyz);
	vec3 bitang = normalize(TexV.xyz);
	mat3 toTangentSpace = mat3(
		tang.x, bitang.x, FaceNormal.x,
		tang.y, bitang.y, FaceNormal.y,
		tang.z, bitang.z, FaceNormal.z );
	vec4 normalmap = 2 * texture(NormalMap, TexCoords) - vec4(1.0);
	vec2 offset = -normalmap.w * 0.03 * (toTangentSpace * ViewDir).xy;
	normalmap = 2 * texture(NormalMap, TexCoords+offset) - vec4(1.0);
	vec3 normal = normalmap.xyz;
	vec3 LightColor = AmbientLight;
	vec3 camReflectDir = reflect(toTangentSpace * ViewDir, normal);
	for (int i = 0; i < NumLights; i++) {
		vec3 lightDir = toTangentSpace * normalize(Light[i].Position - FragPosition);
		if (lightDir.z < 0) {
			continue;
		}
		float dirMod = dot(normal, lightDir); // -1 to 1
		dirMod = max(dirMod, 0);
		float lightDist = distance(Light[i].Position, FragPosition);
		float distMod = (Light[i].Power - lightDist) / Light[i].Power;
		distMod = pow(clamp(distMod, 0.0, 1.0), 2);
		if (DiffuseLighting) {
			LightColor += 1.5 * Light[i].Color * dirMod * distMod;
		}
		if (SpecularLighting) {
			LightColor += 1.5 * Light[i].Color * pow(max(dot(lightDir, camReflectDir), 0.0), 10) * distMod * 0.8;
		}
	}
	if (uint(SurfFlags & uint(1)) != uint(0)) { // FIXME: stupid casts
		LightColor = vec3(1.0);
	}
	LightColor = min(LightColor, vec3(1.0, 1.0, 1.0));

	if (DrawMode == 1) {
		FragColor.rgb = LightColor;
	} else if (DrawMode == 2) {
		FragColor.rgb = texture(DiffuseMap, TexCoords).rgb;
	} else if (DrawMode == 3) {
		FragColor.rgb = vec3(texture(NormalMap, TexCoords).a);
	} else if (DrawMode == 4) {
		FragColor.rgb = texture(NormalMap, TexCoords).rgb;
	} else {
		FragColor.rgb = texture(DiffuseMap, TexCoords+offset).rgb;
		FragColor.rgb *= LightColor;
	}
	FragColor.a = 1.0;
}
")

prog = GL.Program([GL.VERTEX_SHADER   => vertex_shader_src,
                   GL.FRAGMENT_SHADER => fragment_shader_src])

uModel = GL.Uniform(prog, "ModelMatrix")
uView = GL.Uniform(prog, "ViewMatrix")
uProj = GL.Uniform(prog, "ProjMatrix")

uCamPos = GL.Uniform(prog, "CameraPosition")

uNormal = GL.Uniform(prog, "FaceNormal")

uSurfFlags = GL.Uniform(prog, "SurfFlags")
uTexU = GL.Uniform(prog, "TexU")
uTexV = GL.Uniform(prog, "TexV")
uTexW = GL.Uniform(prog, "TexW")
uTexH = GL.Uniform(prog, "TexH")

aPosition = GL.GetAttribLocation(prog, "VertexPosition")

vao = GL.GenVertexArray()
GL.BindVertexArray(vao)

positionBuf = GL.GenBuffer()
GL.BindBuffer(GL.ARRAY_BUFFER, positionBuf)
GL.BufferData(GL.ARRAY_BUFFER, bsp.vertices, GL.STATIC_DRAW)
GL.EnableVertexAttribArray(aPosition)
GL.VertexAttribPointer(aPosition, 3, GL.FLOAT, false, 0, 0)
GL.BindBuffer(GL.ARRAY_BUFFER, 0)

GL.BindVertexArray(0)

frames = 0
tic()

function rotationMatrix{T<:Real}(eyeDir::AbstractVector{T}, upDir::AbstractVector{T})
	rightDir = cross(eyeDir, upDir)
	rightDir /= norm(rightDir)
	upDir = cross(rightDir, eyeDir)

	Xx = rightDir.e1
	Xy = rightDir.e2
	Xz = rightDir.e3

	Yx = upDir.e1
	Yy = upDir.e2
	Yz = upDir.e3

	Zx = -eyeDir.e1
	Zy = -eyeDir.e2
	Zz = -eyeDir.e3

	Matrix4x4{Float32}(
	Vector4{Float32}(Xx, Yx, Zx,  0),
	Vector4{Float32}(Xy, Yy, Zy,  0),
	Vector4{Float32}(Xz, Yz, Zz,  0),
	Vector4{Float32}( 0,  0,  0,  1))
end

uDrawMode = GL.Uniform(prog, "DrawMode")
uAmbient = GL.Uniform(prog, "AmbientLight")
uDiffuse = GL.Uniform(prog, "DiffuseLighting")
uSpecular = GL.Uniform(prog, "SpecularLighting")
numLightsUniform = GL.Uniform(prog, "NumLights")
lightUniforms = Array(GL.Uniform, 0)
for i = 0:maxLights-1
	light = string("Light[", i, "].")
	push!(lightUniforms, GL.Uniform(prog, string(light, "Position")))
	push!(lightUniforms, GL.Uniform(prog, string(light, "Color")))
	push!(lightUniforms, GL.Uniform(prog, string(light, "Power")))
end

GL.UseProgram(prog)
write(GL.Uniform(prog, "DiffuseMap"), int32(0))
write(GL.Uniform(prog, "NormalMap"), int32(1))

light1 = Player.State()
Player.movedir!(light1, GL.Vec3(1,0,0), true)
light1_pow = float32(0)

draw_mode = 0
ambient_lighting_on = true
diffuse_lighting_on = true
specular_lighting_on = true
wireframe_only = false

toggle_ambient_light() = global ambient_lighting_on = !ambient_lighting_on
toggle_diffuse_light() = global diffuse_lighting_on = !diffuse_lighting_on
toggle_specular_light() = global specular_lighting_on = !specular_lighting_on
toggle_wireframe() = global wireframe_only = !wireframe_only

function fire(apply::Bool)
	if apply
		light1.position = Player.self.position
		light1.orientation = Player.self.orientation
		light1.speed = 0
		global light1_pow = float32(300)
	else
		light1.speed = 300
	end
end

bind(GLFW.MOUSE_BUTTON_LEFT, fire)
bind(GLFW.MOUSE_BUTTON_RIGHT, in_grab)
bind(GLFW.KEY_ESC, in_release)

# WASD in Dvorak
bind(',', forward)
bind('A', moveleft)
bind('O', back)
bind('E', moveright)
bind(' ', moveup)
bind(GLFW.KEY_LCTRL, movedown)
bind(GLFW.KEY_LSHIFT, speed)

bind(GLFW.KEY_UP, lookup)
bind(GLFW.KEY_DOWN, lookdown)
bind(GLFW.KEY_LEFT, left)
bind(GLFW.KEY_RIGHT, right)

draw_mode_0() = global draw_mode = 0
draw_mode_1() = global draw_mode = 1
draw_mode_2() = global draw_mode = 2
draw_mode_3() = global draw_mode = 3
draw_mode_4() = global draw_mode = 4

bind('0', draw_mode_0)
bind('1', draw_mode_1)
bind('2', draw_mode_2)
bind('3', draw_mode_3)
bind('4', draw_mode_4)

bind('6', toggle_ambient_light)
bind('7', toggle_diffuse_light)
bind('8', toggle_specular_light)

bind('[', toggle_wireframe)

GLFW.SetKeyCallback(Input.event)
GLFW.SetMouseButtonCallback(Input.event)
GLFW.SetMousePosCallback(Input.look_event)
GLFW.SetMouseWheelCallback(Input.wheel_event)

while GLFW.GetWindowParam(GLFW.OPENED)

	Player.move!(Player.self)
	Player.move!(light1)
	eyedir, updir, rightdir = sphere2cartesian(Player.self.orientation)
	rotMat = rotationMatrix(eyedir, updir)

	GL.Clear(GL.COLOR_BUFFER_BIT | GL.DEPTH_BUFFER_BIT)

	write(uModel, modelMatrix)
	transMat = translationMatrix(-Player.self.position)
	viewMat = rotMat * transMat
	write(uView, viewMat)
	write(uProj, projMatrix)

	write(uCamPos, Player.self.position)

	write(uDrawMode, int32(draw_mode))
	if wireframe_only
		write(uAmbient, GL.Vec3(0.2, 0.2, 0.2))
	elseif ambient_lighting_on
		write(uAmbient, GL.Vec3(0.075, 0.075, 0.075))
	else
		write(uAmbient, GL.Vec3(0, 0, 0))
	end
	write(uDiffuse, diffuse_lighting_on)
	write(uSpecular, specular_lighting_on)

	write(lightUniforms[1], light1.position)
	write(lightUniforms[2], GL.Vec3(1.0, 1.0, 0.0))
	write(lightUniforms[3], light1_pow)

	GL.BindVertexArray(vao)

	for face = search(bsp, Player.self.position).faces
		write(uNormal, face.normal)

		write(uSurfFlags, face.texture.flags)
		write(uTexU, face.u_axis)
		write(uTexV, face.v_axis)
		write(uTexW, face.texture.width)
		write(uTexH, face.texture.height)

		write(numLightsUniform, int32(length(face.lights)+1))
		i = 4
		for light = face.lights
			write(lightUniforms[i], light.origin); i += 1
			write(lightUniforms[i], light.color); i += 1
			write(lightUniforms[i], light.power); i += 1
		end

		GL.ActiveTexture(GL.TEXTURE0)
		GL.BindTexture(GL.TEXTURE_2D, face.texture.diffuse)
		GL.ActiveTexture(GL.TEXTURE1)
		GL.BindTexture(GL.TEXTURE_2D, face.texture.normal)

		draw = GL.TRIANGLES
		if wireframe_only
			draw = GL.LINE_LOOP
		end
		GL.BindBuffer(GL.ELEMENT_ARRAY_BUFFER, face.ibo);
		GL.DrawElements(draw, length(face.indices), GL.UNSIGNED_SHORT, 0)
		GL.BindBuffer(GL.ELEMENT_ARRAY_BUFFER, 0);

		GL.GetError()
	end

	GL.BindVertexArray(0)

	GLFW.SwapBuffers()
	frames += 1
end

seconds = toc()
println(frames / seconds, " FPS")

GLFW.Terminate()


text:
# good projects are tested via project_examples.jl
const MissingMilestoneException = ProjectManagement.MissingMilestoneException
const BadMilestoneDurationException = ProjectManagement.BadMilestoneDurationException
const PreStartTasksException = ProjectManagement.PreStartTasksException
const PostFinishTasksException = ProjectManagement.PostFinishTasksException
const CyclicPathException = ProjectManagement.CyclicPathException
const DisconnectedPathException = ProjectManagement.DisconnectedPathException
const UnlinkedTasksException = ProjectManagement.UnlinkedTasksException
const UnspecifiedTasksException = ProjectManagement.UnspecifiedTasksException

@testset "bad projects" begin
    @test_throws MissingMilestoneException Project(
        (start=0, a=1),
        [:start => :a,]
    )

    @test_throws MissingMilestoneException Project(
        (a=1, finish=0,),
        [:a => :finish,]
    )

    @test_throws BadMilestoneDurationException Project(
        (start=1, finish=0,),
        [:start => :finish,]
    )

    @test_throws BadMilestoneDurationException Project(
        (start=0, finish=(0,0),),
        [:start => :finish,]
    )

    @test_throws PreStartTasksException Project(
        (start=0, a=1, finish=0,),
        [:a => :start, :start => :finish,]
    )

    @test_throws PostFinishTasksException Project(
        (start=0, a=1, finish=0,),
        [:start => :finish, :finish=>:a]
    )

    @test_throws CyclicPathException Project(
        (start=0, a=1, b=1, finish=0,),
        [:start => :a, :a => :b, :b  => :a, :a => :finish,]
    )

    @test_throws DisconnectedPathException Project(
        (start=0, a=1, b=1, finish=0,),
        [:start => :a, :a => :b, :a => :finish,]
    )

    @test_throws UnspecifiedTasksException Project(
        (start=0, a=1, finish=0,),
        [:start => :a, :a => :b, :b => :finish, :a => :finish,]
    )

    @test_throws UnlinkedTasksException Project(
        (start=0, a=1, b=1, finish=0,),
        [:start => :a, :a => :finish,]
    )
end

@testset "Error messages" begin
    msg = sprint(showerror, BadMilestoneDurationException(:start, 2))
    @test occursin("start", msg)
    @test occursin("2", msg)
end


text:
# Activating Julia environment
using Pkg
Pkg.activate(".")

# Packages
using JuMP, Gurobi, Random, Distributions, LinearAlgebra, DataFrames, Plots

n = 10 # Number of facilities
m = 50 # Number of customers

# Generating random data (please don't change the seeds.)
facilities = 0.6.*rand(MersenneTwister(5), n,2) .+ 0.2;
customers = rand(MersenneTwister(2), m, 2);
c = [LinearAlgebra.norm(customers[i, :] .- facilities[j, :])[1] for j=1:n, i=1:m];
f = rand(MersenneTwister(3), n)*1 .+ 5;
s = rand(MersenneTwister(4), n)*2 .+ 15;
d = rand(MersenneTwister(5), m)*0.5 .+ 0.75

# P matrix
R_D = 0.25
P = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P = (P .>= 0.2*exp(-1/R_D .* R_D)) .* P

""" Plots the solution of the facility location model.
Blue circles are active facilities with different capacities.
Orange plus signs are other potential facility locations.
Rays describe connections between facilities and demand nodes.
"""
function plot_solution(model, x, y, cost = nothing)
    plt = scatter(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s .* value.(x))
    scatter!(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s, markershape = :+)
    for i=1:n
        for j=1:m
            if value(y[i,j]) >= 1e-10
                plot!([customers[j, 1], facilities[i,1]], [customers[j,2], facilities[i,2]], linewidth = value(y[i,j]), legend=false)
            end
        end
    end
    if cost == nothing
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
                title = "Total cost: $(round(objective_value(model), sigdigits=5))")
    else
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
        title = "Total cost: $(round(cost,sigdigits=5))")
    end
    println("Facility cost: $(value(sum(f[j] * x[j] for j = 1:n)))")
    println("Transportation cost: $(value(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)))")
    return plt
end


text:
using ConstructiveGeometry
using Test
using StaticArrays
using LinearAlgebra

G = ConstructiveGeometry

@testset "2d" begin#««
nva(s::G.PolygonXor) = (length.(G.paths(s)), G.Shapes.area(s.poly))
r1 = G.square(2.,1)
r2 = G.square(1.,2)
r3 = G.square(3.,3)
r4 = [4,4] + G.square(1,1)
@test nva(mesh(union(r1,r2))) == ([6], 3)
@test nva(mesh(union(r1,r3))) == ([4], 9)
@test nva(mesh(union(r1,r4))) == ([4,4], 3)
@test nva(mesh(intersect(r1,r2))) == ([4], 1)
@test nva(mesh(intersect(r2,r1))) == ([4], 1)
@test nva(mesh(intersect(r1,r3))) == ([4], 2)
@test nva(mesh(intersect(r3,r1))) == ([4], 2)
@test nva(mesh(intersect(r1,r4))) == ([], 0)
@test nva(mesh(setdiff(r1,r2))) == ([4], 1)
@test nva(mesh(setdiff(r3,r1))) == ([6], 7)
end#»»
# @testset "Extrusion" begin #««1
# # using ConstructiveGeometry: nvertices, nfaces
# # C = vertices(circle(3.),(precision=.01,accuracy=1,symmetry=1))
# # c = [Point(20*cos(i),20*sin(i)) for i in 0:.1:π]; c=[c;[Point(0.,-1.)]]
# # @test (path_extrude(c, C)) != 0
# #
# d=setdiff(square(15.), translate([1,1])*square(8.))
# # e=linear_extrude(10.)*d
# # m=mesh(e)
# # @test nvertices(m) == 16
# # @test nfaces(m) == 32
# end
@testset "Elementary operations on objects" begin #««1
s = square(1)
@testset "Primitives" begin #««2
@test s == square([1,1])
@test s == square(1,1)
end
@testset "Operations" begin #««2
@test union(s, union()) === s
@test length(G.children(union(s, union(s)))) == 2
end
@testset "Transforms" begin #««2
@test 2s == scale(2, s)
@test scale(2)*s == scale(2, s)
@test scale(2)*[s] == scale(2, s)
@test color("red")*s == color("red", s)
end
end
# @testset "Convex hull" begin #««1
# using ConstructiveGeometry: convex_hull, convex_hull_list
# P(x...) = Point(Float64.(x)...)
# CH = convex_hull([P(0,0,0),P(0,0,10),P(10,0,0),P(0,10,0),P(1,1,1),P(1,0,0),])
# @test Set(CH[1]) == Set([P(0,0,0),P(0,0,10),P(10,0,0),P(0,10,0),])
# @test length(CH[2]) ==4
# @test convex_hull_list([ P(0., 0), P(1., 1), P(-1., 1), P(.2, 2.), P(0, .8), P(-.2, .8), ]) == [1,2,4,3]
# # @test convex_hull_list(ConstructiveGeometry.rows(SA[0. 0;1. 1;-1. 1;.2 2.;0 .8;-.2 .8])) == [1,2,4,3]
# @test convex_hull_list([
# 	P(-2.627798062316817, 1.075268817204301),
# 	P(-0.5030257403564974, -1.720430107526882),
# 	P(0.7927283156659947, 2.7956989247311825),
# 	P(0.0, 2.396978520135108),
# 	P(0.0, 0.03278003249806397),
# 	]) == [2,5,3,1]
# @test convex_hull_list([
# 	P(-2.150537634408602, 1.3494700327417308),
# 	P(-0.4301075268817205, -2.097580910437773),
# 	P(2.3655913978494625, 0.04739817471346019),
# 	P(0.0, 0.0),
# 	P(2.3038140933536018, 0.0),
# 	P(0.0, 0.7294358146330306),
# 	])== [2,5,3,6,1]
# end
# @testset "Triangulation" begin #««1
# using ConstructiveGeometry: LibTriangle, triangulate, triangulate_loop, identify_polygons, PolygonXor
# @test length(triangulate_loop(Point{2,Float64}.([(0,0),(1,0),(1,1),(0,1)]))) == 2
# square_with_hole = PolygonXor([[5.0, 5.0], [0.0, 5.0], [0.0, 0.0], [5.0, 0.0]],[[1.0, 1.0], [1.0, 3.0], [3.0, 3.0], [3.0, 1.0]])
# @test identify_polygons(square_with_hole) == [1, -1]
# @test length(triangulate(square_with_hole)) == 8
# @test CG.ladder_triangles(5,4,0,10) ==
# 	[(0,10,1),(1,10,11),(1,11,2),(2,11,12),(2,12,3),(3,12,13),(3,13,4)]
# end
# @testset "Basic geometry" begin#««1
# v = [[-1,0],[0,-1],[1,0],[0,1]]
# m = [CG.circular_sign(i,j) for i in v, j in v]
# @test sign.(m) == [0 1 1 1; -1 0 1 1; -1 -1 0 1; -1 -1 -1 0]
# end
# @testset "Triangle intersection" begin#««1
# TI=CG.TriangleIntersections
# Pos=TI.Constants
# pts(m) = ([Float64.(m[i,:]) for i in 1:size(m,1)]...,)
# function check_inter1(s, t, n)#««
# 	it = inter(s,t)
# 	return length(it)==n &&
# 		all(u[1] ==  pttype(p, s) && u[2] == pttype(p,t) for (p,u) in it)
# end#»»
# function check_inter(s1,t1,n)#««
# 	perm = (([1],),
# 		([1,2],[2,1]),
# 		([1,2,3],[2,3,1],[3,1,2],[1,3,2],[2,1,3],[3,2,1]))
# 	s = (s1[p] for p in perm[length(s1)])
# 	if length(s1) == 2
# 		t = (t1[p] for p in perm[3][1:3])
# 	else
# 		t = (t1[p] for p in perm[length(t1)])
# 	end
# 	return all(check_inter1(s2,t2,n) for s2 in s, t2 in t)
# end#»»
# function pttype(a,(u,v)::NTuple{2},ε=1e-8)#««
# 	abs(TI.det2(u,v,a)) > ε && return -1
# 	t = dot(a-u,v-u)/dot(v-u,v-u)
# 	t < 0-ε && return Pos.invalid
# 	t > 1+ε && return Pos.invalid
# 	t ≤ 0+ε && return Pos.vertex1
# 	t ≥ 1-ε && return Pos.vertex2
# 	return Pos.edge12
# end#»»
# function pttype(a,(p,q,r)::NTuple{3,<:StaticVector{2}},ε=1e-8)#««
# 	(d1,d2,d3) = (TI.det2(q,r,a),TI.det2(r,p,a),TI.det2(p,q,a))
# 	s(x) = (x > ε) ? '+' : (x<-ε) ? '-' : '0'
# 	m = s(d1)s(d2)s(d3)
# 	@assert '+' ∈ m "bad point type ($m)for $a in ($p, $q, $r)"
# 	m == "+++" && return Pos.interior
# 	m == "0++" && return Pos.edge23
# 	m == "+0+" && return Pos.edge31
# 	m == "++0" && return Pos.edge12
# 	m == "+00" && return Pos.vertex1
# 	m == "0+0" && return Pos.vertex2
# 	m == "00+" && return Pos.vertex3
# 	return Pos.invalid
# end#»»
# @inline inter(s::NTuple{2}, t::NTuple{3}) = TI.inter_segment2_triangle2(s,t)
# @inline inter(s::NTuple{3,<:StaticVector{3}}, t::NTuple{3}) = TI.inter(s,t)
# t2=pts(SA[0 0;10 0;0 10])
#
# # vertices:
# @test check_inter(pts(SA[-1 1;10 0]), t2, 2)
# @test check_inter(pts(SA[-1 1;0 10]), t2, 1)
# @test check_inter(pts(SA[-5 5;0 0]), t2, 1)
# # edges:
# @test check_inter(pts(SA[-1 1;2 0]), t2, 2)
# @test check_inter(pts(SA[-1 1;0 2]), t2, 1)
# @test check_inter(pts(SA[-1 1;7 3]), t2, 2)
# # contact with vertices:
# @test check_inter(pts(SA[-5 5;1 -1]), t2, 1)
# # inside:
# @test check_inter(pts(SA[-1 1;1 1]), t2, 2)
# # outside:
# @test check_inter(pts(SA[-1 1;2 10]), t2, 2)
# @test check_inter(pts(SA[-1 5;2 20]), t2, 1)
# @test check_inter(pts(SA[-1 5;2 30]), t2, 0)
# @test check_inter(pts(SA[-1 5;10 -3]), t2, 2)
# function pttype(a,(p,q,r)::NTuple{3,<:StaticVector{3}},ε=1e-8)#««
# 	normal2 = cross(q-p,r-p)
# 	if abs(dot(normal2, p-a)) > ε
# 		return -1
# 	end
# 	proj = TI.Projector(normal2)
# 	(a1, p1, q1, r1) = proj.((a,p,q,r))
# 	return pttype(a1, (p1,q1,r1), ε)
# end#»»
# t3 = pts(SA[0 0 0;5 0 0;0 5 0])
# function dotest(c, a, b, n)#««
# 	# triangle on plane [x+y=c],
# 	# crossing [z=0] at (3-a,a,0) and (3-b,b,0)
# 	# (edges 3 and 2)
# 	t1 = pts(SA[c-a a 1; c-a a -1; c-2b+a 2b-a  -1])
#
# 	check_inter(t3,t1,n)
# 	check_inter(t1,t3,n)
# end#»»
# @test check_inter(pts(SA[2. 0 0;2 2 0;1 1 1]),
# 		pts(SA[3. 2 0;1 2 0;2 1 1]), 2)
# 	# touch configuration:
# @test dotest(0, -2, -1, 0)
# @test dotest(0, -2, 0,  1)
# @test dotest(0, -2, 1, 1)
# @test dotest(0, 0, 1, 1)
# @test dotest(0, 1, 2, 0)
#
# 	# general position:
# @test dotest(3, -2,-1, 0)
# @test dotest(3, -1, 0, 1)
# @test dotest(3, -1, 2, 2)
# @test dotest(3, -1, 3, 2)
# @test dotest(3, -1, 4, 2)
#
# @test dotest(3, 0, 2, 2)
# @test dotest(3, 0, 3, 2)
# @test dotest(3, 0, 4, 2)
#
# @test dotest(3, 1, 2, 2)
# @test dotest(3, 1, 3, 2)
# @test dotest(3, 1, 4, 2)
#
# @test dotest(3, 3, 4, 1)
# @test dotest(3, 5, 4, 0)
#
# 	# arrow configuration:
# @test check_inter(pts(SA[-1 1 1;-1 1 -1;2 2 0]), t3, 2)
#
# 	# border configuration:
# @test check_inter(pts(SA[0 0 0;5 0 0;0 0 1]), t3, 2)
# @test check_inter(pts(SA[-1 1 0; 2 2 0;0 0 1]), t3, 2)
# @test check_inter(pts(SA[1 0 0;3 0 0;0 0 1]), t3, 2)
#
# 	# coplanar
# @test check_inter(pts(SA[1 1 0;2 1 0;1 2 0]), t3, 3)
# @test check_inter(pts(SA[-1 2 0;2 -1 0;6 6 0]), t3, 6)
# @test check_inter(pts(SA[28. 0 0;12 0 0;20 0 -8]),
# 		pts(SA[25. 0 0;20 0 -5;15 0 0]), 2)
# end
# @testset "Surfaces" begin #««1
# using ConstructiveGeometry: Surface, merge, mesh
# # using ConstructiveGeometry: nvertices, nfaces
# v=[[0,-1],[1,0],[0,1],[-1,0]]
# for j in eachindex(v), i in 1:j-1
# 	@test ConstructiveGeometry.circular_lt(v[i], v[j])
# end
# nvf(s) = (CG.nvertices(s), CG.nfaces(s))
# p=Surface([[-1.,0,0],[1,0,0],[0,1,0],[0,0,1]],
#          [[3,2,1],[4,1,2],[4,3,1],[4,2,3]])
# @test nvf(mesh(2p \ p)) == (8, 12)
# # @test connected_components([:a, :b, :c, :d, :e], [[1,2],[1,3],[4,5]]) ==
# # 	[([:a, :b, :c], [[1,2],[1,3]]),
# # 	 ([:d, :e], [[1,2]]) ]
# function pyramid(t=[0,0,0], n=0)
# 	points = ([ t, t+[2.,0,0], t+[2,2,0], t+[0,2,0], t+[1,1,1]])
#   faces = [[4,3,2],[2,1,4],[1,2,5],[2,3,5],[3,4,5],[4,1,5]]
# 	faces1 = [ f .+ n for f in faces ]
# 	return surface(points, faces1)
# end
# p1 = pyramid()
# p2 = pyramid([1,0,0])
# u12 = mesh(p1 ∪ p2)
# i12 = mesh(p1 ∩ p2)
# d12 = mesh(p1 \ p2)
# @test nvf(u12) == (14, 24)
# @test nvf(i12) == (8, 12)
# @test nvf(d12) == (8, 12)
#
# function tetrahedron(;center=[0,0,0], radius=1)
# 	return surface([center+radius*[1,0,-.5], center+radius*[-1,0,-.5],
# 		center+radius*[0,1,.5],center+radius*[0,-1,.5]],
# 		[[1,3,4],[1,2,3],[1,4,2],[2,4,3]])
# end
#
# t1 = tetrahedron()
# t2 = 2t1
# t3 = [3,0,0]+t1
# @test nvf(mesh(t1∪t2)) == (4,4)
# @test nvf(mesh(t1∪t3)) == (8,8)
# @test nvf(mesh(t1∩t2)) == (4,4)
# @test nvf(mesh(t1∩t3)) == (0,0)
# @test nvf(mesh(t1\t1)) == (0,0)
# @test nvf(mesh(t1\t3)) == (4,4)
# @test nvf(mesh(t1\t2)) == (0,0) # t1 ⊂ t2
# @test nvf(mesh(t2\t1)) == (8,8)
#
# end
# @testset "Difference of rotate_extrude()" begin# ««1
# m1=mesh(rotate_extrude(10,[2,0]+circle(1)))
# m2=mesh(rotate_extrude(11,[2,0]+circle(.5)))
# @debug "###### difference ######"
# m3=mesh(m1\m2)
# @test CG.nfaces(m3) == 2*CG.nvertices(m3)
# end
# #»»1
# vim: noet ts=2 fmr=««,»»

# @testset "Spatial sorting" begin#««
# include("../src/SpatialSorting.jl")
# # Bounding box««
# struct BoundingBox{N,T}
# 	min::SVector{N,T}
# 	max::SVector{N,T}
# end
#
# BoundingBox{N}(min::AbstractVector, max::AbstractVector) where{N} =
# 	BoundingBox{N,promote_type(eltype(min), eltype(max))}(min, max)
# BoundingBox(min::StaticVector{N}, max::StaticVector{N}) where{N} =
# 	BoundingBox{N}(min, max)
# function BoundingBox(min::AbstractVector, max::AbstractVector)
# 	length(min) == length(max) || throw(DimensionMismatch(
# 	"min and max must have same length ($(length(min)),$(length(max)))"))
# 	return BoundingBox{length(min)}(min, max)
# end
#
# @inline Base.merge(box1::BoundingBox{N}, box2::BoundingBox{N}) where{N} =
# 	BoundingBox{N}(min.(box1.min, box2.min), max.(box1.max, box2.max))
# @inline Base.:∩(box1::BoundingBox{N}, box2::BoundingBox{N}) where{N} =
# 	BoundingBox{N}(max.(box1.min, box2.min), min.(box1.max, box2.max))
# @inline Base.isempty(box::BoundingBox) = any(box.min .> box.max)
# @inline SpatialSorting.position(box::BoundingBox) = box.min + box.max
#
# function rbb(int1, int2, n)
# 	v = rand(int1, n)
# 	w = rand(int2, n)
# 	return BoundingBox(v, v+w)
# end
# boxes=[ rbb(1:500, 1:10, 3) for _ in 1:3000 ]
#
# function validate(boxes)
# 	int = extrema.(SpatialSorting.intersections(boxes))
# 	for i = 1:length(boxes), j = 1:i-1
# 		a = SpatialSorting.intersects(boxes[i], boxes[j])
# 		b = (j,i) ∈ int
# 		(a ≠ b) && return false
# 	end
# 	return true
# end
#
# @test validate(boxes)
#
# #»»
# end#»»


text:
module ContextTracking

using Dates: now
using DataStructures: Stack
using DocStringExtensions
using ExprTools: combinedef, splitdef

export Context
export context
export @ctx, @memo

include("types.jl")
include("context.jl")
include("singleton.jl")
include("trace.jl")

end # module


text:
include("../mesh/mesh_clear.jl")

@testset "Valid boundry box" begin

end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
using CasesLH, Test

function case_spec_test()
    @testset "CaseSpec" begin
        c1 = CaseSpec(:a);
        @test !has_mods(c1);
        @test !has_xp_mods(c1);
        @test !has_modifier(c1, :a);
        @test isequal(c1, CaseSpec(:a));
        @test isequal(make_case_spec(:a), c1);
        @test isequal(make_string(c1), "[a][Base]");
        @test isequal(case_fn(c1), "a_XP_Base");

        c2 = CaseSpec([:a, :b, :c])
        @test has_mods(c2);
        @test !has_xp_mods(c2);
        @test isequal(get_mods(c2), [Modifier(:b), Modifier(:c)])
        @test has_modifier(c2, :b)
        @test isequal(make_case_spec(c2), c2)
        @test isequal(make_string(c2), "[a_b_c][Base]");

        c3 = CaseSpec([:a, :c, :b]);
        @test isequal(make_case_spec(c2), make_case_spec(c3))
        @test length(make_string(c2)) > length(make_string(c1)) > 0

        c4 = CaseSpec(:a, [:b, :c]);
        @test !has_mods(c4);
        @test has_xp_mods(c4);
        @test isequal(make_string(c4), "[a][b_c]");
        @test isequal(case_fn(c4),  "a_XP_b_c");

        c5 = CaseSpec(:a, [(:b, 1), :c]);
        @test isequal(make_string(c5), "[a][b1_c]");

        arg = ([:a, :b], [(:b, 1), :c]);
        c6 = CaseSpec(arg);
        @test isequal(make_string(c6), "[a_b][b1_c]");
        @test isequal(c6, make_case_spec(arg));

        c7 = CaseSpec(arg);
        c7a = replace_xp_mods(c7, [:d, :e]);
        @test isequal(make_string(c7a), "[a_b][d_e]");

        arg = (:a, [(:b, 1), :c]);
        c6a = CaseSpec(arg);
        @test isequal(make_string(c6a),  make_string(c5));

        c7 = CaseSpec(arg);
        c7b = replace_mods(c7, [:d, :e]);
        @test isequal(make_string(c7b), "[a_d_e][b1_c]");

        # strV = string_vector(c1; prefixStr = "-");
        # @test length(strV) == 1
        # @test strV[1] == "-a"

        # strV = string_vector(c3; prefixStr = "-");
        # @test length(strV) == 3
        # @test isequal(strV,  ["-a", "-b", "-c"])
    end
end

function modify_test()
    @testset "Modify" begin
        xpModV = [:two, :three];
        c1 = CaseSpec(:test, xpModV);
        @test has_xp_mods(c1);
        c2 = remove_xp_mods(c1);
        @test !has_xp_mods(c2);

        c3 = add_modifier(c1, :test2);
        @test isequal(c3, CaseSpec([:test, :test2], xpModV));

        c1 = CaseSpec([:test1, :test2], xpModV);
        c2 = add_modifier(c1, :test11);
        @test isequal(c2, CaseSpec([:test1, :test11, :test2], xpModV));

        c1 = CaseSpec(:test);
        c2 = add_modifier(c1, (:test, 10));
        @test isequal(c2, CaseSpec([:test, (:test, 10)]));
    end
end

function find_main_mod_test()
    @testset "Find main mod" begin
        c1 = CaseSpec([:test, :one]);
        @test isnothing(find_main_mod(c1, :two));
        c2 = add_modifier(c1, :two);
        m = find_main_mod(c2, :two);
        @test isequal(m, Modifier(:two));

        c3 = add_modifier(c1, (:two, 2));
        m = find_main_mod(c3, :two);
        @test isequal(m, Modifier((:two, 2)));
    end
end

@testset "CaseSpec" begin
    case_spec_test();
    modify_test();
    find_main_mod_test();
end

# ----------


text:
"""
    nearneighbor(cmd0::String="", arg1=[]; kwargs...)

Reads arbitrarily located (x,y,z[,w]) triples [quadruplets] and uses a nearest neighbor algorithm
to assign an average value to each node that have one or more points within a radius centered on the node.
The average value is computed as a weighted mean of the nearest point from each sector inside the search
radius. The weighting function used is w(r) = 1 / (1 + d ^ 2), where d = 3 * r / search_radius and r is
distance from the node. This weight is modulated by the weights of the observation points [if supplied].

Full option list at [`nearneighbor`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html)

Parameters
----------

- **I** : **inc** : -- Str or Number --

	*x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#i)
- **N** : **sectors** : -- Number or Str --

    The circular area centered on each node is divided into `sectors` sectors.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#n)
- $(GMT.opt_R)
- **S** : **search_radius** : -- Number --

    Sets the search_radius that determines which data points are considered close to a node.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#s)

- **E** : **empty** : -- Bool or [] --

    Set the value assigned to empty nodes when G is set [NaN].
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#e)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = nearneighbor(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#g)
- $(GMT.opt_V)
- **W** : **weights** : -- Bool or [] --

    Input data have a 4th column containing observation point weights.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/nearneighbor.html#w)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_n)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function nearneighbor(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("nearneighbor", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd, opt_R = parse_R("", d)
	cmd, opt_J = parse_J(cmd, d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_n(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'E', d, [:E :empty])
    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'N', d, [:N :ids])
	cmd = add_opt(cmd, 'S', d, [:S :search_radius])
	cmd = add_opt(cmd, 'W', d, [:Z :weights])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "nearneighbor", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
nearneighbor(arg1=[], cmd0::String=""; kw...) = nearneighbor(cmd0, arg1; kw...)


text:
export BrickletAccelerometerAcceleration
struct BrickletAccelerometerAcceleration
    x::Integer
    y::Integer
    z::Integer
end

export BrickletAccelerometerAccelerationCallbackThreshold
struct BrickletAccelerometerAccelerationCallbackThreshold
    option::Char
    min_x::Integer
    max_x::Integer
    min_y::Integer
    max_y::Integer
    min_z::Integer
    max_z::Integer
end

export BrickletAccelerometerConfiguration
struct BrickletAccelerometerConfiguration
    data_rate::Integer
    full_scale::Integer
    filter_bandwidth::Integer
end

export BrickletAccelerometerIdentity
struct BrickletAccelerometerIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAccelerometer
"""
Measures acceleration in three axis
"""
mutable struct BrickletAccelerometer <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAccelerometer(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_accelerometer")
        deviceInternal = package.BrickletAccelerometer(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the acceleration in x, y and z direction. The values
are given in gₙ/1000 (1gₙ = 9.80665m/s²). The range is
configured with :func:`Set Configuration`.

If you want to get the acceleration periodically, it is recommended
to use the :cb:`Acceleration` callback and set the period with
:func:`Set Acceleration Callback Period`.
"""
function get_acceleration(device::BrickletAccelerometer)
    return device.deviceInternal.get_acceleration()
end

export set_acceleration_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Acceleration` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Acceleration` callback is only triggered if the acceleration has
changed since the last triggering.
"""
function set_acceleration_callback_period(device::BrickletAccelerometer, period)
    device.deviceInternal.set_acceleration_callback_period(period)
end

export get_acceleration_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Acceleration Callback Period`.
"""
function get_acceleration_callback_period(device::BrickletAccelerometer)
    return device.deviceInternal.get_acceleration_callback_period()
end

export set_acceleration_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Acceleration Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the acceleration is *outside* the min and max values"
 "'i'",    "Callback is triggered when the acceleration is *inside* the min and max values"
 "'<'",    "Callback is triggered when the acceleration is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the acceleration is greater than the min value (max is ignored)"
"""
function set_acceleration_callback_threshold(device::BrickletAccelerometer, option, min_x, max_x, min_y, max_y, min_z, max_z)
    device.deviceInternal.set_acceleration_callback_threshold(option, min_x, max_x, min_y, max_y, min_z, max_z)
end

export get_acceleration_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Acceleration Callback Threshold`.
"""
function get_acceleration_callback_threshold(device::BrickletAccelerometer)
    return device.deviceInternal.get_acceleration_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Acceleration Reached`

is triggered, if the threshold

* :func:`Set Acceleration Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletAccelerometer, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletAccelerometer)
    return device.deviceInternal.get_debounce_period()
end

export get_temperature
"""
    $(SIGNATURES)

Returns the temperature of the accelerometer.
"""
function get_temperature(device::BrickletAccelerometer)
    return device.deviceInternal.get_temperature()
end

export set_configuration
"""
    $(SIGNATURES)

Configures the data rate, full scale range and filter bandwidth.
Possible values are:

* Data rate of 0Hz to 1600Hz.
* Full scale range of ±2gₙ up to ±16gₙ.
* Filter bandwidth between 50Hz and 800Hz.

Decreasing data rate or full scale range will also decrease the noise on
the data.
"""
function set_configuration(device::BrickletAccelerometer, data_rate, full_scale, filter_bandwidth)
    device.deviceInternal.set_configuration(data_rate, full_scale, filter_bandwidth)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletAccelerometer)
    return device.deviceInternal.get_configuration()
end

export led_on
"""
    $(SIGNATURES)

Enables the LED on the Bricklet.
"""
function led_on(device::BrickletAccelerometer)
    device.deviceInternal.led_on()
end

export led_off
"""
    $(SIGNATURES)

Disables the LED on the Bricklet.
"""
function led_off(device::BrickletAccelerometer)
    device.deviceInternal.led_off()
end

export is_led_on
"""
    $(SIGNATURES)

Returns *true* if the LED is enabled, *false* otherwise.
"""
function is_led_on(device::BrickletAccelerometer)
    return device.deviceInternal.is_led_on()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAccelerometer)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAccelerometer, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
module Groups

export
    # Group types
    Group,

    OrthogonalGroup, O3,
    SymmetricGroup, Sn,

    # General group methods
    eye,
    rand,
    nearest,
    *,
    inv,
    defrep, defrep!,
    irrep,
    real_irrep,

    # Specific group methods

        # Orthogonal group

        wigner_d̃,
        wigner_D̃,
        wigner_d,
        wigner_D,
        wigner_d̃_matrix,
        wigner_D̃_matrix,
        wigner_d_matrix,
        wigner_D_matrix,

        irrep_O3!,
        irrep_O3,

        # Symmetric group

        standard_tableaux,
        yor

include("types.jl")

include("orthogonal_group/O3.jl")
include("orthogonal_group/O3_irreps.jl")
include("orthogonal_group/O3_real_irreps.jl")

include("symmetric_group/Sn.jl")
include("symmetric_group/Sn_irreps.jl")

end


text:
# indexing

## utilities

using Base.Cartesian

## logical indexing

Base.getindex(xs::AnyCuArray, bools::AbstractArray{Bool}) = getindex(xs, CuArray(bools))

function Base.getindex(xs::AnyCuArray{T}, bools::AnyCuArray{Bool}) where {T}
  bools = reshape(bools, prod(size(bools)))
  indices = cumsum(bools)  # unique indices for elements that are true

  n = @allowscalar indices[end]  # number that are true
  ys = CuArray{T}(undef, n)

  if n > 0
    function kernel(ys::CuDeviceArray{T}, xs::CuDeviceArray{T}, bools, indices)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if i <= length(xs) && bools[i]
            b = indices[i]   # new position
            ys[b] = xs[i]
        end

        return
    end

    kernel = @cuda name="logical_getindex" launch=false kernel(ys, xs, bools, indices)
    config = launch_configuration(kernel.fun)
    threads = Base.min(length(indices), config.threads)
    blocks = cld(length(indices), threads)
    kernel(ys, xs, bools, indices; threads=threads, blocks=blocks)
  end

  unsafe_free!(indices)

  return ys
end

## find*

function Base.findall(bools::AnyCuArray{Bool})
    I = keytype(bools)
    indices = cumsum(reshape(bools, prod(size(bools))))

    n = @allowscalar indices[end]
    ys = CuArray{I}(undef, n)

    if n > 0
        function kernel(ys::CuDeviceArray, bools, indices)
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if i <= length(bools) && bools[i]
                i′ = CartesianIndices(bools)[i]
                b = indices[i]   # new position
                ys[b] = i′
            end

            return
        end

        kernel = @cuda name="findall" launch=false kernel(ys, bools, indices)
        config = launch_configuration(kernel.fun)
        threads = Base.min(length(indices), config.threads)
        blocks = cld(length(indices), threads)
        kernel(ys, bools, indices; threads=threads, blocks=blocks)
    end

    unsafe_free!(indices)

    return ys
end

function Base.findall(f::Function, A::AnyCuArray)
    bools = map(f, A)
    ys = findall(bools)
    unsafe_free!(bools)
    return ys
end

function Base.findfirst(testf::Function, xs::AnyCuArray)
    I = keytype(xs)

    y = CuArray([typemax(Int)])

    function kernel(y::CuDeviceArray, xs::CuDeviceArray)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if i <= length(xs) && testf(xs[i])
            @atomic y[1] = Base.min(y[1], i)
        end

        return
    end

    kernel = @cuda name="findfirst" launch=false kernel(y, xs)
    config = launch_configuration(kernel.fun)
    threads = Base.min(length(xs), config.threads)
    blocks = cld(length(xs), threads)
    kernel(y, xs; threads=threads, blocks=blocks)

    first_i = @allowscalar y[1]
    return first_i == typemax(Int) ? nothing : keys(xs)[first_i]
end

Base.findfirst(xs::AnyCuArray{Bool}) = findfirst(identity, xs)

function Base.findmin(a::AnyCuArray; dims=:)
    if dims == Colon()
        m = minimum(a)
        i = findfirst(x->x==m, a)
        return m,i
    else
        minima = minimum(a; dims=dims)
        i = findfirstval(minima, a)
        return minima,i
    end
end

function Base.findmax(a::AnyCuArray; dims=:)
    if dims == Colon()
        m = maximum(a)
        i = findfirst(x->x==m, a)
        return m,i
    else
        maxima = maximum(a; dims=dims)
        i = findfirstval(maxima, a)
        return maxima,i
    end
end

function findfirstval(vals::AnyCuArray, xs::AnyCuArray)
    ## find the first matching element

    # NOTE: this kernel performs global atomic operations for the sake of simplicity.
    #       if this turns out to be a bottleneck, we will need to cache in local memory.
    #       that requires the dimension-under-reduction to be iterated in first order.
    #       this can be done by splitting the iteration domain eagerly; see the
    #       accumulate kernel for an example, or git history from before this comment.

    indices = fill(typemax(Int), size(vals))

    function kernel(xs, vals, indices)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        R = CartesianIndices(xs)

        if i <= length(R)
            I = R[i]
            Jmax = last(CartesianIndices(vals))
            J = Base.min(I, Jmax)

            @inbounds if xs[I] == vals[J]
                I′ = LinearIndices(xs)[I]      # atomic_min only works with integers
                J′ = LinearIndices(indices)[J] # FIXME: @atomic doesn't handle array ref with CartesianIndices
                @atomic indices[J′] = Base.min(indices[J′], I′)
            end
        end

        return
    end

    kernel = @cuda launch=false kernel(xs, vals, indices)
    config = launch_configuration(kernel.fun)
    threads = Base.min(length(xs), config.threads)
    blocks = cld(length(xs), threads)
    kernel(xs, vals, indices; threads=threads, blocks=blocks)

    ## convert the linear indices to an appropriate type

    kt = keytype(xs)

    if kt == Int
        return indices
    else
        indices′ = CuArray{kt}(undef, size(indices))
        broadcast!(indices′, indices, Ref(keys(xs))) do index, keys
            keys[index]
        end

        return indices′
    end
end


text:
using CSV, GLM, Plots; pyplot()

df = CSV.read("../data/IQalc.csv")
groupA = df[df.Group .== "A", :]
groupB = df[df.Group .== "B", :]
groupC = df[df.Group .== "C", :]

model  = fit(LinearModel, @formula(AlcConsumption ~ IQ), df)
modelA = fit(LinearModel, @formula(AlcConsumption ~ IQ), groupA)
modelB = fit(LinearModel, @formula(AlcConsumption ~ IQ), groupB)
modelC = fit(LinearModel, @formula(AlcConsumption ~ IQ), groupC)

pred(x)  = coef(model)'  * [1, x]
predA(x) = coef(modelA)' * [1, x]
predB(x) = coef(modelB)' * [1, x]
predC(x) = coef(modelC)' * [1, x]

xlims = collect(extrema(df.IQ))

p1 = scatter(df.IQ, df.AlcConsumption, c=:black, msw=0, ma=0.2, label="")
     plot!(xlims, pred.(xlims), c=:black, label="All data")

p2 = scatter(groupA.IQ, groupA.AlcConsumption, c=:blue, msw=0, ma=0.2, label="")
     scatter!(groupB.IQ, groupB.AlcConsumption, c=:red, msw=0, ma=0.2, label="")
     scatter!(groupC.IQ, groupC.AlcConsumption, c=:green,msw=0, ma=0.2, label="")
     plot!(xlims, predA.(xlims), c=:blue, label="Group A")
     plot!(xlims, predB.(xlims), c=:red, label="Group B")
     plot!(xlims, predC.(xlims), c=:green, label="Group C")

plot(p1, p2, xlims=(xlims), ylims=(0,1),
	xlabel="IQ", ylabel="Alcohol Metric", size=(800,400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Define components for simulation and logging machines
# -------------------------------------------------------
# I choose Sim.jl to be self sufficient and not to depend
# on StateMachines.jl.
#

"supertype for state machines in `Sim.jl`"
abstract type SEngine end

"supertype for states"
abstract type SState end

"supertype for events"
abstract type SEvent end

"a state machine is undefined (after creation)"
struct Undefined <: SState end

"a state machine is idle"
struct Idle <: SState end

"a state machine is empty"
struct Empty <: SState end

"a state machine is busy"
struct Busy <: SState end

"a state machine is halted"
struct Halted <: SState end

"`Init(info)`: Init event with some info."
struct Init <: SEvent
    info::Any
end

"`Setup(vars::Array{Symbol,1}, scope::Module)`: setup a logger with some info."
struct Setup <: SEvent
    vars::Array{Symbol,1}
    scope::Module
end

"`Switch(to)`: switch to some other mode"
struct Switch <: SEvent
    to
end

"`Log()`: record command for logging"
struct Log <: SEvent end

"`Step()`: command"
struct Step <: SEvent end

"`Run()`: command"
struct Run <: SEvent
    duration::Float64
end

"`Start()`: command"
struct Start <: SEvent end

"`Stop()`: command"
struct Stop <: SEvent end

"`Resume()`: command"
struct Resume <: SEvent end

"`Clear()`: command"
struct Clear <: SEvent end

"""
    step!(A::SEngine, q::SState, σ::SEvent)

Default transition for clock and logger.

This is called if no otherwise defined transition occurs.

# Arguments
- `A::SEngine`: state machine for which a transition is called
- `q::SState`:  state of the state machine
- `σ::SEvent`:  event, triggering the transition
"""
function step!(A::SEngine, q::SState, σ::SEvent)
    println(stderr, "Warning: undefined transition ",
            "$(typeof(A)), ::$(typeof(q)), ::$(typeof(σ)))")
end


text:
#! /usr/bin/env julia

using Documenter

using GeometryBasics

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

# Copy the README to serve as the homepage
cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"))

makedocs(format=Documenter.HTML(), sitename="GeometryBasics.jl",
         pages=["index.md", "primitives.md", "rectangles.md", "polygons.md", "meshes.md",
                "decomposition.md", "distancefields.md", "metadata.md", "api.md"],
         modules=[GeometryBasics])

deploydocs(repo="github.com/JuliaGeometry/GeometryBasics.jl.git", push_preview=true)


text:
#=
   From BLOG example/mixture-of-gaussian.blog
   """
   A model for mixture of Gaussian distribution.
   The latent variable \verb|z| determines the mixture component of \verb|x|.
   The program asks the possible value of \verb|z| given the evidence of \verb|x| being 0.2.

   @author leili
   @since 2014-06-16
   """

  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           z    0.5481    0.4977     0.0050    0.0047   8954.7436    1.0004     1230.7234
       xpost    0.0560    1.1241     0.0112    0.0108   8908.5957    0.9999     1224.3809

   Cf ~/webppl/mixture_of_gaussian.wppl

=#

using Turing, StatsPlots

include("jl_utils.jl")

@model function mixture_of_gaussian(x=0.2)
    z ~ flip(0.5)
    x ~ z == true ? Normal(0.5, 1.0) : Normal(-0.5, 1.0)

    # Posterior prediction of x
    xpost ~ z == true ? Normal(0.5, 1.0) : Normal(-0.5, 1.0)
end

x=0.2
model = mixture_of_gaussian(x)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export Psrc,
       ModSpt,
       InitModSpt,
       WriteModSpt,
       InfoModSpts,
       MultiStepForward_MOR,
       ModSpts2Spts,
       SrcProj,
       AddPsrc!,
       QFQ,
       lowRankFd,
       ReadRDFD,
       QRbase,
       QR,
       randMix,
       formBaseMatrix,
       formStateMatrix,
       WriteBase

include("randQR.jl")

max_stars_repo_name: WneleiGao/ElasticWave


text:
using DustExtinction
using Test

# =============================================================================
# ccm89

# NOTE: Test is only to precision of 0.015 because there is a discrepancy
# of 0.014 for the B band wavelength of unknown origin (and up to 0.002 in
# other bands).
#
# Note that a and b can be obtained with:
# b = ccm89(wave, 0.)
# a = ccm89(wave, 1.) - b
#
# These differ from the values tablulated in the original paper.
# Could be due to floating point errors in the original paper?

# U, B, V, R, I, J, H, K band effective wavelengths from CCM '89 table 3
x_inv_microns = [2.78; 2.27; 1.82; 1.43; 1.11; 0.80; 0.63; 0.46]
wave = 1.e4 ./ x_inv_microns

# A(lambda)/A(V) for R_V = 3.1 from Table 3 of CCM '89
ref_values = [1.569, 1.337, 1.000, 0.751, 0.479, 0.282, 0.190, 0.114]

for i=1:length(wave)
    @test ccm89(wave[i], 3.1) ≈ ref_values[i] rtol=0.016
end

# =============================================================================
# od94

# NOTE: The tabulated values go to 0.001, but the test is only for matching
# at the 0.005 level, because there is currently a discrepancy up to 0.0047
# of unknown origin.

# Tests od94() at Rv = 3.1 against the widely used values tabulated in
# Schlegel, Finkbeiner and Davis (1998)
# http://adsabs.harvard.edu/abs/1998ApJ...500..525S

# This is tested by evaluating the extinction curve at a (given)
# effective wavelength, since these effective wavelengths:
# "... represent(s) that wavelength on the extinction curve
# with the same extinction as the full passband."

# The test does not include UKIRT L' (which, at 3.8 microns) is
# beyond the range of wavelengths allowed by the function
# or the APM b_J filter which is defined in a non-standard way.

# The SFD98 tabulated values go to 1e-3, so we should be able to match at
# that level.

wave = [3372., 4404., 5428., 6509., 8090.,
        3683., 4393., 5519., 6602., 8046.,
        12660., 16732., 22152.,
        5244., 6707., 7985., 9055.,
        6993.,
        3502., 4676., 4127.,
        4861., 5479.,
        3546., 4925., 6335., 7799., 9294.,
        3047., 4711., 5498.,
        6042., 7068., 8066.,
        4814., 6571., 8183.]

ref_values = [1.664, 1.321, 1.015, 0.819, 0.594,
              1.521, 1.324, 0.992, 0.807, 0.601,
              0.276, 0.176, 0.112,
              1.065, 0.793, 0.610, 0.472,
              0.755,
              1.602, 1.240, 1.394,
              1.182, 1.004,
              1.579, 1.161, 0.843, 0.639, 0.453,
              1.791, 1.229, 0.996,
              0.885, 0.746, 0.597,
              1.197, 0.811, 0.580]

for i=1:length(wave)
    @test od94(wave[i], 3.1) ≈ ref_values[i] rtol=0.0051
end

# =============================================================================
# SFD98 dust maps

# refebv obtained using http://irsa.ipac.caltech.edu/applications/DUST/
# and manually inserting the following lines and reading off the values for
# SFD (1998) reference pixel:
#
# 0. 0. gal
# 90. 0. gal
# 180. 0. gal
# 270. 0. gal
# 5.729577951308233 5.729577951308233 gal
# 5.729577951308233 -5.729577951308233 gal
#
# Note that values are not expected to agree exactly because (1) IRSA
# seems to not do linear interpolation and (2) IRSA seems to use a
# resampled map (pixel values reported by IRSA don't seem to match any
# pixel values in original maps).

refcoords = [(0., 0.),
             (pi/2., 0.),
             (pi, 0.),
             (3pi/2., 0.),
             (0.1, 0.1),
             (0.1, -0.1)]

refebv = [100.0270,
          2.6185,
          1.4182,
          3.4194,
          0.7949,
          0.5680]

if haskey(ENV, "SFD98_DIR")
    dustmap = SFD98Map()
    for i=1:length(refcoords)
        l, b = refcoords[i]
        @test ebv_galactic(dustmap, l, b) ≈ refebv[i] rtol=0.02
    end
else
    println("Skipping SFD98Map test because \$SFD98_DIR not defined.")
end

println("Tests passed.")