text:
abstract type AbstractThermoSolver end
abstract type AbstractHelmholtzSolver <: AbstractThermoSolver end
#volume solver
#with P and T, finds a volume that satisfies pressure(v,T) = P
#without falling on the false loop.
#overload this for cubics and consistent EOS
function volume_solver(model::AbstractHelmholtzModel,P,T,x,v0 = nothing;no_pts = 7,all_volumes = false)
    p(z) = core_pressure(model,z,T,x)
    fp(z) = p(z)-P
    dfp(z) = ForwardDiff.derivative(fp,z)
    if isnothing(v0)
    min_v = dot(covolumes(model),x)
    max_v = 40*T/P #approx 5 times ideal gas
    #this is to be sure that (min_v,max_v) is a bracketing interval
    while fp(max_v) > 0
        max_v *= 2
    end
    while fp(min_v) < 0
        min_v *=0.5
    end
    vv = find_zeros(fp,min_v,max_v,no_pts=no_pts)
    if all_volumes == true
        return vv
    end
    if length(vv) == 1
        return vv[1]
    else length(vv) >= 1
        v1 = vv[1]
        v2 = vv[end]
        g1 = core_gibbs(model,P,v1,T,x)
        g2 = core_gibbs(model,P,v2,T,x)
        if g1 < g2
            return v2
        else
            return v1
        end
    end
else
    try
        return find_zero(fp,v0)
    catch
        return volume_solver(model,P,T,x)
    end
end
end

function core_mol_volume(model::AbstractHelmholtzModel,P0,T0,x0)
    v = volume_solver(model,P0,T0,x0)
    return v
end

function core_mol_volume(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    v = volume_solver(model,P0,T0,x0,v0)
    return v
end
function core_mol_density(model::AbstractHelmholtzModel,P0,T0,x0)
    v = volume_solver(model,P0,T0,x0)
    return inv(v)
end

function core_mol_density(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    return inv(v)
end

function core_mass_volume(model::AbstractHelmholtzModel,P0,T0,x0)
    v = volume_solver(model,P0,T0,x0)
    k = dot(molecular_weight(model),x0)*1e-03 #kg/mol
    return v/k
end

function core_mass_volume(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    v = volume_solver(model,P0,T0,x0,v0)
    k = dot(molecular_weight(model),x0)*1e-03 #kg/mol
    return v/k
end

function core_mass_density(model::AbstractHelmholtzModel,P0,T0,x0)
    v = volume_solver(model,P0,T0,x0)
    k = dot(molecular_weight(model),x0)*1e-03 #kg/mol
    return k/v
end

function core_mass_density(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    v = volume_solver(model,P0,T0,x0,v0)
    k = dot(molecular_weight(model),x0)*1e-03
    return k/v
end

function mol_volume(model::AbstractHelmholtzModel,P0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    return core_mol_volume(model,P1,T1,x0)*1.0u"m^3/mol"
end

function mol_volume(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    v1 = _transform_v(v0,molecular_weight(model),x0)
    return core_mol_volume(model,P1,v1,T1,x0)*1.0u"m^3/mol"
end

function mol_density(model::AbstractHelmholtzModel,P0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    return core_mol_density(model,P1,T1,x0)*1.0u"mol/m^3"
end

function mol_density(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    v1 = _transform_v(v0,molecular_weight(model),x0)
    return core_mol_density(model,P1,v1,T1,x0)*1.0u"mol/m^3"
end

function mass_volume(model::AbstractHelmholtzModel,P0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    return core_mass_volume(model,P1,T1,x0)*1.0u"m^3/kg"
end

function mass_volume(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    v1 = _transform_v(v0,molecular_weight(model),x0)
    return core_mass_volume(model,P1,v1,T1,x0)*1.0u"m^3/kg"
end

function mass_density(model::AbstractHelmholtzModel,P0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    return core_mass_density(model,P1,T1,x0)*1.0u"kg/m^3"
end

function mass_density(model::AbstractHelmholtzModel,P0,v0,T0,x0)
    P1 = _transform_P(P0)
    T1 = _transform_T(T0)
    v1 = _transform_v(v0,molecular_weight(model),x0)
    return core_mass_density(model,P1,v1,T1,x0)*1.0u"kg/m^3"
end

function pt_flash(method::T,
    model::T2,
    P0,
    T0,
    x0,
    options=nothing) where T <: AbstractHelmholtzSolver where T2 <: AbstractHelmholtzModel
    throw(error("this method is not implemented."))
end

function pt_flash(
    model::T2,
    P0::Unitful.Pressure,
    T0::Unitful.Temperature,
    x0,
    method::T,
    options=nothing) where T <: AbstractHelmholtzSolver where T2 <: AbstractHelmholtzModel
    return core_pt_flash(method,
    model,
    _transform_P(P0),
    _transform_T(T0),
    x0,
    options)
end

function pt_flash(model::AbstractHelmholtzModel,
    phase0::HelmholtzPhase,
    method::AbstractHelmholtzSolver,options)
    P0 = core_pressure(model,phase0)
    T0 = temperature(phase0)
    x0 = mol_fraction(model,phase0)
    return core_pt_flash(model,P0,T0,x0,method,options)
end

#spinodal solver
#solves the equation dp/dv = 0 and gives the outer values as a tuple
function spinodal_solver(model::AbstractHelmholtzModel,T,x,predictor::Symbol=:PengRobinson)
    fp(z) = core_pressure(model,z,T,x)
    dfp(z) = core_dpdv(model,z,T,x)
    d2p(z) = ForwardDiff.derivative(dfp,z)
    min_v = dot(covolumes(model),x)
    P0 = dot(pressure_predictor(model,T,predictor),x) #by raoult, sum(xiP0i) = P
    max_v = 40*T/P0 # 0.01 atm
    while dfp(max_v) > 0
        max_v *= 2
    end
    while dfp(min_v) > 0
        min_v *=0.5
    end
    #this is to be sure that (min_v,max_v) is a bracketing interval
    dv = find_zeros(dfp,min_v,max_v,no_pts=21)
    if length(dv) > 2
        return (dv[1],dv[end],:success)
    elseif length(dv) == 0
        return (NaN,NaN,:zero_points)
    elseif length(dv) == 1
        return (dv[1],dv[1],:one_point)
    elseif length(dv) == 2
        if (d2p(dv[1]) > 0.0) && (d2p(dv[2]) < 0.0) #v1 checked
            return (dv[1],dv[2],:success)
        end
        return (dv[1],dv[2],:two_points)
    end
end

#pressure_predictor:
#tries to predict a saturation pressure given only the saturation temperature
#TODO: include a trait to let the model choose his own predictor?
#let the user build his own predictors?
#an ideal predictor would be an antoine equation, for example
pressure_predictor(model,T,pred::Symbol) = pressure_predictor(model,T,Val(pred))
function pressure_predictor(model,T,pred::Val{:PengRobinson})
    Tc = critical_temperature(model)
    Tr = T ./ Tc
    m(ω) = 0.37464 + 1.54226*ω - 0.26992*ω^2
    _α(ω,_Tr) = (1+m(ω)*(1-sqrt(_Tr)))^2
    α =_α.(acentric_factor(model),Tr)
    Pc = critical_pressure(model)
    aa =-2.605272488488440
    bb = -9.017571450539830
    cc = -19.896014683288000
    dd = 0.579677284391001
    ee = 22.501011849124900
    ff = -0.041440165126830
    xx = α ./ Tr
    @. xx = (aa + cc*xx + ee*xx^2)/(1 + bb*xx + dd*xx^2 + ff*xx^3)
    xx .= exp.(xx) .* Tr
    return xx .* Pc
end

function core_pure_t_flash(model,T,x,options=nothing)
    if (Tc = dot(x,critical_temperature(model))) ≈ T
        return helmholtz_phase(model, dot(critical_volume(model),x), T, x)
    elseif T > Tc
        throw(error("the phase is supercritical at temperature = $T K"))
    end
    p_pred = dot(x,pressure_predictor(model,T,:PengRobinson))
    P = p_pred
    p(z) = core_pressure(model,z,T,x)
    fp(z) = p(z)-P
    min_v = dot(covolumes(model),x)
    max_v = 40*T/p_pred #approx 5 times ideal gas
    #this is to be sure that (min_v,max_v) is a bracketing interval
    while fp(max_v) > 0
        max_v *= 2
    end
    while fp(min_v) < 0
        min_v *=0.5
    end
    vv = find_zeros(fp,min_v,max_v,no_pts=21)
    #@show P
    if length(vv) <= 1
        throw(error("the phase is stable at temperature = $T K"))
    elseif length(vv) == 2
        throw(error("proper initial points not found at temperature = $T K"))
    elseif length(vv) > 2
       v1 = vv[1]
       v2 = vv[end]
        p1 = p(v1)
        p2 = p(v2)
        px = p_pred
        A(z) = core_helmholtz(model,z,T,x)
        v1old = 0
        v2old = 0
        for i = 1:20
            if i > 1
                v1old = v1
                v2old = v2
            end
            _v1 = Threads.@spawn volume_solver(model,px,T,x,0.9*v1)
            _v2 = Threads.@spawn volume_solver(model,px,T,x,1.1*v2)

            v1 = fetch(_v1)
            v2 = fetch(_v2)
            if abs(v1-v1old)/v1 < 1e-15 && i > 1
                #println("v1 condition")
                break
            elseif abs(v2-v2old)/v2 < 1e-15 && i > 1
                #println("v2 condition")
                break
            end
            pold = px
            px = (A(v1)-A(v2))/(v2-v1)
            #@show px
            if abs(pold-px) < 1e-10*px
                #println("P condition")
                break
            end

        end
        phase1 = helmholtz_phase(model,v1,T,x)
        phase2 = helmholtz_phase(model,v2,T,x)
        return [phase1,phase2]
    end
end

function core_pure_p_flash(model,P,x,options=nothing)
    if (Pc = dot(x,critical_pressure(model))) ≈ P
        return helmholtz_phase(model, dot(critical_volume(model),x), T, x)
    elseif P > Pc
        throw(error("the phase is supercritical at temperature = $T K"))
    end
    Tc = dot(x,critical_temperature(model))
    _t_pred0(T) =dot(x,pressure_predictor(model,T,:PengRobinson)) - P
    t_pred = find_zero(_t_pred0,Tc)
    vv = volume_solver(model,P,t_pred,x,all_volumes=true,no_pts=21)
    if length(vv) <= 1
        throw(error("the phase is stable at temperature = $T K"))
    elseif length(vv) == 2
        throw(error("proper initial points not found at temperature = $T K"))
    elseif length(vv) > 2
       v1 = vv[1]
       v2 = vv[end]
        px = P
        _A(z,T) = core_helmholtz(model,z,T,x)
        v1old = 0.0
        v2old = 0.0
        Told = -100.0
        Tx = t_pred
        for i = 1:20
            if i > 1
                v1old = v1
                v2old = v2
            end
            A = z-> _A(z,Tx)
            _v1 = Threads.@spawn volume_solver(model,P,Tx,x,0.9*v1)
            _v2 = Threads.@spawn volume_solver(model,P,Tx,x,1.1*v2)
            v1 = fetch(_v1)
            v2 = fetch(_v2)
            if abs(v1-v1old)/v1 < 1e-15 && i > 1
                #println("v1 condition")
                break
            elseif abs(v2-v2old)/v2 < 1e-15 && i > 1
                #println("v2 condition")
                break
            end
            _px(T) = (_A(v1,T)-_A(v2,T))/(v2-v1) - P
            Told = Tx
            Tx = find_zero(_px,Tx)
            if abs(Told-Tx) < 1e-10*Tx
                #println("P condition")
                break
            end
        end

        phase1 = helmholtz_phase(model,v1,Tx,x)
        phase2 = helmholtz_phase(model,v2,Tx,x)
        return [phase1,phase2]
    end
end

#function predict_saturation_pressure(model::AbstractHelmholtzModel,T0,x0,method)
function molar_volume_rackett(model::AbstractHelmholtzModel,T,x0)
    Pc = critical_pressure(model)
    Tc = critical_temperature(model)
    vc = critical_volume(model)
    Zc = core_compressibility_factor.(model,Pc,vc,Tc,Ref(1.0))
   # return vc
    return R_GAS .* Tc ./ Pc .* Zc .^ (1.0 .+ (1.0 .- T./ Tc).^(2/7) )
end

function pure_t_flash(
    model::T2,
    T0::Unitful.Temperature,
    x0,
    options=nothing) where T <: AbstractHelmholtzSolver where T2 <: AbstractHelmholtzModel
    return core_pure_t_flash_2(
    model,
    _transform_T(T0),
    x0,
    options)
end

function pure_t_flash(
    model::T2,
    T0::Unitful.Temperature,
    x0 = [1.0],
    options=nothing) where T <: AbstractHelmholtzSolver where T2 <: AbstractHelmholtzModel
    return core_pure_t_flash_2(
    model,
    _transform_T(T0),
    x0,
    options)
end


text:
module ClineSimulation

export binomial_migration,
       migration,
       correction,
       initialfreqs

using Distributions: Binomial,  pdf

include("migration.jl")
include("initialfrequencies.jl")

end # module

ext: jl


text:
__precompile__()

module EgyptianFractions

using Base.Iterators: repeated

export efgreedy, efoddgreedy, efharmonic, engelexpand, efengel

function _prep(r::Rational)
  @assert(denominator(r) != 0, "denominator must be > 0")
  rem = abs(big(r))
  ef = BigInt[]
  if rem ≥ 1
    f = floor(rem)
    rem -= f
    append!(ef, collect(repeated(1, Int(f))))
  end
  return (rem, ef)
end

function _greedyloop!(ef::Vector{BigInt}, r::Rational{BigInt})
  c = ceil(1//r)
  push!(ef, c)
  return r - 1//c
end

"""
Performs a greedy (Fibonacci–Sylvester) expansion of the given rational number into a sum of fractions of the form `1/a_1 + 1/a_2 + ...`

This function returns only the denominators of the expansion (i.e., only `[a_1, a_2, ...]`).  If the given rational number `r` satisfies `r ≥ 1`, the first `n` elements of the expansion will be 1, where `n = floor(r)`.  If the given rational number `r` satisfies `r < 0`, all returned denominators will be less than 0, so that the relationship

    r == sum(1 .// efgreedy(r))

is always true.
"""
function efgreedy(r::Rational; nmax::Int = typemax(Int))
  (rem, ef) = _prep(r)
  i = nmax
  while rem != 0 && i > 0
    rem = _greedyloop!(ef, rem)
    i -= 1
  end
  if r < 0
    ef .*= -1
  end
  return ef
end

efgreedy(r::Real; nmax::Int = typemax(Int)) = efgreedy(Rational(big(r)), nmax=nmax)

function _oddgreedyloop!(ef::Vector{BigInt}, r::Rational{BigInt})
  c = ceil(1//r)
  if iseven(numerator(c))
    c += 1
  end
  push!(ef, c)
  return r - 1//c
end

"""
Performs a greedy (Fibonacci–Sylvester) expansion of the given rational number into a sum of fractions of the form `1/a_1 + 1/a_2 + ...`, but only using odd denominators.

It can be shown that, for any rational number `x/y` where `y` is odd, you can write this as a finite sum of fractions where all denominators are odd.  This method sometimes produces an expansion with fewer elements than the traditional greedy algorithm.

This function will throw an exception if the denominator of the rational number supplied is even.

This function returns only the denominators of the expansion (i.e., only `[a_1, a_2, ...]`).  If the given rational number `r` satisfies `r ≥ 1`, the first `n` elements of the expansion will be 1, where `n = floor(r)`.  If the given rational number `r` satisfies `r < 0`, all returned denominators will be less than 0, so that the relationship

    r == sum(1 .// efoddgreedy(r))

is always true.
"""
function efoddgreedy(r::Rational; nmax::Int = typemax(Int))
  @assert isodd(denominator(r)) "denominator of rational ($(denominator(r))) must be odd"
  (rem, ef) = _prep(r)
  i = nmax
  while rem != 0 && i > 0
    rem = _oddgreedyloop!(ef, rem)
    i -= 1
  end
  if r < 0
    ef .*= -1
  end
  return ef
end

efoddgreedy(r::Real; nmax::Int = typemax(Int)) = efoddgreedy(Rational(big(r)), nmax=nmax)

"""
Performs a harmonic expansion of the given rational number into a sum of fractions of the form `1/2 + 1/3 + ...`, and concludes the remainder of the expansion using `efgreedy` if the given rational number cannot be represented as a sum in the harmonic sequence.

If the second argument is specified, the expansion will begin using that value as the first denominator.  In other words, while `efharmonic(r)` will return an array beginning `[2, 3, ...]`, `efharmonic(r, 3)` will return an array beginning `[3, 4, ...]`, and so on.

This function will throw an exception if the second argument `f ≤ 1`.

This function returns only the denominators of the expansion (i.e., only `[a_1, a_2, ...]`).  If the given rational number `r` satisfies `r ≥ 1`, the first `n` elements of the expansion will be 1, where `n = floor(r)`.  If the given rational number `r` satisfies `r < 0`, all returned denominators will be less than 0, so that the relationship

    r == sum(1 .// efharmonic(r))

is always true.
"""
function efharmonic(r::Rational, first::Int = 2; nmax::Int = typemax(Int))
  @assert(first ≥ 2, "harmonic series must start at 2 or greater ($first given)")
  (rem, ef) = _prep(r)
  s = Rational{BigInt}(0)
  i = big(first)
  j = nmax
  while s ≤ rem && j > 0
    hh = 1//i
    if s + hh > rem
      break
    end
    push!(ef, i)
    s += hh
    i += 1
    j -= 1
  end
  rem -= s
  while rem != 0 && j > 0
    rem = _greedyloop!(ef, rem)
    j -= 1
  end
  if r < 0
    ef .*= -1
  end
  return ef
end

efharmonic(r::Real, first::Int = 2; nmax::Int = typemax(Int)) = efharmonic(Rational(big(r)), first, nmax=nmax)

function _engelloop!(ef::Vector{BigInt}, r::Rational{BigInt})
  c = ceil(1//r)
  push!(ef, c)
  return r * c - 1
end

"""
Performs an Engel expansion of the given rational number into a sum of fractions of the form `1/a + 1/(a*b) + 1/(a*b*c)...`.

This function returns only the unique denominators of the expansion (i.e., only `[a, b, ...]`).  If the given rational number `r` satisfies `r ≥ 1`, the first `n` elements of the expansion will be 1, where `n = floor(r)`.
"""
function engelexpand(r::Rational; nmax::Int = typemax(Int))
  (rem, ef) = _prep(r)
  i = nmax
  while rem != 0 && i > 0
    rem = _engelloop!(ef, rem)
    i -= 1
  end
  # Convention:  the Engle expansion of a negative number will have ef[1] < 0.
  # This way, all elements of cumprod(ef) will be negative.
  if !isempty(ef) && r < 0
    ef[1] *= -1
  end
  return ef
end

engelexpand(r::Real; nmax::Int = typemax(Int)) = engelexpand(Rational(big(r)), nmax=nmax)

"""
Performs an Engel expansion of the given rational number into a sum of fractions of the form `1/a + 1/(a*b) + 1/(a*b*c)...`.

This function returns only the denominators of the expansion (i.e., only `[a, a*b, ...]`).  If the given rational number `r` satisfies `r ≥ 1`, the first `n` elements of the expansion will be 1, where `n = floor(r)`.  If the given rational number `r` satisfies `r < 0`, all returned denominators will be less than 0, so that the relationship

    r == sum(1 .// efengel(r))

is always true.
"""
function efengel(r::Rational; nmax::Int = typemax(Int))
  ef = engelexpand(r, nmax=nmax)
  cumprod!(ef, ef)
  return ef
end

efengel(r::Real; nmax::Int = typemax(Int)) = efengel(Rational(big(r)), nmax=nmax)

end # module


text:
using SnoopCompile

# using runtests:
@snoopi_bot "Juno"


text:
using LinearAlgebra
using Test

using UnitBlockAngular

@testset "UnitBlockAngular" begin

Aref = [
    1.0 0.0 1.1 1.2 1.3 1.4;
    0.0 1.0 2.1 2.2 2.3 2.4;
    1.0 0.0 3.1 3.2 3.3 3.4;
    0.0 0.0 1.0 0.0 1.0 0.0;
    0.0 0.0 0.0 1.0 0.0 1.0
]

m0, n0, n, R = 3, 2, 4, 2
B0 = Aref[1:3, 1:2]
B  = Aref[1:3, 3:end]
blocks = [1, 2, 1, 2]

A = UnitBlockAngularMatrix{Float64}(3, 2, 4, 2);
A.B0 .= B0
A.B  .= B
A.blockidx .= blocks

@test size(A) == (m0+R, n0+n)
@test A == Aref

A = UnitBlockAngularMatrix(B0, B, R, blocks)
@test size(A) == (m0+R, n0+n)
@test A == Aref

x, xref = ones(n0 + n), ones(n0 + n)
y, yref = zeros(m0 + R), zeros(m0 + R)

@testset "Matrix-vector" begin
    mul!(yref, Aref, xref)
    mul!(y, A, x)
    @test y ≈ yref

    mul!(yref, Aref, xref, 100.0, 10.0)
    mul!(y, A, x, 100.0, 10.0)
    @test y ≈ yref

    y .= 1
    yref .= 1
    mul!(xref, Aref', yref)
    mul!(x, A', y)
    @test x ≈ xref

    y .= 1
    yref .= 1
    mul!(xref, Aref', yref, 100.0, 10.0)
    mul!(x, A', y, 100.0, 10.0)
    @test y ≈ yref
end

end  # testset

# Factorization
include("kkt_solver.jl")


text:
function linear_QNet(input_size::T, hidden_size::T, output_size::T) where {T<:Integer}
    # Creating the model
    model = Chain(
		# First linear layer
		Dense(input_size => hidden_size, relu),

        # Second linear layer
		Dense(hidden_size => output_size)
    )

    return model
end


text:
using StructuredOptimization
using BenchmarkTools
using JuMP
using SCS
using ECOS

function set_up(S::Int) #S scales problem

	n = S
	m = div(S,5)
	SNR = 5

	srand(123)

	A = sprandn(m,n,5/n)
	println("n = $n, m = $m, nnz(A) = $(countnz(A))")
	x0 = zeros(n)
	x0[randperm(n)[1:div(n,4)+1]] = randn(div(n,4)+1)

	y = A*x0
	y += 10^(-SNR/10)*sqrt(var(y))*randn(length(y))
	lambda = 0.01*norm(A'*y,Inf)
	x = StructuredOptimization.Variable(n)
	@minimize ls(A*x-y)+lambda*norm(x,1) with ZeroFPR(verbose = 0, tol =1e-12)
	setup = A, y, lambda, n
	setup_JuMP = create_JuMP_model(A, y, lambda, n)
	return ~x, setup, setup_JuMP
end

function solve_problem(slv::S, A, y, lambda, n, M, xJ) where {S <: StructuredOptimization.ForwardBackwardSolver}
	x = StructuredOptimization.Variable(n)
	it, = @minimize ls(A*x-y)+lambda*norm(x,1) with slv
	return ~x, it, :UserLimit
end

function create_JuMP_model(A, y, lambda, n)
	M = Model()
	@variables M begin
		x[1:n]
		t[1:n]
		w
	end
	@objective(M,Min,[0.5;lambda*ones(n)]'*[w;t])
	@constraint(M, soc, norm( [1-w;2*(A*x-y)] ) <= 1+w)
	@constraint(M,  x .<= t)
	@constraint(M, -t .<= x)
	return M, x
end

function solve_problem(slv::S, A, y, lambda, n, M, xJ) where {S <: MathProgBase.SolverInterface.AbstractMathProgSolver}
	setsolver(M, slv)
	status = JuMP.solve(M)
	return getvalue(xJ), 0, status
end

function benchmark_LASSO(slv, eps_max, nvar)

	slv = eval(parse(slv))

	verbose, samples, seconds = 0, 5, 2*60
	xopt, setup, setup_JuMP = set_up(nvar)

	if slv == PG
		if nvars == 100000 #100000
			iterations = 1:div(100000,6):100000
		elseif nvars == 10000 #10000
			iterations = 1:div(50000,6):50000
		elseif nvars == 1000 #7875
			iterations = 1:div(7857,6):7857
		end
	elseif slv == FPG
		if nvars == 100000 #28941
			iterations = 1:div(28941,6):28941
		elseif nvars == 10000 #12534
			iterations = 1:div(13000,6):13000
		elseif nvars == 1000 #4027
			iterations = 1:div(4027,6):4027
		end
	elseif slv == ZeroFPR
		if nvars == 100000 #1019
			iterations = 1:div(1019,6):1019
		elseif nvars == 10000 #640
			iterations = 1:div(640,6):640
		elseif nvars == 1000 #163
			iterations = 1:div(163,6):163
		end
	elseif slv == PANOC
		if nvars == 100000 #1019
			iterations = 1:div(2*1019,6):2*1019
		elseif nvars == 10000 #640
			iterations = 1:div(2*640,6):2*640
		elseif nvars == 1000 #163
			iterations = 1:div(2*163,6):2*163
		end
	elseif slv == ECOSSolver
		if nvars == 100000 #104
			iterations = 1:div(104,6):104
		elseif nvars == 10000 #30
			iterations = 1:div(35,6):35
		elseif nvars == 1000 #23
			iterations = 1:div(25,6):25
		end
	elseif slv == SCSSolver
		if nvars == 100000
			itSCS = 2680 #with default (low) tolerance
			iterations = 1:div(itSCS,3):6*div(itSCS,3)
		elseif nvars == 10000
			itSCS = 420  #with default (low) tolerance
			iterations = 1:div(itSCS,3):8*div(itSCS,3)
		elseif nvars == 1000
			itSCS = 180  #with default (low) tolerance
			iterations = 1:div(itSCS,3):6*div(itSCS,3)
		end
	end

	#iterations = [100000]
	err = zeros(length(iterations))
	t   = zeros(length(iterations))
	its = zeros(length(iterations))
	status = [:NotSolved]

	for i in eachindex(iterations)

		if slv == SCSSolver
			solver = slv(verbose = verbose, max_iters = iterations[i], eps = 1e-20)
		elseif slv == ECOSSolver
			solver = slv(verbose = verbose, maxit = iterations[i])
		else
			solver = slv(verbose = verbose, maxit = iterations[i], tol = 1e-20)
		end

		b = @benchmarkable((x,it,status) = solve_problem(solver, setup..., setup_JuMP...),
				   setup = (
					setup  = deepcopy($setup);
					setup_JuMP  = deepcopy($setup_JuMP);
					solver = deepcopy($solver);
					x = nothing;
					status = nothing;
					it = 0;
					),
				   teardown = (
					   $err[$i] = norm(x-$xopt)/norm($xopt);
					   $status[1] = status;
					   ),
				   evals = 1, samples = samples, seconds = seconds)

		results = run(b)
		t[i] = time(median(results))
		its[i] = iterations[i]
		if err[i] <= eps_max || status[1] == :Optimal
			t = t[1:i]
			err = err[1:i]
			its = its[1:i]
			break
		end
	end

	return t, err, its
end
BLAS.set_num_threads(5)

nvars = 100000
solvers = [
#	   "SCSSolver",
#  	   "ECOSSolver",
	   "PG",
	   "FPG",
	   "ZeroFPR",
	   "PANOC"
	   ]
T, ERR = [],[]
for slv in solvers
	t, err = benchmark_LASSO(slv, 1e-8, nvars)
	push!(T,t)
	push!(ERR,err)
end

using PyPlot
using DataFrames, CSV
save_stuff = true
dirpath = "/home/nantonel/Proximal_Gradient_Algorithms/fig"
#mkpath("data/")
figure()
for i in eachindex(solvers)
	plot(1e-9.*T[i], log10.(ERR[i]),  label = solvers[i],":*")
    if save_stuff
        CSV.write("$dirpath/lasso_$(nvars)_$(solvers[i]).cvs",
                  DataFrame(q = log10.(ERR[i]), t = 1e-9.*T[i]) )
    end
end
legend()
ylabel("time")
xlabel("solution quality")


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f01568ca-a480-11eb-20c5-039d52fd7c70
begin
	using StatsPlots
	using PlutoUI
end

# ╔═╡ 72bc5d01-1a22-4837-af46-db56c02406e6
md"# Joint, Marginal, and Conditional Densities"

# ╔═╡ a010879a-571a-4098-8c7d-5880044c5e6b
md"Consider the **joint** density

$$f(x, y) = 4x^2y + 2y^5, \quad \text{ for } y \in [0, 1] \text{ and } x \in [0, 1].$$

We can visualize this as"

# ╔═╡ 008cc8d7-8a36-4643-a282-c68fb849da55
md"""The **marginal** density with respect to to $y$ is give by

$$f(y) = \frac{4}{3}y + 2y^5,$$

and the **conditional** density of $y$ given $x$ is

$$g(y|x) = \frac{4x^2y + 2y^5}{2x^2 + 1/3}.$$

Given some value of $x$, the conditional density $g(y|x)$ will change, but since the marginal density $f(y)$ has no dependence on $x$ it does not change as $x$ changes.  For the marginal density $f(y)$, we've already integrated over $x$ to remove the dependence on it.
"""

# ╔═╡ 4c36ebbe-8175-43f5-84f3-5a5aa046223e
md"Notice that at $x = 0.6$, $g(y|x) \approx f(y)$, and for all other values of $x$ the densities are clearly different."

# ╔═╡ 369567b4-b4a7-4eae-b5ea-43236f7375da
md"Adjust the value of $x$ with the slider below."

# ╔═╡ ac44d408-3ee0-4bce-aa68-92cb786934de
md"``x=`` $(@bind x Slider(0:0.1:1, default = 0.8, show_value = true))"

# ╔═╡ 47866c5f-b961-4db0-a5cd-0eb26ccf3359
md"To calculate the probabilities above, you need to perform integration, since these are continuous distributions:

$$\mathbb{P}[Y \in [0.2, 0.3] ] = \int_{[0.2, 0.3]} f(y) dy$$

and

$$\mathbb{P}[Y\in [0.2, 0.3] | X = x] = \int_{[0.2, 0.3]}  g(y|x) dy$$
"

# ╔═╡ a9b83b43-38f5-42ff-aa3c-b345ecf7fc3e
y = 0:0.05:1;

# ╔═╡ ed32a057-c61d-4dae-ad51-9b8a0b1ecd43
plot(y, y, (x, y) -> 4x^2 * y + 2y^5,
	xlabel = "x", ylabel = "y", zlabel = "f(x, y)",	t = :wireframe)

# ╔═╡ 393f86db-3259-45f9-84da-6205d5d2c928
f(y) = 4y/3 + 2y^5;

# ╔═╡ 60ab4f05-5819-4b48-9265-c531d4359153
g(y) = (4x^2 * y + 2y^5) / (2x^2 + 1/3);

# ╔═╡ 4398444f-1708-4755-8104-12e206b6640e
begin
	yfillmin = 0.2
	yfillmax = 0.3
	yfill = 0.2:0.05:0.3
end;

# ╔═╡ fc5259cb-1c3c-476f-868c-7ef244d9a7a7
begin
	F(y) = (2y^2 + y^6) / 3
	G(y) = (2x^2 * y^2 + y^6 / 3) / (2x^2 + 1/3)
end;

# ╔═╡ 2fc3983e-2b6c-42c6-8e44-b6178ad2e42b
begin
	vF = round(F(yfillmax) - F(yfillmin); digits = 3)
	vG = round(G(yfillmax) - G(yfillmin); digits = 3)
end;

# ╔═╡ 9e438242-5426-4dbc-af3f-c9f07b15b154
md"The probabilities ``\mathbb{P}[Y \in [0.2, 0.3] ] =`` $vF and ``\mathbb{P}[Y\in [0.2, 0.3] | X = x] =`` $vG are shown in the plot below.  The two shades of grey are meant to aid visualization."

# ╔═╡ c7a8192e-a2f9-4192-92e7-fbe316e946ac
begin
	plot(y, f.(y), label = "f(y)", color = :blue, ylabel = "Density", xlabel = "y")
	plot!(y, g.(y), label = "g(y|x = $(x))", color = :orange, legend = :topleft)
	if vF > vG
		plot!(yfill, f.(yfill), fill = (0, :grey90), label = "", seriescolor = :blue)
		plot!(yfill, g.(yfill), fill = (0, :grey70), label = "", seriescolor = :orange)
	else
		plot!(yfill, g.(yfill), fill = (0, :grey90), label = "", seriescolor = :orange)
		plot!(yfill, f.(yfill), fill = (0, :grey70), label = "", seriescolor = :blue)
	end
end

# ╔═╡ Cell order:
# ╟─72bc5d01-1a22-4837-af46-db56c02406e6
# ╟─a010879a-571a-4098-8c7d-5880044c5e6b
# ╟─ed32a057-c61d-4dae-ad51-9b8a0b1ecd43
# ╟─008cc8d7-8a36-4643-a282-c68fb849da55
# ╟─4c36ebbe-8175-43f5-84f3-5a5aa046223e
# ╟─9e438242-5426-4dbc-af3f-c9f07b15b154
# ╟─369567b4-b4a7-4eae-b5ea-43236f7375da
# ╟─ac44d408-3ee0-4bce-aa68-92cb786934de
# ╟─c7a8192e-a2f9-4192-92e7-fbe316e946ac
# ╟─47866c5f-b961-4db0-a5cd-0eb26ccf3359
# ╟─a9b83b43-38f5-42ff-aa3c-b345ecf7fc3e
# ╟─393f86db-3259-45f9-84da-6205d5d2c928
# ╟─60ab4f05-5819-4b48-9265-c531d4359153
# ╟─4398444f-1708-4755-8104-12e206b6640e
# ╟─fc5259cb-1c3c-476f-868c-7ef244d9a7a7
# ╟─2fc3983e-2b6c-42c6-8e44-b6178ad2e42b
# ╟─f01568ca-a480-11eb-20c5-039d52fd7c70

max_stars_repo_name: tim-au/pluto

max_issues_repo_name: tim-au/pluto

max_forks_repo_name: tim-au/pluto


text:
# SynchronySDK integration file

using SynchronySDK
using DocStringExtensions # temporary while $(SIGNATURES) is in use in this file

const Syncr = SynchronySDK

"""
    SyncrSLAM

An object definition containing the require variables to leverage the server side SLAM solution per user, robot, and session.
"""
mutable struct SyncrSLAM
  robotId::AbstractString
  sessionId::AbstractString
  syncrconf::Union{Void, SynchronySDK.SynchronyConfig}
  robot
  session

  # Constructor
  SyncrSLAM(
    robotId::AbstractString,
    sessionId::AbstractString,
    syncrconf::Union{Void, SynchronySDK.SynchronyConfig};
    robot=nothing,
    session=nothing
   ) = new(
        robotId,
        sessionId,
        syncrconf,
        robot,
        session
      )
end

"""
$(SIGNATURES)

Get  Synchrony configuration, default filepath location assumed as `~/Documents/synchronyConfig.json`.
"""
function loadSyncrConfig(;
            filepath::AS=joinpath(ENV["HOME"],"Documents","synchronyConfig.json")
         ) where {AS <: AbstractString}
  #
  println(" - Retrieving Synchrony Configuration...")
  configFile = open(filepath)
  configData = JSON.parse(readstring(configFile))
  close(configFile)
  Unmarshal.unmarshal(SynchronyConfig, configData) # synchronyConfig =
end

"""
$(SIGNATURES)

Confirm that the robot already exists, create if it doesn't.
"""
function syncrRobot(
           synchronyConfig::SynchronyConfig,
           robotId::AS
         ) where {AS <: AbstractString}
  #
  println(" - Creating or retrieving robot '$robotId'...")
  robot = nothing
  if(SynchronySDK.isRobotExisting(synchronyConfig, robotId))
      println(" -- Robot '$robotId' already exists, retrieving it...")
      robot = getRobot(synchronyConfig, robotId)
  else
      # Create a new one
      println(" -- Robot '$robotId' doesn't exist, creating it...")
      newRobot = RobotRequest(robotId, "My New Bot", "Description of my neat robot", "Active")
      robot = addRobot(synchronyConfig, newRobot)
  end
  println(robot)
  robot
end

"""
$(SIGNATURES)

Get sessions, if it already exists, add to it.
"""
function syncrSession(
            synchronyConfig::SynchronyConfig,
            robotId::AS,
            sessionId::AS
         ) where {AS<: AbstractString}
  #
  println(" - Creating or retrieving data session '$sessionId' for robot...")
  session = nothing
  if(SynchronySDK.isSessionExisting(synchronyConfig, robotId, sessionId))
      println(" -- Session '$sessionId' already exists for robot '$robotId', retrieving it...")
      session = getSession(synchronyConfig, robotId, sessionId)
  else
      # Create a new one
      println(" -- Session '$sessionId' doesn't exist for robot '$robotId', creating it...")
      newSessionRequest = SessionDetailsRequest(sessionId, "Submersible vehicle.", "Pose3", false)
      session = addSession(synchronyConfig, robotId, newSessionRequest)
  end
  println(session)
  session
end

"""
$(SIGNATURES)

Intialize the `sslaml` object using configuration file defined in `syncrconfpath`.
"""
function initialize!(sslaml::SyncrSLAM;
            syncrconfpath::AS=joinpath(ENV["HOME"],"Documents","synchronyConfig.json")
         ) where {AS <: AbstractString}
  # 1. Get a Synchrony configuration
  sslaml.syncrconf = loadSyncrConfig(filepath=syncrconfpath)

  # 2. Confirm that the robot already exists, create if it doesn't.
  sslaml.robot = syncrRobot(sslaml.syncrconf, sslaml.robotId)

  # Make sure that this session is not already populated
  if isSessionExisting(slam_client.syncrconf, sslaml.robotId, sslaml.sessionId)
      warn("There is already a session named '$sessionId' for robot '$robotId'. This example will fail if it tries to add duplicate nodes. We strongly recommend providing a new session name.")
      print(" Should we delete it? [Y/N] - ")
      if lowercase(strip(readline(STDIN))) == "y"
          deleteSession(sslaml.syncrconf, sslaml.robotId, sslaml.sessionId)
          println(" -- Session '$sessionId' deleted!")
      else
          warn("Okay, but this may break the example! Caveat Userus!")
      end
  end

  # 3. Create or retrieve the session.
  sslaml.session = syncrSession(sslaml.syncrconf, sslaml.robotId, sslaml.sessionId)

  nothing
end


text:
using DiscreteMarkovChains
using Test

@testset "Expected Failures" begin
    # Different lengths
    @test_throws(
        ErrorException(
            "The state space, [1, 2], and "*
            "transition matrix should be the same size."
        ),
        DiscreteMarkovChain([1, 2], [1]),
    )

    # Repeated elements in the state space
    @test_throws(
        ErrorException("The state space, [1, 1], must have unique elements."),
        DiscreteMarkovChain([1, 1], [1 0; 0 1]),
    )

    # Not row stochastic
    @test_throws(
        ErrorException(
            "The transition matrix, [0.5 0.6; 1.0 0.0], "*
            "should be row-stochastic (each row must sum up to 1)."
        ),
        DiscreteMarkovChain([1, 2], [0.5 0.6; 1.0 0.0]),
    )
end

@testset "Communication Classes" begin
    # Empty test
    T = [][:,:]
    X = DiscreteMarkovChain([], T)
    @test periodicities(X) == ([], [], [])

    state_space = [1, 2, 3]
    T = [
        0 1 0;
        1 0 0;
        0 0 1;
    ]
    X = DiscreteMarkovChain(state_space, T)
    @test periodicities(X) == ([[1, 2], [3]], [true, true], [2, 1])

    T = [
        0 5 5 0 0;
        0 0 0 10 0;
        5 0 5 0 0;
        0 10 0 0 0;
        0 3 0 3 4;
    ]/10
    X = DiscreteMarkovChain(0:4, T)
    classes, recurrence, periods = periodicities(X)
    @test classes == [[1, 3], [0, 2], [4]]
    @test recurrence == [true, false, false]
    @test periods == [2, 1, 1]

    T = [
        0 0 0 10 0 0;
        5 0 5 0 0 0;
        0 4 0 0 0 6;
        10 0 0 0 0 0;
        0 10 0 0 0 0;
        0 0 0 5 5 0;
    ]/10
    X = DiscreteMarkovChain(0:5, T)
    classes, recurrence, periods = periodicities(X)
    @test classes == [[0, 3], [1, 2, 5, 4]]
    @test recurrence == [true, false]
    @test periods == [2, 2]

    T = [
        2 0 0 3 0 0 3 2 0 0;
        0 10 0 0 0 0 0 0 0 0;
        0 2 2 0 0 0 0 0 3 3;
        0 0 0 3 0 0 6 1 0 0;
        0 0 0 0 5 5 0 0 0 0;
        0 0 0 0 0 10 0 0 0 0;
        4 0 0 5 0 0 1 0 0 0;
        2 0 0 4 0 0 2 2 0 0;
        3 0 1 0 0 0 0 0 4 2;
        0 0 4 0 0 0 0 0 3 3;
    ]/10
    X = DiscreteMarkovChain(0:9, T)
    classes, recurrence, periods = periodicities(X)
    @test classes == [[0, 3, 6, 7], [1], [2, 8, 9], [5], [4]]
    @test recurrence == [true, true, false, true, false]
    @test periods == [1, 1, 1, 1, 1]
end

@testset "Decompose and Canonical Form" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test decompose(X) == ([], T, T, T)

    T = [
        1 0 0;
        1//3 1//3 1//3;
        0 1//4 3//4;
    ]
    X = DiscreteMarkovChain(["Sunny", "Cloudy", "Rainy"], T)
    @test decompose(X) == (
        ["Sunny", "Cloudy", "Rainy"],
        T[1:1, 1:1], T[2:3, 1:1], T[2:3, 2:3],
    )
    @test canonical_form(X) == (["Sunny", "Cloudy", "Rainy"], T)

    T = [
        1//4 3//4 0;
        1//3 1//3 1//3;
        0 1//4 3//4;
    ]
    X = DiscreteMarkovChain([1, 2, 3], T)
    @test decompose(X) == ([1, 2, 3], T, Array{Any}(undef, 0, 3), Array{Any}(undef, 0, 0))
    @test canonical_form(X) == ([1, 2, 3], T)

    T = [
        1 0 0 0 0;
        1//2 0 1//2 0 0;
        0 1//2 0 1//2 0;
        0 0 1//2 0 1//2;
        0 0 0 0 1;
    ]
    X = DiscreteMarkovChain(0:4, T)
    states, A, B, C = decompose(X)
    @test states == [0, 4, 1, 2, 3]
    @test A == [1 0; 0 1]
    @test B == [1//2 0; 0 0; 0 1//2]
    @test C == [0 1//2 0; 1//2 0 1//2; 0 1//2 0]
    states, new_matrix = canonical_form(X)
    @test states == [0, 4, 1, 2, 3]
    desired_new_matrix = [
        1 0 0 0 0;
        0 1 0 0 0;
        1//2 0 0 1//2 0;
        0 0 1//2 0 1//2;
        0 1//2 0 1//2 0;
    ]
    @test new_matrix == desired_new_matrix
end

@testset "Regular, Ergodic, Absorbing" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test !is_regular(X)
    @test !is_ergodic(X)
    @test !is_absorbing(X)

    T = [
        0 4 0 0 0;
        1 0 3 0 0;
        0 2 0 2 0;
        0 0 3 0 1;
        0 0 0 4 0;
    ]/4.0
    X = DiscreteMarkovChain(1:5, T)
    @test !is_regular(X)
    @test is_ergodic(X)
    @test !is_absorbing(X)

    T = [
        0 1;
        1 0;
    ]
    X = DiscreteMarkovChain(1:2, T)
    @test !is_regular(X)
    @test is_ergodic(X)
    @test !is_absorbing(X)

    T = [
        2 1 1;
        2 0 2;
        1 1 2;
    ]/4
    X = DiscreteMarkovChain(1:3, T)
    @test is_regular(X)
    @test is_ergodic(X)
    @test !is_absorbing(X)

    T = [
        1 1;
        1 1;
    ]/2
    X = DiscreteMarkovChain(1:2, T)
    @test is_regular(X)
    @test is_ergodic(X)
    @test !is_absorbing(X)
end

@testset "Stationary Distribution" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test stationary_distribution(X) == Any[]

    T = [
        1//5 2//5 2//5;
        1//10 1//2 2//5;
        3//5 3//10 1//10;
    ]
    X = DiscreteMarkovChain(1:3, T)
    @test stationary_distribution(X) == [11//39, 16//39, 4//13]

    T = [
        0 1 0 0 0;
        .25 0 .75 0 0;
        0 .5 0 .5 0;
        0 0 .75 0 .25;
        0 0 0 1 0;
    ]
    X = DiscreteMarkovChain(0:4, T)
    @test stationary_distribution(X) ≈ [.0625, .2500, .3750, .2500, .0625]
end

@testset "Fundamental Matrix" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test fundamental_matrix(X) == Array{Any}(undef, 0, 0)

    T = [
        1 0 0 0 0;
        1//2 0 1//2 0 0;
        0 1//2 0 1//2 0;
        0 0 1//2 0 1//2;
        0 0 0 0 1;
    ]
    X = DiscreteMarkovChain([0, 1, 2, 3, 4], T)
    desired_new_matrix = [
        3//2 1 1//2;
        1 2 1;
        1//2 1 3//2;
    ]
    @test fundamental_matrix(X) == desired_new_matrix
end

@testset "Expected Time To Absorption" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test mean_time_to_absorption(X) == Any[]

    T = [
        1 0 0 0 0;
        1//2 0 1//2 0 0;
        0 1//2 0 1//2 0;
        0 0 1//2 0 1//2;
        0 0 0 0 1;
    ]
    X = DiscreteMarkovChain([0, 1, 2, 3, 4], T)
    @test mean_time_to_absorption(X) == [2, 3, 2]

    T = [
        1 0 0 0 0;
        0 1 0 0 0;
        0 0 0 1 0;
        1//4 1//4 1//4 0 1//4;
        1//2 0 0 1//2 0;
    ]
    X = DiscreteMarkovChain([0, 1, 2, 3, 4], T)
    @test mean_time_to_absorption(X) == [12//5, 7//5, 6//5]
end

@testset "Exit Probability" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test exit_probabilities(X) == Array{Any}(undef, 0, 0)

    T = [
        1 0 0 0 0;
        0 1 0 0 0;
        0 0 0 1 0;
        1//4 1//4 1//4 0 1//4;
        1//2 0 0 1//2 0;
    ]
    X = DiscreteMarkovChain([0, 1, 2, 3, 4], T)
    desired_new_matrix = [
        3//5 2//5;
        3//5 2//5;
        4//5 1//5;
    ]
    @test exit_probabilities(X) == desired_new_matrix

    T = [
        2//5 3//5 0 0 0 0;
        9//10 1//10 0 0 0 0;
        0 0 1 0 0 0;
        0 0 0 0 1 0;
        0 1//4 1//4 1//4 0 1//4;
        0 1//2 0 0 1//2 0;
    ]
    X = DiscreteMarkovChain([0, 1, 2, 3, 4, 5], T)
    desired_new_matrix = [
        0 6//10 4//10;
        0 6//10 4//10;
        0 8//10 2//10;
    ]
    @test exit_probabilities(X) == desired_new_matrix
end

@testset "First Passage" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test first_passage_probabilities(X, 1) == Array{Any}(undef, 0, 0)

    T = [
        2//10 4//10 4//10;
        3//10 3//10 4//10;
        5//10 4//10 1//10;
    ]
    X = DiscreteMarkovChain(["A", "B", "C"], T)
    @test first_passage_probabilities(X, 1) == T

    @test first_passage_probabilities(X, 1, "A", "C") == 4//10
    @test first_passage_probabilities(X, 2, "A", "C") == 24//100
    @test first_passage_probabilities(X, 3, "A", "C") == 144//1000
    @test first_passage_probabilities(X, 4, "A", "C") == 864//10000

    @test first_passage_probabilities(X, 1, "A", "A") == 2//10
    @test first_passage_probabilities(X, 2, "A", "A") == 32//100
    @test first_passage_probabilities(X, 3, "A", "A") == 184//1000
    @test first_passage_probabilities(X, 4, "A", "A") == 1152//10000
end

@testset "Mean Recurrence Time" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test mean_recurrence_time(X) == Any[]

    T = [
        0 1 0;
        0 0 1;
        1 0 0;
    ]
    X = DiscreteMarkovChain(T)
    @test mean_recurrence_time(X) ≈ [3, 3, 3]

    T = [
        0 1 0 0 0;
        .25 0 .75 0 0;
        0 .5 0 .5 0;
        0 0 .75 0 .25;
        0 0 0 1 0;
    ]
    X = DiscreteMarkovChain(0:4, T)
    @test mean_recurrence_time(X) ≈ [16.0000, 4.0000, 8/3, 4.0000, 16.000]
end

@testset "Mean First Passage Time" begin
    # Empty test
    T = Array{Any}(undef, 0, 0)
    X = DiscreteMarkovChain([], T)
    @test mean_first_passage_time(X) == Array{Any}(undef, 0, 0)

    T = [
        1//2 1//4 1//4;
        1//2 0 1//2;
        1//4 1//4 1//2;
    ]
    X = DiscreteMarkovChain(["A", "B", "C"], T)
    desired_new_matrix = [
        0 4 10//3;
        8//3 0 8//3;
        10//3 4 0;
    ]
    @test mean_first_passage_time(X) == desired_new_matrix

    T = [
        0 1 0 0 0;
        .25 0 .75 0 0;
        0 .5 0 .5 0;
        0 0 .75 0 .25;
        0 0 0 1 0;
    ]
    X = DiscreteMarkovChain(0:4, T)
    desired_new_matrix = [
        0 1 8/3 19/3 64/3;
        15 0 5/3 16/3 61/3;
        56/3 11/3 0 11/3 56/3;
        61/3 16/3 5/3 0 15;
        64/3 19/3 8/3 1 0;
    ]
    @test mean_first_passage_time(X) ≈ desired_new_matrix
end


text:
abstract type PDXObject end

mutable struct PDFormXObject <: PDXObject
    doc::PDDoc
    cosXObj::CosIndirectObject{CosStream}
    matrix::Matrix{Float32}
    bbox::CDRect{Float32}
    fonts::Dict{CosName, PDFont}
    xobjs::Dict{CosName, PDXObject}
    content_objects::PDPageObjectGroup
    function PDFormXObject(doc::PDDoc, cosxobj::CosIndirectObject{CosStream})
        mat = get(cosxobj, cn"Matrix")
        box = get(cosxobj, cn"BBox")
        @assert box !== CosNull "Invalid Form XObject without bounding box"
        matrix = mat === CosNull ?
            [1f0 0f0 0f0; 0f0 1f0 0f0; 0f0 0f0 1f0] :
            hcat(reshape((get.(CosFloat.(get(mat)))), (2, 3))', [0f0, 0f0, 1f0])
        bbox = CDRect{Float32}(CDRect(box))
        fonts = Dict{CosName, PDFont}()
        xobjs = Dict{CosName, PDXObject}()
        new(doc, cosxobj, matrix, bbox, fonts, xobjs, PDPageObjectGroup())
    end
end

mutable struct PDImageXObject <: PDXObject
    doc::PDDoc
    obj::CosIndirectObject{CosStream}
end

mutable struct PDDefaultXObject <: PDXObject
    doc::PDDoc
    obj::CosObject
end

function createPDXObject(doc::PDDoc, cosstm::CosObject)
    otype = get(cosstm, cn"Type")
    @assert otype === cn"XObject" || otype === CosNull
    subtype = get(cosstm, cn"Subtype")
    subtype === cn"Form"  && return PDFormXObject(doc, cosstm)
    subtype === cn"Image" && return PDImageXObject(doc, cosstm)
    return PDDefaultXObject(doc, cosobj)
end

function find_resource(xobj::PDFormXObject,
                       restype::CosName,
                       resname::CosName)
    cosdoc = xobj.doc.cosDoc
    resref = get(xobj.cosXObj, cn"Resources")
    resref === CosNull && return CosNull
    resources = cosDocGetObject(cosdoc, resref)
    resources === CosNull && return CosNull
    ress = cosDocGetObject(cosdoc, resources, restype)
    ress === CosNull && return CosNull
    res = cosDocGetObject(cosdoc, ress, resname)
    return res
end

get_font(xobj::PDXObject, fontname::CosName) =
    get!(xobj.fonts, fontname,
         get_pd_font!(xobj.doc, find_resource(xobj, cn"Font", fontname)))

get_xobject(xobj::PDXObject, xobjname::CosName) =
    get!(xobj.xobjs, xobjname,
         get_pd_xobject!(xobj.doc,
                         find_resource(xobj, cn"XObject", xobjname)))

function load_content_objects(xobj::PDFormXObject)
    stm = xobj.cosXObj
    bufstm = decode(stm)
    # try
        load_objects(xobj.content_objects, bufstm)
    # finally
    #    util_close(bufstm)
    # end
    return nothing
end

Do(xobj::PDDefaultXObject, state::GState) = nothing

Do(xobj::PDImageXObject, state::GState) = nothing

function Do(xobj::PDFormXObject, state::GState)
    isempty(xobj.content_objects) && load_content_objects(xobj)
    isempty(xobj.content_objects) && return state
    xstate = GState{:PDFIO}()
    ctm = state[:CTM]
    nctm = xobj.matrix*ctm
    xstate[:CTM] = nctm
    xstate[:source] = xobj
    xstate[:text_layout] = state[:text_layout]
    xstate[:h_profile] = state[:h_profile]
    evalContent!(xobj.content_objects, xstate)
    return state
end


text:
using DataDeps

register(DataDep(
  "wall-following-robot-navigation-24",
  "http://archive.ics.uci.edu/ml/datasets/Wall-Following+Robot+Navigation+Data",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00194/sensor_readings_24.data",
  "6e9728dedef8b35e31ba0be11aea85aa53bd248ab33652777df81b95a4cc8c00",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=25,
        feature_indices=1:24,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{L}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {L <: PowerSystems.ElectricLoad, D <: AbstractControllableLoadForm, S <: PM.AbstractPowerFormulation}

    dev_set = [a.second for a in args if a.first == :devices]

    isempty(dev_set) ? devices = [d for d in sys.loads if (d.available == true && !isa(d,PowerSystems.StaticLoad))] : devices = dev_set[1]

    p_cl = activepowervariables(m, devices, sys.time_periods);

    varnetinjectiterate!(netinjection.var_active, p_cl, sys.time_periods, devices)

    activepower(m, devices, category_formulation, system_formulation, sys.time_periods)

    cost = variablecost(m, devices, category_formulation, system_formulation)

    add_to_cost!(m, cost)

end

function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{L}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {L <: PowerSystems.ElectricLoad, D <: AbstractControllableLoadForm, S <: AbstractACPowerModel}

    dev_set = [d for d in sys.loads if (d.available == true && !isa(d,PowerSystems.StaticLoad))]

    constructdevice!(m, netinjection, category, category_formulation, PM.AbstractPowerFormulation, sys; devices = dev_set)

    q_cl = reactivepowervariables(m, dev_set, sys.time_periods);

    varnetinjectiterate!(netinjection.var_reactive, q_cl, sys.time_periods, dev_set)

    reactivepower(m, dev_set, category_formulation, system_formulation, sys.time_periods)

end


text:
# ---
# title: 712. Minimum ASCII Delete Sum for Two Strings
# id: problem712
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-ascii-delete-sum-for-two-strings/description/>
# hidden: true
# ---
#
# Given two strings `s1, s2`, find the lowest ASCII sum of deleted characters to
# make two strings equal.
#
# **Example 1:**
#
#
#
#     Input: s1 = "sea", s2 = "eat"
#     Output: 231
#     Explanation: Deleting "s" from "sea" adds the ASCII value of "s" (115) to the sum.
#     Deleting "t" from "eat" adds 116 to the sum.
#     At the end, both strings are equal, and 115 + 116 = 231 is the minimum sum possible to achieve this.
#
#
# **Example 2:**
#
#
#
#     Input: s1 = "delete", s2 = "leet"
#     Output: 403
#     Explanation: Deleting "dee" from "delete" to turn the string into "let",
#     adds 100[d]+101[e]+101[e] to the sum.  Deleting "e" from "leet" adds 101[e] to the sum.
#     At the end, both strings are equal to "let", and the answer is 100+101+101+101 = 403.
#     If instead we turned both strings into "lee" or "eet", we would get answers of 433 or 417, which are higher.
#
#
# **Note:**
#
# * `0 < s1.length, s2.length <= 1000`.
# * All elements of each string will have an ASCII value in `[97, 122]`.
#
#
## @lc code=start
using LeetCode

function minimum_delete_sum(s1::String, s2::String)
    m, n = length(s1) + 1, length(s2) + 1
    dp = fill(0, m, n)
    dp[2:end, 1] .= cumsum(codeunits(s1))
    dp[1, 2:end] .= cumsum(codeunits(s2))
    for i in 2:m, j in 2:n
        dp[i, j] = (s1[i - 1] == s2[j - 1]) ? (dp[i - 1, j - 1]) :
                   min(dp[i - 1, j] + Int(s1[i - 1]), dp[i, j - 1] + Int(s2[j - 1]))
    end
    return dp[m, n]
end
## @lc code=end


text:
#=
Popis
  Funkcia vrati podgraf zo zadaneho grafu na zaklade indexov.

Parametre
  verts: pole indexov vrcholov ::Array{Int64, 1}
  g: datova struktura grafu :: GenericGraph

Navratova hodnota
  sub_graph: podgraf grafu g :: GenericGraph

Poznamky
  Funkcia vytvori novy podgraf zo zadaneho grafu.
=#

function get_subgraph(verts::Array{Int64, 1}, g::GenericGraph)

  sub_graph = Graphs.graph(Graphs.ExVertex[], Graphs.ExEdge{Graphs.ExVertex}[], is_directed=false)

  for v in sort(verts)
    add_vertex!(sub_graph, get_vertex(v, g))
  end

  for e in edges(g)
    if in(source(e).index, verts) && in(target(e).index, verts)

      edge_sub = add_edge!(sub_graph, get_vertex(source(e).index, g), get_vertex(target(e).index, g))

      for (att,value) in e.attributes
        edge_sub.attributes[att] = value
      end

    end
  end

  return sub_graph

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
using Gtk
params = IgusRobotParams()
rob = IgusRobot("igusRob",params)

if ask_dialog("You are trying to start the Igus robot hardware test. Please ensure that the robot is turned on, is able to reference itself and can move at least 50mm in the positive direction. Please continue only if the robot is safe to move!","Cancel","Start test")
    if isReferenced(rob)
        @test ask_dialog("The robot says it is already referenced. This usually means that the robot has not been powered off since the last reference drive. Is this the case?")
    end

    @test state(rob)==:INIT
    @test dof(rob)==1
    @test namedPositions(rob)["origin"]==[0]u"mm"
    @test collect(keys(namedPositions(rob))) == ["origin"]

    setup(rob)
    @test state(rob)==:DISABLED

    @test_throws RobotStateError moveAbs(rob, 1u"mm")
    @test_throws RobotStateError moveAbs(rob, [1]u"mm")
    @test_throws RobotStateError setup(rob)

    enable(rob)
    @test state(rob)==:READY
    moveRel(rob, [10]u"mm")
    @test ask_dialog("Did the robot move 10mm in the positive direction?")
    moveRel(rob, 10u"mm")
    @test ask_dialog("Did the robot move 10mm in the positive direction?")

    @test_throws RobotAxisRangeError moveRel(rob,10u"m")
    @test_throws RobotDOFError moveRel(rob,[10u"mm",10u"mm"])
    @test_throws RobotDOFError moveAbs(rob,[10u"mm",10u"mm"])

    doReferenceDrive(rob)
    @test isReferenced(rob)

    moveAbs(rob, [20]u"mm")
    @test ask_dialog("Did the robot move to x=20mm?")
    @test getPosition(rob)==[20]u"mm"

    moveAbs(rob, 50u"mm")
    @test ask_dialog("Did the robot move to x=50mm?")
    @test getPosition(rob)==[50]u"mm"

    teachNamedPosition(rob, "pos1")
    @test issetequal(keys(namedPositions(rob)), ["origin", "pos1"])
    @test namedPositions(rob)["pos1"] == [50]u"mm"
    @test_throws RobotTeachError teachNamedPosition(rob, "pos1")

    moveAbs(rob, [10u"mm"],20u"mm/s")
    @test ask_dialog("Did the robot move faster to x=10mm?")
    @test getPosition(rob)==[10]u"mm"

    gotoPos(rob, "pos1")
    @test ask_dialog("Did the robot move to x=50mm?")
    @test getPosition(rob)==[50]u"mm"

    gotoPos(rob, "origin", 20u"mm/s")
    @test ask_dialog("Did the robot move faster to x=0mm?")

    moveAbs(rob, [50u"mm"],[20u"mm/s"])
    @test ask_dialog("Did the robot move faster to x=20mm?")
    @test getPosition(rob)==[50]u"mm"

    reset(rob)
    @test state(rob)==:INIT
    setup(rob)
    enable(rob)

    @test isReferenced(rob) #resetting of the igus robot should not lose the referencing

    disable(rob)
    @test ask_dialog("Did everything seem normal?")
else
    @test_broken 0
end


text:
module BinaryPlatformTests
import ..Pkg # ensure we are using the correct Pkg

using Test, Pkg.BinaryPlatforms
import Pkg.BinaryPlatforms: platform_name

# The platform we're running on
const platform = @inferred Platform platform_key_abi()

# This is a compatability test; once we've fully migrated away from Pkg.BinaryPlatforms
# to the new Base.BinaryPlatforms module, we can throw away the shim definitions in
# `BinaryPlatforms_compat.jl` and drop these tests.
@testset "Compat - PlatformNames" begin
    # Ensure the platform type constructors are well behaved
    @testset "Platform constructors" begin
        @test_throws ArgumentError Linux(:not_a_platform)
        @test_throws ArgumentError Linux(:x86_64; libc=:crazy_libc)
        @test_throws ArgumentError Linux(:x86_64; libc=:glibc, call_abi=:crazy_abi)
        @test_throws ArgumentError Linux(:x86_64; libc=:glibc, call_abi=:eabihf)
        @test_throws ArgumentError Linux(:armv7l; libc=:glibc, call_abi=:kekeke)
        @test_throws ArgumentError MacOS(:i686)
        @test_throws ArgumentError MacOS(:x86_64; libc=:glibc)
        @test_throws ArgumentError MacOS(:x86_64; call_abi=:eabihf)
        @test_throws ArgumentError Windows(:x86_64; libc=:glibc)
        @test_throws ArgumentError Windows(:x86_64; call_abi=:eabihf)
        @test_throws ArgumentError FreeBSD(:not_a_platform)
        @test_throws ArgumentError FreeBSD(:x86_64; libc=:crazy_libc)
        @test_throws ArgumentError FreeBSD(:x86_64; call_abi=:crazy_abi)
        @test_throws ArgumentError FreeBSD(:x86_64; call_abi=:eabihf)

        # Test copy constructor
        cabi = CompilerABI(;
            libgfortran_version=v"3",
            libstdcxx_version=v"3.4.18",
            cxxstring_abi=:cxx03,
        )
        cabi2 = CompilerABI(cabi; cxxstring_abi=:cxx11)
        @test libgfortran_version(cabi) == libgfortran_version(cabi2)
        @test libstdcxx_version(cabi) == libstdcxx_version(cabi2)
        @test cxxstring_abi(cabi) != cxxstring_abi(cabi2)

        # Explicitly test that we can pass arguments to UnknownPlatform,
        # and it doesn't do anything.
        @test UnknownPlatform(:riscv; libc=:fuschia_libc) == UnknownPlatform()
    end

    @testset "Platform properties" begin
        # Test that we can get the name of various platforms
        for T in (Linux, MacOS, Windows, FreeBSD)
            @test endswith(lowercase(string(T)), lowercase(platform_name(T(:x86_64))))
        end

        # Test that we can get the arch of various platforms
        @test arch(Linux(:aarch64; libc=:musl)) == :aarch64
        @test arch(Windows(:i686)) == :i686
        @test arch(FreeBSD(:amd64)) == :x86_64
        @test arch(FreeBSD(:i386)) == :i686
        @test arch(UnknownPlatform(:ppc64le)) == nothing

        # Test that our platform_dlext stuff works
        @test platform_dlext(Linux(:x86_64)) == platform_dlext(Linux(:i686))
        @test platform_dlext(Windows(:x86_64)) == platform_dlext(Windows(:i686))
        @test platform_dlext(MacOS()) != platform_dlext(Linux(:armv7l))
        @test platform_dlext(FreeBSD(:x86_64)) == platform_dlext(Linux(:x86_64))
        @test platform_dlext() == platform_dlext(platform)

        @test wordsize(Linux(:i686)) == wordsize(Linux(:armv7l)) == 32
        @test wordsize(MacOS()) == wordsize(Linux(:aarch64)) == 64
        @test wordsize(FreeBSD(:x86_64)) == wordsize(Linux(:powerpc64le)) == 64

        @test call_abi(Linux(:x86_64)) === nothing
        @test call_abi(Linux(:armv6l)) == :eabihf
        @test call_abi(Linux(:armv7l; call_abi=:eabihf)) == :eabihf
        @test call_abi(UnknownPlatform(;call_abi=:eabihf)) === nothing

        @test triplet(Windows(:i686)) == "i686-w64-mingw32"
        @test triplet(Linux(:x86_64; libc=:musl)) == "x86_64-linux-musl"
        @test triplet(Linux(:armv7l; libc=:musl)) == "armv7l-linux-musleabihf"
        @test triplet(Linux(:armv6l; libc=:musl, call_abi=:eabihf)) == "armv6l-linux-musleabihf"
        @test triplet(Linux(:x86_64)) == "x86_64-linux-gnu"
        @test triplet(Linux(:armv6l)) == "armv6l-linux-gnueabihf"
        @test triplet(MacOS()) == "x86_64-apple-darwin14"
        @test triplet(FreeBSD(:x86_64)) == "x86_64-unknown-freebsd11.1"
        @test triplet(FreeBSD(:i686)) == "i686-unknown-freebsd11.1"
    end

    @testset "Valid DL paths" begin
        # Test some valid dynamic library paths
        @test valid_dl_path("libfoo.so.1.2.3", Linux(:x86_64))
        @test valid_dl_path("libfoo.1.2.3.so", Linux(:x86_64))
        @test valid_dl_path("libfoo-1.2.3.dll", Windows(:x86_64))
        @test valid_dl_path("libfoo.1.2.3.dylib", MacOS())
        @test !valid_dl_path("libfoo.dylib", Linux(:x86_64))
        @test !valid_dl_path("libfoo.so", Windows(:x86_64))
        @test !valid_dl_path("libfoo.dll", MacOS())
        @test !valid_dl_path("libfoo.so.1.2.3.", Linux(:x86_64))
        @test !valid_dl_path("libfoo.so.1.2a.3", Linux(:x86_64))
    end

    @testset "platforms_match()" begin
        # Just do a quick combinatorial sweep for completeness' sake for platform matching
        for libgfortran_version in (nothing, v"3", v"5"),
            libstdcxx_version in (nothing, v"3.4.18", v"3.4.26"),
            cxxstring_abi in (nothing, :cxx03, :cxx11)

            cabi = CompilerABI(;
                libgfortran_version=libgfortran_version,
                libstdcxx_version=libstdcxx_version,
                cxxstring_abi=cxxstring_abi,
            )
            @test platforms_match(Linux(:x86_64), Linux(:x86_64, compiler_abi=cabi))
            @test platforms_match(Linux(:x86_64, compiler_abi=cabi), Linux(:x86_64))

            # Also test auto-string-parsing
            @test platforms_match(triplet(Linux(:x86_64)), Linux(:x86_64, compiler_abi=cabi))
            @test platforms_match(Linux(:x86_64), triplet(Linux(:x86_64, compiler_abi=cabi)))
        end

        # Ensure many of these things do NOT match
        @test !platforms_match(Linux(:x86_64), Linux(:i686))
        @test !platforms_match(Linux(:x86_64), Windows(:x86_64))
        @test !platforms_match(Linux(:x86_64), MacOS())
        @test !platforms_match(Linux(:x86_64), UnknownPlatform())

        # Make some explicitly non-matching cabi's
        base_cabi = CompilerABI(;
            libgfortran_version=v"5",
            cxxstring_abi=:cxx11,
        )
        for arch in (:x86_64, :i686, :aarch64, :armv6l, :armv7l),
            cabi in (
                CompilerABI(libgfortran_version=v"3"),
                CompilerABI(cxxstring_abi=:cxx03),
                CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11),
                CompilerABI(libgfortran_version=v"3", cxxstring_abi=:cxx03),
            )

            @test !platforms_match(Linux(arch, compiler_abi=base_cabi), Linux(arch, compiler_abi=cabi))
        end
    end

    @testset "Sys.is* overloading" begin
        # Test that we can indeed ask if something is linux or windows, etc...
        @test Sys.islinux(Linux(:aarch64))
        @test !Sys.islinux(Windows(:x86_64))
        @test Sys.iswindows(Windows(:i686))
        @test !Sys.iswindows(Linux(:x86_64))
        @test Sys.isapple(MacOS())
        @test !Sys.isapple(Linux(:powerpc64le))
        @test Sys.isbsd(MacOS())
        @test Sys.isbsd(FreeBSD(:x86_64))
        @test !Sys.isbsd(Linux(:powerpc64le; libc=:musl))
    end
end

@testset "select_platform" begin
    platforms = Dict(
        # Typical binning test
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "linux4",
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4")) => "linux7",
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"5")) => "linux8",

        # Ambiguity test
        Linux(:aarch64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "linux4",
        Linux(:aarch64, compiler_abi=CompilerABI(libgfortran_version=v"3", libstdcxx_version=v"3.4.18")) => "linux5",

        MacOS(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "mac4",
        Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)) => "win",
    )

    @test select_platform(platforms, Linux(:x86_64)) == "linux8"
    @test select_platform(platforms, Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4"))) == "linux7"

    # Ambiguity test
    @test select_platform(platforms, Linux(:aarch64)) == "linux5"
    @test select_platform(platforms, Linux(:aarch64; compiler_abi=CompilerABI(libgfortran_version=v"3"))) == "linux5"
    @test select_platform(platforms, Linux(:aarch64; compiler_abi=CompilerABI(libgfortran_version=v"4"))) == nothing

    @test select_platform(platforms, MacOS(:x86_64)) == "mac4"
    @test select_platform(platforms, MacOS(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4"))) == nothing

    @test select_platform(platforms, Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11))) == "win"
    @test select_platform(platforms, Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx03))) == nothing

    # Poor little guy
    @test select_platform(platforms, FreeBSD(:x86_64)) == nothing
end

end # module


text:
macro parse_opt_field_param(struct_variable, field_name, field_type)
    return quote
        if $(esc(field_name)) != nothing
            $(esc(struct_variable)).$(field_name) = $(field_type)($(esc(field_name)))
        end
    end
end

macro parse_opt_field_param(struct_variable, field_name)
    return quote
        if $(esc(field_name)) != nothing
            $(esc(struct_variable)).$(field_name) = $(esc(field_name))
        end
    end
end


text:
module TransposeMultithreading

using KernelAbstractions.Extras: @unroll
using SIMD: Vec, vload, vstore
using StaticArrays: MArray, MVector
using MyGEMM.Types: TilingParams
using MyGEMM.Utils: packA!, packB!

export multiply_transpose!

"""
  multiply_transpose!(C, A, B, params)

Performs `C += A * B'` with packing and multithreading option.
"""
@inline function multiply_transpose!(C, A, B, params::TilingParams)
  FT  = eltype(C)

  # Get the matrix params
  m, n = size(C)
  k = size(B, 2)

  # Check the matrix params
  @assert size(C) == (m, n)
  @assert size(A) == (m, k)
  @assert size(B) == (n, k)

  # Get tiling parameters
  mr = params.mr
  nr = params.nr
  vl = params.vl
  mc = params.mc
  nc = params.nc
  kc = params.kc

  # Since `B_packed` lives in the L3 cache we need to rescale it
  nc = cld(cld(min(n, nc), Threads.nthreads()), nr) * nr
  params = TilingParams{mc, nc, kc, mr, nr, vl}()

  # We assume that the matrices are even divisible by `nr` and `mr`,
  # that the blocking perfectly partitions the matrix,
  # and, vectors are perfectly divided by the chosen vector unit
  @assert mod(m, mr)  == 0
  @assert mod(n, nr)  == 0
  @assert mod(mc, mr) == 0
  @assert mod(nc, nr) == 0
  @assert mod(mr, vl) == 0

  # Prepare tuples for packing
  A_packed = ntuple(j -> Array{FT}(undef, mc * kc), Threads.nthreads())
  B_packed = ntuple(j -> Array{FT}(undef, kc * nc), Threads.nthreads())

  # Execute C = C + A * B'
  loop5!(C, A, B, A_packed, B_packed, params)
end

@inline function loop5!(C, A, B, A_packed, B_packed, params::TilingParams)
  n = size(C, 2)
  nc = params.nc
  @inbounds Threads.@threads for j1 = 1:nc:n
    j2 = min(n, j1 + nc - 1)
    C_j = @view C[:, j1:j2]
    B_j = @view B[j1:j2, :]
    loop4!(
      C_j,
      A,
      B_j,
      A_packed[Threads.threadid()],
      B_packed[Threads.threadid()],
      params
    )
  end
end

@inline function loop4!(C, A, B, A_packed, B_packed, params::TilingParams)
  k = size(A, 2)
  kc = params.kc
  mr = params.mr
  nr = params.nr
  @inbounds for p1 = 1:kc:k
    p2 = min(k, p1 + kc - 1)
    A_p = @view A[:, p1:p2]
    B_p = @view B[:, p1:p2]
    real_kc = (p2 - p1 + 1)
    kc_nr = real_kc * nr
    kc_mr = real_kc * mr
    packB!(B_packed, B_p, params)
    loop3!(C, A_p, A_packed, B_packed, kc_nr, kc_mr, params)
  end
end

@inline function loop3!(C, A, A_packed, B_packed, kc_nr, kc_mr, params::TilingParams)
  m = size(C, 1)
  mc = params.mc
  @inbounds for i1 = 1:mc:m
    i2 = min(m, i1 + mc - 1)
    C_i = @view C[i1:i2, :]
    A_i = @view A[i1:i2, :]
    packA!(A_packed, A_i, params)
    loop2!(C_i, A_packed, B_packed, kc_nr, kc_mr, params)
  end
end

@inline function loop2!(C, A_packed, B_packed, kc_nr, kc_mr, params::TilingParams)
  n = size(C, 2)
  nr = params.nr
  @inbounds for (t, j1) = enumerate(1:nr:n)
    j2 = min(n, j1 + nr - 1)
    C_j = @view C[:, j1:j2]
    tile_rng = (1 + kc_nr * (t - 1)) : (t * kc_nr)
    B_packed_j = @view B_packed[tile_rng]
    loop1!(C_j, A_packed, B_packed_j, kc_mr, params)
  end
end

@inline function loop1!(C, A_packed, B_packed, kc_mr, params::TilingParams)
  m = size(C, 1)
  mr = params.mr
  @inbounds for (t, i1) = enumerate(1:mr:m)
    i2 = min(m, i1 + mr - 1)
    C_i = @view C[i1:i2, :]
    tile_rng = (1 + kc_mr * (t - 1)) : (t * kc_mr)
    A_packed_i = @view A_packed[tile_rng]
    microknl!(C_i, A_packed_i, B_packed, params)
  end
end

@inline function microknl!(C, A, B, params::TilingParams)
  # Get relevant tile params
  m = size(parent(C), 1)
  mr = params.mr
  nr = params.nr
  vl = params.vl
  mr_vl = div(mr, vl)
  pend = div(size(B, 1), nr)

  # Set relevant types
  FT = eltype(C)
  VT = Vec{vl, FT}

  # Initialize vector storage
  C_tile = MArray{Tuple{mr_vl, nr}, VT}(undef)
  A_panel = MVector{mr_vl, VT}(undef)

  # Load the columns of the microtile of C
  @inbounds @unroll for j = 1:nr
    @unroll for i = 1:mr_vl
      offset = ((i - 1) * vl + (j - 1) * m) * sizeof(FT)
      C_tile[i, j] = vload(VT, pointer(C) + offset, nothing, Val(true))
    end
  end

  # Do the rank-one updates
  @inbounds @unroll for p = 1:pend
    # Load pieces of `A_panel`
    @unroll for i = 1:mr_vl
      offset = ((p - 1) * mr + (i - 1) * vl) * sizeof(FT)
      A_panel[i] = vload(VT, pointer(A) + offset, nothing, Val(true))
    end

    # Do outer product updates
    @unroll for j = 1:nr
      β = B[j + nr * (p - 1)]
      @unroll for i = 1:mr_vl
        C_tile[i, j] = muladd(A_panel[i], β, C_tile[i, j])
      end
    end
  end

  # Write back from vector storage
  @inbounds @unroll for j = 1:nr
    @unroll for i = 1:mr_vl
      offset = ((i - 1) * vl + (j - 1) * m) * sizeof(FT)
      vstore(C_tile[i, j], pointer(C) + offset, nothing, Val(true), Val(true))
    end
  end
end

end # module

max_stars_repo_name: bischtob/MyGEMM

max_issues_repo_name: bischtob/MyGEMM

max_forks_repo_name: bischtob/MyGEMM


text:
isdefined(Base, :__precompile__) && __precompile__()

module IntArrays

export IntArray, IntVector, IntMatrix,
    radixsort,
    radixsort!

import Base:
    convert,
    call,
    linearindexing,
    getindex,
    setindex!,
    size,
    length,
    sizeof,
    similar,
    fill!,
    copy!,
    resize!,
    push!,
    pop!,
    append!,
    reverse!

include("buffer.jl")
include("array.jl")
include("vector.jl")
include("matrix.jl")

end # module


text:
using Jchemo, LinearAlgebra

n = 50 ; p = 7
Xtrain = rand(n, p)
m = 3
Xtest = rand(m, p)

nlv = 5
fm = pcasvd(Xtrain; nlv = nlv) ;
#fm = pcaeigen(Xtrain; nlv = nlv) ;
pnames(fm)
fm.T
fm.T' * fm.T
fm.P' * fm.P

Jchemo.transform(fm, Xtest)

res = Jchemo.summary(fm, Xtrain) ;
pnames(res)
res.explvar
res.contr_var
res.coord_var
res.cor_circle

# Weighted PCA
w = collect(1:n)
nlv = 5
fm = pcasvd(Xtrain, w; nlv = nlv) ;
D = Diagonal(mweights(w)) ;
fm.T' * D * fm.T
fm.P' * fm.P


text:
using BinaryBuilder, Pkg

name = "METIS"
version = v"4.0.3"

# Collection of sources required to build METIS
sources = [
    ArchiveSource("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/OLD/metis-$(version).tar.gz",
                  "5efa35de80703c1b2c4d0de080fafbcf4e0d363a21149a1ad2f96e0144841a55"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == *-mingw* ]]; then
  COPTIONS="${COPTIONS} -D__VC__"  # to resolve missing srand48/drand48 symbols
fi

# build libmetis.a
cd $WORKSPACE/srcdir/metis-*
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
  atomic_patch -p1 ${f}
done
cd Lib
make -j${nproc} COPTIONS="${COPTIONS} -fPIC"
cd ..

# We copy the .a files into ${prefix}/lib since the main purpose is to link them in other builds.
# Specifically this is in a separate location than the typical location for libraries on Windows.
mkdir -p ${prefix}/lib
mv libmetis.a ${prefix}/lib
mkdir -p ${prefix}/include
cp Lib/metis.h ${prefix}/include
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    FileProduct("lib/libmetis.a", :libmetis_a)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"6")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
export ImmutablePolynomial

"""
    ImmutablePolynomial{T<:Number}(coeffs::AbstractVector{T}, var=:x)

Construct an immutable (static) polynomial from its coefficients `a`,
lowest order first, optionally in terms of the given variable `x`
where `x` can be a character, symbol, or string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct this through
`ImmutablePolynomial((a_0, a_1, ..., a_n))`.

The usual arithmetic operators are overloaded to work with polynomials
as well as with combinations of polynomials and scalars. However,
operations involving two polynomials of different variables causes an
error, though for `+` and `*` operations, constant polynomials are
treated as having no variable. (This adds a runtime check, but is useful when working with matrices of polynomials.)

As the coefficient size is a compile-time constant, immutable polynomials can take advantage of faster polynomial evaluation provided by `evalpoly` from Julia 1.4.

    # Examples

```jldoctest
julia> using  Polynomials

julia> ImmutablePolynomial((1, 0, 3, 4))
ImmutablePolynomial(1 + 3*x^2 + 4*x^3)

julia> ImmutablePolynomial((1, 2, 3), :s)
ImmutablePolynomial(1 + 2*s + 3*s^2)

julia> one(ImmutablePolynomial)
ImmutablePolynomial(1.0)
```

!!! note
    This was modeled after https://github.com/tkoolen/StaticUnivariatePolynomials.jl by @tkoolen.

"""
struct ImmutablePolynomial{N, T <: Number} <: StandardBasisPolynomial{T}
    coeffs::NTuple{N, T}
    var::Symbol
    function ImmutablePolynomial{N,T}(coeffs::NTuple{N,T}, var::Symbol=:x) where {N, T <: Number}
        new{N,T}(coeffs, var)
    end
    function ImmutablePolynomial{N,T}(coeffs::NTuple{M,S}, var::Symbol=:x) where {N, M, T, S<: Number}
        M > N && throw(ArgumentError("coeffs too big for the specified type"))
        if M < N
            cs = NTuple{N,T}(i <= M ? T(coeffs[i]) : zero(T) for i in 1:N)
        else
            cs = NTuple{N,T}(T(c) for c in coeffs)
        end
        new{N,T}(cs, var)
    end
    function ImmutablePolynomial{N,T}(coeffs::AbstractVector{S}, var::Symbol=:x) where {N, T <: Number, S}
        M = findlast(!iszero, coeffs)
        M == nothing && return zero(ImmutablePolynomial{N,T})
        if M < N
            cs = NTuple{N,T}(i <= M ? T(coeffs[i]) : zero(T) for i in 1:N)
        else
            cs = NTuple{N,T}(T(c) for c in coeffs)
        end
        new{N,T}(cs, var)
    end
end

@register1 ImmutablePolynomial

## promote N,M case
function Base.promote(p::ImmutablePolynomial{N,T}, q::ImmutablePolynomial{M,S}) where {N,T,M,S}
    NN, R = max(N,M), promote_type(T,S)
    ImmutablePolynomial{NN,R}(p.coeffs, p.var), ImmutablePolynomial{NN,R}(q.coeffs, q.var)
end
## Need defaults when N not specified
# tuple
ImmutablePolynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T} =
     ImmutablePolynomial{N,T}(coeffs, Symbol(var))

# vector
function ImmutablePolynomial(coeffs::AbstractVector{T}, var::SymbolLike = :x) where{N,T}
    M = length(coeffs)
    M == 0 && return zero(ImmutablePolynomial{1, T}, var)
    ImmutablePolynomial{M,T}(NTuple{M,T}(x for x in coeffs), Symbol(var))
end

# number
ImmutablePolynomial(n::T, var::SymbolLike = :x) where {T <: Number} =
     ImmutablePolynomial{1,T}(NTuple{1,T}(n), Symbol(var))

# variable
ImmutablePolynomial(var::SymbolLike = :x)  = variable(ImmutablePolynomial{2, Int}, var)

# Convenience; pass tuple to Polynomial
# Not documented, not sure this is a good idea as P(...)::P is not true...
Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T} =
    ImmutablePolynomial(coeffs, var)
function Polynomial{T}(coeffs::NTuple{N,S}, var::SymbolLike = :x) where{N,T,S}
    ImmutablePolynomial{N,T}(T.(coeffs), var)
end

##
## ----
##

# overrides from common.jl due to coeffs possibly being padded, coeffs being non mutable, ...
Base.promote_rule(::Type{<:ImmutablePolynomial{N, T}}, ::Type{<:ImmutablePolynomial{M,S}}) where {N,T,M,S} =
    ImmutablePolynomial{max(N,M), promote_type(T, S)}

Base.convert(::Type{<:ImmutablePolynomial{N, T}}, p::ImmutablePolynomial{M,S}) where {N,T,M,S} =
    ImmutablePolynomial{N,T}(p.coeffs,  p.var)

Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p), p.var)

function Base.hash(p::ImmutablePolynomial{N,T}, h::UInt) where {N,T}
    n = findlast(!iszero, coeffs(p))
    n == nothing && return hash(p.var, hash(NTuple{0,T}(),h))
    hash(p.var, hash(coeffs(p)[1:n], h))
end

Base.one(::Type{ImmutablePolynomial{N,T}},
         var=:x) where {N,T}  = ImmutablePolynomial(tuple((i==1 ? one(T) : zero(T) for i in  1:N)...), var)
Base.one(::Type{ImmutablePolynomial{N}}, var=:x) where {N}  = one(ImmutablePolynomial{N,Int}, var)
Base.zero(::Type{ImmutablePolynomial{N,T}}, var=:x) where {N,T} = ImmutablePolynomial(zeros(T, N), var)
Base.zero(::Type{ImmutablePolynomial{N}}, var=:x) where {N}  = zero(ImmutablePolynomial{N,Int}, var)

function degree(p::ImmutablePolynomial{N,T}) where {N, T}
    n = findlast(!iszero, coeffs(p))
    n == nothing ? -1 : n-1
end
isconstant(p::ImmutablePolynomial{0}) = true
isconstant(p::ImmutablePolynomial{1}) = true
function  isconstant(p::ImmutablePolynomial{N})  where {N}
    for i in 2:length(p.coeffs)
        !iszero(p.coeffs[i]) && return false
    end
    return true
end
for op in [:isequal, :(==)]
    @eval function Base.$op(p1::ImmutablePolynomial{N,T}, p2::ImmutablePolynomial{M,S}) where {N,T,M,S}
        (p1.var == p2.var) || return false
        p1s, p2s = coeffs(p1), coeffs(p2)
        (N == M  && $op(p1s,p2s)) &&  return  true
        n1 = findlast(!iszero, p1s) # now trim out zeros
        n2 = findlast(!iszero, p2s)
        (n1 == nothing && n2 == nothing) && return true
        (n1 == nothing || n2  == nothing) && return false
        $op(p1s[1:n1],p2s[1:n2]) &&  return true
        false
    end
end

function Base.isapprox(p1::ImmutablePolynomial{N,T},
                       p2::ImmutablePolynomial{M,S};
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {N,T,M,S}
    p1.var == p2.var || error("p1 and p2 must have same var")
    NN = max(N,M)
    for i in 1:NN-1
        isapprox(p1[i],p2[i], rtol=rtol, atol=atol) || return false
    end
    true
end

function Base.isapprox(p1::ImmutablePolynomial{N,T}, n::S;
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {N, T,S <: Number}
    return isapprox(p1, ImmutablePolynomial(n, p1.var))
end

function Base.chop(p::ImmutablePolynomial{N,T};
              rtol::Real = Base.rtoldefault(real(T)),
              atol::Real = 0)  where {N, T}
    cs = coeffs(p)
    for i in N:-1:1
        if !isapprox(cs[i], zero(T), rtol=rtol, atol=atol)
            return ImmutablePolynomial{i,T}(cs[1:i], p.var)
        end
    end
    zero(ImmutablePolynomial{0,T})
end

function Base.truncate(p::ImmutablePolynomial{N,T};
                       rtol::Real = Base.rtoldefault(real(T)),
                       atol::Real = 0)  where {N, T}
    q = chop(p, rtol=rtol, atol=atol)
    iszero(q) && return  q
    cs = coeffs(q)
    thresh = maximum(abs,cs) * rtol + atol
    ImmutablePolynomial{length(cs), T}(map(c->abs(c) <= thresh ? zero(T) : c, coeffs(q)), p.var)
end

# no in-place chop! and truncate!
chop!(p::ImmutablePolynomial; kwargs...) =  chop(p; kwargs...)
truncate!(p::ImmutablePolynomial; kwargs...) =  truncate(p; kwargs...)

LinearAlgebra.conj(p::P) where {P <: ImmutablePolynomial} = P(conj([aᵢ for aᵢ in coeffs(p)]))

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

(p::ImmutablePolynomial{N, T})(x::S) where {N, T,S} = evalpoly(x, coeffs(p))

function Base.:+(p1::ImmutablePolynomial{N,T}, p2::ImmutablePolynomial{M,S}) where {N,T,M,S}

    isconstant(p1) && return p2 + p1[0]
    isconstant(p2) && return p1 + p2[0]
    p1.var != p2.var && error("Polynomials must have same variable")

    R = Base.promote_op(+, T,S)
    NN = max(N, M)
    P = ImmutablePolynomial{NN,R}
    ⊕(convert(P,p1), convert(P,p2))
end

function ⊕(p1::ImmutablePolynomial{N,T}, p2::ImmutablePolynomial{N,T}) where {N,T}
    cs = NTuple{N, T}(p1[i] + p2[i] for i in 0:N-1)
    return ImmutablePolynomial{N, T}(cs, p1.var)
end

function Base.:+(p::ImmutablePolynomial{N, T}, c::S) where {N, T, S<:Number}
    R = promote_type(T,S)
    as = NTuple{N,R}(i == 1 ? ai + c : ai for (i,ai) in enumerate(p.coeffs))
    ImmutablePolynomial(as, p.var)
end

function Base.:*(p1::ImmutablePolynomial{N,T}, p2::ImmutablePolynomial{M,S}) where {N,T,M,S}
    isconstant(p1) && return p2 * p1[0]
    isconstant(p2) && return p1 * p2[0]
    p1.var != p2.var && error("Polynomials must have same variable")
    p1 ⊗ p2
end

## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl
@generated function ⊗(p1::ImmutablePolynomial{N,T}, p2::ImmutablePolynomial{M,S}) where {N,T,M,S}
    R = eltype(one(T)*one(S))
    P = M + N - 1
    exprs = Any[nothing for i = 1 : P]
    for i in 1 : N
        for j in 1 : M
            k = i + j - 1
            if exprs[k] === nothing
                exprs[k] = :(p1.coeffs[$i] * p2.coeffs[$j])
            else
                exprs[k] = :(muladd(p1.coeffs[$i], p2.coeffs[$j], $(exprs[k])))
            end
        end
    end
    return quote
        Base.@_inline_meta
        ImmutablePolynomial{$(max(0,P)), $(R)}(tuple($(exprs...)), p1.var)
    end
end

function Base.:*(p::ImmutablePolynomial{N,T}, c::S) where {N,T,S <: Number}
    R = eltype(one(T)*one(S))
    cs = NTuple{N,R}(p[i]*c for i in eachindex(p))
    return ImmutablePolynomial{N,R}(cs, p.var)
end

Base.:-(p::ImmutablePolynomial{N,T}) where {N,T} = ImmutablePolynomial(NTuple{N,T}(-pi for pi in p.coeffs), p.var)

Base.to_power_type(p::ImmutablePolynomial{N,T}) where {N,T} = p


text:
using Documenter
using SDtoolbox

makedocs(
    sitename = "SDtoolbox",
    format = Documenter.HTML(prettyurls = false),
    modules = [SDtoolbox],
    pages = [
        "Home" => "index.md",
        "Postshock" => "postshock.md",
        "ZND" => "znd.md",
        "CV" => "cv.md",
        "Cell size calculation" => "cell_size.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/daharn/SDtoolbox.jl"
)


text:
export HPoint, getz, dist, midpoint, RandomHPoint, between, polar

"""
`HPoint(z::Complex)` creates a new point in the hyperbolic plane.
The argument `z` must have absolute value less than 1.

`HPoint(r,theta)` creates a new point with polar coordinates
`(r,theta)`. See also: `polar`.
"""
struct HPoint <: HObject
    z::Complex{Float64}
    attr::Dict{Symbol,Any}
    function HPoint(z::Complex)
        if _mag(z) >=  1
            throw(DomainError(z, "absolute value is too large"))
        end
        P = new(z,Dict{Symbol,Any}())
        set_color(P)
        set_radius(P)
        return P
    end
end

HPoint(P::HPoint) = HPoint(getz(P))  # copy constructor
HPoint(z::Number) = HPoint(Complex(z))
HPoint() = HPoint(0)

function show(io::IO,P::HPoint)
    r,theta = polar(P)
    print(io,"HPoint($r, $theta)")
end

"""
`RandomHPoint()` generates a point at random in the hyperbolic plane.
"""
RandomHPoint() = HPoint(randn(), 2*pi*rand())

"""
`getz(P::HPoint)` returns the point (complex number) in the
interior of the unit disc that represents `P`.
"""
getz(P::HPoint) = P.z

"""
`dist(P,Q)` gives the distance betwen two points in the
hyperbolic plane. If `Q` is omitted, give the distance
from `P` to `HPoint(0)`.
"""
function dist(P::HPoint, Q::HPoint)
    a = getz(P)
    b = getz(Q)
    delta = 2 * _mag(a-b)/(1-_mag(a))/(1-_mag(b))
    return acosh(1+delta)
end

dist(P::HPoint) = dist(P, HPoint(0))

HPoint(r::Real, theta::Real) = HPoint( solve_dist(r) * exp(im*theta) )

"""
`polar(P::HPoint)` gives the polar coordinates of `P`
"""
function polar(P::HPoint)
    r = dist(P)
    theta = angle(getz(P))
    return (r,theta)
end

(==)(P::HPoint,Q::HPoint) = _mag(getz(P)-getz(Q)) <= THRESHOLD*eps(1.0)

function _dist(t::Real)
    delta = 2*t*t/(1-t*t)
    return acosh(1+delta)
end

"""
`solve_dist(d)` is the inverse of `_dist()`
"""
function solve_dist(d::Real)
    ex = exp(d)
    return (ex-1)/(ex+1)
end

"""
`midpoint(p,q)` finds the mid point of the line segment
from `p` to `q`. Also `midpoint(L::HSegment)`.
"""
function midpoint(p::HPoint, q::HPoint)::HPoint
    if p==q
        return p
    end
    d = dist(p,q)
    f = move2xplus(p,q)

    t = solve_dist(d/2)

    r = HPoint(t)
    h = inv(f)
    return h(r)
end

"""
`between(a,b,c)` determines if the hyperbolic point `b`
lies on the segment from `a` to `c`.
"""
function between(a::HPoint, b::HPoint, c::HPoint)::Bool
    if a==b || b==c
        return true
    end
    f = move2xplus(a,c)
    aa = getz(f(a))
    bb = getz(f(b))
    cc = getz(f(c))

    if abs(imag(bb)) > THRESHOLD * eps(1.0)
        return false
    end

    if real(aa)-THRESHOLD*eps(1.0) <= real(bb) <= real(cc)+THRESHOLD*eps(1.0)
        return true
    end
    return false
end


text:
module HTMLWriterTests

using Compat.Test
using Compat

import Documenter.Writers.HTMLWriter: jsescape, generate_version_file

@testset "HTMLWriter" begin
    @test jsescape("abc123") == "abc123"
    @test jsescape("▶αβγ") == "▶αβγ"
    @test jsescape("") == ""

    @test jsescape("a\nb") == "a\\nb"
    @test jsescape("\r\n") == "\\r\\n"
    @test jsescape("\\") == "\\\\"

    @test jsescape("\"'") == "\\\"\\'"

    # Ref: #639
    @test jsescape("\u2028") == "\\u2028"
    @test jsescape("\u2029") == "\\u2029"
    @test jsescape("policy to
 delete.") == "policy to\\u2028 delete."

    mktempdir() do tmpdir
        versions = ["stable", "latest", "release-0.2", "release-0.1", "v0.2.6", "v0.1.1", "v0.1.0"]
        cd(tmpdir) do
            mkdir("foobar")
            for version in versions
                mkdir(version)
            end
        end

        generate_version_file(tmpdir)

        versions_file = joinpath(tmpdir, "versions.js")
        @test isfile(versions_file)
        contents = String(read(versions_file))
        @test !occursin("foobar", contents) # only specific directories end up in the versions file
        # let's make sure they're in the right order -- they should be sorted in the output file
        last = 0:0
        for version in versions
            this = Compat.findfirst(version, contents)
            @test this !== nothing
            @test first(last) < first(this)
            last = this
        end
    end
end

end