text:
function gsreg(
    equation::String;
    data=nothing,
    datanames=nothing,
    intercept=INTERCEPT_DEFAULT,
    outsample=OUTSAMPLE_DEFAULT,
    criteria=CRITERIA_DEFAULT,
    ttest=TTEST_DEFAULT,
    method=METHOD_DEFAULT,
    vectoroperation=VECTOR_OPERATION_DEFAULT,
    modelavg=MODEL_AVG_DEFAULT,
    residualtest=RESIDUAL_TEST_DEFAULT,
    time=TIME_DEFAULT,
    summary=SUMMARY_DEFAULT,
    csv=CSV_DEFAULT,
    resultscsv=CSV_DEFAULT,
    orderresults=ORDER_RESULTS_DEFAULT,
    onmessage=ON_MESSAGE_DEFAULT,
    parallel=PARALLEL_DEFAULT
    )

    return gsreg(
        equation,
        data,
        datanames=datanames,
        intercept=intercept,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        method=method,
        vectoroperation=vectoroperation,
        modelavg=modelavg,
        residualtest=residualtest,
        time=time,
        summary=summary,
        resultscsv=resultscsv,
        csv=csv,
        orderresults=orderresults,
        onmessage=onmessage,
        parallel=parallel
    )
end

function gsreg(
    equation::String,
    data;
    datanames=nothing,
    intercept=INTERCEPT_DEFAULT,
    outsample=OUTSAMPLE_DEFAULT,
    criteria=CRITERIA_DEFAULT,
    ttest=TTEST_DEFAULT,
    method=METHOD_DEFAULT,
    vectoroperation=VECTOR_OPERATION_DEFAULT,
    modelavg=MODEL_AVG_DEFAULT,
    residualtest=RESIDUAL_TEST_DEFAULT,
    time=TIME_DEFAULT,
    summary=SUMMARY_DEFAULT,
    csv=CSV_DEFAULT,
    resultscsv=CSV_DEFAULT,
    orderresults=ORDER_RESULTS_DEFAULT,
    onmessage=ON_MESSAGE_DEFAULT,
    parallel=PARALLEL_DEFAULT
    )

    equation = equation_str_to_strarr(equation)

    return gsreg(
        equation,
        data,
        datanames=datanames,
        intercept=intercept,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        method=method,
        vectoroperation=vectoroperation,
        modelavg=modelavg,
        residualtest=residualtest,
        time=time,
        summary=summary,
        resultscsv=resultscsv,
        csv=csv,
        orderresults=orderresults,
        onmessage=onmessage,
        parallel=parallel
    )
end

function gsreg(
    equation::Array{String},
    data;
    datanames=nothing,
    intercept=INTERCEPT_DEFAULT,
    outsample=OUTSAMPLE_DEFAULT,
    criteria=CRITERIA_DEFAULT,
    ttest=TTEST_DEFAULT,
    method=METHOD_DEFAULT,
    vectoroperation=VECTOR_OPERATION_DEFAULT,
    modelavg=MODEL_AVG_DEFAULT,
    residualtest=RESIDUAL_TEST_DEFAULT,
    time=TIME_DEFAULT,
    summary=SUMMARY_DEFAULT,
    csv=CSV_DEFAULT,
    resultscsv=CSV_DEFAULT,
    orderresults=ORDER_RESULTS_DEFAULT,
    onmessage=ON_MESSAGE_DEFAULT,
    parallel=PARALLEL_DEFAULT
    )

    datanames = get_datanames(data, datanames)
    equation = equation_strarr_to_symarr(equation, datanames)

    if isempty(equation)
        error(VARIABLES_NOT_OR_VALID_OR_NOT_DEFINED)
    end

    return gsreg(
        equation,
        data,
        datanames=datanames,
        intercept=intercept,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        method=method,
        vectoroperation=vectoroperation,
        modelavg=modelavg,
        residualtest=residualtest,
        time=time,
        summary=summary,
        resultscsv=resultscsv,
        csv=csv,
        orderresults=orderresults,
        onmessage=onmessage,
        parallel=parallel
    )
end

function gsreg(
    equation::Array{Symbol},
    data;
    datanames=nothing,
    intercept=INTERCEPT_DEFAULT,
    outsample=OUTSAMPLE_DEFAULT,
    criteria=CRITERIA_DEFAULT,
    ttest=TTEST_DEFAULT,
    method=METHOD_DEFAULT,
    vectoroperation=VECTOR_OPERATION_DEFAULT,
    modelavg=MODEL_AVG_DEFAULT,
    residualtest=RESIDUAL_TEST_DEFAULT,
    time=TIME_DEFAULT,
    summary=SUMMARY_DEFAULT,
    csv=CSV_DEFAULT,
    resultscsv=CSV_DEFAULT,
    orderresults=ORDER_RESULTS_DEFAULT,
    onmessage=ON_MESSAGE_DEFAULT,
    parallel=PARALLEL_DEFAULT
    )

    datanames = datanames_strarr_to_symarr!(datanames)
    depvar = equation[1]
    expvars = equation[2:end]
    data = convert_if_is_tuple_to_array(data, datanames)

    if time != nothing
        data = sort_data_by_time(data, time, datanames)
    end

    data = filter_data_valid_columns(data, depvar, expvars, datanames)
    data = filter_rows_with_empty_values(data)
    data = convert_if_is_dataframe_to_array(data)

    return gsreg(
        equation,
        data,
        datanames;
        intercept=intercept,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        method=method,
        vectoroperation=vectoroperation,
        modelavg=modelavg,
        residualtest=residualtest,
        time=time,
        summary=summary,
        resultscsv=resultscsv,
        csv=csv,
        orderresults=orderresults,
        onmessage=onmessage,
        parallel=parallel
    )
end

function gsreg(
    equation::Array{Symbol},
    data,
    datanames::Array;
    intercept=INTERCEPT_DEFAULT,
    outsample=OUTSAMPLE_DEFAULT,
    criteria=CRITERIA_DEFAULT,
    ttest=TTEST_DEFAULT,
    method=METHOD_DEFAULT,
    vectoroperation=VECTOR_OPERATION_DEFAULT,
    modelavg=MODEL_AVG_DEFAULT,
    residualtest=RESIDUAL_TEST_DEFAULT,
    time=TIME_DEFAULT,
    summary=SUMMARY_DEFAULT,
    csv=CSV_DEFAULT,
    resultscsv=CSV_DEFAULT,
    orderresults=ORDER_RESULTS_DEFAULT,
    onmessage=ON_MESSAGE_DEFAULT,
    parallel=PARALLEL_DEFAULT
    )

    if method == "precise"
        datatype = Float64
    elseif method == "fast"
        datatype = Float32
    else
        error(METHOD_INVALID)
    end

    if !isa(data, Array{datatype})
        data = convert(Array{datatype}, data)
    end

    if outsample != OUTSAMPLE_DEFAULT
        if outsample < 0
            error(OUTSAMPLE_LOWER_VALUE)
        elseif size(data, 1) - outsample < INSAMPLE_MIN_SIZE + size(data, 2) - 1
            error(OUTSAMPLE_HIGHER_VALUE)
        end
    end

    if outsample == false && :rmseout in criteria
        error(OUTSAMPLE_MISMATCH)
    end

    if criteria == CRITERIA_DEFAULT
        if outsample != OUTSAMPLE_DEFAULT
            criteria = CRITERIA_DEFAULT_OUTSAMPLE
        else
            criteria = CRITERIA_DEFAULT_INSAMPLE
        end
    end

    if time != nothing && time ∉ datanames
        error(TIME_VARIABLE_INEXISTENT)
    end

    if resultscsv != csv
        if resultscsv != CSV_DEFAULT && csv != CSV_DEFAULT
            error(CSV_DUPLICATED_PARAMETERS)
        elseif csv != CSV_DEFAULT
            resultscsv = csv
        end
    end

    # TODO: Is this been used?
    if parallel != nothing
        if parallel > nworkers()
            error("Number of parallel workers can not exceed available cores. Use addprocs()")
        end

        if parallel < 1
            error("Number of workers can not be less than 1")
        end
    end

    if size(data, 1) < size(equation[2:end], 1) + 1
        error(NO_ENOUGH_OBSERVATIONS)
    end

    if !in_vector(equation, datanames)
        error(SELECTED_VARIABLES_DOES_NOT_EXISTS)
    end

    depvar = equation[1]
    expvars = equation[2:end]

    expvars_data = data[ 1:end, 2:end ]
    corrmatrix = cor(expvars_data)
    s = size(corrmatrix, 1)
    corrminusiden = corrmatrix - Matrix{Float64}(I, s,s)
    maxcorr = maximum(corrminusiden)

    if maxcorr > 0.999
        error("Your explanatory variables are not linearly independent -perfect multicollinearity among covariates-")
    end

    result = gsreg(
        depvar,
        expvars,
        data,
        datanames=datanames,
        intercept=intercept,
        outsample=outsample,
        criteria=criteria,
        ttest=ttest,
        vectoroperation=vectoroperation,
        modelavg=modelavg,
        residualtest=residualtest,
        time=time,
        summary=summary,
        datatype=datatype,
        orderresults=orderresults,
        onmessage=onmessage,
        parallel=parallel
    )
    if resultscsv != nothing
        export_csv(resultscsv, result)
    end
    return result
end


text:
"""
    find_ratbasis(list::AbstractVector{<:StaticVector{N,<:Rational{T}}}) where {N,T}

Given a list of `N`-dimensional rational vectors, return a `N×N` matrix whose columns
are independent vectors from the list.
If no such `N` vectors exist, the returned matrix will not be invertible.
"""
function find_ratbasis(list::AbstractVector{<:StaticVector{N,<:Rational{T}}}) where {N,T}
    U = widen(T)
    n = length(list)
    ratbasis = zero(SizedMatrix{N,N,Rational{U}})
    i = 1
    @inbounds while i <= n
        l = list[i]
        if !iszero(l)
            ratbasis[:,1] .= l
            i += 1
            break
        end
        i += 1
    end
    j = 2
    while i <= n && j <= N
        l = list[i]
        if iszero(l)
            i += 1
            continue
        end
        ratbasis[:,j] .= l
        if rank(@view ratbasis[:,1:j]) == j
            j += 1
        end
        i += 1
    end
    return ratbasis
end

# Specialized version for N = 3
function find_ratbasis(list::AbstractVector{<:StaticVector{3,<:Rational{T}}}) where T
    U = widen(T)
    n = length(list)
    ratbasis = zero(SizedMatrix{3,3,Rational{U}})
    # We start by computing a basis for the i first vectors of v where i (≤ n) is minimal
    i = 1
    @inbounds while i <= n
        if !iszero(list[i])
            ratbasis[:,3] .= list[i]
            i += 1
            break
        end
        i += 1
    end
    @inbounds while i <= n
        v = list[i]
        valid = false
        # v is valid if it is linearly independant of ratbasis[:,3]
        if iszero(v[1])
            valid = (!iszero(ratbasis[1,3]) & (!iszero(v[2]) | !iszero(v[3]))) ||
                    (v[2] * ratbasis[3,3] != v[3] * ratbasis[2,3])
        else
            x = ratbasis[1,3] // v[1]
            valid = ratbasis[2,3] != x * v[2] || ratbasis[3,3] != x * v[3]
        end
        if valid
            ratbasis[:,2] .= v
            i += 1
            break
        end
        i += 1
    end
    d = zero(Rational{U})
    @inbounds while i <= n
        v = list[i]
        ratbasis[:,1] .= v
        d = det(ratbasis)
        iszero(d) || break
        i += 1
    end
    # At this point ratbasis is a basis in which the i-1 first coordinates are
    # [0, 0, 1] followed by a (possibly empty) series of [0, 0, x], then
    # [0, 1, 0] followed by a (possibly empty) series of [0, x, y], then [1, 0, 0]
    @inbounds if d < 0
        ratbasis[:,1] .= .-ratbasis[:,1] # Ensure a direct transformation
    end
    return SMatrix{3,3,Rational{U}}(ratbasis)
end

# Specialized version for N = 2
function find_ratbasis(list::AbstractVector{<:StaticVector{2,<:Rational{T}}}) where T
    U = widen(T)
    n = length(list)
    ratbasis = zero(SizedMatrix{2,2,Rational{U}})
    i = 1
    @inbounds while i <= n
        if !iszero(list[i])
            ratbasis[:,1] .= list[i]
            i += 1
            break
        end
        i += 1
    end
    d = zero(Rational{U})
    @inbounds while i <= n
        v = list[i]
        ratbasis[:,2] .= v
        d = det(ratbasis)
        iszero(d) || break
        i += 1
    end
    @inbounds if d < 0
        ratbasis[:,1] .= .-ratbasis[:,1]
    end
    return SMatrix{2,2,Rational{U}}(ratbasis)
end

"""
    normal_basis_rational(list::AbstractVector{<:StaticVector{N,<:Rational{T}}}) where {N,T}

Given a list of `N`-dimensional rational vectors, return a basis for the space
spanned by integer combinations of these vectors.

This basis is deterministically computed from the input.

It should depend only on the spanned space, not on the exact input (although this
assertion should be considered experimental for now).
"""
function normal_basis_rational(list::AbstractVector{<:StaticVector{N,<:Rational{T}}}) where {N,T}
    U = widen(T)
    n = length(list)
    ratbasis = find_ratbasis(list)
    # ratbasis should always be invertible in our setting: if not, it means that
    # we missed the fact that the dimensionality of the graph was strictly lower
    # than N, which should have been detected much earlier.
    # In this case, the definition of invratbasis will cause a failure.

    lcms = SizedVector{N,T}(ones(T, N))
    # per-column lcm of the vectors expressed in ratbasis
    expressed_in_ratbasis = Vector{SVector{N,Rational{T}}}(undef, n)
    invratbasis = inv(ratbasis) # See comment above in case of error here.
    @inbounds for j in 1:n
        expressed_in_ratbasis[j] = invratbasis * list[j]
        lcms .= lcm.(lcms, denominator.(expressed_in_ratbasis[j]))
    end

    expressed_in_basis = [MVector{N,Int}(undef) for _ in 1:n]
    @inbounds for j in 1:n
        x = expressed_in_ratbasis[j]
        expressed_in_basis[j] .= div.(numerator.(x).*lcms, denominator.(x))
    end
    #= lcms contains the lcm of all denominators of the coordinates of translations
       expressed in ratbasis.
       We now work in an artificial basis in which the coordinates are all integers.
       The integer coordinates are stored in expressed_in_basis.
       expressed_in_basis is now going to be reduced to its normal form.
    =#

    intbasis, _ = PeriodicGraphs.normal_basis(expressed_in_basis)

    # We finally go back from the integer normal form to a rational basis that expresses
    # all the original vectors with integer coordinates
    return SMatrix{N,N,Rational{U}}(ratbasis * [@inbounds(intbasis[i,j] // lcms[i]) for i in 1:N, j in 1:N])

    # newcoords = [SVector{N,Int}(normalization * x) for x in intcoords]
    # @toggleassert (newbasis,) .* newcoords == list

    # return newbasis, newcoords
end

# Specialized version for N = 1
function normal_basis_rational(list::AbstractVector{<:StaticVector{1,<:Rational{T}}}) where T
    nzl = [x[] for x in list if !iszero(x[])]
    den = lcm(denominator.(nzl))
    S = widen(T)
    num = minimum(S[abs(S(den*x)) for x in nzl])
    return SMatrix{1,1,Rational{S}}((num // den,))
end

ext: jl


text:
# draw histgram of sensor data

module HistogramMean
    using Plots, DataFrames, CSV, Statistics
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df_200_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        a = sum(df_200_mm.lidar)
        println("Sum of sensor data = $a")

        b = length(df_200_mm.lidar)
        println("Length of sensor data = $b")

        mean_1 = a / b
        mean_2 = mean(df_200_mm.lidar)
        println("Mean(Sum / Length) = $mean_1")
        println("Mean(mean()) = $mean_2")

        bin_min_max = maximum(df_200_mm.lidar) - minimum(df_200_mm.lidar)
        histogram(df_200_mm.lidar, bins=bin_min_max, color=:orange,
                  label="histogram")
        plot!([mean_1], st=:vline, color=:red,
              ylim=(0, 5000), linewidth=5, label="mean")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/freq_dist/histogram/histogram_mean.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using CUDA

@inline function CheckBlockBounds(offset, block, grids)
        x::Int32 = offset.x + block.x
        if x < 1 || x > grids
                return true
        end
        x = offset.y + block.y
        if x < 1 || x > grids
                return true
        end
        x = offset.z + block.z
        if x < 1 || x > grids
                return true
        end
        return false
end
@inline function BoundsDistance(offset, query, dist)

    dx::Float32 = (offset.x == 1)  ? (1/grids - mod(query[1], 1/grids)) : 0
    dx = (offset.x == -1) ? (mod(query[1], 1/grids)) : dx
    dx = CUDA.pow(dx,Int32(2))
    dy::Float32 = (offset.y == 1)  ? (1/grids - mod(query[2], 1/grids)) : 0
    dy = (offset.y == -1) ? (mod(query[2], 1/grids)) : dy
    dy = CUDA.pow(dy,Int32(2))
    dz::Float32 = (offset.z == 1)  ? (1/grids - mod(query[3], 1/grids)) : 0
    dz = (offset.z == -1) ? (mod(query[3], 1/grids)) : dz
    dz = CUDA.pow(dz,Int32(2))
    boundsDist::Float32 = CUDA.sqrt(dx+dy+dz)
    local_skip = boundsDist > dist

    return local_skip

end
function cuGridKnnWithSkip(Points, Queries,
        Pointsperblock, Queriesperblock,
        IntPointsperblock, IntQueriesperblock,
        Distances, Neighbours,
        dimensions, grids, offset)

    if CheckBlockBounds(offset, blockIdx(), grids)
        return nothing
    end
    q_bid::Int32 = ((blockIdx().z - 1)*gridDim().y*gridDim().x
            + (blockIdx().y - 1)*gridDim().x
            + blockIdx().x)
    p_bid::Int32 = ((blockIdx().z + offset.z- 1)*gridDim().y*gridDim().x
            + (blockIdx().y + offset.y - 1)*gridDim().x
            + blockIdx().x + offset.x)
    tid::Int32 = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride::Int32 = (blockDim().x)*(blockDim().y)

    SharedPoints = @cuDynamicSharedMem(Float32, (dimensions, stride))
    SharedQueries = @cuDynamicSharedMem(Float32, (dimensions, stride),
                                dimensions*stride*sizeof(Float32))
    SharedSkip = @cuDynamicSharedMem(Bool, (blockDim().y),
                                2*dimensions*stride*sizeof(Float32))
    query = @view SharedQueries[:, tid]

    startPoints::Int32 = 0
    (threadIdx().x == 1)  && (@inbounds  startPoints = IntPointsperblock[p_bid])
    startPoints = shfl_sync(0xffffffff, startPoints, 1);

    startQueries::Int32 = 0
    (threadIdx().x == 1) && (@inbounds startQueries = IntQueriesperblock[q_bid])
    startQueries = shfl_sync(0xffffffff, startQueries, 1);

    totalPoints::Int32 = 0
    (threadIdx().x == 1) && (@inbounds totalPoints = Pointsperblock[p_bid])
    totalPoints = shfl_sync(0xffffffff, totalPoints, 1);

    totalQueries::Int32 = 0
    (threadIdx().x == 1)  && (@inbounds totalQueries = Queriesperblock[q_bid])
    totalQueries = shfl_sync(0xffffffff, totalQueries, 1);

    dist::Float32 = 0
    nb::Int32 = 0
    local_skip::Bool = true
    global_skip::Bool = true
    for q::Int32 = 0:stride:(totalQueries - 1)
        ## Read queries
        qIndex::Int32 = startQueries + q + tid
        if tid + q <= totalQueries
                for d::Int32 = 1:dimensions
                        @inbounds query[d] = Queries[qIndex, d]
                end
                @inbounds dist = Distances[qIndex]
                @inbounds nb = Neighbours[qIndex]

                dx::Float32 = (offset.x == 1)  ? (1/grids - mod(query[1], 1/grids)) : 0
                dx = (offset.x == -1) ? (mod(query[1], 1/grids)) : dx
                dx = CUDA.pow(dx,Int32(2))
                dy::Float32 = (offset.y == 1)  ? (1/grids - mod(query[2], 1/grids)) : 0
                dy = (offset.y == -1) ? (mod(query[2], 1/grids)) : dy
                dy = CUDA.pow(dy,Int32(2))
                dz::Float32 = (offset.z == 1)  ? (1/grids - mod(query[3], 1/grids)) : 0
                dz = (offset.z == -1) ? (mod(query[3], 1/grids)) : dz
                dz = CUDA.pow(dz,Int32(2))
                boundsDist::Float32 = CUDA.sqrt(dx+dy+dz)
                local_skip = boundsDist > dist
        end

        local_skip = vote_all(local_skip)
        (threadIdx().x==1) && (SharedSkip[threadIdx().y] = local_skip)
        sync_threads()
        if threadIdx().x==1
            for w = 1:blockDim().y
                global_skip = global_skip && SharedSkip[w]
            end
        end
        global_skip = shfl_sync(FULL_MASK, global_skip, 1)
        global_skip && continue
        for p::Int32 = 0:stride:(totalPoints - 1)
            ## Read points
            sync_threads()
            if tid + p <= totalPoints
                for d::Int32 = 1:dimensions
                    @inbounds SharedPoints[d, tid] = Points[startPoints + p + tid, d]
                end
            end
            sync_threads()
            ## Calculate and compare distance
            local_skip && continue
            bounds::Int32 = CUDA.min(stride, totalPoints-p)
            for i::Int32 = 1:bounds
                tempdist::Float32 = 0
                for d::Int32 = 1:dimensions
                    @inbounds tempdist += CUDA.pow(query[d]
                                        -SharedPoints[d, (i+tid-2)%bounds+1],2)
                end
                tempdist = CUDA.sqrt(tempdist)
                if tempdist < dist
                    dist = tempdist
                    nb = startPoints + p + (i+tid-2)%bounds+1
                end
            end
        end
        if tid + q <= totalQueries
            @inbounds Distances[qIndex] = dist
            @inbounds Neighbours[qIndex] = nb
        end
    end
    return nothing
end

function cuda_knn_with_skip(OrderedPoints, OrderedQueries,
    PointsPerBlock, QueriesPerBlock,
    IntegralPointsPerBlock, IntegralQueriesPerBlock,
    numOfPoints, numOfQueries, numOfGrids, dimensions)

    ## Data transfer
    devPoints = CuArray(OrderedPoints)
    devQueries = CuArray(OrderedQueries)
    devPointsPerBlock = CuArray(PointsPerBlock)
    devQueriesPerBlock = CuArray(QueriesPerBlock)
    devIntegralPointsPerBlock = CuArray(IntegralPointsPerBlock)
    devIntegralQueriesPerBlock = CuArray(IntegralQueriesPerBlock)
    devRes = CUDA.fill(Float32(100), numOfQueries)
    devNeighbours = CuArray(zeros(Int32, numOfQueries))
    ## Config
    dimensions = Int32(dimensions)
    grids = Int32(numOfGrids)
    offset = (x=Int32(0),y=Int32(0),z=Int32(0))
    thread_groups = 2
    shmem=2*thread_groups*32*dimensions*sizeof(Float32) + thread_groups*sizeof(Int32)
    ## Running for main block.
    @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
      cuGridKnnWithSkip(devPoints, devQueries,
        devPointsPerBlock, devQueriesPerBlock,
        devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
        devRes, devNeighbours, dimensions, grids, offset))
    #Running for neighbour blocks.
    for x::Int32 = -1:1
        for y::Int32 = -1:1
            for z::Int32 = -1:1
                if x == 0 && y == 0 && z == 0
                    continue
                end
                offset = (x=x,y=y,z=z)
                @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
                  cuGridKnnWithSkip(devPoints, devQueries,
                    devPointsPerBlock, devQueriesPerBlock,
                    devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
                    devRes, devNeighbours, dimensions, grids, offset))
            end
        end
    end

    return Array(devNeighbours), Array(devRes)

end


text:
module Macros
import Base.@time

macro __LOCATION__(); return QuoteNode(__source__); end

dump(
    @__LOCATION__()
)

time() = println("macros")

@time time()

end

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
@testset "maha" begin
    rng = MersenneTwister(123456)
    D_in = 3
    v1 = rand(rng, D_in)
    v2 = rand(rng, D_in)

    U = UpperTriangular(rand(rng, 3, 3))
    P = Matrix(Cholesky(U, 'U', 0))
    @assert isposdef(P)

    k = @test_deprecated MahalanobisKernel(; P=P)
    @test k isa TransformedKernel{SqExponentialKernel,<:LinearTransform}
    @test k.transform.A ≈ sqrt(2) .* U
    @test k(v1, v2) ≈ exp(-sqmahalanobis(v1, v2, P))

    # Standardised tests.
    @testset "ColVecs" begin
        x0 = ColVecs(randn(D_in, 3))
        x1 = ColVecs(randn(D_in, 3))
        x2 = ColVecs(randn(D_in, 2))
        TestUtils.test_interface(k, x0, x1, x2)
    end
    @testset "RowVecs" begin
        x0 = RowVecs(randn(3, D_in))
        x1 = RowVecs(randn(3, D_in))
        x2 = RowVecs(randn(2, D_in))
        TestUtils.test_interface(k, x0, x1, x2)
    end
end


text:
"Returns the time-scale for advection on a regular grid across a single grid cell
 for ShallowWaterModel."

import Oceananigans.Utils: cell_advection_timescale

function shallow_water_cell_advection_timescale(uh, vh, h, grid, formulation)

    Δxmin = minimum(grid.Δxᶜᵃᵃ)
    Δymin = minimum(grid.Δyᵃᶜᵃ)

    uhmax = maximum(abs, uh)
    vhmax = maximum(abs, vh)
    hmin  = minimum(abs,  h)

    return min(Δxmin / uhmax, Δymin / vhmax) * hmin
end

function shallow_water_cell_advection_timescale(u, v, h, grid, ::VectorInvariantFormulation)

    Δxmin = minimum(grid.Δxᶜᵃᵃ)
    Δymin = minimum(grid.Δyᵃᶜᵃ)

    umax = maximum(abs, u)
    vmax = maximum(abs, v)

    return min(Δxmin / umax, Δymin / vmax)
end

cell_advection_timescale(model::ShallowWaterModel) =
    shallow_water_cell_advection_timescale(
        model.solution[1].data.parent,
        model.solution[2].data.parent,
        model.solution.h.data.parent,
        model.grid, model.formulation
        )


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept;
#
# import typedb.client.api.concept.thing.Thing;
# import typedb.client.api.concept.type.AttributeType;
# import typedb.client.api.concept.type.EntityType;
# import typedb.client.api.concept.type.RelationType;
# import typedb.client.api.concept.type.ThingType;
#
# import javax.annotation.CheckReturnValue;
# import javax.annotation.Nullable;
#
# public interface ConceptManager {
#
#     @CheckReturnValue
#     ThingType getRootThingType();
#
#     @CheckReturnValue
#     EntityType getRootEntityType();
#
#     @CheckReturnValue
#     RelationType getRootRelationType();
#
#     @CheckReturnValue
#     AttributeType getRootAttributeType();
#
#     @Nullable
#     @CheckReturnValue
#     ThingType getThingType(String label);
#
#     @Nullable
#     @CheckReturnValue
#     Thing getThing(String iid);
#
#     @Nullable
#     @CheckReturnValue
#     EntityType getEntityType(String label);
#
#     EntityType putEntityType(String label);
#
#     @Nullable
#     @CheckReturnValue
#     RelationType getRelationType(String label);
#
#     RelationType putRelationType(String label);
#
#     @Nullable
#     @CheckReturnValue
#     AttributeType getAttributeType(String label);
#
#     AttributeType putAttributeType(String label, AttributeType.ValueType valueType);
# }


text:
if VERSION <= v"1.0.5"
    # https://github.com/JuliaLang/julia/pull/29442
    _oneunit(::CartesianIndex{N}) where {N} = _oneunit(CartesianIndex{N})
    _oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N)))
else
    _oneunit = Base.oneunit
end


text:
using SSDPClient
using Test

str = ssdpquery("uuid", timeoutsecs = 15)
@test contains(str, "USN")

#str = ssdpquery("Roku", timeoutsecs = 30)
# @test contains(str, "Roku")


text:
using JSON3, StructTypes, BenchmarkTools

mutable struct T
    a::Int
    b::Float64
    c::String
    null::Nothing
    maybe::Union{Int, Nothing}
    T() = new()
end

JSON3.StructType(::Type{T}) = JSON3.Mutable()
StructTypes.StructType(::Type{T}) = StructTypes.Mutable()

@btime JSON3.read("{}", T)

@btime JSON3.read(b"""
{
    "a": 1,
    "b": 3.14,
    "c": "hey",
    "null": null,
    "maybe": 10
}""", T)

# master:   4.114 μs (62 allocations: 6.70 KiB)
# 0.1.13:    3.784 μs (16 allocations: 400 bytes)
# new:   950.759 ns (11 allocations: 432 bytes)


text:
using CompScienceMeshes
using BEAST
using Base.Test

m = meshrectangle(1.0, 1.0, 0.5)
b = meshsegment(1.0, 0.3, 3)

edges = skeleton(m, 1)

vps = cellpairs(m, edges)
rwg = raviartthomas(m, vps)

nt = ntrace(rwg, b)

#nt = BEAST.ntrace2(rwg, b)

n = 0
Σ = geometry(nt)
on_bnd = overlap_gpredicate(b)
for i in eachindex(rwg.fns)
    length(nt.fns[i]) == 1 || continue
    n += 1
    c = nt.fns[i][1].cellid
    edge = chart(Σ, Σ.faces[c])
    @test on_bnd(edge)
    @test nt.fns[i][1].refid == 1
    @test nt.fns[i][1].coeff == 2.0
end

@test n == 2

## test the scalar trace of a lgrange basis
using CompScienceMeshes
using BEAST
using Base.Test
m = meshrectangle(1.0, 1.0, 0.5)
b = meshsegment(1.0, 0.3, 3)

X = lagrangec0d1(m,b)

@test numfunctions(X) == 2
@test length(X.fns[1]) == 3
@test length(X.fns[2]) == 6

Y = strace(X,b)
A = find(length(f) for f in Y.fns)

@test length(A) == 1
@test length(Y.fns[1]) == 2
@test Y.fns[1][1].coeff ≈ 1
@test Y.fns[1][2].coeff ≈ 1


text:
"""
`merge`

Syntax:
    `merge_11 <varlist> using <DataFrame> , [options]`
    `merge_m1 <varlist> using <DataFrame> , [options]`
    `merge_1m <varlist> using <DataFrame> , [options]`

Description

Performs 1:1, 1:m, and m:1 merges.

Examples:
```julia
people = DataFrame(ID = [20, 40], Name = ["John Doe", "Jane Doe"])
jobs = DataFrame(ID = [20, 20, 60], Job = ["Lawyer", "Economist", "Astronaut"])
Douglass.set_active_df(:people)
d"merge_1m :ID using jobs"
```

Notable differences to Stata:
  - Note that the commands are e.g. `merge_m1` not `merge m:1`. There is an underscore and the colon is missing.
  - The commands currently do not produce any variable `_merge` and no summary statistics.
  - m:m merges are currently not allowed.

Notes:
    - Julia's DataFrames.jl has SQL-type joins, of which the merges here are but a small subset-- so much more
      powerful than this here.

"""
macro merge_11(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            merge_11 <varlist> using <DataFrame> , [options]
    """)
end

macro merge_11(t::Symbol,
    by::Nothing,
    sort::Nothing,
    arguments::Vector{Symbol},
    filter::Nothing,
    use::String,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@merge!($t, :one_to_one, $arguments, $use, $options )
        end
    )
end

macro merge_m1(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            merge_m1 <varlist> using <DataFrame> , [options]
    """)
end
macro merge_m1(t::Symbol,
    by::Nothing,
    sort::Nothing,
    arguments::Vector{Symbol},
    filter::Nothing,
    use::String,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@merge!($t, :m_to_one, $arguments, $use, $options )
        end
    )
end

macro merge_1m(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            merge_1m <varlist> using <DataFrame> , [options]
    """)
end
macro merge_1m(t::Symbol,
    by::Nothing,
    sort::Nothing,
    arguments::Vector{Symbol},
    filter::Nothing,
    use::String,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@merge!($t, :one_to_m, $arguments, $use, $options )
        end
    )
end

macro merge_mm(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            merge_mm <varlist> using <DataFrame> , [options]
    """)
end
macro merge_mm(t::Symbol,
    by::Nothing,
    sort::Nothing,
    arguments::Vector{Symbol},
    filter::Nothing,
    use::String,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@merge!($t, :m_to_m, $arguments, $use, $options )
        end
    )
end

# merge <type> <keys> using <rhs> , [options]
# Note: this could also be a function, really
macro merge!(t::Symbol, type::QuoteNode, keys::Vector{Symbol}, rhs::String, options::Union{Dict{String,Any}, Nothing})

    # parse rhs
    ex_rhs = Meta.parse(rhs)

    isa(ex_rhs, Symbol) || error("Merge operations expect a non-quoted symbol after `using`.")

    return esc(
        quote
            # make sure keys are present in both master and using
            # assert_vars_present($t, $keys)
            # assert_vars_present($rhs, $keys)
            for k in $keys
                (k ∈ propertynames($t)) || error("Variable $(k) not present in active DataFrame.")
            end
            for k in $keys
                (k ∈ propertynames($ex_rhs)) || error("Variable $(k) not present in RHS DataFrame.")
            end

            # outer = 1:1 or m:1
            # outer with lhs and rhs switched: 1:m

            # whenever the merge has a '1' in Stata, the keys must be uniquely identifying observations
            # this needs to be checked
            if ($type == :one_to_one) || ($type == :m_to_one)
                # check that we are restraining things correctly
                ($type == :one_to_one) && (Douglass.unique_obs($t, $keys) || error("Keys are not uniquely identifying observations in master DataFrame."))
                Douglass.unique_obs($ex_rhs, $keys) || error("Keys are not uniquely identifying observations in using DataFrame.")
                # do the merge
                $t = outerjoin($t, $ex_rhs, on = $keys, makeunique = false, indicator=nothing, validate=(false,false))

            elseif ($type == :one_to_m)
                Douglass.unique_obs($t, $keys) || error("Keys are not uniquely identifying observations in master DataFrame.")
                # do the merge
                $t = outerjoin($t, $ex_rhs, on = $keys, makeunique = false, indicator=nothing, validate=(false,false))
            elseif ($type == :m_to_m)
                error("m:m mergers are not allowed.")
            else
                error("Invalid merge type.")
            end
            $t
        end
    )
end

# This is just a helper to get the syntax right.
"""
`merge`

Syntax:
    `merge_11 <varlist> using <DataFrame> , [options]`
    `merge_m1 <varlist> using <DataFrame> , [options]`
    `merge_1m <varlist> using <DataFrame> , [options]`

"""
macro merge(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        `merge` does not exist.
        For merges, the syntax is:
            `merge_11 <varlist> using <DataFrame> , [options]`
            `merge_m1 <varlist> using <DataFrame> , [options]`
            `merge_1m <varlist> using <DataFrame> , [options]`
    """)
end


text:
const uri_pages = ["index", "home", "search", "main", "post", "homepage", "category", "register", "login", "faq", "about", "terms","privacy", "author"]
const uri_paths = ["app", "main", "wp-content", "search", "category", "tag", "categories", "tags", "blog", "posts", "list", "explore"]
const uri_extensions = [".html", ".html", ".html", ".htm", ".htm", ".php", ".php", ".jsp", ".asp"]

pixels()::Int = rand(1:1024)
user_name()::String = bothify(executor(data["faker"]["internet"]["user_name_formats"]))
domain_suffix()::String = executor(data["faker"]["internet"]["domain_suffix"])
safe_email(items::Tuple=(user_name(),"@example.", domain_suffix()))::String = reduce(string, items)
free_email_domain()::String = executor(data["faker"]["internet"]["free_email"])
free_email(items::Tuple=(user_name(),"@",free_email_domain()))::String = reduce(string, items)
domain_word()::String = lowercase(join(map( x -> length(x)>0 ? x[1] : "" , split(Faker.company(), [',', ' ','-'])),""))
domain_name(items::Tuple=(domain_word(),".",domain_suffix()))::String = reduce(string, items)
company_email(items::Tuple=(user_name(),"@",domain_name()))::String = reduce(string, items)

email()::String = executor(data["faker"]["internet"]["email_formats"])
url()::String = executor(data["faker"]["internet"]["url_formats"])

ipv4()::String = string(rand( 172:255),".",rand( 0:255),".",rand(0:255),".",rand(0:255))
ipv6()::String = join(map( x -> string(rand(4100:65535),base=16) ,1:8), ":")
mac_address()::String = join(map( x -> string(Base.hex(unsigned(rand( 0x00:0xff)),2,false)) ,1:6), ":")
uri_page()::String = rand(uri_pages)

uri_path(;deep::Int=0)::String =
    (if deep == 0
        deep = rand(1:3)
    end;
    return join(map( x -> string(rand(uri_paths)) , 0:deep), "/") )

uri_extension()::String = rand(uri_extensions)
uri()::String = executor(data["faker"]["internet"]["uri_formats"])
image_url()::String = executor(data["faker"]["internet"]["image_placeholder_services"])


text:
using DataFrames
using Base.Test

include("$(Pkg.dir())/MultiDimEquations/src/MultiDimEquations.jl")

# TEST 1: Testing both defVars()and the @meq macro using a single IndexedTable
df = CSV.read(IOBuffer("""
reg	prod	var	value
us	banana	production	10
us	banana	transfCoef	0.6
us	banana	trValues	2
us	apples	production	7
us	apples	transfCoef	0.7
us	apples	trValues	5
us	juice	production	NA
us	juice	transfCoef	NA
us	juice	trValues	NA
eu	banana	production	5
eu	banana	transfCoef	0.7
eu	banana	trValues	1
eu	apples	production	8
eu	apples	transfCoef	0.8
eu	apples	trValues	4
eu	juice	production	NA
eu	juice	transfCoef	NA
eu	juice	trValues    NA
"""), delim=" ", ignorerepeated=true)
variables =  vcat(unique(DataFrames.dropmissing(df)[:var]),["consumption"])
#defVars(variables,df;dfName="df",varNameCol="var", valueCol="value")
data = defVars(variables, df, tableName="data", varNameCol="var", valueCol="value")
products = ["banana","apples","juice"]
primPr   = products[1:2]
secPr    = [products[3]]
reg      = ["us","eu"]
# equivalent to [production!(sum(trValues_(r,pp) * transfCoef_(r,pp)  for pp in primPr), r, sp) for r in reg, sp in secPr]
@meq production!(r in reg, sp in secPr)   = sum(trValues_(r,pp) * transfCoef_(r,pp)  for pp in primPr)
@meq consumption!(r in reg, pp in primPr) = production_(r,pp) - trValues_(r,pp)
@meq consumption!(r in reg, sp in secPr)  = production_(r, sp)
totalConsumption = sum(consumption_(r,p) for r in reg, p in products)

@test totalConsumption == 26.6

# # TEST 2: Testing both the "old" defVarsDf()and the @meq macro using a single DataFrame as base
# df2 = wsv"""
# reg	prod	var	value
# us	banana	production	10
# us	banana	transfCoef	0.6
# us	banana	trValues	2
# us	apples	production	7
# us	apples	transfCoef	0.7
# us	apples	trValues	5
# us	juice	production	NA
# us	juice	transfCoef	NA
# us	juice	trValues	NA
# eu	banana	production	5
# eu	banana	transfCoef	0.7
# eu	banana	trValues	1
# eu	apples	production	8
# eu	apples	transfCoef	0.8
# eu	apples	trValues	4
# eu	juice	production	NA
# eu	juice	transfCoef	NA
# eu	juice	trValues    NA
# """
# variables =  vcat(unique(dropna(df[:var])),["consumption"])
# defVarsDf(variables,df2;dfName="df2",varNameCol="var", valueCol="value")
# products = ["banana","apples","juice"]
# primPr   = products[1:2]
# secPr    = [products[3]]
# reg      = ["us","eu"]
# # equivalent to [production!(sum(trValues_(r,pp) * transfCoef_(r,pp)  for pp in primPr), r, sp) for r in reg, sp in secPr]
# @meq production!(r in reg, sp in secPr)   = sum(trValues_(r,pp) * transfCoef_(r,pp)  for pp in primPr)
# @meq consumption!(r in reg, pp in primPr) = production_(r,pp) - trValues_(r,pp)
# @meq consumption!(r in reg, sp in secPr)  = production_(r, sp)
# totalConsumption = sum(consumption_(r,p) for r in reg, p in products)
#
# @test totalConsumption == 26.6

# TEST : Testing the @meq macro with individual IndexedTables
a = IndexedTable([["a","a","b","b"],[1,2,1,2]]...,[1,2,3,4])
dim1 = ["a","b"]
dim2 = [1,2]

@meq a[d1 in dim1,2] = a[d1,1]+3
tot = sum(a[d1,d2] for d1 in dim1, d2 in dim2)
@test tot == 14

# Test n: Fake test, this should not pass
# @test 1 == 2


text:
# arcs, circles, ellipses, curves, pie, sector, bezier

"""
    circle(x, y, r, action=:none)

Make a circle of radius `r` centered at `x`/`y`.

`action` is one of the actions applied by `do_action`, defaulting to `:none`. You can
also use `ellipse()` to draw circles and place them by their centerpoint.
"""
function circle(x, y, r, action=:none)
    if action != :path
        newpath()
    end
    Cairo.arc(get_current_cr(), x, y, r, 0, 2pi)
    do_action(action)
end

"""
    circle(pt, r, action=:none)

Make a circle centered at `pt`.
"""
circle(centerpoint::Point, r, action=:none) =
    circle(centerpoint.x, centerpoint.y, r, action)

"""
    circle(pt1::Point, pt2::Point, action=:none)

Make a circle that passes through two points that define the diameter:
"""
function circle(pt1::Point, pt2::Point, action=:none)
    center = midpoint(pt1, pt2)
    radius = distance(pt1, pt2)/2
    circle(center, radius, action)
end

"""
    circle(pt1::Point, pt2::Point, pt3::Point, action=:none)

Make a circle that passes through three points.
"""
function circle(pt1::Point, pt2::Point, pt3::Point, action=:none)
    center = midpoint(pt1, pt2)
    radius = distance(pt1, pt2)/2
    circle(center3pts(pt1, pt2, pt3)..., action)
end

"""
    center3pts(a::Point, b::Point, c::Point)

Find the radius and center point for three points lying on a circle.

returns `(centerpoint, radius)` of a circle.

If there's no such circle, the function returns `(Point(0, 0), 0)`.

If two of the points are the same, use `circle(pt1, pt2)` instead.
"""
function center3pts(p1::Point, p2::Point, p3::Point)
    norm2(p::Point) = (p.x)^2+(p.y)^2

    α1 = norm2(p3-p2)*(norm2(p2-p1)+norm2(p1-p3)-norm2(p3-p2))
    α2 = norm2(p1-p3)*(norm2(p3-p2)+norm2(p2-p1)-norm2(p1-p3))
    α3 = norm2(p2-p1)*(norm2(p1-p3)+norm2(p3-p2)-norm2(p2-p1))

    if α1+α2+α3 ≠ 0.0
        c = (α1*p1+α2*p2+α3*p3)/(α1+α2+α3)
        r = √(norm2(p1-c))
        return c, r
    else
        @warn "center3pts(): There are no circles which pass through $p1, $p2 and $p3."
        return (Point(0, 0), 0)
    end
end

"""
    ellipse(xc, yc, w, h, action=:none)

Make an ellipse, centered at `xc/yc`, fitting in a box of width `w` and height `h`.
"""
function ellipse(xc, yc, w, h, action=:none)
    x  = xc - w/2
    y  = yc - h/2
    # kappa = 4.0 * (sqrt(2.0) - 1.0) / 3.0
    kappa = 0.5522847498307936
    ox = (w / 2) * kappa  # control point offset horizontal
    oy = (h / 2) * kappa  # control point offset vertical
    xe = x + w            # x-end
    ye = y + h            # y-end
    xm = x + w / 2        # x-middle
    ym = y + h / 2        # y-middle
    move(x, ym)
    curve(x, ym - oy, xm - ox, y, xm, y)
    curve(xm + ox, y, xe, ym - oy, xe, ym)
    curve(xe, ym + oy, xm + ox, ye, xm, ye)
    curve(xm - ox, ye, x, ym + oy, x, ym)
    do_action(action)
end

"""
    ellipse(cpt, w, h, action=:none)

Make an ellipse, centered at point `c`, with width `w`, and height `h`.
"""
ellipse(c::Point, w, h, action=:none) = ellipse(c.x, c.y, w, h, action)

"""
    squircle(center::Point, hradius, vradius, action=:none;
        rt = 0.5, stepby = pi/40, vertices=false)

Make a squircle or superellipse (basically a rectangle with rounded corners).
Specify the center position, horizontal radius (distance from center to a side),
and vertical radius (distance from center to top or bottom):

The root (`rt`) option defaults to 0.5, and gives an intermediate shape. Values
less than 0.5 make the shape more rectangular. Values above make the shape more
round. The horizontal and vertical radii can be different.
"""
function squircle(center::Point, hradius, vradius, action=:none;
                  rt = 0.5,
                  vertices=false,
                  stepby = pi/40,
                  reversepath=false)
    points = Point[]
    for theta in 0:stepby:2pi
        xpos = center.x + ^(abs(cos(theta)), rt) * hradius * sign(cos(theta))
        ypos = center.y + ^(abs(sin(theta)), rt) * vradius * sign(sin(theta))
        push!(points, Point(xpos, ypos))
    end
    result = reversepath ? reverse(points) : points
    if !vertices
        poly(points, action, close=true, reversepath=reversepath)
    end
    return result
end

"""
    arc(xc, yc, radius, angle1, angle2, action=:none)

Add an arc to the current path from `angle1` to `angle2` going clockwise, centered
at xc, yc.

Angles are defined relative to the x-axis, positive clockwise.
"""
function arc(xc, yc, radius, angle1, angle2, action=:none)
    Cairo.arc(get_current_cr(), xc, yc, radius, angle1, angle2)
    do_action(action)
end

"""
    arc(centerpoint::Point, radius, angle1, angle2, action=:none)

Add an arc to the current path from `angle1` to `angle2` going clockwise, centered
at `centerpoint`.
"""
arc(centerpoint::Point, radius, angle1, angle2, action=:none) =
    arc(centerpoint.x, centerpoint.y, radius, angle1, angle2, action)

"""
    carc(xc, yc, radius, angle1, angle2, action=:none)

Add an arc to the current path from `angle1` to `angle2` going counterclockwise,
centered at `xc`/`yc`.

Angles are defined relative to the x-axis, positive clockwise.
"""
function carc(xc, yc, radius, angle1, angle2, action=:none)
    Cairo.arc_negative(get_current_cr(), xc, yc, radius, angle1, angle2)
    do_action(action)
end

"""
    carc(centerpoint::Point, radius, angle1, angle2, action=:none)

Add an arc centered at `centerpoint` to the current path from `angle1` to
`angle2`, going counterclockwise.
"""
carc(centerpoint::Point, radius, angle1, angle2, action=:none) =
    carc(centerpoint.x, centerpoint.y, radius, angle1, angle2, action)

"""
      arc2r(c1::Point, p2::Point, p3::Point, action=:none)

Add a circular arc centered at `c1` that starts at `p2` and ends at `p3`, going clockwise,
to the current path.

`c1`-`p2` really determines the radius. If `p3` doesn't lie on the circular path,
 it will be used only as an indication of the arc's length, rather than its position.
"""
function arc2r(c1::Point, p2::Point, p3::Point, action=:none)
    r = distance(c1, p2)
    startangle = atan(p2.y - c1.y, p2.x - c1.x)
    endangle   = atan(p3.y - c1.y, p3.x - c1.x)
    if endangle < startangle
        endangle = mod2pi(endangle + 2pi)
    end
    arc(c1, r, startangle, endangle, action)
end

"""
    carc2r(c1::Point, p2::Point, p3::Point, action=:none)

Add a circular arc centered at `c1` that starts at `p2` and ends at `p3`,
going counterclockwise, to the current path.

`c1`-`p2` really determines the radius. If `p3` doesn't lie on the circular
path, it will be used only as an indication of the arc's length, rather than its position.
"""
function carc2r(c1::Point, p2::Point, p3::Point, action=:none)
    r = distance(c1, p2)
    startangle = atan(p2.y - c1.y, p2.x - c1.x)
    endangle   = atan(p3.y - c1.y, p3.x - c1.x)
    if startangle < endangle
        startangle = mod2pi(startangle + 2pi)
    end
    carc(c1, r, startangle, endangle, action)
end

"""
    isarcclockwise(c::Point, A::Point, B::Point)

Return `true` if an arc centered at `c` going from `A` to `B` is clockwise.

If `c`, `A`, and `B` are collinear, then a hemispherical arc could be
either clockwise or not.
"""
function isarcclockwise(c::Point, A::Point, B::Point)
    a = A - c
    b = B - c
    return crossproduct(a, b) > 0
end

"""
    sector(centerpoint::Point, innerradius, outerradius, startangle, endangle, action:none)

Draw an annular sector centered at `centerpoint`.

"""
function sector(centerpoint::Point, innerradius::Real, outerradius::Real, startangle::Real,
                endangle::Real, action::Symbol=:none)
    (innerradius > outerradius) && throw(DomainError(outerradius, "outer radius must be larger than inner radius $(innerradius)"))
    gsave()
    translate(centerpoint)
    newpath()
    move(innerradius * cos(startangle), innerradius * sin(startangle))
    line(outerradius * cos(startangle), outerradius * sin(startangle))
    arc(0, 0, outerradius, startangle, endangle, :none)
    line(innerradius * cos(endangle), innerradius * sin(endangle))
    carc(0, 0, innerradius, endangle, startangle, :none)
    closepath()
    grestore()
    do_action(action)
end

"""
    sector(innerradius::Real, outerradius::Real, startangle::Real, endangle::Real,
       action::Symbol=:none)

Draw an annular sector centered at the origin.
"""
sector(innerradius::Real, outerradius::Real, startangle::Real, endangle::Real,
       action::Symbol=:none) =
    sector(O, innerradius, outerradius, startangle, endangle, action)

"""
    sector(centerpoint::Point, innerradius, outerradius, startangle, endangle,
        cornerradius, action:none)

Draw an annular sector with rounded corners, basically a bent sausage shape,
centered at `centerpoint`.

TODO: The results aren't 100% accurate at the moment. There are small
discontinuities where the curves join.

The cornerradius is reduced from the supplied value if neceesary to prevent overshoots.
"""
function sector(centerpoint::Point, innerradius::Real, outerradius::Real, startangle::Real,
                endangle::Real, cornerradius::Real, action::Symbol=:none)
    (innerradius > outerradius) && throw(DomainError(outerradius, "outer radius must be larger than inner radius $(innerradius)"))
    gsave()
    translate(centerpoint)
    # some work is done using polar coords to calculate the points

    # attempts to prevent pathological cases
    cornerradius = min(cornerradius, abs(outerradius-innerradius)/2)
    if endangle < startangle
        endangle = mod2pi(endangle + 2pi)
    end

    # TODO reduce given corner radius to prevent messes when spanning angle too small
    # 4 is a magic number
    while abs(endangle - startangle) < 4.0(atan(cornerradius, innerradius))
        cornerradius *= 0.75
    end

    # first inner corner
    alpha1 = asin(cornerradius/(innerradius+cornerradius))
    p1p2center = (innerradius + cornerradius, startangle + alpha1)
    p1 = (innerradius, startangle + alpha1)
    p2 = (innerradius + cornerradius, startangle)
    # first outer
    alpha2 = asin(cornerradius/(outerradius-cornerradius))
    p3p4center = (outerradius - cornerradius, startangle + alpha2)
    p3 = (outerradius - cornerradius, startangle)
    p4 = (outerradius, startangle + alpha2)
    # last outer
    p5p6center = (outerradius - cornerradius, endangle - alpha2)
    p5 = (outerradius, endangle - alpha2)
    p6 = (outerradius - cornerradius, endangle)
    # last inner
    p7p8center = (innerradius + cornerradius, endangle - alpha1)
    p7 = (innerradius + cornerradius, endangle)
    p8 = (innerradius, endangle - alpha1)

    # make path
    move(@polar(p1))
    newpath()
    # inner corner
    arc(@polar(p1p2center), cornerradius, slope(@polar(p1p2center), @polar(p1)),
        slope(@polar(p1p2center), @polar(p2)), :none)
    line(@polar(p3))
    # outer corner
    arc(@polar(p3p4center), cornerradius, slope(@polar(p3p4center), @polar(p3)),
        slope(@polar(p3p4center), @polar(p4)), :none)
    # outside arc
    arc(O, outerradius, slope(O, @polar(p4)), slope(O, @polar(p5)), :none)
    # last outside corner
    arc(@polar(p5p6center), cornerradius, slope(@polar(p5p6center), @polar(p5)),
        slope(@polar(p5p6center), @polar(p6)), :none)
    line(@polar(p7))
    # last inner corner
    arc(@polar(p7p8center), cornerradius, slope(@polar(p7p8center), @polar(p7)),
        slope(@polar(p7p8center), @polar(p8)), :none)
    s1, s2 = slope(O, @polar(p8)), slope(O, @polar(p1))
    if s1 < s2
        s2 = mod2pi(s2 + 2pi)
    end
    carc(O, innerradius, s1, s2, :none)
    closepath()
    do_action(action)
    grestore()
end

"""
    sector(innerradius::Real, outerradius::Real, startangle::Real, endangle::Real,
       cornerradius::Real, action::Symbol=:none)

Draw an annular sector with rounded corners, centered at the current origin.
"""
sector(innerradius::Real, outerradius::Real, startangle::Real, endangle::Real,
       cornerradius::Real, action::Symbol=:none) =
    sector(O, innerradius, outerradius, startangle, endangle, cornerradius, action)

"""
    pie(x, y, radius, startangle, endangle, action=:none)

Draw a pie shape centered at `x`/`y`. Angles start at the positive x-axis and
are measured clockwise.
"""
function pie(x, y, radius, startangle, endangle, action=:none)
    gsave()
    translate(x, y)
    newpath()
    move(0, 0)
    line(radius * cos(startangle), radius * sin(startangle))
    arc(0, 0, radius, startangle, endangle, :none)
    closepath()
    grestore()
    do_action(action)
end

"""
    pie(centerpoint, radius, startangle, endangle, action=:none)

Draw a pie shape centered at `centerpoint`.

Angles start at the positive x-axis and are measured clockwise.
"""
pie(centerpoint::Point, radius, startangle, endangle, action) =
    pie(centerpoint.x, centerpoint.y, radius, startangle, endangle, action)

"""
    pie(radius, startangle, endangle, action=:none)

Draw a pie shape centered at the origin

"""
pie(radius, startangle, endangle, action=:none) = pie(O, radius, startangle, endangle, action)

"""
    curve(x1, y1, x2, y2, x3, y3)
    curve(p1, p2, p3)

Add a Bézier curve.

The spline starts at the current position, finishing at `x3/y3` (`p3`),
following two control points `x1/y1` (`p1`) and `x2/y2` (`p2`).
"""
curve(x1, y1, x2, y2, x3, y3) = Cairo.curve_to(get_current_cr(), x1, y1, x2, y2, x3, y3)
curve(pt1, pt2, pt3)          = curve(pt1.x, pt1.y, pt2.x, pt2.y, pt3.x, pt3.y)

"""
    circlepath(center::Point, radius, action=:none;
        reversepath=false,
        kappa = 0.5522847498307936)

Draw a circle using Bézier curves.

The magic value, `kappa`, is `4.0 * (sqrt(2.0) - 1.0) / 3.0`.
"""
function circlepath(center::Point, radius, action=:none;
    reversepath=false,
    kappa = 0.5522847498307936)
    function northtoeast(center::Point, radius, kappa)
        curve(center.x + (radius * kappa), center.y + radius, center.x + radius,
        center.y + (radius * kappa), center.x + radius, center.y )
    end

    function easttosouth(center::Point, radius, kappa)
        curve(center.x + radius, center.y - (radius * kappa), center.x + (radius * kappa),
        center.y - radius, center.x, center.y - radius)
    end

    function southtowest(center::Point, radius, kappa)
        curve(center.x - (radius * kappa), center.y - radius, center.x - radius,
        center.y - (radius * kappa), center.x - radius, center.y)
    end

    function westtonorth(center::Point, radius, kappa)
        curve(center.x - radius, center.y + (radius * kappa), center.x - (radius * kappa),
        center.y + radius, center.x, center.y + radius)
    end

    function northtowest(center::Point, radius, kappa)
        curve(center.x - (radius * kappa), center.y + radius, center.x - radius,
        center.y + (radius * kappa), center.x - radius, center.y )
    end

    function westtosouth(center::Point, radius, kappa)
        curve(center.x - radius, center.y - (radius * kappa), center.x - (radius * kappa),
        center.y - radius, center.x, center.y - radius )
    end

    function southtoeast(center::Point, radius, kappa)
        curve(center.x + (radius * kappa), center.y - radius, center.x + radius,
        center.y - (radius * kappa), center.x + radius, center.y)
    end

    function easttonorth(center::Point, radius, kappa)
        curve(center.x + radius, center.y + (radius * kappa), center.x + (radius * kappa),
        center.y + radius, center.x, center.y + radius )
    end

    move(center.x, center.y + radius)
    if !reversepath
        northtoeast(center, radius, kappa)
        easttosouth(center, radius, kappa)
        southtowest(center, radius, kappa)
        westtonorth(center, radius, kappa)
    else
        northtowest(center, radius, kappa)
        westtosouth(center, radius, kappa)
        southtoeast(center, radius, kappa)
        easttonorth(center, radius, kappa)
    end
    do_action(action)
end

"""
    ellipse(focus1::Point, focus2::Point, k, action=:none;
            stepvalue=pi/100,
            vertices=false,
            reversepath=false)

Build a polygon approximation to an ellipse, given two points and a distance, `k`, which is
the sum of the distances to the focii of any points on the ellipse (or the shortest length
of string required to go from one focus to the perimeter and on to the other focus).
"""
function ellipse(focus1::Point, focus2::Point, k, action=:none;
                 stepvalue=pi/100,
                 vertices=false,
                 reversepath=false)
    a = k/2  # a = ellipse's major axis, the widest part
    cpoint = midpoint(focus1, focus2)
    dc = distance(focus1, cpoint)
    b = sqrt(abs(a^2 - dc^2)) # minor axis, hopefuly not 0
    phi = slope(focus1, focus2) # angle between the major axis and the x-axis
    points = Point[]
    drawing = false
    for t in 0:stepvalue:2pi
        xt = cpoint.x + a * cos(t) * cos(phi) - b * sin(t) * sin(phi)
        yt = cpoint.y + a * cos(t) * sin(phi) + b * sin(t) * cos(phi)
        push!(points, Point(xt, yt))
    end
    vertices ? points : poly(points, action, close=true, reversepath=reversepath)
end

"""
    ellipse(focus1::Point, focus2::Point, pt::Point, action=:none;
            stepvalue=pi/100,
            vertices=false,
            reversepath=false)

Build a polygon approximation to an ellipse, given two points and a point somewhere on the
ellipse.
"""
function ellipse(focus1::Point, focus2::Point, pt::Point, action=:none;
                 stepvalue=pi/100,
                 vertices=false,
                 reversepath=false)
    k = distance(focus1, pt) + distance(focus2, pt)
    ellipse(focus1, focus2, k, action, stepvalue=stepvalue,
        vertices=vertices, reversepath=reversepath)
end

"""
    hypotrochoid(R, r, d, action=:none;
            stepby=0.01,
            period=0.0,
            vertices=false)

Make a hypotrochoid with short line segments. (Like a Spirograph.) The curve is traced by a
point attached to a circle of radius `r` rolling around the inside  of a fixed circle of
radius `R`, where the point is a distance `d` from  the center of the interior circle.
Things get interesting if you supply non-integral values.

Special cases include the hypocycloid, if `d` = `r`, and an ellipse, if `R` = `2r`.

`stepby`, the angular step value, controls the amount of detail, ie the smoothness of the
polygon,

If `period` is not supplied, or 0, the lowest period is calculated for you.

The function can return a polygon (a list of points), or draw the points directly using
the supplied `action`. If the points are drawn, the function returns a tuple showing how
many points were drawn and what the period was (as a multiple of `pi`).
"""
function hypotrochoid(R, r, d, action=:none;
                      close    = true,
                      stepby   = 0.01,
                      period   = 0.0,
                      vertices = false)
    function nextposition(t)
        x = (R - r) * cos(t) + (d * cos(((R - r)/r) * t))
        y = (R - r) * sin(t) - (d * sin(((R - r)/r) * t))
        return Point(x, y)
    end
    # try to calculate the period exactly
    if isapprox(period, 0)
        period = 2pi * (r/gcd(convert(Int, floor(R)), convert(Int, floor(r))))
    end
    counter=1
    points=Point[]
    for t = 0:stepby:period
        push!(points, nextposition(t))
    end
    # don't repeat end point if it's more or less the same as the start point
    isapprox(points[1], points[end]) && pop!(points)
    vertices ? points : (poly(points, action, close=close); (length(points), period/pi))
end

"""
    epitrochoid(R, r, d, action=:none;
            stepby=0.01,
            period=0,
            vertices=false)

Make a epitrochoid with short line segments. (Like a Spirograph.) The curve is traced by a
point attached to a circle of radius `r` rolling around the outside of a fixed circle of
radius `R`, where the point is a distance `d` from the center of the circle.
Things get interesting if you supply non-integral values.

`stepby`, the angular step value, controls the amount of detail, ie the smoothness of the
polygon.

If `period` is not supplied, or 0, the lowest period is calculated for you.

The function can return a polygon (a list of points), or draw the points directly using
the supplied `action`. If the points are drawn, the function returns a tuple showing how
many points were drawn and what the period was (as a multiple of `pi`).
"""
function epitrochoid(R, r, d, action=:none;
                     close    = true,
                     stepby   = 0.01,
                     period   = 0,
                     vertices = false)
    function nextposition(t)
        x = (R + r) * cos(t) - (d * cos(((R - r)/r) * t))
        y = (R + r) * sin(t) - (d * sin(((R - r)/r) * t))
        return Point(x, y)
    end
    # try to calculate the period exactly
    if isapprox(period, 0)
        period = 2pi * (r/gcd(convert(Int, floor(R)), convert(Int, floor(r))))
    end
    counter=1
    points=Point[]
    for t = 0:stepby:period
        push!(points, nextposition(t))
    end
    # don't repeat end point if it's more or less the same as the start point
    isapprox(points[1], points[end]) && pop!(points)
    vertices ? points : (poly(points, action, close=close); (length(points), period/pi))
end

"""
    spiral(a, b, action::Symbol=:none;
                     stepby = 0.01,
                     period = 4pi,
                     vertices = false,
                     log=false)

Make a spiral. The two primary parameters `a` and `b` determine the start radius, and the
tightness.

For linear spirals (`log=false`), `b` values are:

    lituus: -2

    hyperbolic spiral: -1

    Archimedes' spiral: 1

    Fermat's spiral: 2

For logarithmic spirals (`log=true`):

    golden spiral: b = ln(phi)/ (pi/2) (about 0.30)

Values of `b` around 0.1 produce tighter, staircase-like spirals.
"""
function spiral(a, b, action::Symbol=:none;
                 stepby = 0.01,
                 period = 4pi,
                 vertices = false,
                 log=false)
    function nextpositionlog(t)
        ebt = exp(t * b)
        if abs(ebt) < 10^8 # arbitrary cutoff to avoid NaNs and Infs
            x = a * ebt * cos(t)
            y = a * ebt * sin(t)
            return Point(x, y)
        else
            return nothing
        end
    end
    function nextpositionlin(t)
        tpn = t ^ (1/b)
        if tpn < 10^8
            x = a * tpn * cos(t)
            y = a * tpn * sin(t)
            return Point(x, y)
        else
            return nothing
        end
    end
    log ? nextpos=nextpositionlog : nextpos = nextpositionlin
    points = Point[]
    for t = 0:stepby:period
        pt = nextpos(t)
        if isa(nextpos(t), Point)
            push!(points, pt)
        end
    end
    if vertices == false
        poly(points, action) # no close by default :)
    end
    return points
end

"""
    intersection2circles(pt1, r1, pt2, r2)

Find the area of intersection between two circles, the first centered at `pt1` with radius
`r1`, the second centered at `pt2` with radius `r2`.

If one circle is entirely within another, that circle's area is returned.
"""
function intersection2circles(pt1, rad1, pt2, rad2)
    # via casey and jùlio on slack
    # squared radii
    rr1, rr2 = rad1 * rad1, rad2 * rad2
    d = distance(pt1, pt2)

    # trivial cases
    d ≥ rad1 + rad2      && return 0.0
    d ≤ abs(rad2 - rad1) && return π * min(rr1, rr2)

    # First center point to the middle line
    a_distancecenterfirst = (rr1 - rr2 + (d^2)) / (2d)

    # Second centre point to the middle line
    b_distancecentersecond = d - a_distancecenterfirst

    # Half of the middle line
    h_height = sqrt(rr1 - a_distancecenterfirst^2)

    # central angle for the first circle
    alpha = mod2pi(atan(h_height, a_distancecenterfirst) * 2.0 + 2π)

    # Central angle for the second circle
    beta = mod2pi(atan(h_height, b_distancecentersecond) * 2.0 + 2π)

    #  Area of the first circular segment
    A1 = rr1 / 2.0 * (alpha - sin(alpha))
    # Area of the second circular segment
    A2 = rr2 / 2.0 * (beta - sin(beta))
    return A1 + A2
end

"""
    intersectioncirclecircle(cp1, r1, cp2, r2)

Find the two points where two circles intersect, if they do. The first circle is centered
at `cp1` with radius `r1`, and the second is centered at `cp1` with radius `r1`.

Returns

    (flag, ip1, ip2)

where `flag` is a Boolean `true` if the circles intersect at the points `ip1` and `ip2`. If
the circles don't intersect at all, or one is completely inside the other, `flag` is `false`
and the points are both Point(0, 0).

Use `intersection2circles()` to find the area of two overlapping circles.

In the pure world of maths, it must be possible that two circles 'kissing' only have a
single intersection point. At present, this unromantic function reports that two kissing
circles have no intersection points.
"""
function intersectioncirclecircle(cp1, r1, cp2, r2)
    r1squared = r1^2
    r2squared = r2^2
    d = distance(cp1, cp2)
    if d > (r2 + r1) # circles do not overlap
        return (false, O, O)
    elseif (d <= abs(r1 - r2)) && (r1 >= r2)
        # second circle is completely inside first circle
        return (false, O, O)
    elseif (d <= abs(r1 - r2)) && (r1 < r2)
        # first circle is completely inside second circle
        return (false, O, O)
    end
    a = (r1squared - r2squared + d^2) / 2d
    h = sqrt(r1squared - a^2)
    p0 = cp1 + a * (cp2 - cp1)/d
    p3 = Point(
            p0.x + h * (cp2.y - cp1.y )/d,
            p0.y - h * (cp2.x - cp1.x )/d)
    p4 = Point(
            p0.x - h * (cp2.y - cp1.y )/d,
            p0.y + h * (cp2.x - cp1.x )/d)
    return (true, p3, p4)
end

"""
    circlepointtangent(through::Point, radius, targetcenter::Point, targetradius)

Find the centers of up to two circles of radius `radius` that pass through point
`through` and are tangential to a circle that has radius `targetradius` and
center `targetcenter`.

This function returns a tuple:

* (0, O, O)      - no circles exist

* (1, pt1, O)    - 1 circle exists, centered at pt1

* (2, pt1, pt2)  - 2 circles exist, with centers at pt1 and pt2

(The O are just dummy points so that three values are always returned.)
"""
function circlepointtangent(through::Point, radius, targetcenter::Point, targetradius)
    distx = targetcenter.x - through.x
    disty = targetcenter.y - through.y
    dsq = distance(through, targetcenter)^2
    if isless(dsq, 10e-6) # coincident
        return (0, O, O)
    else
        sqinv=0.5/dsq
        s = dsq - ((2radius + targetradius) * targetradius)
        root = 4(radius^2) * dsq - s^2
        s *= sqinv
        if isless(dsq, 0.0) # no center possible
            return (0, O, O)
        else
            if isless(root, 10e-6) # only one circle possible
                x = through.x + distx * s
                y = through.y + disty * s
                if isless(abs(distance(through, Point(x, y)) - radius), 10e-6)
                    return (1, Point(x, y), O)
                else
                    return (0, O, O)
                end
            else # two circles are possible
                root = sqrt(root) * sqinv
                xconst = through.x + distx * s
                yconst = through.y + disty * s
                xvar = disty * root
                yvar = distx * root
            return (2, Point(xconst - xvar, yconst + yvar), Point(xconst + xvar, yconst - yvar))
            end
        end
    end
end

"""
    circletangent2circles(radius, circle1center::Point, circle1radius, circle2center::Point, circle2radius)

Find the centers of up to two circles of radius `radius` that are tangent to the
two circles defined by `circle1...` and `circle2...`. These two circles can
overlap, but one can't be inside the other.

* (0, O, O)      - no such circles exist

* (1, pt1, O)    - 1 circle exists, centered at pt1

* (2, pt1, pt2)  - 2 circles exist, with centers at pt1 and pt2

(The O are just dummy points so that three values are always returned.)
"""
function circletangent2circles(radius, circle1center::Point, circle1radius, circle2center::Point, circle2radius)
    modradius1 = radius + circle1radius
    modradius2 = circle2radius - circle1radius
    return circlepointtangent(circle1center, modradius1, circle2center, modradius2)
end

"""
    arc2sagitta(p1::Point, p2::Point, s, action=:none)

Make a clockwise arc starting at `p1` and ending at `p2` that reaches a height of `s`, the sagitta, at the middle. Might append to current path...

Return tuple of the center point and the radius of the arc.
"""
function arc2sagitta(p1::Point, p2::Point, s, action=:none)
    if isapprox(s, 0.0)
        throw(error("Height of arc $s should be greater than 0.0"))
    end
    l = distance(p1, p2)/2
    r = (s^2 + l^2) / 2s
    flag, ip1, ip2 = intersectioncirclecircle(p1, r, p2, r)
    if flag
        if r <= s
            arc2r(ip1, p1, p2, action)
            result = (ip1, r)
        else
            arc2r(ip2, p1, p2, action)
            result = (ip2, r)
        end
    else
        result = (O, 0.0)
    end
    return result
end

"""
    carc2sagitta(p1::Point, p2::Point, s, action=:none)

Make a counterclockwise arc starting at `p1` and ending at `p2` that reaches a height of `s`, the sagitta, at the middle. Might append to current path...

Return tuple of center point and radius of arc.
"""
function carc2sagitta(p1::Point, p2::Point, s, action=:none)
    if isapprox(s, 0.0)
        throw(error("Height of arc $s should be greater than 0.0"))
    end
    l = distance(p1, p2)/2
    r = (s^2 + l^2) / 2s
    flag, ip1, ip2 = intersectioncirclecircle(p1, r, p2, r)
    if flag
        if r <= s
            carc2r(ip2, p1, p2, action)
            result = (ip1, r)
        else
            carc2r(ip1, p1, p2, action)
            result = (ip2, r)
        end
    else
        result = (O, 0.0)
    end
    return result
end

"""
    circlecircleoutertangents(cpt1::Point, r1, cpt2::Point, r2)

Return four points, `p1`, `p2, `p3`, `p4`, where a line
through `p1` and `p2`, and a line through `p3` and `p4`, form the outer
tangents to the circles defined by `cpt1/r1` and `cpt2/r2`.

Returns four identical points (`O`) if one of the circles lies inside the other.
"""
function circlecircleoutertangents(cpt1::Point, r1, cpt2::Point, r2)
    dx = cpt2.x - cpt1.x
    dy = cpt2.y - cpt1.y
    dist = distance(cpt1, cpt2)

    if dist <= abs(r2 - r1)
        # the circles overlap, no tangents
        @debug "circlecircleoutertangents: circles overlap"
        return O, O, O, O
    end

    # rotation from the x-axis.
    angle1 = atan(dy, dx)

    # relative angle of the normals
    angle2 = acos((r1 - r2)/dist)

    pt1 = Point(cpt1.x + r1 * cos(angle1 + angle2),
        cpt1.y + r1 * sin(angle1 + angle2))

    pt2 = Point(cpt2.x + r2 * cos(angle1 + angle2),
        cpt2.y + r2 * sin(angle1 + angle2))

    pt3 = Point(cpt1.x + r1 * cos(angle1 - angle2),
        cpt1.y + r1 * sin(angle1 - angle2))

    pt4 = Point(cpt2.x + r2 * cos(angle1 - angle2),
        cpt2.y + r2 * sin(angle1 - angle2))

    return pt1, pt2, pt3, pt4
end

"""
    pointcircletangent(point::Point, circlecenter::Point, circleradius)

Find the two points on a circle that lie on tangent lines
passing through an external point.

If both points are O, the external point is inside the
circle, and the result is `(O, O)`.
"""
function pointcircletangent(point::Point, circlecenter::Point, circleradius)
    d = distance(circlecenter, point)
    if d <= circleradius # isinside
        @debug "pointcircletangent: external point is inside circle"
        # for type stability, return two points anyway?
        return O, O
    end
    dx, dy = point.x - circlecenter.x, point.y - circlecenter.y
    dxr, dyr = -dy, dx # rotate
    ρ = circleradius/d
    ad = ρ^2
    bd = ρ * sqrt(1 - ρ^2)
    tangentpoint1x = circlecenter.x + (ad * dx) + (bd * dxr)
    tangentpoint1y = circlecenter.y + (ad * dy) + (bd * dyr)
    tangentpoint2x = circlecenter.x + (ad * dx) - (bd * dxr)
    tangentpoint2y = circlecenter.y + (ad * dy) - (bd * dyr)
    return Point(tangentpoint1x, tangentpoint1y), Point(tangentpoint2x, tangentpoint2y)
end

"""
circlecircleinnertangents(circle1center::Point, circle1radius, circle2center::Point, circle2radius)

Find the inner tangents of two circles. These are tangent lines that
cross as they skim past one circle and touch the other.

Returns the four points: tangentpoint1 on circle 1,
tangentpoint1 on circle2, tangentpoint2 on circle 1,
tangentpoint2 on circle2.

Returns `(O, O, O, O)` if inner tangents can't be found (eg when the circles overlap).

Use `circlecircleoutertangents()` to find the outer tangents.
"""
function circlecircleinnertangents(circle1center::Point, circle1radius, circle2center::Point, circle2radius)
    # point where the inner tangents will cross
    ip = (circle2center * circle1radius + circle1center * circle2radius) / (circle1radius + circle2radius)
    d = distance(circle1center, circle2center)
    # circles overlap?
    d <= (circle1radius + circle2radius) && return O, O, O, O
    tp1, tp2 = pointcircletangent(ip, circle1center, circle1radius)
    tp3, tp4 = pointcircletangent(ip, circle2center, circle2radius)
    # double check tp1 and tp3 are collinear, just in case the order was wrong
    if slope(tp1, ip) ≈ slope(ip, tp3)
        return tp1, tp3, tp2, tp4
    elseif slope(tp1, ip) ≈ slope(ip, tp4)
        return tp1, tp4, tp2, tp3
    end
end

"""
    ellipseinquad(qgon, action=:none)

Calculate a Bézier-based ellipse that fits inside the
quadrilateral `qgon`, an array of with at least four Points,
then apply `action`.

Returns `ellipsecenter, ellipsesemimajor, ellipsesemiminor,
ellipseangle`:

`ellipsecenter` the ellipse center

`ellipsesemimajor` ellipse semimajor axis

`ellipsesemiminor` ellipse semiminor axis

`ellipseangle` ellipse rotation

The function returns `O, 0, 0, 0` if a suitable ellipse
can't be found. (The qgon is probably not a convex polygon.)

### References

http://faculty.mae.carleton.ca/John_Hayes/Papers/InscribingEllipse.pdf
"""
function ellipseinquad(qgon, action=:none)
    p1, p2, p3, p4 = qgon
    length(unique([p1, p2, p3, p4])) != 4 && throw(error("ellipseinquad() the 4 points must be different: $(qgon)"))
    if !ispolyclockwise(qgon)
        reverse!(qgon)
    end
    if !ispolyconvex(qgon)
        return O, 0, 0, 0
    end
    r = [p1.x, p1.y, 1]
    s = [p2.x, p2.y, 1]
    t = [p3.x, p3.y, 1]
    u = [p4.x, p4.y, 1]

    # Create the 'projective collineation' matrix
    # whatever that is
    e1 = s[1] * t[1] * u[2] - r[1] * t[1] * u[2] - s[1] * t[2] * u[1] + r[1] * t[2] * u[1] - r[1] * s[2] * u[1] + r[2] * s[1] * u[1] + r[1] * s[2] * t[1] - r[2] * s[1] * t[1]
    e2 = r[1] * t[1] * u[2] - r[1] * s[1] * u[2] - s[1] * t[2] * u[1] + s[2] * t[1] * u[1] - r[2] * t[1] * u[1] + r[2] * s[1] * u[1] + r[1] * s[1] * t[2] - r[1] * s[2] * t[1]
    e3 = s[1] * t[1] * u[2] - r[1] * s[1] * u[2] - r[1] * t[2] * u[1] - s[2] * t[1] * u[1] + r[2] * t[1] * u[1] + r[1] * s[2] * u[1] + r[1] * s[1] * t[2] - r[2] * s[1] * t[1]
    e4 = s[2] * t[1] * u[2] - r[2] * t[1] * u[2] - r[1] * s[2] * u[2] + r[2] * s[1] * u[2] - s[2] * t[2] * u[1] + r[2] * t[2] * u[1] + r[1] * s[2] * t[2] - r[2] * s[1] * t[2]
    e5 = -s[1] * t[2] * u[2] + r[1] * t[2] * u[2] + s[2] * t[1] * u[2] - r[1] * s[2] * u[2] - r[2] * t[2] * u[1] + r[2] * s[2] * u[1] + r[2] * s[1] * t[2] - r[2] * s[2] * t[1]
    e6 = s[1] * t[2] * u[2] - r[1] * t[2] * u[2] + r[2] * t[1] * u[2] - r[2] * s[1] * u[2] - s[2] * t[2] * u[1] + r[2] * s[2] * u[1] + r[1] * s[2] * t[2] - r[2] * s[2] * t[1]
    e7 = s[1] * u[2] - r[1] * u[2] - s[2] * u[1] + r[2] * u[1] - s[1] * t[2] + r[1] * t[2] + s[2] * t[1] - r[2] * t[1]
    e8 = t[1] * u[2] - s[1] * u[2] - t[2] * u[1] + s[2] * u[1] + r[1] * t[2] - r[2] * t[1] - r[1] * s[2] + r[2] * s[1]
    e9 = t[1] * u[2] - r[1] * u[2] - t[2] * u[1] + r[2] * u[1] + s[1] * t[2] - s[2] * t[1] + r[1] * s[2] - r[2] * s[1]
    mat = [e1 e2 e3 ; e4 e5 e6 ; e7 e8 e9]
    matinv = inv(mat)
    #  coefficients of the ellipse equation ax2 + 2bxy + cy2 + 2dx + 2fy + g = 0
    a = matinv[1,1] * matinv[1,1] + matinv[2,1] * matinv[2,1] - matinv[3,1] * matinv[3,1]
    b = matinv[1,1] * matinv[1,2] + matinv[2,1] * matinv[2,2] - matinv[3,1] * matinv[3,2]
    c = matinv[1,2] * matinv[1,2] + matinv[2,2] * matinv[2,2] - matinv[3,2] * matinv[3,2]
    d = matinv[1,1] * matinv[1,3] + matinv[2,1] * matinv[2,3] - matinv[3,1] * matinv[3,3]
    f = matinv[1,2] * matinv[1,3] + matinv[2,2] * matinv[2,3] - matinv[3,2] * matinv[3,3]
    g = matinv[1,3] * matinv[1,3] + matinv[2,3] * matinv[2,3] - matinv[3,3] * matinv[3,3]

    #  ellipse centre
    b2ac = b*b - a*c
    ellipsecenter = Point((c*d - b*f) / b2ac, (a*f - b*d) / b2ac)

    #  axes
    k = 2(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g) / b2ac
    ac24b2 = sqrt((a - c) * (a - c) + 4b*b)
    ellipsesemimajor = sqrt(abs(k / abs( ac24b2 - a - c)))
    ellipsesemiminor = sqrt(abs(k / abs(-ac24b2 - a - c)))

    #  rotation angle
    if b == 0
        ellipseangle = (a < c) ? 0 : π/2
    else
        ellipseangle = (a < c) ? 0.5 * atan(2b / (a-c)) : π/2 + 0.5 * atan(2b / (a-c))
    end

    # finally
    @layer begin
        translate(ellipsecenter)
        rotate(ellipseangle)
        ellipse(O, 2ellipsesemimajor, 2ellipsesemiminor, action)
    end
    return ellipsecenter, ellipsesemimajor, ellipsesemiminor, ellipseangle
end

"""
    crescent(pos, innerradius, outeradius, action=nothing;
        vertices=false,
        reversepath=false,
        steps = 30)

Create a crescent-shaped polygon, aligned with the current x-axis.
If the inner radius is 0, you'll get a semicircle.

See also ```crescent(pos1, innerradius, pos2, outeradius...)```.

## Examples

Create a filled crescent shape with outer radius of 200, inner radius of 130.
```
crescent(O, 130, 200, :fill)
```

Create a stroked crescent shape; the inner radius of 0 produces a semicircle.

```
crescent(O, 0, 200, :stroke)
```
"""
function crescent(pos::Point, innerradius::Real, outerradius::Real, action=nothing;
        vertices=false,
        reversepath=false,
        steps = 30)
    Δ = π/steps
    ptlist = Point[]
    move(pos + Point(innerradius * sin(0), outerradius * cos(0)))
    push!(ptlist, currentpoint())
    for i in 1:steps
        pt = pos + Point(innerradius * sin(Δ * i), outerradius * cos(Δ * i))
        push!(ptlist, pt)
    end
    for i in 1:steps
        pt = pos + Point(outerradius * sin(π - (Δ * i)), outerradius * cos(π - (Δ * i)))
        push!(ptlist, pt)
    end
    if !vertices
        poly(ptlist, action, close=true, reversepath=reversepath)
    end
    return ptlist
end

"""
crescent(cp1, r1, cp2, r2, action=nothing;
            vertices=false,
            reversepath=false)

Create a crescent-shaped polygon, aligned with the current
x-axis, by finding the intersection of two circles.
The two center positions should be different.

See also ```crescent(point, innerradius, outeradius...)```.

## Examples

Create a filled crescent shape from two circles.

```
crescent(O, 100, O + (60, 0), 150, :fill)
```
"""
function crescent(cp1::Point, r1::Real, cp2::Point, r2::Real, action=nothing;
        vertices=false,
        reversepath=false)
    flag, ip1, ip2 = intersectioncirclecircle(cp1, r1, cp2, r2)
    if flag
        move(ip2)
        carc(cp2, r2, slope(cp2, ip2), slope(cp2, ip1), :path)
        arc(cp1, r1, slope(cp1, ip1), slope(cp1, ip2), :path)
        ptlist = pathtopoly()[1]
        if !vertices
            poly(ptlist, action, close=true, reversepath=reversepath)
        end
        return ptlist
    else
        @info "crescent(): the two arcs won't intersect "
    end
end
# eof


text:
##
# Searching
##########################################################

function search_part_idx(file, terms::Array, mapto_url::Bool, limit::Int)
    results = search_part_idx(file, terms, mapto_url)
    if limit > 0
        while length(results) > limit
            pop!(results)
        end
    end
    results
end

function search_part_idx(file, terms::Array, mapto_url::Bool)
    results = search_part_idx(file, terms)

    if mapto_url
        id_map_file = joinpath(id_to_doc_location, basename(file))
        id_to_doc = cached_get(id_map_file, ()->as_deserialized(openable(id_map_file)))
        result_docs = map(id->get(id_to_doc, id, ""), results)
        return Set(filter(x->!isempty(x), result_docs))
    else
        return Set(results)
    end
end

function search_part_idx(file, terms::Array)
    part_idx = cached_get(file, ()->as_deserialized(file))
    results = IntSet()
    for term in terms
        union!(results, get(part_idx, term, []))
    end
    results
end

function search_index(terms::String; mapto_url::Bool=true, limit::Int=10)
    sd = StringDocument(terms)
    td = preprocess(sd)
    terms = tokens(td)
    terms = filter(tok->!isempty(tok), terms)

    # The Blocks library provides many ways to divide work among processors.
    # Here, we call it on the part_idx_location directory, which tells it to
    # divide up the files in that directory.
    # The call localpart(master_idx) returns a list of files to be handled
    # by the calling processor.
    master_idx = cached_get(part_idx_location, ()->Block(openable(part_idx_location), false, 2))

    results = @parallel union for i in 1:nworkers()
        local_files = Any[]
        for b in localpart(master_idx)
            append!(local_files, b)
        end
        search_results = map(file->search_part_idx(file, terms, mapto_url, limit), local_files)
        isempty(search_results) ? Any[] : reduce(union, search_results)
    end
    if limit > 0
        while length(results) > limit
            pop!(results)
        end
    end
    collect(results)
end


text:
# Sensors and Actuators B 149 (2010) 189-193

# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Define beam
    λ = LinRange(400,1000,500) # wavelength range [nm]
    θ = [0.] # angle of incidence [degrees]
    beam = PlaneWave(λ,θ)
    # Define layers with d=1/4 by default
    p1 = 0.86
    p2 = 0.54
    n = [RIdb.air(beam.λ) RIdb.silicon(beam.λ)]
    l0 = LayerTMMO(RIdb.air(beam.λ)) # incident medium
    l1 = LayerTMMO(RI.looyenga([p1 1-p1],n); type=:OT)
    l2 = LayerTMMO(RI.looyenga([p2 1-p2],n); type=:OT)
    l3 = LayerTMMO(RIdb.glass(beam.λ./1e3)) # emergent medium
    layers = vec([l0 l1 l2 l1 l2 l1 l2 l1 l2 repeat([l1], 1, 10) l2 l1 l2 l1 l2 l1 l2 l1 l3])
    # Reference wavelength
    λ0 = 730.0
    # call main script
    sol = tmm_optics(beam, layers; λ0=λ0, emfflag=true, h=10)
    return sol
end

sol = main()

# plot the R, T and A spectra
plot(Spectrum1D(),
    sol.beam.λ, [sol.Spectra.Rp, sol.Spectra.Tp, 1.0.-(sol.Spectra.Rp.+sol.Spectra.Tp)],
    label=["Reflectance" "Transmittance" "Absorbance"],
    line=([:solid :dash :dashdot]),
    xlims=(sol.beam.λ[1], sol.Beam.λ[end]),
    yaxis=("Transmittance", (0.,1.0)),
)
gui()

plot(EMF2D(), sol.beam.λ, sol.Misc.ℓ, sol.Field.emfp[:,1,:])
gui()

# plot the refractive index profile
plot(RIprofile(), sol)
gui()


text:
# This file was generated, do not modify it. # hide
using Random
Random.seed!(1) # hide
@show randn(2)


text:
# example of a function that have a result type
function sum(x, y)::Float64
    return x + y
end

# two use cases
a = sum(1, 2)
b = sum(1.2, 2)

println(a)
println(b)

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
"""
    importASCII(file::AbstractString [, outfile])

Imports an ESRI Ascii grid to a .grd file and returns a new VerySimpleRaster
object. Specify `outfile` to keep the newly created file.
"""
function importASCII(x::AbstractString, outfile = "")
    file = open(x, "r")
    nc = parse(Int, match(r"NCOLS (.+)", readline(file)).captures[1])
    nr = parse(Int, match(r"NROWS (.+)", readline(file)).captures[1])
    xll = parse(Float64, match(r"XLLCORNER (.+)", readline(file)).captures[1])
    yll = parse(Float64, match(r"YLLCORNER (.+)", readline(file)).captures[1])
    cell = parse(Float64, match(r"CELLSIZE (.+)", readline(file)).captures[1])
    NA = parse(Float64, match(r"NODATA_value (.+)", readline(file)).captures[1])
    tmp = Vector{Float64}(undef, nc)

    outfile == "" && (outfile = tempname())
    outfile[end-3:end] ∈ (".gri", "grd") && (outfile = outfile[1:end-4])
    IO = open(outfile*".gri", "w")

    for row in nr:-1:1
        tmp .= parse.(Float64, split(readline(file), " "))
        for r in tmp
            write(IO, r)
        end
    end
    close(file)
    close(IO)
    mat = Mmap.mmap(outfile*".gri", Matrix{Float64}, (nc, nr))
    vsr = VerySimpleRaster(mat, NA, xll, yll, cell, cell, "", outfile*".gri")
    write_grd_header(outfile*".grd", vsr)
    vsr
end


text:
using CUDAdrv

"""
    empiricalbandwidth(nbytes=2*1024^3; devicenumber=0, numtests=10)

Compute the emperical bandwidth in GB/s of a CUDA device using `nbytes` of
memory.  The device to test can be slected with `devicenumber` and the
bandwidth is an average of `ntests`.
"""
function empiricalbandwidth(nbytes=1024^2; devicenumber=0, ntests=10, pinned=false)
    dev = CuDevice(devicenumber)
    ctx = CuContext(dev)
    #stm = CuStream()
    stm = CuStream(CUDAdrv.STREAM_NON_BLOCKING)

    d = rand(Char, nbytes)
    a = pinned ? Mem.alloc(Mem.Host, nbytes) : pointer(d)
    b = Mem.alloc(Mem.Device, nbytes)

    Mem.copy!(a, b, nbytes)
    Mem.copy!(b, a, nbytes)

    t0, t1 = CuEvent(), CuEvent()
    record(t0, stm)
    for n = 1:ntests
        Mem.copy!(a, b, nbytes, async=true, stream=stm)
    end
    record(t1, stm)
    synchronize(t1)
    t_dtoh = elapsed(t0, t1)
    bandwidth_dtoh = nbytes*ntests/(t_dtoh*1e9)

    t0, t1 = CuEvent(), CuEvent()
    record(t0, stm)
    for n = 1:ntests
        Mem.copy!(b, a, nbytes, async=true, stream=stm)
    end
    record(t1, stm)
    synchronize(t1)
    t_htod = elapsed(t0, t1)
    bandwidth_htod = nbytes*ntests/(t_htod*1e9)

    pinned && Mem.free(a)
    Mem.free(b)

    (bandwidth_dtoh, bandwidth_htod)
end


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

Reads (length,azimuth) pairs and plot a windrose diagram.

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

Parameters
----------

- $(GMT.opt_J)
- **A** : **sector** : -- Str or number --

	Gives the sector width in degrees for sector and rose diagram.
	[`-A`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#a)
- $(GMT.opt_B)
- **C** : **vectors** : -- Str --

	Plot vectors showing the principal directions given in the mode_file file.
	[`-C`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#c)
- **D** : **shift** : -- Bool or [] --

	Shift sectors so that they are centered on the bin interval (e.g., first sector is centered on 0 degrees).
	[`-D`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#d)
- **F** : **no_scale** : -- Bool or [] --

	Do not draw the scale length bar [Default plots scale in lower right corner].
	[`-F`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#f)
- **G** : **fill** : -- Number or Str --

	Selects shade, color or pattern for filling the sectors [Default is no fill].
	[`-G`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#g)
- **I** : **inquire** : -- Bool or [] --

	Inquire. Computes statistics needed to specify a useful -R. No plot is generated.
	[`-I`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#i)
- **L** : **pen** : -- Number or Str --

	Draw bar outline using the specified pen thickness. [Default is no outline].
	[`-L`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#l)
- **M** : -- Bool or [] --

	Used with -C to modify vector parameters.
	[`-M`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#m)
- $(GMT.opt_P)
- **Q** : **alpha** : -- Str or [] --

	Sets the confidence level used to determine if the mean resultant is significant.
	[`-Q`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#q)
- $(GMT.opt_R)
- **S** : **radius** : -- Bool or [] --

	Specifies radius of plotted circle (append a unit from c|i|p).
	[`-S`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#s)
- **T** : -- Bool or [] --

	Specifies that the input data are orientation data (i.e., have a 180 degree ambiguity)
	instead of true 0-360 degree directions [Default].
	[`-T`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#t)
- **W** : **pen** : -- Str or tuple --

	Set pen attributes for sector outline or rose plot. [Default is no outline].
	[`-W`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#w)
- **Z** : **scale** : -- Str --

	Multiply the data radii by scale.
	[`-Z`](http://gmt.soest.hawaii.edu/doc/latest/psrose.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function rose(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	arg2 = []		# May be needed if GMTcpt type is sent in via C
	N_args = isempty_(arg1) ? 0 : 1

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

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, "")
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di)

	cmd = add_opt(cmd, 'A', d, [:A :sector])
	cmd = add_opt(cmd, 'C', d, [:C :vectors])
	cmd = add_opt(cmd, 'D', d, [:D :shift])
	cmd = add_opt(cmd, 'F', d, [:F :no_scale])
	cmd = add_opt(cmd, 'G', d, [:G :fill])
	cmd = add_opt(cmd, 'I', d, [:I :inquire])
	cmd = cmd * opt_pen(d, 'L', [:L :pen])
	cmd = add_opt(cmd, 'M', d, [:M])
	cmd = add_opt(cmd, 'Q', d, [:Q :alpha])
	cmd = add_opt(cmd, 'S', d, [:S :radius])
	cmd = add_opt(cmd, 'T', d, [:T :radius])
	cmd = cmd * opt_pen(d, 'W', [:W :pen])
	cmd = add_opt(cmd, 'Z', d, [:Z :scale])

	cmd = finish_PS(d, cmd, output, K, O)
	return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "psrose", arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
rose!(cmd0::String="", arg1=[]; K=true, O=true,  first=false, kw...) =
	rose(cmd0, arg1; K=K, O=O,  first=first, kw...)

rose(arg1=[], cmd0::String=""; K=false, O=false,  first=true, kw...) =
	rose(cmd0, arg1; K=K, O=O,  first=first, kw...)

rose!(arg1=[], cmd0::String=""; K=true, O=true,  first=false, kw...) =
	rose(cmd0, arg1; K=K, O=O,  first=first, kw...)

psrose  = rose 			# Alias
psrose! = rose!			# Alias


text:
# recoupling function
# two options for the type:
#   - either all arguments are Float64 or energy variables are complex
function change_basis_3from1(σ1, cosθ1, ϕ1, cosθ23, ϕ23,
        m1sq, m2sq, m3sq, s)
        # calculate σ3 in (23) frame
        m1 = sqrt(m1sq); sqrts = sqrt(s)
        sqrtλσ1 = sqrt((sqrts-m1)^2-σ1)*sqrt((sqrts+m1)^2-σ1)/(2*sqrt(σ1));
        σ3 = m1sq + m2sq +
        # 2*(  E2 * E1 -
            2*( (σ1 + m2sq - m3sq)/(2*sqrt(σ1)) * (s-m1sq-σ1)/(2*sqrt(σ1)) -
            # |p2|*cos(theta23) * (-|p1|)
            sqrt(Kallen(σ1, m2sq, m3sq)/(4*σ1))*cosθ23 * (-sqrtλσ1) );
        # caluclate p3 in (23) frame
        p23bu = sqrt(Kallen(σ1, m2sq, m3sq)/(4*σ1));
        p3_in23 = [(σ1 + m3sq - m2sq)/sqrt(4*σ1),
                   -p23bu*sqrt(1-cosθ23^2)*cos(ϕ23),
                   -p23bu*sqrt(1-cosθ23^2)*sin(ϕ23),
                   -p23bu*cosθ23];
        # boost to lab frame from (23) frame
        γ1 = (s + σ1 - m1sq)/sqrt(4*s*σ1);
        β1 = sqrt(1.0-1.0/γ1^2);
        p3_b = [γ1*(p3_in23[1]+β1*p3_in23[4]),
                p3_in23[2],
                p3_in23[3],
                γ1*(β1*p3_in23[1]+p3_in23[4])];

        ct1 = cosθ1; st1 = sqrt(1.0-cosθ1^2); cp1 = cos(ϕ1); sp1 = sin(ϕ1);
        # Rz(phi1) * Ry(theta1) * p3_boost
        p3_rot = [p3_b[1],
                  cp1*ct1* p3_b[2] + (-sp1)* p3_b[3] + cp1*st1* p3_b[4],
                  sp1*ct1* p3_b[2] +   cp1 * p3_b[3] + sp1*st1* p3_b[4],
                     -st1* p3_b[2] +     0 * p3_b[3] +     ct1* p3_b[4]];

        cosθ3 = -p3_rot[4]/sqrt(p3_b[2]^2+p3_b[3]^2+p3_b[4]^2);
        ϕ3 = (p3_rot[2] != zero(p3_rot[2])) ? atan(-p3_rot[3], -p3_rot[2]) : rand()*one(p3_rot[2]);

        cosθ12_n = m2sq + m3sq + 2* (σ3+m2sq-m1sq)/(2*sqrt(σ3)) * (s-m3sq-σ3)/(2*sqrt(σ3)) - σ1;
        m3 = sqrt(m3sq)
        sqrtλσ3 = sqrt((sqrts-m3)^2-σ3)*sqrt((sqrts+m3)^2-σ3)
        cosθ12_d = 2 * sqrt(Kallen(σ3, m1sq, m2sq)/(4*(σ3))) * sqrtλσ3 / (2*sqrt(σ3));
        cosθ12 = cosθ12_d ≈ 0.0+0.0im ? 2.0*rand()-1.00*one(σ1) : cosθ12_n/cosθ12_d;

        n1 = [0.,
             -sqrt(1.0-cosθ1^2)*cos(ϕ1),
             -sqrt(1.0-cosθ1^2)*sin(ϕ1),
             -cosθ1];
        ct3 = cosθ3; st3 = sqrt(1.0-cosθ3^2); cp3 = cos(ϕ3); sp3 = sin(ϕ3);
        # Rz(phi23) * Ry(theta23) * p3_boost
        n1_rot = [n1[1],
                  cp3*ct3* n1[2] + sp3*ct3* n1[3] + (-st3)* n1[4],
                   (-sp3)* n1[2] +    cp3 * n1[3] +   0.0 * n1[4],
                  cp3*st3* n1[2] + sp3*st3* n1[3] +   ct3 * n1[4]];
        # println(n1_rot[3], " ", n1_rot[2])
        ϕ12 = (n1_rot[2] != zero(n1_rot[2])) ? atan(n1_rot[3], n1_rot[2]) : rand()*one(n1_rot[2]);
        return σ3,cosθ3,ϕ3,cosθ12,ϕ12
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay