text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

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

function drawcurve!(canvassize::Int, curve::Spline{P,S,N,M}, numpoints::Int; linewidth = 0.5, curvecolor = RGB(1, 0, 1), controlpointcolor = RGB(0, 0, 0), controlpointsize = 5, controlpolygoncolor = RGB(0, 0, 0)) where {P,S,N,M}
    step = 1.0 / numpoints

    point1 = point(curve, 0.0) .* canvassize
    # println("point1 $point1")
    point1 = [point1[1], canvassize - point1[2]] # flip y because of the coordinate system SVG uses

    Luxor.setline(linewidth)
    Luxor.sethue(curvecolor)

    for u in step:step:1.0
        point2 = point(curve, u) .* canvassize
        point2 = [point2[1], canvassize - point2[2]]

        Luxor.line(Luxor.Point(point1...), Luxor.Point(point2...), :fillstroke)
        point1 = copy(point2)
    end

    Luxor.sethue(controlpointcolor)
    # draw control points
    controlpoints = euclideancontrolpoints(curve)
    for point in controlpoints
        pt = point .* canvassize
        pt = [pt[1], canvassize - pt[2]]
        Luxor.circle(Luxor.Point(pt...), controlpointsize, :fill)
    end

    Luxor.sethue(controlpolygoncolor)
    for i in 1:(size(controlpoints)[1] - 1)
        pt1 = controlpoints[i] .* canvassize
        pt1 = [pt1[1], canvassize - pt1[2]]
        pt2 = controlpoints[i + 1] .* canvassize
        pt2 = [pt2[1], canvassize - pt2[2]]
        Luxor.line(Luxor.Point(pt1...), Luxor.Point(pt2...), :fillstroke)
    end
end

function drawcurves(curves::Vararg{Spline{P,S,N,M}}; numpoints::Int = 200, canvassize::Int = 2000) where {P,S,N,M}
    canvas = Luxor.Drawing(canvassize, canvassize)
    Luxor.background("white")
    drawcurve!(canvassize, curves[1], numpoints, linewidth = 5, curvecolor = RGB(0, 1, 1))
    for curve in curves[2:end]
        drawcurve!(canvassize, curve, numpoints, linewidth = 1, controlpointcolor = RGB(1, 0, 0), controlpointsize = 3, controlpolygoncolor = RGB(0, 1, 0))
    end
    Luxor.finish()
    Luxor.preview()
end

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

# all functions follow the pattern draw(obj) and draw!(scene, obj) where the first case draws the object in a blank scene and displays it, and
# the second case draws the object in an existing scene, draw!(obj) can also be used to draw the object in the current scene

global current_main_scene = nothing
global current_layout_scene = nothing
global current_3d_scene = nothing
global current_mode = nothing           # modes:    nothing, :default  -> Original Vis beheviour
                                        #           :pluto             -> support pluto notebooks
                                        #           :docs              -> support documenter figures

# added the following 2 functions to allow us to hack the drawing mechanisim while in a pluto notebook
set_current_main_scene(scene) = (global current_main_scene = scene)
set_current_3d_scene(lscene) = (global current_3d_scene = lscene)

get_current_mode() = begin global current_mode; return current_mode end
set_current_mode(mode) = (global current_mode = mode)

show(image) = imshow(image)
display(scene = current_main_scene) = begin
    global current_mode;
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
    Makie.display(scene)
end

"""
    scene(resolution = (1000, 1000))

Create a new Makie scene with the given resolution including control buttons.
"""
function scene(resolution = (1000, 1000))
    @assert resolution[1] > 0 && resolution[2] > 0

    scene, layout = Makie.layoutscene(resolution = resolution)
    global current_main_scene = scene
    global current_layout_scene = layout
    lscene = layout[1, 1] = Makie.LScene(scene, scenekw = (camera = Makie.cam3d_cad!, axis_type = Makie.axis3d!, raw = false))
    global current_3d_scene = lscene

    # in these modes we want to skip the creation of the utility buttons as these modes are not interactive
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene, lscene
    end

    threedbutton = Makie.Button(scene, label = "3D", buttoncolor = RGB(0.8, 0.8, 0.8), height = 40, width = 80)
    twodxbutton = Makie.Button(scene, label = "2D-x", buttoncolor = RGB(0.8, 0.8, 0.8), height = 40, width = 80)
    twodybutton = Makie.Button(scene, label = "2D-y", buttoncolor = RGB(0.8, 0.8, 0.8), height = 40, width = 80)
    savebutton = Makie.Button(scene, label = "Screenshot", buttoncolor = RGB(0.8, 0.8, 0.8), height = 40, width = 160)

    Makie.on(threedbutton.clicks) do nclicks
        make3d(lscene)
        yield()
    end

    Makie.on(twodybutton.clicks) do nclicks
        make2dy(lscene)
        yield()
    end

    Makie.on(twodxbutton.clicks) do nclicks
        make2dx(lscene)
        yield()
    end

    Makie.on(savebutton.clicks) do nclicks
        Vis.save("screenshot.png")
        yield()
    end

    layout[2, 1] = Makie.grid!(hcat(threedbutton, twodxbutton, twodybutton, savebutton), tellwidth = false, tellheight = true)
    return scene, lscene
end

function make3d(scene::Makie.LScene = current_3d_scene)
    s = scene.scene
    # use 3d camera
    Makie.cam3d_cad!(s)
    # reset scene rotation
    s.transformation.rotation[] = Makie.Quaternion(0.0, 0.0, 0.0, 1.0)
    # show all the axis ticks
    s[Makie.OldAxis].attributes.showticks[] = (true, true, true)
    # reset tick and axis label rotation and position
    rot = (Makie.Quaternion(0.0, 0.0, -0.7071067811865476, -0.7071067811865475), Makie.Quaternion(0.0, 0.0, 1.0, 0.0), Makie.Quaternion(0.0, 0.7071067811865475, 0.7071067811865476, 0.0))
    s[Makie.OldAxis].attributes.ticks.rotation[] = rot
    s[Makie.OldAxis].attributes.names.rotation[] = rot
    s[Makie.OldAxis].attributes.ticks.align = ((:left, :center), (:right, :center), (:right, :center))
    s[Makie.OldAxis].attributes.names.align = ((:left, :center), (:right, :center), (:right, :center))
    # reset scene limits to automatic
    s.limits[] = Makie.Automatic()
    Makie.update!(s)
end

function make2dy(scene::Makie.LScene = current_3d_scene)
    s = scene.scene
    # use 2d camera
    Makie.cam2d!(s)

    scene_transform = Makie.qrotation(Makie.Vec3f0(0, 1, 0), 0.5pi)
    scene_transform_inv = Makie.qrotation(Makie.Vec3f0(0, 1, 0), -0.5pi)    # to use with the ticks and names

    # set rotation to look onto yz plane
    s.transformation.rotation[] = scene_transform
    # hide x ticks

    # there is a bug in Makie 0.14.2 which cause an exception setting the X showticks to false.
    # we work wround it by making sure the labels we want to turn off are ortogonal to the view direction
    # s[Makie.OldAxis].attributes.showticks[] = (false, true, true)
    s[Makie.OldAxis].attributes.showticks[] = (true, true, true)

    # set tick and axis label rotation and position
    s[Makie.OldAxis].attributes.ticks.rotation[] = (0.0, scene_transform_inv, scene_transform_inv)
    s[Makie.OldAxis].attributes.names.rotation[] = s[Makie.OldAxis].attributes.ticks.rotation[]
    s[Makie.OldAxis].attributes.ticks.align = ((:right, :center), (:right, :center), (:center, :right))
    s[Makie.OldAxis].attributes.names.align = ((:left, :center), (:left, :center), (:center, :left))
    # update the scene limits automatically to get true reference values
    s.limits[] = Makie.Automatic()
    Makie.update_limits!(s)
    # manually set the scene limits to draw the axes correctly
    o, w = Makie.origin(s.data_limits[]), Makie.widths(s.data_limits[])
    s.limits[] = Makie.FRect3D((1000.0f0, o[2], o[3]), (w[2], w[2], w[3]))
    # set the eye (i.e. light) position to behind the camera
    s.camera.eyeposition[] = (0, 0, -100)
    Makie.update!(s)
end

function make2dx(scene::Makie.LScene = current_3d_scene)
    s = scene.scene
    # use 2d camera
    Makie.cam2d!(s)

    scene_transform= Makie.qrotation(Makie.Vec3f0(0, 0, 1), 0.5pi) * Makie.qrotation(Makie.Vec3f0(1, 0, 0), 0.5pi)
    scene_transform_inv=Makie.qrotation(Makie.Vec3f0(1, 0, 0), -0.5pi) * Makie.qrotation(Makie.Vec3f0(0, 0, 1), -0.5pi)

    # set rotation to look onto yz plane
    s.transformation.rotation[] = scene_transform
    # hide y ticks

    # there is a bug in Makie 0.14.2 which cause an exception setting the X showticks to false.
    # we work wround it by making sure the labels we want to turn off are ortogonal to the view direction
    s[Makie.OldAxis].attributes.showticks[] = (true, false, true)

    # set tick and axis label rotation and position
    s[Makie.OldAxis].attributes.ticks.rotation[] = (scene_transform_inv, 0.0, scene_transform_inv)
    s[Makie.OldAxis].attributes.names.rotation[] = s[Makie.OldAxis].attributes.ticks.rotation[]
    s[Makie.OldAxis].attributes.ticks.align = ((:right, :center), (:right, :center), (:center, :center))
    s[Makie.OldAxis].attributes.names.align = ((:left, :center), (:right, :center), (:center, :center))
    # update the scene limits automatically to get true reference values
    s.limits[] = Makie.Automatic()
    Makie.update_limits!(s)
    # manually set the scene limits to draw the axes correctly
    o, w = Makie.origin(s.data_limits[]), Makie.widths(s.data_limits[])
    s.limits[] = Makie.FRect3D((o[1], -1000.0f0, o[3]), (w[1], w[1], w[3]))
    # set the eye (i.e. light) position to behind the camera
    s.camera.eyeposition[] = (0, 0, -100)
    Makie.update!(s)
end

"""
    draw(ob; resolution = (1000, 1000), kwargs...)

Draw an object in a new scene.
`kwargs` depends on the object type.
"""
function draw(ob; resolution = (1000, 1000), kwargs...)
    scene, lscene = Vis.scene(resolution)
    draw!(lscene, ob; kwargs...)
    display(scene)

    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

"""
    draw!([scene = currentscene], ob; kwargs...)

Draw an object in an existing scene.
`kwargs` depends on the object type.
"""
function draw!(ob; kwargs...)
    if current_3d_scene === nothing
        scene, lscene = Vis.scene()
    else
        scene = current_main_scene
        lscene = current_3d_scene
    end
    draw!(lscene, ob; kwargs...)
    display(scene)

    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

"""
    save(path::String)

Save the current Makie scene to an image file.
"""
function save(path::String)
    # save closes the window so just display it again as a work-around
    # for some reason the size isn't maintained automatically so we just reset it manually
    size = Makie.size(current_main_scene)
    Makie.save(path, current_3d_scene.scene)
    Makie.resize!(current_main_scene, size)
    display(current_main_scene)
end
function save(::Nothing) end

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

function λtoRGB(λ::T, gamma::T = 0.8) where {T<:Real}
    wavelength = λ * 1000 # λ is in um, need in nm
    if (wavelength >= 380 && wavelength <= 440)
        attenuation = 0.3 + 0.7 * (wavelength - 380) / (440 - 380)
        R = ((-(wavelength - 440) / (440 - 380)) * attenuation)^gamma
        G = 0.0
        B = (1.0 * attenuation)^gamma
    elseif (wavelength >= 440 && wavelength <= 490)
        R = 0.0
        G = ((wavelength - 440) / (490 - 440))^gamma
        B = 1.0
    elseif (wavelength >= 490 && wavelength <= 510)
        R = 0.0
        G = 1.0
        B = (-(wavelength - 510) / (510 - 490))^gamma
    elseif (wavelength >= 510 && wavelength <= 580)
        R = ((wavelength - 510) / (580 - 510))^gamma
        G = 1.0
        B = 0.0
    elseif (wavelength >= 580 && wavelength <= 645)
        R = 1.0
        G = (-(wavelength - 645) / (645 - 580))^gamma
        B = 0.0
    elseif (wavelength >= 645 && wavelength <= 750)
        attenuation = 0.3 + 0.7 * (750 - wavelength) / (750 - 645)
        R = (1.0 * attenuation)^gamma
        G = 0.0
        B = 0.0
    else
        R = 0.0
        G = 0.0
        B = 0.0
    end
    return RGB(R, G, B)
end

indexedcolor(i::Int) = ColorSchemes.hsv[rem(i / (2.1 * π), 1.0)]
indexedcolor2(i::Int) = ColorSchemes.hsv[1.0 - rem(i / (2.1 * π), 1.0)] .* 0.5

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

## displaying imported meshes

Base.:*(a::Transform, p::GeometryBasics.PointMeta) = a * p.main
Base.:*(a::Real, p::GeometryBasics.PointMeta{N,S}) where {S<:Real,N} = GeometryBasics.Point{N,S}((a * SVector{N,S}(p))...)
Base.:*(a::Transform, p::GeometryBasics.Point{N,S}) where {S<:Real,N} = GeometryBasics.Point{N,S}((a.rotation * SVector{N,S}(p) + a.translation)...)

function draw!(scene::Makie.LScene, ob::AbstractString; color = :gray, linewidth = 3, shaded::Bool = true, wireframe::Bool = false, transform::Transform{Float64} = identitytransform(Float64), scale::Float64 = 1.0, kwargs...)
    if any(endswith(lowercase(ob), x) for x in [".obj", "ply", ".2dm", ".off", ".stl"])
        meshdata = FileIO.load(ob)
        if transform != identitytransform(Float64) || scale != 1.0
            coords = [transform * (scale * p) for p in GeometryBasics.coordinates(meshdata)]
            meshdata = GeometryBasics.Mesh(coords, GeometryBasics.faces(meshdata))
        end
        Makie.mesh!(GeometryBasics.normal_mesh(meshdata); kwargs..., color = color, shading = shaded, visible = shaded)
        if wireframe
            if shaded
                Makie.wireframe!(scene[end][1], color = (:black, 0.1), linewidth = linewidth)
            else
                Makie.wireframe!(scene[end][1], color = color, linewidth = linewidth)
            end
        end
    else
        @error "Unsupported file type"
    end
end

## GEOMETRY

"""
    draw!(scene::Makie.LScene, surf::Surface{T}; numdivisions = 20, normals = false, normalcolor = :blue, kwargs...)

Transforms `surf` into a mesh using [`makemesh`](@ref) and draws the result.
`normals` of the surface can be drawn at evenly sampled points with provided `normalcolor`.
`numdivisions` determines the resolution with which the mesh is triangulated.
`kwargs` is passed on to the [`TriangleMesh`](@ref) drawing function.
"""
function draw!(scene::Makie.LScene, surf::Surface{T}; numdivisions::Int = 30, normals::Bool = false, normalcolor = :blue, kwargs...) where {T<:Real}
    mesh = makemesh(surf, numdivisions)
    if nothing === mesh
        return
    end
    draw!(scene, mesh; kwargs..., normals = false)
    if normals
        ndirs = Makie.Point3f0.(samplesurface(surf, normal, numdivisions ÷ 10))
        norigins = Makie.Point3f0.(samplesurface(surf, point, numdivisions ÷ 10))
        Makie.arrows!(scene, norigins, ndirs, arrowsize = 0.2, arrowcolor = normalcolor, linecolor = normalcolor, linewidth = 2)
    end
end

"""
    draw!(scene::Makie.LScene, tmesh::TriangleMesh{T}; linewidth = 3, shaded = true, wireframe = false, color = :orange, normals = false, normalcolor = :blue, transparency = false, kwargs...)

Draw a [`TriangleMesh`](@ref), optionially with a visible `wireframe`. `kwargs` are passed on to [`Makie.mesh`](http://makie.juliaplots.org/stable/plotting_functions.html#mesh).
"""
function draw!(scene::Makie.LScene, tmesh::TriangleMesh{T}; linewidth = 3, shaded::Bool = true, wireframe::Bool = false, color = :orange, normals::Bool = false, normalcolor = :blue, transparency::Bool = false, kwargs...) where {T<:Real}
    points, indices = makiemesh(tmesh)
    if length(points) > 0 && length(indices) > 0
        Makie.mesh!(scene, points, indices; kwargs..., color = color, shading = shaded, transparency = transparency, visible = shaded)
        if wireframe
            mesh = scene.scene[end][1]
            if shaded
                Makie.wireframe!(scene, mesh, color = (:black, 0.1), linewidth = linewidth)
            else
                Makie.wireframe!(scene, mesh, color = color, linewidth = linewidth)
            end
        end
    end
    if normals
        @warn "Normals being drawn from triangulated mesh, precision may be low"
        norigins = [Makie.Point3f0(centroid(t)) for t in tmesh.triangles[1:10:end]]
        ndirs = [Makie.Point3f0(normal(t)) for t in tmesh.triangles[1:10:end]]
        if length(norigins) > 0
            Makie.arrows!(scene, norigins, ndirs, arrowsize = 0.2, arrowcolor = normalcolor, linecolor = normalcolor, linewidth = 2)
        end
    end
end

"""
    draw!(scene::Makie.LScene, meshes::Vararg{S}; colors::Bool = false, kwargs...) where {T<:Real,S<:Union{TriangleMesh{T},Surface{T}}}

Draw a series of [`TriangleMesh`](@ref) or [`Surface`](@ref) objects, if `colors` is true then each mesh will be colored automatically with a diverse series of colors.
`kwargs` are is passed on to the drawing function for each element.
"""
function draw!(scene::Makie.LScene, meshes::Vararg{S}; colors::Bool = false, kwargs...) where {T<:Real,S<:Union{TriangleMesh{T},Surface{T}}}
    for i in 1:length(meshes)
        if colors
            col = indexedcolor2(i)
        else
            col = :orange
        end
        draw!(scene, meshes[i]; kwargs..., color = col)
    end
end

function draw(meshes::Vararg{S}; kwargs...) where {T<:Real,S<:Union{TriangleMesh{T},Surface{T}}}
    scene, lscene = Vis.scene()
    draw!(lscene, meshes...; kwargs...)
    Makie.display(scene)
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

"""
    draw!(scene::Makie.LScene, csg::Union{CSGTree,CSGGenerator}; numdivisions::Int = 20, kwargs...)

Convert a CSG object ([`CSGTree`](@ref) or [`CSGGenerator`](@ref)) to a mesh using [`makemesh`](@ref) with resolution set by `numdivisions` and draw the resulting [`TriangleMesh`](@ref).
"""
draw!(scene::Makie.LScene, csg::CSGTree{T}; numdivisions::Int = 30, kwargs...) where {T<:Real} = draw!(scene, makemesh(csg, numdivisions); kwargs...)
draw!(scene::Makie.LScene, csg::CSGGenerator{T}; kwargs...) where {T<:Real} = draw!(scene, csg(); kwargs...)

"""
    draw!(scene::Makie.LScene, bbox::BoundingBox{T}; kwargs...)

Draw a [`BoundingBox`](@ref) as a wireframe, ie series of lines.
"""
function draw!(scene::Makie.LScene, bbox::BoundingBox{T}; kwargs...) where {T<:Real}
    p1 = SVector{3,T}(bbox.xmin, bbox.ymin, bbox.zmin)
    p2 = SVector{3,T}(bbox.xmin, bbox.ymax, bbox.zmin)
    p3 = SVector{3,T}(bbox.xmin, bbox.ymax, bbox.zmax)
    p4 = SVector{3,T}(bbox.xmin, bbox.ymin, bbox.zmax)
    p5 = SVector{3,T}(bbox.xmax, bbox.ymin, bbox.zmin)
    p6 = SVector{3,T}(bbox.xmax, bbox.ymax, bbox.zmin)
    p7 = SVector{3,T}(bbox.xmax, bbox.ymax, bbox.zmax)
    p8 = SVector{3,T}(bbox.xmax, bbox.ymin, bbox.zmax)
    Makie.linesegments!(scene, [p1, p2, p2, p3, p3, p4, p4, p1, p1, p5, p2, p6, p3, p7, p4, p8, p5, p6, p6, p7, p7, p8, p8, p5]; kwargs...)
end

## OPTICS

"""
    draw!(scene::Makie.LScene, ass::LensAssembly; kwargs...)

Draw each element in a [`LensAssembly`](@ref), with each element automatically colored differently.
"""
function draw!(scene::Makie.LScene, ass::LensAssembly{T}; kwargs...) where {T<:Real}
    for (i, e) in enumerate(elements(ass))
        draw!(scene, e; kwargs..., color = indexedcolor2(i))
    end
end

"""
    draw!(scene::Makie.LScene, sys::AbstractOpticalSystem; kwargs...)

Draw each element in the lens assembly of an [`AbstractOpticalSystem`](@ref), with each element automatically colored differently, as well as the detector of the system.
"""
function draw!(scene::Makie.LScene, sys::CSGOpticalSystem{T}; kwargs...) where {T<:Real}
    draw!(scene, sys.assembly; kwargs...)
    draw!(scene, sys.detector; kwargs...)
end

draw!(scene::Makie.LScene, sys::AxisymmetricOpticalSystem{T}; kwargs...) where {T<:Real} = draw!(scene, sys.system; kwargs...)

onlydetectorrays(system::Q, tracevalue::LensTrace{T,3}) where {T<:Real,Q<:AbstractOpticalSystem{T}} = onsurface(detector(system), point(tracevalue))

## RAY GEN
"""
    drawtracerays(system::Q; raygenerator::S = Source(transform = Transform.translation(0.0,0.0,10.0), origins = Origins.RectGrid(10.0,10.0,25,25),directions = Constant(0.0,0.0,-1.0)), test::Bool = false, trackallrays::Bool = false, colorbysourcenum::Bool = false, colorbynhits::Bool = false, rayfilter::Union{Nothing,Function} = onlydetectorrays, kwargs...)

Displays a model of the optical system.
`raygenerator` is an iterator that generates rays.
If `trackallrays` is true then ray paths from the emitter will be displayed otherwise just the final rays that intersect the image detector will be shown, not the entire ray path.
`colorbysourcenum` and `colorbynhits` will color rays accordingly, otherwise rays will be colored according to their wavelength.

By default only ray paths that eventually intersect the detector surface are displayed. If you want to display all ray paths set `rayfilter = nothing`.

Also `drawtracerays!` to add to an existing scene, with `drawsys` and `drawgen` to specify whether `system` and `raygenerator` should be drawn respectively.
"""
function drawtracerays(system::Q; raygenerator::S = Source(transform = translation(0.0,0.0,10.0), origins = Origins.RectGrid(10.0,10.0,25,25),directions = Constant(0.0,0.0,-1.0)), test::Bool = false, trackallrays::Bool = false, colorbysourcenum::Bool = false, colorbynhits::Bool = false, rayfilter::Union{Nothing,Function} = onlydetectorrays, verbose::Bool = false, resolution::Tuple{Int,Int} = (1000, 1000), kwargs...) where {T<:Real,Q<:AbstractOpticalSystem{T},S<:AbstractRayGenerator{T}}
    verbose && println("Drawing System...")
    s, ls = Vis.scene(resolution)

    drawtracerays!(ls, system, raygenerator = raygenerator, test = test, colorbysourcenum = colorbysourcenum, colorbynhits = colorbynhits, rayfilter = rayfilter, trackallrays = trackallrays, verbose = verbose, drawsys = true, drawgen = true; kwargs...)

    display(s)
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return s
    end
end

drawtracerays!(system::Q; kwargs...) where {T<:Real,Q<:AbstractOpticalSystem{T}} = drawtracerays!(current_3d_scene, system; kwargs...)

function drawtracerays!(scene::Makie.LScene, system::Q; raygenerator::S = Source(transform = translation(0.0,0.0,10.0), origins = Origins.RectGrid(10.0,10.0,25,25),directions = Constant(0.0,0.0,-1.0)), test::Bool = false, trackallrays::Bool = false, colorbysourcenum::Bool = false, colorbynhits::Bool = false, rayfilter::Union{Nothing,Function} = onlydetectorrays, verbose::Bool = false, drawsys::Bool = false, drawgen::Bool = false, kwargs...) where {T<:Real,Q<:AbstractOpticalSystem{T},S<:AbstractRayGenerator{T}}
    raylines = Vector{LensTrace{T,3}}(undef, 0)

    drawgen && draw!(scene, raygenerator, norays = true; kwargs...)
    drawsys && draw!(scene, system; kwargs...)

    verbose && println("Tracing...")
    for (i, r) in enumerate(raygenerator)
        if i % 1000 == 0 && verbose
            print("\r $i / $(length(raygenerator))")
        end
        allrays = Vector{LensTrace{T,3}}(undef, 0)
        if trackallrays
            res = trace(system, r, trackrays = allrays, test = test)
        else
            res = trace(system, r, test = test)
        end

        if trackallrays && !isempty(allrays)
            if rayfilter === nothing || rayfilter(system, allrays[end])
                # filter on the trace that is hitting the detector
                for r in allrays
                    push!(raylines, r)
                end
            end
        elseif res !== nothing
            if rayfilter === nothing || rayfilter(system, res)
                push!(raylines, res)
            end
        end
    end
    verbose && print("\r")

    verbose && println("Drawing Rays...")
    draw!(scene, raylines, colorbysourcenum = colorbysourcenum, colorbynhits = colorbynhits; kwargs...)
end

"""
    drawtraceimage(system::Q; raygenerator::S = Source(transform = translation(0.0,0.0,10.0), origins = Origins.RectGrid(10.0,10.0,25,25),directions = Constant(0.0,0.0,-1.0)), test::Bool = false)

Traces rays from `raygenerator` through `system` and shows and returns the detector image. `verbose` will print progress updates.
"""
function drawtraceimage(system::Q; raygenerator::S = Source(transform = translation(0.0,0.0,10.0), origins = Origins.RectGrid(10.0,10.0,25,25),directions = Constant(0.0,0.0,-1.0)), test::Bool = false, verbose::Bool = false) where {T<:Real,Q<:AbstractOpticalSystem{T},S<:AbstractRayGenerator{T}}
    resetdetector!(system)
    start_time = time()
    for (i, r) in enumerate(raygenerator)
        if i % 1000 == 0 && verbose
            dif = round(time() - start_time, digits = 1)
            left = round((time() - start_time) / i * (length(raygenerator) - i), digits = 1)
            print("\rTraced: $i / $(length(raygenerator))     Elapsed: $(dif)s     Left: $(left)s")
        end
        res = trace(system, r, test = test)
    end
    verbose && print("\r")
    tot = round(time() - start_time, digits = 1)
    verbose && println("Completed in $(tot)s")
    show(detectorimage(system))
    return detectorimage(system)
end

## RAYS, LINES AND POINTS

"""
    draw!(scene::Makie.LScene, rays::AbstractVector{<:AbstractRay{T,N}}; kwargs...)

Draw a vector of [`Ray`](@ref) or [`OpticalRay`](@ref) objects.
"""
function draw!(scene::Makie.LScene, rays::AbstractVector{<:AbstractRay{T,N}}; kwargs...) where {T<:Real,N}
    for r in rays
        draw!(scene, r; kwargs...)
    end
end

"""
    draw!(scene::Makie.LScene, traces::AbstractVector{LensTrace{T,N}}; kwargs...)

Draw a vector of [`LensTrace`](@ref) objects.
"""
function draw!(scene::Makie.LScene, traces::AbstractVector{LensTrace{T,N}}; kwargs...) where {T<:Real,N}
    for t in traces
        draw!(scene, t; kwargs...)
    end
end

"""
    draw!(scene::Makie.LScene, trace::LensTrace{T,N}; colorbysourcenum::Bool = false, colorbynhits::Bool = false, kwargs...)

Draw a [`LensTrace`](@ref) as a line which can be colored automatically by its `sourcenum` or `nhits` attributes.
The alpha is determined by the `power` attribute of `trace`.
"""
function draw!(scene::Makie.LScene, trace::LensTrace{T,N}; colorbysourcenum::Bool = false, colorbynhits::Bool = false, kwargs...) where {T<:Real,N}
    if colorbysourcenum
        color = indexedcolor(sourcenum(trace))
    elseif colorbynhits
        color = indexedcolor(nhits(trace))
    else
        color = λtoRGB(wavelength(trace))
    end
    draw!(scene, (origin(ray(trace)), point(intersection(trace))); kwargs..., color = RGBA(color.r, color.g, color.b, sqrt(power(trace))), transparency = true)
end

"""
    draw!(scene::Makie.LScene, ray::OpticalRay{T,N}; colorbysourcenum::Bool = false, colorbynhits::Bool = false, kwargs...)

Draw an [`OpticalRay`](@ref) which can be colored automatically by its `sourcenum` or `nhits` attributes.
The alpha of the ray is determined by the `power` attribute of `ray`.
`kwargs` are passed to `draw!(scene, ray::Ray)`.
"""
function draw!(scene::Makie.LScene, r::OpticalRay{T,N}; colorbysourcenum::Bool = false, colorbynhits::Bool = false, kwargs...) where {T<:Real,N}
    if colorbysourcenum
        color = indexedcolor(sourcenum(r))
    elseif colorbynhits
        color = indexedcolor(nhits(r))
    else
        color = λtoRGB(wavelength(r))
    end
    draw!(scene, ray(r); kwargs..., color = RGBA(color.r, color.g, color.b, sqrt(power(r))), transparency = true, rayscale = power(r))
end

"""
    draw!(scene::Makie.LScene, ray::Ray{T,N}; color = :yellow, rayscale = 1.0, kwargs...)

Draw a [`Ray`](@ref) in a given `color` optionally scaling the size using `rayscale`.
`kwargs` are passed to [`Makie.arrows`](http://makie.juliaplots.org/stable/plotting_functions.html#arrows).
"""
function draw!(scene::Makie.LScene, ray::AbstractRay{T,N}; color = :yellow, rayscale = 1.0, kwargs...) where {T<:Real,N}
    arrow_size = min(0.05, rayscale * 0.05)
    Makie.arrows!(scene, [Makie.Point3f0(origin(ray))], [Makie.Point3f0(rayscale * direction(ray))]; kwargs..., arrowsize = arrow_size, arrowcolor = color, linecolor = color, linewidth=arrow_size * 0.5)
end

"""
    draw!(scene::Makie.LScene, du::DisjointUnion{T}; kwargs...)

Draw each [`Interval`](@ref) in a [`DisjointUnion`](@ref).
"""
function draw!(scene::Makie.LScene, du::DisjointUnion{T}; kwargs...) where {T<:Real}
    draw!(scene, intervals(du); kwargs...)
end

"""
    draw!(scene::Makie.LScene, intervals::AbstractVector{Interval{T}}; kwargs...)

Draw a vector of [`Interval`](@ref)s.
"""
function draw!(scene::Makie.LScene, intervals::AbstractVector{Interval{T}}; kwargs...) where {T<:Real}
    for i in intervals
        draw!(scene, i; kwargs...)
    end
end

"""
    draw!(scene::Makie.LScene, interval::Interval{T}; kwargs...)

Draw an [`Interval`](@ref) as a line with circles at each [`Intersection`](@ref) point.
"""
function draw!(scene::Makie.LScene, interval::Interval{T}; kwargs...) where {T<:Real}
    if !(interval isa EmptyInterval)
        l = lower(interval)
        u = upper(interval)
        if !(l isa RayOrigin)
            draw!(scene, l)
        else
            @warn "Negative half space, can't draw ray origin"
        end
        if !(u isa Infinity)
            draw!(scene, u)
        else
            @warn "Positive half space, can't draw end point"
        end
        if !(l isa RayOrigin) && !(u isa Infinity)
            draw!(scene, (point(l), point(u)); kwargs...)
        end
    end
end

"""
    draw!(scene::Makie.LScene, intersection::Intersection; normal::Bool = false, kwargs...)

Draw an [`Intersection`](@ref) as a circle, optionally showing the surface normal at the point.
"""
function draw!(scene::Makie.LScene, intersection::Intersection; normal::Bool = false, kwargs...)
    draw!(scene, point(intersection); kwargs...)
    if normal
        draw!(scene, Ray(point(intersection), normal(intersection)); kwargs...)
    end
end

"""
    draw!(scene::Makie.LScene, lines::AbstractVector{Tuple{AbstractVector{T},AbstractVector{T}}}; kwargs...)

Draw a vector of lines.
"""
function draw!(scene::Makie.LScene, lines::AbstractVector{Tuple{P,P}}; kwargs...) where {T<:Real,P<:AbstractVector{T}}
    for l in lines
        draw!(scene, l; kwargs...)
    end
end

"""
    draw!(scene::Makie.LScene, line::Tuple{AbstractVector{T},AbstractVector{T}}; color = :yellow, kwargs...)

Draw a line between two points, `kwargs` are passed to [`Makie.linesegments`](http://makie.juliaplots.org/stable/plotting_functions.html#linesegments).
"""
function draw!(scene::Makie.LScene, line::Tuple{P,P}; color = :yellow, kwargs...) where {T<:Real,P<:AbstractVector{T}}
    Makie.linesegments!(scene, [line[1], line[2]]; kwargs..., color = color)
end

"""
    draw!(s::Makie.LScene, point::AbstractVector{T}; kwargs...)

Draw a single point, `kwargs` are passed to `draw!(scene, points::AbstractVector{AbstractVector{T}})`.
"""
function draw!(s::Makie.LScene, point::AbstractVector{T}; kwargs...) where {T<:Real}
    draw!(s, [point]; kwargs...)
end

"""
    draw!(scene::Makie.LScene, points::AbstractVector{AbstractVector{T}}; markersize = 20, color = :black, kwargs...)

Draw a vector of points.
`kwargs` are passed to [`Makie.scatter`](http://makie.juliaplots.org/stable/plotting_functions.html#scatter).
"""
function draw!(scene::Makie.LScene, points::AbstractVector{P}; markersize = 20, color = :black, kwargs...) where {T<:Real,P<:AbstractVector{T}}
    Makie.scatter!(scene, points, markersize = markersize, color = color, strokewidth = 0; kwargs...)
end

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

function plotOPD!(sys::AxisymmetricOpticalSystem{T}; label = nothing, color = nothing, collimated::Bool = true, axis::Int = 1, samples::Int = 5, wavelength::T = 0.55, sourcepos::SVector{3,T} = SVector{3,T}(0.0, 0.0, 10.0), kwargs...) where {T<:Real}
    sysdiam = semidiameter(sys)
    rays = Vector{OpticalRay{T,3}}(undef, 0)
    positions = Vector{T}(undef, 0)
    if axis == 2
        if collimated
            dir = -sourcepos
            for v in 0:samples
                s = (2 * (v / samples) - 1.0) * sysdiam
                p = SVector{3,T}(0, s, 0)
                push!(positions, s)
                push!(rays, OpticalRay(p - dir, dir, 1.0, wavelength))
            end
        else
            for v in 0:samples
                s = (2 * (v / samples) - 1.0) * sysdiam
                p = SVector{3,T}(0, s, 0)
                dir = p - sourcepos
                push!(positions, s)
                push!(rays, OpticalRay(sourcepos, dir, 1.0, wavelength))
            end
        end
    elseif axis == 1
        if collimated
            dir = -sourcepos
            for v in 0:samples
                s = (2 * (v / samples) - 1.0) * sysdiam
                p = SVector{3,T}(s, 0, 0)
                push!(positions, s)
                push!(rays, OpticalRay(p - dir, dir, 1.0, wavelength))
            end
        else
            for v in 0:samples
                s = (2 * (v / samples) - 1.0) * sysdiam
                p = SVector{3,T}(s, 0, 0)
                dir = p - sourcepos
                push!(positions, s)
                push!(rays, OpticalRay(sourcepos, dir, 1.0, wavelength))
            end
        end
    end
    raygenerator = RayListSource(rays)
    chiefray = OpticalRay(sourcepos, -sourcepos, 1.0, wavelength)
    chiefres = trace(sys, chiefray, test = true)
    xs = Vector{T}()
    opls = Vector{T}()
    for (i, r) in enumerate(raygenerator)
        lt = trace(sys, r, test = true)
        if !(nothing === lt)
            push!(xs, positions[i + 1])
            push!(opls, (pathlength(lt) - pathlength(chiefres)) / (wavelength / 1000)) # wavelength is um while OPD is mm
        end
    end
    p = Plots.plot!(xs, opls, color = nothing === color ? λtoRGB(wavelength) : color, label = label)
    Plots.xlabel!("Relative position on entrance pupil")
    Plots.ylabel!("OPD (waves)")
    return p
end

"""
    spotdiag(sys::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; size = (500, 500), kwargs...)

Plot a spot diagram for an arbitrary [`CSGOpticalSystem`](@ref) and [`OpticalRayGenerator`](@ref).
All rays from `raygenerator` will be traced through `sys` and their intersection location on the detector plotted.

Also `spotdiag!` of the same arguments to add to an existing plot.
"""
function spotdiag(sys::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; kwargs...) where {T<:Real}
    Plots.plot()
    return spotdiag!(sys, raygenerator; kwargs...)
end

"""
    spotdiag(sys::AxisymmetricOpticalSystem{T}; size = (500, 500), hexapolar::Bool = true, collimated::Bool = true, samples::Int = 5, wavelength::T = 0.55, sourceangle::T = zero(T), sourcepos::SVector{3,T} = SVector{3,T}(0.0, 0.0, 10.0), kwargs...)

Plot a spot diagram for an [`AxisymmetricOpticalSystem`](@ref), rays are distributed across the entrance pupil of the system either in a hexapolar of rectangular grid pattern depending on `hexapolar`.

The input rays can be `collimated`, in which case the `sourceangle` parameter determines their direction.
Otherwise rays are treated as coming from a point source at `sourcepos`.

Also `spotdiag!` of the same arguments to add to an existing plot.
"""
function spotdiag(sys::AxisymmetricOpticalSystem{T}; kwargs...) where {T<:Real}
    Plots.plot()
    return spotdiag!(sys; kwargs...)
end

function spotdiag!(sys::AxisymmetricOpticalSystem{T}; color = nothing, hexapolar::Bool = true, collimated::Bool = true, samples::Int = 5, wavelength::T = 0.55, sourceangle::T = zero(T), sourcepos::SVector{3,T} = SVector{3,T}(0.0, 0.0, 10.0), kwargs...) where {T<:Real}
    if hexapolar
        source = HexapolarField(sys, samples = samples, collimated = collimated, wavelength = wavelength, sourcepos = sourcepos, sourceangle = sourceangle)
    else
        source = GridField(sys, samples = samples, collimated = collimated, wavelength = wavelength, sourcepos = sourcepos, sourceangle = sourceangle)
    end
    return spotdiag!(sys.system, source, color = color; kwargs...)
end

function spotdiag!(sys::CSGOpticalSystem{T}, raygenerator::OpticalRayGenerator{T}; color = nothing, label = nothing, size = (500, 500), kwargs...) where {T<:Real}
    xs = Vector{T}()
    ys = Vector{T}()
    det = detector(sys)
    for (i, r) in enumerate(raygenerator)
        lt = trace(sys, r, test = true)
        if !(nothing === lt)
            push!(xs, dot(point(lt) - centroid(det), det.uvec))
            push!(ys, dot(point(lt) - centroid(det), det.vvec))
        end
    end
    # relative to centroid
    cen = sum(xs) / length(xs), sum(ys) / length(ys)
    xs .-= cen[1]
    ys .-= cen[2]
    p = Plots.scatter!(xs, ys, color = (color === nothing) ? λtoRGB(wavelength(raygenerator[0])) : color, label = label, marker = :+, markerstrokewidth = 0, markerstrokealpha = 0, size = size, aspect_ratio = :equal)
    Plots.xlabel!(p, "x (mm)")
    Plots.ylabel!(p, "y (mm)")
    return p
end

spotdiaggrid(sys::AxisymmetricOpticalSystem{T}, raygenerators::Vector{<:OpticalRayGenerator{T}}; kwargs...) where {T<:Real} = spotdiaggrid(sys.system, raygenerators; kwargs...)

function spotdiaggrid(sys::CSGOpticalSystem{T}, raygenerators::Vector{<:OpticalRayGenerator{T}}; colorbynum::Bool = false, kwargs...) where {T<:Real}
    plots = []
    x = Int(floor(sqrt(length(raygenerators))))
    y = length(raygenerators) ÷ x
    size = 2000 / max(x, y)
    xmin = 0
    xmax = 0
    ymin = 0
    ymax = 0
    for (i, r) in enumerate(raygenerators)
        if colorbynum
            p = spotdiag(sys, r; color = indexedcolor2(i), kwargs..., size = (size, size))
        else
            p = spotdiag(sys, r; kwargs..., size = (size, size))
        end
        xmin = min(min(Plots.xlims(p)...), xmin)
        xmax = max(max(Plots.xlims(p)...), xmax)
        ymin = min(min(Plots.ylims(p)...), ymin)
        ymax = max(max(Plots.ylims(p)...), ymax)
        push!(plots, p)
    end
    # set all plots to the same range
    for p in plots
        Plots.xlims!(p, xmin, xmax)
        Plots.ylims!(p, ymin, ymax)
    end
    return Plots.plot(plots..., layout = Plots.grid(x, y), size = (y * size, x * size))
end

"""
    surfacesag(object::Union{CSGTree{T},Surface{T}}, resolution::Tuple{Int,Int}, halfsizes::Tuple{T,T}; offset::T = T(10), position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 10.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0))

Calculates and displays the surface sag of an arbitrary [`Surface`](@ref) or [`CSGTree`](@ref).

Rays are shot in a grid of size defined by `resolution` across a arectangular area defined by `halfsizes`.
This rectangle is centred at `postion` with normal along `direction` and rotation defined by `rotationvec`.
`offset` is subtracted from the sag measurements to provide values relative to the appropriate zero level.
"""
function surfacesag(object::Union{CSGTree{T},Surface{T}}, resolution::Tuple{Int,Int}, halfsizes::Tuple{T,T}; offset::T = T(10), position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 10.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real}
    direction = normalize(direction)
    if abs(dot(rotationvec, direction)) == one(T)
        rotationvec = SVector{3,T}(1.0, 0.0, 0.0)
    end
    uvec = normalize(cross(normalize(rotationvec), direction))
    vvec = normalize(cross(direction, uvec))
    image = Array{T,2}(undef, resolution)
    image .= NaN
    Threads.@threads for x in 0:(resolution[1] - 1)
        for y in 0:(resolution[2] - 1)
            u = 2 * (x / (resolution[1] - 1)) - 1.0
            v = 2 * (y / (resolution[2] - 1)) - 1.0
            pos = position + (halfsizes[1] * u * uvec) + (halfsizes[2] * v * vvec)
            r = Ray(pos, direction)
            int = closestintersection(surfaceintersection(object, r), false)
            if int !== nothing
                int = int::Intersection{T,3}
                sag = α(int) - offset
                image[y + 1, x + 1] = -sag
            end
        end
    end
    p = Plots.heatmap((-halfsizes[1]):(2 * halfsizes[1] / resolution[1]):halfsizes[1], (-halfsizes[2]):(2 * halfsizes[2] / resolution[2]):halfsizes[2], image, c = :jet1, aspect_ratio = :equal, size = (1000, 1000))
    Plots.xlims!(p, -halfsizes[1], halfsizes[1])
    Plots.ylims!(p, -halfsizes[2], halfsizes[2])
    return p
end

"""
    eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::Transform{T} = identitytransform(T))

Visualise the images formed when tracing `assembly` with a human eye for an evenly sampled `sample_points_x × sample_points_y` grid in the angular range of eyeball rotations `-eye_rotation_x:eye_rotation_x` and `-eye_rotation_y:eye_rotation_y` in each dimension respectively.
`resolution` is the size of the detector image (necessarily square).

The eye must be positioned appropriately relative to the system using `eye_transform`, this should transform the eye to the correct position and orientation when at 0 rotation in both dimensions.
By default the eye is directed along the positive z-axis with the vertex of the cornea at the origin.

The result is displayed as a 4D image - the image seen by the eye is shown in 2D as normal with sliders to vary eye rotation in x and y.
The idea being that the whole image should be visible for all rotations in the range.
"""
function eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::Transform{T} = identitytransform(T)) where {T<:Real}
    out_image = zeros(Float32, resolution, resolution, sample_points_y, sample_points_x)
    xrange = sample_points_x > 1 ? range(-eye_rotation_x, stop = eye_rotation_x, length = sample_points_x) : [0.0]
    yrange = sample_points_y > 1 ? range(-eye_rotation_y, stop = eye_rotation_y, length = sample_points_y) : [0.0]

    for (xi, erx) in enumerate(xrange)
        for (yi, ery) in enumerate(yrange)
            print("Position: ($xi, $yi) / ($sample_points_x, $sample_points_y)\r")
            r = OpticSim.rotmatd(T, ery, erx, 0.0)
            ov = OpticSim.rotate(eye_transform, SVector{3,T}(0.0, 0.0, 13.0))
            t = r * ov - ov
            sys = HumanEye.ModelEye(assembly, pupil_radius = pupil_radius, detpixels = resolution, transform = eye_transform * Transform(r, t))
            # Vis.draw(sys)
            # Vis.draw!((eye_transform.translation - ov - SVector(0.0, 20.0, 0.0), eye_transform.translation - ov + SVector(0.0, 20.0, 0.0)))
            # Vis.draw!((eye_transform.translation - ov - SVector(20.0, 0.0, 0.0), eye_transform.translation - ov + SVector(20.0, 0.0, 0.0)))
            # Vis.draw!(eye_transform.translation, markersize = 500.0)
            out_image[:, :, yi, xi] = traceMT(sys, raygen, printprog = false)
        end
    end
    # this is a 4D image - ImageView makes sliders to visualise the dimensions >2!
    Vis.show(out_image)
end

"""
    eyebox_eval_planar(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eyebox::Rectangle{T}, sample_points_x::Int, sample_points_y::Int, vsize::T; pupil_radius::T = T(2.0), resolution::Int = 512)

Visualise the images formed when tracing `assembly` for multiple pupil positions within a planar `eyebox`.
Any angles which are present in a circle radius `pupil_radius` around each sampling point on an even `sample_points_x × sample_points_y` grid on the `eyebox` are added to the sub-image at that grid point.

A paraxial lens focal length 1mm is placed at the eyebox and a detector of size `vsize × vsize`mm placed 1mm behind it.
The normal of the detector rectangle should point towards the system (and away from the fake detector).
Any rays which miss the detector are ignored.

The pupil is always fully contained in the eyebox, i.e., the extreme sample position in u would be `eyebox.halfsizeu - pupil_radius`, for example.

The result is displayed as a 4D image - each sub-image is shown as normal with sliders to vary eye rotation in x and y.
The idea being that the whole FoV should be visible for all rotations in the range.
"""
function eyebox_eval_planar(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eyebox::Rectangle{T}, sample_points_x::Int, sample_points_y::Int, vsize::T; pupil_radius::T = T(2.0), resolution::Int = 512) where {T<:Real}
    sys = CSGOpticalSystem(assembly, eyebox, 1, 1)
    hits = tracehitsMT(sys, raygen)
    println("Calculating results...")
    det = Rectangle(vsize, vsize, normal(eyebox), centroid(eyebox) - normal(eyebox))
    if length(hits) > 0
        out_image = zeros(Float32, resolution, resolution, sample_points_y, sample_points_x)
        for res in hits
            dir = direction(ray(res))
            u, v = uv(res)
            x = u * eyebox.halfsizeu
            y = v * eyebox.halfsizev

            raydirection = normalize((centroid(eyebox) - point(res)) + direction(ray(res)) / dot(-normal(eyebox), direction(ray(res))))
            focal_plane_hit = surfaceintersection(det, Ray(point(res), raydirection))
            if focal_plane_hit isa EmptyInterval
                continue
            end
            px, py = OpticSim.uvtopix(det, uv(OpticSim.halfspaceintersection(focal_plane_hit)), (resolution, resolution))

            hsu = eyebox.halfsizeu - pupil_radius
            hsv = eyebox.halfsizev - pupil_radius
            xrange = sample_points_x > 1 ? range(-hsu, stop = hsu, length = sample_points_x) : [0.0]
            yrange = sample_points_y > 1 ? range(-hsv, stop = hsv, length = sample_points_y) : [0.0]

            for (xi, posx) in enumerate(xrange)
                for (yi, posy) in enumerate(yrange)
                    if (x - posx)^2 + (y - posy)^2 < pupil_radius^2
                        out_image[py, px, yi, xi] += sourcepower(ray(res))
                    end
                end
            end
        end
        # this is a 4D image - ImageView makes sliders to visualise the dimensions >2!
        Vis.show(out_image)
    end
end


text:
#########################################################################
#
# rdiff differentiation function
# copy from c:/Users/awf/.julia/v0.4/ReverseDiffSource
#
#########################################################################

#########################################################################
#
#   Expression to graph conversion
#
#########################################################################
#
###########  Parameterized type to ease AST exploration  ############
#  type ExH{H}
#    head::Symbol
#    args::Vector
#    typ::Any
#  end
#  toExH(ex::Expr) = ExH{ex.head}(ex.head, ex.args, ex.typ)
#  toExpr(ex::ExH) = Expr(ex.head, ex.args...)
#
#  typealias ExEqual    ExH{:(=)}
#  typealias ExDColon   ExH{:(::)}
#  typealias ExColon    ExH{:(:)}
#  typealias ExPEqual   ExH{:(+=)}
#  typealias ExMEqual   ExH{:(-=)}
#  typealias ExTEqual   ExH{:(*=)}
#  typealias ExTrans    ExH{symbol("'")}
#  typealias ExCall     ExH{:call}
#  typealias ExBlock    ExH{:block}
#  typealias ExLine     ExH{:line}
#  typealias ExVcat     ExH{:vcat}
#  typealias ExVect     ExH{:vect}
#  typealias ExCell1d   ExH{:cell1d}
#  typealias ExCell     ExH{:cell1d}
#  typealias ExFor      ExH{:for}
#  typealias ExRef      ExH{:ref}
#  typealias ExIf       ExH{:if}
#  typealias ExComp     ExH{:comparison}
#  typealias ExDot      ExH{:.}
#  typealias ExTuple    ExH{:tuple}
#  typealias ExReturn   ExH{:return}
#  typealias ExBody     ExH{:body}
#  typealias ExQuote    ExH{:QuoteNode}
#

#  s     : expression to convert
#  svars : vars set since the toplevel graph (helps separate globals / locals)
function awftograph(s, evalmod=Main, svars=Any[])

    explore(ex::Any)       = error("[awftograph] unmanaged type $ex ($(typeof(ex)))")
    explore(ex::Expr)      = explore(toExH(ex))
    explore(ex::ExH)       = error("[awftograph] unmanaged expr type $(ex.head) in ($ex)")

    explore(ex::ExLine)         = nothing     # remove line info
    explore(ex::LineNumberNode) = nothing     # remove line info
    explore(ex::QuoteNode)      = addnode!(g, NConst(ex.value))  # consider as constant

    explore(ex::ExReturn)  = explore(ex.args[1]) # focus on returned statement

    explore(ex::ExVcat)    = explore(Expr(:call, :vcat, ex.args...) )  # translate to vcat() call, and explore
    explore(ex::ExVect)    = explore(Expr(:call, :vcat, ex.args...) )  # translate to vcat() call, and explore
    explore(ex::ExCell1d)  = explore(Expr(:call, :(Base.cell_1d), ex.args...) )  # translate to cell_1d() call, and explore
    explore(ex::ExTrans)   = explore(Expr(:call, :transpose, ex.args[1]) )  # translate to transpose() and explore
    explore(ex::ExColon)   = explore(Expr(:call, :colon, ex.args...) )  # translate to colon() and explore
    explore(ex::ExTuple)   = explore(Expr(:call, :tuple, ex.args...) )  # translate to tuple() and explore

    explore(ex::ExPEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :+, args[1], args[2])) ) )
    explore(ex::ExMEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :-, args[1], args[2])) ) )
    explore(ex::ExTEqual)  = (args = ex.args ; explore( Expr(:(=), args[1], Expr(:call, :*, args[1], args[2])) ) )

    explore(ex::Real)      = addnode!(g, NConst(ex))

    explore(ex::ExBlock)   = map( explore, ex.args )[end]
    explore(ex::ExBody)    = map( explore, ex.args )[end]

    explore(ex::ExComp)    = addnode!(g, NComp(ex.args[2], [explore(ex.args[1]), explore(ex.args[3])]))

    # explore(ex::ExDot)     = addnode!(g, NDot(ex.args[2],     [ explore(ex.args[1]) ]))
    explore(ex::ExDot)     = explore(Expr(:call, :getfield, ex.args...))
    explore(ex::ExRef)     = explore(Expr(:call, :getindex, ex.args...))

    function explore(ex::Symbol)
        hassym(g.seti, ex)       && return getnode(g.seti, ex)
        hassym(g.exti, ex)       && return getnode(g.exti, ex)

        nn = addnode!(g, NExt(ex))    # create external node for this var
        g.exti[nn] = ex
        return nn
    end

    function explore(ex::ExCall)
        sf = ex.args[1]
        if sf == :getindex
            nv = explore(ex.args[2])
            ps = indexspec(nv, ex.args[3:end])
            return addnode!(g, NRef(:getidx, vcat([nv], ps)))

        elseif sf == :setindex!
            isa(ex.args[2], Symbol) ||
                error("[awftograph] setindex! only allowed on variables, $(ex.args[2]) found")

            nv  = explore(ex.args[2]) # node whose subpart is assigned
            ps  = indexspec(nv, ex.args[4:end])
            rhn = addnode!(g, NSRef(:setidx,
                                    [ nv,                               # var modified in pos #1
                                      explore(ex.args[3]),              # value affected in pos #2
                                      ps...] ))                         # dims

            rhn.precedence = filter(n -> nv in n.parents && n != rhn, g.nodes)
            g.seti[rhn] = ex.args[2]

            return rhn

        elseif sf == :getfield
            return addnode!(g, NDot(ex.args[3], [ explore(ex.args[2]) ]))

        elseif sf == :setfield!
            isa(ex.args[2], Symbol) ||
                error("[awftograph] setfield! only allowed on variables, $(ex.args[2]) found")

            nv  = explore(ex.args[2]) # node whose subpart is assigned
            rhn = addnode!(g, NSDot(ex.args[3],
                                    [ nv,                               # var modified in pos #1
                                      explore(ex.args[4])]))            # value affected in pos #2

            rhn.precedence = filter(n -> nv in n.parents && n != rhn, g.nodes)
            g.seti[rhn] = ex.args[2]

            return rhn

        else
            return addnode!(g, NCall(  :call,
                                        map(explore, ex.args[1:end]) ))
        end
    end

    function explore(ex::ExEqual)
        lhs = ex.args[1]

        if isSymbol(lhs)  # x = ....
            lhss = lhs

            # set before ? call explore
            if lhss in union(svars, collect(syms(g.seti)))
                vn = explore(lhss)
                rhn  = addnode!(g, NSRef(:setidx,
                                         [ vn,                     # var modified in pos #1
                                           explore(ex.args[2]) ])) # value affected in pos #2
                rhn.precedence = filter(n -> vn in n.parents && n != rhn, g.nodes)

            else # never set before ? assume it is created here
                rhn = explore(ex.args[2])

                # we test if RHS has already a symbol
                # if it does, to avoid loosing it, we create an NIn node
                if hasnode(g.seti, rhn)
                    rhn = addnode!(g, NIn(lhss, [rhn]))
                end
            end

        elseif isRef(lhs)   # x[i] = ....
            lhss = lhs.args[1]
            rhn = explore( Expr(:call, :setindex!, lhss, ex.args[2], lhs.args[2:end]...) )
            # vn = explore(lhss) # node whose subpart is assigned
            # rhn  = addnode!(g, NSRef(:setidx, [ vn,    # var modified in pos #1
            #                                     explore(ex.args[2]), # value in pos #2
            #                                     map(explore, lhs.args[2:end])] ))  # indexing starting at #3
            # rhn.precedence = filter(n -> vn in n.parents && n != rhn, g.nodes)

        elseif isDot(lhs)   # x.field = ....
            lhss = lhs.args[1]
            rhn = explore( Expr(:call, :setfield!, lhss, lhs.args[2], ex.args[2]) )

            # vn = explore(lhss) # node whose subpart is assigned
            # rhn  = addnode!(g, NSDot(lhs.args[2], [ vn, explore(ex.args[2])] ))
            # rhn.precedence = filter(n -> vn in n.parents && n != rhn, g.nodes)

        else
            error("[awftograph] $(toExpr(ex)) not allowed on LHS of assigment")
        end

        g.seti[rhn] = lhss

        return nothing
    end

    function explore(ex::ExFor)
        is = ex.args[1].args[1]
        isa(is, Symbol) ||
            error("[awftograph] for loop using several indexes : $is ")

        # explore the index range
        nir = explore(ex.args[1].args[2])

        # explore the for block as a separate graph
        nsvars = union(svars, collect(syms(g.seti)))
        g2 = awftograph(ex.args[2], evalmod, nsvars)

        # create "for" node
        nf = addnode!(g, NFor( Any[ is, g2 ] ))
        nf.parents = [nir]  # first parent is indexing range fo the loop

        # create onodes (node in parent graph) for each exti
        for (k, sym) in g2.exti.kv
            sym==is  && continue # loop index should be excluded
            pn = explore(sym)  # look in setmap, externals or create it
            g2.exto[pn] = sym
            push!(nf.parents, pn) # mark as parent of for loop
        end

        # create onodes and 'Nin' nodes for each seti
        #  will be restricted to variables that are defined in parent
        #   (others are assumed to be local to the loop)
        for (k, sym) in g2.seti.kv
            if sym in nsvars && sym != is # only for variables set in parent scope
                pn = explore(sym)                   # create node if needed
                rn = addnode!(g, NIn(sym, [nf]))    # exit node for this var in this graph
                g.seti[rn] = sym                    # signal we're setting the var
                g2.seto[rn] = sym

                append!(nf.precedence, filter(n -> pn in n.parents && n != nf, g.nodes))

                # create corresponding exti if it's not already done
                if !hassym(g2.exto, sym)
                    g2.exto[pn] = sym
                    push!(nf.parents, pn) # mark as parent of for loop
                end
            end
        end
    end

    #### translates ':' and 'end' special symbols in getindex / setindex!
    function indexspec(nv, as)
        p  = ExNode[]
        for (i,na) in enumerate(as)
            if length(as)==1 # single dimension
                ns = addgraph!(:( length(x) ), g,  Dict(:x => nv) )
            else # several dimensions
                ns = addgraph!(:( size(x, $i) ), g,  Dict(:x => nv) )
            end

            na==:(:) && (na = Expr(:(:), 1, :end) )  # replace (:) with (1:end)

            # explore the dimension expression
            nsvars = union(svars, collect(syms(g.seti)))
            ng = awftograph(na, evalmod, nsvars)

            # find mappings for onodes, including :end
            vmap = Dict()
            for (k, sym) in ng.exti.kv
                vmap[sym] = sym == :end ? ns : explore(sym)
            end

            nd = addgraph!(ng, g, vmap)
            push!(p, nd)
        end
        p
    end

    #  top level graph
    g = ExGraph()

    exitnode = explore(s)
    # exitnode = nothing if only variable assigments in expression
    #          = ExNode of last calc otherwise

    # id is 'nothing' for unnassigned last statement
    exitnode!=nothing && ( g.seti[exitnode] = nothing )

    # Resolve external symbols that are Functions, DataTypes or Modules
    # and turn them into constants
    for en in filter(n -> isa(n, NExt) & !in(n.main, svars) , keys(g.exti))
        if isdefined(evalmod, en.main)  # is it defined
            tv = evalmod.eval(en.main)
            isa(tv, TypeConstructor) && error("[awftograph] TypeConstructors not supported: $ex $(tv), use DataTypes")
            if isa(tv, DataType) || isa(tv, Module) || isa(tv, Function)
                delete!(g.exti, en)
                nc = addnode!(g, NConst( tv ))
                fusenodes(g, nc, en)
            end
        end
    end

    g
end

##########  function version   ##############

function awfrdiff(f::Function, sig0::Tuple; order::Int=1, evalmod=Main, debug=false, allorders=true)
    display("HELLO")
    sig = map( typeof, sig0 )
    fs = methods(f, sig)
    length(fs) == 0 && error("no function '$f' found for signature $sig")
    length(fs) > 1  && error("several functions $f found for signature $sig")  # is that possible ?

    fdef  = fs[1].func.code
    fcode = Base.uncompressed_ast(fdef)
    fargs = fcode.args[1]  # function parameters

    cargs = [ (fargs[i], sig0[i]) for i in 1:length(sig0) ]
    dex = awfrdiff(fcode.args[3]; order=order, evalmod=evalmod, debug=debug,
                allorders=allorders, cargs...)

    # Note : new function is created in the same module as original function
    myf = fdef.module.eval( :( $(Expr(:tuple, fargs...)) -> $dex ) )
end

######### expression version   ################
# TODO : break this huge function in smaller blocks

function awfrdiff(ex; outsym=nothing, order::Int=1, evalmod=Main, debug=false, allorders=true, params...)

    length(params) >= 1 || error("There should be at least one parameter specified, none found")

    order <= 1 ||
    length(params) == 1 || error("Only one param allowed for order >= 2")

    order <= 1 ||
    isa(params[1][2], Vector) ||
    isa(params[1][2], Real)   || error("Param should be a real or vector for order >= 2")

    paramsym    = Symbol[ e[1] for e in params]
    paramvalues = [ e[2] for e in params]
    parval      = Dict(zip(paramsym, paramvalues))

    g = awftograph(ex, evalmod)

    hassym(g.seti, outsym) ||
    error("can't find output var $( outsym==nothing ? "" : outsym)")

    # reduce to variable of interest
    g.seti = NSMap([getnode(g.seti, outsym)], [ outsym ])

    g |> splitnary! |> prune! |> simplify!
    calc!(g, params=parval, emod=evalmod)

    ov = getnode(g.seti, outsym).val
    isa(ov, Real) || error("output var should be a Real, $(typeof(ov)) found")

    voi = Any[ outsym ]

    if order == 1
        dg = reversegraph(g, getnode(g.seti, outsym), paramsym)
        append!(g.nodes, dg.nodes)

        for p in paramsym
            nn = getnode(dg.seti, dprefix(p))  # find the exit node of deriv of p
            ns = newvar("_dv")
            g.seti[nn] = ns
            push!(voi, ns)
        end

        g |> splitnary! |> prune! |> simplify!

    elseif order > 1 && isa(paramvalues[1], Real)
        for i in 1:order
            dg = reversegraph(g, getnode(g.seti, voi[i]), paramsym)
            append!(g.nodes, dg.nodes)
            nn = collect(nodes(dg.seti))[1]  # only a single node produced
            ns = newvar("_dv")
            g.seti[nn] = ns
            push!(voi, ns)

            g |> splitnary! |> prune! |> simplify!

            calc!(g, params=parval, emod=evalmod)
        end

    elseif order > 1 && isa(paramvalues[1], Vector)
        # do first order as usual
        dg = reversegraph(g, getnode(g.seti, outsym), paramsym)
        append!(g.nodes, dg.nodes)
        ns = newvar(:_dv)
        g.seti[ collect(nodes(dg.seti))[1] ] = ns
        push!(voi, ns)

        g |> splitnary! |> prune! |> simplify!

        # now order 2 to n
        for i in 2:order
            # launch derivation on a single value of the preceding
            #   derivation vector
            no = getnode(g.seti, voi[i])
            si = newvar(:_idx)
            ni = addnode!(g, NExt(si))
            ns = addnode!(g, NRef(:getidx, [ no, ni ]))

            calc!(g, params=Dict(zip([paramsym; si], [paramvalues; 1])), emod=evalmod)
            dg = reversegraph(g, ns, paramsym)

            #### We will now wrap dg in a loop scanning all the elements of 'no'
            # first create ext nodes to make dg a complete subgraph
            dg2 = ExNode[]
            nmap = Dict()
            for n in dg.nodes  # n = dg.nodes[2]
                for (j, np) in enumerate(n.parents)  # j,np = 1, n.parents[1]
                    if haskey(nmap, np) # already remapped
                        n.parents[j] = nmap[np]

                    elseif np == ni # it's the loop index
                        nn = NExt(si)
                        push!(dg2, nn)
                        dg.exti[nn] = si
                        n.parents[j] = nn
                        nmap[np] = nn

                    elseif np == ns # it's the selected element of the deriv vector
                        # create 'no' ref if needed
                        if !haskey(nmap, no)
                            sn = newvar()
                            nn = NExt(sn)
                            push!(dg2, nn)
                            dg.exti[nn] = sn
                            dg.exto[no] = sn
                            nmap[no] = nn
                        end

                        nn = NRef(:getidx, [ nmap[no], nmap[ni] ])
                        push!(dg2, nn)
                        nmap[ns] = nn

                    elseif !(np in dg.nodes) # it's not in dg (but in g)
                        sn = newvar()
                        nn = NExt(sn)
                        push!(dg2, nn)
                        dg.exti[nn] = sn
                        dg.exto[np] = sn
                        n.parents[j] = nn
                        nmap[np] = nn

                    end
                end

                # update onodes in for loops
                if isa(n, NFor)
                    g2 = n.main[2]
                    for (o,s) in g2.exto
                        if haskey(nmap, o)
                            g2.exto[ nmap[o] ] = s  # replace
                        end
                    end
                end
            end
            append!(dg.nodes, dg2)
            # dg |> prune! |> simplify!

            # create for loop node
            nf = addnode!(g, NFor(Any[ si, dg ] ) )

            # create param size node
            nsz = addgraph!( :( length( x ) ), g,  Dict( :x => getnode(g.exti, paramsym[1]) ) )

            # create (n-1)th derivative size node
            ndsz = addgraph!( :( sz ^ $(i-1) ), g, Dict( :sz => nsz ) )

            # create index range node
            nid = addgraph!( :( 1:dsz ),  g,  Dict( :dsz => ndsz ) )
            push!(nf.parents, nid)

            # pass size node inside subgraph
            sst = newvar()
            inst = addnode!(dg, NExt(sst))
            dg.exti[inst] = sst
            dg.exto[nsz]  = sst
            push!(nf.parents, nsz)

            # create result node (alloc in parent graph)
            nsa = addgraph!( :( zeros( $( Expr(:tuple, [:sz for j in 1:i]...) ) ) ),
                            g,  Dict( :sz => nsz ) )
            ssa = newvar()
            insa = addnode!(dg, NExt(ssa))
            dg.exti[insa] = ssa
            dg.exto[nsa]  = ssa
            push!(nf.parents, nsa)

            # create result node update (in subgraph)
            nres = addgraph!( :( res[ ((sidx-1)*st+1):(sidx*st) ] = dx ; res ), dg,
                              Dict(:res  => insa,
                                          :sidx => nmap[ni],
                                          :st   => inst,
                                          :dx   => collect(dg.seti)[1][1] ) )
            dg.seti = NSMap([nres], [ssa])

            # create exit node for result
            nex = addnode!(g, NIn(ssa, [nf]))
            dg.seto = NSMap([nex], [ssa])

            # update parents of for loop
            append!( nf.parents, setdiff(collect( nodes(dg.exto)), nf.parents[2:end]) )

            ns = newvar(:_dv)
            g.seti[nex] = ns
            push!(voi, ns)

            g |> splitnary! |> prune! |> simplify!
        end

    end

    if !allorders  # only keep the last derivative
        voi = [voi[end]]
    end

    if length(voi) > 1  # create tuple if multiple variables
        voin = map( s -> getnode(g.seti, s), voi )
        nf = addnode!(g, NConst(tuple))
        exitnode = addnode!(g, NCall(:call, [nf, voin...]))
    else
        exitnode = getnode(g.seti, voi[1])
    end
    g.seti = NSMap( [exitnode], [nothing])  # make this the only exitnode of interest

    g |> splitnary! |> prune! |> simplify!

    resetvar()
    debug ? g : tocode(g)
end

ext: jl

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
using OffsetArrays,Test
using LegendrePolynomials

@testset "Pl array" begin
	x = 2rand() - 1
    P = Pl(x,lmax=5)
    @test P[0]≈1
    @test P[1]≈x
    @test P[2]≈(3x^2-1)/2
    @test P[3]≈(5x^3-3x)/2
    @test P[4]≈(35x^4-30x^2+3)/8
    @test P[5]≈(63x^5-70x^3+15x)/8
end

@testset "dPl array" begin
	x = 2rand() - 1
    P = dPl(x,lmax=5)
    @test P[0]≈0
    @test P[1]≈1
    @test P[2]≈3x
    @test P[3]≈(-3+15x^2)/2
    @test P[4]≈1/8*(-60x + 140x^3)
    @test P[5]≈1/8*(15 - 210x^2 + 315x^4)
end

@testset "d2Pl array" begin
	x = 2rand() - 1
    P = d2Pl(x,lmax=5)
    @test P[0]≈0
    @test P[1]≈0
    @test P[2]≈3
    @test P[3]≈15x
    @test P[4]≈1/8*(-60 + 420x^2)
    @test P[5]≈1/8*(-420x + 1260x^3)
end

@testset "allmodes vs one" begin
    x = 2rand() - 1
    lmax = 10
    P,dP,d2P = Pl_dPl_d2Pl(x,lmax=lmax)
    @testset "Pl" begin
        for l=0:lmax
            @test P[l]==Pl(x,l)
        end
    end
    @testset "dPl" begin
        for l=0:lmax
            @test dP[l]==dPl(x,l)
        end
    end
    @testset "d2Pl" begin
        for l=0:lmax
            @test d2P[l]==d2Pl(x,l)
        end
    end
end

@testset "boundary" begin

    lmax = 10
    P1,dP1 = Pl_dPl(1,lmax=lmax)
    Pm1,dPm1 = Pl_dPl(-1,lmax=lmax)

    @testset "Pl" begin
        for l in axes(P1,1)
            @test P1[l] == 1
        end
        for l in axes(Pm1,1)
            @test Pm1[l] == (-1)^l
        end
    end
    @testset "dPl" begin
        for l in axes(dP1,1)
            @test dP1[l] == l*(l+1)/2
        end
        for l in axes(dPm1,1)
            @test dPm1[l] == (-1)^(l+1) * l*(l+1)/2
        end
    end
end

@testset "Askey-Gasper" begin
    x=0.5
    for lmax=1:100:1_000
        @test sum(Pl(x,lmax=lmax)) >= 0
    end
end


text:
using Pkg: develop, instantiate, PackageSpec
develop(PackageSpec(path=pwd()))
using LightQuery

instantiate()
using Documenter: deploydocs, makedocs

makedocs(sitename = "LightQuery.jl", modules = [LightQuery], doctest = false)
deploydocs(repo = "github.com/bramtayl/LightQuery.jl.git")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module ResourceManagementClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_DeploymentOperationsApi.jl")
include("api_DeploymentsApi.jl")
include("api_OperationsApi.jl")
include("api_ProvidersApi.jl")
include("api_ResourceGroupsApi.jl")
include("api_ResourcesApi.jl")
include("api_TagsApi.jl")

# export models
export convert, Alias
export convert, AliasPath
export convert, AliasPattern
export convert, BasicDependency
export convert, CloudError
export convert, DebugSetting
export convert, Dependency
export convert, Deployment
export convert, DeploymentExportResult
export convert, DeploymentExtended
export convert, DeploymentExtendedFilter
export convert, DeploymentListResult
export convert, DeploymentOperation
export convert, DeploymentOperationProperties
export convert, DeploymentOperationsListResult
export convert, DeploymentProperties
export convert, DeploymentPropertiesExtended
export convert, DeploymentValidateResult
export convert, DeploymentWhatIf
export convert, DeploymentWhatIfProperties
export convert, DeploymentWhatIfSettings
export convert, ErrorAdditionalInfo
export convert, ErrorResponse
export convert, ExportTemplateRequest
export convert, GenericResource
export convert, GenericResourceExpanded
export convert, GenericResourceFilter
export convert, HttpMessage
export convert, Identity
export convert, IdentityUserAssignedIdentities
export convert, OnErrorDeployment
export convert, OnErrorDeploymentExtended
export convert, Operation
export convert, OperationDisplay
export convert, OperationListResult
export convert, ParametersLink
export convert, Plan
export convert, Provider
export convert, ProviderListResult
export convert, ProviderResourceType
export convert, Resource
export convert, ResourceGroup
export convert, ResourceGroupExportResult
export convert, ResourceGroupFilter
export convert, ResourceGroupListResult
export convert, ResourceGroupPatchable
export convert, ResourceGroupProperties
export convert, ResourceListResult
export convert, ResourceProviderOperationDisplayProperties
export convert, ResourceReference
export convert, ResourcesMoveInfo
export convert, ScopedDeployment
export convert, ScopedDeploymentWhatIf
export convert, Sku
export convert, StatusMessage
export convert, SubResource
export convert, TagCount
export convert, TagDetails
export convert, TagValue
export convert, Tags
export convert, TagsListResult
export convert, TagsPatchResource
export convert, TagsResource
export convert, TargetResource
export convert, TemplateHashResult
export convert, TemplateLink
export convert, WhatIfChange
export convert, WhatIfOperationProperties
export convert, WhatIfOperationResult
export convert, WhatIfPropertyChange

# export operations
export convert, DeploymentOperationsApi, DeploymentsApi, OperationsApi, ProvidersApi, ResourceGroupsApi, ResourcesApi, TagsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
mapch(f, x, y) = mapc(f, x, y)
mapch(f, x::C, y::C) where C <: Union{HSV, HSL, HSI} =
    C(f(Hue(x), Hue(y)), f(x.s, y.s), f(comp3(x), comp3(y)))
mapch(f, x::C, y::C) where C <: Union{LCHab, LCHuv} =
    C(f(x.l, y.l), f(x.c, y.c), f(Hue(x), Hue(y)))

mapca(fc, a, x::C) where C <: TransparentColorN{2} = C(fc(comp1(x)), a)
mapca(fc, a, x::C) where C <: TransparentColorN{3} = C(fc(comp1(x)), fc(comp2(x)), a)
mapca(fc, a, x::C) where C <: TransparentColorN{4} = C(fc(comp1(x)), fc(comp2(x)), fc(comp3(x)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{2} =
    C(fc(comp1(x), comp1(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{3} =
    C(fc(comp1(x), comp1(y)), fc(comp2(x), comp2(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4} =
    C(fc(comp1(x), comp1(y)), fc(comp2(x), comp2(y)), fc(comp3(x), comp3(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4, <:Union{HSV, HSL, HSI}} =
    C(fc(Hue(x), Hue(y)), fc(x.s, y.s), fc(comp3(x), comp3(y)), a)
mapca(fc, a, x::C, y::C) where C <: TransparentColorN{4, <:Union{LCHab, LCHuv}} =
    C(fc(x.l, y.l), fc(x.c, y.c), fc(Hue(x), Hue(y)))

# complement
_n(v::T) where T = oneunit(T) - v

# linear interpolation
_w(v1::T, v2::T, w) where T = _w(v1, _n(w), v2, w)

function _w(h1::Hue{T}, h2::Hue{T}, w) where T
    h1a = ifelse(h1.isgray, h2.angle, h1.angle)
    h2a = ifelse(h2.isgray, h1.angle, h2.angle)
    d0 = h2a - h1a
    d = ifelse(abs(d0) > T(180), d0 - copysign(T(360), d0), d0)
    a = muladd(d, w, h1a)
    convert(T, ifelse(a > T(360), a - T(360), ifelse(a < T(0), a + T(360), a)))
end

_w(v1::T, w1, v2::T, w2) where T = convert(T, muladd(w1, v1, w2 * v2))

function _w(v1::T, w1::T, v2::T, w2::T) where T <: Union{N0f8, N0f16}
    r1, rw1 = reinterpret(v1), reinterpret(w1)
    r2, rw2 = reinterpret(v2), reinterpret(w2)
    m = widemul(r1, rw1) + widemul(r2, rw2)
    if T === N0f8
        z = (m + ((m + 0x80) >> 0x8) + 0x80) >> 0x8
    else
        z = (m + ((m + 0x8000) >> 0x10) + 0x8000) >> 0x10
    end
    reinterpret(T, unsafe_trunc(typeof(r1), z))
end

function _w(h1::Hue{T}, w1, h2::Hue{T}, w2) where T
    w = w1 + w2
    w == oneunit(w) && return _w(h1, h2, w2)
    w == zero(w) && return zero(T)
    _w(h1, h2, w2 / w) * w
end

_w_safe(v1::T, w1, v2::T, w2) where T = min(oneunit(T), _w(v1, w1, v2, w2))
_w_safe(v1::T, w1, v2::T, w2) where T <: FixedPoint =
    convert(T, min(float(typemax(T)), muladd(w1, float(v1), w2 * float(v2))))

_comp(op::CompositeOperation{:clear}, c1, c2) = mapc(v1 -> zero(v1), c1)

_comp(op::CompositeOperation{:copy}, c1, c2) = c2

_comp(op::CompositeOperation{:destination}, c1, c2) = c1

function _comp(op::typeof(CompositeSourceOver), c1, c2)
    k1 = mul(alpha(c1), _n(alpha(c2)))
    k2 = alpha(c2)
    a = k1 + k2
    k1a = a == zero(a) ? a : k1 / a
    k2a = a == zero(a) ? a : oneunit(a) - k1a
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a, c1, c2)
end

_comp(op::typeof(CompositeDestinationOver), c1, c2) = _comp(CompositeSourceOver, c2, c1)

function _comp(op::typeof(CompositeSourceIn), c1, c2)
    a = mul(alpha(c1), alpha(c2))
    ifelse(a == zero(a), mapca(v -> a, a, c2), mapca(v -> v, a, c2))
end

_comp(op::typeof(CompositeDestinationIn), c1, c2) = _comp(CompositeSourceIn, c2, c1)

function _comp(op::typeof(CompositeSourceOut), c1, c2)
    a = mul(_n(alpha(c1)), alpha(c2))
    ifelse(a == zero(a), mapca(v -> a, a, c2), mapca(v -> v, a, c2))
end

_comp(op::typeof(CompositeDestinationOut), c1, c2) = _comp(CompositeSourceOut, c2, c1)

function _comp(op::typeof(CompositeSourceAtop), c1, c2)
    a1, a2 = alpha(c1), alpha(c2)
    k1a = a1 == zero(a1) ? a1 : _n(a2)
    k2a = a1 == zero(a1) ? a1 : a2
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a1, c1, c2)
end

_comp(op::typeof(CompositeDestinationAtop), c1, c2) = _comp(CompositeSourceAtop, c2, c1)

function _comp(op::CompositeOperation{:xor}, c1, c2)
    m = mul(alpha(c1), alpha(c2))
    k1, k2 = alpha(c1) - m, alpha(c2) - m
    a = k1 + k2
    k1a = a == zero(a) ? a : k1 / a
    k2a = a == zero(a) ? a : oneunit(a) - k1a
    mapca((v1, v2) -> _w(v1, k1a, v2, k2a), a, c1, c2)
end

function _comp(op::CompositeOperation{:lighter}, c1, c2)
    k1, k2 = float(alpha(c1)), float(alpha(c2))
    a = k1 + k2
    ac = min(oneunit(a), a)
    k1a = a == zero(a) ? a : k1 / ac
    k2a = a == zero(a) ? a : k2 / ac
    mapca((v1, v2) -> _w_safe(v1, k1a, v2, k2a), ac, c1, c2)
end

mix_alpha(opacity::Nothing, a) = a
mix_alpha(opacity::Real, a) = opacity * a

@inline (mode::BlendMode)(c1, c2; opacity=nothing, op=CompositeSourceOver) =
    _blend_c(mode, c1, c2, opacity, op)

@inline (op::CompositeOperation)(c1, c2; opacity=nothing, mode=BlendNormal) =
    _blend_c(mode, c1, c2, opacity, op)

"""
    blend(c1, c2; mode=BlendNormal, opacity=1, op=CompositeSourceOver)

Create the mixed color of two colors `c1` and `c2`. The `c1` means the backdrop
color and the `c2` means the source color.

`mode` specifies the blend mode, e.g. [`BlendMultiply`](@ref).

`opacity` modifies the source (i.e. `c2`-side) alpha by multiplication.

`op` specifies the composite operations, e.g. [`CompositeSourceAtop`](@ref).

The return type is the same as `c1`.

# Examples
```jldoctest
julia> blend(RGB(1, 0.5, 0), RGB(0, 0.5, 1), mode=BlendLighten)
RGB{Float64}(1.0,0.5,1.0)
```
"""
@inline blend(c1, c2; mode::BlendMode=BlendNormal, opacity=nothing, op=CompositeSourceOver) =
    _blend_c(mode, c1, c2, opacity, op)

@inline _blend_c(mode::BlendMode, c1::Color, c2::Color, opacity, op) =
    _blend_cc(mode, c1, convert(typeof(c1), c2), opacity, op)

@inline _blend_c(mode::BlendMode, c1::Color, c2::TransparentColor, opacity, op) =
    _blend_cc(mode, c1, convert(typeof(c1), color(c2)), mix_alpha(opacity, alpha(c2)), op)

@inline _blend_c(mode::BlendMode, c1::TransparentColor, c2::Color, opacity, op) =
    _blend_tc(mode, c1, convert(color_type(c1), c2), opacity, op)

@inline _blend_c(mode::BlendMode, c1::TransparentColor, c2::TransparentColor, opacity, op) =
    _blend_tc(mode, c1, convert(color_type(c1), color(c2)), mix_alpha(opacity, alpha(c2)), op)

# without opacity
@inline _blend_cc(mode::BlendMode, c1, c2, ::Nothing, ::OpaqueInSourceOutOperation) =
    _blend(mode, c1, c2)

_blend_cc(mode::BlendMode, c1, c2, ::Nothing, ::OpaqueInDestinationOutOperation) = c1

@inline _blend_cc(mode::BlendMode, c1, c2, ::Nothing, op::CompositeOperation{:lighter}) =
    mapc((v1, v2) -> _w_safe(v1, oneunit(v1), v2, oneunit(v2)), c1, _blend(mode, c1, c2))

@inline function _blend_tc(mode::BlendMode, c1, c2, ::Nothing, op)
    cm = mapch((v1, v2) -> _w(v1, v2, alpha(c1)), c2, _blend(mode, color(c1), c2))
    _comp(op, c1, typeof(c1)(cm))
end

# with opacity
@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeSourceOver)) =
    mapch((v1, v2) -> _w(v1, v2, opacity), c1, _blend(mode, c1, c2))

@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeSourceAtop)) =
    _blend_cc(mode, c1, c2, opacity, CompositeSourceOver)

_blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeDestination)) = c1

_blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::typeof(CompositeDestinationOver)) = c1

@inline _blend_cc(mode::BlendMode, c1, c2, opacity::Real, ::CompositeOperation{:lighter}) =
    mapc((v1, v2) -> _w_safe(v1, oneunit(opacity), v2, opacity), c1, _blend(mode, c1, c2))

@inline function _blend_tc(mode::BlendMode, c1, c2, opacity::Real, op)
    cm = mapch((v1, v2) -> _w(v1, v2, alpha(c1)), c2, _blend(mode, color(c1), c2))
    _comp(op, c1, typeof(c1)(cm, opacity))
end

_blend(::BlendMode{:normal}, c1::C, c2::C) where C <: Color = c2

_blend(::BlendMode{:multiply}, c1::C, c2::C) where C <: Color = mapch(mul, c1, c2)

mul(v1, v2) = v1 * v2
@fastmath function mul(v1::N, v2::N) where N <: Normed
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    m = widemul(rv1, rv2) / convert(floattype(N), r1)
    reinterpret(N, unsafe_trunc(typeof(r1), round(m)))
end
function mul(v1::N0f8, v2::N0f8)
    m = UInt16(reinterpret(v1)) * UInt16(reinterpret(v2))
    z = m + ((m + 0x80) >> 0x8) + 0x80
    reinterpret(N0f8, unsafe_trunc(UInt8, z >> 0x8))
end
function mul(v1::N0f16, v2::N0f16)
    m = UInt32(reinterpret(v1)) * UInt32(reinterpret(v2))
    z = m + ((m + 0x8000) >> 0x10) + 0x8000
    reinterpret(N0f16, unsafe_trunc(UInt16, z >> 0x10))
end

function _blend(::BlendMode{:multiply}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = c1.l * c2.l / T(100)
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

function mul_lab(v1::T, v2::T) where T
    t, s = T(128), T(1/128)
    v1r = min(t - v1, t + v1)
    v2r = t - flipsign(v2, v1)
    flipsign(muladd(v1r, -s * v2r, t), v1)
end

_blend(::BlendMode{:screen}, c1::C, c2::C) where C <: Color =
    mapch((v1, v2) -> _n(mul(_n(v1), _n(v2))), c1, c2)

function _blend(::BlendMode{:screen}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = muladd(T(100) - c1.l, (c2.l - T(100)) / T(100), T(100))
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

_blend(::BlendMode{:overlay}, c1::C, c2::C) where C <: Color =
    _blend(BlendHardLight, c2, c1)

function _blend(::BlendMode{:overlay}, c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l1r = min(c1.l, T(100) - c1.l)
    mr = (l1r + l1r) * ifelse(c1.l == l1r, c2.l, T(100) - c2.l) / T(100)
    l = ifelse(c1.l == l1r, mr, T(100) - mr)
    C(l, mul_lab(comp2(c1), comp2(c2)), mul_lab(comp3(c1), comp3(c2)))
end

_blend(::BlendMode{:darken}, c1::C, c2::C) where C <: Color = mapch(min, c1, c2)

_blend(::BlendMode{:lighten}, c1::C, c2::C) where C <: Color = mapch(max, c1, c2)

_blend(::typeof(BlendColorDodge), c1::C, c2::C) where C <: Color = mapch(dodge, c1, c2)

function dodge(v1, v2)
    v1 == zero(v1) && return v1
    v2 == oneunit(v2) && return v2
    div01(v1, _n(v2))
end

div01(v1, v2) = convert(typeof(v1), min(1, v1 / v2))
@fastmath function div01(v1::N0f8, v2::N0f8)
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float32(rv1) * r1 / Float32(rv2)
    reinterpret(N0f8, unsafe_trunc(typeof(r1), min(Float32(r1), round(d))))
end
@fastmath function div01(v1::N0f16, v2::N0f16)
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float64(rv1) * r1 / Float64(rv2)
    reinterpret(N0f16, unsafe_trunc(typeof(r1), min(Float64(r1), round(d))))
end

_blend(::typeof(BlendColorBurn), c1::C, c2::C) where C <: Color = mapch(burn, c1, c2)

# Note that "color-burn" tends to cause the loss of significance.
function burn(v1, v2)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    _n(div01(_n(v1), v2))
end
burn(v1::Float32, v2::Float32) = Float32(burn(Float64(v1), Float64(v2)))
@fastmath function burn(v1::N0f8, v2::N0f8)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float32(r1) - (Float32(r1) - Float32(rv1)) * r1 / Float32(rv2)
    reinterpret(N0f8, unsafe_trunc(typeof(r1), max(0.0f0, round(d))))
end
@fastmath function burn(v1::N0f16, v2::N0f16)
    v1 == oneunit(v1) && return v1
    v2 == zero(v2) && return v2
    rv1, rv2, r1 = reinterpret(v1), reinterpret(v2), reinterpret(oneunit(v2))
    d = Float64(r1) - (Float64(r1) - Float64(rv1)) * r1 / Float64(rv2)
    reinterpret(N0f16, unsafe_trunc(typeof(r1), max(0.0, round(d))))
end

_blend(::typeof(BlendHardLight), c1::C, c2::C) where C <: Color = mapch(hardlight, c1, c2)

function hardlight(v1, v2)
    v2r = min(v2, _n(v2))
    mr = mul(v2r + v2r, ifelse(v2 == v2r, v1, _n(v1)))
    ifelse(v2 == v2r, mr, _n(mr))
end

_blend(::typeof(BlendHardLight), c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}} =
    _blend(BlendOverlay, c2, c1)

_blend(::typeof(BlendSoftLight), c1::C, c2::C) where C <: Color = mapch(softlight, c1, c2)

function softlight(v1, v2)
    v2r = min(v2, _n(v2))
    nv2r2 = _n(v2r + v2r)
    if v2r == v2 # if v2 <= 0.5
        return mul(v1, _n(mul(nv2r2, _n(v1))))
    else
        m = mul(v1, _n(v2))
        @fastmath d(v) = v <= 0.25 ? muladd(4v - 3, v, 1) * 4v : sqrt(v)
        return m + m + mul(nv2r2, convert(typeof(v1), d(v1)))
    end
end
softlight(v1::X, v2::X) where X <: FixedPoint = softlight(float(v1), float(v2)) % X

function _blend(::typeof(BlendSoftLight), c1::C, c2::C) where {T, C <: Union{Lab{T}, Luv{T}}}
    l = softlight(c1.l / T(100), c2.l / T(100))
    a = softlight(muladd(comp2(c1), T(1/256), T(0.5)),
                  muladd(comp2(c2), T(1/256), T(0.5)))
    b = softlight(muladd(comp3(c1), T(1/256), T(0.5)),
                  muladd(comp3(c2), T(1/256), T(0.5)))
    C(l * T(100), muladd(a, T(256), T(-128)) , muladd(b, T(256), T(-128)))
end

_blend(::BlendMode{:difference}, c1::C, c2::C) where C <: Color = mapch(difference, c1, c2)

difference(v1, v2) = max(v1, v2) - min(v1, v2)

_blend(::BlendMode{:exclusion}, c1::C, c2::C) where C <: Color = mapch(exclusion, c1, c2)

function exclusion(v1, v2)
    m = mul(v1, v2)
    (v1 - m) + (v2 - m)
end
function exclusion(v1::N0f8, v2::N0f8)
    rv1, rv2 = UInt16(reinterpret(v1)), UInt16(reinterpret(v2))
    m = rv1 * rv2
    z = muladd(rv1 + rv2, 0xFF, -(m + m))
    reinterpret(N0f8, unsafe_trunc(UInt8, (z + ((z + 0x80) >> 0x8) + 0x80) >> 0x8))
end
function exclusion(v1::N0f16, v2::N0f16)
    rv1, rv2 = UInt32(reinterpret(v1)), UInt32(reinterpret(v2))
    m = rv1 * rv2
    z = muladd(rv1 + rv2, 0xFFFF, -(m + m))
    reinterpret(N0f16, unsafe_trunc(UInt16, (z + ((z + 0x8000) >> 0x10) + 0x8000) >> 0x10))
end

lum100(c::AbstractRGB) = lum100(red(c), green(c), blue(c))
lum100(r::T, g::T, b::T) where T = muladd(T(59), g, muladd(T(11), b, T(30) * r))
function lum100(r::T, g::T, b::T) where T <: Normed
    F = floattype(T)
    rr = reinterpret(r)
    rg = reinterpret(g)
    rb = reinterpret(b)
    ri = reinterpret(oneunit(T))
    muladd(F(59), rg, muladd(F(11), rb, F(30) * rr)) / F(ri)
end
function lum100(r::T, g::T, b::T) where T <: Normed{UInt8}
    F = floattype(T)
    rr = reinterpret(r)
    rg = reinterpret(g)
    rb = reinterpret(b)
    unsafe_trunc(Int32, 0x268cf3 * rr + 0x4bd0f0 * rg + 0xe229d * rb) * 4.65661322835f-8
end
function setlum(c::C, l::AbstractFloat) where {T, C <: AbstractRGB{T}}
    lc = lum100(c)
    r, g, b = red(c), green(c), blue(c)
    n = min(r, g, b)
    x = max(r, g, b)

    d = (lc - l) / 100
    if n < d
        # k = l / (lc - 100n)
        # cr, cg, cb = map(c -> c * k - n * k, (r, g, b))
        if n == r
            cr = zero(T)
            cb = (b - n) * l / muladd(T(-70), r - g, T(11) * (b - g))
            cg = max(cr, muladd(T(-11), cb, l) * T(1/59))
        elseif n == g
            cg = zero(T)
            cb = (b - n) * l / muladd(T(30), r - g, T(11) * (b - g))
            cr = max(cg, muladd(T(-11), cb, l) * T(1/30))
        else
            cb = zero(T)
            cr = (r - n) * l / muladd(T(30), r - g, T(-89) * (b - g))
            cg = max(cb, muladd(T(-30), cr, l) * T(1/59))
        end
        return C(cr, cg, cb)
    elseif d < x - oneunit(x)
        # k1 = 1 / (100x - lc)
        # kl = l / (100x - lc)
        # mapc(v -> x * kl - v * kl + (100v - lc) * k1, c)
        if x == r
            k1 = muladd(T(70), r - g, T(-11) * (b - g))
            k2 = muladd(T(-30), r - g, T(89) * (b - g))
            cr = oneunit(T)
            cb = muladd(x - b, l, k2) / k1
            cg = min(cr, muladd(T(-11), cb, l - T(30)) * T(1/59))
        elseif x == g
            k1 = muladd(T(-30), r - g, T(-11) * (b - g))
            k2 = muladd(T(-30), r - g, T(89) * (b - g))
            cg = oneunit(T)
            cb = muladd(x - b, l, k2) / k1
            cr = min(cg, muladd(T(-11), cb, l - T(59)) * T(1/30))
        else
            k1 = muladd(T(-30), r - g, T(89) * (b - g))
            k2 = muladd(T(70), r - g, T(-11) * (b - g))
            cb = oneunit(T)
            cr = muladd(x - r, l, k2) / k1
            cg = min(cb, muladd(T(-30), cr, l - T(11)) * T(1/59))
        end
        return C(cr, cg, cb)
    end
    return mapc(v -> v - d, c)
end
function setlum(c::C, l::AbstractFloat) where {T <: Normed, C <: AbstractRGB{T}}
    convert(C, setlum(convert(RGB{floattype(T)}, c), l))
end

sat(c::AbstractRGB) = sat(red(c), green(c), blue(c))
sat(r::T, g::T, b::T) where T = max(r, g, b) - min(r, g, b)
function setsat(c::C, s) where C <: AbstractRGB
    r, g, b = float(red(c)), float(green(c)), float(blue(c))
    cmin_rg, cmax_rg = minmax(r, g)
    cmin = min(cmin_rg, b)
    cmax = max(cmax_rg, b)
    cmax == cmin && return C(zero(eltype(C)))

    cmid = max(cmin_rg, min(cmax_rg, b))

    cmid2 = ((cmid - cmin) * float(s)) / (cmax - cmin)
    cmax2 = float(s)

    r2 = ifelse(cmax == r, cmax2, ifelse(cmin == r, zero(r), cmid2))
    g2 = ifelse(cmax == g, cmax2, ifelse(cmin == g, zero(g), cmid2))
    b2 = ifelse(cmax == b, cmax2, ifelse(cmin == b, zero(b), cmid2))

    C(r2, g2, b2)
end

_blend(m::M, c1::C, c2::C) where {C <: Color, M <: NonSeparableBlendMode} =
    convert(C, _blend(m, convert(RGB, c1), convert(RGB, c2)))

_blend(::BlendMode{:hue}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(setsat(c2, sat(c1)), lum100(c1))

function _blend(::BlendMode{:hue}, c1::C, c2::C) where C <: Union{Lab, Luv}
    if abs(comp2(c2)) >= abs(comp3(c2))
        a = copysign(max(abs(comp2(c1)), abs(comp3(c1))), comp2(c2))
        b = ifelse(a == zero(a), a, a * comp3(c2) / comp2(c2))
    else
        b = copysign(max(abs(comp2(c1)), abs(comp3(c1))), comp3(c2))
        a = ifelse(b == zero(b), b, b * comp2(c2) / comp3(c2))
    end
    C(c1.l, a, b)
end

_blend(::BlendMode{:hue}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c2.h, c1.s, comp3(c1))

_blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(setsat(c1, sat(c2)), lum100(c1))

function _blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: Union{Lab, Luv}
    if abs(comp2(c1)) >= abs(comp3(c1))
        a = copysign(max(abs(comp2(c2)), abs(comp3(c2))), comp2(c1))
        b = ifelse(a == zero(a), a, a * comp3(c1) / comp2(c1))
    else
        b = copysign(max(abs(comp2(c2)), abs(comp3(c2))), comp3(c1))
        a = ifelse(b == zero(b), b, b * comp2(c1) / comp3(c1))
    end
    C(c1.l, a, b)
end

_blend(::BlendMode{:saturation}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c1.h, c2.s, comp3(c1))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(c2, lum100(c1))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: Union{Lab, Luv} =
    C(c1.l, comp2(c2), comp3(c2))

_blend(::BlendMode{:color}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c2.h, c2.s, comp3(c1))

_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: AbstractRGB =
    setlum(c1, lum100(c2))

_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: Union{Lab, Luv} =
    C(c2.l, comp2(c1), comp3(c1))

_blend(::BlendMode{:luminosity}, c1::C, c2::C) where C <: Union{HSV, HSL, HSI} =
    C(c1.h, c1.s, comp3(c2))


text:
using EBayes
using Documenter

makedocs(
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
        "Reference" => ["Empirical Bayes Samples" => "ebayes_samples.md",
                        "Example Priors" => "prior_distributions.md"],
    ],
    repo="https://github.com/nignatiadis/EBayes.jl/blob/{commit}{path}#L{line}",
    sitename="EBayes.jl",
    authors="Nikos Ignatiadis <nikos.ignatiadis01@gmail.com>"
)

deploydocs(;
    repo="github.com/nignatiadis/EBayes.jl",
)


text:
using CSV
using Distributions


text:
using AmbientForcing
using NetworkDynamics, Distributions
using LightGraphs
using OrdinaryDiffEq

const N_plastic = 10 # number of nodes
k = 4  # average degree
g = barabasi_albert(N_plastic, k)

function kuramoto_plastic_edge!(de, e, v_s, v_d, p, t)
    de[1] =  e[2] * sin(v_s[1] - v_d[1] + α) / N_plastic - e[1]
    de[2] = - ϵ * (sin(v_s[1] - v_d[1] + β) + e[2])

    nothing
end

function kuramoto_plastic_vertex!(dv, v, edges, p, t)
    dv .= 0
    for e in edges
        dv .-= e[1]
    end
end

const ϵ = 0.1
const α = .2π
const β = -.95π

plasticvertex = ODEVertex(f! = kuramoto_plastic_vertex!, dim =1)
mass_matrix_plasticedge = zeros(2,2)
mass_matrix_plasticedge[2,2] = 1. # First variables is set to 0

plasticedge = ODEEdge(f! = kuramoto_plastic_edge!, dim=2, sym=[:e, :de], coupling=:undirected,mass_matrix = mass_matrix_plasticedge);
kuramoto_plastic! = network_dynamics(plasticvertex, plasticedge, g)

x0_plastic = rand(106)
g_nd = constraint_equations(kuramoto_plastic!)
sum(g_nd(x0_plastic))

x0_plastic = zeros(106)
g_nd = constraint_equations(kuramoto_plastic!)
g_nd(x0_plastic)

Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform)
z_new = ambient_forcing(kuramoto_plastic!, x0_plastic, 2.0, Frand)

sum(g_nd(z_new))

idx = idx_exclusive(kuramoto_plastic!, ["e_22", "de_22"])
Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform, idx)
z_new = ambient_forcing(kuramoto_plastic!, x0_plastic, 2.0, Frand)

sum(g_nd(z_new))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# lid.jl
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

println("[START] starting training LID Model")
@time bkgmodel, fextractor, model = lid_train("data/text.tsv");
println("[END]   finished training LID Model")

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

println("[START] starting Test Set 1")

server = text_rest_server(fextractor, model)

@async run(server, host=ip"127.0.0.1", port=8000)

ret = Requests.get("http://localhost:8000/lid")
res = Requests.json(ret)
println(res)
@test res["resource"] == "lid"
@test length(res["all_languages"]) == 18

ret = Requests.put("http://localhost:8000/lid"; data = "hello world")
res = Requests.json(ret)
println(res)
@test res["resource"] == "lid"
@test res["text"] == "hello world"
@test res["language"] == "en"

ret = Requests.get("http://localhost:8000/blah")
res = Requests.json(ret)
println(res)
@test res["resource"] == 404

# close(server)
close(server.http.sock)

println("[END]   finished Test Set 1")

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

println("[START] starting saving LID Model")
@time lid_save("data/test.jld", bkgmodel, model)
println("[END]   finished saving LID Model")

println("[START] starting loading LID Model")
@time bkgmodel, fextractor, model = lid_load("data/test.jld")
println("[END]   finished loading LID Model")

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

println("[START] starting Test Set 2")

server = text_rest_server(fextractor, model)

@async run(server, host=ip"127.0.0.1", port=8000)

ret = Requests.get("http://localhost:8000/lid")
res = Requests.json(ret)
println(res)
@test res["resource"] == "lid"
@test length(res["all_languages"]) == 18

ret = Requests.put("http://localhost:8000/lid"; data = "hello world")
res = Requests.json(ret)
println(res)
@test res["resource"] == "lid"
@test res["text"] == "hello world"
@test res["language"] == "en"

ret = Requests.get("http://localhost:8000/blah")
res = Requests.json(ret)
println(res)
@test res["resource"] == 404

# close(server)
close(server.http.sock)

println("[END]   finished Test Set 2")

max_forks_repo_name: trevorlewis/TEXT-Language-REST


text:
import MLBase
import NumericalIntegration
import StatsBase

"""
"""
function auprc(
        ytrue::AbstractVector{<:Integer},
        yscore::AbstractVector{<:Real},
        )
    allprecisions, allrecalls, allthresholds = prcurve(
        ytrue,
        yscore,
        )
    #
    permutation = sortperm(allthresholds; rev=true)
    allprecisions = allprecisions[permutation]
    allrecalls = allrecalls[permutation]
    allthresholds = allthresholds[permutation]
    #
    x = allrecalls
    y = allprecisions
    areaunderprcurve = trapz(x, y)
    return areaunderprcurve
end