text:
import FLOWFarm; const ff = FLOWFarm

turbine_x = [0.0,0.0]
turbine_y = [0.0,0.0]
turbine_z = [0.0,0.0]
rotor_diameter = [40.0,40.0]
hub_height = [90.0,90.0]
ai = 1.0/3.0
turbine_ai = [ai,ai]
ct = 0.7
turbine_ct = [ct,ct]
yaw = 0.0
turbine_yaw = [yaw,yaw]
ambient_ti = 0.1
turbine_local_ti = [ambient_ti,ambient_ti]

sorted_turbine_index = [1]

wind_speed = 8.1
air_density = 1.1716  # kg/m^3
winddirections = [270.0*pi/180.0]
windspeeds = [wind_speed]
windprobabilities = [1.0]
measurementheights = [hub_height[1]]
shearexponent = 0.15
wtvelocities = [wind_speed]
ambient_tis = [ambient_ti]
wind_shear_model = ff.PowerLawWindShear(shearexponent)
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)

alpha = 0.1
wakedeficitmodel = ff.JensenTopHat(alpha)
horizontal_spread_rate = alpha
wakedeflectionmodel = ff.JiminezYawDeflection(horizontal_spread_rate)
wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)

# loc = [7.0*rotor_diameter, 0.0, hub_height]
# cp = 0.8
# generator_efficiency = 0.944
#
# wind_speed = 8.1
# air_density = 1.1716  # kg/m^3
#
#
#
#
# winddirections = [270.0*pi/180.0]
# windspeeds = [wind_speed]
# windprobabilities = [1.0]
# measurementheights = [hub_height]
# shearexponent = 0.15
# turbine_inflow_velcities = [wind_speed]
# ambient_tis = [ambient_ti]
# nturbines = 1
#
# ct_model = ff.ThrustModelConstantCt(ct)
# power_model = ff.PowerModelConstantCp(cp)
# wind_shear_model = [ff.PowerLawWindShear(shearexponent)]
#
# cut_in_speed = 0.0
# cut_out_speed = 25.0
# rated_speed = 12.0
# rated_power = 1.0176371581904552e6
#
# turbine1 = ff.TurbineDefinition(1, [rotor_diameter], [hub_height], [cut_in_speed], [rated_speed], [cut_out_speed], [rated_power], [generator_efficiency], ct_model, power_model)
# turbine_definitions = [turbine1]
# turbine_definition_ids = [1]
#
#
# windfarm = ff.WindFarm(turbine_x, turbine_y, turbine_z, turbine_definition_ids, turbine_definitions)
# windfarmstate = ff.SingleWindFarmState(1, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai, turbine_inflow_velcities, [0.0], [ambient_ti],sorted_turbine_index)
# windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)
#
#
# horizontal_spread_rate = alpha
# wakedeflectionmodel = ff.JiminezYawDeflection(horizontal_spread_rate)
# wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
# localtimodel = ff.LocalTIModelNoLocalTI()
#
#
# pd1 = ff.WindFarmProblemDescription(windfarm, windresource, [windfarmstate])


text:
#
using LinearSolve

"""
 plug into OrdinaryDiffEq -
    - write all functionality of DefaultLinSolve,
      any of the stuff from DiffEqBase/linear_nonlinear.jl
    - wrap Krylov.jl, IterativeSolvers.jl
"""

A = rand(5, 5)
b = rand(5)
prob = LinearProblem(A, b)

# Factorization
A * solve(prob, LUFactorization() ) ≈ b
A * solve(prob, QRFactorization() ) ≈ b
A * solve(prob, SVDFactorization()) ≈ b
A * solve(prob, KrylovJL(A, b))     ≈ b

ext: jl


text:
function test_dict_key_change()

    # create a Dict with a single entry keyed by `x`
    d = Dict()
    x = []
    d[x] = 1

    # push a number to the `x` array and test if the key is still found in `d`
    found = Bool[]
    for i in 1:100
        push!(x, i)
        push!(found, haskey(d, x))
    end

    println(sum(found), " out of 100 tests are found")
end


text:
#' # AdventOfCode2017
#'
#' [![Build Status](https://travis-ci.org/ellisvalentiner/AdventOfCode2017.jl.svg?branch=master)](https://travis-ci.org/ellisvalentiner/AdventOfCode2017.jl)
#'
#' [![Coverage Status](https://coveralls.io/repos/ellisvalentiner/AdventOfCode2017.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/ellisvalentiner/AdventOfCode2017.jl?branch=master)
#'
#' [![codecov.io](http://codecov.io/github/ellisvalentiner/AdventOfCode2017.jl/coverage.svg?branch=master)](http://codecov.io/github/ellisvalentiner/AdventOfCode2017.jl?branch=master)

#+ echo=false
using AdventOfCode2017
using Markdown
import Markdown: MD, Table
import Statistics: median

table = ["Day" "Elapsed (s)" "Allocated (MiB)" "Garbage Collection (s)"]
for puzzle in filter(x-> x≠ :AdventOfCode2017, names(AdventOfCode2017))
    m = match(r"[0-9]+", String(puzzle))
    day = string(parse(Int, m.match))
    @info "Puzzle $day"
    @eval begin
        val, t, bytes, gctime, memallocs = @timed $puzzle()
        t = round(t, digits=4)
        kib = round(bytes/1024^2, digits=4)
        gctime = round(gctime, digits=4)
        row = [$day t kib gctime]
        table = vcat(table, row)
    end
end
table = vcat(table, ["**Median**" round.(median(table[2:end, 2:end], dims=1), digits=4)])
table = vcat(table, ["**Total**" round.(sum(table[2:end, 2:end], dims=1), digits=4)])
performance = MD(Table(Any[map(x->Any[x], table[i, :]) for i in 1:size(table, 1)], Symbol[:l, :r, :r, :r]))

#' ## Performance

#+ echo=false
#+ results="as_is"
print(performance)


text:
using NumericalAnalysis: Polynomial
using Test

x = [1.0, 1.3, 1.6, 1.9, 2.2]
y = [0.7651977, 0.6200860, 0.4554022, 0.2818186, 0.1103623]

Polynomial.NCSpline([0,1,2,3], [1, exp(1), exp(2), exp(3)])
Polynomial.CCSpline([0,1,2,3], [1, exp(1), exp(2), exp(3)],[1, exp(3)])
@testset "Polynomial test" begin
    @test abs(Polynomial.Lagrange(x,y)(1.5) - 0.5118199)<=0.0000001
    @test abs(Polynomial.NDDF(x,y)(1.5) - 0.5118199)<=0.0000001
    @test abs(abs(Polynomial.Neville(x,y,1.5) - 0.5118199)<=0.0000001)
end


text:
module Essential

include("ensemblegenemap.jl")
include("idMap.jl")

function essentialGenesHugo()
    ensembletohugo = EnsembleGeneMap.ensembleToHugo()
    essentialgenesHugo = String[]

    for line in readlines("./data/essential_9606_all_gene_status.txt")
        fields = chomp(line)
        line_parts = split(fields)
        ensembleGene = line_parts[1]
        if haskey(ensembletohugo, ensembleGene)
            push!(essentialgenesHugo, ensembletohugo[ensembleGene])
        end
    end

    return essentialgenesHugo
end

function getNodes(dbidfile)
    essentialgenesHugo = essentialGenesHugo()
    genetonode = DbIdNameMapping.geneToNodes(dbidfile)

    essentialnodes = String[]

    for hugogene in essentialgenesHugo
        if haskey(genetonode, hugogene)
            for node in genetonode[hugogene]
                push!(essentialnodes, node)
            end
        end
    end

    return essentialnodes
end

function getGenes(nodes, dbidfile)
    return collect(intersect(Set(nodes), Set(getNodes(dbidfile))))
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
n_states(p::SubHuntPOMDP) = p.size^4*4*2+1

function state_index(p::SubHuntPOMDP, s::SubState)
    if s == END_KILL
        return n_states(p)
    else
        o = s.own
        t = s.target
        return sub2ind((p.size, p.size, p.size, p.size, 4, 2), o[1], o[2], t[1], t[2], s.goal, Int(s.aware)+1)
    end
end

function state_from_index(p::SubHuntPOMDP, s::Int)
    if s == n_states(p)
        return END_KILL
    else
        is = ind2sub((p.size, p.size, p.size, p.size, 4, 2), s)
        o = SVector(is[1], is[2])
        t = SVector(is[3], is[4])
        return SubState(o, t, is[5], is[6]==2)
    end
end

# the state space is represented by the POMDP itself
# rand(rng::AbstractRNG, p::SubHuntPOMDP)

states(p::SubHuntPOMDP) = p
iterator(p::SubHuntPOMDP) = p
Base.start(p::SubHuntPOMDP) = 1
Base.done(p::SubHuntPOMDP, i::Int) = i > n_states(p)
Base.next(p::SubHuntPOMDP, i::Int) = (state_from_index(p, i), i+1)
Base.eltype(p::SubHuntPOMDP) = SubState
Base.length(p::SubHuntPOMDP) = n_states(p)

struct SubHuntInitDist
    p::SubHuntPOMDP
end
initial_state_distribution(p::SubHuntPOMDP) = SubHuntInitDist(p)

function rand(rng::AbstractRNG, d::SubHuntInitDist)
    mid = round(Int, d.p.size/2)
    own = SVector(mid, mid)
    goal = rand(rng, 1:4)
    goal_dir = DIR[goal]
    start = rand(rng, 1:d.p.size)
    if goal_dir[1] == 0
        target = SVector(start, goal_dir[2] > 0 ? 1 : d.p.size)
    else
        target = SVector(goal_dir[1] > 0 ? 1 : d.p.size, start)
    end
    return SubState(own, target, goal, false)
end

sampletype(::Type{SubHuntInitDist}) = SubState

actions(p::SubHuntPOMDP) = 1:6
n_actions(p::SubHuntPOMDP) = 6


text:
abstract type Lattice{T<:Real} end
abstract type InfiniteLattice{T<:Real} <: Lattice{T} end
abstract type FiniteLattice{T<:Real,L<:InfiniteLattice{T}} <: Lattice{T} end

"""returns a Tuple = (2D coordinates of emitter in the lattice,emitter)"""
emitter(a::L, x::Int, y::Int) where {L<:InfiniteLattice} = (a.origin + x * a.e1 + y * a.e2, a.emitter)

emitterpitchx(a::L) where {L<:Lattice} = norm(a.e₁)
emitterpitchy(a::L) where {L<:Lattice} = norm(a.e₂)

struct SizedLattice{T<:Real,L<:InfiniteLattice{T}} <: FiniteLattice{T,L}
    lattice::L
    e₁emitters::Int
    e₂emitters::Int

    function SizedLattice(lattice::L, e₁emitters::Int, e₂emitters::Int) where {T<:Real,L<:InfiniteLattice{T}}
        @assert e₁emitters >= T(0)
        @assert e₂emitters >= T(0)
        return new{T,L}(lattice, e₁emitters, e₂emitters)
    end
end

emitterpitchx(a::SizedLattice) = emitterpitchx(norm(a.lattice.e₁))
emitterpitchy(a::SizedLattice) = emitterpitchy(norm(a.lattice.e₂))

function emitter(a::FiniteLattice{T,L}, x::Int, y::Int) where {T<:Real,L<:InfiniteLattice{T}}
    @assert T(0) <= x <= a.e₁emitters
    @assert T(0) <= y <= a.e₂emitters
    return emitter(a.lattice, x, y) #because L is of type InfiteLattice and not FiniteLattice this call will not cause infinite recursion.
end

struct PrimitiveLattice{S<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real} <: InfiniteLattice{T}
    origin::SVector{2,T}
    e₁::SVector{2,T}
    e2::SVector{2,T}
    emitter::PlanarEmitter{S,P,T}

    PrimitiveLattice(e₁::SVector{2,T}, e₂::SVector{2,T}, emitter::PlanarEmitter{S,P,T}; origin = SVector{2,T}(T(0), T(0))) where {S<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real} = new{S,P,T}(e₁, e₂, emitter, origin)
end

hexagonalemitterlattice(pitch::T, emitter::PlanarEmitter{S,P,T}; origin = SVector{2,T}(T(0), T(0))) where {S<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real} = PrimitiveLattice{S,P,T}(origin, SVector{2,T}(emitterpitch, T(0)), SVector{2,T}(emitterpitch / T(2), emitterpitch * sqrt(T(3)) / T(2)), emitter)

squareemitterlattice(emitterpitch::T, emitter::PlanarEmitter{S,P,T}) where {S<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real} = PrimitiveLattice{S,P,T}(SVector{2,T}(emitterpitch, T(0)), SVector{2,T}(T(0), emitterpitch), emitter)

#This needs work. We probably want some kind of hierarchical lattice with superemitter and subemitter functions.
struct RGBPixelLattice{R<:AbstractSpectrum,G<:AbstractSpectrum,B<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real} <: InfiniteLattice{T}
    red::PrimitiveLattice{R,P,T}
    green::PrimitiveLattice{G,P,T}
    blue::PrimitiveLattice{B,P,T}

    function RGBPixelLattice(::Type{R}, ::Type{G}, ::Type{B}, ::Type{P}, subemitterwidth::T, rgbemitterpitch::T) where {R<:AbstractSpectrum,G<:AbstractSpectrum,B<:AbstractSpectrum,P<:AbstractAngularPowerDistribution,T<:Real}
        @assert 3 * subemitterwidth <= rgbemitterpitch

        subemitterpitch = rgbemitterpitch / T(3.0)
        e₁ = SVector{2,T}(rgbemitterpitch, T(0))
        e2 = SVector{2,T}(T(0), rgbemitterpitch)

        return new{R,G,B,P,T}(PrimitiveLattice(e1, e2, PlanarEmitter{R,P,T}(subemitterwidth, rgbemitterpitch)), PrimitiveLattice(e1, e2, PlanarEmitter{G,P,T}(subemitterwidth, rgbemitterpitch), origin = SVector{2,T}(subemitterpitch, T(0.0))), PrimitiveLattice(e1, e2, PlanarEmitter{B,P,T}(subemitterwidth, rgbemitterpitch), origin = SVector{2,T}(2 * subemitterpitch, T(0.0))))
    end
end


text:
mutable struct TreeLayout{T,F}
    rows::Vector{Vector{Vector{NamedTuple{(:node,:position),Tuple{T,Ref{Point{2,F}}}}}}}
    x_min::Vector{F}
    x_max::Vector{F}
    y_min::Vector{F}
    y_max::Vector{F}
    x_margin::F
    y_margin::F
end
TreeLayout(root::T; pos::Point{2,F}=zero(Point2f), x_margin=0.1,y_margin=0.1, height=(_)->1.0, width=(_)->1.0) where {T,F} =
    TreeLayout{T,F}([[[(node=root,position=Ref(pos))]]],F[pos[1]-width(root)/2],F[pos[1]+width(root)/2],F[pos[2]],F[pos[2]+height(root)],x_margin,y_margin)

function add_right_branch!(tree::TreeLayout{T,F}, branch::TreeLayout{T,F}; keep_centered=true) where {T,F}
    # add dummy root to right branch
    pushfirst!(branch.x_min, Inf)
    pushfirst!(branch.x_max, -Inf)
    pushfirst!(branch.y_min, zero(F))
    pushfirst!(branch.y_max, zero(F))

    # get all arrays to the same size
    new_length = maximum(length.((
                                    tree.x_min, tree.x_max, tree.y_min, tree.y_max,
                                    branch.x_min, branch.x_max, branch.y_min, branch.y_max
                                )))
    append!(tree.x_min, fill(Inf, new_length-length(tree.x_min)))
    append!(tree.x_max, fill(-Inf, new_length-length(tree.x_max)))
    append!(tree.y_min, fill(0, new_length-length(tree.y_min)))
    append!(tree.y_max, fill(0, new_length-length(tree.y_max)))
    append!(branch.x_min, fill(Inf, new_length-length(branch.x_min)))
    append!(branch.x_max, fill(-Inf, new_length-length(branch.x_max)))
    append!(branch.y_min, fill(0, new_length-length(branch.y_min)))
    append!(branch.y_max, fill(0, new_length-length(branch.y_max)))

    append!(tree.rows, [[NamedTuple{(:node,:position),Tuple{T,Point{2,F}}}[]] for i in 1:(new_length-length(tree.rows))])

    # compute x-offset between the trees
    tree.x_margin = max(tree.x_margin, branch.x_margin)
    tree.y_margin = max(tree.y_margin, branch.y_margin)
    x_offset = max(maximum(tree.x_max .+ tree.x_margin .- branch.x_min),zero(F))

    # compute new "hull" of the combined tree
    tree.x_min .= min.(tree.x_min, branch.x_min.+x_offset)
    tree.x_max .= max.(tree.x_max, branch.x_max.+x_offset)

    # raise row
    min_height = max.(tree.y_max .- tree.y_min, branch.y_max .- branch.y_min)

    # shift branch nodes to right positions and integrate into tree
    tmp_y = tree.y_max[1] + tree.y_margin
    for (i,(tree_row,branch_row,row_height)) in enumerate(zip(tree.rows[2:end], branch.rows, min_height[2:end]))
        for clique in branch_row
            for node in clique
                node.position[] = Point{2,F}(node.position[][1]+x_offset, tmp_y)
            end
        end

        # on the lowest level, cliques are merged
        if i==1
            append!(tree_row[1], branch_row[1])
        else
            append!(tree_row, branch_row)
        end
        tree.y_min[i+1] = tmp_y
        tree.y_max[i+1] = tmp_y+row_height

        tmp_y = tree.y_max[i+1] + tree.y_margin
    end

    # compute new common root position
    new_root_x = (tree.x_min[2] + tree.x_max[2])/2
    root_width = tree.x_max[1] - tree.x_min[1]
    tree.x_min[1] = new_root_x - root_width/2
    tree.x_max[1] = new_root_x + root_width/2
    tree.rows[1][1][1].position[] = Point{2,F}(new_root_x,tree.rows[1][1][1].position[][2])

    #maybe re-center the entire tree
    if keep_centered
        tree.x_min .-= new_root_x
        tree.x_max .-= new_root_x
        for row in tree.rows
            for clique in row
                for node in clique
                   node.position[] -= Point{2,F}(new_root_x,zero(F))
                end
            end
        end
    end
    nothing
end

function layout(dendrite::DendriteSegment; pos=zero(Point2f), params...)
    tree = TreeLayout(dendrite;pos, params...)

    for child in dendrite.children
        child_branch = layout(child; params...)
        add_right_branch!(tree,child_branch)
    end
    return tree
end

(anglepoint(center::P, angle, radius) where {P}) = center + P(cos(angle),sin(angle))*radius

function excSynapse(foot::Point{2,F}, dist, radius; csegs=20, dir=:left) where F
    res = Point{2,F}[]
    if dir == :left
        arc = anglepoint.(Ref(foot+Point{2,F}(-dist,0)), LinRange(pi/4, pi*7/4, csegs), radius)
        push!(res, foot+Point{2,F}(0,sqrt(2)*radius/2), arc..., foot-Point{2,F}(0,sqrt(2)*radius/2))
    else
        arc = anglepoint.(Ref(foot+Point{2,F}( dist,0)), LinRange(3pi/4, -3pi/4, csegs), radius)
        push!(res, foot+Point{2,F}(0,sqrt(2)*radius/2), arc..., foot-Point{2,F}(0,sqrt(2)*radius/2))
    end
    return res
end

function drawDendrite(bottom::Point{2,F}, width, height, cr; csegs=20, soma=false) where F
    cr = min(width/2-eps(F), height/2-eps(F), cr)

    # inner corners
    ictl = bottom .+ Point2(cr-width/2, -cr+height)
    ictr = bottom .+ Point2(-cr+width/2, -cr+height)
    icbl = bottom .+ Point2(cr-width/2, cr)
    icbr = bottom .+ Point2(-cr+width/2, cr)

    cstr = anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs), cr)
    cstl = anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs), cr)

    res = Point{2,F}[]
    append!(res, cstr, cstl)
    if soma
        csb = anglepoint.(Ref(bottom), LinRange(3pi/4, 9pi/4, 3*csegs), sqrt(2)*width/2)
        append!(res, csb)
        return res
    else
        csbl = anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs), cr)
        csbr = anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs), cr)
        append!(res, csbl, csbr)
        return res
    end
end

@recipe(NeuronPlot, neuron, offset) do scene
    Attributes(
        portside=:both,
        route_to_bottom=true,
        dend_width=0.1f0,
        syn_radius=0.025f0,
        dend_margin=0.1f0,
        row_margin=0.1f0,
        syn_margin=0.1f0,
        syn_stem=0.03f0,
        connector_kwargs=Dict(:linewidth=>5, :color=>:black),
        inputs_kwargs=Dict(:linewidth=>3, :color=>:black),
        dendrite_kwargs=Dict(:strokewidth=>2, :color=>:silver, :strokecolor=>:black),
        spine_kwargs=Dict(:strokewidth=>2, :strokecolor=>:black),
        csegs=20
    )
end

function Makie.plot!(neuronplot::NeuronPlot)
    F = Float32
    neuron = neuronplot[1]
    offset = neuronplot[2]

    (;
        portside, route_to_bottom, dend_width, syn_radius, dend_margin,
        row_margin, syn_margin, syn_stem, connector_kwargs, inputs_kwargs,
        dendrite_kwargs, spine_kwargs, csegs
    ) = neuronplot.attributes

    # prepare outputs
    all_ports = Dict{Symbol,Point{2,F}}()

    #hidedecorations!(ax)
    #neuronplot.aspect[]=DataAspect()

    connector_lines = Observable(Point{2,F}[])
    dendrites = Observable(Dict{Symbol,Vector{Point{2,F}}}())
    spines = Observable(Dict{Tuple{Symbol,Symbol},Vector{Point{2,F}}}())
    inputs = Observable(Dict{Symbol,Vector{Point{2,F}}}())
    y_extent = Observable(zeros(F,2))
    x_extent = Observable(zeros(F,2))

    onany(neuron,
        portside, route_to_bottom, dend_width, syn_radius, dend_margin,
        row_margin, syn_margin, syn_stem, csegs
    ) do neuron,
        portside, route_to_bottom, dend_width, syn_radius, dend_margin,
        row_margin, syn_margin, syn_stem, csegs

        empty!(connector_lines[])
        empty!(dendrites[])
        empty!(spines[])
        empty!(inputs[])

        # do initial layout to figure out horizontal placement
        tree = layout(neuron.dendrite; pos=Point{2,F}(zero(F),√(2)*dend_width/2), height=_->one(F), width=_->dend_width, x_margin=dend_margin, y_margin=row_margin)

        all_inputs = unique([
            flatten(node.node.exc_inputs for row in tree.rows for clique in row for node in clique)...,
            flatten(node.node.inh_inputs for row in tree.rows for clique in row for node in clique)...
        ])
        all_input_ys = Dict(name=>F[] for name in all_inputs)
        merge!(inputs[], Dict(name=>Point{2,F}[] for name in all_inputs))

        # figure out on which side to draw the symbol wires
        input_side = if portside == :left
            Dict(name => :left for name in all_inputs)
        elseif portside == :right
            Dict(name => :right for name in all_inputs)
        elseif portside == :both
            res = Dict{Symbol,Symbol}()
            # compute center of mass of synapse x-positions for each input symbol
            center_of_mass = Dict{Symbol,F}()
            for row in tree.rows
                for clique in row
                    for node in clique
                        for inp in node.node.exc_inputs
                            center_of_mass[inp] = get(center_of_mass, inp, zero(F)) + node.position[][1]
                        end
                        for inp in node.node.inh_inputs
                            center_of_mass[inp] = get(center_of_mass, inp, zero(F)) + node.position[][1]
                        end
                    end
                end
            end

            # assign side based on center of mass
            for name in all_inputs
                res[name] =  center_of_mass[name] <= 0 ? :left : :right
            end
            res
        else
            ArgumentError("'portside' must be either ':left', ':right' or ':both', not '$(portside)'.")
        end

        isleft(name) = input_side[name]==:left
        all_inputs_left = filter(isleft, all_inputs)
        all_inputs_right = filter(!isleft, all_inputs)

        # figure out which inputs can coexist on the left and right side
        row_inputs_left_padded = Vector{Union{Nothing,Symbol}}[]
        row_inputs_right_padded = Vector{Union{Nothing,Symbol}}[]
        for row in tree.rows
            from_left = Dict{Symbol,F}()
            from_right = Dict{Symbol,F}()
            # find rightmost and leftmost extent of left and right inputs, respectively
            for clique in row
                for node in clique
                    node_x = node.position[][1]
                    for input in node.node.exc_inputs
                        if isleft(input)
                            from_left[input] = max(node_x, get(from_left, input, typemin(F)))
                        else
                            from_right[input] = min(node_x, get(from_right, input, typemax(F)))
                        end
                    end
                    for input in node.node.inh_inputs
                        if isleft(input)
                            from_left[input] = max(node_x,get(from_left,input,typemin(F)))
                        else
                            from_right[input] = min(node_x,get(from_right,input,typemax(F)))
                        end
                    end
                end
            end

            # sort left and right inputs to conform with the overall order
            row_inputs = unique([
                flatten(node.node.exc_inputs for clique in row for node in clique)...,
                flatten(node.node.inh_inputs for clique in row for node in clique)...
            ])
            row_inputs_left = filter(isleft, row_inputs)
            row_inputs_right = filter(!isleft, row_inputs)
            order_left = sortperm(collect(indexin(row_inputs_left, all_inputs_left)))
            order_right = sortperm(collect(indexin(row_inputs_right, all_inputs_right)))
            permute!(row_inputs_left, order_left)
            permute!(row_inputs_right, order_right)

            # greedy algorithm: alternatingly place left and right if possible
            left=true
            order = NTuple{2,Union{Nothing,Symbol}}[]
            while length(row_inputs_left) + length(row_inputs_right) > 0
                next_left = isempty(row_inputs_left) ? nothing  : first(row_inputs_left)
                next_right = isempty(row_inputs_right) ? nothing  : first(row_inputs_right)

                if !isnothing(next_left) && !isnothing(next_right) && from_left[next_left] + syn_radius + syn_stem ≤ from_right[next_right]
                    # both inputs fit
                    push!(order,(next_left,next_right))
                    popfirst!(row_inputs_left)
                    popfirst!(row_inputs_right)
                elseif left && !isnothing(next_left)
                    push!(order,(next_left,nothing))
                    popfirst!(row_inputs_left)
                else
                    push!(order,(nothing,next_right))
                    popfirst!(row_inputs_right)
                end
                left != left
            end
            rl,rr = zip(order...)
            push!(row_inputs_left_padded,collect(rl))
            push!(row_inputs_right_padded,collect(rr))
        end

        # store for each node how tall the corresponding row is
        node_height = Dict(
            node.node.name => (length(row_inputs_left_padded[i])+1)*syn_margin
            for (i,row) in enumerate(tree.rows) for clique in row for node in clique
        )

        # re-layout with correct assigment of synapses to left and right side
        tree = layout(neuron.dendrite; pos=Point{2,F}(zero(F),√(2)*dend_width/2), height=node->node_height[node.name], width=x->dend_width, x_margin=dend_margin, y_margin=row_margin)

        # if we route the signals all the way to the bottom, they should be staggered
        if route_to_bottom
            input_x = Dict(
                (name=>minimum(tree.x_min) - dend_margin - syn_margin/2 * i for (i,name) in enumerate(all_inputs_left))...,
                (name=>maximum(tree.x_max) + dend_margin + syn_margin/2 * i for (i,name) in enumerate(all_inputs_right))...
            )
        else
            input_x = Dict(
                (name=>minimum(tree.x_min) - dend_margin for name in all_inputs_left)...,
                (name=>maximum(tree.x_max) + dend_margin for name in all_inputs_right)...
            )
        end

        x_extent[] .= extrema([collect(values(input_x)); tree.x_min; tree.x_max; zero(F)])
        y_extent[] .= (zero(F), tree.y_max[end])
        notify(x_extent)
        notify(y_extent)

        # compute, for each input, which other inputs (on the same side) it's vertical connector intersects with
        input_intersections = Dict(name=>F[] for name in all_inputs)

        # iterate all rows
        for (i,row) in enumerate(tree.rows)
            # compute y-position for all inputs
            input_y = Dict(
                (v=>k*syn_margin + tree.y_min[i] + syn_margin/2 for (k,v) in enumerate(row_inputs_left_padded[i]) if !isnothing(v))...,
                (v=>k*syn_margin + tree.y_min[i] + syn_margin/2 for (k,v) in enumerate(row_inputs_right_padded[i]) if !isnothing(v))...
            )
            mergewith!(union, all_input_ys, input_y)

            # iterate all cliques in the row
            for clique in row
                if isempty(clique)
                    continue
                end

                # target point
                pt5 = (first(clique).position[] + last(clique).position[])/2 - Point{2,F}(zero(F),tree.y_margin[])
                # joint
                pt4 = pt5 + Point{2,F}(zero(F),tree.y_margin[]/3)

                # iterate all nodes in the clique
                for (j, node) in enumerate(clique)
                    # start points
                    pt1 = node.position[]
                    node_x = pt1[1]

                    # knee points
                    pt2 = pt1 - Point{2,F}(zero(F), tree.y_margin[] / 3)
                    # foot points
                    pt3 = (pt2 + pt5) / 2

                    # add leg (if not root of the tree)
                    if i > 1
                        push!(connector_lines[], pt1, pt2, pt3, pt4, Point{2,F}(NaN, NaN))
                    end

                    # sort synapses by y-position
                    order = sortperm([input_y[inp] for inp in [node.node.exc_inputs; node.node.inh_inputs]])
                    node_inputs = [node.node.exc_inputs; node.node.inh_inputs][order]

                    # draw synapses
                    for inp_name in node_inputs
                        foot = Point{2,F}(node_x + (isleft(inp_name) ? -dend_width / 2 : dend_width / 2), input_y[inp_name] - syn_margin / 2)

                        # add the horizontal connector
                        pushfirst!(inputs[][inp_name],
                            Point{2,F}(input_x[inp_name], input_y[inp_name]),
                            foot + Point{2,F}((isleft(inp_name) ? -1 : 1) * (syn_stem + syn_margin / 2), syn_margin / 2),
                            foot + Point{2,F}((isleft(inp_name) ? -1 : 1) * (syn_stem + syn_radius * sqrt(2) / 2), syn_radius * sqrt(2) / 2),
                            Point{2,F}(NaN, NaN)
                        )

                        spines[][(node.node.name, inp_name)] = excSynapse(foot, syn_stem, syn_radius; csegs, dir=input_side[inp_name])
                    end

                    dendrites[][node.node.name] = drawDendrite(pt1, dend_width, tree.y_max[i] - tree.y_min[i], F(Inf); csegs, soma=i == 1)
                end

                if i>1
                    # add vertical line
                    push!(connector_lines[], pt4, pt5, Point{2,F}(NaN,NaN))
                end
            end
        end

        # complete input axons
        if route_to_bottom
            potential_intersections = F[]
            for inp_name in reverse(all_inputs_left)
                input_intersections = sort(unique(filter(<(maximum(all_input_ys[inp_name])),potential_intersections)))
                arcs = [anglepoint.(Ref(Point{2,F}(input_x[inp_name], inter)), LinRange(3pi/2,pi/2,csegs), syn_margin/4) for inter in input_intersections]
                port = Point{2,F}(input_x[inp_name], 0)
                pushfirst!(inputs[][inp_name], port, flatten(arcs)...)
                all_ports[inp_name] = port

                append!(potential_intersections, all_input_ys[inp_name])
            end
            potential_intersections = F[]
            for inp_name in reverse(all_inputs_right)
                input_intersections = sort(unique(filter(<(maximum(all_input_ys[inp_name])),potential_intersections)))
                arcs = [anglepoint.(Ref(Point{2,F}(input_x[inp_name], inter)), LinRange(3pi/2,pi/2,csegs), syn_margin/4) for inter in input_intersections]
                port = Point{2,F}(input_x[inp_name], 0)
                pushfirst!(inputs[][inp_name], port, flatten(arcs)...)
                all_ports[inp_name] = port

                append!(potential_intersections, all_input_ys[inp_name])
            end
        end
    end

    # force update
    neuron[] = neuron[]

    neuronplot[:x_extent] = x_extent
    neuronplot[:y_extent] = y_extent

    #draw clique connectors
    lines!(neuronplot, @lift($(connector_lines).+Ref($(offset))); connector_kwargs[]...)

    # draw spines
    if !isempty(spines[])
        col = pop!(spine_kwargs[], :color, :gray)
        neuronplot[:spines_color] = @lift(Dict(name=>col for name in keys($(spines))))
        neuronplot[:spines]=poly!(neuronplot, @lift(Vector{Point{2,F}}[val.+ Ref($(offset)) for val in values($(spines))]); spine_kwargs[]...,
            color=@lift([$(neuronplot[:spines_color])[key] for key in keys($(spines))]), spine_kwargs[]...)
        neuronplot[:spine_ids]=@lift(keys($(spines)))
    else
        neuronplot[:spines_color]   = Nothing[]
        neuronplot[:spines]         = Nothing[]
        neuronplot[:spine_ids]      = Nothing[]
    end

    # draw branches
    col = pop!(dendrite_kwargs[], :color, :gray)
    neuronplot[:dendrites_color] = @lift(Dict(name=>col for name in keys($(dendrites))))
    neuronplot[:dendrites]=poly!(neuronplot, @lift([val.+ Ref($(offset)) for val in values($(dendrites))]);
        color=@lift([$(neuronplot[:dendrites_color])[key] for key in keys($(dendrites))]), dendrite_kwargs[]...)
    neuronplot[:dendrite_ids]=@lift(keys($(dendrites)))

    #draw incoming axons
    if !isempty(inputs[])
        col = pop!(inputs_kwargs[], :color, :gray)
        neuronplot[:inputs_color] = @lift(Dict(name=>col for name in keys($(inputs))))
        neuronplot[:axons]=lines!(neuronplot, @lift(Point{2,F}[([val.+ Ref($offset); [Point{2,F}(NaN,NaN)]] for val in values($inputs))...;]);
            color=@lift([(fill($(neuronplot[:inputs_color])[key], length(val)+1) for (key,val) in pairs($(inputs)))...;]), inputs_kwargs[]...)
        neuronplot[:inputs_ids]=@lift(collect(zip(keys($(inputs)),length.(values($(inputs))))))
        neuronplot[:inputs_ports] = @lift(Dict(key=>first(value)+$(offset) for (key, value) in pairs($(inputs))))
    else
        neuronplot[:inputs_color]   = Nothing[]
        neuronplot[:axons]          = Nothing[]
        neuronplot[:inputs_ids]     = Nothing[]
        neuronplot[:inputs_ports]   = Nothing[]
    end

    neuronplot[:soma] = offset

    return neuronplot
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

# service methods for TypeDB
const _TypeDB_methods = MethodDescriptor[
        MethodDescriptor("databases_contains", 1, CoreDatabaseManager_Contains_Req, CoreDatabaseManager_Contains_Res),
        MethodDescriptor("databases_create", 2, CoreDatabaseManager_Create_Req, CoreDatabaseManager_Create_Res),
        MethodDescriptor("databases_all", 3, CoreDatabaseManager_All_Req, CoreDatabaseManager_All_Res),
        MethodDescriptor("database_schema", 4, CoreDatabase_Schema_Req, CoreDatabase_Schema_Res),
        MethodDescriptor("database_delete", 5, CoreDatabase_Delete_Req, CoreDatabase_Delete_Res),
        MethodDescriptor("session_open", 6, Session_Open_Req, Session_Open_Res),
        MethodDescriptor("session_close", 7, Session_Close_Req, Session_Close_Res),
        MethodDescriptor("session_pulse", 8, Session_Pulse_Req, Session_Pulse_Res),
        MethodDescriptor("transaction", 9, Channel{Transaction_Client}, Channel{Transaction_Server})
    ] # const _TypeDB_methods
const _TypeDB_desc = ServiceDescriptor("typedb.protocol.TypeDB", 1, _TypeDB_methods)

TypeDB(impl::Module) = ProtoService(_TypeDB_desc, impl)

mutable struct TypeDBStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    TypeDBStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_TypeDB_desc, channel))
end # mutable struct TypeDBStub

mutable struct TypeDBBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    TypeDBBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_TypeDB_desc, channel))
end # mutable struct TypeDBBlockingStub

databases_contains(stub::TypeDBStub, controller::ProtoRpcController, inp::CoreDatabaseManager_Contains_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[1], controller, inp, done)
databases_contains(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::CoreDatabaseManager_Contains_Req) = call_method(stub.impl, _TypeDB_methods[1], controller, inp)

databases_create(stub::TypeDBStub, controller::ProtoRpcController, inp::CoreDatabaseManager_Create_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[2], controller, inp, done)
databases_create(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::CoreDatabaseManager_Create_Req) = call_method(stub.impl, _TypeDB_methods[2], controller, inp)

databases_all(stub::TypeDBStub, controller::ProtoRpcController, inp::CoreDatabaseManager_All_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[3], controller, inp, done)
databases_all(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::CoreDatabaseManager_All_Req) = call_method(stub.impl, _TypeDB_methods[3], controller, inp)

database_schema(stub::TypeDBStub, controller::ProtoRpcController, inp::CoreDatabase_Schema_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[4], controller, inp, done)
database_schema(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::CoreDatabase_Schema_Req) = call_method(stub.impl, _TypeDB_methods[4], controller, inp)

database_delete(stub::TypeDBStub, controller::ProtoRpcController, inp::CoreDatabase_Delete_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[5], controller, inp, done)
database_delete(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::CoreDatabase_Delete_Req) = call_method(stub.impl, _TypeDB_methods[5], controller, inp)

session_open(stub::TypeDBStub, controller::ProtoRpcController, inp::Session_Open_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[6], controller, inp, done)
session_open(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::Session_Open_Req) = call_method(stub.impl, _TypeDB_methods[6], controller, inp)

session_close(stub::TypeDBStub, controller::ProtoRpcController, inp::Session_Close_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[7], controller, inp, done)
session_close(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::Session_Close_Req) = call_method(stub.impl, _TypeDB_methods[7], controller, inp)

session_pulse(stub::TypeDBStub, controller::ProtoRpcController, inp::Session_Pulse_Req, done::Function) = call_method(stub.impl, _TypeDB_methods[8], controller, inp, done)
session_pulse(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::Session_Pulse_Req) = call_method(stub.impl, _TypeDB_methods[8], controller, inp)

transaction(stub::TypeDBStub, controller::ProtoRpcController, inp::Channel{Transaction_Client}, done::Function) = call_method(stub.impl, _TypeDB_methods[9], controller, inp, done)
transaction(stub::TypeDBBlockingStub, controller::ProtoRpcController, inp::Channel{Transaction_Client}) = call_method(stub.impl, _TypeDB_methods[9], controller, inp)

export TypeDB, TypeDBStub, TypeDBBlockingStub, databases_contains, databases_create, databases_all, database_schema, database_delete, session_open, session_close, session_pulse, transaction


text:
module RRRMCTest

using RRRMC
using Test
using Random

function gen_timeout_hook(t = 1.0)
    t += time()
    return (args...) -> (time() ≤ t)
end

function checkenergy_hook(it, X, C, acc, E)
    @test E ≈ RRRMC.energy(X, C) atol=1e-11
    return true
end

function checkenergy_hook_EO(it, X, C, E, Emin)
    @test E ≈ RRRMC.energy(X, C) atol=1e-11
    return true
end

function test()

    Random.seed!(8426732438942) # get reproducible results...

    graphs = [
        RRRMC.GraphTwoSpin(),
        RRRMC.GraphThreeSpin(),
        RRRMC.GraphFields(10),
        RRRMC.GraphFields(10, (-1,1)),
        RRRMC.GraphFields(10, (-1.5,0.5)),
        RRRMC.GraphFieldsNormalDiscretized(10, (-1,0,1)),
        RRRMC.GraphFieldsNormalDiscretized(10, (-1.5,0.0,1.5)),
        RRRMC.GraphIsing1D(10),
        RRRMC.GraphPSpin3(39, 5),
        RRRMC.GraphRRG(10, 3),
        RRRMC.GraphRRG(10, 3, (-1,0,1)),
        RRRMC.GraphRRG(10, 3, (-1.0,0.0,1.0)),
        RRRMC.GraphRRG(10, 3, (-1.0,0,1)),
        RRRMC.GraphRRG(10, 3, (-1//1,0//1,1//1)),
        RRRMC.GraphRRGNormalDiscretized(10, 3, (-1,0,1)),
        RRRMC.GraphRRGNormalDiscretized(10, 3, (-1.0,0.0,1.0)),
        RRRMC.GraphRRGNormalDiscretized(10, 3, (-1.0,0,1)),
        RRRMC.GraphRRGNormalDiscretized(10, 3, (-1//1,0//1,1//1)),
        RRRMC.GraphRRGNormal(10, 3),
        RRRMC.GraphEA(2, 3),
        RRRMC.GraphEA(2, 3, (-1,0,1)),
        RRRMC.GraphEA(2, 3, (-1.0,0.0,1.0)),
        RRRMC.GraphEA(2, 3, (-1.0,0,1)),
        RRRMC.GraphEA(2, 3, (-1//1,0//1,1//1)),
        RRRMC.GraphEANormalDiscretized(2, 3, (-1,0,1)),
        RRRMC.GraphEANormalDiscretized(2, 3, (-1.0,0.0,1.0)),
        RRRMC.GraphEANormalDiscretized(2, 3, (-1.0,0,1)),
        RRRMC.GraphEANormalDiscretized(2, 3, (-1//1,0//1,1//1)),
        RRRMC.GraphEANormal(2, 3),
        RRRMC.GraphEA(3, 2),
        RRRMC.GraphEA(3, 2, (-1,0,1)),
        RRRMC.GraphEA(3, 2, (-1.0,0.0,1.0)),
        RRRMC.GraphEA(3, 2, (-1.0,0,1)),
        RRRMC.GraphEA(3, 2, (-1//1,0//1,1//1)),
        RRRMC.GraphEANormalDiscretized(3, 2, (-1,0,1)),
        RRRMC.GraphEANormalDiscretized(3, 2, (-1.0,0.0,1.0)),
        RRRMC.GraphEANormalDiscretized(3, 2, (-1.0,0,1)),
        RRRMC.GraphEANormalDiscretized(3, 2, (-1//1,0//1,1//1)),
        RRRMC.GraphEANormal(3, 2),
        RRRMC.GraphSK(10),
        RRRMC.GraphSKNormal(10),
        RRRMC.GraphSAT(10, 3, 4.2),
        RRRMC.GraphPercLinear(101, 30),
        RRRMC.GraphPercStep(101, 30),
        RRRMC.GraphPercXEntr(101, 30, 1.5),
        RRRMC.GraphCommStep(25, 5, 30),
        RRRMC.GraphCommStep(21, 5, 30, fc=true),
        RRRMC.GraphCommReLU(24, 6, 30),
        RRRMC.GraphCommReLU(22, 6, 30, fc=true),
        RRRMC.GraphCommQu(24, 6, 30),
        RRRMC.GraphCommQu(22, 6, 30, fc=true),
        RRRMC.GraphQuant(10, 8, 0.5, 2.0, RRRMC.GraphEmpty, 10),
        RRRMC.GraphQuant(10, 8, 0.5, 2.0, RRRMC.GraphSK, RRRMC.SK.gen_J(10)),
        RRRMC.GraphQuant(10, 8, 0.5, 2.0, RRRMC.GraphSKNormal, RRRMC.SK.gen_J_gauss(10)),
        RRRMC.GraphQuant(3, 8, 0.5, 2.0, RRRMC.GraphThreeSpin),
        RRRMC.GraphQPercLinearT(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphQPercStepT(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphQCommStepT(25, 5, 30, 5, 0.5, 2.0),
        RRRMC.GraphQCommStepT(21, 5, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphQCommReLUT(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphQCommReLUT(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphQCommQuT(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphQCommQuT(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphRobustEnsemble(10, 8, 1.5, 2.0, RRRMC.GraphEmpty, 10),
        RRRMC.GraphRobustEnsemble(10, 8, 1.5, 2.0, RRRMC.GraphSK, RRRMC.SK.gen_J(10)),
        RRRMC.GraphRobustEnsemble(10, 8, 1.5, 2.0, RRRMC.GraphSKNormal, RRRMC.SK.gen_J_gauss(10)),
        RRRMC.GraphRobustEnsemble(3, 8, 1.5, 2.0, RRRMC.GraphThreeSpin),
        RRRMC.GraphPercLinearRE(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphPercStepRE(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommStepRE(25, 5, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommStepRE(21, 5, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphCommReLURE(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommReLURE(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphCommQuRE(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommQuRE(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphLocalEntropy(10, 8, 1.5, 2.0, RRRMC.GraphEmpty, 10),
        RRRMC.GraphLocalEntropy(10, 8, 1.5, 2.0, RRRMC.GraphSKNormal, RRRMC.SK.gen_J_gauss(10)),
        RRRMC.GraphLocalEntropy(3, 8, 1.5, 2.0, RRRMC.GraphThreeSpin),
        RRRMC.GraphRobustEnsemble(20, 4, 1.5, 2.0, RRRMC.GraphQuant, 5, 4, 0.5, 2.0, RRRMC.GraphSK, RRRMC.SK.gen_J(5)),
        RRRMC.GraphPercLinearLE(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphPercStepLE(101, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommStepLE(25, 5, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommStepLE(21, 5, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphCommReLULE(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommReLULE(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphCommQuLE(24, 6, 30, 5, 0.5, 2.0),
        RRRMC.GraphCommQuLE(22, 6, 30, 5, 0.5, 2.0, fc=true),
        RRRMC.GraphMixed(RRRMC.GraphRRGNormal(18, 3),
                         RRRMC.GraphPSpin3(18, 5)),
        RRRMC.GraphMixed(RRRMC.GraphRRG(10, 3),
                         RRRMC.GraphSK(10)),
        RRRMC.GraphMixed(RRRMC.GraphSK(31),
                         RRRMC.GraphPercStep(31, 30),
                         RRRMC.GraphSAT(31, 3, 4.2)),
        RRRMC.GraphAddFields(randn(10), RRRMC.GraphRRGNormal(10, 3)),
        RRRMC.GraphAddSubFields(randn(10), RRRMC.GraphRRGNormal(10, 3)),
       ]

    β = 2.0
    iters = 10_000
    st = 100
    τ = 1.3
    samples = iters ÷ st
    quiet = true

    function abbrev(X, n)
        s = string(X)
        l = length(s)
        return l > n ? chop(s, tail = l-n) : rpad(s, n)
    end

    @testset "$(abbrev(X, 60))" for X in graphs
        # @show X
        E, C = standardMC(X, β, iters, step=st, quiet=quiet)
        E, C = standardMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)
        E, C = standardMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=gen_timeout_hook())

        E, C = bklMC(X, β, iters, step=st, quiet=quiet)
        E, C = bklMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)
        E, C = bklMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=gen_timeout_hook())

        E, C = wtmMC(X, β, samples, step=Float64(st), quiet=quiet)
        E, C = wtmMC(X, β, samples, step=Float64(st), quiet=quiet, C0=C, hook=checkenergy_hook)
        E, C = wtmMC(X, β, samples, step=Float64(st), quiet=quiet, C0=C, hook=gen_timeout_hook())

        E, C = rrrMC(X, β, iters, step=st, quiet=quiet)
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=gen_timeout_hook())
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, staged_thr=0.0, hook=checkenergy_hook)
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, staged_thr=0.0)
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, staged_thr=1.0, hook=checkenergy_hook)
        E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, staged_thr=1.0)

        C, Emin, Cmin, itmin = extremal_opt(X, τ, iters, step=st, quiet=quiet)
        C, Emin, Cmin, itmin = extremal_opt(X, τ, iters, step=st, quiet=quiet, C0=C)
        C, Emin, Cmin, itmin = extremal_opt(X, τ, iters, step=st, quiet=quiet, hook=checkenergy_hook_EO)
        C, Emin, Cmin, itmin = extremal_opt(X, τ, iters, step=st, quiet=quiet, hook=gen_timeout_hook())

        if X isa RRRMC.DoubleGraph
            X0 = RRRMC.inner_graph(X)
            E, C = bklMC(X0, β, iters, step=st, quiet=quiet)
            E, C = bklMC(X0, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)

            E, C = rrrMC(X, β, iters, step=st, quiet=quiet)
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, hook=gen_timeout_hook())
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, staged_thr=0.0, hook=checkenergy_hook)
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, staged_thr=0.0)
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, staged_thr=1.0, hook=checkenergy_hook)
            E, C = rrrMC(X, β, iters, step=st, quiet=quiet, C0=C, staged_thr=1.0)

            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet)
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, C0=C, hook=checkenergy_hook)
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, C0=C, hook=gen_timeout_hook())
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, staged_thr=0.0, hook=checkenergy_hook)
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, C0=C, staged_thr=0.0)
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, staged_thr=1.0, hook=checkenergy_hook)
            E, C = rrrMC(X0, β, iters, step=st, quiet=quiet, C0=C, staged_thr=1.0)

            C, Emin, Cmin, itmin = extremal_opt(X0, τ, iters, step=st, quiet=quiet)
            C, Emin, Cmin, itmin = extremal_opt(X0, τ, iters, step=st, quiet=quiet, C0=C)
            C, Emin, Cmin, itmin = extremal_opt(X0, τ, iters, step=st, quiet=quiet, hook=checkenergy_hook_EO)
            C, Emin, Cmin, itmin = extremal_opt(X0, τ, iters, step=st, quiet=quiet, hook=gen_timeout_hook())
        end
    end
end

test()

end # module


text:
# ---
# title: 1391. Check if There is a Valid Path in a Grid
# id: problem1391
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/check-if-there-is-a-valid-path-in-a-grid/description/>
# hidden: true
# ---
#
# Given a _m_ x _n_ `grid`. Each cell of the `grid` represents a street. The
# street of `grid[i][j]` can be:
#
#   * **1** which means a street connecting the left cell and the right cell.
#   * **2** which means a street connecting the upper cell and the lower cell.
#   * **3**  which means a street connecting the left cell and the lower cell.
#   * **4** which means a street connecting the right cell and the lower cell.
#   * **5** which means a street connecting the left cell and the upper cell.
#   * **6** which means a street connecting the right cell and the upper cell.
#
# ![](https://assets.leetcode.com/uploads/2020/03/05/main.png)
#
# You will initially start at the street of the upper-left cell `(0,0)`. A valid
# path in the grid is a path which starts from the upper left cell `(0,0)` and
# ends at the bottom-right cell `(m - 1, n - 1)`. **The path should only follow
# the streets**.
#
# **Notice** that you are **not allowed** to change any street.
#
# Return _true_  if there is a valid path in the grid or _false_ otherwise.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/03/05/e1.png)
#
#
#
#     Input: grid = [[2,4,3],[6,5,2]]
#     Output: true
#     Explanation: As shown you can start at cell (0, 0) and visit all the cells of the grid to reach (m - 1, n - 1).
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/03/05/e2.png)
#
#
#
#     Input: grid = [[1,2,1],[1,2,1]]
#     Output: false
#     Explanation: As shown you the street at cell (0, 0) is not connected with any street of any other cell and you will get stuck at cell (0, 0)
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[1,1,2]]
#     Output: false
#     Explanation: You will get stuck at cell (0, 1) and you cannot reach cell (0, 2).
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[1,1,1,1,1,1,3]]
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: grid = [[2],[2],[2],[2],[2],[2],[6]]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 300`
#   * `1 <= grid[i][j] <= 6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
abstract type AbstractDOI <: ShortCode end
struct EmDOI <: AbstractDOI
    doi
    highlight # Author to highlight when displaying
end
struct DOI <: AbstractDOI
    doi
end
struct ShortDOI <: AbstractDOI
    shortdoi
    ShortDOI(doi::String) = length(doi) > 10 ? new(shortdoi(DOI(doi))) : new(doi)
end
EmDOI(doi::String) = EmDOI(doi, "")
ShortDOI(doi::AbstractDOI) = ShortDOI(shortdoi(doi))

function Base.getproperty(obj::AbstractDOI, sym::Symbol)
    if sym == :doi
        return getdoi(obj)                                  # different for ShortDOI
    elseif sym in [:author, :title, :source_title, :page]   # string types
        return fetch_metadata(obj)[sym]
    elseif sym in [:year, :volume, :citation_count, :issue] # integer types
        return parse(Int, fetch_metadata(obj)[sym])
    elseif sym == :reference                                # DOI type
        return split(fetch_metadata(obj)[sym], ";") .|> x->DOI(replace(x, " "=>""))
    else # fallback to getfield
        return getfield(obj, sym)
    end
end
getdoi(obj::AbstractDOI) = getfield(obj, :doi)
getdoi(obj::ShortDOI) = expand(obj)

function Base.show(io::IO, ::MIME"text/plain", doi::AbstractDOI)
    print(io, join((doi.author, doi.title, string(doi.year), shortdoi(doi)), " "))
end

function Base.show(io::IO, ::MIME"text/html", doi::AbstractDOI)
    print(io, "
$(emph_author(doi)) $(doi.title), $(doi.source_title) ($(doi.year))
     $(shortdoi(doi)), cited by $(doi.citation_count)
")
end

function Base.show(io::IO, ::MIME"text/html", dois::Array{T} where T<:AbstractDOI)
    print(io, "
")
    for doi in dois
        print(io, "
$(emph_author(doi))) $(doi.title), $(doi.source_title) ($(doi.year))
        $(shortdoi(doi)), cited by $(doi.citation_count)
")
    end
    print(io, "
")
end

@memoize function fetch_metadata(doi::AbstractDOI) fetch_metadata(doi.doi) end
@memoize function fetch_metadata(doi)
    r = HTTP.get("https://opencitations.net/index/api/v1/metadata/$(doi)")
    rj = JSON3.read(r.body)
    return rj[1]
end

@memoize shortdoi(doi::AbstractDOI) = fetch_shortdoi(doi).ShortDOI

@memoize function fetch_shortdoi(doi::AbstractDOI)
    fetch_shortdoi(doi.doi)
end

@memoize function fetch_shortdoi(doi::String)
    r = HTTP.get("https://shortdoi.org/$(doi)?format=json")
    rj = JSON3.read(r.body)
end
"""
    expand(doi::ShortDOI)

Get full DOI from doi.org
"""
@memoize function expand(doi::ShortDOI)
    r = HTTP.get("https://doi.org/api/handles/$(doi.shortdoi)")
    return JSON3.read(r.body)[:values][2][:data][:value]
end

"""
    emph_author(authors, author="", em="b")

Add emphasis to selected author display (e.g. for CV use)
"""
function emph_author(doi::AbstractDOI)
    emph_author(doi.author)
end
function emph_author(doi::EmDOI)
    emph_author(doi.author, doi.highlight)
end
function emph_author(authors, author="", em="b")
    if length(author) > 2
        names = split(author)
        lnfirst = names[end] * ", " * names[1]
        short = names[end] * ", " * names[1][1] * "."
        return replace(replace(replace(authors, author => "<$em>"*author*"</$em>"),short=>"<$em>"*short*"</$em>"),lnfirst => "<$em>" * lnfirst * "</$em>")
    else
        return authors
    end
end


text:
using CUDA
using NESSie
using ImplicitArrays
using Test

using CuNESSie
using CuNESSie: CuPosition, CuPositionVector, CuTriangle, CuTriangleVector, PositionVector,
    TriangleVector, laplacepot_single, laplacepot_double, regularyukawapot_single,
    regularyukawapot_double

@testset "CuNESSie.jl" begin
    @testset "device" begin
        @testset "math.jl" begin include("device/math.jl") end
        @testset "laplace.jl" begin include("device/laplace.jl") end
        @testset "yukawa.jl" begin include("device/yukawa.jl") end
    end
end


text:
using Formatting: printfmt
using LinearAlgebra: Diagonal, Hermitian, eigen, norm, tr, diag, dot

using PyCall: pyimport
pyscf = pyimport("pyscf")
mp = pyimport("pyscf.mp")   # Had to import mp alone ??!
cc = pyimport("pyscf.cc") # Had to import mp alone ??!

import PySCF: scf_rhf, mymp2, short_mp2, ccsd

#########################################################################
# Benchmark own HF, MP₂, and CCSD against PySCF
function benchmark_h2o_sto3g()
   mol = pyscf.gto.Mole()
   mol.build(verbose = 0,
             basis = "sto-3g",
             #atom = pyscf_atom_from_xyz("../data/h2.xyz"),
             #atom = pyscf_atom_from_xyz("../data/acetaldehyde.xyz"),
             atom = """8 0  0.    0
                       1 0 -0.757 0.587
                       1 0  0.757 0.587""",
            )

   mf = pyscf.scf.RHF(mol)

   rhf = scf_rhf(mol,verbose=false)
   e_rhf_me, e, C = rhf.Emol, rhf.e, rhf.C
   e_rhf_pyscf    = mf.kernel()

   #mf.analyze(verbose=3,ncol=10, digits=9)

   println( "           E(HF)     ")
   println( "        -------------")
   printfmt(" ME:    {1:13.8f}  \n", e_rhf_me)
   printfmt(" PySCF: {1:13.8f}  \n", e_rhf_pyscf)

   eri   = mol.intor("cint2e_sph", aosym="s8")
   println()
   e_mp2_me1   = @time mymp2(e,C,eri,mol,alg=:ao2mo_noddy) # Naive loop  ao2mo with my MP₂
   e_mp2_me2   = @time mymp2(e,C,eri,mol,alg=:ao2mo_smart) # Smart loop  ao2mo with my MP₂
   e_mp2_me3   = @time short_mp2(e,C,eri,mol)              # @tensor     ao2mo with my MP₂
   e_mp2_pyscf = @time mp.MP2(mf).kernel()[1]              # pure PySCF MP₂ solution
   println()

   println( "           Ecorr(MP₂)        E(MP₂)   ")
   println( "          -----------   -------------")
   printfmt(" ME1:   {1:13.8f}   {2:13.8f}  \n", e_mp2_me1,  e_rhf_me   +e_mp2_me1)
   printfmt(" ME2:   {1:13.8f}   {2:13.8f}  \n", e_mp2_me2,  e_rhf_me   +e_mp2_me1)
   printfmt(" ME3:   {1:13.8f}   {2:13.8f}  \n", e_mp2_me3,  e_rhf_me   +e_mp2_me1)
   printfmt(" PySCF: {1:13.8f}   {2:13.8f}  \n", e_mp2_pyscf,e_rhf_pyscf+e_mp2_pyscf)
   println()

   @time e_ccsd_me, emp2 = ccsd(e,C,eri,mol,verbose=false)
   @time e_ccsd_pyscf    = cc.ccsd.CC(mf).kernel()[1]
   emp2 ≈ e_mp2_me3 || @warn "E_MP2 from CC is not accurate"
   println()

   println("            E(HF)          Ecorr(MP2)           E(MP2)         Ecorr(CCSD)  ")
   println("        --------------   --------------   ---------------   ---------------")
   printfmt(" ME:    {1:13.10f}    {2:13.10f}    {3:13.10f}    {4:13.10f}\n",
                    e_rhf_me,      e_mp2_me3,     e_rhf_me+e_mp2_me3,    e_ccsd_me)
   printfmt(" PySCF: {1:13.10f}    {2:13.10f}    {3:13.10f}    {4:13.10f}\n",
                     e_rhf_pyscf,  e_mp2_pyscf,e_rhf_pyscf+e_mp2_pyscf,e_ccsd_pyscf)
end

benchmark_h2o_sto3g()


text:
# Code to extract parameters from unlabeled parameter vector

# model assumes model parameters = [ Period, K, k, h, M0 ] for each planet followed by each RV offset and jitter
const num_param_per_planet = 5
const num_jitters = 0   # WARNING: Setting this to zero removes jitter from parameter list, setting it to one leads to errors
pl_offset(plid::Integer) = (plid-1)*num_param_per_planet
num_planets(theta::Vector) = floor(Int64,(length(theta)-num_jitters)//num_param_per_planet)  # WARNING: Assumes num_rvoffsets<num_param_per_planet
num_obs_offsets(theta::Vector) = length(theta)-num_planets(theta)*num_param_per_planet-num_jitters
obs_offset(theta::Vector,obsid::Integer = 1) = num_planets(theta)*num_param_per_planet+obsid
jitter_offset(theta::Vector,obsid::Integer = 1) = length(theta)-(num_jitters-1)

# constants defining unit system/scale for modified Jeffrys priors
const P0 = 1.0  # units of days
const K0 = 1.0  # units of m/s
const Jitter0 = 1.0  # units of m/s

set_period(theta::Vector, P; plid::Integer = 1) = theta[1+pl_offset(plid)] = log1p(P/P0)
set_amplitude(theta::Vector, K; plid::Integer = 1) = theta[2+pl_offset(plid)] = log1p(K/K0)
set_ecosw(theta::Vector, h; plid::Integer = 1) = theta[3+pl_offset(plid)] = h
set_esinw(theta::Vector, k; plid::Integer = 1) = theta[4+pl_offset(plid)] = k
set_w_plus_mean_anomaly_at_t0(theta::Vector, wpM0; plid::Integer = 1) = theta[5+pl_offset(plid)] = wpM0
set_rvoffset(theta::Vector, C; obsid::Integer = 1) = theta[obs_offset(theta,obsid)] = C
set_jitter(theta::Vector, jitter; obsid::Integer = 1) = theta[jitter_offset(theta,obsid)] = log1p(jitter/Jitter0)
function set_ew(theta::Vector, e, w; plid::Integer = 1)
  set_ecosw(theta,e*cos(w),plid=plid)
  set_esinw(theta,e*sin(w),plid=plid)
  return
end
function set_ewM0(theta::Vector, e, w, M0; plid::Integer = 1)
  set_ecosw(theta,e*cos(w),plid=plid)
  set_esinw(theta,e*sin(w),plid=plid)
  set_w_plus_mean_anomaly_at_t0(theta,w+M0,plid=plid)
  return
end
function set_PKewM0(theta::Vector, P, K, e, w, M0; plid::Integer = 1)
  set_period(theta,P,plid=plid)
  set_amplitude(theta,K,plid=plid)
  set_ecosw(theta,e*cos(w),plid=plid)
  set_esinw(theta,e*sin(w),plid=plid)
  set_w_plus_mean_anomaly_at_t0(theta,w+M0,plid=plid)
  return
end

extract_period(theta::Vector; plid::Integer = 1) = P0*(exp(theta[1+pl_offset(plid)])-1)
extract_amplitude(theta::Vector; plid::Integer = 1) = K0*(exp(theta[2+pl_offset(plid)])-1)
extract_ecosw(theta::Vector; plid::Integer = 1) = theta[3+pl_offset(plid)]
extract_esinw(theta::Vector; plid::Integer = 1) = theta[4+pl_offset(plid)]
extract_w(theta::Vector; plid::Integer = 1) = atan(extract_esinw(theta,plid), extract_ecosw(theta,plid))
extract_w_plus_mean_anomaly_at_t0(theta::Vector; plid::Integer = 1) = theta[5+pl_offset(plid)]
extract_mean_anomaly_at_t0(theta::Vector; plid::Integer = 1) = extract_w_plus_mean_anomaly_at_t0(theta,plid)-extract_w(theta,plid)
extract_rvoffset(theta::Vector; obsid::Integer = 1) = theta[obs_offset(theta,obsid)]
extract_jitter(theta::Vector; obsid::Integer = 1) = Jitter0*(exp(theta[jitter_offset(theta,obsid)])-1)
function extract_ewM0(theta::Vector; plid::Integer = 1)
  h = extract_ecosw(theta,plid=plid)
  k = extract_esinw(theta,plid=plid)
  wpM0 = extract_w_plus_mean_anomaly_at_t0(theta,plid=plid)
  e = sqrt(h*h+k*k)
  w = atan2(k,h)
  M0 = wpM0-w
  return (e,w,M0)
end

function extract_PKhkM(theta::Vector; plid::Integer = 1)
  local P,K,ecosw,esinw,M0
  P = extract_period(theta, plid=plid)
  K = extract_amplitude(theta, plid=plid)
  ecosw = extract_ecosw(theta, plid=plid)
  esinw = extract_esinw(theta, plid=plid)
  M0 = extract_w_plus_mean_anomaly_at_t0(theta, plid=plid) - atan2(esinw,ecosw)
  return (P,K,ecosw,esinw,M0)
end

function extract_PKewM(theta::Vector; plid::Integer = 1)
  local P,K,e,w,M0
  P = extract_period(theta, plid=plid)
  K = extract_amplitude(theta, plid=plid)
  (e,w,M0) = extract_ewM0(theta, plid=plid)
  return (P,K,e,w,M0)
end

# check if parameters are valid, i.e., e<1
function is_valid(p::Vector)
  local num_pl,j,k,ecc_sq
  num_pl = num_planets(p)
  @assert length(num_pl) >= 1
  for plid in 1:num_pl
    h = extract_ecosw(p,plid=plid)
    k = extract_esinw(p,plid=plid)
    ecc_sq = h*h+k*k
    if !(ecc_sq<1.0)
      return false
    end
  end
  return true
end

export num_planets, num_obs_offsets
export set_period, set_amplitude, set_ecosw, set_esinw, set_w_plus_mean_anomaly_at_t0, set_mean_anomaly_at_t0, set_rvoffset, set_jitter, set_ew, set_ewM0
export extract_period, extract_amplitude, extract_ecosw, extract_esinw, extract_w_plus_mean_anomaly_at_t0, extract_mean_anomaly_at_t0, extract_ewM0, extract_PKhkM, extract_PKewM, extract_rvoffset, extract_jitter
export extract_PKhkM, isvalid


text:
## Describes how dimensions of interpretation work

mutable struct Dimension
    startline::UInt
    current_line::UInt
    endline::UInt
    namespace::Namespace
    name::String
    busy::Bool
    isritual::Bool
end
Dimension(cl::I where I<:Integer, el::I where I<:Integer, ns::Namespace, name::String) = Dimension(cl, cl, el, ns, name, false, false)
Ritual(cl::I where I<:Integer, el::I where I<:Integer, ns::Namespace, name::String) = Dimension(cl, cl, el, ns, name, false, true)

function dimension_forward!(dim::Dimension)::Nothing
    try
        refline::Ref{UInt} = Ref{UInt}(dim.current_line)
        aetherlang_eval!(lines[dim.current_line], dim.namespace, refline)
        dim.current_line = refline[]+1
        if (dim.current_line > dim.endline && dim.isritual)
            dim.current_line = dim.startline;
        elseif dim.current_line > dim.endline || dim.current_line < dim.startline
            global tofinish
            push!(tofinish, dim.name)
        end
    catch e
        aetherlang_print_exception(e, lines[dim.current_line], dim.name) # also exits
    end
    dim.busy = false; # is set to true in Interpreter.jl
    nothing
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
# Rational Quadratic Isotropic Covariance Function

"""
    RQIso <: Isotropic{SqEuclidean}

Isotropic Rational Quadratic kernel (covariance)
```math
k(x,x') = σ²(1 + (x - x')ᵀ(x - x')/(2αℓ²))^{-α}
```
with length scale ``ℓ``, signal standard deviation ``σ``, and shape parameter ``α``.
"""
mutable struct RQIso <: Isotropic{SqEuclidean}
    "Squared length scale"
    ℓ2::Float64
    "Signal variance"
    σ2::Float64
    "Shape parameter"
    α::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        RQIso(ll:Float64, lσ::Float64, lα::Float64)

    Create `RQIso` with length scale `exp(ll)`, signal standard deviation `exp(lσ)`, and
    shape parameter `exp(lα)`.
    """
    RQIso(ll::Float64, lσ::Float64, lα::Float64) =
        new(exp(2 * ll), exp(2 * lσ), exp(lα), [])
end

function set_params!(rq::RQIso, hyp::VecF64)
    length(hyp) == 3 || throw(ArgumentError("Rational Quadratic function has three parameters"))
    rq.ℓ2, rq.σ2, rq.α = exp(2 * hyp[1]), exp(2 * hyp[2]), exp(hyp[3])
end

get_params(rq::RQIso) = Float64[log(rq.ℓ2) / 2, log(rq.σ2) / 2, log(rq.α)]
get_param_names(rq::RQIso) = [:ll, :lσ, :lα]
num_params(rq::RQIso) = 3

Statistics.cov(rq::RQIso, r::Float64) = rq.σ2 * (1 + r / (2 * rq.α * rq.ℓ2))^(-rq.α)

@inline dk_dll(rq::RQIso, r::Float64) =
    (s = r / rq.ℓ2; rq.σ2 * s * (1 + s / (2 * rq.α))^(-rq.α - 1)) # dK_d(log ℓ)dK_dℓ
@inline function dk_dlα(rq::RQIso, r::Float64)
    s = r / rq.ℓ2
    part = 1 + s / (2 * rq.α)
    rq.σ2 * part^(-rq.α) * (s / (2 * part) - rq.α * log(part))  # dK_d(log α)
end
@inline function dk_dθp(rq::RQIso, r::Float64, p::Int)
    if p==1
        return dk_dll(rq, r)
    elseif p==2
        return dk_dlσ(rq, r)
    elseif p==3
        return dk_dlα(rq, r)
    else
        return NaN
    end
end


text:
cDevice = Cint(0)
cChannel = Cint(0)
hzFreq = Cdouble(0)
szDeviceName = zeros(Cchar, 32)
szSN = zeros(Cchar, 32)
fIsInUse = Bool(0)
hdwf = HDWF(0)
szError = zeros(Cchar, 512)

# detect connected all supported devices
if FDwfEnum(enumfilterAll, Ref(cDevice)) == 0
    FDwfGetLastErrorMsg(szError)
    error("FDwfEnum: $szError")
end
# list information about each device
@info "Found $cDevice devices:"
for i in 0:cDevice-1
    # we use 0 based indexing
    FDwfEnumDeviceName(i, szDeviceName);
    FDwfEnumSN(i, szSN);
    println("Device: $(i+1) name: $szDeviceName $szSN")
    # before opening, check if the device isn�t already opened by other application, like: WaveForms
    FDwfEnumDeviceIsOpened(i, Ref(fIsInUse));
    if fIsInUse != 1
        if FDwfDeviceOpen(i, Ref(hdwf)) != 1
            FDwfGetLastErrorMsg(szError)
            println("FDwfDeviceOpen: $szError")
            break
        end
        FDwfAnalogInChannelCount(hdwf, Ref(cChannel))
        FDwfAnalogInFrequencyInfo(hdwf, C_NULL, Ref(hzFreq))
        println("number of analog input channels: $cChannel maximum freq.: $hzFreq Hz")
        FDwfDeviceClose(hdwf)
        hdwf = hdwfNone
    end
end
# before application exit make sure to close all opened devices by this process
FDwfDeviceCloseAll()


text:
# plot first 5000 time points (before stabilization)
include("../Model_scripts/Functions.jl");
# plot the number of infected & uninfected hosts
plot_population_numbers(
    hcat(labels...),
    "by=0.1\n \nNumber of infected and uninfected hosts",
    100,
    5000,
    "$figure_directory/supplemental/supp_1a.png"
)
# plot the average weighted R0
plot_reproductive_rate(6, 5000, "$figure_directory/supplemental/supp_1c.png")
# plot the average weighted V0
plot_vertical_reproductive_ratio(1, 5000, "$figure_directory/supplemental/supp_1e.png")
# plot the weighted average mortality rate through time
plot_mortality(1, 5000, "$figure_directory/supplemental/supp_1g.png")
# plot the evenness through time
plot_evenness(1, 5000, "$figure_directory/supplemental/supp_1i.png")


text:
@testset "dijkstra (dense)" begin
    begin
        g = AdjacencyMatrix(5)
        setedge!(g, 1, 2, 10)
        setedge!(g, 1, 3, 3)
        setedge!(g, 1, 4, 20)
        setedge!(g, 2, 4, 5)
        setedge!(g, 3, 2, 2)
        setedge!(g, 4, 5, 11)
        setedge!(g, 3, 5, 15)
    end
    @test dijkstra(g, 1) == [0, 5, 3, 10, 18]
end

@testset "dijkstra (sparse)" begin
    begin
        g = AdjacencyList(5)
        setedge!(g, 1, 2, 10)
        setedge!(g, 1, 3, 3)
        setedge!(g, 1, 4, 20)
        setedge!(g, 2, 4, 5)
        setedge!(g, 3, 2, 2)
        setedge!(g, 4, 5, 11)
        setedge!(g, 3, 5, 15)
    end
    @test dijkstra(g, 1) == [0, 5, 3, 10, 18]
end

@testset "jarnik (dense)" begin
    begin
        g = AdjacencyMatrix(6)
            setedge!(g, 1, 3, 7)
            setedge!(g, 3, 1, 7)
            setedge!(g, 1, 5, 9)
            setedge!(g, 5, 1, 9)
            setedge!(g, 3, 2, 5)
            setedge!(g, 2, 3, 5)
            setedge!(g, 3, 4, 1)
            setedge!(g, 4, 3, 1)
            setedge!(g, 3, 6, 2)
            setedge!(g, 6, 3, 2)
            setedge!(g, 2, 6, 6)
            setedge!(g, 6, 2, 6)
            setedge!(g, 4, 6, 2)
            setedge!(g, 6, 4, 2)
            setedge!(g, 5, 6, 1)
            setedge!(g, 6, 5, 1)
        result = AdjacencyMatrix(6)
            setedge!(result, 1, 3, 7)
            setedge!(result, 3, 1, 7)
            setedge!(result, 3, 2, 5)
            setedge!(result, 2, 3, 5)
            setedge!(result, 3, 4, 1)
            setedge!(result, 4, 3, 1)
            setedge!(result, 3, 6, 2)
            setedge!(result, 6, 3, 2)
            setedge!(result, 5, 6, 1)
            setedge!(result, 6, 5, 1)
    end
    @test jarnik(g, 1) == result
end


text:
#=
polylagrange:
- Julia version:
- Author: ymocquar
- Date: 2019-11-25
=#
include("polyexp.jl")
function getpolylagrange(k::Int64, j::Int64, N::DataType)
    @assert k <= j "_getpolylagrange(k=$k,j=$j) k must be less or equal to j"
    @assert N <: Signed "the type $N must be an Integer"
    result = Polynomial([one(Complex{Rational{N}})])
    for l = 0:j
        if l != k
            result *= Polynomial([l // 1, 1 // 1]) / (l - k)
        end
    end
    return result
end
function interpolate(tab, order, value, N::DataType)
    T = (N == BigInt) ? BigFloat : Float64
    res = zeros(Complex{T}, size(tab[1]))
    for i = 0:order
        res .+= getpolylagrange(i, order, N)(-value) * tab[i+1]
    end
    return res
end
function getpolylagrange(t::Vector{T}, k, j) where {T<:Number}
    result = Polynomial([one(T)])
    for l = 0:j
        if l != k
            result *= Polynomial([-t[l+1], one(T)]) / (t[k+1] - t[l+1])
        end
    end
    return result
end
function interpolate(tab_time::Vector{T}, tab, order, value) where {T<:Number}
    res = zeros(Complex{T}, size(tab[1]))
    for i = 0:order
        res .+= getpolylagrange(tab_time, i, order)(value) * tab[i+1]
    end
    return res
end
interpolate(tab, order, value) = interpolate(tab, order, value, BigInt)
struct CoefExpABRational
    tab_coef::Any
    function CoefExpABRational(
        order::Int64,
        epsilon::AbstractFloat,
        list_tau,
        dt::AbstractFloat,
    )
        n_tau = size(list_tau, 1)
        T = typeof(epsilon)
        tab_coef = zeros(Complex{T}, n_tau, order + 1, order + 1)
        N = T == BigFloat ? BigInt : Int64
        epsilon = rationalize(N, epsilon, tol = epsilon * 10 * Base.eps(T))
        dt = rationalize(N, dt, tol = dt * 10 * Base.eps(T))
        list_tau = rationalize.(N, list_tau)
        pol_x = Polynomial([0 // 1, 1 // dt])
        for j = 0:order
            for k = 0:j
                res = view(tab_coef, :, k + 1, j + 1)
                pol = getpolylagrange(k, j, N)
                pol2 = pol(pol_x)
                for ind = 1:n_tau
                    ell = list_tau[ind]
                    pol3 = undef
                    pol_int = if ell == 0
                        # in this case the exponentiel value is always 1
                        Polynomials.integrate(pol2)
                    else
                        pol3 = PolyExp(pol2, im * ell / epsilon, -im * ell * dt / epsilon)
                        Polynomials.integrate(pol3)
                    end
                    res[ind] = pol_int(dt) - pol_int(0)
                end
            end
        end
        return new(tab_coef)
    end
end
struct CoefExpAB
    tab_coef::Any
    tab_coef_neg::Any
    function CoefExpAB(order::Int64, epsilon::AbstractFloat, n_tau, dt)
        T = typeof(epsilon)
        N = T == BigFloat ? BigInt : Int64
        new_prec = precision(BigFloat)
        new_prec += T == BigFloat ? order * 16 : 0
        setprecision(BigFloat, new_prec) do
            list_tau = [collect(0:n_tau/2-1); collect(-n_tau/2:-1)]
            T2 = BigFloat
            epsilon = T2(epsilon)
            dt = T2(dt)
            tab_coef = zeros(Complex{T2}, n_tau, order + 1, order + 1)
            pol_x = Polynomial([0, 1 / dt])
            for j = 0:order
                for k = 0:j
                    res = view(tab_coef, :, k + 1, j + 1)
                    pol = getpolylagrange(k, j, N)
                    pol2 = pol(pol_x)
                    for ind = 1:n_tau
                        ell = list_tau[ind]
                        pol_int = if ell == 0
                            # in this case the exponentiel value is always 1
                            Polynomials.integrate(pol2)
                        else
                            Polynomials.integrate(
                                PolyExp(pol2, im * ell / epsilon, -im * ell * dt / epsilon),
                            )
                        end
                        res[ind] = pol_int(dt) - pol_int(0)
                    end
                end
            end # end of for j=....
        end # end of setprecision(...)
        # conversion to the new precision
        tab_coef = T.(real(tab_coef)) + im * T.(imag(tab_coef))
        tab_coef_neg = -conj(tab_coef)
        return new(tab_coef, tab_coef_neg)
    end
end


text:
function mean(a)
    sum = 0
    count = length(a)
    for x in a
        sum = sum + x
    end

    m = sum / count
    #println(sum)
    println()

    m
end

b = Array(Int32, 100000000)
rand!(1:20, b)

@time mean(b)

max_stars_repo_name: dwdii/DataAcqMgmt

max_issues_repo_name: dwdii/DataAcqMgmt

max_forks_repo_name: dwdii/DataAcqMgmt


text:
using Diversity
using EcoBase
using DataFrames

"""
    jostalpha(proportions::AbstractMatrix, qs)

Calculates Jost's alpha diversity of a series of columns representing
independent community counts, for a series of orders, repesented as a
vector of qs. This is just the naive-community ecosystem diversity
divided by the naive-community beta diversity.

# Arguments:

- `proportions` relative proportions of different individuals / species
                in population (vector, or matrix where columns are
                for individual sub-communities)

- `qs` single number or vector of orders of diversity measurement

# Returns:

- DataFrame of diversities
"""
function jostalpha(proportions::AbstractMatrix, qs)
    md = metacommunityDiversity(RawAlpha(Metacommunity(proportions)), qs)
    md[!,:diversity] ./= qD(reshape(mapslices(sum, proportions, dims=(1,)),
                            size(proportions)[2]), qs)
    md[!,:measure] .= "JostAlpha"
    return md
end

function jostalpha(asm::EcoBase.AbstractAssemblage, qs)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return jostalpha(occurrences(asm), qs)
end

"""
    jostbeta(proportions::AbstractMatrix, qs)

Calculates Jost's beta diversity of a series of columns representing
independent community counts, for a series of orders, repesented as a
vector of qs. This is just the naive gamma diversity divided by
Jost's alpha diversity

# Arguments:

- `proportions` relative proportions of different individuals / species
                in population (vector, or matrix where columns are
                for individual sub-communities)

- `qs` single number or vector of orders of diversity measurement

# Returns:

- DataFrame of diversities
"""
function jostbeta(proportions::AbstractMatrix, qs)
    md = metacommunityDiversity(Gamma(Metacommunity(proportions)), qs)
    ja = jostalpha(proportions, qs)
    j = innerjoin(md, ja, on=[:q, :type_level, :type_name,
                         :partition_level, :partition_name, :div_type],
             makeunique=true)
    j[!,:diversity] ./= j[!,:diversity_1]
    j[!,:measure] .= "JostBeta"
    select!(j, Not([:diversity_1, :measure_1]))
    return j
end

function jostbeta(asm::EcoBase.AbstractAssemblage, qs)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return jostbeta(occurrences(asm), qs)
end


text:
"""
The PreAmp supertype corresponds to the charge sensitive preamplifier component
"""
abstract type PreAmp end

"""
The GenericPreAmp is a dummy PreAmp model that accounts for
    effects such as decay and rise time, offset, noise and gain.

This dummy model involves no real electronics response, only gain.

This struct is currently mutable because in the case of noise modelling based on data,
    the gain has to match the one in the baselines extracted from data,
    so in that case the PreAmp is initialized with gain = 0, and the gain
    is calculated later. I assume this is temporary, as the user should know which
    gain / max_e the data was produced with and give it as simulation input.
"""
@with_kw mutable struct GenericPreAmp <: PreAmp
    "PreAmp exp decay time"
    τ_decay::typeof(1.0*μs_unit) = T(5)*u"μs"

    "PreAmp rise time"
    τ_rise::typeof(1.0*ns_unit) = T(15)*u"ns"

    "Preamp offset in keV"
    offset::typeof(1.0*energy_unit) = 0u"keV";

    "Maximum DAQ energy"
    max_e::typeof(1.0*energy_unit) = 10_000u"keV" # == typemax(UInt16)

    "Preamp gain"
    # If offset = 0, it means gain has to be inferred (from the data baseline offset)
    # How does it end up being Float64?
    gain = offset == 0u"keV" ? 0 : typemax(UInt16) / ((max_e+offset)/uconvert(u"keV", germanium_ionization_energy))

    "Preamp Gaussian noise"
    noise_σ::typeof(1.0*energy_unit)= 0u"keV"
end

"""
CC2 circuit model (ToDo)
"""
struct CC2 <: PreAmp end

"""
    GenericPreAmp(sim_conf)

PropDict -> GenericPreAmp

Construct a GenericPreAmp instance based on given simulation configuration <sim_conf>
"""
function GenericPreAmp(sim_conf::PropDict)
    GenericPreAmp(
        τ_decay=T(sim_conf.setup.preamp.t_decay)*u"μs",
        τ_rise=T(sim_conf.setup.preamp.t_rise)*u"ns",
        offset = T(sim_conf.setup.preamp.offset)u"keV",
        max_e = T(sim_conf.setup.preamp.max_e)u"keV",
        noise_σ = haskey(sim_conf.setup.preamp, :noise_sigma) ? T(sim_conf.setup.preamp.noise_sigma)u"keV" : 0u"keV"
    )
end

"""
    PreAmp(sim_conf)

PropDict -> <PreAmp>

Construct a PreAmp supertype instance based on given simulation configuration <sim_conf>.
Returned type depends on the settings in <sim_conf>.
Currently only GenericPreAmp is available.
"""
function PreAmp(sim_conf::PropDict)
    GenericPreAmp(sim_conf)
end

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

"""
    simulate(wf, preamp)

RDWaveform, GenericPreAmp -> RDWaveform

Simulate the effecs of the preamp on the waveform.
"""
function simulate(wf::RDWaveform, preamp::GenericPreAmp)
    ## rise and decay time
    csa_filter = RadiationDetectorDSP.simple_csa_response_filter(
        preamp.τ_rise / step(wf.time),
        uconvert(u"ns", preamp.τ_decay) / step(wf.time))

    wf_preamp = RDWaveform(wf.time, filt(csa_filter, wf.value))

    ## noise
    wf_preamp = simulate_noise(wf_preamp, preamp)

    ## offset
    # wf values are in eV but without u"eV" attached
    offset = ustrip(uconvert(u"eV", preamp.offset))
    wf_preamp = RDWaveform(wf_preamp.time, wf_preamp.value .+ offset)

    ## gain
    wf_preamp = RDWaveform(wf_preamp.time, wf_preamp.value * preamp.gain)

    ## what is this step?
    # normalizing by germanium ionization energy?
    # basically converting to number of charge carriers?
    # but after the gain?
    # I guess it's like a dummy MeV -> ADC (calibration curve?)
    wf_preamp = RDWaveform(wf_preamp.time, wf_preamp.value ./ ustrip(germanium_ionization_energy))

    wf_preamp
end


text:
@safetestset "Assets functionality" begin

  @safetestset "Assets paths" begin
    using Genie, Genie.Assets

    Genie.config.base_path = "/proxy/8000"

    @test Genie.Assets.asset_path("s.js") == "/proxy/8000/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", host = "") == "/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", host = "/") == "/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", host = "//") == "/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", path = "/") == "/proxy/8000/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", path = "//") == "/proxy/8000/assets/js/s.js"
    @test Genie.Assets.asset_path("s.js", path = "foo") == "/proxy/8000/assets/js/foo/s.js"
    @test Genie.Assets.asset_path("s.js", path = "foo/bar") == "/proxy/8000/assets/js/foo/bar/s.js"
    @test Genie.Assets.asset_path("s.js", path = "foo/bar/baz") == "/proxy/8000/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", path = "/foo/bar/baz") == "/proxy/8000/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", path = "/foo/bar/baz/") == "/proxy/8000/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", host = "abc", path = "/foo/bar/baz/") == "/abc/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", host = "abc/def", path = "/foo/bar/baz/") == "/abc/def/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", host = "/abc/def", path = "/foo/bar/baz/") == "/abc/def/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js", host = "/abc/def/", path = "/foo/bar/baz/") == "/abc/def/assets/js/foo/bar/baz/s.js"

    Genie.config.base_path = "/proxy/8000/"
    @test Genie.Assets.asset_path("s.js", path = "/foo/bar/baz/") == "/proxy/8000/assets/js/foo/bar/baz/s.js"

    Genie.config.base_path = "proxy/8000"
    @test Genie.Assets.asset_path("s.js", path = "/foo/bar/baz/") == "/proxy/8000/assets/js/foo/bar/baz/s.js"
    @test Genie.Assets.asset_path("s.js") == "/proxy/8000/assets/js/s.js"
  end;

end;


text:
# Autogenerated wrapper script for FFTW_jll for x86_64-apple-darwin-libgfortran5
export libfftw3, libfftw3_mpi, libfftw3f, libfftw3f_mpi

using MPItrampoline_jll
using MicrosoftMPI_jll
JLLWrappers.@generate_wrapper_header("FFTW")
JLLWrappers.@declare_library_product(libfftw3, "@rpath/libfftw3.3.dylib")
JLLWrappers.@declare_library_product(libfftw3_mpi, "@rpath/libfftw3_mpi.3.dylib")
JLLWrappers.@declare_library_product(libfftw3f, "@rpath/libfftw3f.3.dylib")
JLLWrappers.@declare_library_product(libfftw3f_mpi, "@rpath/libfftw3f_mpi.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(MPItrampoline_jll, MicrosoftMPI_jll)
    JLLWrappers.@init_library_product(
        libfftw3,
        "lib/libfftw3.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3_mpi,
        "lib/libfftw3_mpi.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3f,
        "lib/libfftw3f.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfftw3f_mpi,
        "lib/libfftw3f_mpi.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ArgParse,KUparser
include("conll17.jl")
macro tm(_x) :(if o[:fast]; $(esc(_x)); else; info("$(now()) "*$(string(_x))); $(esc(_x)); end) end
macro msg(x); :(if !o[:fast]; info($x); end); end

# Some default data
grctxt = "/mnt/ai/data/nlp/conll17/ud-treebanks-conll2017/UD_Ancient_Greek/grc-ud-dev.txt"
grcudp = "/mnt/ai/data/nlp/conll17/UDPipe/udpipe-ud-2.0-conll17-170315/models/ancient_greek-ud-2.0-conll17-170315.udpipe"

function main(args="")
    isa(args, AbstractString) && (args=split(args))
    s = ArgParseSettings()
    s.description="Random parser"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        # Parsing options
        ("--parse"; default=grctxt; help="file in plain text format to be parsed")
        ("--udpipe"; default=grcudp; help="UDpipe model file")
        ("--ptype"; default="ArcEagerR1"; help="Parser type")
        ("--fast"; action=:store_true; help="fewer messages")
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
    end

    # Process args
    o = parse_args(args, s; as_symbols=true)
    @msg(s.description)
    @msg(string("opts=",[(k,v) for (k,v) in o]...))
    if o[:seed] > 0; setseed(o[:seed]); end
    ptype = eval(parse(o[:ptype]))

    # Parse
    if o[:parse] != nothing
        # we need to test with predicted toks and tags from txt data.
        txtfile = tempname()
        @tm run(pipeline(`udpipe --tokenize --tag $(o[:udpipe]) $(o[:parse])`, txtfile))
        @tm global corpus = readconllu(txtfile)
        @tm parses = rparse(ptype, corpus)
        @tm writeconllu(corpus, parses, txtfile)
        rm(txtfile)
    end
end

if basename(PROGRAM_FILE)==basename(@__FILE__); main(ARGS); end


text:
#!/usr/bin/env julia
# separate file because I can't get an equivalent of @__FILE__ to work in IJulia
# TODO file an issue about that and replace this line with a link to it
function clone_or_pull(url; pkg=false)
    dir = basename(url)
    if pkg
        name = splitext(dir)[1]
        if isdir(Pkg.dir(name))
            run(setenv(`git pull`; dir=Pkg.dir(name)))
        else
            Pkg.clone(string("https://", url))
            symlink(Pkg.dir(name), dir)
            # for my own pushing convenience because ssh is annoying
            run(setenv(`git remote add tkelman https://tkelman@$url`; dir=dir))
        end
    else
        if isdir(dir)
            run(setenv(`git pull`; dir=dir))
        else
            run(`git clone https://$url`)
            # for my own pushing convenience because ssh is annoying
            run(setenv(`git remote add tkelman https://tkelman@$url`; dir=dir))
        end
    end
end
cd(dirname(dirname(@__FILE__))) do
    clone_or_pull("github.com/tkelman/Rmath-julia")
    clone_or_pull("github.com/tkelman/Rmath.jl"; pkg=true)
end