text:
"""
    SixDOFAeroEuler
    SixDOFAeroEuler(x::AbstractVector)
    SixDOFAeroEuler(dss::DSState)

Six degrees of freedom dynamic system using TAS, α, β for velocity representation and Euler
angles for attitude.

Flat Earth hypothesis is applied and Earth reference frame is considered inertial.

It is considered that the aircraft xb-zb plane is a plane of symmetry so that Jxy and Jyz
cross-product of inertia are zero and will not be taken into account.

# Fields
- `x::SVector{13, T}`: state vector.
  - [tas (m/s), α (rad), β (rad), ϕ (rad), θ (rad), ψ (rad), p (rad/s), q (rad/s), r (rad/s),
    x (m), y (m), z (m), pow (%)].
"""
struct SixDOFAeroEuler{T}<:DSState{T}
    x::SVector{13, T}
end

SixDOFAeroEuler(x::AbstractVector) = SixDOFAeroEuler(SVector{13, eltype(x)}(x))

SixDOFAeroEuler(dss::DSState) = SixDOFAeroEuler(
    [
        get_tasαβ(dss)...,
        get_euler_angles(dss)[3:-1:1]...,
        get_ang_vel_body(dss)...,
        get_earth_position(dss)...,
        get_engine_power(dss),
    ]
)

get_x_names(dss::SixDOFAeroEuler) = [:tas, :α, :β, :ϕ, :θ, :ψ, :p, :q, :r, :x, :y, :z, :pow]

# Mandatory getters for DSState
get_earth_position(dss::SixDOFAeroEuler) = get_x(dss)[10:12]
get_height(dss::SixDOFAeroEuler) = -get_x(dss)[12]

get_euler_angles(dss::SixDOFAeroEuler) = get_x(dss)[6:-1:4]

get_tas(dss::SixDOFAeroEuler) = get_x(dss)[1]
get_α(dss::SixDOFAeroEuler) = get_x(dss)[2]
get_β(dss::SixDOFAeroEuler) = get_x(dss)[3]
get_tasαβ(dss::SixDOFAeroEuler) = get_x(dss)[1:3]
get_body_velocity(dss::SixDOFAeroEuler) = wind2body(get_tas(dss), 0, 0, get_α(dss), get_β(dss))
get_horizon_velocity(dss::SixDOFAeroEuler) = body2horizon(
    get_body_velocity(dss)..., get_euler_angles(dss)...
)

get_ang_vel_body(dss::SixDOFAeroEuler) = get_x(dss)[7:9]
get_euler_angles_rates(dss::SixDOFAeroEuler) = pqr_2_ψθϕ_dot(
    get_ang_vel_body(dss)..., get_euler_angles(dss)[2:3]...
)

get_engine_power(dss::SixDOFAeroEuler) = get_x(dss)[13]

# Mandatory getters for DSStateDot
get_tas_dot(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[1]
get_α_dot(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[2]
get_β_dot(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[3]
get_tasαβ_dot(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[1:3]
get_accel_body(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = tasαβ_dot_to_uvw_dot(
    get_tasαβ(dssd)..., get_tasαβ_dot(dssd)...
)

get_ang_accel_body(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[7:9]
get_engine_power_dot(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T} = get_xdot(dssd)[13]

function state_eqs(dss::SixDOFAeroEuler, time, mass, inertia, forces, moments, h, pow_dot)

    xdot = sixdof_aero_earth_euler_fixed_mass(
         time,
         get_x(dss),
         mass,
         inertia,
         forces,
         moments,
         h
    )

    xdot = [xdot..., pow_dot]

    return DSStateDot(dss, xdot)
end

"""
    sixdof_aero_earth_euler_fixed_mass!(x_dot, time, x, mass, inertia, forces, moments, h)

Mutating version of [six_dof_aero_euler_fixed_mass](@ref).
"""
function sixdof_aero_earth_euler_fixed_mass!(
    x_dot, time, x, mass, inertia, forces, moments, h,
    )
    # C     x(1)  -> tas (m/s)
    # C     x(2)  -> α (rad)
    # C     x(3)  -> β (rad)
    # C     x(4)  -> ϕ (rad)
    # C     x(5)  -> θ (rad)
    # C     x(6)  -> ψ (rad)
    # C     x(7)  -> p (rad/s)
    # C     x(8)  -> q (rad/s)
    # C     x(9)  -> r (rad/s)
    # C     x(10) -> North (m)
    # C     x(11) -> East (m)
    # C     x(12) -> Altitude (m)

    # Assign state
    tas = x[1]
    α = x[2]
    β = x[3]
    ϕ = x[4]
    θ = x[5]
    ψ = x[6]
    p = x[7]
    q = x[8]
    r = x[9]

    # Unpack forces
    Fx, Fy, Fz = forces
    # Unpack moments
    L, M, N = moments
    # Unpack angular momentum contributions
    hx, hy, hz = h
    # Unpack inertia
    Ixx, Iyy, Izz = inertia[1, 1], inertia[2, 2], inertia[3, 3]
    Ixz = inertia[1, 3]

    # Get ready for state equations
    u, v, w = wind2body(tas, 0, 0, α, β)

    sψ, cψ = sin(ψ), cos(ψ)
    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    # Force equations
    udot = r * v - q * w + Fx / mass
    vdot = p * w - r * u + Fy / mass
    wdot = q * u - p * v + Fz / mass

    x_dot[1], x_dot[2], x_dot[3] = uvw_dot_to_tasαβ_dot(u, v, w, udot, vdot, wdot)

    # Kinematics
    x_dot[6], x_dot[5], x_dot[4] = pqr_2_ψθϕ_dot(p, q, r, θ, ϕ)

    # Moments
    pq = p * q
    qr = q * r

    rhy_qhz = (r * hy - q * hz)
    qhx_phy = (q * hx - p * hy)

    # If inertia is constant this terms are constant too.
    # TODO: think about passing them as arguments to improve speed.
    IxzS = Ixz^2
    xpq = Ixz * (Ixx - Iyy + Izz)
    gam = Ixx * Izz - IxzS
    xqr = Izz * (Izz - Iyy) + IxzS
    zpq = (Ixx - Iyy) * Ixx + IxzS
    ypr = Izz - Ixx

    x_dot[7] = (xpq * pq - xqr * qr + Izz * (L + rhy_qhz) + Ixz * (N + qhx_phy)) / gam
    x_dot[8] = (ypr * p * r - Ixz * (p^2 - r^2) + M - r * hx + p * hz) / Iyy
    x_dot[9] = (zpq * pq - xpq * qr + Ixz * (L + rhy_qhz) + Ixx * (N + qhx_phy)) / gam

    # Navigation
    t1 = sϕ * cψ
    t2 = cϕ * sθ
    t3 = sϕ * sψ
    s1 = cθ * cψ
    s2 = cθ * sψ
    s3 = t1 * sθ - cϕ * sψ
    s4 = t3 * sθ + cϕ * cψ
    s5 = sϕ * cθ
    s6 = t2 * cψ + t3
    s7 = t2 * sψ - t1
    s8 = cϕ * cθ

    x_dot[10] = u * s1 + v * s3 + w * s6  # North speed
    x_dot[11] = u * s2 + v * s4 + w * s7  # East speed
    x_dot[12] = -u * sθ + v * s5 + w * s8  # Vertical speed (possitive downwards)

end

"""
    sixdof_aero_earth_euler_fixed_mass(time, x, mass, inertia, forces, moments, h)

State equations for `SixDOFAeroEuler` dynamic system.

# Arguments

- `time::Number`: time (s).
- `x::AbstractArray`: state vector.
  - [tas (m/s), α (rad), β (rad), ϕ (rad), θ (rad), ψ (rad), p (rad/s), q (rad/s), r (rad/s),
    x (m), y (m), z (m), pow (%)].
- `mass::Number`: mass (kg).
- `inertia::AbstractMatrix`: 3x3 inertia tensor (kg·m^2).
- `forces::AbstractVector`: body axis forces [Fx, Fy, Fz] (N).
- `moments::AbstractVector`: body axis moments (L, M, N) (N·m).
- `h::AbstractVector`:  Additional angular momentum contributions such as those coming from
 spinning rotors (kg·m²/s).

# Returns
- `x_dot::Array`: time derivative of the state vector.

# Notes

It is considered that the aircraft xb-zb plane is a plane of symmetry so that Jxy and Jyz
cross-product of inertia are zero and will not be taken into account.

# See also

[six_dof_aero_euler_fixed_mass!](@ref)
"""
function sixdof_aero_earth_euler_fixed_mass(time, x, mass, inertia, forces, moments, h)
    x_dot = Array{eltype(x)}(undef, 12)
    sixdof_aero_earth_euler_fixed_mass!(x_dot, time, x, mass, inertia, forces, moments, h)
    return x_dot
end

function get_trimmer_cost(dssd::DSStateDot{S, N, T}) where {S<:SixDOFAeroEuler, N, T}
    x_dot = get_xdot(dssd)
    return [x_dot[1:3]..., x_dot[7:9]...]
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT, BAT.Logging
using Compat.Test
using Distributions, PDMats, StatsBase

@testset "mcmc_rand" begin
    @testset "rand" begin
        set_log_level!(BAT, LOG_WARNING)

        mvec = [-0.3, 0.3]
        cmat = [1.0 1.5; 1.5 4.0]
        Σ = @inferred PDMat(cmat)
        mv_dist = MvNormal(mvec, Σ)
        density = @inferred MvDistDensity(mv_dist)
        bounds = @inferred HyperRectBounds([-5, -8], [5, 8], reflective_bounds)
        nsamples_per_chain = 2000
        nchains = 4

        algorithmMW = @inferred MetropolisHastings()
        @test BAT.mcmc_compatible(algorithmMW, GenericProposalDist(mv_dist), NoParamBounds(2))
        samples, sampleids, stats = @inferred rand(
            MCMCSpec(algorithmMW, density, bounds),
            nsamples_per_chain,
            nchains,
            max_time = Inf,
            granularity = 1
        )

        @test length(samples) == length(sampleids)
        @test length(samples) == nchains * nsamples_per_chain
        @test samples.params[:, findmax(samples.log_value)[2]] == stats.mode

        cov_samples = cov(samples.params, FrequencyWeights(samples.weight), 2; corrected=true)
        mean_samples = mean(Array(samples.params), FrequencyWeights(samples.weight), 2)

        @test isapprox(mean_samples, mvec; atol = 0.2)
        @test isapprox(cov_samples, cmat; atol = 0.5)

        algorithmPW = @inferred MetropolisHastings(MHAccRejProbWeights())
        samples, sampleids, stats = @inferred rand(
            MCMCSpec(algorithmPW, mv_dist, bounds),
            nsamples_per_chain,
            nchains,
            max_time = Inf,
            granularity = 1
        )

        @test length(samples) == length(sampleids)
        @test samples.params[:, findmax(samples.log_value)[2]] == stats.mode

        cov_samples = cov(samples.params, FrequencyWeights(samples.weight), 2; corrected=true)
        mean_samples = mean(Array(samples.params), FrequencyWeights(samples.weight), 2)

        @test isapprox(mean_samples, mvec; atol = 0.2)
        @test isapprox(cov_samples, cmat; atol = 0.5)

        algorithmFW = @inferred MetropolisHastings(MHPosteriorFractionWeights())

        samples, sampleids, stats = @inferred rand(
            MCMCSpec(algorithmFW, density, bounds),
            nsamples_per_chain,
            nchains,
            max_time = Inf,
            granularity = 1
        )

        @test length(samples) == length(sampleids)
        @test samples.params[:, findmax(samples.log_value)[2]] == stats.mode

        cov_samples = cov(samples.params, FrequencyWeights(samples.weight), 2; corrected=true)
        mean_samples = mean(Array(samples.params), FrequencyWeights(samples.weight), 2)

        @test isapprox(mean_samples, mvec; atol = 0.2)
        @test isapprox(cov_samples, cmat; atol = 0.5)

        algorithmDS = @inferred DirectSampling()
        @test BAT.mcmc_compatible(algorithmDS, GenericProposalDist(mv_dist), NoParamBounds(2))
        samples, sampleids, stats = @inferred rand(
            MCMCSpec(algorithmDS, density, bounds),
            nsamples_per_chain,
            nchains,
            max_time = Inf,
            granularity = 1
        )

        cov_samples = cov(samples.params, FrequencyWeights(samples.weight), 2; corrected=true)
        mean_samples = mean(Array(samples.params), FrequencyWeights(samples.weight), 2)

        @test isapprox(mean_samples, mvec; atol = 0.2)
        @test isapprox(cov_samples, cmat; atol = 0.5)
    end
end

ext: jl


text:
using Random
using Distributed

struct SerialIterator
    len::Int
    get_example::Function
    batchsize::Int
    indices::Vector
    function SerialIterator(dataset::Dataset, batchsize::Int; shuffle=true)
        indices=Vector(1:dataset.len)
        if shuffle
            shuffle!(indices)
        end
        _get_example = i -> get_example(dataset, i)
        new(length(dataset), _get_example, batchsize, indices)
    end
end

function Base.iterate(diter::SerialIterator, state=(1, 0))
    idx_start, count = state
    if idx_start + diter.batchsize > diter.len
        return nothing
    else
        indices = diter.indices[idx_start:idx_start+diter.batchsize-1]
        element = diter.get_example.(indices)
        return (element, (idx_start + diter.batchsize, count + 1))
    end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
load("flatten.jl")

code = quote
    A = B.*C + D[j,i]
    dest[2i, 2j] = A
    dest2[] = 2A
end

symbols = Dict{Symbol,SymEntry}()
flat_code = {}
receive = ex->append!(flat_code, {ex})
context = Context(symbols, receive)

value = flatten(context, code)

print_code(flat_code)
println("value = $value")
print_context(context)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
@doc raw"""
    translate_forces_moments(forces_1, moments_1, p1, p2)

Calculate equivalent moment at point 2 (p2) given the forces (forces_1)  and moments
(moments_1) at point 1 (p1).

`` r = p_1 - p_2 ``

``M_{2} = M_{1} + r \times F_{1}``
"""
function translate_forces_moments(forces_1, moments_1, p1, p2)
    r21 = p1 - p2
    moments_2 = moments_1 + r21 × forces_1
    return moments_2
end

@doc raw"""
    steiner_inertia(cg, inertia_g, mass, p2)

Calculate the inertia tensor of a rigid solid at point 2 (p2), given the inertia tensor
(inertia_g) at the center of gravity (cg) and the mass of the system.

`` r = p_2 - cg ``

``I_{2} = I_{cg} + m (r^T · r I - r · r^T) ``
"""
function steiner_inertia(cg, inertia1, mass, p2)
    r = p2 - cg
    inertia2 = inertia1 + mass * ((r' * r) * I - r*r')
    return inertia2
end

"""
    coordinated_turn_bank(ψ_dot, α, β, tas, γ[, g])

Calculate roll angle (ϕ) [rad] for a given turn rate, angle of attack, angle of sideslip,
tas and flight path angle in the absence of wind for a coordinated turn bank.

Imposes sum of forces along y body axis equal to zero.

# Arguments
- `ψ_dot`: turn rate (rad/s).
- `α`: angle of attack (rad).
- `β`: angle of sideslip (rad).
- `tas`: true air speed (m/s).
- `γ`: flight path angle (rad).
- `g`: gravity (m/s²). Optional. Default value is `gD`.
"""
function coordinated_turn_bank(ψ_dot, α, β, tas, γ, g=gD)
    G = ψ_dot * tas / g

    if abs(γ) < 1e-8
        ϕ = G * cos(β) / (cos(α) - G * sin(α) * sin(β))
        ϕ = atan(ϕ)
    else
        a = 1 - G * tan(α) * sin(β)
        b = sin(γ) / cos(β)
        c = 1 + G^2 * cos(β)^2

        sq = sqrt(c * (1 - b^2) + G^2 * sin(β)^2)

        num = (a - b^2) + b * tan(α) * sq
        den = a ^ 2 - b^2 * (1 + c * tan(α)^2)

        ϕ = atan(G * cos(β) / cos(α) * num / den)
    end
    return ϕ
end


text:
function simd_loop_example_from_manual(x, y, z)
    s = zero(eltype(z))
    n = min(length(x),length(y),length(z))
    @simd for i in 1:n
        @inbounds begin
            z[i] = x[i]-y[i]
            s += z[i]*z[i]
        end
    end
    s
end

function simd_loop_with_multiple_reductions(x, y, z)
    # Use non-zero initial value to make sure reduction values include it.
    (s,t) = (one(eltype(x)),one(eltype(y)))
    @simd for i in 1:length(z)
        @inbounds begin
            s += x[i]
            t += 2*y[i]
            s += z[i]   # Two reductions go into s
        end
    end
    (s,t)
end

for T in {Int32,Int64,Float32,Float64}
   # Try various lengths to make sure "remainder loop" works
   for n in {0,1,2,3,4,255,256,257}
        # Dataset chosen so that results will be exact with only 24 bits of mantissa
        a = convert(Array{T},[2*j+1 for j in 1:n])
        b = convert(Array{T},[3*j+2 for j in 1:n])
        c = convert(Array{T},[5*j+3 for j in 1:n])
        s = simd_loop_example_from_manual(a,b,c)

        @test a==[2*j+1 for j in 1:n]
        @test b==[3*j+2 for j in 1:n]
        @test c==[-j-1 for j in 1:n]
        @test s==sum(c.*c)
        (s,t) = simd_loop_with_multiple_reductions(a,b,c)
        @test s==sum(a)+sum(c)+1
        @test t==2*sum(b)+1
    end
end

# Test that scope rules match regular for
let j=4
    # Use existing local variable.
    @simd for j=1:0 end
    @test j==4
    @simd for j=1:3 end
    @test j==3

    # Use global variable
    global simd_glob = 4
    @simd for simd_glob=1:0 end
    @test simd_glob==4
    @simd for simd_glob=1:3 end
    @test simd_glob==3

    # Index that is local to loop
    @simd for simd_loop_local=1:0 end
    simd_loop_local_present = true
    try
        simd_loop_local += 1
    catch
        simd_loop_local_present = false
    end
    @test !simd_loop_local_present
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function istriustrict(A::Matrix)
    m, n = size(A)
    for j = 1:min(n,m-1), i = j:m
        if A[i,j] != 0.
            return false
        end
    end
    return true
end

function istrilstrict(A::Matrix)
    m, n = size(A)
    for j = 2:n, i = 1:min(j,m)
        if A[i,j] != 0.
            return false
        end
    end
    return true
end

function simd_scale!(x, a)
    @simd for i=1:length(x)
        @inbounds x[i] *= a
    end
    nothing
end

function simd_copy!(x, y)
    @assert length(x) == length(y)
    @simd for i=1:length(y)
        @inbounds y[i] = x[i]
    end
    nothing
end

function simd_copy_scale!(a, x, y)
    @assert length(x) == length(y)
    @simd for i=1:length(y)
        @inbounds y[i] = a*x[i]
    end
    nothing
end

"Copy the first dimension of a 2D array y into a 1D array x."
function simd_copy_xy_first!(x, y, j)
    @assert length(x) == size(y, 1)
    @simd for i=1:length(x)
        @inbounds x[i] = y[i,j]
    end
    nothing
end

"Copy the first dimension of a 3D array y into a 1D array x."
function simd_copy_xy_first!(x, y, j, k)
    @assert length(x) == size(y, 1)
    @simd for i=1:length(x)
        @inbounds x[i] = y[i,j,k]
    end
    nothing
end

"Copy a 1D array x into the first dimension of a 2D array y."
function simd_copy_yx_first!(x, y, j)
    @assert length(x) == size(y, 1)
    @simd for i=1:size(y, 1)
        @inbounds y[i,j] = x[i]
    end
    nothing
end

"Copy a 1D array x into the first dimension of a 3D array y."
function simd_copy_yx_first!(x, y, j, k)
    @assert length(x) == size(y, 1)
    @simd for i=1:size(y, 1)
        @inbounds y[i,j,k] = x[i]
    end
    nothing
end

function simd_copy_yx_first_last!(x, y, j)
    @assert size(x, 1) == size(y, 1)
    @assert size(x, 2) == size(y, 3)
    @simd for k=1:size(y, 3)
        @simd for i=1:size(y, 1)
            @inbounds y[i,j,k] = x[i,k]
        end
    end
    nothing
end

function simd_xpy!(x, y)
    @assert length(x) == length(y)
    @simd for i=1:length(y)
        @inbounds y[i] += x[i]
    end
    nothing
end

function simd_axpy!(a, x, y)
    @assert length(x) == length(y)
    @simd for i=1:length(y)
        @inbounds y[i] += a*x[i]
    end
    nothing
end

function simd_wxpy!(w, x, y)
    @assert length(x) == length(y) == length(w)
    @simd for i=1:length(w)
        @inbounds w[i] = x[i] + y[i]
    end
    nothing
end

function simd_waxpy!(w, a, x, y)
    @assert length(x) == length(y) == length(w)
    simd_copy!(y, w)
    simd_axpy!(a, x, w)
    nothing
end

function simd_mult!(w, X, y)
    @assert length(w) == size(X, 1)
    @assert length(y) == size(X, 2)
    for i=1:length(w)
        w[i] = 0.
        @simd for j=1:length(y)
            @inbounds w[i] += X[i,j] * y[j]
        end
    end
    nothing
end

function simd_mult!(w, a, X, y)
    @assert length(w) == size(X, 1)
    @assert length(y) == size(X, 2)
    @simd for i=1:length(w)
        w[i] = 0
        @simd for j=1:length(y)
            @inbounds w[i] += a * X[i,j] * y[j]
        end
    end
    nothing
end


text:
# Bipartite matching, using the Ford–Fulkerson algorithm (with DFS). (Algorithms
# with better asymptotic running time exist, but they may not be more efficient
# in practice, and they tend to be quite a bit more complex.) `X` should be a
# matrix where `X[i, g] ≠ 0` indicates an edge between `i` and `j`. I.e., this
# may be a boolean matrix, or any other numeric matrix where a matching along
# nonzero values (e.g., using the valuation matrix directly in the case of MNW)
# is desired. The function returns an iterator over pairs `(i => g)`,
# indicating that `i` is matched with `g`.
function bipartite_matching(X)

    n, m = size(X)

    mate, seen = zeros(Int, m), falses(m)

    function match!(i)
        for g = 1:m
            (iszero(X[i, g]) || seen[g]) && continue
            seen[g] = true
            if mate[g] == 0 || match!(mate[g])
                mate[g] = i
                return true
            end
        end
        return false
    end

    for i = 1:n
        seen .= false
        match!(i)
    end

    return ((i => g) for (g, i) in pairs(mate) if i ≠ 0)

end

bipartite_matching(X::Profile) = bipartite_matching(matrix(X))


text:
using Test
using DomainGistry

@testset "Domain Name Generation Fxn" begin
    @test DomainGistry.generate("example") == ["example.com" ,"example.edu" ,"example.net" ,"example.org" ,"example.site","example.co"  ,"example.io"  ,"example.ai"  ,"example.app" ,"example.ca"  ,"example.uk"  ,"example.ua"  ,"example.us"  ,"example.ru"  ,"example.ch"]
    @test DomainGistry.generate("example","common") == ["example.com" ,"example.edu" ,"example.net" ,"example.org" ,"example.site","example.co"  ,"example.io"  ,"example.ai"  ,"example.app" ,"example.ca"  ,"example.uk"  ,"example.ua"  ,"example.us"  ,"example.ru"  ,"example.ch"]
end


text:
export S²
export ComplexLine
export Cartesian
export Spherical
export Geographic

"""
    Represents a point in 2-sphere.
"""
abstract type S² end

"""
    Represents a point in the complex line.

field: z.
"""
struct ComplexLine <: S²
    z::Complex
end

"""
    Represents a point in the Cartesian coordinate system.

fields: x, y and z.
"""
struct Cartesian <: S²
    r::ℝ³
    Cartesian(r::ℝ³) = new(r)
    Cartesian(x::Real, y::Real, z::Real) = Cartesian(ℝ³(x, y, z))
    Cartesian(a::Array) = begin
        @assert(length(a) == 3, "The length of the input vector must be exactly 3.")
        Cartesian(a...)
    end
end

"""
    Represents a point in the spherical coordinate system.

fields: r, ϕ and θ.
"""
struct Spherical <: S²
    r::Float64
    ϕ::Float64
    θ::Float64
    Spherical(r::Real, ϕ::Real, θ::Real) = begin
        @assert(r ≥ 0, "r must be non-negative.")
        @assert(-π ≤ ϕ ≤ π, "ϕ must be in [-π, π].")
        @assert(0 ≤ θ ≤ π, "θ must be in [0, π].")
        new(float(r), float(ϕ), float(θ))
    end
    Spherical(a::Array) = begin
        @assert(length(a) == 3, "The length of the input vector must be exactly 3.")
        Spherical(a...)
    end
end

"""
    Represents a point in the geographic coordinate system.

fields: r, ϕ and θ.
"""
struct Geographic <: S²
    r::Float64
    ϕ::Float64
    θ::Float64
    Geographic(r::Real, ϕ::Real, θ::Real) = begin
        @assert(r ≥ 0, "r must be non-negative: $r")
        #@assert(-π ≤ ϕ ≤ π, "ϕ must be in [-π, π]: $ϕ")
        #@assert(-π/2 ≤ θ ≤ π/2, "θ must be in [-π/2, π/2]: $θ")
        new(float(r), float(ϕ), float(θ))
    end
    Geographic(a::Array) = begin
        @assert(length(a) == 2, "The length of the input vector must be exactly 2.")
        Geographic(a...)
    end
end

Base.real(cl::ComplexLine) = real(cl.z)
Base.imag(cl::ComplexLine) = imag(cl.z)
Base.abs(cl::ComplexLine) = abs(cl.z)
Base.angle(cl::ComplexLine) = angle(cl.z)
Base.conj(cl::ComplexLine) = ComplexLine(conj(cl.z))

Base.vec(cl::ComplexLine) = [cl.z; conj(cl).z]
Base.vec(c::Cartesian) = vec(c.r)
Base.vec(s::Spherical) = [s.r; s.ϕ; s.θ]
Base.vec(g::Geographic) = [g.r; g.ϕ; g.θ]

Base.isapprox(c1::Cartesian, c2::Cartesian) = isapprox(vec(c1), vec(c2))
Base.isapprox(s1::S², s2::S²) = isapprox(Cartesian(s1), Cartesian(s2))
Base.isapprox(a::Array{<:S²,N} where N,
              b::Array{<:S²,N} where N;
              atol::Float64 = TOLERANCE) = begin
    for (elementa, elementb) in zip(a, b)
        if isapprox(vec(elementa), vec(elementb), atol = atol) == false
            return false
        end
    end
    return true
end

ℝ³(c::Cartesian) = c.r
ℝ³(s::S²) = Cartesian(s).r

norm(c::Cartesian) = norm(c.r)

ComplexLine(c::Cartesian) = begin
    x, y, z = vec(c)
    ComplexLine(x / (1 - z) + im * y / (1 - z))
end
ComplexLine(r::ℝ³) = ComplexLine(Cartesian(r))
ComplexLine(s::Spherical) = ComplexLine(Cartesian(s))
ComplexLine(g::Geographic) = ComplexLine(Cartesian(g))

Cartesian(cl::ComplexLine) = begin
    x, y = real(cl), imag(cl)
    d = x^2 + y^2 + 1
    Cartesian(2x / d, 2y / d, (d - 2) / d)
end
Cartesian(s::Spherical) = begin
    r, ϕ, θ = vec(s)
    x = r * sin(θ) * cos(ϕ)
    y = r * sin(θ) * sin(ϕ)
    z = r * cos(θ)
    Cartesian(x, y, z)
end
Cartesian(g::Geographic) = begin
    r, ϕ, θ = vec(g)
    x = r * cos(θ) * cos(ϕ)
    y = r * cos(θ) * sin(ϕ)
    z = r * sin(θ)
    Cartesian(x, y, z)
end

Spherical(c::Cartesian) = begin
    x, y, z = vec(c)
    r = norm(c)
    Spherical(r, atan(y, x), acos(z / r))
end
Spherical(cl::ComplexLine) = Spherical(Cartesian(cl))
Spherical(g::Geographic) = Spherical(g.r, g.ϕ, pi / 2 - g.θ)
Spherical(r::ℝ³) = Spherical(Cartesian(r))

Geographic(g::Geographic) = g
Geographic(s::Spherical) = Geographic(s.r, s.ϕ, pi / 2 - s.θ)
Geographic(c::Cartesian) = begin
    x, y, z = vec(c)
    r = norm(c)
    Geographic(r, atan(y, x), asin(z / r))
end
Geographic(r::ComplexLine) = Geographic(Cartesian(r))
Geographic(r::ℝ³) = Geographic(Cartesian(r))

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module RootHitFiles

import Base, Tables

using Mmap: mmap
using Parsers

export RootHitFile

"""
    RootHitFile(file::Union{IOStream, AbstractString})

represents a `.root.hits` file, which can be iterated or read to yield events, which are
tables of hits.
"""
struct RootHitFile
    stream::IOBuffer
end

RootHitFile(stream::IOStream) = RootHitFile(IOBuffer(mmap(stream)))
RootHitFile(path::AbstractString) = occursin(r".root.hits$", path) ? RootHitFile((open(path))) :
    throw(ArgumentError("$path is not a .root.hits file"))

EventTuple = NamedTuple{
    (:eventnum, :primcount, :pos, :E, :time, :particleID, :trkID, :trkparentID, :volumeID),
    Tuple{
        Int32, Int32, Vector{NTuple{3, Float32}}, Vector{Float32}, Vector{Float32},
        Vector{Int32}, Vector{Int32}, Vector{Int32}, Vector{String}
    }
}

function Base.read(f::RootHitFile)
    eventnum  = Parsers.parse(Int32, f.stream)
    hitcount  = Parsers.parse(Int32, f.stream)
    primcount = Parsers.parse(Int32, f.stream)

    # skip newline
    skip(f.stream, 1)

    pos         = Vector{NTuple{3, Float32}}(undef, hitcount)
    E           = Vector{            Float32}(undef, hitcount)
    time        = Vector{            Float32}(undef, hitcount)
    particleID  = Vector{              Int32}(undef, hitcount)
    trkID       = Vector{              Int32}(undef, hitcount)
    trkparentID = Vector{              Int32}(undef, hitcount)
    volumeID    = Vector{             String}(undef, hitcount)

    @inbounds for i in 1:hitcount
        pos[i] = (
            Parsers.parse(Float32, f.stream),
            Parsers.parse(Float32, f.stream),
            Parsers.parse(Float32, f.stream)
        )
        E[i]           = Parsers.parse(Float32, f.stream)
        time[i]        = Parsers.parse(Float32, f.stream)
        particleID[i]  = Parsers.parse(Int32, f.stream)
        trkID[i]       = Parsers.parse(Int32, f.stream)
        trkparentID[i] = Parsers.parse(Int32, f.stream)
        volumeID[i]    = String(readuntil(f.stream, UInt8('\n')))
    end

    return EventTuple((
        eventnum, primcount, pos, E, time,
        particleID, trkID, trkparentID, volumeID
    ))
end

Base.eof(f::RootHitFile) = eof(f.stream)

function Base.iterate(f::RootHitFile, state = nothing)
    eof(f) && return nothing
    return read(f), nothing
end

Base.IteratorSize(::Type{RootHitFile}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{RootHitFile}) = Base.HasEltype()
Base.eltype(::RootHitFile) = EventTuple

Tables.isrowtable(::Type{RootHitFile}) = true
Tables.schema(f::RootHitFile) = Tables.Schema(eltype(f))

end # RootHitFiles


text:
export prepare, getedge

using PyCall

# https://github.com/dhoegh/Hawk.jl/blob/master/src/Hawk.jl

filename = abspath(joinpath(dirname(@__FILE__), "hed.py"))

@pyimport imp
(path, name) = dirname(filename), basename(filename)
(name, ext) = rsplit(name, '.')

(file, filename, data) = imp.find_module(name, [path])
conv = imp.load_module(name, file, filename, data)

struct EdgeModel
    predictor
end

function prepare(model_path=Pkg.dir("HED","src","HED_pretrained_bsds.npy"); gpu=false)::EdgeModel
    println("model_path ", model_path)
    predictor =conv[:prepare](model_path)
    EdgeModel(predictor)
end

function getedge(m::EdgeModel, image_path, output=nothing)
    conv[:getedge](m.predictor, image_path, output)
end


text:
using Plots, Test

@testset "Showaxis" begin
    for value in Plots._allShowaxisArgs
        @test plot(1:5, showaxis = value)[1][:yaxis][:showaxis] isa Bool
    end
    @test plot(1:5, showaxis = :y)[1][:yaxis][:showaxis] == true
    @test plot(1:5, showaxis = :y)[1][:xaxis][:showaxis] == false
end

@testset "Magic axis" begin
    @test plot(1, axis = nothing)[1][:xaxis][:ticks] == []
    @test plot(1, axis = nothing)[1][:yaxis][:ticks] == []
end # testset

@testset "Categorical ticks" begin
    p1 = plot('A':'M', 1:13)
    p2 = plot('A':'Z', 1:26)
    p3 = plot('A':'Z', 1:26, ticks = :all)
    @test Plots.get_ticks(p1[1], p1[1][:xaxis])[2] == string.('A':'M')
    @test Plots.get_ticks(p2[1], p2[1][:xaxis])[2] == string.('C':3:'Z')
    @test Plots.get_ticks(p3[1], p3[1][:xaxis])[2] == string.('A':'Z')
end

@testset "Ticks getter functions" begin
    ticks1 = ([1, 2, 3], ("a", "b", "c"))
    ticks2 = ([4, 5], ("e", "f"))
    p1 = plot(1:5, 1:5, 1:5, xticks = ticks1, yticks = ticks1, zticks = ticks1)
    p2 = plot(1:5, 1:5, 1:5, xticks = ticks2, yticks = ticks2, zticks = ticks2)
    p = plot(p1, p2)
    @test xticks(p) == yticks(p) == zticks(p) == [ticks1, ticks2]
    @test xticks(p[1]) == yticks(p[1]) == zticks(p[1]) == ticks1
end

@testset "Axis limits" begin
    pl = plot(1:5, xlims = :symmetric, widen = false)
    @test Plots.xlims(pl) == (-5, 5)

    pl = plot(1:3)
    @test Plots.xlims(pl) == Plots.widen(1, 3)

    pl = plot([1.05, 2.0, 2.95], ylims = :round)
    @test Plots.ylims(pl) == (1, 3)

    pl = plot(1:3, xlims = (1, 5))
    @test Plots.xlims(pl) == (1, 5)

    pl = plot(1:3, xlims = (1, 5), widen = true)
    @test Plots.xlims(pl) == Plots.widen(1, 5)
end

@testset "3D Axis" begin
    ql = quiver([1, 2], [2, 1], [3, 4], quiver = ([1, -1], [0, 0], [1, -0.5]), arrow = true)
    @test ql[1][:projection] == "3d"
end

@testset "twinx" begin
    pl = plot(1:10, margin = 2Plots.cm)
    twpl = twinx(pl)
    pl! = plot!(twinx(), -(1:10))
    @test twpl[:right_margin] == 2Plots.cm
    @test twpl[:left_margin] == 2Plots.cm
    @test twpl[:top_margin] == 2Plots.cm
    @test twpl[:bottom_margin] == 2Plots.cm
end

@testset "aliases" begin
    p = plot(1:2, xl = "x label")
    @test p[1][:xaxis][:guide] === "x label"
    p = plot(1:2, xrange = (0, 3))
    @test xlims(p) === (0, 3)
    p = plot(1:2, xtick = [1.25, 1.5, 1.75])
    @test p[1][:xaxis][:ticks] == [1.25, 1.5, 1.75]
    p = plot(1:2, xlabelfontsize = 4)
    @test p[1][:xaxis][:guidefontsize] == 4
    p = plot(1:2, xgα = .07)
    @test p[1][:xaxis][:gridalpha] ≈ .07
    p = plot(1:2, xgridls = :dashdot)
    @test p[1][:xaxis][:gridstyle] === :dashdot
    p = plot(1:2, xgridcolor = :red)
    @test p[1][:xaxis][:foreground_color_grid] === RGBA{Float64}(1.,0.,0.,1.)
    p = plot(1:2, xminorgridcolor = :red)
    @test p[1][:xaxis][:foreground_color_minor_grid] === RGBA{Float64}(1.,0.,0.,1.)
    p = plot(1:2, xgrid_lw = .01)
    @test p[1][:xaxis][:gridlinewidth] ≈ .01
    p = plot(1:2, xminorgrid_lw = .01)
    @test p[1][:xaxis][:minorgridlinewidth] ≈ .01
    p = plot(1:2, xtickor = :out)
    @test p[1][:xaxis][:tick_direction] === :out
end


text:
"""
`KernelProduct(kernels::Array{Kernel})`
Create a multiplication of kernels.
One can also use the operator `*`
```
k1 = SqExponentialKernel()
k2 = LinearKernel()
k = KernelProduct([k1,k2])
kernelmatrix(k,X) == kernelmatrix(k1,X).*kernelmatrix(k2,X)
kernelmatrix(k,X) == kernelmatrix(k1*k2,X)
```
"""
struct KernelProduct{T,Tr} <: Kernel{T,Tr}
    kernels::Vector{Kernel}
end

function KernelProduct(kernels::AbstractVector{<:Kernel})
    KernelProduct{eltype(kernels),Transform}(kernels)
end

params(k::KernelProduct) = params.(k.kernels)
opt_params(k::KernelProduct) = opt_params.(k.kernels)
duplicate(k::KernelProduct,θ) = KernelProduct(duplicate.(k.kernels,θ))

Base.:*(k1::Kernel,k2::Kernel) = KernelProduct([k1,k2])
Base.:*(k1::KernelProduct,k2::KernelProduct) = KernelProduct(vcat(k1.kernels,k2.kernels)) #TODO Add test
Base.:*(k::Kernel,kp::KernelProduct) = KernelProduct(vcat(k,kp.kernels))
Base.:*(kp::KernelProduct,k::Kernel) = KernelProduct(vcat(kp.kernels,k))

Base.length(k::KernelProduct) = length(k.kernels)
metric(k::KernelProduct) = metric.(k.kernels) #TODO Add test
transform(k::KernelProduct) = transform.(k.kernels) #TODO Add test
transform(k::KernelProduct,x::AbstractVecOrMat) = transform.(k.kernels,[x]) #TODO Add test
transform(k::KernelProduct,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim) #TODO Add test

hadamard(x,y) = x.*y

function kernelmatrix(
    κ::KernelProduct,
    X::AbstractMatrix;
    obsdim::Int=defaultobs)
    reduce(hadamard,kernelmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
end

function kernelmatrix(
    κ::KernelProduct,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int=defaultobs)
    reduce(hadamard,_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
end

function kerneldiagmatrix(
    κ::KernelProduct,
    X::AbstractMatrix;
    obsdim::Int=defaultobs) #TODO Add test
    reduce(hadamard,kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
end


text:
"""
Create a geometry object of the appropriate type from it's well known
binary (WKB) representation.

### Parameters
* **pabyData**: pointer to the input BLOB data.
* **hSRS**: handle to the spatial reference to be assigned to the created
    geometry object. This may be `NULL`.
* **phGeometry**: the newly created geometry object will be assigned to
    the indicated handle on return. This will be NULL in case of failure.
    If not NULL, *phGeometry should be freed with OGR_G_DestroyGeometry()
    after use.
* **nBytes**: the number of bytes of data available in pabyData, or -1 if
    it is not known, but assumed to be sufficient.
### Returns
OGRERR_NONE if all goes well, otherwise any of OGRERR_NOT_ENOUGH_DATA,
OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or OGRERR_CORRUPT_DATA may be returned.
"""
function fromWKB(data)
    geom = Ptr{GDAL.OGRGeometryH}(0)
    result = GDAL.createfromwkb(data, C_NULL, geom, -1)
    if result != GDAL.OGRERR_NONE
        error("Failed to create geometry from WKB")
    end
    geom
end

function fromWKB(data, spatialref::SpatialRef)
    geom = Ptr{GDAL.OGRGeometryH}(0)
    result = GDAL.createfromwkb(data, spatialref, geom, -1)
    if result != GDAL.OGRERR_NONE
        error("Failed to create geometry from WKB")
    end
    geom
end

"""
Create a geometry object of the appropriate type from it's well known text
(WKT) representation.

### Parameters
* **ppszData**: input zero terminated string containing WKT representation
    of the geometry to be created. The pointer is updated to point just beyond
    that last character consumed.
* **hSRS**: handle to the spatial reference to be assigned to the created
    geometry object. This may be `NULL`.
* **phGeometry**: the newly created geometry object will be assigned to the
    indicated handle on return. This will be NULL if the method fails. If not
    `NULL`, *phGeometry should be freed with OGR_G_DestroyGeometry() after use.

### Returns
OGRERR_NONE if all goes well, otherwise any of OGRERR_NOT_ENOUGH_DATA,
OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or OGRERR_CORRUPT_DATA may be returned.
"""
function fromWKT(data)
    geom = Ptr{GDAL.OGRGeometryH}(0)
    result = GDAL.createfromwkt(data, C_NULL, geom)
    if result != GDAL.OGRERR_NONE
        error("Failed to create geometry from WKT")
    end
    geom
end

function fromWKT(data, spatialref::SpatialRef)
    geom = Ptr{GDAL.OGRGeometryH}(0)
    result = GDAL.createfromwkt(data, spatialref, geom)
    if result != GDAL.OGRERR_NONE
        error("Failed to create geometry from WKT")
    end
    geom
end

#"""
#     OGR_G_CreateFromFgf(unsigned char * pabyData,
#                         OGRSpatialReferenceH hSRS,
#                         OGRGeometryH * phGeometry,
#                         int nBytes,
#                         int * pnBytesConsumed) -> OGRErr
# """
# function createfromfgf(arg1,arg2::Ptr{OGRSpatialReferenceH},
#     arg3,arg4::Integer,arg5)
#     ccall((:OGR_G_CreateFromFgf,libgdal),OGRErr,(Ptr{Cuchar},
# Ptr{OGRSpatialReferenceH},Ptr{OGRGeometryH},Cint,Ptr{Cint}),
# arg1,arg2,arg3,arg4,arg5)
# end

"Destroy geometry object."
destroy(geom::Geometry) = GDAL.destroygeometry(geom)

"Create an empty geometry of desired type."
creategeom(geomtype::GDAL.OGRwkbGeometryType) = GDAL.creategeometry(geomtype)

function creategeom(f::Function, args...)
    geom = creategeom(args...)
    try
        f(geom)
    finally
        destroy(geom)
    end
end

# """
#     OGR_G_ApproximateArcAngles(double dfCenterX,
#                                double dfCenterY,
#                                double dfZ,
#                                double dfPrimaryRadius,
#                                double dfSecondaryRadius,
#                                double dfRotation,
#                                double dfStartAngle,
#                                double dfEndAngle,
#                                double dfMaxAngleStepSizeDegrees) -> OGRGeometryH
# Stroke arc to linestring.
# ### Parameters
# * **dfCenterX**: center X
# * **dfCenterY**: center Y
# * **dfZ**: center Z
# * **dfPrimaryRadius**: X radius of ellipse.
# * **dfSecondaryRadius**: Y radius of ellipse.
# * **dfRotation**: rotation of the ellipse clockwise.
# * **dfStartAngle**: angle to first point on arc (clockwise of X-positive)
# * **dfEndAngle**: angle to last point on arc (clockwise of X-positive)
# * **dfMaxAngleStepSizeDegrees**: the largest step in degrees along the arc, zero to use the default setting.
# ### Returns
# OGRLineString geometry representing an approximation of the arc.
# """
# function approximatearcangles(dfCenterX::Real,dfCenterY::Real,dfZ::Real,dfPrimaryRadius::Real,dfSecondaryAxis::Real,dfRotation::Real,dfStartAngle::Real,dfEndAngle::Real,dfMaxAngleStepSizeDegrees::Real)
#     checknull(ccall((:OGR_G_ApproximateArcAngles,libgdal),Ptr{OGRGeometryH},(Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble),dfCenterX,dfCenterY,dfZ,dfPrimaryRadius,dfSecondaryAxis,dfRotation,dfStartAngle,dfEndAngle,dfMaxAngleStepSizeDegrees))
# end

# """
#     OGR_G_ForceToPolygon(OGRGeometryH hGeom) -> OGRGeometryH
# Convert to polygon.
# ### Parameters
# * **hGeom**: handle to the geometry to convert (ownership surrendered).
# ### Returns
# the converted geometry (ownership to caller).
# """
# function forcetopolygon(arg1::Ptr{OGRGeometryH})
#     checknull(ccall((:OGR_G_ForceToPolygon,libgdal),Ptr{OGRGeometryH},(Ptr{OGRGeometryH},),arg1))
# end

# """
#     OGR_G_ForceToLineString(OGRGeometryH hGeom) -> OGRGeometryH
# Convert to line string.
# ### Parameters
# * **hGeom**: handle to the geometry to convert (ownership surrendered).
# ### Returns
# the converted geometry (ownership to caller).
# """
# function forcetolinestring(arg1::Ptr{OGRGeometryH})
#     checknull(ccall((:OGR_G_ForceToLineString,libgdal),Ptr{OGRGeometryH},(Ptr{OGRGeometryH},),arg1))
# end

# """
#     OGR_G_ForceToMultiPolygon(OGRGeometryH hGeom) -> OGRGeometryH
# Convert to multipolygon.
# ### Parameters
# * **hGeom**: handle to the geometry to convert (ownership surrendered).
# ### Returns
# the converted geometry (ownership to caller).
# """
# function forcetomultipolygon(arg1::Ptr{OGRGeometryH})
#     checknull(ccall((:OGR_G_ForceToMultiPolygon,libgdal),Ptr{OGRGeometryH},(Ptr{OGRGeometryH},),arg1))
# end

# """
#     OGR_G_ForceToMultiPoint(OGRGeometryH hGeom) -> OGRGeometryH
# Convert to multipoint.
# ### Parameters
# * **hGeom**: handle to the geometry to convert (ownership surrendered).
# ### Returns
# the converted geometry (ownership to caller).
# """
# function forcetomultipoint(arg1::Ptr{OGRGeometryH})
#     checknull(ccall((:OGR_G_ForceToMultiPoint,libgdal),Ptr{OGRGeometryH},(Ptr{OGRGeometryH},),arg1))
# end

# """
#     OGR_G_ForceToMultiLineString(OGRGeometryH hGeom) -> OGRGeometryH
# Convert to multilinestring.
# ### Parameters
# * **hGeom**: handle to the geometry to convert (ownership surrendered).
# ### Returns
# the converted geometry (ownership to caller).
# """
# function forcetomultilinestring(arg1::Ptr{OGRGeometryH})
#     checknull(ccall((:OGR_G_ForceToMultiLineString,libgdal),Ptr{OGRGeometryH},(Ptr{OGRGeometryH},),arg1))
# end

"""
Convert to another geometry type.

### Parameters
* **hGeom**: the input geometry - ownership is passed to the method.
* **eTargetType**: target output geometry type.
* **papszOptions**: options as a null-terminated list of strings or NULL.
"""
forceto(geom::Geometry, targettype::GDAL.OGRwkbGeometryType) =
    GDAL.forceto(geom, targettype, C_NULL)

"""
Get the dimension of this geometry.

### Returns
0 for points, 1 for lines and 2 for surfaces.
"""
getdim(geom::Geometry) = GDAL.getdim(geom)

"""
Get the dimension of the coordinates in this geometry.

### Returns
in practice this will return 2 or 3. It can also return 0 in the case of
an empty point.
"""
getcoorddim(geom::Geometry) = GDAL.getcoordinatedimension(geom)

"Set the coordinate dimension."
setcoorddim(geom::Geometry, dim::Integer) = GDAL.setcoordinatedimension(geom, dim)

"Make a copy of this object, with the spatial reference of the original."
clone(geom::Geometry) = GDAL.clone(geom)

# """
#     OGR_G_GetEnvelope(OGRGeometryH hGeom,
#                       OGREnvelope * psEnvelope) -> void
# Computes and returns the bounding envelope for this geometry in the passed psEnvelope structure.
# ### Parameters
# * **hGeom**: handle of the geometry to get envelope from.
# * **psEnvelope**: the structure in which to place the results.
# """
# function getenvelope(arg1::Ptr{OGRGeometryH},arg2)
#     ccall((:OGR_G_GetEnvelope,libgdal),Void,(Ptr{OGRGeometryH},Ptr{OGREnvelope}),arg1,arg2)
# end

# """
#     OGR_G_GetEnvelope3D(OGRGeometryH hGeom,
#                         OGREnvelope3D * psEnvelope) -> void
# Computes and returns the bounding envelope (3D) for this geometry in the passed psEnvelope structure.
# ### Parameters
# * **hGeom**: handle of the geometry to get envelope from.
# * **psEnvelope**: the structure in which to place the results.
# """
# function getenvelope3d(arg1::Ptr{OGRGeometryH},arg2)
#     ccall((:OGR_G_GetEnvelope3D,libgdal),Void,(Ptr{OGRGeometryH},Ptr{OGREnvelope3D}),arg1,arg2)
# end

"""
Assign geometry from well known binary data.

### Parameters
* **hGeom**: handle on the geometry to assign the well know binary data to.
* **pabyData**: the binary input data.
* **nSize**: the size of pabyData in bytes, or zero if not known.

### Returns
OGRERR_NONE if all goes well, otherwise any of OGRERR_NOT_ENOUGH_DATA,
OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or OGRERR_CORRUPT_DATA may be returned.
"""
function fromWKB(data, nsize::Integer = 0)
    geom = Ptr{GDAL.OGRGeometryH}()
    result = GDAL.importfromwkb(geom, data, nsize)
    if result != OGRERR_NONE
        error("Failed to import geometry from WKB")
    end
    geom
end

"""
Convert a geometry well known binary format.

### Parameters
* **hGeom**: handle on the geometry to convert to a well know binary data from.
* **eOrder**: One of wkbXDR or [wkbNDR] indicating MSB or LSB byte order resp.
* **pabyDstBuffer**: a buffer into which the binary representation is written.
    This buffer must be at least OGR_G_WkbSize() byte in size.
### Returns
Currently OGRERR_NONE is always returned.
"""
function toWKB(geom::Geometry, eorder::GDAL.OGRwkbByteOrder = GDAL.wkbNDR)
    buffer = Array(Cuchar, wkbsize(geom))
    GDAL.exporttowkb(geom, eorder, pointer(buffer))
    buffer
end

"""
    OGR_G_ExportToIsoWkb(OGRGeometryH hGeom,
                         OGRwkbByteOrder eOrder,
                         unsigned char * pabyDstBuffer) -> OGRErr
Convert a geometry into SFSQL 1.2 / ISO SQL/MM Part 3 well known binary format.
### Parameters
* **hGeom**: handle on the geometry to convert to a well know binary data from.
* **eOrder**: One of wkbXDR or [wkbNDR] indicating MSB or LSB byte order resp.
* **pabyDstBuffer**: a buffer into which the binary representation is written.
    This buffer must be at least OGR_G_WkbSize() byte in size.
### Returns
Currently OGRERR_NONE is always returned.
"""
function toISOWKB(geom::Geometry, eorder::GDAL.OGRwkbByteOrder = GDAL.wkbNDR)
    buffer = Array(Cuchar, wkbsize(geom))
    GDAL.exporttoisowkb(geom, eorder, pointer(buffer))
    buffer
end

"Returns size (in bytes) of related binary representation."
wkbsize(geom::Geometry) = GDAL.wkbsize(geom)

# """
#     OGR_G_ImportFromWkt(OGRGeometryH hGeom,
#                         char ** ppszSrcText) -> OGRErr
# Assign geometry from well known text data.
# ### Parameters
# * **hGeom**: handle on the geometry to assign well know text data to.
# * **ppszSrcText**: pointer to a pointer to the source text. The pointer is updated to pointer after the consumed text.
# ### Returns
# OGRERR_NONE if all goes well, otherwise any of OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or OGRERR_CORRUPT_DATA may be returned.
# """
# function importfromwkt(arg1::Ptr{OGRGeometryH},arg2)
#     ccall((:OGR_G_ImportFromWkt,libgdal),OGRErr,(Ptr{OGRGeometryH},Ptr{Cstring}),arg1,arg2)
# end

# """
#     OGR_G_ExportToWkt(OGRGeometryH hGeom,
#                       char ** ppszSrcText) -> OGRErr
# Convert a geometry into well known text format.
# ### Parameters
# * **hGeom**: handle on the geometry to convert to a text format from.
# * **ppszSrcText**: a text buffer is allocated by the program, and assigned to the passed pointer. After use, *ppszDstText should be freed with OGRFree().
# ### Returns
# Currently OGRERR_NONE is always returned.
# """
# function exporttowkt(arg1::Ptr{OGRGeometryH},arg2)
#     ccall((:OGR_G_ExportToWkt,libgdal),OGRErr,(Ptr{OGRGeometryH},Ptr{Cstring}),arg1,arg2)
# end

# """
#     OGR_G_ExportToIsoWkt(OGRGeometryH hGeom,
#                          char ** ppszSrcText) -> OGRErr
# Convert a geometry into SFSQL 1.2 / ISO SQL/MM Part 3 well known text format.
# ### Parameters
# * **hGeom**: handle on the geometry to convert to a text format from.
# * **ppszSrcText**: a text buffer is allocated by the program, and assigned to the passed pointer. After use, *ppszDstText should be freed with OGRFree().
# ### Returns
# Currently OGRERR_NONE is always returned.
# """
# function exporttoisowkt(arg1::Ptr{OGRGeometryH},arg2)
#     ccall((:OGR_G_ExportToIsoWkt,libgdal),OGRErr,(Ptr{OGRGeometryH},Ptr{Cstring}),arg1,arg2)
# end

"Fetch geometry type code"
getgeomtype(geom::Geometry) = GDAL.getgeometrytype(geom)

"Fetch WKT name for geometry type."
getgeomname(geom::Geometry) = GDAL.getgeometryname(geom)

# """
#     OGR_G_DumpReadable(OGRGeometryH hGeom,
#                        FILE * fp,
#                        const char * pszPrefix) -> void
# Dump geometry in well known text format to indicated output file.
# ### Parameters
# * **hGeom**: handle on the geometry to dump.
# * **fp**: the text file to write the geometry to.
# * **pszPrefix**: the prefix to put on each line of output.
# """
# function dumpreadable(arg1::Ptr{OGRGeometryH},arg2,arg3)
#     ccall((:OGR_G_DumpReadable,libgdal),Void,(Ptr{OGRGeometryH},Ptr{FILE},Cstring),arg1,arg2,arg3)
# end

"Convert geometry to strictly 2D."
flattento2d(geom::Geometry) = GDAL.flattento2d(geom)

"Force rings to be closed."
closerings(geom::Geometry) = GDAL.closerings(geom)

"Create geometry from GML."
fromGML(data) = GDAL.createfromgml(data)

"Convert a geometry into GML format."
toGML(geom::Geometry) = GDAL.exporttogml(geom)

# """
#     OGR_G_ExportToGMLEx(OGRGeometryH hGeometry,
#                         char ** papszOptions) -> char *
# Convert a geometry into GML format.
# ### Parameters
# * **hGeometry**: handle to the geometry.
# * **papszOptions**: NULL-terminated list of options.
# ### Returns
# A GML fragment or NULL in case of error.
# """
# function exporttogmlex(arg1::Ptr{OGRGeometryH},papszOptions)
#     bytestring(ccall((:OGR_G_ExportToGMLEx,libgdal),Cstring,(Ptr{OGRGeometryH},Ptr{Cstring}),arg1,papszOptions))
# end

# """
#     OGR_G_CreateFromGMLTree(const CPLXMLNode * psTree) -> OGRGeometryH
# """
# function createfromgmltree(arg1)
#     checknull(ccall((:OGR_G_CreateFromGMLTree,libgdal),Ptr{OGRGeometryH},(Ptr{CPLXMLNode},),arg1))
# end

# """
#     OGR_G_ExportToGMLTree(OGRGeometryH hGeometry) -> CPLXMLNode *
# """
# function exporttogmltree(arg1::Ptr{OGRGeometryH})
#     ccall((:OGR_G_ExportToGMLTree,libgdal),Ptr{CPLXMLNode},(Ptr{OGRGeometryH},),arg1)
# end

# """
#     OGR_G_ExportEnvelopeToGMLTree(OGRGeometryH hGeometry) -> CPLXMLNode *
# """
# function exportenvelopetogmltree(arg1::Ptr{OGRGeometryH})
#     ccall((:OGR_G_ExportEnvelopeToGMLTree,libgdal),Ptr{CPLXMLNode},(Ptr{OGRGeometryH},),arg1)
# end

"""
Convert a geometry into KML format.

### Parameters
* **geom**: the geometry to be converted.
* **altitudemode**: value to write in altitudeMode element, or NULL.
"""
toKML(geom::Geometry) = GDAL.exporttokml(geom, C_NULL)
toKML(geom::Geometry, altitudemode) = GDAL.exporttokml(geom, altitude)

"Convert a geometry into GeoJSON format."
toJSON(geom::Geometry) = GDAL.exporttojson(geom)

# """
#     OGR_G_ExportToJsonEx(OGRGeometryH hGeometry,
#                          char ** papszOptions) -> char *
# Convert a geometry into GeoJSON format.
# ### Parameters
# * **hGeometry**: handle to the geometry.
# * **papszOptions**: a null terminated list of options.
# ### Returns
# A GeoJSON fragment or NULL in case of error.
# """
# function exporttojsonex(arg1::Ptr{OGRGeometryH},papszOptions)
#     bytestring(ccall((:OGR_G_ExportToJsonEx,libgdal),Cstring,(Ptr{OGRGeometryH},Ptr{Cstring}),arg1,papszOptions))
# end

"""Create a geometry object of the appropriate type from its
GeoJSON representation"""
fromJSON(data) = GDAL.creategeometryfromjson(data)

"Assign spatial reference to this object."
setspatialref(geom::Geometry, spatialref::SpatialRef) =
    GDAL.assignspatialreference(geom, spatialref)

"Returns spatial reference system for geometry."
getspatialref(geom::Geometry) = GDAL.getspatialreference(geom)

"""
    OGR_G_Transform(OGRGeometryH hGeom,
                    OGRCoordinateTransformationH hTransform) -> OGRErr
Apply arbitrary coordinate transformation to geometry.

### Parameters
* **hGeom**: handle on the geometry to apply the transform to.
* **hTransform**: handle on the transformation to apply.

### Returns
OGRERR_NONE on success or an error code.
"""
function transform(geom::Geometry, coordtransform::CoordTransform)
    result = GDAL.transform(geom, coordtransform)
    if result != GDAL.OGRERR_NONE
        error("Failed to transform geometry")
    end
    geom
end

"Transform geometry to new spatial reference system."
function transformto(geom::Geometry, spatialref::SpatialRef)
    result = GDAL.transformto(geom, spatialref)
    if result != GDAL.OGRERR_NONE
        error("Failed to transform geometry to the new SRS")
    end
end

"""
Compute a simplified geometry.

### Parameters
* **geom**: the geometry.
* **tol**: the distance tolerance for the simplification.
"""
simplify(geom::Geometry, tol::Real) = GDAL.simplify(geom, tol)

"""
Simplify the geometry while preserving topology.

### Parameters
* **geom**: the geometry.
* **tol**: the distance tolerance for the simplification.
"""
simplifypreservetopology(geom::Geometry, tol::Real) =
    GDAL.simplifypreservetopology(geom, tol)

"""
Return a Delaunay triangulation of the vertices of the geometry.

### Parameters
* **geom**: the geometry.
* **tol**: optional snapping tolerance to use for improved robustness
* **onlyedges**: if TRUE, will return a MULTILINESTRING, otherwise it
    will return a GEOMETRYCOLLECTION containing triangular POLYGONs.
"""
delaunaytriangulation(geom::Geometry, tol::Real, onlyedges::Bool) =
    GDAL.delaunaytriangulation(geom, tol, onlyedges)

"""
Modify the geometry such it has no segment longer than the given distance.

### Parameters
* **geom**: the geometry to segmentize
* **maxLength**: the maximum distance between 2 points after segmentization
"""
segmentize(geom::Geometry, maxlength::Real) =
    GDAL.segmentize(geom, maxlength)

"Returns whether the features intersect"
intersects(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.intersects(geom1, geom2))

"Returns TRUE if the geometries are equivalent."
equals(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.equals(geom1, geom2))

"Returns TRUE if the geometries are disjoint."
disjoint(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.disjoint(geom1, geom2))

"Returns TRUE if the geometries are touching."
touches(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.touches(geom1, geom2))

"Returns TRUE if the geometries are crossing."
crosses(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.crosses(geom1, geom2))

"Returns TRUE if geom1 is contained within geom2."
within(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.within(geom1, geom2))

"Returns TRUE if geom1 contains geom2."
contains(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.contains(geom1, geom2))

"Returns TRUE if the geometries overlap."
overlaps(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.overlaps(geom1, geom2))

"Returns the boundary of the geometry."
boundary(geom1::Geometry, geom2::Geometry) = GDAL.boundary(geom1, geom2)

"Returns the convex hull of the geometry."
convexhull(geom::Geometry) = GDAL.convexhull(geom)

"""
Compute buffer of geometry.

### Parameters
* **geom**: the geometry.
* **dist**: the buffer distance to be applied. Should be expressed into the
    same unit as the coordinates of the geometry.
* **quadsegs**: the number of segments used to approximate a 90 degree
    (quadrant) of curvature.
"""
buffer(geom::Geometry, dist::Real, quadsegs::Integer=8) =
    GDAL.buffer(geom, dist, quadsegs)

"""
Compute the intersection of the geometries.

### Returns
a new geometry representing the union or NULL if an error occurs.
"""
intersection(geom1::Geometry, geom2::Geometry) =
    GDAL.intersection(geom1, geom2)

"""
Compute the union of the geometries.

### Returns
a new geometry representing the union or NULL if an error occurs.
"""
union(geom1::Geometry, geom2::Geometry) =
    GDAL.union(geom1, geom2)

"""
Compute the union of the geometry using cascading.

### Returns
a new geometry representing the union or NULL if an error occurs.
"""
unioncascaded(geom::Geometry) = GDAL.unioncascaded(geom)

"""
Returns a point guaranteed to lie on the surface or NULL if an error occurred.
"""
pointonsurface(geom::Geometry) = GDAL.pointonsurface(geom)

"""
Compute the difference of the geometries.

### Returns
a new geometry representing the difference or NULL if the difference is empty
or an error occurs.
"""
difference(geom1::Geometry, geom2::Geometry) =
    GDAL.difference(geom1, geom2)

"""
Compute the symmetric difference of the geometries.

### Returns
a new geometry representing the symmetric difference or NULL if the difference
is empty or an error occurs.
"""
symdifference(geom1::Geometry, geom2::Geometry) =
    GDAL.symdifference(geom1, geom2)

"""
Compute distance between two geometries.

### Returns
the distance between the geometries or -1 if an error occurs.
"""
distance(geom1::Geometry, geom2::Geometry) =
    GDAL.distance(geom1, geom2)

"""
Compute length of a geometry.

### Returns
the length or 0.0 for unsupported geometry types.
"""
geomlength(geom::Geometry) = GDAL.length(geom)

"""
Compute geometry area.

### Returns
the area or 0.0 for unsupported geometry types.
"""
geomarea(geom::Geometry) = GDAL.area(geom)

"Compute the geometry centroid."
function centroid!(geom::Geometry, centroid::Geometry)
    result = GDAL.centroid(geom, centroid)
    if result != GDAL.OGRERR_NONE
        error("Failed to compute the geometry centroid")
    end
end

"""
Fetch point at given distance along curve.

### Parameters
* **geom**: curve geometry.
* **distance**: distance along the curve at which to sample position.
    This distance should be between zero and geomlength() for this curve.
"""
pointalongline(geom::Geometry, distance::Real) =
    GDAL.value(geom, distance)

"Clear geometry information."
empty!(geom::Geometry) = GDAL.empty(geom)

"Returns TRUE if the geometry has no points, otherwise FALSE."
isempty(geom::Geometry) = Bool(GDAL.isempty(geom))

"Returns TRUE if the geometry is valid, otherwise FALSE."
isvalid(geom::Geometry) = Bool(GDAL.isvalid(geom))

"Returns TRUE if the geometry is simple, otherwise FALSE."
issimple(geom::Geometry) = Bool(GDAL.issimple(geom))

"Returns TRUE if the geometry is a ring, otherwise FALSE."
isring(geom::Geometry) = Bool(GDAL.isring(geom))

"""
Polygonizes a set of sparse edges.

### Returns
a handle to a newly allocated geometry now owned by the caller.
"""
polygonize(geom::Geometry) = GDAL.polygonize(geom)

equal(geom1::Geometry, geom2::Geometry) =
    Bool(GDAL.equal(geom1, geom2))

"Fetch number of points from a geometry."
npoint(geom::Geometry) = GDAL.getpointcount(geom)

# """
#     OGR_G_GetPoints(OGRGeometryH hGeom,
#                     void * pabyX,
#                     int nXStride,
#                     void * pabyY,
#                     int nYStride,
#                     void * pabyZ,
#                     int nZStride) -> int
# Returns all points of line string.
# ### Parameters
# * **hGeom**: handle to the geometry from which to get the coordinates.
# * **pabyX**: a buffer of at least (sizeof(double) * nXStride * nPointCount) bytes, may be NULL.
# * **nXStride**: the number of bytes between 2 elements of pabyX.
# * **pabyY**: a buffer of at least (sizeof(double) * nYStride * nPointCount) bytes, may be NULL.
# * **nYStride**: the number of bytes between 2 elements of pabyY.
# * **pabyZ**: a buffer of at last size (sizeof(double) * nZStride * nPointCount) bytes, may be NULL.
# * **nZStride**: the number of bytes between 2 elements of pabyZ.
# ### Returns
# the number of points
# """
# function getpoints(hGeom::Ptr{OGRGeometryH},pabyX,nXStride::Integer,pabyY,nYStride::Integer,pabyZ,nZStride::Integer)
#     ccall((:OGR_G_GetPoints,libgdal),Cint,(Ptr{OGRGeometryH},Ptr{Void},Cint,Ptr{Void},Cint,Ptr{Void},Cint),hGeom,pabyX,nXStride,pabyY,nYStride,pabyZ,nZStride)
# end

"Fetch the x coordinate of a point from a geometry, at index i."
getx(geom::Geometry, i::Integer) = GDAL.getx(geom, i)

"Fetch the y coordinate of a point from a geometry, at index i."
gety(geom::Geometry, i::Integer) = GDAL.gety(geom, i)

"Fetch the z coordinate of a point from a geometry, at index i."
getz(geom::Geometry, i::Integer) = GDAL.getz(geom, i)

"""
Fetch a point in line string or a point geometry, at index i.

### Parameters
* **i**: the vertex to fetch, from 0 to getNumPoints()-1, zero for a point.
"""
getpoint!(geom::Geometry, i::Integer, x, y, z) =
    GDAL.getpoint(geom, i, x, y, z)

function getpoint(geom::Geometry, i::Integer)
    x = Ref{Cdouble}()
    y = Ref{Cdouble}()
    z = Ref{Cdouble}()
    getpoint!(geom, i, x, y, z)
    (x[], y[], z[])
end

"""
Set number of points in a geometry.

### Parameters
* **geom**: the geometry.
* **n**: the new number of points for geometry.
"""
setpointcount(geom::Geometry, n::Integer) = GDAL.setpointcount(geom, n)

"""
Set the location of a vertex in a point or linestring geometry.

### Parameters
* **geom**: handle to the geometry to add a vertex to.
* **i**: the index of the vertex to assign (zero based) or zero for a point.
* **x**: input X coordinate to assign.
* **y**: input Y coordinate to assign.
* **z**: input Z coordinate to assign (defaults to zero).
"""
setpoint(geom::Geometry, i::Integer, x, y, z) =
    GDAL.setpoint(geom, i, x, y, z)

"""
Set the location of a vertex in a point or linestring geometry.

### Parameters
* **geom**: handle to the geometry to add a vertex to.
* **i**: the index of the vertex to assign (zero based) or zero for a point.
* **x**: input X coordinate to assign.
* **y**: input Y coordinate to assign.
"""
setpoint(geom::Geometry, i::Integer, x, y) =
    GDAL.setpoint_2d(geom, i, x, y)

"""
Add a point to a geometry (line string or point).

### Parameters
* **geom**: the geometry to add a point to.
* **x**: x coordinate of point to add.
* **y**: y coordinate of point to add.
* **z**: z coordinate of point to add.
"""
addpoint(geom::Geometry, x, y, z) = GDAL.addpoint(geom, x, y, z)

"""
Add a point to a geometry (line string or point).

### Parameters
* **geom**: the geometry to add a point to.
* **x**: x coordinate of point to add.
* **y**: y coordinate of point to add.
"""
addpoint(geom::Geometry, x, y) = GDAL.addpoint_2d(geom, x, y)

# """
#     OGR_G_SetPoints(OGRGeometryH hGeom,
#                     int nPointsIn,
#                     void * pabyX,
#                     int nXStride,
#                     void * pabyY,
#                     int nYStride,
#                     void * pabyZ,
#                     int nZStride) -> void
# Assign all points in a point or a line string geometry.
# ### Parameters
# * **hGeom**: handle to the geometry to set the coordinates.
# * **nPointsIn**: number of points being passed in padfX and padfY.
# * **pabyX**: list of X coordinates (double values) of points being assigned.
# * **nXStride**: the number of bytes between 2 elements of pabyX.
# * **pabyY**: list of Y coordinates (double values) of points being assigned.
# * **nYStride**: the number of bytes between 2 elements of pabyY.
# * **pabyZ**: list of Z coordinates (double values) of points being assigned (defaults to NULL for 2D objects).
# * **nZStride**: the number of bytes between 2 elements of pabyZ.
# """
# function setpoints(hGeom::Ptr{OGRGeometryH},nPointsIn::Integer,pabyX,nXStride::Integer,pabyY,nYStride::Integer,pabyZ,nZStride::Integer)
#     ccall((:OGR_G_SetPoints,libgdal),Void,(Ptr{OGRGeometryH},Cint,Ptr{Void},Cint,Ptr{Void},Cint,Ptr{Void},Cint),hGeom,nPointsIn,pabyX,nXStride,pabyY,nYStride,pabyZ,nZStride)
# end

"The number of elements in a geometry or number of geometries in container."
ngeom(geom::Geometry) = GDAL.getgeometrycount(geom)

"""
Fetch geometry from a geometry container.

### Parameters
* **geom**: the geometry container from which to get a geometry from.
* **i**: index of the geometry to fetch, between 0 and getNumGeometries() - 1.
"""
fetchgeom(geom::Geometry, i::Integer) =
    GDAL.getgeometryref(geom, i)

"""
Add a geometry to a geometry container.

### Returns
OGRERR_NONE if successful, or OGRERR_UNSUPPORTED_GEOMETRY_TYPE if the geometry
type is illegal for the type of existing geometry.
"""
function addgeom(geom::Geometry, subgeom::Geometry)
    result = GDAL.addgeometry(geom, subgeom)
    if result != GDAL.OGRERR_NONE
        error("Failed to add geometry. The geometry type could be illegal")
    end
end

"""
Add a geometry directly to an existing geometry container.

### Parameters
* **geom**: existing geometry.
* **subgeom**: geometry to add to the existing geometry.
"""
function addgeomdirectly(geom::Geometry, subgeom::Geometry)
    result = GDAL.addgeometrydirectly(geom, subgeom)
    if result != GDAL.OGRERR_NONE
        error("Failed to add geometry. The geometry type could be illegal")
    end
end

"""
Remove a geometry from an exiting geometry container.

### Parameters
* **geom**: the existing geometry to delete from.
* **i**: the index of the geometry to delete. A value of -1 is a special flag
    meaning that all geometries should be removed.
* **todelete**: if TRUE the geometry will be destroyed, otherwise it will not.
    The default is TRUE as the existing geometry is considered to own the
    geometries in it.
"""
function removegeom(geom::Geometry, i::Integer, todelete::Bool = true)
    result = GDAL.removegeometry(geom, i, todelete)
    if result != GDAL.OGRERR_NONE
        error("Failed to remove geometry. The index could be out of range")
    end
end

function removeallgeoms(geom::Geometry, todelete::Bool = true)
    result = GDAL.removegeometry(geom, -1, todelete)
    if result != GDAL.OGRERR_NONE
        error("Failed to remove geometries. The index could be out of range")
    end
end

"""
Returns if this geometry is or has curve geometry.

### Parameters
* **geom**: the geometry to operate on.
* **nonlinear**: set it to TRUE to check if the geometry is or contains a
    CIRCULARSTRING.
"""
hascurvegeom(geom::Geometry, nonlinear::Bool) =
    Bool(GDAL.hascurvegeometry(geom, nonlinear))

"""
Return, possibly approximate, linear version of this geometry.

### Parameters
* **geom**: the geometry to operate on.
* **stepsize**: the largest step in degrees along the arc, zero to use the
    default setting.
* **options**: options as a null-terminated list of strings or NULL.
    See OGRGeometryFactory::curveToLineString() for valid options.
"""
getlineargeom(geom::Geometry, stepsize::Real=0.0) =
    GDAL.getlineargeometry(geom, stepsize, C_NULL)

"Return curve version of this geometry."
getcurvegeom(geom::Geometry) = GDAL.getcurvegeometry(geom, C_NULL)

"""
Build a ring from a bunch of arcs.

### Parameters
* **lines**: handle to an OGRGeometryCollection (or OGRMultiLineString)
    containing the line string geometries to be built into rings.
* **besteffort**: not yet implemented???.
* **autoclose**: indicates if the ring should be close when first and last
    points of the ring are the same.
* **tol**: whether two arcs are considered close enough to be joined.
"""
function polygonfromedges(lines::Geometry, besteffort::Bool, autoclose::Bool,
                          tol::Real)
    peErr = Ref{GDAL.OGRErr}()
    result = GDAL.buildpolygonfromedges(lines, besteffort,
                                        autoclose, tol, peErr)
    if peErr[] != OGRERR_NONE
        error("Failed to build polygon from edges.")
    end
    result
end

# """
#     OGRSetNonLinearGeometriesEnabledFlag(int bFlag) -> void
# Set flag to enable/disable returning non-linear geometries in the C API.
# ### Parameters
# * **bFlag**: TRUE if non-linear geometries might be returned (default value). FALSE to ask for non-linear geometries to be approximated as linear geometries.
# ### Returns
# a point or NULL.
# """
# function setnonlineargeometriesenabledflag(bFlag::Integer)
#     ccall((:OGRSetNonLinearGeometriesEnabledFlag,libgdal),Void,(Cint,),bFlag)
# end

# """
#     OGRGetNonLinearGeometriesEnabledFlag(void) -> int
# Get flag to enable/disable returning non-linear geometries in the C API.
# """
# function getnonlineargeometriesenabledflag()
#     ccall((:OGRGetNonLinearGeometriesEnabledFlag,libgdal),Cint,())
# end


text:
"""
    API_REGISTRY_NAME::Symbol

The name of a global constant that holds a list `Vector{Symbol}` of public names
in each module.
"""
const API_REGISTRY_NAME = gensym(:API_REGISTRY_NAME)

"""
    @public name₁ name₂ … nameₙ
    @public @macroname

Declare public API names for the current module.

The second form `@public @macroname` is equivalent to `@public var"@macroname"`.

# Extended help

Note that where this macro is invoked is important. Consider:

```julia
module A1
    using PublicAPI: @public
    @public B, C
    module B
        f() = nothing
    end
    module C
        using PublicAPI: @public
        using ..B: f
        @public f
    end
end
```

and

```julia
module A2
    using PublicAPI: @public
    @public B, C
    module B
        using PublicAPI: @public
        f() = nothing
        @public f
    end
    module C
        using ..B: f
    end
end
```

The fully-qualified names `A1.C.f` and `A2.B.f` are public but `A1.B.f` and
`A2.C.f` are private.
"""
macro public(name::Symbol, othernames::Symbol...)
    return public_impl(__module__, [name, othernames...])
end

macro public(macrocall::Expr)
    if macrocall.head !== :macrocall
        error("Unsupported syntax: $macrocall")
    elseif length(macrocall.args) != 2
        error("Expected single macro name as in `@public @macro`; got: $macrocall")
    end
    return public_impl(__module__, Symbol[macrocall.args[1]])
end

function public_impl(__module__::Module, allnames::Vector{Symbol})
    @gensym namespace
    expr = quote
        # Allocate the "registry" (just a `Vector{Symbol}`) if it hasn't been
        # defined yet:
        if $Base.:!($Base.@isdefined $API_REGISTRY_NAME)
            const $API_REGISTRY_NAME = $Symbol[]
        end

        # Register the API names:
        $union!($API_REGISTRY_NAME, $(QuoteNode(allnames)))

        # Check that these names are defined after the module is finalized
        module $namespace
        __init__() = $check_public_names($__module__, $(QuoteNode(allnames)))
        end
    end
    @assert expr.head === :block
    return esc(Expr(:toplevel, expr.args...))
end

function check_public_names(__module__::Module, allnames::Vector{Symbol})
    undefined = filter(allnames) do n
        try
            getproperty(__module__, n)
            false
        catch
            true
        end
    end
    if !isempty(undefined)
        are = length(undefined) == 1 ? "is" : "are"
        println(stderr, "ERROR: Following API $are declared but not defined")
        for n in undefined
            println(stderr, "  $__module__.$n")
        end
        error("PublicAPI: $(length(undefined)) undefined public API")
    end
end


text:
using DifferentialForms

using StaticArrays
using Test

@testset "bitsign" begin
    for b in false:true
        @test abs(bitsign(b)) === 1
        @test signbit(bitsign(b)) === b
        @test bitsign(Int(b)) === bitsign(b)
    end
    for n in 1:100
        i = Int(rand(Int8))
        j = Int(rand(Int8))
        @test bitsign(i) * bitsign(j) == bitsign(i + j)
    end
end

@testset "sort_perm (StaticArray)" begin
    @test sort_perm(SVector{0,Int}()) == (SVector{0,Int}(), 0)
    @test sort_perm(SVector(1)) == (SVector(1), 0)
    @test sort_perm(SVector(1, 2)) == (SVector(1, 2), 0)
    @test sort_perm(SVector(2, 1)) == (SVector(1, 2), 1)
    @test sort_perm(SVector(1, 2, 3)) == (SVector(1, 2, 3), 0)
    @test sort_perm(SVector(1, 3, 2)) == (SVector(1, 2, 3), 1)
    @test sort_perm(SVector(2, 1, 3)) == (SVector(1, 2, 3), 1)
    @test sort_perm(SVector(2, 3, 1)) == (SVector(1, 2, 3), 2)
    @test sort_perm(SVector(3, 1, 2)) == (SVector(1, 2, 3), 2)
    @test sort_perm(SVector(3, 2, 1)) == (SVector(1, 2, 3), 3)

    for N in 0:10
        xs = rand(SVector{N,Int})
        ys, s = sort_perm(xs)
        @test issorted(ys)
        @test s >= 0
    end

    for N in 2:10
        xs = SVector{N}(1:N)
        n = rand(0:10)
        # permute n times
        for i in 1:n
            j = rand(1:(N - 1))
            xsj = xs[j]
            xsj1 = xs[j + 1]
            xs = setindex(xs, xsj1, j)
            xs = setindex(xs, xsj, j + 1)
        end
        ys, s = sort_perm(xs)
        @test issorted(ys)
        @test iseven(s - n)
    end
end

@testset "sort_perm (AbstractVector)" begin
    @test sort_perm(Vector{Int}()) == (Vector{Int}(), 0)
    @test sort_perm([1]) == ([1], 0)
    @test sort_perm([1, 2]) == ([1, 2], 0)
    @test sort_perm([2, 1]) == ([1, 2], 1)
    @test sort_perm([1, 2, 3]) == ([1, 2, 3], 0)
    @test sort_perm([1, 3, 2]) == ([1, 2, 3], 1)
    @test sort_perm([2, 1, 3]) == ([1, 2, 3], 1)
    @test sort_perm([2, 3, 1]) == ([1, 2, 3], 2)
    @test sort_perm([3, 1, 2]) == ([1, 2, 3], 2)
    @test sort_perm([3, 2, 1]) == ([1, 2, 3], 3)

    for N in 0:10
        xs = rand(Int, N)
        ys, s = sort_perm(xs)
        @test issorted(ys)
        @test s >= 0
    end

    for N in 2:10
        xs = collect(1:N)
        n = rand(0:10)
        # permute n times
        for i in 1:n
            j = rand(1:(N - 1))
            xs[j], xs[j + 1] = xs[j + 1], xs[j]
        end
        ys, s = sort_perm(xs)
        @test issorted(ys)
        @test iseven(s - n)
    end
end


text:
using Test
using QuantumLattices.Prerequisites.SimpleTrees

@simpletree struct ATree end
@simpletree struct BTree end {<:AbstractString,::Int}

@testset "AbstractSimpleTree" begin
    tree=ATree{String,Int}()

    @test tree|>eltype==Pair{String,Int}
    @test tree|>keytype==String
    @test tree|>valtype==Int
    @test tree|>root==nothing

    push!(tree,"L0",0)
    push!(tree,"L0","L1-1",2)
    push!(tree,"L0","L1-2",3)
    push!(tree,"L1-1","L2-1",4)
    push!(tree,"L1-1","L2-2",5)
    push!(tree,"L1-2","L2-3",6)
    push!(tree,"L1-2","L2-4",7)
    tree["L0"]=1

    @test tree|>root=="L0"
    @test tree|>length==7
    @test keys(tree,simpletreedepth)|>collect==["L0","L1-1","L2-1","L2-2","L1-2","L2-3","L2-4"]
    @test keys(tree,simpletreewidth)|>collect==["L0","L1-1","L1-2","L2-1","L2-2","L2-3","L2-4"]
    @test values(tree,simpletreedepth)|>collect==[1,2,4,5,3,6,7]
    @test values(tree,simpletreewidth)|>collect==[1,2,3,4,5,6,7]
    @test pairs(tree,simpletreedepth)|>collect==[("L0",1),("L1-1",2),("L2-1",4),("L2-2",5),("L1-2",3),("L2-3",6),("L2-4",7)]
    @test pairs(tree,simpletreewidth)|>collect==[("L0",1),("L1-1",2),("L1-2",3),("L2-1",4),("L2-2",5),("L2-3",6),("L2-4",7)]
    @test [isleaf(tree,node) for node in keys(tree,simpletreewidth)]==[false,false,false,true,true,true,true]
    @test [level(tree,node) for node in keys(tree,simpletreewidth)]==[1,2,2,3,3,3,3]
    @test [ancestor(tree,"L2-1",i) for i=0:2]==["L2-1","L1-1","L0"]
    @test [ancestor(tree,"L2-3",i) for i=0:2]==["L2-3","L1-2","L0"]
    @test descendants(tree,"L0",0)==["L0"]
    @test descendants(tree,"L0",1)==["L1-1","L1-2"]
    @test descendants(tree,"L0",2)==["L2-1","L2-2","L2-3","L2-4"]
    @test siblings(tree,"L0")==[]
    @test siblings(tree,"L1-1")==["L1-2"]
    @test siblings(tree,"L2-1")==["L2-2"]
    @test leaves(tree)==["L2-1","L2-2","L2-3","L2-4"]

    backup=deepcopy(tree)

    sub=subtree(tree,"L0")
    empty!(tree)
    append!(tree,sub)
    @test tree==backup
    @test isequal(tree,backup)

    sub=subtree(tree,"L1-1")
    delete!(tree,"L1-1")
    append!(tree,"L0",sub)
    move!(tree,"L1-2","L0")
    @test tree==backup
    @test isequal(tree,backup)
end

@testset "SimpleTree" begin
    tree=SimpleTree{String,Int}()
    @test tree|>eltype==tree|>typeof|>eltype==Pair{String,Int}
    @test tree|>keytype==tree|>typeof|>keytype==String
    @test tree|>valtype==tree|>typeof|>valtype==Int
    @test tree|>root==nothing
end

@testset "@simpletree" begin
    @test BTree|>valtype==Int
    @test BTree{String}|>keytype==String
    @test BTree{String}|>valtype==Int
    @test BTree{String}().TREECORE==SimpleTreeCore{String,Int}()
end


text:
module HiGHS

import HiGHS_jll: libhighs

include("wrapper/ctypes.jl")
include("wrapper/libhighs_common.jl")
include("wrapper/libhighs_api.jl")

include("MOI_wrapper.jl")

# HiGHS exports all `Highs_xxx` symbols. If you don't want all of these symbols
# in your environment, then use `import HiGHS` instead of `using HiGHS`.

for sym in names(@__MODULE__, all = true)
    if startswith(string(sym), "Highs_")
        @eval export $sym
    end
end

end # module HiGHS


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function contains(ustr, s)

itsin = false;
if occursin(s, string(ustr))
    itsin = true;
end

return itsin;

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
# testing ERB filterbank design
fb = make_erb_filterbank(16000, 23, 150)
fb_matlab = readcsv(open(joinpath(dirname(@__FILE__), "data", "ERB_filter_coeffs.csv")))
@test length(fb.filters) == size(fb_matlab, 1)
for k=1:length(fb.filters)
	sos = fb.filters[k]
	biquads = sos.biquads
	g = sos.g
	for b=1:4
		@test_approx_eq biquads[b].b0 fb_matlab[k,1]
		@test_approx_eq biquads[b].b2 fb_matlab[k,6]
		@test_approx_eq biquads[b].a1 fb_matlab[k,8]
		@test_approx_eq biquads[b].a2 fb_matlab[k,9]
	end
	@test_approx_eq biquads[1].b1 fb_matlab[k,2]
	@test_approx_eq biquads[2].b1 fb_matlab[k,3]
	@test_approx_eq biquads[3].b1 fb_matlab[k,4]
	@test_approx_eq biquads[4].b1 fb_matlab[k,5]
	@test_approx_eq 1/g fb_matlab[k,10]
end

# testing ERB filterbank response
x = zeros(10000)
x[1] = 1.0
y = filt(fb, x)
y_matlab = readcsv(open(joinpath(dirname(@__FILE__), "data", "ERB_filter_response.csv")))
@test_approx_eq y y_matlab'


text:
#############
# Repo Type #
#############

type Repo <: GitHubType
    name::Nullable{GitHubString}
    full_name::Nullable{GitHubString}
    description::Nullable{GitHubString}
    language::Nullable{GitHubString}
    default_branch::Nullable{GitHubString}
    owner::Nullable{Owner}
    parent::Nullable{Repo}
    source::Nullable{Repo}
    id::Nullable{Int}
    size::Nullable{Int}
    subscribers_count::Nullable{Int}
    forks_count::Nullable{Int}
    stargazers_count::Nullable{Int}
    watchers_count::Nullable{Int}
    open_issues_count::Nullable{Int}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    homepage::Nullable{HttpCommon.URI}
    pushed_at::Nullable{Dates.DateTime}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    has_issues::Nullable{Bool}
    has_wiki::Nullable{Bool}
    has_downloads::Nullable{Bool}
    has_pages::Nullable{Bool}
    private::Nullable{Bool}
    fork::Nullable{Bool}
    permissions::Nullable{Dict}
end

Repo(data::Dict) = json2github(Repo, data)
Repo(full_name::AbstractString) = Repo(Dict("full_name" => full_name))

namefield(repo::Repo) = repo.full_name

###############
# API Methods #
###############

# repos #
#-------#

function repo(repo_obj; options...)
    result = gh_get_json("/repos/$(name(repo_obj))"; options...)
    return Repo(result)
end

# forks #
#-------#

function forks(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/forks"; options...)
    return map(Repo, results), page_data
end

function create_fork(repo; options...)
    result = gh_post_json("/repos/$(name(repo))/forks"; options...)
    return Repo(result)
end

# contributors/collaborators #
#----------------------------#

function contributors(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/contributors"; options...)
    results = [Dict("contributor" => Owner(i), "contributions" => i["contributions"]) for i in results]
    return results, page_data
end

function collaborators(repo; options...)
    results, page_data = gh_get_json("/repos/$(name(repo))/collaborators"; options...)
    return map(Owner, results), page_data
end

function iscollaborator(repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    r = gh_get(path; handle_error = false, options...)
    r.status == 204 && return true
    r.status == 404 && return false
    handle_response_error(r)  # 404 is not an error in this case
    return false
end

function add_collaborator(repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    return gh_put(path; options...)
end

function remove_collaborator(repo, user; options...)
    path = "/repos/$(name(repo))/collaborators/$(name(user))"
    return gh_delete(path; options...)
end

# stats #
#-------#

function stats(repo, stat, attempts = 3; options...)
    path = "/repos/$(name(repo))/stats/$(name(stat))"
    local r
    for a in 1:attempts
        r = gh_get(path; handle_error = false, options...)
        r.status == 200 && return r
        sleep(2.0)
    end
    return r
end


text:
using BlockArrays
using ControlSystems
using LinearAlgebra
using PredictiveControl
using Test

const MPC = PredictiveControl

# Some utilities for testing
include( "../../src/utilities.jl" )
include( "../testUtils.jl" )

# Create a sample system
A = [1.0 0.0 0.0 0.0;
     1.0 2.0 0.0 0.0;
     0.0 0.0 3.0 0.0;
     0.0 0.0 1.0 4.0]
B = [1.0 0.0;
     0.0 0.0;
     0.0 1.0;
     0.0 0.0]
C = to_matrix( Float64, I, 4 )

sys = StateSpace( A, B, C, 0, 0.1 )

Q = diagm( [1.0, 2.0, 3.0, 4.0] )
R = diagm( [1.0, 2.0] )

N  = 5
nx = ControlSystems.nstates( sys )
nu = ControlSystems.ninputs( sys )

###################################################################
# Test the matrices of an uncontrolled system
###################################################################
clqr = ConstrainedTimeInvariantLQR( sys, N, Q, R, :Q )

# Form a propagation matrix to test
H = MPC.hessian( clqr, PrimalFullCondensing )
J = MPC.linearcoefficients( clqr, PrimalFullCondensing )

@test blocksize( H ) == (N, N)
@test size( H ) == (nu*N, nu*N)

@test blocksize( J ) == (N, 1)
@test size( J ) == (nu*N, nx)

# Make sure the blocks have the proper values for the system
for i = 1:N
#    b = view( mat, Block( i, 1 ) )

#    @test size( b ) == (nx, nx)
#    @test b == sys.A^i
end

###################################################################
# Test the matrices of a controlled system
###################################################################
clqr = ConstrainedTimeInvariantLQR( sys, N, Q, R, :Q, K = :dlqr )

# Form a propagation matrix to test
H = MPC.hessian( clqr, PrimalFullCondensing )
J = MPC.linearcoefficients( clqr, PrimalFullCondensing )

@test blocksize( H ) == (N, N)
@test size( H ) == (nu*N, nu*N)

@test blocksize( J ) == (N, 1)
@test size( J ) == (nu*N, nx)

# Make sure the blocks have the proper values for the system
for i = 1:N
#    b = view( mat, Block( i, 1 ) )

#    @test size( b ) == (nx, nx)
#    @test b == sys.A^i
end


text:
#!/usr/bin/env julia
using Covid19Inc, PGFPlotsX
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
p = Covid19Inc.covid_pop(724, :i10o90, "2020-04-13") |> Covid19Inc.covidpl
pgfsave("../../images/724_covid19_2020-04-13.svg", p)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
abstract type TracerBC end

"""
    ImpermeableTracer() :: TracerBC

No tracer diffusion across boundary
"""
struct ImpermeableTracer <: TracerBC end
function atmos_tracer_boundary_state!(
    nf,
    bc_tracer::ImpermeableTracer,
    atmos,
    args...,
)
    nothing
end
function atmos_tracer_normal_boundary_flux_second_order!(
    nf,
    bc_tracer::ImpermeableTracer,
    atmos,
    args...,
)
    nothing
end


text:
function genericPDGSBP(x, maxiter, burnin, seed, filename = "/results")
  ### set seed and create a folder in the "savelocation" path to store results
  srand(seed)
  savelocation = string(pwd(), filename, "/seed$seed")
  mkpath(savelocation)

  for i in 1:1size(x)[1]
    writedlm(string(savelocation, "/x$i.txt"), x[i])
  end
  n = [length(x[j]) for j in 1:1:size(x)[1]]

  ### initialization
  NSAMPLES = length(n)
  epsilons = eye(Int64, NSAMPLES)
  idx = Int(maxiter/1000)
  times_elapsed = Array(Float64, idx)

  ### initialize clustering and geometric-slice variables
  ### each collection of data x[j] has its corresponding clustering variables d[j] and slice var N[j]
  # the delta[j] indicate the infinite mixture within the finite mixture each observation came from
  d = cell(NSAMPLES)
  N = cell(NSAMPLES)
  delta = cell(NSAMPLES)
  for j in 1:1:NSAMPLES
    d[j] = rand(DiscreteUniform(1,5),n[j])
    delta[j] = rand(Multinomial(1, 1/NSAMPLES * ones(NSAMPLES)), n[j])'
    N[j] = copy(d[j])
  end

  ### set up prior hyperparameters of the geometric stick breaking base measure G_0
  ### here we set as G_0 the independent NormalxGamma measure
  ### prior hyperparameter of the mean component
  mu0 = 0.0
  tau0 = 0.001
  sigma0 = sqrt(1.0 / tau0)

  ### prior hyperparameters of the precision components
  at = 0.001
  bt = 0.001

  ### initialize mixture selection probabilities and Dirichlet hyperprior alpha
  p = (1/NSAMPLES) * ones(NSAMPLES, NSAMPLES)
  alphajl = ones(NSAMPLES, NSAMPLES)
  P = cell(maxiter)

  ### Initialize geometric probability and hyperparameters of the conjugate beta prior
  z = 0.5
  az = 1.0
  bz = 1.0

  ### construct a symmetric matrix of random measures with pairwise dependence
  ### each element of the matrix is a Geometric stick breaking mixture
  Q = Matrix(NSAMPLES, NSAMPLES)
  M = maximum([maximum(N[r]) for r in 1:1:NSAMPLES])
  for j in 1:1:NSAMPLES
    for i in 1:1:j
      Q[i, j] = GSB(zeros(M), ones(M), z, zeros(M), zeros(M), zeros(maxiter))
    end
  end
  Q = Symmetric(Q)

  ### store clusters and samples from predictives in arrays
  fclusters = Array(Int64, NSAMPLES, maxiter)
  pred = Array(Float64, NSAMPLES, maxiter - burnin)

  println("Starting Gibbs...")
  start_time = tic()

  for its in 1:1:maxiter
    M = maximum([maximum(N[r]) for r in 1:1:NSAMPLES])

    for l in 1:1:NSAMPLES
      for j in 1:1:l
        Q[j,l].mu = zeros(M)
        Q[j,l].w = zeros(M)
        Q[j,l].W = zeros(M)
        for k in 1:1:M
          Q[j,l].w[k] = Q[j,l].z * (1 - Q[j,l].z)^(k - 1)

          ### sampling of locations for the idiosyncratic
          # sampling of the means
          ksij, mj, gj = 0.0, 0.0, 0.0
          if j .== l
            for i in 1:1:n[j]
              if (d[j][i] .== k) & (delta[j][i,l] .== 1)
                mj += 1
                ksij += x[j][i]
              end
            end
          else # sample the locations for the common parts between densities fⱼ and fl
            for i in 1:1:n[j]
              if (d[j][i] .== k) & (delta[j][i,l] .== 1)
                mj += 1
                ksij += x[j][i]
              end
            end
            for i in 1:1:n[l]
              if (d[l][i] .== k) & (delta[l][i,j] .== 1)
                mj += 1
                ksij += x[l][i]
              end
            end
          end # end if i==j
          meanstar = Q[j,l].lambda[k] * ksij / (tau0 + Q[j,l].lambda[k] * mj)
          varstar = 1.0 / (tau0 + Q[j,l].lambda[k] * mj)
          Q[j,l].mu[k] = rand(Normal(meanstar, varstar^0.5))

          ### sampling of location for the idiosyncratic parts
          # sampling of precisions
          if j .== l
            for i in 1:1:n[j]
              if (d[j][i] .== k) & (delta[j][i,l] .== 1)
                gj += (x[j][i] - Q[j,l].mu[k])^2
              end
            end
          else # sample locations for the common parts
            for i in 1:1:n[j]
              if (d[j][i] .== k) & (delta[j][i,l] .== 1)
                gj += (x[j][i] - Q[j,l].mu[k])^2
              end
            end
            for i in 1:1:n[l]
              if (d[l][i] .== k) & (delta[l][i,j] .== 1)
                gj += (x[l][i] - Q[j,l].mu[k])^2
              end
            end
          end # if j==l
          alphastar = at + 0.5 * mj
          betastar = bt + 0.5 * gj
          Q[j,l].lambda[k] = rand(Gamma(alphastar, betastar^(-1)))

        end # k ind
      end
    end

    for j in 1:1:NSAMPLES
      d[j], delta[j] = sampleClustAllocs(Q[j,:], x[j], p[j,:], N[j]) # sample clustering variables
      fclusters[j, its] = uniqueClusters(d[j], delta[j]) # unique clusters on j density at its iteration
      N[j] = sampleGeomSlice(Q[j,:], d[j], delta[j])  # sample geometric slice variables
      p[j, :] = sampleSelectionProbabilities(alphajl[j, :], delta[j]) # update the mixture selection probabilities
    end
    P[its] = copy(p)

    if its > burnin
      for j in 1:1:NSAMPLES
        ### density estimation
        pred[j, its - burnin] = samplePredictive(Q[j,:], p[j,:], mu0, sigma0, at, bt) # sample from the predictive density
      end
    end

    ### update geometric probability
    for l in 1:1:NSAMPLES
      for j in 1:1:l
        if j == l # for the idiosyncratic part
          Q[j,l].z = rand(Beta(az + 2 * sum(delta[j][:, l].==1), bz + sum(N[j][delta[j][:,l].==1]) - sum(delta[j][:, l].==1)))
        else  # and for the common parts
          Q[j,l].z = rand(Beta(az + 2 * (sum(delta[j][:, l].==1) + sum(delta[l][:, j].==1)),
             bz + (sum(N[j][delta[j][:,l].==1]) + sum(N[l][delta[l][:,j].==1])) - (sum(delta[j][:, l].==1) + sum(delta[l][:,j].==1))))
        end
        Q[j,l].Z[its] = Q[j,l].z
      end
    end

    ### adjust array of precision for next iteration
    if maximum([maximum(N[r]) for r in 1:1:NSAMPLES]) > M
      for l in 1:1:NSAMPLES
        for j in 1:1:l
          Q[j, l].lambda = [Q[j, l].lambda;ones(maximum([maximum(N[r]) for r in 1:1:NSAMPLES]) - M)]
        end
      end
    end

    if its%1000 == 0
      end_time = toc()
      times_elapsed[Int(its/1000)] = end_time
      println("MCMC iterations: $its", " in ", end_time, " seconds")
      start_time = tic()
    end

  end ### end gibbs
  writedlm(string(savelocation, "/predictives.txt"), pred')
  writedlm(string(savelocation, "/clusters.txt"), fclusters')
  writedlm(string(savelocation, "/times.txt"), times_elapsed)
  writedlm(string(savelocation, "/Prob.txt"), P)

  for l in 1:1:NSAMPLES
    for j in 1:1:l
      writedlm(string(savelocation, "/L$l$j.txt"), Q[j,l].Z)
    end
  end

  return pred, P, fclusters, times_elapsed
end

max_stars_repo_name: cmerkatas/PDGSBP

max_issues_repo_name: cmerkatas/PDGSBP

max_forks_repo_name: cmerkatas/PDGSBP


text:
@inline function _line_interp_weights(x, x1, x2)
  lx = x2 - x1
  dx = (x - x1)/lx
  return (1 - dx, dx)
end

@inline function _tri_interp_weights(x, x1, x2, y, y1, y2)
  lx = (x2 - x1)
  dx = (x - x1)/lx

  ly = (y2 - y1)
  dy = (y - y1)/ly

  return (1 - dx, -dy + dx, dy)
end

@inline function _square_interp_weights(x, x1, x2, y, y1, y2)
  lx = (x2 - x1)
  dx = (x - x1)/lx

  ly = (y2 - y1)
  dy = (y - y1)/ly

  w = ((1 - dx) * (1 - dy),
       (1 - dx) * dy,
       dx * (1 - dy),
       dx * dy)

  return w
end

@inline function _line_interp(x, x1, x2, f1, f2)
  w = _line_interp_weights(x, x1, x2)
  return w[1] * f1 + w[2] * f2
end

@inline function _line_interp!(f, x, x1, x2, f1, f2)
  w = _line_interp_weights(x, x1, x2)
  return f .= w[1] .* f1 .+ w[2] .* f2
end

@inline function _tri_interp(x, x1, x2, y, y1, y2, v1, v2, v3)
  w = _tri_interp_weights(x, x1, x2, y, y1, y2)
  return w[1] * v1 + w[2] * v2 + w[3] * v3
end

@inline function _tri_interp!(f, x, x1, x2, y, y1, y2, v1, v2, v3)
  w = _tri_interp_weights(x, x1, x2, y, y1, y2)
  return f .= w[1] .* v1 .+ w[2] .* v2 .+ w[3] .* v3
end

@inline function _square_interp(x, x1, x2, y, y1, y2, f11, f12, f21, f22)
  w = _square_interp_weights(x, x1, x2, y, y1, y2)
  return w[1] * f11 + w[2] * f12 + w[3] * f21 + w[4] * f22
end

@inline function _square_interp!(f, x, x1, x2, y, y1, y2, f11, f12, f21, f22)
  w = _square_interp_weights(x, x1, x2, y, y1, y2)
  return f .= w[1] .* f11 .+ w[2] .* f12 .+ w[3] .* f21 .+ w[4] .* f22
end


text:
failed_to_open_str = "Arduino SerialPort is not connected. Call open(ard::Arduino) first"

arduino_buffer_size = 1020
msg_block_size = 256

mutable struct Arduino # SerialSocket
    portname::String
    baudrate::Int64
    sp::Union{Nothing, SerialPort}

    # Buffer for dumping read bytes into
    read_buffer::MVector{arduino_buffer_size, UInt8}

    # Buffer for encoded messages
    msg_out_buffer::MVector{msg_block_size, UInt8}
    msg_out_length::Int64

    # Buffer for decoded messages
    msg_in_buffer::MVector{msg_block_size, UInt8}
    msg_in_length::Int64

    # ProtoBuf IO stream
    proto_iob::IOBuffer

    function Arduino(portname::String, baudrate::Int)
        sp = nothing

        read_buffer = @MVector zeros(UInt8, arduino_buffer_size)

        msg_out_buffer = @SVector zeros(UInt8, msg_block_size)
        msg_out_length = 0

        msg_in_buffer = @SVector zeros(UInt8, msg_block_size)
        msg_in_length = 0

        proto_iob = IOBuffer(zeros(UInt8, msg_block_size), read=true, write=true,
                            #  maxsize=msg_block_size,
                             append=false
                             )

        return new(portname, baudrate, sp,
                   read_buffer,
                   msg_out_buffer, msg_out_length,
                   msg_in_buffer, msg_in_length,
                   proto_iob)
    end
end

####################################################################
#                        Serial port I/O
####################################################################
function Base.close(ard::Arduino)::Nothing
    if ard.sp === nothing
        return
    end
    if isopen(ard)
        close(ard.sp)
    end
end

function Base.isopen(ard::Arduino)::Bool
    if ard.sp === nothing
        return false
    end

    return isopen(ard.sp)
end

function Base.open(ard::Arduino)::SerialPort
    close(ard)
    ard.sp = LibSerialPort.open(ard.portname, ard.baudrate)
    return ard.sp
end

function Base.bytesavailable(ard::Arduino)
    isopen(ard) || error(failed_to_open_str)

    return bytesavailable(ard.sp)
end

function messageBytesAvailable(ard::Arduino)
    isopen(ard) || error(failed_to_open_str)

    if bytesavailable(ard) > 2
        return bytesavailable(ard) - 2
    else
        return bytesavailable(ard)
    end
end

# function Base.read!(ard::Arduino, array::AbstractVector{UInt8})
function read_into_buffer!(ard::Arduino)
    isopen(ard) || error(failed_to_open_str)

    if Base.bytesavailable(ard) <= arduino_buffer_size
        bytes_read = Base.bytesavailable(ard)
    else
        bytes_read = Base.bytesavailable(arduino_buffer_size)
    end

    readbytes!(ard.sp, ard.read_buffer, bytes_read)
    return bytes_read
end

function retrieve_encoded_msg(ard::Arduino)
    isopen(ard) || error(failed_to_open_str)

    port_buffered_bytes = Base.bytesavailable(ard)

    if msg_block_size < port_buffered_bytes
        readbytes!(ard.sp, ard.read_buffer, port_buffered_bytes)

        last_msg_last_byte = findfirst(x -> x == 0x00, ard.read_buffer)
        next_msg_last_byte = findnext(x -> x == 0x00, ard.read_buffer, last_msg_last_byte + 1)

        if next_msg_last_byte !== nothing
            encoded_msg = @view ard.read_buffer[last_msg_last_byte+1:next_msg_last_byte]

            return encoded_msg
        end
    end

    return nothing
end

function Base.readuntil(ard::Arduino, delim)
    read_into_buffer!(ard)
    last_ind = findfirst(x -> x==delim, ard.read_buffer)

    if last_ind isa Int64
        return view(ard.read_buffer, 1:last_ind)
    else
        return ard.read_buffer
    end
end

function Base.write(ard::Arduino, x)
    isopen(ard) || error(failed_to_open_str)

    return write(ard.sp, x)
end


text:
module TestAqua

import Setfield
using Aqua
using BangBang
using InitialValues
using Test

Aqua.test_all(InitialValues)

@testset "Test ambiguities with `BangBang`" begin
    if v"1.6.0-DEV.816" <= VERSION < v"1.6.0-DEV.875"
        # Maybe remove this branch?
        @warn "Ignoring ambiguities from `Base` to workaround JuliaLang/julia#36962"
        packages = [InitialValues, BangBang]
    else
        packages = [InitialValues, BangBang, Base]
    end
    Aqua.test_ambiguities(packages; exclude = [Base.get, Setfield.set, Setfield.modify])
end

@testset "Compare test/Project.toml and test/environments/main/Project.toml" begin
    @test Text(read(joinpath(@__DIR__, "Project.toml"), String)) ==
          Text(read(joinpath(@__DIR__, "environments", "main", "Project.toml"), String))
end

end  # module


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allapprox(v1,v2,tol)
    flag = length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=tol) for i = 1:length(v1)])
end

vbasis = TensorProductBasis(2,1)
sbasis = TensorProductBasis(1,1)
vquad = tensor_product_quadrature(2,2)
squad = tensor_product_quadrature(1,2)
facequads = repeat([squad],4)
facemaps = HDGElasticity.reference_cell_facemaps(2)
normals = HDGElasticity.reference_normals()
lambda = 1.0
mu = 2.0
Dhalf = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
facescale = HDGElasticity.face_determinant_jacobian(cellmap)
stabilization = 0.1

locop = HDGElasticity.local_operator(vbasis,vquad,facequads,facemaps,
    Dhalf,stabilization,cellmap)
facelhops = HDGElasticity.local_hybrid_operator(vbasis,sbasis,facequads,facemaps,
    normals,Dhalf,stabilization,cellmap)
lhop = hcat(facelhops...)

tcomp = [1.0,0.0]
hybloc1 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,
    vbasis,facequads[1],facemaps[1],normals[1],tcomp,Dhalf,
    stabilization,facescale[1])
HHt1 = stabilization*HDGElasticity.HHop(sbasis,squad,tcomp,facescale[1])
ucomp = [0.0,-1.0]
HHu1 = HDGElasticity.HHop(sbasis,squad,ucomp,facescale[1])

hybloc2 = facelhops[2]'
R2 = -HDGElasticity.linear_form(facescale[2]*[1.0,0.0],sbasis,squad)
HH2 = stabilization*HDGElasticity.HHop(sbasis,facequads[2],facescale[2])

hybloc3 = facelhops[3]'
HH3 = stabilization*HDGElasticity.HHop(sbasis,facequads[3],facescale[3])

tcomp = [0.,1.]
hybloc4 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,
    vbasis,facequads[4],facemaps[4],normals[4],tcomp,Dhalf,
    stabilization,facescale[4])
HHt4 = stabilization*HDGElasticity.HHop(sbasis,squad,tcomp,facescale[4])
ucomp = [-1.0,0.0]
HHu4 = HDGElasticity.HHop(sbasis,squad,ucomp,facescale[4])

hlop = vcat(hybloc1,hybloc2,hybloc3,hybloc4)

hhop = zeros(16,16)
hhop[1:4,1:4] = HHt1 - HHu1
hhop[5:8,5:8] = HH2
hhop[9:12,9:12] = HH3
hhop[13:16,13:16] = HHt4 - HHu4

cellop = hlop*(locop\lhop) - hhop
rhs = zeros(16)
rhs[5:8] = R2

sol = cellop\rhs
Hsol = reshape(sol,2,:)

locsol = locop\(lhop*sol)
L = -Dhalf*reshape(locsol[1:12],3,:)
U = reshape(locsol[13:20],2,:)

E = (lambda+2mu) - lambda^2/(lambda+2mu)
e11 = 1/E
e22 = -lambda/(lambda+2mu)*e11
u1 = e11*2
u2 = e22*1

testHsol = zeros(2,8)
testHsol[1,[2,3,4,6]] .= u1
testHsol[2,[4,5,6,8]] .= u2

testL = zeros(3,4)
testL[1,:] .= 1.0

testU = zeros(2,4)
testU[2,2] = u2
testU[1,3] = u1
testU[1,4] = u1
testU[2,4] = u2

@test allapprox(testL,L,1e-10)
@test allapprox(testU,U,1e-10)
@test allapprox(Hsol,testHsol,1e-10)

hybloc1 = zeros(4,20)
HH1 = HDGElasticity.HHop(sbasis,facequads[1],facescale[1])

hybloc2 = facelhops[2]'
R2 = -HDGElasticity.linear_form(facescale[2]*[0.0,1.0],sbasis,facequads[2])
HH2 = stabilization*HDGElasticity.HHop(sbasis,facequads[2],facescale[2])

hybloc3 = facelhops[3]'
R3 = -HDGElasticity.linear_form(facescale[3]*[1.0,0.0],sbasis,facequads[3])
HH3 = stabilization*HDGElasticity.HHop(sbasis,facequads[3],facescale[3])

hybloc4 = facelhops[4]'
R4 = -HDGElasticity.linear_form(facescale[4]*[0.0,-1.0],sbasis,facequads[4])
HH4 = stabilization*HDGElasticity.HHop(sbasis,facequads[4],facescale[4])

hlop = vcat(hybloc1,hybloc2,hybloc3,hybloc4)

hhop = zeros(16,16)
hhop[1:4,1:4] = -HH1
hhop[5:8,5:8] = HH2
hhop[9:12,9:12] = HH3
hhop[13:16,13:16] = HH4

cellop = hlop*(locop\lhop) - hhop
rhs = zeros(16)
rhs[5:8] = R2
rhs[9:12] = R3
rhs[13:16] = R4

sol = cellop\rhs
Hsol = reshape(sol,2,:)

locsol = locop\(lhop*sol)
L = -Dhalf*reshape(locsol[1:12],3,:)
U = reshape(locsol[13:20],2,:)

e12 = 1.0/(2mu)
gamma = 2*e12

testL = zeros(3,4)
testL[3,:] .= 1.0

testU = zeros(2,4)
testU[1,2] = gamma
testU[1,4] = gamma

testH = zeros(2,8)
testH[1,4] = gamma
testH[1,5] = gamma
testH[1,6] = gamma
testH[1,8] = gamma

@test allapprox(Hsol,testH,1e-10)
@test allapprox(L,testL,1e-10)
@test allapprox(U,testU,1e-10)