text:
using SparseArrays
using Statistics
using LinearAlgebra
using Arpack
"""
    (pm,v) = partition_fiedler(W; method, v)

Partition the vertices of a graph according to the Fiedler vector

### Input Arguments
* `W::SparseMatrixCSC{Float64,Int}`: the edge weight matrix
* `method::Symbol`: the parition method to be used (:L or :Lrw or else ...)
* `v::Vector{Float64}`: the Fiedler vector supplied if `v` is a null vector

### Output Arguments
* `pm::Vector{Int}`: a vector of 1's and -1's
* `v::Vector{Float64}`: the Fiedler vector used for graph partitioning

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito)
Translated and revised by Naoki Saito, Feb. 9, 2017
"""
function partition_fiedler(W::SparseMatrixCSC{Float64,Int};
                           method::Symbol = :Lrw,
                           v::Vector{Float64} =Vector{Float64}(undef, 0))
#
# Easy case: the Fiedler vector is provided
#
    if !isempty(v)
        (pm, v) = partition_fiedler_pm(v)
        if size(W, 1) == length(v)
            # recovering the algebraic connectivity a(G)
            # Note: regardless of v coming from L or Lrw, the eigenvalue
            # computation below is fine: v'*L*v = a(G) or arw(G) because
            # for L: v'*L*v = \lambda v'*v = \lambda
            # for Lrw: vrw'*L*vrw = \lambda vrw'*D*vrw = \lambda
            # Note also, the result of the quadratic form is scalar
            # mathematically, but in julia, it's 1x1 array; so need to
            # convert it to a real scalar via `val = val[1]`.
            val = v' * (Diagonal(vec(sum(W, dims = 1))) - W) * v; val = val[1]
            pm = partition_fiedler_troubleshooting(pm, v, W, val)
        else
            pm = partition_fiedler_troubleshooting(pm)
        end
        return pm, v
    end

#
# Now, we need to compute the Fiedler vector here. Preliminaries first.
#
    N = size(W, 1)
    eigs_flag = 0

    sigma = eps()
    cutoff = 128 # this value could be changed...

    # handle the case when there are 2 nodes
    if N == 2
        pm = [1; -1]
        v = pm / sqrt(2)
        return pm, v
    end

    # Use L (unnormalized) if specified or edge weights are very small
    if method == :L || minimum(sum(W, dims = 1)) < 10^3 * sigma
        if N > cutoff # for a relative large W
            v0 = ones(N) / sqrt(N)
            try
                # MATLAB, [v,val,eigs_flag] = eigs(diag(sum(W))-W,2,sigma,opts);
                val, vtmp = eigs(sparse(Diagonal(vec(sum(W, dims = 1)))) - W,
                                 nev = 2, sigma = sigma, v0 = v0)
                # val = real.(val)
                # vtmp = real.(vtmp)
                if isa(val[1], Complex{Float64})
                    eigs_flag = 1
                end
            catch emsg
                @warn("Exception in eigs(L) occurred: ", emsg)
                eigs_flag = 2
            end
            if eigs_flag == 0   # no problem in `eigs` is detected.
                val, ind = findmax(val) # val is set to be a scalar.
                v = vec(vtmp[:, ind])   # This is the Fiedler vector!
                # In MATLAB, val was returned as a diagonal matrix, so we had to
                # do all these manipulations while julia, val is simply a vector.
                # [~,ind] = max(diag(val));
                # v = v(:,ind);
                # val = val(ind,ind);
            end
        end
        if N <= cutoff || eigs_flag != 0 # if W is small or eigs had a problem,
                                         # then use full svd.
            vtmp, val = svd(Matrix(Diagonal(vec(sum(W, dims = 1))) - W)) # v <-> U, val <-> S
            # Diagonal(vec(sum(W,1))) is Matrix{Float64} while W is SparseMatrix
            # But the results of the subtraction becomes Matrix{Float64}.
            # Also be careful here! val[end-2] == val[end-1] can happen.
            # If that is the case, vtmp[:,end-2] could also be the Fiedler vec.
            v = vec(vtmp[:, end - 1]) # v is set to be a vector.
            val = val[end - 1]       # val is set to a scalar.
            # In MATLAB, `full` was necessary to convert a sparse matrix format
            # to a usual matrix format as follows
            # [v,val,~] = svd(full(diag(sum(W))-W));
        end
    elseif method == :Lrw # Otherwise, use L_rw, which is normally preferred.
        if N > cutoff           # for a relatively large W
            v0 = ones(N) / sqrt(N)
            try
                # MATLAB: [v,val,eigs_flag] = ...
                #           eigs(diag(sum(W))-W,diag(sum(W)),2,sigma,opts);
                # This is L*v = \lambda*D*v case.
                temp = sparse(Diagonal(vec(sum(W, dims = 1))))
                val, vtmp = eigs(temp - W, temp,
                                 nev = 2, sigma = sigma, v0 = v0)
                # val = real.(val)
                # vtmp = real.(vtmp)
                if isa(val[1], Complex{Float64})
                    eigs_flag = 1
                end
            catch emsg
                @warn("Exception in eigs(Lrw) occurred: ", emsg)
                eigs_flag = 2
            end

            if eigs_flag == 0
                val, ind = findmax(val) # val is set to be a scalar.
                v = vtmp[:, ind]
                # MATLAB: v = (full(sum(W,2)).^(-0.5)) .* v
                # v = vec((sum(W, dims = 2).^(-0.5)) .* v) # This is the Fiedler vector!
            end
        end
        if N <= cutoff || eigs_flag != 0 # if W is small or eigs had a problem,
                                         # then use full svd
            colsumW = vec(sum(W, dims = 1))
            D = Diagonal(colsumW)
            D2 = Diagonal(colsumW.^(-0.5))
            vtmp, val = svd(Matrix(D2 * (D - W) * D2)) # SVD of Lsym
            # MATLAB: [v,val,~] = svd(full(...
            # bsxfun(@times,bsxfun(@times,full(sum(W,2)).^(-0.5),diag(sum(W))-W),
            #  full(sum(W,1)).^(-0.5)) ) );
            #
            # SVD is ordered in the order of nonincreasing singular values,
            # so, (end-1) entry corresponds to the algebraic connectivity.
            # Also be careful here! val[end-2] == val[end-1] can happen.
            # If that is the case, vtmp[:,end-2] could also be the Fiedler vec.
            val = val[end - 1]   # val is set to be a scalar.
            v = vtmp[:, end - 1] # v is set to be a vector.
            v = vec((sum(W, dims = 2).^(-0.5)) .* v) # This is the Fiedler vector of Lrw!
        end
    else                        # Unknown method is specified.
        error("Graph partitioning method :", method, " is not recognized!")
    end

    # partition using the Fiedler vector and troubleshoot any potential issues
    pm, v = partition_fiedler_pm(v)
    pm = partition_fiedler_troubleshooting(pm, v, W, val)
    return pm, v
end # of partition_fiedler

"""
    (pm, v) = partition_fiedler_pm(v::Vector{Float64})

    Partition an input Graph based on the Fiedler vector `v`

### Input Argument
* `v::Vector{Float64}`: the input Fiedler vector

### Output Arguments
* `pm::Vector{Int}`: the partition info
* `v::Vector{Float64}`: the output Fiedler vector if requested

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito)
Translated and revised by Naoki Saito, Feb. 10, 2017
"""
function partition_fiedler_pm(v::Vector{Float64})
    # define a tolerance: if abs(x) < tol, then we view x = 0.
    tol = 10^3 * eps()

    # set the first nonzero element to be positive
    row = 1
    while abs(v[row]) < tol && row < length(v)
        row += 1
    end
    if v[row] < 0
        v = -v
    end

    # assign each point to either region 1 or region -1, and assign any zero
    # entries to the smaller region
    if sum(v .>= tol) > sum(v .<= -tol) # more (+) than (-) entries
        pm = 2*(v .>= tol) .- 1
    else
        pm = 2*(v .<= -tol) .- 1
    end

    # make sure the first point is assigned to region 1 (not -1)
    if pm[1] < 0
        pm = -pm
    end
    return pm, v
end # of partition_fiedler_pm

"""
    pm = partition_fiedler_troubleshooting(pm,v,W,val)

Troubleshoot potential issues with the partitioning

### Input Arguments
* `pm::Vector{Int}`: an input partition info (+1 or -1) of each node
* `v::Vector{Float64}`: the Fiedler vector of the graph
* `W::SparseMatrixCSC{Float64,Int}`: an input edge weight matrix
* `val::Float64`: the algebraic connectivity (i.e., ``\\lambda_2``)

### Ouput Arguments
* `pm::Vector{Int}`: a final partition info vector

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito)
Translated and revised by Naoki Saito, Feb. 10, 2017
"""
function partition_fiedler_troubleshooting(pm::Vector{Int},v::Vector{Float64},
                                           W::SparseMatrixCSC{Float64,Int},
                                           val::Float64)
    # Initial setup
    N = length(pm)
    tol = 10^3 * eps()

    # If pm vector indicates no partition (single sign) or ambiguities via
    # 0 entries, then we'll trouble shoot.
    if sum(pm .< 0) == 0 || sum(pm .> 0) == 0 || sum(abs.(pm)) < N
        # Case 1: an input graph is not connected (i.e., alg. conn < tol)
        if val < tol
            pm = 2 .* (abs.(v) .> tol) .- 1
            while sum(pm .< 0) == 0 || sum(pm .> 0) == 0
                tol = 10 * tol    # relax the tolerance
                pm = 2 .* (abs.(v) .> tol) .- 1
                if tol > 1
                    pm[1:Int(ceil(N / 2))] .= 1
                    pm[(Int(ceil(N / 2)) + 1):N] .= -1
                end
            end
        # Case 2: it is connected, but something funny happened
        else
            pm = 2 .* (v .>= mean(v)) .- 1
            if sum(abs.(pm)) < N
                # assign the near-zero points based on the values of v
                # at their neighbor nodes
                pm0 = find(pm .== 0)
                pm[pm0] = (W[pm0, :] * v .> tol) - (W[pm0, :] * v .< -tol)
                # assign any remaining zeros to the group with fewer members
                pm[find(pm .== 0)] = (sum(pm .> 0) - sum(pm .< 0)) .>= 0
            end
        end
        # if one region has no points after all the above processing
        if sum(pm .< 0) == 0 || sum(pm .> 0) == 0
            pm[1:Int(ceil(N / 2))] .= 1
            pm[Int(ceil(N / 2) + 1):N] .= -1
        end
    end

# make sure that the first point is assigned as a 1
    if pm[1] < 0
        pm = -pm
    end
    return pm
end # of partition_fiedler_troubleshooting

"""
    pm = partition_fiedler_troubleshooting(pm)

Troubleshoot potential issues with the partitioning

### Input Arguments
* `pm::Vector{Int}`: an input partition info (+1 or -1) of each node

### Ouput Arguments
* `pm::Vector{Int}`: a final partition info vector

Copyright 2015 The Regents of the University of California

Implemented by Jeff Irion (Adviser: Dr. Naoki Saito)
Translated and revised by Naoki Saito, Feb. 10, 2017
"""
function partition_fiedler_troubleshooting(pm::Vector{Int})
    N = length(pm)
    if sum(pm .< 0) == 0 || sum(pm .> 0) == 0
        pm[1:Int(ceil(N / 2))] = 1
        pm[Int(ceil(N / 2) + 1):N] = -1
    end
    # make sure that the first point is assigned as a 1
    if pm[1] < 0
        pm = -pm
    end
    return pm
end


text:
using simplexpop
using Test

# Make these the same order as their includes in simplexpop.jl.
@testset "sort.jl" begin include("test_sort.jl") end
@testset "hilbert.jl" begin include("test_hilbert.jl") end
@testset "assign_weights.jl" begin include("test_assign_weights.jl") end

ext: jl


text:
using Base.Test

@testset "No common tests" begin
    nothing
end


text:
## Max length of array is 2^16
## A sharder is (currenty) (subject to change due to unhandled reflections) a function that takes a RayResults struct and returns a RGB value
# DefaultShaderArray[1] = background
const DefaultShaderArray = [
    (x::RayResults)->RGB(0.35,0.35,0.35),
    #(x::RayResults)->(RGB(0.8,0.8,0.8)*10*SignedDistanceFields.clamp01(1/norm(x.Position))),
    (x::RayResults)->(RGB(x.Grad[1]/2+0.5,x.Grad[2]/2+0.5,x.Grad[3]/2+0.5)),
    (x::RayResults)->RGB(0.22,0.545,0.133)*(abs(x.Position[1]/(5*sqrt(2)))),
    (x::RayResults)->RGB(0.584,0.345,0.689)*(abs(x.Position[1]/(5*sqrt(2)))),
    (x::RayResults)->RGB(0.796,0.235,0.2)*(abs(x.Position[1]/(5*sqrt(2)))),
    (x::RayResults)->RGB(0.35,0.35,0.55)]

## CSGBuilding

function render!(frame, node::CSGNode, cam, x, y, shaders)
    for i=1:x ## loop able to produce more than 1 gigapixel per hour per core on my laptop if fed correctly
        for j=1:y
            a = los(cam, i, j)
            temp = raymarch(a[1], a[2], node, 1000.0)
            frame[i,j] = shaders[temp.Shader](temp)
        end
    end
    frame
end

function render(node::CSGNode, sizet, camf::Function, shaders)
    @unpack x, y = sizet
    #Höhe Breite - x
    #Heigth Width - y
    cam = camf(x, y)
    #frame = Array{RayResults}(undef,x,y);
    myPic = Array{RGB}(undef,x,y);

    render!(myPic, node, cam, x, y, shaders)
    return myPic
end

## SignedDistanceFields

function render!(frame, scene, cam, x, y, shaders)
    for i=1:x ## loop able to produce more than 1 gigapixel per hour per core on my laptop if fed correctly
        for j=1:y
            a = los(cam, i, j)
            #frame[i,j]
            temp = raymarch(a[1], a[2], scene, 1000.0)
            ## upto 3x over head due to excessive copying
            ## Stencils on this ^ would be fun :)
            frame[i,j] = shaders[temp.Shader](temp)
            ##TODO:reflections have no place in this model
            ## some improvement possible with the shader
        end
    end
    frame
end

function render(scene, sizet, camf::Function, shaders)
    @unpack x, y = sizet
    #Höhe Breite - x
    #Heigth Width - y
    cam = camf(x, y)
    frame = Array{RGB}(undef,x,y);
    ## render
    render!(frame, scene, cam, x, y, shaders)
    return frame
end


text:
# This is an example of two synaptically-coupled single-compartment neurons with Hodgkin-Huxley dynamics.
# It's written for illustrative purposes and uses the high-level interface of ModelingToolkit.jl as it already exists today.

using DifferentialEquations, ModelingToolkit, Plots, OrdinaryDiffEq

# Channel gating functions
αm(V) = (0.01*(V + 55))/(1 - exp(-(V + 55)/10))
βn(V) = 0.125 * exp(-(V + 65)/80)
αm(V) = (0.1*(V + 40))/(1 - exp(-(V + 40)/10))
βm(V) = 4*exp(-(V + 65)/18)
αh(V) = 0.07*exp(-(V+65)/20)
βh(V) = 1/(1 + exp(-(V + 35)/10))

# Synaptic transmitter release where
# Tmax = 1 mM Vthreshold = 1 mV Kp = 1 mV
# T(Vpre) = 1 / (1 + exp(-(Vpre - 2)/5))

# Synaptic transmitter release via heaviside step at Vthreshold
heaviside(Vpre) = ifelse(Vpre > 1.0, 1.0, 0.0)

# Steady-state channel gating
m∞(V) = αm(V)/(αm(V) + βm(V))
h∞(V) = αh(V)/(αh(V) + βh(V))
n∞(V) = αn(V)/(αn(V) + βn(V))

@parameters Ie[1:2] t gl gna gk gs El Ena Ek Es Cm ar ad
@variables V[1:2](t) m[1:2](t) h[1:2](t) n[1:2](t) s[1:2](t)
@derivatives D'~t
@register heaviside(V)
# HH equations
revs = reverse(s)
eqs = [(@. D(V) ~ (Ie - gl*(V - El) - gna*m^3*h*(V - Ena) - gk*n^4*(V - Ek) - gs*revs*(V - Es))/Cm)...,
       (@. D(m) ~ αm(V) * (1 - m) - βm(V)*m)...,
       (@. D(h) ~ αh(V) * (1 - h) - βh(V)*h)...,
       (@. D(n) ~ αn(V) * (1 - n) - βn(V)*n)...,
       (@. D(s) ~ ar*heaviside(V)*(1-s) - ad*s)...]

sys = ODESystem(eqs, t, [V..., m..., h..., n..., s...], [Ie..., gl, gna, gk, gs, El, Ena, Ek, Es, Cm, ar, ad])

Vinit = -65.0
radius = 0.0025 # given in cm
area = 4*pi*(radius)^2 # area in cm²

# Set applied current to zero initially
u0 = [(V .=> Vinit)...
      (m .=> m∞(Vinit))...
      (h .=> h∞(Vinit))...
      (n .=> n∞(Vinit))...
      (s .=> 0.0)...]

p = [Ie[1] => 0.8e-3 # tonic excitability to force first neuron to fire repetitively
     Ie[2] => 0.0
     gl  => 0.3*area
     gna => 120.0*area
     gk  => 36.0*area
     gs  => 0.025*area
     El  => -54.4
     Ena => 50.0
     Ek  => -77.0
     Es  => 0.0
     Cm  => 1.0e-3*area
     ar  => 1.1
     ad  => 0.19]

tspan = (0.0, 250.0)
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, KenCarp47())
plot(sol, vars=(0,1))
plot!(sol, vars=(0,2))

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation


text:
@with_kw mutable struct CrosswalkDriver <: DriverModel{LonAccelDirection}
    a::LonAccelDirection = LonAccelDirection(0., 1)
    navigator::RouteFollowingIDM = RouteFollowingIDM()
    crosswalk::Lane = Lane()
    conflict_lanes::Vector{Lane{Float64}} = Lane{Float64}[]
    intersection_entrances::Vector{Lane{Float64}} = Lane{Float64}[]
    ped_model::ConstantPedestrian = ConstantPedestrian()
    ped_start::Float64 = 4.0 # assume a 5meter buffer before a pedestrian reaches the road
    stop_delta::Float64 = 0.7
    accel_tol::Float64 = 0.1
    d_tol::Float64 = 0.5

    #states
    yield::Bool = true
    priority::Bool = false
    stop::Bool = false
    wait_list::Vector{Int64} = Int64[]

    # transition
    clear::Bool = false

    debug::Bool = false
end

function AutomotiveDrivingModels.reset_hidden_state!(model::CrosswalkDriver)
    model.priority = false
    model.stop = false
    model.clear = false
    model
end

AutomotiveDrivingModels.get_name(::CrosswalkDriver) = "CrosswalkDriver"
Base.rand(model::CrosswalkDriver) = model.a

function AutomotiveDrivingModels.observe!(model::CrosswalkDriver, scene::Scene, roadway::Roadway, egoid::Int)
    ego = scene[findfirst(egoid, scene)]
    AutomotiveDrivingModels.observe!(model.navigator, scene, roadway, egoid) # set the direction
    dir = model.navigator.dir
    a_lon =0.
    a_lon_idm = model.navigator.a
    dist_to_cw = get_distance_to_crosswalk(model, ego, roadway, -model.stop_delta)
    # if !model.yield
    #     a_lon = a_lon_idm
    #     # println("veh ", egoid, " not yielding to crosswalk ", model.crosswalk.tag)
    # else
    #     if !model.priority && !model.stop # reach stop line
    #         a_lon = min(a_lon_idm, stop_at_dist(model, ego, dist_to_cw))
    #     elseif !model.priority && model.stop # wait
    #         a_lon = -model.navigator.d_max # negative to make sure v stays 0
    #     else
    #         a_lon = a_lon_idm # just idm
    #     end
    # end
    # if !model.stop
    #     update_stop!(model, ego, roadway, dist_to_cw)
    # end
    # grow_wait_list!(model, scene, roadway, egoid)
    # ungrow_wait_list!(model, scene, roadway, egoid)
    # update_priority!(model, scene, roadway, egoid)

    if !model.yield
        a_lon = a_lon_idm
    else
        grow_wait_list!(model, scene, roadway, egoid)
        ungrow_wait_list!(model, scene, roadway, egoid)
        update_priority!(model, scene, roadway, egoid)
        if !model.priority
            a_lon = min(a_lon_idm, AutomotivePOMDPs.stop_at_dist(model, ego, dist_to_cw))
            # println("dist to cw $dist_to_cw, crosswalk $(model.crosswalk.tag), $(get_lane(roadway, car) ∈ model.intersection_entrances)")
            # println("priority $(model.priority), crossing $(is_crossing(ped, model.crosswalk, model.conflict_lanes, roadway)), a_lon $a_lon, v $(car.v)")
        else
            a_lon = a_lon_idm # just idm
        end
    end

    if model.debug
        println("ID ", egoid, " priority ", model.priority, " stop ", model.stop, " a ", a_lon, " wait list ", model.wait_list, " dist_to_cw ", dist_to_cw)
    end
    model.a = LonAccelDirection(a_lon, dir)
end

"""
Check if all the pedestrian have crossed
"""
function update_priority!(model::CrosswalkDriver, scene::Scene, roadway::Roadway, egoid::Int)
    ego = scene[findfirst(egoid, scene)]
    model.priority = has_passed(model, ego, roadway) || isempty(model.wait_list)
end

function has_passed(model::CrosswalkDriver, ego::Vehicle, roadway::Roadway)
    cw_length = get_end(model.crosswalk)
    # cw_width = model.crosswalk.width
    cw_center = posg(Frenet(model.crosswalk, cw_length/2, 0.), roadway)
    cw_to_car = get_front(ego) - cw_center
    car_vec = get_front(ego) - ego.state.posG
    has_passed = dot(cw_to_car, car_vec) > 0.
    return has_passed
end

function has_passed(model::CrosswalkDriver, ego::VehicleState, roadway::Roadway)
    cw_length = get_end(model.crosswalk)
    cw_center = posg(Frenet(model.crosswalk, cw_length/2), roadway)
    cw_to_car = ego.posG - cw_center
    car_vec = polar(1., ego.posG.θ) - ego.posG + ego.posG
    has_passed = dot(cw_to_car, car_vec) > 0.
    return has_passed
end

function grow_wait_list!(model::CrosswalkDriver, scene::Scene, roadway::Roadway, egoid::Int)
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN &&
           is_crossing(veh, model.crosswalk, model.conflict_lanes, scene, roadway) &&
           !(veh.id ∈ model.wait_list)
            push!(model.wait_list, veh.id)
        end
    end
end

function ungrow_wait_list!(model::CrosswalkDriver, scene::Scene, roadway::Roadway, egoid::Int)
    n = length(model.wait_list)
    to_remove = Int64[]
    for i=1:n
        pedid = model.wait_list[i]
        ped_ind = findfirst(pedid, scene)
        if ped_ind == nothing
            push!(to_remove, i)
            continue
        end
        ped = scene[ped_ind]
        if !is_crossing(ped, model.crosswalk, model.conflict_lanes, scene, roadway)
            push!(to_remove, i)
        end
    end
    deleteat!(model.wait_list, to_remove)
end

function is_crossing(ped::Vehicle, crosswalk::Lane{T}, conflict_lanes::Vector{Lane{T}}, scene::Scene, roadway::Roadway) where T
    # check if the pedestrian is in the conflict zone
    ped_lane = get_lane(roadway, ped)
    if ped_lane.tag != crosswalk.tag
        return false
    end
    for lane in conflict_lanes
        ped_f = Frenet(proj(ped.state.posG, lane, roadway, move_along_curves=false), roadway)
        if abs(ped_f.t) < lane.width/2 && get_lane(roadway, ped).tag == crosswalk.tag
            return true
        end
    end
    # at this point, the pedestrian is not on the road
    # check if the pedestrian is going to cross or not
    if AutomotivePOMDPs.direction_from_center(ped, crosswalk) > 0. && get_lane(roadway, ped).tag == crosswalk.tag && is_free(crosswalk, scene, roadway) #&& !isapprox(ped.state.v, 0.)
        return true
    end
    return false
end

function is_free(crosswalk::Lane, scene::Scene, roadway::Roadway)
    #check that no vehicle is stopped on the crosswalk
    for veh in scene
        if veh.def.class == AgentClass.PEDESTRIAN
            continue
        end
        posG = veh.state.posG
        cw_posF = Frenet(proj(posG, crosswalk, roadway,move_along_curves=false), roadway)
        lane = get_lane(roadway, veh)
        t_car = cw_posF.t
        if abs(cw_posF.t) < crosswalk.width/2 + veh.def.length/2*sin(cw_posF.ϕ)
            return false
        end
    end
    return true
end

function get_distance_to_crosswalk(model::CrosswalkDriver, veh::VehicleState, roadway::Roadway, delta::Float64 = 0.)
    lane = get_lane(roadway, veh)
    Δs = 0.
    # if lane ∈ model.conflict_lanes
    #     cw_length = get_end(model.crosswalk)
    #     cw_center = posg(Frenet(model.crosswalk, cw_length/2), roadway)
    #     collision_point = VecSE2(cw_center.x, veh.posG.y)
    #     cw_proj = Frenet(collision_point, lane, roadway)
    #     Δs = cw_proj.s - veh.posF.s + delta - model.crosswalk.width
    # elseif lane ∈ model.intersection_entrances # stop before the intersection to not block traffic
    #     Δs = get_end(lane) - veh.posF.s + delta - model.crosswalk.width
    # end
    if lane ∈ model.intersection_entrances # stop before the intersection to not block traffic
        Δs = get_end(lane) - veh.posF.s + delta - model.crosswalk.width
    elseif lane ∈ model.conflict_lanes
        cw_length = get_end(model.crosswalk)
        cw_start = model.crosswalk.curve[1].pos
        cw_center = posg(Frenet(model.crosswalk, cw_length/2), roadway)
        cw_to_car = veh.posG - cw_center
        car_vec = polar(1., veh.posG.θ) - veh.posG + veh.posG
        d = dot(car_vec, cw_to_car)
        Δs = d + delta - model.crosswalk.width/2
    end
    return Δs
end

function get_distance_to_crosswalk(model::CrosswalkDriver, veh::Vehicle, roadway::Roadway, delta::Float64 = 0.)
    return get_distance_to_crosswalk(model, veh.state, roadway, delta)
end

function update_stop!(model::CrosswalkDriver, veh::Vehicle, roadway::Roadway, dist_to_cw::Float64)
    update_stop!(model, veh.state, roadway, dist_to_cw)
end

function update_stop!(model::CrosswalkDriver, veh::VehicleState, roadway::Roadway, dist_to_cw::Float64)
    if veh.v ≈ 0. && isapprox(dist_to_cw, 0, atol=1.)
        model.stop = true
    end
end

# # croswalk driver
# @with_kw mutable struct CrosswalkDriver <: DriverModel{LonAccelDirection}
#     a::LonAccelDirection = LonAccelDirection(0., 1)
#     navigator::RouteFollowingIDM = RouteFollowingIDM()
#     crosswalk::Lane = Lane()
#     intersect::Lane = Lane()
#     ped_model::ConstantPedestrian = ConstantPedestrian()
#     ped_start::Float64 = 4.0 # assume a 5meter buffer before a pedestrian reaches the road
#     stop_delta::Float64 = 0.0
#     accel_tol::Float64 = 0.1
#     d_tol::Float64 = 0.5
#
#     #states
#     go::Bool = false
#     stop::Bool = false
#
#     # transition
#     clear::Bool = false
#
#     debug::Bool = false
# end
#
# AutomotiveDrivingModels.get_name(::CrosswalkDriver) = "CrosswalkDriver"
# Base.rand(model::CrosswalkDriver) = model.a
#
# function AutomotiveDrivingModels.reset_hidden_state!(model::CrosswalkDriver)
#     model.go = false
#     model.stop = false
#     model
# end
#
# function AutomotiveDrivingModels.observe!(model::CrosswalkDriver, scene::Scene, roadway::Roadway, egoid::Int)
#     ego = scene[findfirst(egoid, scene)]
#     AutomotiveDrivingModels.observe!(model.navigator, scene, roadway, egoid) # set the direction
#     dir = model.navigator.dir
#     a_lon =0.
#     a_lon_idm = model.navigator.a
#     if model.go
#         a_lon = a_lon_idm
#     else
#         update_clear!(model, ego, scene, roadway)
#         if model.clear
#             model.go = true
#         else # brake!
#             model.clear = false
#             d = dist_to_stop(model, ego, roadway)
#             if model.debug
#                 println("braking! d = ", d, "a = ", stop_at_dist(model, ego, d))
#             end
#             a_lon = min(a_lon_idm, stop_at_dist(model, ego, d))
#             if ego.state.v ≈ 0. && d < model.d_tol
#                 a_lon = 0.
#             end
#         end
#     end
#     if model.debug
#         println("veh ", egoid, " go ", model.go, " stop ", model.stop, " clear ", model.clear)
#     end
#     model.a = LonAccelDirection(a_lon, dir)
# end
#
#
# function update_clear!(model::CrosswalkDriver, ego::Vehicle, scene::Scene, roadway::Roadway)
#     # check how many pedestrians are present
#     peds = Vehicle[]
#     for veh in scene
#         if veh.def.class == AgentClass.PEDESTRIAN && get_lane(roadway, veh) == model.crosswalk
#             push!(peds, veh)
#         end
#     end
#     cw_length = get_end(model.crosswalk)
#     cw_center = posg(Frenet(model.crosswalk, cw_length/2), roadway)
#     collision_point = VecSE2(cw_center.x, ego.state.posG.y)
#     # set to the maximum value
#     buffer_ped = VehicleState(Frenet(model.crosswalk, model.ped_start), roadway,
#                               model.ped_model.v_desired)
#     # println("compute time ped")
#     min_ped_ttc = compute_time_of_approach(model, collision_point, buffer_ped, roadway, -model.intersect.width/2)
#     if model.debug
#         println("buffer ped ttc ", min_ped_ttc)
#     end
#     ego_ttc = compute_time_of_approach(model, collision_point, ego.state, roadway,
#                                        +model.crosswalk.width/2 + ego.def.length/2)
#     if ego.state.v ≈ 0.
#         ego_ttc = sqrt((model.crosswalk.width + ego.def.length)/model.navigator.a_max/2) # it is stopped
#     end
#     if model.debug
#         println("ego ttc ", ego_ttc)
#     end
#     has_crossed = -0.5*model.intersect.width/model.ped_model.v_desired
#     if model.debug
#         println("has_crossed ttc", has_crossed)
#     end
#     for ped in peds
#         ttc = compute_time_of_approach(model, collision_point, ped.state, roadway, -model.intersect.width/2)
#         if model.debug
#             println("ped $(ped.id) ttc ", ttc)
#         end
#         if has_crossed < ttc < min_ped_ttc
#             min_ped_ttc = ttc
#         end
#     end
#     if model.debug
#         println("min ped ttc ", min_ped_ttc)
#     end
#     if has_crossed < min_ped_ttc < ego_ttc
#         model.clear = false
#     else
#         model.clear = true
#     end
# end
#
# """
# compute the time to reach the projection of cw on the ego lane
# """
# function compute_time_of_approach(model::CrosswalkDriver, cw::VecSE2, ego::VehicleState, roadway::Roadway, delta::Float64=0.)
#     lane = get_lane(roadway, ego)
#     cw_proj = Frenet(cw, lane, roadway)
#     Δs = cw_proj.s - ego.posF.s + delta
#     # println("lane ", lane)
#     # println("cw_proj ", cw_proj.s)
#     # println("delta ", delta)
#     # println("ego posF", ego.posF.s)
#     # println("ds ", Δs)
#     v = max(ego.v, model.ped_model.v_desired)
#     time_of_approach = Δs/v
# end
#
# function dist_to_stop(model::CrosswalkDriver, ego::Vehicle, roadway::Roadway)
#     cw_length = get_end(model.crosswalk)
#     cw_center = posg(Frenet(model.crosswalk, cw_length/2), roadway)
#     collision_point = VecSE2(cw_center.x, ego.state.posG.y)
#     intersect_pt = Frenet(collision_point, model.intersect, roadway)
#     return intersect_pt.s - model.crosswalk.width/2 - ego.state.posF.s - ego.def.length/2 - model.stop_delta
# end
#
# function is_crossing(ped::Vehicle, crosswalk::Lane, roadway::Roadway)
# end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
@inline getdiagvalue(A::TronSparseMatrixCSC, i) = A.diag_vals[i]
@inline getdiagvalue(A::TronDenseMatrix, i) = A.vals[i,i]

function dssyax(n, A::TronDenseMatrix, x, y)
    zero = 0.0

    @inbounds for i=1:n
        y[i] = A.vals[i,i]*x[i]
    end

    @inbounds for j=1:n
        rowsum = zero
        for i=j+1:n
            rowsum += A.vals[i,j]*x[i]
            y[i] += A.vals[i,j]*x[j]
        end
        y[j] += rowsum
    end

    return
end

"""
Subroutine dssyax

This subroutine computes the matrix-vector product y = A*x,
where A is a symmetric matrix with the strict lower triangular
part in compressed column storage.
"""
dssyax(A::TronSparseMatrixCSC, x, y) = dssyax(A.n, A, x, y)
function dssyax(n, A::TronSparseMatrixCSC, x, y)
    zero = 0.0

    @inbounds for i in 1:n
        y[i] = A.diag_vals[i]*x[i]
    end

    @inbounds for j in 1:n
        rowsum = zero
        for i in A.colptr[j]:A.colptr[j+1]-1
            row = A.rowval[i]
            rowsum += A.tril_vals[i]*x[row]
            @inbounds y[row] += A.tril_vals[i]*x[j]
        end
        y[j] += rowsum
    end

    return
end

function reorder!(B::TronDenseMatrix, A::TronDenseMatrix, indfree, nfree, iwa)
    nnz = 0
    @inbounds for j=1:nfree
        jfree = indfree[j]
        B.vals[j,j] = A.vals[jfree,jfree]
        for i=jfree+1:A.n
            if iwa[i] > 0
                nnz += 1
                B.vals[iwa[i],j] = A.vals[i,jfree]
            end
        end
    end
    B.n = nfree
    return nnz
end

# Update matrix B inplace
function reorder!(B::TronSparseMatrixCSC, A::TronSparseMatrixCSC, indfree, nfree, iwa)
    B.colptr[1] = 1
    nnz = 0
    @inbounds for j=1:nfree
        jfree = indfree[j]
        B.diag_vals[j] = A.diag_vals[jfree]
        for ip = A.colptr[jfree]:A.colptr[jfree+1]-1
            if iwa[A.rowval[ip]] > 0
                nnz = nnz + 1
                B.rowval[nnz] = iwa[A.rowval[ip]]
                B.tril_vals[nnz] = A.tril_vals[ip]
            end
        end
        B.colptr[j+1] = nnz + 1
    end
    return nnz
end


text:
using MultivariateStats
using LinearAlgebra
import Statistics: mean
using Test

@testset "MDS" begin

    d = 3
    n = 10
    X0 = randn(d, n)
    G0 = X0'X0
    D0 = MultivariateStats.pairwise((x,y)->norm(x-y), eachcol(X0), symmetric=true)

    ## conversion between dmat and gram

    @assert issymmetric(D0)
    @assert issymmetric(G0)

    D = gram2dmat(G0)
    @test issymmetric(D)
    @test D ≈ D0

    G = dmat2gram(D0)
    @test issymmetric(G)
    @test gram2dmat(G) ≈ D0

    ## classical MDS
    @test_throws UndefKeywordError fit(MDS, X0, maxoutdim=3)   # disambiguate distance matrix from data matrix
    M = fit(MDS, X0, maxoutdim=3, distances=false)
    @test indim(M) == d
    @test outdim(M) == 3
    @test size(projection(M)) == (n,3)
    @test length(eigvals(M)) == 3
    @test stress(M) ≈ 0.0 atol = 1e-10

    X = transform(M)
    @test size(X) == (3,n)
    @test MultivariateStats.pairwise((x,y)->norm(x-y), eachcol(X0), symmetric=true) ≈ D0

    @test_throws DimensionMismatch transform(M, rand(d+1))
    y = transform(M, X0[:, 1])
    @test X[:, 1] ≈ y

    # use only distance matrix
    M = fit(MDS, D0, maxoutdim=3, distances=true)
    @test isnan(indim(M))
    @test outdim(M) == 3
    @test stress(M) ≈ 0.0 atol = 1e-10

    X = transform(M)
    @test size(X) == (3,n)
    @test MultivariateStats.pairwise((x,y)->norm(x-y), eachcol(X0), symmetric=true) ≈ D0

    @test_throws AssertionError transform(M, X0[:, 1])
    @test_throws DimensionMismatch transform(M, rand(d+1); distances = true)
    d = MultivariateStats.pairwise((x,y)->norm(x-y), X0, X0[:,2]) #|> vec
    y = transform(M, d, distances=true)
    @test X[:, 2] ≈ y

    #Test MDS embeddings in dimensions >= number of points
    M = fit(MDS, [0. 1.; 1. 0.], maxoutdim=2, distances=true)
    @test outdim(M) == 2
    @test transform(M) == [-0.5 0.5; 0 0]

    M = fit(MDS, [0. 1.; 1. 0.], maxoutdim=3, distances=true)
    @test outdim(M) == 3
    @test transform(M) == [-0.5 0.5; 0 0; 0 0]

    #10 - dmat2gram produces negative definite matrix
    let D = [1.0 0.5 0.3181818181818182 0.38095238095238093 0.6111111111111112 0.36363636363636365 0.3333333333333333 0.4444444444444444 0.45454545454545453 0.38095238095238093
            0.5 1.0 0.3333333333333333 0.4166666666666667 0.32 0.37037037037037035 0.35 0.38461538461538464 0.38461538461538464 0.4
            0.3181818181818182 0.3333333333333333 1.0 0.38095238095238093 0.125 0.22727272727272727 0.47058823529411764 0.35294117647058826 0.2608695652173913 0.21052631578947367
            0.38095238095238093 0.4166666666666667 0.38095238095238093 1.0 0.17391304347826086 0.2962962962962963 0.38095238095238093 0.45454545454545453 0.3076923076923077 0.5454545454545454
            0.6111111111111112 0.32 0.125 0.17391304347826086 1.0 0.5 0.15 0.391304347826087 0.38095238095238093 0.3333333333333333
            0.36363636363636365 0.37037037037037035 0.22727272727272727 0.2962962962962963 0.5 1.0 0.2857142857142857 0.44 0.4 0.3181818181818182
            0.3333333333333333 0.35 0.47058823529411764 0.38095238095238093 0.15 0.2857142857142857 1.0 0.42105263157894735 0.45 0.23529411764705882
            0.4444444444444444 0.38461538461538464 0.35294117647058826 0.45454545454545453 0.391304347826087 0.44 0.42105263157894735 1.0 0.5416666666666666 0.5555555555555556
            0.45454545454545453 0.38461538461538464 0.2608695652173913 0.3076923076923077 0.38095238095238093 0.4 0.45 0.5416666666666666 1.0 0.4
            0.38095238095238093 0.4 0.21052631578947367 0.5454545454545454 0.3333333333333333 0.3181818181818182 0.23529411764705882 0.5555555555555556 0.4 1.0],
        Xt = [-0.27529104101488666 0.006134513718202863 0.33298809606740326 0.2608994458893664 -0.46185275796909575 -0.23734315039370618 0.29972782027671513 0.03827901455843394 -0.04096713097883363 0.07742518984640051
            -0.08177061420820278 -0.0044504235228030225 -0.3271919093638943 0.28206254638779243 -0.0954706915166714 -0.07137742126520012 -0.30754933764853587 0.18582658369448027 -0.03715307349750036 0.45707434094053534]',
        M = transform(fit(MDS, D, maxoutdim=2, distances=true))'
        @test M ≈ Xt .* cis.(angle.(sum(conj.(Xt) .* M, dims=1)))
    end

    #10 - test degenerate problem
    @test transform(fit(MDS, zeros(10, 10), maxoutdim=3, distances=false)) == zeros(3, 10)

    # out-of-sample
    D = [0 1 2 1;
         1 0 1 2;
         2 1 0 1;
         1 2 1 0.0f32]

    M = fit(MDS, sqrt.(D), maxoutdim=2, distances=true)
    X = transform(M)
    @test D ≈ MultivariateStats.pairwise((x,y)->sum(abs2, x-y), eachcol(X), symmetric=true)
    @test eltype(X) == Float32

    a = Float32[0.5, 0.5, 0.5, 0.5]
    A = vcat(hcat(D, a), hcat(a', zeros(Float32, 1, 1)))
    M⁺ = fit(MDS, sqrt.(A), maxoutdim=2, distances=true)
    X⁺ = transform(M⁺)
    @test A ≈ MultivariateStats.pairwise((x,y)->sum(abs2, x-y), eachcol(X⁺), symmetric=true)

    y = transform(M, a, distances=true)
    Y = [X y]
    @test A ≈ MultivariateStats.pairwise((x,y)->sum(abs2, x-y), eachcol(Y), symmetric=true)
    @test eltype(Y) == Float32

    # different input types
    d = 3
    X = randn(Float64, d, 10)
    XX = convert.(Float32, X)

    y = randn(Float64, d)
    yy = convert.(Float32, y)

    M = fit(MDS, X, maxoutdim=3, distances=false)
    MM = fit(MDS, XX, maxoutdim=3, distances=false)

    # test that mixing types doesn't error
    transform(M, yy)
    transform(MM, y)

    # type stability
    for func in (projection, eigvals, stress)
        @test eltype(func(M)) == Float64
        @test eltype(func(MM)) == Float32
    end
end


text:
using Conda
using PyCall
# Conda.add("pcraster"; channel="conda-forge")

pcr = pyimport("pcraster")
np = pyimport("numpy")

function spread_pcr(points, initial, friction; res=2, nodata=NaN)
    pcr.setclone(size(points)..., res, 0, 0)
    ppoints = pcr.numpy2pcr(pcr.Nominal, points, nodata)
    @time presult = pcr.spread(ppoints, initial, friction)
    transpose(pcr.pcr2numpy(presult, nodata))
end


text:
include("electrical_components.jl")

@named source = ConstantVoltage(V=10.0)
@named resistor = Resistor(R=1.0)
@named inductor1 = Inductor(L=1.0e-2)
@named inductor2 = Inductor(L=2.0e-2)
@named ground = Ground()

eqs = [
       connect(source.p, resistor.p)
       connect(resistor.n, inductor1.p)
       connect(inductor1.n, inductor2.p)
       connect(source.n, inductor2.n)
       connect(inductor2.n, ground.g)
      ]

@named ll_model = ODESystem(eqs, t)
ll_model = compose(ll_model, [source, resistor, inductor1, inductor2, ground])

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
function construct_service!(canonical_model::CanonicalModel,
                            service::Type{SD},
                            service_formulation::Type{SV},
                            system_formulation::Type{S},
                            sys::PSY.System;
                            kwargs...) where {SD<:PSY.Service,
                                              SV<:AbstractServiceFormulation,
                                              S<: PM.AbstractPowerModel}

    return

end

#=

##################
This code still need to be rewritten for the new infrastructure in PowerSimulations
##################

function get_devices(sys::PSY.System, device::Type{PSY.ThermalGen})
    return sys.generators.thermal
end
function get_devices(sys::PSY.System, device::Type{PSY.RenewableGen})
    return sys.generators.renewable
end
function get_devices(sys::PSY.System, device::Type{PSY.HydroGen})
    return sys.generators.hydro
end
function get_devices(sys::PSY.System, device::Type{PSY.PSY.ElectricLoad})
    return sys.loads
end

function construct_service!(m::JuMP.AbstractModel, service::PSY.StaticReserve, device_formulation::Type{RampLimitedReserve}, devices::Array{NamedTuple{(:device, :formulation), Tuple{DataType, DataType}}}, sys::PSY.System; kwargs...)

    dev_set = Array{NamedTuple{(:device, :formulation), Tuple{PSY.Device, DataType}}}([])

    for device in devices
        if device != nothing
            D = get_devices(sys, device.device)
            for d in D
                if d in service.contributingdevices
                    push!(dev_set, (device=d, formulation=device.formulation))
                end
            end
        end
    end

    if !isempty(dev_set)

        p_rsv =  reservevariables(m, dev_set, sys.time_periods)
        m =  reserves(m, dev_set, service, sys.time_periods)

    end

    return m

end
=#


text:
function laplacian4!(v, u, dx)
    v[begin] = (u[end] + u[begin+1] - 2u[begin])/dx^2
    @. @views v[begin+1:end-1] = (u[begin:end-2] + u[begin+2:end] - 2u[begin+1:end-1])/dx^2
    v[end] = (u[end-1] + u[begin] - 2u[1end])/dx^2
    return
end

function heateq4(u0, dx, tmax, N=200)
    t = 0:dx:tmax
    dt = step(t)
    u = similar(u0, length(u0), length(t)+1)
    u[:, 1] = u0
    v = similar(u0)
    N = 100
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian4!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

using StaticArrays

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
m = 200
noise1, noise2 = 0.3randn(m), 0.3randn(m)
f(x) = SVector{m}((1 .+ noise1)*sin(x) .- (1 .+ noise2)*cos(2x))
u0 = f.(x);

tmax = 1.0
t, u = heateq4(u0, dx, tmax);

write("sol_t.txt", repr(t))
write("sol_u.txt", repr(u))

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
    SysLED(i::Int32)
Type representing the system LEDs on the BeagleBone. The LEDs are indexed by
i ∈ [1,2,3,4].
"""
type SysLED <: IO_Object
    i::Int32
    basedir::String
    filestream::IOStream
    function SysLED(i::Int32)
        i ∉ [1,2,3,4] && error("Invalid SysLED index: $i")

        # Export system for testing
        basedir = export_led(i)

        # open filestream
        brightness_filestream = open("$(basedir)/beaglebone:green:usr$(i-1)/brightness","r+")
        return new(i, basedir, brightness_filestream)
    end
end

"""
    write!(led::SysLED, val::Bool, debug::Bool=false)
Turns the LED 'SysLed' on/off for val = true/false respectively.
"""
function write!(led::SysLED, entry::String, debug::Bool=false)
    debug && return
    entry ∉ ["0", "1"] && error("Invalid SysLED entry $(entry), valid options are 0 and 1 ::String")
    seekstart(led.filestream)
    write(led.filestream, "$entry\n")
    flush(led.filestream)
end

"""
    l = read(led::SysLED, debug::Bool=false)
Reads the current brightness value from the LED 'SysLED'.
"""
function read(led::SysLED, debug::Bool=false)
    debug && return
    seekstart(led.filestream)
    return readline(led.filestream)
end

"""
    teardown(led::SysLED, debug::Bool=false)
Closes all open filestreams for the SysLED 'led'.
"""
function teardown(led::SysLED, debug::Bool=false)
  debug && return
  close(led.filestream)

  if isdefined(:RUNNING_TESTS)
    # Remove the dummy file system for testing
    try
      #println("$(led.basedir)/beaglebone:green:usr$(led.i-1)")
      rm("$(led.basedir)/beaglebone:green:usr$(led.i-1)"; recursive=true)
    catch
      error("Could not remove the requested LED testfiles for channel beaglebone:green:usr$(led.i-1).")
    end
  end
end

"""
    export_led(i::Int32, debug::Bool=false)
Exports a dummy filesystem for testing the LED implementation
"""
function export_led(i::Int32, debug::Bool=false)
  debug && return

  if isdefined(:RUNNING_TESTS)
    # Export a dummy file system for testing
    basedir = "$(pwd())/testfilesystem/leds"
    try
      #println("$(basedir)/beaglebone:green:usr$(i-1)")
      mkpath("$(basedir)/beaglebone:green:usr$(i-1)")
    catch
      error("Could not export the LED device for beaglebone:green:usr$(i-1)) for testing as the directory $(basedir)/beaglebone:green:usr$(i-1) already exists.")
    end
    try
      f = open("$(basedir)/beaglebone:green:usr$(i-1)/brightness", "w"); write(f,"0"); close(f);
    catch
      error("Could not open the requested LED testfiles for beaglebone:green:usr$(i-1)/brightness.")
    end
  else
    basedir = "/sys/class/leds"
  end
  return basedir
end

"""
    to_string(led::SysLED, debug::Bool=false)
Generates a string representation of the GPIO device.
"""
function to_string(led::SysLED, debug::Bool=false)
  debug && return
  message = "\nID: $(led.i)\n\nAvailable filestream:\n"
  message = string(message, "  name=$(led.filestream.name) - write/read=$(iswritable(led.filestream))/$(isreadable(led.filestream))\n")
  return message
end


text:
###############################################################################
#
#   arb.jl : Arb real numbers
#
#   Copyright (C) 2015 Tommy Hofmann
#   Copyright (C) 2015 Fredrik Johansson
#
###############################################################################

import Base: ceil

import SpecialFunctions: digamma, zeta, polygamma, erf, erfi, erfc,
                         besselj, bessely, besseli, besselk

export ball, radius, midpoint, contains, contains_zero,
       contains_negative, contains_positive, contains_nonnegative,
       contains_nonpositive, iszero,
       isnonzero, isexact, isint, ispositive, isfinite,
       isnonnegative, isnegative, isnonpositive, add!, mul!,
       sub!, div!, strongequal, prec, overlaps, unique_integer,
       accuracy_bits, trim, ldexp, setunion,
       const_pi, const_e, const_log2, const_log10, const_euler,
       const_catalan, const_khinchin, const_glaisher,
       floor, ceil, hypot, sqrt, rsqrt, sqrt1pm1, root,
       log, log1p, expm1, sin, cos, sinpi, cospi, tan, cot,
       tanpi, cotpi, sinh, cosh, tanh, coth, atan, asin, acos,
       atanh, asinh, acosh, gamma, lgamma, rgamma, digamma, zeta,
       sincos, sincospi, sinhcosh, atan2,
       agm, fac, binom, fib, bernoulli, risingfac, risingfac2, polylog,
       chebyshev_t, chebyshev_t2, chebyshev_u, chebyshev_u2, bell, numpart

###############################################################################
#
#   Basic manipulation
#
###############################################################################

elem_type(::Type{ArbField}) = arb

parent_type(::Type{arb}) = ArbField

doc"""
    base_ring(R::ArbField)
> Returns `Union{}` since an Arb field does not depend on any other ring.
"""
base_ring(R::ArbField) = Union{}

doc"""
    base_ring(x::arb)
> Returns `Union{}` since an Arb field does not depend on any other ring.
"""
base_ring(x::arb) = Union{}

doc"""
    parent(x::arb)
> Return the parent of the given Arb field element.
"""
parent(x::arb) = x.parent

isdomain_type(::Type{arb}) = true

isexact_type(::Type{arb}) = false

doc"""
    zero(R::ArbField)
> Return exact zero in the given Arb field.
"""
zero(R::ArbField) = R(0)

doc"""
    one(R::ArbField)
> Return exact one in the given Arb field.
"""
one(R::ArbField) = R(1)

# TODO: Add hash (and document under arb basic functionality)

doc"""
    accuracy_bits(x::arb)
> Return the relative accuracy of $x$ measured in bits, capped between
> `typemax(Int)` and `-typemax(Int)`.
"""
function accuracy_bits(x::arb)
  return ccall((:arb_rel_accuracy_bits, :libarb), Int, (Ref{arb},), x)
end

function deepcopy_internal(a::arb, dict::ObjectIdDict)
  b = parent(a)()
  ccall((:arb_set, :libarb), Void, (Ref{arb}, Ref{arb}), b, a)
  return b
end

################################################################################
#
#  Conversions
#
################################################################################

doc"""
    convert(::Type{Float64}, x::arb)
> Return the midpoint of $x$ rounded down to a machine double.
"""
function convert(::Type{Float64}, x::arb)
    t = ccall((:arb_mid_ptr, :libarb), Ref{arf_struct}, (Ref{arb}, ), x)
    # 4 == round to nearest
    return ccall((:arf_get_d, :libarb), Float64, (Ref{arf_struct}, Int), t, 4)
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::ArbField)
  print(io, "Real Field with ")
  print(io, prec(x))
  print(io, " bits of precision and error bounds")
end

function show(io::IO, x::arb)
  d = ceil(parent(x).prec * 0.30102999566398119521)
  cstr = ccall((:arb_get_str, :libarb), Ref{UInt8}, (Ref{arb}, Int, UInt),
                                                  x, Int(d), UInt(0))
  print(io, unsafe_string(cstr))
  ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
end

needs_parentheses(x::arb) = false

################################################################################
#
#  Containment
#
################################################################################

doc"""
    overlaps(x::arb, y::arb)
> Returns `true` if any part of the ball $x$ overlaps any part of the ball $y$,
> otherwise return `false`.
"""
function overlaps(x::arb, y::arb)
  r = ccall((:arb_overlaps, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y)
  return Bool(r)
end

#function contains(x::arb, y::arf)
#  r = ccall((:arb_contains_arf, :libarb), Cint, (Ref{arb}, Ref{arf}), x, y)
#  return Bool(r)
#end

doc"""
    contains(x::arb, y::fmpq)
> Returns `true` if the ball $x$ contains the given rational value, otherwise
> return `false`.
"""
function contains(x::arb, y::fmpq)
  r = ccall((:arb_contains_fmpq, :libarb), Cint, (Ref{arb}, Ref{fmpq}), x, y)
  return Bool(r)
end

doc"""
    contains(x::arb, y::fmpz)
> Returns `true` if the ball $x$ contains the given integer value, otherwise
> return `false`.
"""
function contains(x::arb, y::fmpz)
  r = ccall((:arb_contains_fmpz, :libarb), Cint, (Ref{arb}, Ref{fmpz}), x, y)
  return Bool(r)
end

function contains(x::arb, y::Int)
  r = ccall((:arb_contains_si, :libarb), Cint, (Ref{arb}, Int), x, y)
  return Bool(r)
end

doc"""
    contains(x::arb, y::Integer)
> Returns `true` if the ball $x$ contains the given integer value, otherwise
> return `false`.
"""
contains(x::arb, y::Integer) = contains(x, fmpz(y))

doc"""
    contains(x::arb, y::Rational{Integer})
> Returns `true` if the ball $x$ contains the given rational value, otherwise
> return `false`.
"""
contains(x::arb, y::Rational{T}) where {T <: Integer} = contains(x, fmpq(y))

doc"""
    contains(x::arb, y::BigFloat)
> Returns `true` if the ball $x$ contains the given floating point value,
> otherwise return `false`.
"""
function contains(x::arb, y::BigFloat)
  r = ccall((:arb_contains_mpfr, :libarb), Cint,
              (Ref{arb}, Ref{BigFloat}), x, y)
  return Bool(r)
end

doc"""
    contains(x::arb, y::arb)
> Returns `true` if the ball $x$ contains the ball $y$, otherwise return
> `false`.
"""
function contains(x::arb, y::arb)
  r = ccall((:arb_contains, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y)
  return Bool(r)
end

doc"""
    contains_zero(x::arb)
> Returns `true` if the ball $x$ contains zero, otherwise return `false`.
"""
function contains_zero(x::arb)
   r = ccall((:arb_contains_zero, :libarb), Cint, (Ref{arb}, ), x)
   return Bool(r)
end

doc"""
    contains_negative(x::arb)
> Returns `true` if the ball $x$ contains any negative value, otherwise return
> `false`.
"""
function contains_negative(x::arb)
   r = ccall((:arb_contains_negative, :libarb), Cint, (Ref{arb}, ), x)
   return Bool(r)
end

doc"""
    contains_positive(x::arb)
> Returns `true` if the ball $x$ contains any positive value, otherwise return
> `false`.
"""
function contains_positive(x::arb)
   r = ccall((:arb_contains_positive, :libarb), Cint, (Ref{arb}, ), x)
   return Bool(r)
end

doc"""
    contains_nonnegative(x::arb)
> Returns `true` if the ball $x$ contains any nonnegative value, otherwise
> return `false`.
"""
function contains_nonnegative(x::arb)
   r = ccall((:arb_contains_nonnegative, :libarb), Cint, (Ref{arb}, ), x)
   return Bool(r)
end

doc"""
    contains_nonpositive(x::arb)
> Returns `true` if the ball $x$ contains any nonpositive value, otherwise
> return `false`.
"""
function contains_nonpositive(x::arb)
   r = ccall((:arb_contains_nonpositive, :libarb), Cint, (Ref{arb}, ), x)
   return Bool(r)
end

################################################################################
#
#  Comparison
#
################################################################################

doc"""
    isequal(x::arb, y::arb)
> Return `true` if the balls $x$ and $y$ are precisely equal, i.e. have the
> same midpoints and radii.
"""
function isequal(x::arb, y::arb)
  r = ccall((:arb_equal, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y)
  return Bool(r)
end

function ==(x::arb, y::arb)
    return Bool(ccall((:arb_eq, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

function !=(x::arb, y::arb)
    return Bool(ccall((:arb_ne, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

function >(x::arb, y::arb)
    return Bool(ccall((:arb_gt, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

function >=(x::arb, y::arb)
    return Bool(ccall((:arb_ge, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

function <(x::arb, y::arb)
    return Bool(ccall((:arb_lt, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

function <=(x::arb, y::arb)
    return Bool(ccall((:arb_le, :libarb), Cint, (Ref{arb}, Ref{arb}), x, y))
end

==(x::arb, y::Int) = x == arb(y)
!=(x::arb, y::Int) = x != arb(y)
<=(x::arb, y::Int) = x <= arb(y)
>=(x::arb, y::Int) = x >= arb(y)
<(x::arb, y::Int) = x < arb(y)
>(x::arb, y::Int) = x > arb(y)

==(x::Int, y::arb) = arb(x) == y
!=(x::Int, y::arb) = arb(x) != y
<=(x::Int, y::arb) = arb(x) <= y
>=(x::Int, y::arb) = arb(x) >= y
<(x::Int, y::arb) = arb(x) < y
>(x::Int, y::arb) = arb(x) > y

==(x::arb, y::fmpz) = x == arb(y)
!=(x::arb, y::fmpz) = x != arb(y)
<=(x::arb, y::fmpz) = x <= arb(y)
>=(x::arb, y::fmpz) = x >= arb(y)
<(x::arb, y::fmpz) = x < arb(y)
>(x::arb, y::fmpz) = x > arb(y)

==(x::fmpz, y::arb) = arb(x) == y
!=(x::fmpz, y::arb) = arb(x) != y
<=(x::fmpz, y::arb) = arb(x) <= y
>=(x::fmpz, y::arb) = arb(x) >= y
<(x::fmpz, y::arb) = arb(x) < y
>(x::fmpz, y::arb) = arb(x) > y

==(x::arb, y::Integer) = x == fmpz(y)
!=(x::arb, y::Integer) = x != fmpz(y)
<=(x::arb, y::Integer) = x <= fmpz(y)
>=(x::arb, y::Integer) = x >= fmpz(y)
<(x::arb, y::Integer) = x < fmpz(y)
>(x::arb, y::Integer) = x > fmpz(y)

==(x::Integer, y::arb) = fmpz(x) == y
!=(x::Integer, y::arb) = fmpz(x) != y
<=(x::Integer, y::arb) = fmpz(x) <= y
>=(x::Integer, y::arb) = fmpz(x) >= y
<(x::Integer, y::arb) = fmpz(x) < y
>(x::Integer, y::arb) = fmpz(x) > y

==(x::arb, y::Float64) = x == arb(y)
!=(x::arb, y::Float64) = x != arb(y)
<=(x::arb, y::Float64) = x <= arb(y)
>=(x::arb, y::Float64) = x >= arb(y)
<(x::arb, y::Float64) = x < arb(y)
>(x::arb, y::Float64) = x > arb(y)

==(x::Float64, y::arb) = arb(x) == y
!=(x::Float64, y::arb) = arb(x) != y
<=(x::Float64, y::arb) = arb(x) <= y
>=(x::Float64, y::arb) = arb(x) >= y
<(x::Float64, y::arb) = arb(x) < y
>(x::Float64, y::arb) = arb(x) > y

==(x::arb, y::BigFloat) = x == arb(y)
!=(x::arb, y::BigFloat) = x != arb(y)
<=(x::arb, y::BigFloat) = x <= arb(y)
>=(x::arb, y::BigFloat) = x >= arb(y)
<(x::arb, y::BigFloat) = x < arb(y)
>(x::arb, y::BigFloat) = x > arb(y)

==(x::BigFloat, y::arb) = arb(x) == y
!=(x::BigFloat, y::arb) = arb(x) != y
<=(x::BigFloat, y::arb) = arb(x) <= y
>=(x::BigFloat, y::arb) = arb(x) >= y
<(x::BigFloat, y::arb) = arb(x) < y
>(x::BigFloat, y::arb) = arb(x) > y

==(x::arb, y::fmpq) = x == arb(y, prec(parent(x)))
!=(x::arb, y::fmpq) = x != arb(y, prec(parent(x)))
<=(x::arb, y::fmpq) = x <= arb(y, prec(parent(x)))
>=(x::arb, y::fmpq) = x >= arb(y, prec(parent(x)))
<(x::arb, y::fmpq) = x < arb(y, prec(parent(x)))
>(x::arb, y::fmpq) = x > arb(y, prec(parent(x)))

==(x::fmpq, y::arb) = arb(x, prec(parent(y))) == y
!=(x::fmpq, y::arb) = arb(x, prec(parent(y))) != y
<=(x::fmpq, y::arb) = arb(x, prec(parent(y))) <= y
>=(x::fmpq, y::arb) = arb(x, prec(parent(y))) >= y
<(x::fmpq, y::arb) = arb(x, prec(parent(y))) < y
>(x::fmpq, y::arb) = arb(x, prec(parent(y))) > y

==(x::arb, y::Rational{T}) where {T <: Integer} = x == fmpq(y)
!=(x::arb, y::Rational{T}) where {T <: Integer} = x != fmpq(y)
<=(x::arb, y::Rational{T}) where {T <: Integer} = x <= fmpq(y)
>=(x::arb, y::Rational{T}) where {T <: Integer} = x >= fmpq(y)
<(x::arb, y::Rational{T}) where {T <: Integer} = x < fmpq(y)
>(x::arb, y::Rational{T}) where {T <: Integer} = x > fmpq(y)

==(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) == y
!=(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) != y
<=(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) <= y
>=(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) >= y
<(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) < y
>(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) > y

################################################################################
#
#  Predicates
#
################################################################################

function isunit(x::arb)
   !iszero(x)
end

doc"""
    iszero(x::arb)
> Return `true` if $x$ is certainly zero, otherwise return `false`.
"""
function iszero(x::arb)
   return Bool(ccall((:arb_is_zero, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isnonzero(x::arb)
> Return `true` if $x$ is certainly not equal to zero, otherwise return
> `false`.
"""
function isnonzero(x::arb)
   return Bool(ccall((:arb_is_nonzero, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isone(x::arb)
> Return `true` if $x$ is certainly not equal to oneo, otherwise return
> `false`.
"""
function isone(x::arb)
   return Bool(ccall((:arb_is_one, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isfinite(x::arb)
> Return `true` if $x$ is finite, i.e. having finite midpoint and radius,
> otherwise return `false`.
"""
function isfinite(x::arb)
   return Bool(ccall((:arb_is_finite, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isexact(x::arb)
> Return `true` if $x$ is exact, i.e. has zero radius, otherwise return
> `false`.
"""
function isexact(x::arb)
   return Bool(ccall((:arb_is_exact, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isint(x::arb)
> Return `true` if $x$ is an exact integer, otherwise return `false`.
"""
function isint(x::arb)
   return Bool(ccall((:arb_is_int, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    ispositive(x::arb)
> Return `true` if $x$ is certainly positive, otherwise return `false`.
"""
function ispositive(x::arb)
   return Bool(ccall((:arb_is_positive, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isnonnegative(x::arb)
> Return `true` if $x$ is certainly nonnegative, otherwise return `false`.
"""
function isnonnegative(x::arb)
   return Bool(ccall((:arb_is_nonnegative, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isnegative(x::arb)
> Return `true` if $x$ is certainly negative, otherwise return `false`.
"""
function isnegative(x::arb)
   return Bool(ccall((:arb_is_negative, :libarb), Cint, (Ref{arb},), x))
end

doc"""
    isnonpositive(x::arb)
> Return `true` if $x$ is certainly nonpositive, otherwise return `false`.
"""
function isnonpositive(x::arb)
   return Bool(ccall((:arb_is_nonpositive, :libarb), Cint, (Ref{arb},), x))
end

################################################################################
#
#  Parts of numbers
#
################################################################################

doc"""
    ball(mid::arb, rad::arb)
> Constructs an `arb` enclosing the range $[m-|r|, m+|r|]$, given the pair
> $(m, r)$.
"""
function ball(mid::arb, rad::arb)
  z = arb(mid, rad)
  z.parent = parent(mid)
  return z
end

doc"""
    radius(x::arb)
> Return the radius of the ball $x$ as an Arb ball.
"""
function radius(x::arb)
  z = parent(x)()
  ccall((:arb_get_rad_arb, :libarb), Void, (Ref{arb}, Ref{arb}), z, x)
  return z
end

doc"""
    midpoint(x::arb)
> Return the midpoint of the ball $x$ as an Arb ball.
"""
function midpoint(x::arb)
  z = parent(x)()
  ccall((:arb_get_mid_arb, :libarb), Void, (Ref{arb}, Ref{arb}), z, x)
  return z
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::arb)
  z = parent(x)()
  ccall((:arb_neg, :libarb), Void, (Ref{arb}, Ref{arb}), z, x)
  return z
end

################################################################################
#
#  Binary operations
#
################################################################################

for (s,f) in ((:+,"arb_add"), (:*,"arb_mul"), (://, "arb_div"), (:-,"arb_sub"))
  @eval begin
    function ($s)(x::arb, y::arb)
      z = parent(x)()
      ccall(($f, :libarb), Void, (Ref{arb}, Ref{arb}, Ref{arb}, Int),
                           z, x, y, parent(x).prec)
      return z
    end
  end
end

for (f,s) in ((:+, "add"), (:*, "mul"))
  @eval begin
    #function ($f)(x::arb, y::arf)
    #  z = parent(x)()
    #  ccall(($("arb_"*s*"_arf"), :libarb), Void,
    #              (Ref{arb}, Ref{arb}, Ref{arf}, Int),
    #              z, x, y, parent(x).prec)
    #  return z
    #end

    #($f)(x::arf, y::arb) = ($f)(y, x)

    function ($f)(x::arb, y::UInt)
      z = parent(x)()
      ccall(($("arb_"*s*"_ui"), :libarb), Void,
                  (Ref{arb}, Ref{arb}, UInt, Int),
                  z, x, y, parent(x).prec)
      return z
    end

    ($f)(x::UInt, y::arb) = ($f)(y, x)

    function ($f)(x::arb, y::Int)
      z = parent(x)()
      ccall(($("arb_"*s*"_si"), :libarb), Void,
      (Ref{arb}, Ref{arb}, Int, Int), z, x, y, parent(x).prec)
      return z
    end

    ($f)(x::Int, y::arb) = ($f)(y,x)

    function ($f)(x::arb, y::fmpz)
      z = parent(x)()
      ccall(($("arb_"*s*"_fmpz"), :libarb), Void,
                  (Ref{arb}, Ref{arb}, Ref{fmpz}, Int),
                  z, x, y, parent(x).prec)
      return z
    end

    ($f)(x::fmpz, y::arb) = ($f)(y,x)
  end
end

#function -(x::arb, y::arf)
#  z = parent(x)()
#  ccall((:arb_sub_arf, :libarb), Void,
#              (Ref{arb}, Ref{arb}, Ref{arf}, Int), z, x, y, parent(x).prec)
#  return z
#end

#-(x::arf, y::arb) = -(y - x)

function -(x::arb, y::UInt)
  z = parent(x)()
  ccall((:arb_sub_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, y, parent(x).prec)
  return z
end

-(x::UInt, y::arb) = -(y - x)

function -(x::arb, y::Int)
  z = parent(x)()
  ccall((:arb_sub_si, :libarb), Void,
              (Ref{arb}, Ref{arb}, Int, Int), z, x, y, parent(x).prec)
  return z
end

-(x::Int, y::arb) = -(y - x)

function -(x::arb, y::fmpz)
  z = parent(x)()
  ccall((:arb_sub_fmpz, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{fmpz}, Int),
              z, x, y, parent(x).prec)
  return z
end

-(x::fmpz, y::arb) = -(y-x)

+(x::arb, y::Integer) = x + fmpz(y)

-(x::arb, y::Integer) = x - fmpz(y)

*(x::arb, y::Integer) = x*fmpz(y)

//(x::arb, y::Integer) = x//fmpz(y)

+(x::Integer, y::arb) = fmpz(x) + y

-(x::Integer, y::arb) = fmpz(x) - y

*(x::Integer, y::arb) = fmpz(x)*y

//(x::Integer, y::arb) = fmpz(x)//y

#function //(x::arb, y::arf)
#  z = parent(x)()
#  ccall((:arb_div_arf, :libarb), Void,
#              (Ref{arb}, Ref{arb}, Ref{arf}, Int), z, x, y, parent(x).prec)
#  return z
#end

function //(x::arb, y::UInt)
  z = parent(x)()
  ccall((:arb_div_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, y, parent(x).prec)
  return z
end

function //(x::arb, y::Int)
  z = parent(x)()
  ccall((:arb_div_si, :libarb), Void,
              (Ref{arb}, Ref{arb}, Int, Int), z, x, y, parent(x).prec)
  return z
end

function //(x::arb, y::fmpz)
  z = parent(x)()
  ccall((:arb_div_fmpz, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{fmpz}, Int),
              z, x, y, parent(x).prec)
  return z
end

function //(x::UInt, y::arb)
  z = parent(y)()
  ccall((:arb_ui_div, :libarb), Void,
              (Ref{arb}, UInt, Ref{arb}, Int), z, x, y, parent(y).prec)
  return z
end

function //(x::Int, y::arb)
  z = parent(y)()
  t = arb(x)
  ccall((:arb_div, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, t, y, parent(y).prec)
  return z
end

function //(x::fmpz, y::arb)
  z = parent(y)()
  t = arb(x)
  ccall((:arb_div, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, t, y, parent(y).prec)
  return z
end

function ^(x::arb, y::arb)
  z = parent(x)()
  ccall((:arb_pow, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, parent(x).prec)
  return z
end

function ^(x::arb, y::fmpz)
  z = parent(x)()
  ccall((:arb_pow_fmpz, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{fmpz}, Int),
              z, x, y, parent(x).prec)
  return z
end

^(x::arb, y::Integer) = x^fmpz(y)

function ^(x::arb, y::UInt)
  z = parent(x)()
  ccall((:arb_pow_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, y, parent(x).prec)
  return z
end

function ^(x::arb, y::fmpq)
  z = parent(x)()
  ccall((:arb_pow_fmpq, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{fmpq}, Int),
              z, x, y, parent(x).prec)
  return z
end

+(x::fmpq, y::arb) = parent(y)(x) + y
+(x::arb, y::fmpq) = x + parent(x)(y)
-(x::fmpq, y::arb) = parent(y)(x) - y
//(x::arb, y::fmpq) = x//parent(x)(y)
//(x::fmpq, y::arb) = parent(y)(x)//y
-(x::arb, y::fmpq) = x - parent(x)(y)
*(x::fmpq, y::arb) = parent(y)(x) * y
*(x::arb, y::fmpq) = x * parent(x)(y)
^(x::fmpq, y::arb) = parent(y)(x) ^ y

+(x::Float64, y::arb) = parent(y)(x) + y
+(x::arb, y::Float64) = x + parent(x)(y)
-(x::Float64, y::arb) = parent(y)(x) - y
//(x::arb, y::Float64) = x//parent(x)(y)
//(x::Float64, y::arb) = parent(y)(x)//y
-(x::arb, y::Float64) = x - parent(x)(y)
*(x::Float64, y::arb) = parent(y)(x) * y
*(x::arb, y::Float64) = x * parent(x)(y)
^(x::Float64, y::arb) = parent(y)(x) ^ y
^(x::arb, y::Float64) = x ^ parent(x)(y)

+(x::BigFloat, y::arb) = parent(y)(x) + y
+(x::arb, y::BigFloat) = x + parent(x)(y)
-(x::BigFloat, y::arb) = parent(y)(x) - y
//(x::arb, y::BigFloat) = x//parent(x)(y)
//(x::BigFloat, y::arb) = parent(y)(x)//y
-(x::arb, y::BigFloat) = x - parent(x)(y)
*(x::BigFloat, y::arb) = parent(y)(x) * y
*(x::arb, y::BigFloat) = x * parent(x)(y)
^(x::BigFloat, y::arb) = parent(y)(x) ^ y
^(x::arb, y::BigFloat) = x ^ parent(x)(y)

+(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) + y
+(x::arb, y::Rational{T}) where {T <: Integer} = x + fmpq(y)
-(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) - y
-(x::arb, y::Rational{T}) where {T <: Integer} = x - fmpq(y)
//(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x)//y
//(x::arb, y::Rational{T}) where {T <: Integer} = x//fmpq(y)
*(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) * y
*(x::arb, y::Rational{T}) where {T <: Integer} = x * fmpq(y)
^(x::Rational{T}, y::arb) where {T <: Integer} = fmpq(x) ^ y
^(x::arb, y::Rational{T}) where {T <: Integer} = x ^ fmpq(y)

/(x::arb, y::arb) = x // y
/(x::fmpz, y::arb) = x // y
/(x::arb, y::fmpz) = x // y
/(x::Int, y::arb) = x // y
/(x::arb, y::Int) = x // y
/(x::UInt, y::arb) = x // y
/(x::arb, y::UInt) = x // y
/(x::fmpq, y::arb) = x // y
/(x::arb, y::fmpq) = x // y
/(x::Rational{T}, y::arb) where {T <: Integer} = x // y
/(x::arb, y::Rational{T}) where {T <: Integer} = x // y

divexact(x::arb, y::arb) = x // y
divexact(x::fmpz, y::arb) = x // y
divexact(x::arb, y::fmpz) = x // y
divexact(x::Int, y::arb) = x // y
divexact(x::arb, y::Int) = x // y
divexact(x::UInt, y::arb) = x // y
divexact(x::arb, y::UInt) = x // y
divexact(x::fmpq, y::arb) = x // y
divexact(x::arb, y::fmpq) = x // y
divexact(x::Rational{T}, y::arb) where {T <: Integer} = x // y
divexact(x::arb, y::Rational{T}) where {T <: Integer} = x // y

################################################################################
#
#  Absolute value
#
################################################################################

doc"""
    abs(x::arb)
> Return the absolute value of $x$.
"""
function abs(x::arb)
  z = parent(x)()
  ccall((:arb_abs, :libarb), Void, (Ref{arb}, Ref{arb}), z, x)
  return z
end

################################################################################
#
#  Inverse
#
################################################################################

doc"""
    inv(x::arb)
> Return the multiplicative inverse of $x$, i.e. $1/x$.
"""
function inv(x::arb)
  z = parent(x)()
  ccall((:arb_inv, :libarb), Void,
              (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
  return parent(x)(z)
end

################################################################################
#
#  Shifting
#
################################################################################

doc"""
    ldexp(x::arb, y::Int)
> Return $2^yx$. Note that $y$ can be positive, zero or negative.
"""
function ldexp(x::arb, y::Int)
  z = parent(x)()
  ccall((:arb_mul_2exp_si, :libarb), Void,
              (Ref{arb}, Ref{arb}, Int), z, x, y)
  return z
end

doc"""
    ldexp(x::arb, y::fmpz)
> Return $2^yx$. Note that $y$ can be positive, zero or negative.
"""
function ldexp(x::arb, y::fmpz)
  z = parent(x)()
  ccall((:arb_mul_2exp_fmpz, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{fmpz}), z, x, y)
  return z
end

################################################################################
#
#  Miscellaneous
#
################################################################################

doc"""
    trim(x::arb)
> Return an `arb` interval containing $x$ but which may be more economical,
> by rounding off insignificant bits from the midpoint.
"""
function trim(x::arb)
  z = parent(x)()
  ccall((:arb_trim, :libarb), Void, (Ref{arb}, Ref{arb}), z, x)
  return z
end

doc"""
    unique_integer(x::arb)
> Return a pair where the first value is a boolean and the second is an `fmpz`
> integer. The boolean indicates whether the interval $x$ contains a unique
> integer. If this is the case, the second return value is set to this unique
> integer.
"""
function unique_integer(x::arb)
  z = fmpz()
  unique = ccall((:arb_get_unique_fmpz, :libarb), Int,
    (Ref{fmpz}, Ref{arb}), z, x)
  return (unique != 0, z)
end

doc"""
    setunion(x::arb, y::arb)
> Return an `arb` containing the union of the intervals represented by $x$ and
> $y$.
"""
function setunion(x::arb, y::arb)
  z = parent(x)()
  ccall((:arb_union, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, parent(x).prec)
  return z
end

################################################################################
#
#  Constants
#
################################################################################

doc"""
    const_pi(r::ArbField)
> Return $\pi = 3.14159\ldots$ as an element of $r$.
"""
function const_pi(r::ArbField)
  z = r()
  ccall((:arb_const_pi, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_e(r::ArbField)
> Return $e = 2.71828\ldots$ as an element of $r$.
"""
function const_e(r::ArbField)
  z = r()
  ccall((:arb_const_e, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_log2(r::ArbField)
> Return $\log(2) = 0.69314\ldots$ as an element of $r$.
"""
function const_log2(r::ArbField)
  z = r()
  ccall((:arb_const_log2, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_log10(r::ArbField)
> Return $\log(10) = 2.302585\ldots$ as an element of $r$.
"""
function const_log10(r::ArbField)
  z = r()
  ccall((:arb_const_log10, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_euler(r::ArbField)
> Return Euler's constant $\gamma = 0.577215\ldots$ as an element of $r$.
"""
function const_euler(r::ArbField)
  z = r()
  ccall((:arb_const_euler, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_catalan(r::ArbField)
> Return Catalan's constant $C = 0.915965\ldots$ as an element of $r$.
"""
function const_catalan(r::ArbField)
  z = r()
  ccall((:arb_const_catalan, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_khinchin(r::ArbField)
> Return Khinchin's constant $K = 2.685452\ldots$ as an element of $r$.
"""
function const_khinchin(r::ArbField)
  z = r()
  ccall((:arb_const_khinchin, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

doc"""
    const_glaisher(r::ArbField)
> Return Glaisher's constant $A = 1.282427\ldots$ as an element of $r$.
"""
function const_glaisher(r::ArbField)
  z = r()
  ccall((:arb_const_glaisher, :libarb), Void, (Ref{arb}, Int), z, prec(r))
  return z
end

################################################################################
#
#  Real valued functions
#
################################################################################

# real - real functions

doc"""
    floor(x::arb)
> Compute the floor of $x$, i.e. the greatest integer not exceeding $x$, as an
> Arb.
"""
function floor(x::arb)
   z = parent(x)()
   ccall((:arb_floor, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    ceil(x::arb)
> Return the ceiling of $x$, i.e. the least integer not less than $x$, as an
> Arb.
"""
function ceil(x::arb)
   z = parent(x)()
   ccall((:arb_ceil, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sqrt(x::arb)
> Return the square root of $x$.
"""
function sqrt(x::arb)
   z = parent(x)()
   ccall((:arb_sqrt, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    rsqrt(x::arb)
> Return the reciprocal of the square root of $x$, i.e. $1/\sqrt{x}$.
"""
function rsqrt(x::arb)
   z = parent(x)()
   ccall((:arb_rsqrt, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sqrt1pm1(x::arb)
> Return $\sqrt{1+x}-1$, evaluated accurately for small $x$.
"""
function sqrt1pm1(x::arb)
   z = parent(x)()
   ccall((:arb_sqrt1pm1, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    log(x::arb)
> Return the principal branch of the logarithm of $x$.
"""
function log(x::arb)
   z = parent(x)()
   ccall((:arb_log, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    log1p(x::arb)
> Return $\log(1+x)$, evaluated accurately for small $x$.
"""
function log1p(x::arb)
   z = parent(x)()
   ccall((:arb_log1p, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    exp(x::arb)
> Return the exponential of $x$.
"""
function Base.exp(x::arb)
   z = parent(x)()
   ccall((:arb_exp, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    expm1(x::arb)
> Return $\exp(x)-1$, evaluated accurately for small $x$.
"""
function expm1(x::arb)
   z = parent(x)()
   ccall((:arb_expm1, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sin(x::arb)
> Return the sine of $x$.
"""
function sin(x::arb)
   z = parent(x)()
   ccall((:arb_sin, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    cos(x::arb)
> Return the cosine of $x$.
"""
function cos(x::arb)
   z = parent(x)()
   ccall((:arb_cos, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sinpi(x::arb)
> Return the sine of $\pi x$.
"""
function sinpi(x::arb)
   z = parent(x)()
   ccall((:arb_sin_pi, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    cospi(x::arb)
> Return the cosine of $\pi x$.
"""
function cospi(x::arb)
   z = parent(x)()
   ccall((:arb_cos_pi, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    tan(x::arb)
> Return the tangent of $x$.
"""
function tan(x::arb)
   z = parent(x)()
   ccall((:arb_tan, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    cot(x::arb)
> Return the cotangent of $x$.
"""
function cot(x::arb)
   z = parent(x)()
   ccall((:arb_cot, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    tanpi(x::arb)
> Return the tangent of $\pi x$.
"""
function tanpi(x::arb)
   z = parent(x)()
   ccall((:arb_tan_pi, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    cotpi(x::arb)
> Return the cotangent of $\pi x$.
"""
function cotpi(x::arb)
   z = parent(x)()
   ccall((:arb_cot_pi, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sinh(x::arb)
> Return the hyperbolic sine of $x$.
"""
function sinh(x::arb)
   z = parent(x)()
   ccall((:arb_sinh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    cosh(x::arb)
> Return the hyperbolic cosine of $x$.
"""
function cosh(x::arb)
   z = parent(x)()
   ccall((:arb_cosh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    tanh(x::arb)
> Return the hyperbolic tangent of $x$.
"""
function tanh(x::arb)
   z = parent(x)()
   ccall((:arb_tanh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    coth(x::arb)
> Return the hyperbolic cotangent of $x$.
"""
function coth(x::arb)
   z = parent(x)()
   ccall((:arb_coth, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    atan(x::arb)
> Return the arctangent of $x$.
"""
function atan(x::arb)
   z = parent(x)()
   ccall((:arb_atan, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    asin(x::arb)
> Return the arcsine of $x$.
"""
function asin(x::arb)
   z = parent(x)()
   ccall((:arb_asin, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    acos(x::arb)
> Return the arccosine of $x$.
"""
function acos(x::arb)
   z = parent(x)()
   ccall((:arb_acos, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    atanh(x::arb)
> Return the hyperbolic arctangent of $x$.
"""
function atanh(x::arb)
   z = parent(x)()
   ccall((:arb_atanh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    asinh(x::arb)
> Return the hyperbolic arcsine of $x$.
"""
function asinh(x::arb)
   z = parent(x)()
   ccall((:arb_asinh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    acosh(x::arb)
> Return the hyperbolic arccosine of $x$.
"""
function acosh(x::arb)
   z = parent(x)()
   ccall((:arb_acosh, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    gamma(x::arb)
> Return the Gamma function evaluated at $x$.
"""
function gamma(x::arb)
   z = parent(x)()
   ccall((:arb_gamma, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    lgamma(x::arb)
> Return the logarithm of the Gamma function evaluated at $x$.
"""
function lgamma(x::arb)
   z = parent(x)()
   ccall((:arb_lgamma, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    rgamma(x::arb)
> Return the reciprocal of the Gamma function evaluated at $x$.
"""
function rgamma(x::arb)
   z = parent(x)()
   ccall((:arb_rgamma, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    digamma(x::arb)
> Return the  logarithmic derivative of the gamma function evaluated at $x$,
> i.e. $\psi(x)$.
"""
function digamma(x::arb)
   z = parent(x)()
   ccall((:arb_digamma, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    zeta(x::arb)
> Return the Riemann zeta function evaluated at $x$.
"""
function zeta(x::arb)
   z = parent(x)()
   ccall((:arb_zeta, :libarb), Void, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
   return z
end

doc"""
    sincos(x::arb)
> Return a tuple $s, c$ consisting of the sine $s$ and cosine $c$ of $x$.
"""
function sincos(x::arb)
  s = parent(x)()
  c = parent(x)()
  ccall((:arb_sin_cos, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), s, c, x, parent(x).prec)
  return (s, c)
end

doc"""
    sincospi(x::arb)
> Return a tuple $s, c$ consisting of the sine $s$ and cosine $c$ of $\pi x$.
"""
function sincospi(x::arb)
  s = parent(x)()
  c = parent(x)()
  ccall((:arb_sin_cos_pi, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), s, c, x, parent(x).prec)
  return (s, c)
end

doc"""
    sinpi(x::fmpq, r::ArbField)
> Return the sine of $\pi x$ in the given Arb field.
"""
function sinpi(x::fmpq, r::ArbField)
  z = r()
  ccall((:arb_sin_pi_fmpq, :libarb), Void,
        (Ref{arb}, Ref{fmpq}, Int), z, x, prec(r))
  return z
end

doc"""
    cospi(x::fmpq, r::ArbField)
>  Return the cosine of $\pi x$ in the given Arb field.
"""
function cospi(x::fmpq, r::ArbField)
  z = r()
  ccall((:arb_cos_pi_fmpq, :libarb), Void,
        (Ref{arb}, Ref{fmpq}, Int), z, x, prec(r))
  return z
end

doc"""
    sincospi(x::fmpq, r::ArbField)
> Return a tuple $s, c$ consisting of the sine and cosine of $\pi x$ in the
> given Arb field.
"""
function sincospi(x::fmpq, r::ArbField)
  s = r()
  c = r()
  ccall((:arb_sin_cos_pi_fmpq, :libarb), Void,
        (Ref{arb}, Ref{arb}, Ref{fmpq}, Int), s, c, x, prec(r))
  return (s, c)
end

doc"""
    sinhcosh(x::arb)
> Return a tuple $s, c$ consisting of the hyperbolic sine and cosine of $x$.
"""
function sinhcosh(x::arb)
  s = parent(x)()
  c = parent(x)()
  ccall((:arb_sinh_cosh, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), s, c, x, parent(x).prec)
  return (s, c)
end

doc"""
    atan2(x::arb, y::arb)
> Return atan2$(b,a) = \arg(a+bi)$.
"""
function atan2(x::arb, y::arb)
  z = parent(x)()
  ccall((:arb_atan2, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, parent(x).prec)
  return z
end

doc"""
    agm(x::arb, y::arb)
> Return the arithmetic-geometric mean of $x$ and $y$
"""
function agm(x::arb, y::arb)
  z = parent(x)()
  ccall((:arb_agm, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, parent(x).prec)
  return z
end

doc"""
    zeta(s::arb, a::arb)
> Return the Hurwitz zeta function $\zeta(s,a)$.
"""
function zeta(s::arb, a::arb)
  z = parent(s)()
  ccall((:arb_hurwitz_zeta, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, s, a, parent(s).prec)
  return z
end

doc"""
    hypot(x::arb, y::arb)
> Return $\sqrt{x^2 + y^2}$.
"""
function hypot(x::arb, y::arb)
  z = parent(x)()
  ccall((:arb_hypot, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, parent(x).prec)
  return z
end

function root(x::arb, n::UInt)
  z = parent(x)()
  ccall((:arb_root, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, n, parent(x).prec)
  return z
end

doc"""
    root(x::arb, n::Int)
> Return the $n$-th root of $x$. We require $x \geq 0$.
"""
root(x::arb, n::Int) = x < 0 ? throw(DomainError()) : root(x, UInt(n))

doc"""
    fac(x::arb)
> Return the factorial of $x$.
"""
fac(x::arb) = gamma(x+1)

function fac(n::UInt, r::ArbField)
  z = r()
  ccall((:arb_fac_ui, :libarb), Void, (Ref{arb}, UInt, Int), z, n, r.prec)
  return z
end

doc"""
    fac(n::Int, r::ArbField)
> Return the factorial of $n$ in the given Arb field.
"""
fac(n::Int, r::ArbField) = n < 0 ? fac(r(n)) : fac(UInt(n), r)

doc"""
    binom(x::arb, n::UInt)
> Return the binomial coefficient ${x \choose n}$.
"""
function binom(x::arb, n::UInt)
  z = parent(x)()
  ccall((:arb_bin_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, n, parent(x).prec)
  return z
end

doc"""
    binom(n::UInt, k::UInt, r::ArbField)
> Return the binomial coefficient ${n \choose k}$ in the given Arb field.
"""
function binom(n::UInt, k::UInt, r::ArbField)
  z = r()
  ccall((:arb_bin_uiui, :libarb), Void,
              (Ref{arb}, UInt, UInt, Int), z, n, k, r.prec)
  return z
end

doc"""
    fib(n::fmpz, r::ArbField)
> Return the $n$-th Fibonacci number in the given Arb field.
"""
function fib(n::fmpz, r::ArbField)
  z = r()
  ccall((:arb_fib_fmpz, :libarb), Void,
              (Ref{arb}, Ref{fmpz}, Int), z, n, r.prec)
  return z
end

function fib(n::UInt, r::ArbField)
  z = r()
  ccall((:arb_fib_ui, :libarb), Void,
              (Ref{arb}, UInt, Int), z, n, r.prec)
  return z
end

doc"""
    fib(n::Int, r::ArbField)
> Return the $n$-th Fibonacci number in the given Arb field.
"""
fib(n::Int, r::ArbField) = n >= 0 ? fib(UInt(n), r) : fib(fmpz(n), r)

doc"""
    gamma(x::fmpz, r::ArbField)
> Return the Gamma function evaluated at $x$ in the given Arb field.
"""
function gamma(x::fmpz, r::ArbField)
  z = r()
  ccall((:arb_gamma_fmpz, :libarb), Void,
              (Ref{arb}, Ref{fmpz}, Int), z, x, r.prec)
  return z
end

doc"""
    gamma(x::fmpq, r::ArbField)
> Return the Gamma function evaluated at $x$ in the given Arb field.
"""
function gamma(x::fmpq, r::ArbField)
  z = r()
  ccall((:arb_gamma_fmpq, :libarb), Void,
              (Ref{arb}, Ref{fmpq}, Int), z, x, r.prec)
  return z
end

function zeta(n::UInt, r::ArbField)
  z = r()
  ccall((:arb_zeta_ui, :libarb), Void,
              (Ref{arb}, UInt, Int), z, n, r.prec)
  return z
end

doc"""
    zeta(n::Int, r::ArbField)
> Return the Riemann zeta function $\zeta(n)$ as an element of the given Arb
> field.
"""
zeta(n::Int, r::ArbField) = n >= 0 ? zeta(UInt(n), r) : zeta(r(n))

function bernoulli(n::UInt, r::ArbField)
  z = r()
  ccall((:arb_bernoulli_ui, :libarb), Void,
              (Ref{arb}, UInt, Int), z, n, r.prec)
  return z
end

doc"""
    bernoulli(n::Int, r::ArbField)
> Return the $n$-th Bernoulli number as an element of the given Arb field.
"""
bernoulli(n::Int, r::ArbField) = n >= 0 ? bernoulli(UInt(n), r) : throw(DomainError)

function risingfac(x::arb, n::UInt)
  z = parent(x)()
  ccall((:arb_rising_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Int), z, x, n, parent(x).prec)
  return z
end

doc"""
    risingfac(x::arb, n::Int)
> Return the rising factorial $x(x + 1)\ldots (x + n - 1)$ as an Arb.
"""
risingfac(x::arb, n::Int) = n < 0 ? throw(DomainError()) : risingfac(x, UInt(n))

function risingfac(x::fmpq, n::UInt, r::ArbField)
  z = r()
  ccall((:arb_rising_fmpq_ui, :libarb), Void,
              (Ref{arb}, Ref{fmpq}, UInt, Int), z, x, n, r.prec)
  return z
end

doc"""
    risingfac(x::fmpq, n::Int, r::ArbField)
> Return the rising factorial $x(x + 1)\ldots (x + n - 1)$ as an element of the
> given Arb field.
"""
risingfac(x::fmpq, n::Int, r::ArbField) = n < 0 ? throw(DomainError()) : risingfac(x, UInt(n), r)

function risingfac2(x::arb, n::UInt)
  z = parent(x)()
  w = parent(x)()
  ccall((:arb_rising2_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, UInt, Int), z, w, x, n, parent(x).prec)
  return (z, w)
end

doc"""
    risingfac2(x::arb, n::Int)
> Return a tuple containing the rising factorial $x(x + 1)\ldots (x + n - 1)$
> and its derivative.
"""
risingfac2(x::arb, n::Int) = n < 0 ? throw(DomainError()) : risingfac2(x, UInt(n))

doc"""
    polylog(s::arb, a::arb)
> Return the polylogarithm Li$_s(a)$.
"""
function polylog(s::arb, a::arb)
  z = parent(s)()
  ccall((:arb_polylog, :libarb), Void,
              (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, s, a, parent(s).prec)
  return z
end

doc"""
    polylog(s::Int, a::arb)
> Return the polylogarithm Li$_s(a)$.
"""
function polylog(s::Int, a::arb)
  z = parent(a)()
  ccall((:arb_polylog_si, :libarb), Void,
              (Ref{arb}, Int, Ref{arb}, Int), z, s, a, parent(a).prec)
  return z
end

function chebyshev_t(n::UInt, x::arb)
  z = parent(x)()
  ccall((:arb_chebyshev_t_ui, :libarb), Void,
              (Ref{arb}, UInt, Ref{arb}, Int), z, n, x, parent(x).prec)
  return z
end

function chebyshev_u(n::UInt, x::arb)
  z = parent(x)()
  ccall((:arb_chebyshev_u_ui, :libarb), Void,
              (Ref{arb}, UInt, Ref{arb}, Int), z, n, x, parent(x).prec)
  return z
end

function chebyshev_t2(n::UInt, x::arb)
  z = parent(x)()
  w = parent(x)()
  ccall((:arb_chebyshev_t2_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Ref{arb}, Int), z, w, n, x, parent(x).prec)
  return z, w
end

function chebyshev_u2(n::UInt, x::arb)
  z = parent(x)()
  w = parent(x)()
  ccall((:arb_chebyshev_u2_ui, :libarb), Void,
              (Ref{arb}, Ref{arb}, UInt, Ref{arb}, Int), z, w, n, x, parent(x).prec)
  return z, w
end

doc"""
    chebyshev_t(n::Int, x::arb)
> Return the value of the Chebyshev polynomial $T_n(x)$.
"""
chebyshev_t(n::Int, x::arb) = n < 0 ? throw(DomainError()) : chebyshev_t(UInt(n), x)

doc"""
    chebyshev_u(n::Int, x::arb)
> Return the value of the Chebyshev polynomial $U_n(x)$.
"""
chebyshev_u(n::Int, x::arb) = n < 0 ? throw(DomainError()) : chebyshev_u(UInt(n), x)

doc"""
    chebyshev_t2(n::Int, x::arb)
> Return the tuple $(T_{n}(x), T_{n-1}(x))$.
"""
chebyshev_t2(n::Int, x::arb) = n < 0 ? throw(DomainError()) : chebyshev_t2(UInt(n), x)

doc"""
    chebyshev_u2(n::Int, x::arb)
> Return the tuple $(U_{n}(x), U_{n-1}(x))$
"""
chebyshev_u2(n::Int, x::arb) = n < 0 ? throw(DomainError()) : chebyshev_u2(UInt(n), x)

doc"""
    bell(n::fmpz, r::ArbField)
> Return the Bell number $B_n$ as an element of $r$.
"""
function bell(n::fmpz, r::ArbField)
  z = r()
  ccall((:arb_bell_fmpz, :libarb), Void,
              (Ref{arb}, Ref{fmpz}, Int), z, n, r.prec)
  return z
end

doc"""
    bell(n::Int, r::ArbField)
> Return the Bell number $B_n$ as an element of $r$.
"""
bell(n::Int, r::ArbField) = bell(fmpz(n), r)

doc"""
    numpart(n::fmpz, r::ArbField)
> Return the number of partitions $p(n)$ as an element of $r$.
"""
function numpart(n::fmpz, r::ArbField)
  z = r()
  ccall((:arb_partitions_fmpz, :libarb), Void,
              (Ref{arb}, Ref{fmpz}, Int), z, n, r.prec)
  return z
end

doc"""
    numpart(n::fmpz, r::ArbField)
> Return the number of partitions $p(n)$ as an element of $r$.
"""
numpart(n::Int, r::ArbField) = numpart(fmpz(n), r)

################################################################################
#
#  Unsafe operations
#
################################################################################

function zero!(z::arb)
   ccall((:arb_zero, :libarb), Void, (Ref{arb},), z)
   return z
end

for (s,f) in (("add!","arb_add"), ("mul!","arb_mul"), ("div!", "arb_div"),
              ("sub!","arb_sub"))
  @eval begin
    function ($(Symbol(s)))(z::arb, x::arb, y::arb)
      ccall(($f, :libarb), Void, (Ref{arb}, Ref{arb}, Ref{arb}, Int),
                           z, x, y, parent(x).prec)
      return z
    end
  end
end

function addeq!(z::arb, x::arb)
    ccall((:arb_add, :libarb), Void, (Ref{arb}, Ref{arb}, Ref{arb}, Int),
                           z, z, x, parent(x).prec)
    return z
end

################################################################################
#
#  Unsafe setting
#
################################################################################

for (typeofx, passtoc) in ((arb, Ref{arb}), (Ref{arb}, Ref{arb}))
  for (f,t) in (("arb_set_si", Int), ("arb_set_ui", UInt),
                ("arb_set_d", Float64))
    @eval begin
      function _arb_set(x::($typeofx), y::($t))
        ccall(($f, :libarb), Void, (($passtoc), ($t)), x, y)
      end

      function _arb_set(x::($typeofx), y::($t), p::Int)
        _arb_set(x, y)
        ccall((:arb_set_round, :libarb), Void,
                    (($passtoc), ($passtoc), Int), x, x, p)
      end
    end
  end

  @eval begin
    function _arb_set(x::($typeofx), y::fmpz)
      ccall((:arb_set_fmpz, :libarb), Void, (($passtoc), Ref{fmpz}), x, y)
    end

    function _arb_set(x::($typeofx), y::fmpz, p::Int)
      ccall((:arb_set_round_fmpz, :libarb), Void,
                  (($passtoc), Ref{fmpz}, Int), x, y, p)
    end

    function _arb_set(x::($typeofx), y::fmpq, p::Int)
      ccall((:arb_set_fmpq, :libarb), Void,
                  (($passtoc), Ref{fmpq}, Int), x, y, p)
    end

    function _arb_set(x::($typeofx), y::arb)
      ccall((:arb_set, :libarb), Void, (($passtoc), Ref{arb}), x, y)
    end

    function _arb_set(x::($typeofx), y::arb, p::Int)
      ccall((:arb_set_round, :libarb), Void,
                  (($passtoc), Ref{arb}, Int), x, y, p)
    end

    function _arb_set(x::($typeofx), y::AbstractString, p::Int)
      s = string(y)
      err = ccall((:arb_set_str, :libarb), Int32,
                  (($passtoc), Ref{UInt8}, Int), x, s, p)
      err == 0 || error("Invalid real string: $(repr(s))")
    end

    function _arb_set(x::($typeofx), y::BigFloat)
      m = ccall((:arb_mid_ptr, :libarb), Ref{arf_struct},
                  (($passtoc), ), x)
      r = ccall((:arb_rad_ptr, :libarb), Ref{mag_struct},
                  (($passtoc), ), x)
      ccall((:arf_set_mpfr, :libarb), Void,
                  (Ref{arf_struct}, Ref{BigFloat}), m, y)
      ccall((:mag_zero, :libarb), Void, (Ref{mag_struct}, ), r)
    end

    function _arb_set(x::($typeofx), y::BigFloat, p::Int)
      m = ccall((:arb_mid_ptr, :libarb), Ref{arf_struct}, (($passtoc), ), x)
      r = ccall((:arb_rad_ptr, :libarb), Ref{mag_struct}, (($passtoc), ), x)
      ccall((:arf_set_mpfr, :libarb), Void,
                  (Ref{arf_struct}, Ref{BigFloat}), m, y)
      ccall((:mag_zero, :libarb), Void, (Ref{mag_struct}, ), r)
      ccall((:arb_set_round, :libarb), Void,
                  (($passtoc), ($passtoc), Int), x, x, p)
    end
  end
end

################################################################################
#
#  Parent object overloading
#
################################################################################

function (r::ArbField)()
  z = arb()
  z.parent = r
  return z
end

function (r::ArbField)(x::Int)
  z = arb(fmpz(x), r.prec)
  z.parent = r
  return z
end

function (r::ArbField)(x::UInt)
  z = arb(fmpz(x), r.prec)
  z.parent = r
  return z
end

function (r::ArbField)(x::fmpz)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

(r::ArbField)(x::Integer) = r(fmpz(x))

function (r::ArbField)(x::fmpq)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

(r::ArbField)(x::Rational{T}) where {T <: Integer} = r(fmpq(x))

#function call(r::ArbField, x::arf)
#  z = arb(arb(x), r.prec)
#  z.parent = r
#  return z
#end

function (r::ArbField)(x::Float64)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

function (r::ArbField)(x::arb)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

function (r::ArbField)(x::AbstractString)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

function (r::ArbField)(x::Irrational)
  if x == pi
    return const_pi(r)
  elseif x == e
    return const_e(r.prec)
  else
    error("constant not supported")
  end
end

function (r::ArbField)(x::BigFloat)
  z = arb(x, r.prec)
  z.parent = r
  return z
end

################################################################################
#
#  Arb real field constructor
#
################################################################################

# see inner constructor for ArbField


text:
using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff

using LinearAlgebra
using Random
import GLPK

s = rand(1:100)
s = 98
@info "Seed $s"
Random.seed!(s)

n = Int(1e2)
k = 3000

xpi = rand(n * n);
total = sum(xpi);
xpi = reshape(xpi, n, n)
const xp = xpi # ./ total;

# better for memory consumption as we do coordinate-wise ops

function cf(x, xp)
    return LinearAlgebra.norm(x .- xp)^2 / n^2
end

function cgrad!(storage, x, xp)
    return @. storage = 2 * (x - xp) / n^2
end

# initial direction for first vertex
direction_vec = Vector{Float64}(undef, n * n)
randn!(direction_vec)
direction_mat = reshape(direction_vec, n, n)

# BirkhoffPolytopeLMO via Hungarian Method
lmo_native = FrankWolfe.BirkhoffPolytopeLMO()

# BirkhoffPolytopeLMO realized via LP solver
lmo_moi = FrankWolfe.convert_mathopt(lmo_native, GLPK.Optimizer(), dimension=n)

# choose between lmo_native (= Hungarian Method) and lmo_moi (= LP formulation solved with GLPK)
lmo = lmo_native

x00 = FrankWolfe.compute_extreme_point(lmo, direction_mat)

FrankWolfe.benchmark_oracles(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    () -> randn(n, n),
    lmo;
    k=100,
)

# BCG run

x0 = deepcopy(x00)

@time x, v, primal, dual_gap, trajectoryBCG = FrankWolfe.blended_conditional_gradient(
    x -> cf(x, xp),
    (str, x) -> cgrad!(str, x, xp),
    lmo,
    x0,
    max_iteration=k,
    L=100,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    linesearch_tol=1e-7,
    emphasis=FrankWolfe.memory,
    trajectory=true,
    verbose=true,
);

data = [trajectoryBCG]
label = ["BCG"]

FrankWolfe.plot_trajectories(data, label)


text:
if basename(pwd()) == "aoc"
    cd("2017/4")
end

part1(filename::AbstractString) = count(isvalidpassphrase.(eachline(filename)))

function isvalidpassphrase(line::AbstractString)
    words = split(line)
    length(words) == length(unique!(sort!(words)))
end

@assert part1("example.txt") == 2
part1("input.txt")

part2(filename::AbstractString) = count(isvalidpassphrase2.(eachline(filename)))

function isvalidpassphrase2(line::AbstractString)
    words = split(line)
    length(words) == length(unique!(sort!([join(sort!(split(word, ""))) for word in words])))
end

part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
## Hashing

equality_data(g::AbstractFPGroupElement) = (normalform!(g); word(g))

bitget(h::UInt, n::Int) = Bool((h & (1 << n)) >> n)
bitclear(h::UInt, n::Int) = h & ~(1 << n)
bitset(h::UInt, n::Int) = h | (1 << n)
bitset(h::UInt, v::Bool, n::Int) = v ? bitset(h, n) : bitclear(h, n)

# We store hash of a word in field `savedhash` to use it as cheap proxy to
# determine non-equal words. Additionally bits of `savehash` store boolean
# information as follows
# * `savedhash & 1` (the first bit): is word in normal form?
# * `savedhash & 2` (the second bit): is the hash valid?
const __BITFLAGS_MASK = ~(~(UInt(0)) << 2)

isnormalform(g::AbstractFPGroupElement) = bitget(g.savedhash, 0)
_isvalidhash(g::AbstractFPGroupElement) = bitget(g.savedhash, 1)

_setnormalform(h::UInt, v::Bool) = bitset(h, v, 0)
_setvalidhash(h::UInt, v::Bool) = bitset(h, v, 1)

_setnormalform!(g::AbstractFPGroupElement, v::Bool) = g.savedhash = _setnormalform(g.savedhash, v)
_setvalidhash!(g::AbstractFPGroupElement, v::Bool) = g.savedhash = _setvalidhash(g.savedhash, v)

# To update hash use this internal method, possibly only after computing the
# normal form of `g`:
function _update_savedhash!(g::AbstractFPGroupElement, data)
    h = hash(data, hash(parent(g)))
    h = (h << count_ones(__BITFLAGS_MASK)) | (__BITFLAGS_MASK & g.savedhash)
    g.savedhash = _setvalidhash(h, true)
    return g
end

function Base.hash(g::AbstractFPGroupElement, h::UInt)
    _isvalidhash(g) || _update_savedhash!(g, equality_data(g))
    return hash(g.savedhash >> count_ones(__BITFLAGS_MASK), h)
end

function Base.copyto!(res::AbstractFPGroupElement, g::AbstractFPGroupElement)
    @assert parent(res) === parent(g)
    resize!(word(res), length(word(g)))
    copyto!(word(res), word(g))
    res.savedhash = g.savedhash
    return res
end