text:
using RipQP, LLSModels, QPSReader, QuadraticModels, Quadmath
using DelimitedFiles, LinearAlgebra, MatrixMarket, Test

function createQuadraticModel128(qpdata; name = "qp_pb")
  return QuadraticModel(
    convert(Array{Float128}, qpdata.c),
    qpdata.qrows,
    qpdata.qcols,
    convert(Array{Float128}, qpdata.qvals),
    Arows = qpdata.arows,
    Acols = qpdata.acols,
    Avals = convert(Array{Float128}, qpdata.avals),
    lcon = convert(Array{Float128}, qpdata.lcon),
    ucon = convert(Array{Float128}, qpdata.ucon),
    lvar = convert(Array{Float128}, qpdata.lvar),
    uvar = convert(Array{Float128}, qpdata.uvar),
    c0 = Float128(qpdata.c0),
    x0 = zeros(Float128, length(qpdata.c)),
    name = name,
  )
end
qps1 = readqps("QAFIRO.SIF") #lower bounds
qps2 = readqps("HS21.SIF") # low/upp bounds
qps3 = readqps("HS52.SIF") # free bounds
qps4 = readqps("AFIRO.SIF") # LP
qps5 = readqps("HS35MOD.SIF") # fixed variables

@testset "mono_mode" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(presolve = false))
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2))
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats1.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3))
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "multi_mode" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "dynamic_regularization" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic), history = true),
    display = false,
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic)),
    display = false,
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    iconf = InputConfig(sp = K2LDLParams(regul = :dynamic)),
    display = false,
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "centrality_corrections" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(kc = -1), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(kc = 2), display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(kc = 2), display = false)
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "Float128" begin
  qm128_1 = createQuadraticModel128(qps1)
  stats1 = ripqp(
    qm128_1,
    itol = InputTol(ϵ_rb32 = 0.1, ϵ_rb64 = 0.01),
    iconf = InputConfig(mode = :multi, normalize_rtol = false),
    display = false,
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  qm128_2 = createQuadraticModel128(qps2)
  stats2 = ripqp(qm128_2, display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable
end

@testset "refinement" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement = :zoom), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(refinement = :ref), display = false) # automatic centrality corrections computation
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    iconf = InputConfig(mode = :multi, scaling = false, refinement = :multizoom),
    display = false,
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    iconf = InputConfig(mode = :multi, refinement = :multiref),
    display = false,
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "K2_5" begin
  stats1 = ripqp(QuadraticModel(qps1), display = false, iconf = InputConfig(sp = K2_5LDLParams()))
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(sp = K2_5LDLParams(), mode = :multi),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(
    QuadraticModel(qps3),
    display = false,
    iconf = InputConfig(sp = K2_5LDLParams(regul = :dynamic)),
  )
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "IPF" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(solve_method = :IPF, mode = :multi),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), display = false, iconf = InputConfig(solve_method = :IPF))
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(solve_method = :IPF, sp = K2_5LDLParams(), refinement = :zoom),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(solve_method = :IPF, mode = :multi, refinement = :multiref),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(
      solve_method = :IPF,
      mode = :multi,
      refinement = :multizoom,
      sp = K2_5LDLParams(),
    ),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats2 = ripqp(
    QuadraticModel(qps2),
    display = false,
    iconf = InputConfig(
      solve_method = :IPF,
      refinement = :zoom,
      sp = K2LDLParams(regul = :dynamic),
    ),
  )
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable
end

@testset "LLS" begin
  # least-square problem without constraints
  m, n = 20, 10
  A = rand(m, n)
  b = rand(m)
  nls = LLSModel(A, b)
  statsnls = ripqp(nls, itol = InputTol(ϵ_rb = sqrt(eps()), ϵ_rc = sqrt(eps())), display = false)
  x, r = statsnls.solution[1:n], statsnls.solution[(n + 1):end]
  @test norm(x - A \ b) ≤ norm(b) * sqrt(eps())
  @test norm(A * x - b - r) ≤ norm(b) * sqrt(eps())

  # least-square problem with constraints
  A = rand(m, n)
  b = rand(m)
  lcon, ucon = zeros(m), fill(Inf, m)
  C = ones(m, n)
  lvar, uvar = fill(-10.0, n), fill(200.0, n)
  nls = LLSModel(A, b, lvar = lvar, uvar = uvar, C = C, lcon = lcon, ucon = ucon)
  statsnls = ripqp(nls, display = false)
  x, r = statsnls.solution[1:n], statsnls.solution[(n + 1):end]
  @test length(r) == m
  @test norm(A * x - b - r) ≤ norm(b) * sqrt(eps())
end

@testset "matrixwrite" begin
  stats1 = ripqp(
    QuadraticModel(qps1),
    display = false,
    iconf = InputConfig(w = SystemWrite(write = true, name = "test_", kfirst = 4, kgap = 1000)),
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  K = MatrixMarket.mmread("test_K_iter4.mtx")
  rhs_aff = readdlm("test_rhs_iter4_aff.rhs", Float64)[:]
  rhs_cc = readdlm("test_rhs_iter4_cc.rhs", Float64)[:]
  @test size(K, 1) == size(K, 2) == length(rhs_aff) == length(rhs_cc)
end

@testset "KrylovK2" begin
  for kmethod in [:minres, :minres_qlp]
    stats2 = ripqp(
      QuadraticModel(qps2),
      display = true,
      iconf = InputConfig(
        sp = K2KrylovParams(kmethod = kmethod, preconditioner = :Identity),
        history = true,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = true,
      iconf = InputConfig(sp = K2KrylovParams(uplo = :U, kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "KrylovK2_5" begin
  for kmethod in [:minres, :minres_qlp]
    stats1 = ripqp(
      QuadraticModel(qps1),
      display = false,
      iconf = InputConfig(
        sp = K2_5KrylovParams(kmethod = kmethod, preconditioner = :Identity),
        history = true,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats1.objective, -1.59078179, atol = 1e-1)
    @test stats1.status == :acceptable

    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(sp = K2_5KrylovParams(uplo = :U, kmethod = kmethod), solve_method = :IPF),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K2_5KrylovParams(kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "Krylov K3" begin
  for kmethod in [:bicgstab, :usymqr, :qmr]
    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(
        sp = K3KrylovParams(kmethod = kmethod),
        solve_method = :IPF,
        history = true,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K3KrylovParams(uplo = :U, kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "KrylovK3_5" begin
  for kmethod in [:minres, :minres_qlp]
    stats1 = ripqp(
      QuadraticModel(qps1),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(kmethod = kmethod), history = true),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats1.objective, -1.59078179, atol = 1e-1)
    @test stats1.status == :acceptable

    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(uplo = :U, kmethod = kmethod), solve_method = :IPF),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "Krylov K1" begin
  stats4 = ripqp(
    QuadraticModel(qps4),
    display = false,
    iconf = InputConfig(sp = K1KrylovParams(), solve_method = :IPF, history = true),
    itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
  )
  @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
  @test stats4.status == :acceptable

  stats4 = ripqp(
    QuadraticModel(qps4),
    display = false,
    iconf = InputConfig(sp = K1KrylovParams(uplo = :U), solve_method = :PC, history = true),
    itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
  )
  @test isapprox(stats4.objective, -4.6475314286e02, atol = 1e-2)
  @test stats4.status == :acceptable
end

function QuadraticModelMaximize(qps, x0 = zeros(qps.nvar))
  QuadraticModel(
    -qps.c,
    qps.qrows,
    qps.qcols,
    -qps.qvals,
    Arows = qps.arows,
    Acols = qps.acols,
    Avals = qps.avals,
    lcon = qps.lcon,
    ucon = qps.ucon,
    lvar = qps.lvar,
    uvar = qps.uvar,
    c0 = -qps.c0,
    x0 = x0,
    minimize = false,
  )
end

@testset "maximize" begin
  stats1 = ripqp(QuadraticModelMaximize(qps1), display = false)
  @test isapprox(stats1.objective, 1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModelMaximize(qps2), display = false)
  @test isapprox(stats2.objective, 9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModelMaximize(qps3), display = false)
  @test isapprox(stats3.objective, -5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "presolve" begin
  stats5 = ripqp(QuadraticModel(qps5), display = false)
  @test isapprox(stats5.objective, 0.250000001, atol = 1e-2)
  @test stats5.status == :acceptable

  stats5 = ripqp(QuadraticModel(qps5), iconf = InputConfig(sp = K2KrylovParams()), display = false)
  @test isapprox(stats5.objective, 0.250000001, atol = 1e-2)
  @test stats5.status == :acceptable
end


text:
export package_specs

"""
Return a `Vector` of `PackageSpec`s that are needed by the
workspace that's the target of the `URI`."
"""
function package_specs end

function package_specs(string::String)
    add_packages_from(URI(string))
end

function package_specs(uri::URI)
    package_specs = Vector{Pkg.Types.PackageSpec}()
    project = fetch_toml(uri_add_path(uri, "Project.toml"))
    manifest = fetch_toml(uri_add_path(uri, "Manifest.toml"))
    function mfind(name, uuid)
        for (pkg, recs) in manifest
            if name != pkg
                continue
            end
            for rec in recs
                if uuid == rec["uuid"]
                    return rec
                end
            end
        end
        return nothing
    end
    for (name, uuid) in project["deps"]
        rec = mfind(name, uuid)
        if rec == nothing
            throw(ErrorException("Package $name not found in manifest"))
        end
        psargs = Dict(:name => name,
                      :uuid => uuid)
        if haskey(rec, "repo-url")
            psargs[:url] = rec["repo-url"]
        end
        if haskey(rec, "version")
            psargs[:version] = rec["version"]
        end

        push!(package_specs,
              PackageSpec(; psargs...))
    end
    package_specs
end

ext: jl

max_stars_repo_name: MarkNahabedian/WebBasedWorkspace

max_issues_repo_name: MarkNahabedian/WebBasedWorkspace

max_forks_repo_name: MarkNahabedian/WebBasedWorkspace


text:
using DataStructures
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
import Base.cat

View = Plasm.view

# Graphic text font implementation
# Vector definition of printable ASCII codes as one-dimensional LAR models.
# Font design from *Geometric Programming for Computer-Aided Design*, Wiley, 2003.

const ascii32 = ([0.0; 0.0],Array{Int64,1}[[1]])
const ascii33 = ([1.75 1.75 2.0 2.0 1.5 1.5; 1.75 5.5 0.5 1.0 1.0 0.5],Array{Int64,1}[[1,2],[3,4],[4,5],[5,6],[6,3]])
const ascii34 = ([1.0 2.0 2.0 2.0 1.5 1.5 2.0 3.0 3.0 3.0 2.5 2.5; 4.0 5.0 5.5 6.0 6.0 5.5 4.0 5.0 5.5 6.0 6.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3],[7,8],[8,9],[9,10],[10,11],[11,12],[12,9]])
const ascii35 = ([1.0 3.0 1.0 3.0 1.25 1.75 2.25 2.75; 2.5 2.5 3.5 3.5 1.75 4.0 1.75 4.0],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8]])
const ascii36 = ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 2.0 2.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 -0.5 6.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[13,14]])
const ascii37 = ([2.5 2.5 2.0 2.0 2.5 2.5 2.0 2.0 1.5 3.5; 0.0 0.5 0.5 0.0 5.5 6.0 6.0 5.5 5.5 11.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,5],[9,10]])
const ascii38 = ([4.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 5.0 4.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14]])
const ascii39 = ([1.0 2.0 2.0 2.0 1.5 1.5; 4.0 5.0 5.5 6.0 6.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3]])
const ascii40 = ([2.0 1.0 0.5 1.0 2.0; 0.0 1.0 3.0 5.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
const ascii41 = ([2.0 3.0 3.5 3.0 2.0; 0.0 1.0 3.0 5.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
const ascii42 = ([1.0 3.0 2.0 2.0 1.0 3.0 1.0 3.0; 3.0 3.0 2.0 4.0 2.0 4.0 4.0 2.0],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8]])
const ascii43 = ([1.0 3.0 2.0 2.0; 3.0 3.0 2.0 4.0],Array{Int64,1}[[1,2],[3,4]])
const ascii44 = ([1.0 2.0 2.0 2.0 1.5 1.5; -1.0 0.0 0.5 1.0 1.0 0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3]])
const ascii45 = ([1.0 3.0; 3.0 3.0],Array{Int64,1}[[1,2]])
const ascii46 = ([2.0 2.0 1.5 1.5 2.0; 0.0 0.5 0.5 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
const ascii47 = ([1.0 3.0; 0.0 6.0],Array{Int64,1}[[1,2]])
const ascii48 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[4,8]])
const ascii50 = ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 4.0 5.0 6.0 6.0 5.0 4.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
const ascii51 = ([0.0 4.0 2.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 2.0 1.0 0.0 0.0 1.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
const ascii52 = ([4.0 0.0 4.0 4.0; 1.0 1.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii53 = ([4.0 0.0 0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10]])
const ascii54 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 6.0 6.0 5.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11]])
const ascii55 = ([0.0 0.0 4.0 0.0; 5.0 6.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii56 = ([1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 4.0 4.0 3.0 1.0 0.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 4.0 5.0 6.0 6.0 5.0 4.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[6,5],[5,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,6]])
const ascii57 = ([0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 5.0 6.0 6.0 5.0 3.0 2.0 2.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11]])
const ascii58 = ([2.0 2.0 1.5 1.5 2.0 2.0 1.5 1.5; 1.0 1.5 1.5 1.0 3.0 3.5 3.5 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,5]])
const ascii59 = ([2.0 2.0 1.5 1.5 1.0 2.0 2.0 2.0 1.5 1.5; 3.0 3.5 3.5 3.0 -0.5 0.5 1.0 1.5 1.5 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,9],[9,10],[10,7]])
const ascii60 = ([3.0 0.0 3.0; 6.0 3.0 0.0],Array{Int64,1}[[1,2],[2,3]])
const ascii61 = ([1.0 3.0 1.0 3.0; 2.5 2.5 3.5 3.5],Array{Int64,1}[[1,2],[3,4]])
const ascii62 = ([1.0 4.0 1.0; 6.0 3.0 0.0],Array{Int64,1}[[1,2],[2,3]])
const ascii63 = ([2.0 2.0 1.5 1.5 1.75 1.75 3.0 3.0 2.0 1.0 0.0 0.0; 0.0 0.5 0.5 0.0 1.0 2.75 4.0 5.0 6.0 6.0 5.0 4.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12]])
const ascii64 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 2.0 1.0 2.0 3.0 2.0; 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 1.0 2.0 3.0 2.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13]])
const ascii65 = ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 0.0 5.0 6.0 6.0 5.0 0.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[7,8]])
const ascii66 = ([0.0 0.0 3.0 4.0 4.0 3.0 4.0 4.0 3.0 0.0 0.0 3.0; 0.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[11,12]])
const ascii67 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
const ascii68 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 1.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii69 = ([4.0 0.0 0.0 4.0 0.0 3.0; 0.0 0.0 6.0 6.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[5,6]])
const ascii70 = ([0.0 0.0 4.0 0.0 3.0; 0.0 6.0 6.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii71 = ([2.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10]])
const ascii72 = ([0.0 0.0 4.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
const ascii73 = ([2.0 2.0 1.0 3.0 1.0 3.0; 0.0 6.0 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
const ascii74 = ([0.0 1.0 2.0 3.0 3.0 2.0 4.0; 1.0 0.0 0.0 1.0 6.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[6,7]])
const ascii75 = ([4.0 0.0 4.0 0.0 0.0; 6.0 3.0 0.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii76 = ([4.0 0.0 0.0; 0.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3]])
const ascii77 = ([0.0 0.0 2.0 4.0 4.0; 0.0 6.0 4.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
const ascii78 = ([0.0 0.0 4.0 4.0 4.0; 0.0 6.0 2.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii79 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
const ascii80 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii81 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0 1.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[10,11]])
const ascii82 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0 3.0 4.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 2.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[8,9]])
const ascii83 = ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12]])
const ascii84 = ([2.0 2.0 0.0 4.0; 0.0 6.0 6.0 6.0],Array{Int64,1}[[1,2],[3,4]])
const ascii85 = ([0.0 0.0 1.0 3.0 4.0 4.0; 6.0 1.0 0.0 0.0 1.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6]])
const ascii86 = ([0.0 2.0 4.0; 6.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3]])
const ascii87 = ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 6.0 3.0 0.0 3.0 0.0 3.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii88 = ([0.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0],Array{Int64,1}[[1,2],[3,4]])
const ascii89 = ([0.0 2.0 4.0 2.0 2.0; 6.0 2.0 6.0 2.0 0.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii90 = ([0.0 4.0 0.0 4.0; 6.0 6.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii91 = ([2.0 1.0 1.0 2.0; 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii92 = ([1.0 3.0; 6.0 0.0],Array{Int64,1}[[1,2]])
const ascii93 = ([2.0 3.0 3.0 2.0; 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii94 = ([1.0 2.0 3.0; 5.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3]])
const ascii95 = ([1.0 4.0; 0.0 0.0],Array{Int64,1}[[1,2]])
const ascii96 = ([2.0 2.0 3.0 2.5 2.5 2.0; 4.5 5.0 6.0 4.0 4.5 4.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,1],[1,6],[6,4]])
const ascii97 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
const ascii98 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0 1.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10],[10,11]])
const ascii99 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
const ascii100 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0 3.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10],[10,11]])
const ascii101 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
const ascii102 = ([4.0 4.0 3.0 2.0 1.0 1.0 0.0 2.0; 3.0 4.0 5.0 5.0 4.0 0.0 1.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[7,8]])
const ascii103 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 -1.0 -1.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[1,9],[9,10],[10,11],[11,12]])
const ascii104 = ([4.0 4.0 3.0 1.0 0.0 0.0 0.0 1.0; 0.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[6,7],[7,8]])
const ascii105 = ([1.0 3.0 1.0 3.0 2.0 2.0 2.25 2.25 1.75 1.75; 0.0 0.0 3.0 3.0 0.0 3.0 3.75 4.25 4.25 3.75],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8],[8,9],[9,10],[10,7]])
const ascii106 = ([1.0 3.0 2.0 2.0 1.0 0.0 2.25 2.25 1.75 1.75; 3.0 3.0 3.0 0.0 -1.0 0.0 3.75 4.25 4.25 3.75],Array{Int64,1}[[1,2],[3,4],[4,5],[5,6],[7,8],[8,9],[9,10],[10,7]])
const ascii107 = ([0.0 1.0 1.0 0.0 4.0 2.0 1.0 3.0 4.0; 0.0 0.0 3.0 3.0 0.0 0.0 1.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[5,6],[6,7],[7,8],[8,9]])
const ascii108 = ([2.0 2.0 1.0 1.0 3.0; 0.0 5.0 5.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii109 = ([4.0 4.0 2.0 2.0 2.0 0.0 0.0 0.0; 0.0 3.0 2.0 0.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,6],[7,8]])
const ascii110 = ([3.0 3.0 1.0 1.0 1.0; 0.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii111 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
const ascii112 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
const ascii113 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
const ascii114 = ([0.0 2.0 1.0 1.0 1.0 2.0 3.0 4.0; 0.0 0.0 0.0 3.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[3,4],[5,6],[6,7],[7,8]])
const ascii115 = ([0.0 4.0 3.0 1.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 1.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
const ascii116 = ([1.0 3.0 2.0 2.0 2.0 3.0; 0.0 0.0 0.0 5.0 4.0 4.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
const ascii117 = ([0.0 1.0 1.0 2.0 3.0 4.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii118 = ([0.0 1.0 3.0 4.0; 3.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
const ascii119 = ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 3.0 2.0 0.0 2.0 0.0 2.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii120 = ([0.0 1.0 4.0 1.0 4.0; 3.0 3.0 0.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
const ascii121 = ([0.0 1.0 2.5 0.0 1.0 4.0; 3.0 3.0 1.5 0.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,6]])
const ascii122 = ([0.0 0.0 3.0 0.0 3.0 4.0; 2.0 3.0 3.0 0.0 0.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6]])
const ascii123 = ([2.5 2.0 2.0 1.5 2.0 2.0 2.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii124 = ([2.0 2.0; 0.0 5.0],Array{Int64,1}[[1,2]])
const ascii125 = ([1.5 2.0 2.0 2.5 2.0 2.0 1.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
const ascii126 = ([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4]])

"""
	hpcs::Array

"""
const hpcs = [
	ascii32,ascii33,ascii34,ascii35,ascii36,ascii37,ascii38,ascii39,ascii40,ascii41,
	ascii42,ascii43,ascii44,ascii45,ascii46,ascii47,ascii48,ascii49,ascii50,ascii51,
	ascii52,ascii53,ascii54,ascii55,ascii56,ascii57,ascii58,ascii59,ascii60,ascii61,
	ascii62,ascii63,ascii64,ascii65,ascii66,ascii67,ascii68,ascii69,ascii70,ascii71,
	ascii72,ascii73,ascii74,ascii75,ascii76,ascii77,ascii78,ascii79,ascii80,ascii81,
	ascii82,ascii83,ascii84,ascii85,ascii86,ascii87,ascii88,ascii89,ascii90,ascii91,
	ascii92,ascii93,ascii94,ascii95,ascii96,ascii97,ascii98,ascii99,ascii100,ascii101,
	ascii102,ascii103,ascii104,ascii105,ascii106,ascii107,ascii108,ascii109,ascii110,ascii111,
	ascii112,ascii113,ascii114,ascii115,ascii116,ascii117,ascii118,ascii119,ascii120,ascii121,
	ascii122,ascii123,ascii124,ascii125,ascii126 ]

"""
ascii_LAR = DataStructures.OrderedDict{Int,Lar.LAR}()

	ascii_LAR::{Int,Lar.LAR}

*Ordered dictionary* of printable ASCII codes as one-dimensional *LAR models* in 2D.

`Key` is the `integer` ASCII code between 32 and 126.
Font design: *Geometric Programming for Computer-Aided Design*, Wiley, 2003.
# Example
```
julia> Plasm.ascii_LAR[46]
([2.0 2.0 … 1.5 2.0; 0.0 0.5 … 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])

julia> Plasm.ascii_LAR[126]
([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
```
"""
const ascii_LAR = DataStructures.OrderedDict(zip(32:126,Plasm.hpcs))

# Attributes for `text` 2D graphics primitive.
# Reference to GKS ISO Graphics Standard (ISO/IEC 7942-4:1998)

global const  textalignment = "centre" #default value
global const textangle = pi/4 #default value
global const textwidth = 1.0 #default value
global const textheight = 1.0 #default value
global const textspacing = 0.25 #default value
global const fontwidth = 4.0 #default value
global const fontheight = 8.0 #default value
global const fontspacing = 1.0 #default value

"""
	charpols(mystring::String)::Array{LAR,1}

Return the `array` of `LAR` models, for chars in `mystring`.
"""
function charpols(mystring)
    return [ascii_LAR[code] for code in [Int(char) for char in mystring]]
end

"""
	charseq(mystring::String)::Array{Char,1}

# Example
```
julia> Plasm.charseq("PLaSM")
5-element Array{Char,1}:
 'P'
 'L'
 'a'
 'S'
 'M'
```
"""
function charseq(mystring)
	return [char for char in mystring]
end

"""
	text(mystring::String)::LAR

Compute the one-dim *LAR model* drawing the contents of `mystring`

# Example
```
julia> model = Plasm.text("PLaSM")
# output
([0.0 0.0 3.0 4.0 4.0 3.0 0.0 9.0 5.0 5.0 14.0 13.0 11.0 10.0
10.0 11.0 13.0 14.0 14.0 14.0 15.0 16.0 18.0 19.0 19.0 18.0 16.0 15.0 15.0 16.0 18.0
19.0 20.0 20.0 22.0 24.0 24.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 0.0 0.0 6.0 1.0 0.0 0.0 1.0
2.0 3.0 3.0 2.0 0.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 0.0 6.0 4.0
6.0 0.0],
Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[8,9],[9,10],[11,12],[12,13],[13,14],
[14,15],[15,16],[16,17],[17,18],[18,11],[19,20],[21,22],[22,23],[23,24],[24,25],[25,26],
[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[33,34],[34,35],[35,36],[36,37]])

julia> Plasm.view(model)
```
"""
function text(mystring)
	out = comp([ Lar.struct2lar, Lar.Struct, cat, distr,
			cons([ charpols, k(Lar.t(fontspacing+fontwidth,0)) ]),charseq ])(mystring)
	return out
end

"""
	a2a(mat::Matrix)(models::Array{LAR,1})::Struct.body

Local function, service to `textWithAttributes` implementation.
"""
function a2a(mat)
	function a2a0(models)
		assembly = []
		for model in models
			push!( assembly, Lar.Struct([ mat,model ]) )
		end
		assembly
	end
	return a2a0
end

"""
	translate(c::Number)(lar::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function translate(c)
	function translate0(lar)
		xs = lar[1][1,:]
		width = maximum(xs) - minimum(xs)
		apply(Lar.t(width/c,0))(lar)
	end
	return translate0
end

"""
	align(textalignment="centre"::String)(model::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function align(textalignment="centre")
	function align1(model)
		if ( textalignment == "centre" ) out=translate(-2)(model)
		elseif ( textalignment == "right" ) out=translate(-1)(model)
		elseif ( textalignment == "left" ) out=model
		end
		return out
	end
end

"""
 	cat(args)

Redefined locally, as service to `textWithAttributes` implementation.
"""
function cat(args)
	return reduce( (x,y) -> append!(x,y), args; init=[] )
end

"""
	textWithAttributes(textalignment='centre', textangle=0,
		textwidth=1.0, textheight=2.0, textspacing=0.25)(strand::String)::LAR

Partial implementation of the GKS's graphics primitive `text`.

# Example

```
Plasm.view(Plasm.textWithAttributes("left", pi/4)("PLaSM"))
```
"""
function textWithAttributes(textalignment="centre", textangle=0,
							textwidth=1.0, textheight=2.0, textspacing=0.25)
	function textWithAttributes(strand)
		id = x->x
		mat = Lar.s(textwidth/fontwidth,textheight/fontheight)
		comp([
		   apply(Lar.r(textangle)),
		   align(textalignment),
		   Lar.struct2lar,
		   Lar.Struct,
		   cat,
		   distr,
		   cons([ a2a(mat) ∘ charpols,
				k(Lar.t(textwidth+textspacing,0)) ]),
		   charseq ])(strand)
	end
end

"""
	embed(n::Int)(model::LAR)::LAR

Embed a `LAR` `model` of dimension ``d`` in a space ``R^{d+n}``.
The embedding is done by adding ``d`` zero coordinates to each vertex.

# Example

```
julia> square = Lar.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> Plasm.embed(1)(square)
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0; 0.0 0.0 0.0 0.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function embed(n)
	function embed0( larmodel )
		return [larmodel[1]; zeros(n,size(larmodel[1],2))], larmodel[2]
	end
	return embed0
end

"""
	numbering(model::LARmodel)

Display a *wireframe* of `model` cells, of any dimension, with their *ordinal number*.
Different `colors` and size are used for the various dimensional cells.

# Examples

```
model = Lar.cuboidGrid([3,4,2], true)
Plasm.view(Plasm.numbering()(model))

model = Lar.cuboidGrid([10,10], true)
Plasm.view(Plasm.numbering(1.5)(model))
```
"""
function numbering(numberSizeScaling=1)
	p = PyCall.pyimport("pyplasm")
	function numbering0(model)
		V,cells = model
		if size(V,1)==2
			V = Plasm.embed(1)(model)[1]
		end
		wireframe = Plasm.lar2hpc(V,cells[2])
		ns = numberSizeScaling
		gcode = Plasm.textWithAttributes("centre", 0, 0.1ns, 0.2ns, 0.025ns)
		scene = [wireframe]
		for (h,skel) in enumerate(cells)
			  colors = [p["GREEN"], p["YELLOW"], p["CYAN"], p["ORANGE"]]
		 	  nums = []
			  for (k,cell) in enumerate(skel)
				center = sum([V[:,v] for v in cell])/length(cell)
				code = Plasm.embed(1)( gcode(string(k)) )
				scaling = (0.6+0.1h,0.6+0.1h,1)
				push!(nums, Lar.struct2lar( Lar.Struct([
					Lar.t(center...), Lar.s(scaling...), code ]) ))
			end
			hpc = Plasm.lar2hpc(nums)
			push!( scene, p["COLOR"](colors[h])(hpc) )
		end
		p["STRUCT"]( scene )
	end
	return numbering0
end

"""
#	numbering(scaling=0.1)
#		(V::Lar.Points, copEV::Lar.ChainOp, copFE::Lar.ChainOp)::Lar.Hpc

Produce the numbered `Hpc` of `planar_arrangement()` 2D output.
Vertices in `V` are stored by row.

# Example

```julia
using LinearAlgebraicRepresentation
using Plasm, SparseArrays
Lar = LinearAlgebraicRepresentation

V,EV = Lar.randomcuboids(7, 1.0);
V = Plasm.normalize(V,flag=true);
W = convert(Lar.Points, V');
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
V, copEV, copFE = Lar.planar_arrangement(W, cop_EW);

Plasm.view( Plasm.numbering(0.05)((V, copEV, copFE)) )
```
"""
function numbering1(scaling=0.1)
	function numbering0(model::Tuple{Lar.Points,Lar.ChainOp,Lar.ChainOp})
		(V, copEV, copFE) = model
		VV = [[k] for k=1:size(V,1)]
		EV = [findnz(copEV[h,:])[1] for h=1:size(copEV,1)]
		FV = [collect(Set(cat(EV[e] for e in findnz(copFE[i,:])[1]))) for i=1:size(copFE,1)]
		FV = convert(Array{Array{Int64,1},1}, FV)
		model = (convert(Lar.Points, V'), Lar.Cells[VV,EV,FV])
		return Plasm.numbering(scaling)(model)
	end
	return numbering0
end


text:
using DiffEqBayes
using Distributions
using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions
using Plots

gr(fmt=:png)

g1 = @ode_def LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β

r0 = [1.0; 0.0; 0.0]
tspan = (0.0, 30.0)
p = [10.0,28.0,2.66]

prob = ODEProblem(g1,r0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(1,stop=30,length=30))
sig = 0.49
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(3)) for i in 1:length(t)]))

Plots.scatter(t, data[1,:],markersize=4,color=:purple)
Plots.scatter!(t, data[2,:],markersize=4,color=:yellow)
Plots.scatter!(t, data[3,:],markersize=4,color=:black)
plot!(sol)

cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,reltol=1e-5,abstol=1e-5)
plot(sim,vars=(0,1),linealpha=0.4)

cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,reltol=1e-6,abstol=1e-6)
plot(sim,vars=(0,1),linealpha=0.4)

cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,reltol=1e-8,abstol=1e-8)
plot(sim,vars=(0,1),linealpha=0.4)

priors = [Truncated(Normal(10,2),1,15),Truncated(Normal(30,5),1,45),Truncated(Normal(2.5,0.5),1,4)]

@time bayesian_result = stan_inference(prob,t,data,priors;reltol=1e-8,abstol=1e-8,vars=(StanODEData(),InverseGamma(3,2)))

plot_chain(bayesian_result)

@time bayesian_result_turing = turing_inference(prob,Tsit5(),t,data,priors)

plot_chain(bayesian_result_turing)


text:
function temperature_statistics(T)
    mean = 0
    std  = 0
    # Your code goes here!
    return mean, std
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(URIParser.parse_url, (Compat.ASCIIString,))
    precompile(URIParser.parse_authority, (Compat.ASCIIString, Bool,))
    precompile(URIParser.call, (Array{Any, 1}, Type{URIParser.URI}, URIParser.URI,))
    precompile(URIParser.escape_with, (Compat.ASCIIString, Compat.String,))
    precompile(URIParser.call, (Type{URIParser.URI}, Compat.ASCIIString, Compat.ASCIIString, UInt16, Compat.ASCIIString, Compat.ASCIIString, Compat.ASCIIString, Compat.ASCIIString, Bool,))
    precompile(URIParser.is_host_char, (Char,))
    precompile(URIParser.is_url_char, (Char,))
    precompile(URIParser.isnum, (Char,))
    precompile(URIParser.is_mark, (Char,))
    precompile(URIParser.is_userinfo_char, (Char,))
    precompile(URIParser.escape, (Compat.ASCIIString,))
    precompile(URIParser.call, (Type{URIParser.URI}, Compat.ASCIIString,))
    precompile(URIParser.ishex, (Char,))
end


text:
"""
Aggregate fitness reduces a set of individual search objectives to a single
scalar fitness value.
"""
immutable AggregateFitness
  score::Float
  objectives::Vector{Float}
end

type AggregateFitnessScheme <: FitnessScheme
  maximise::Bool
  weights::Vector{Float}

  AggregateFitnessScheme(m::Bool, w::Vector{Float}) = new(m, w)
end

"""
This fitness scheme operates by aggregating the weighted values of each of
the objectives into a single scalar value, which is to be either maximised
or minimised.

**Properties:**

* `weights::Vector{Float}`, a list holding the weights used for each of the
objective values upon aggregation.
* `maximise::Bool`, a flag indicating whether the aggregated fitness value
should be maximised.
"""
aggregate(s::Dict{Any, Any}) =
  # Need to handle weights.
  aggregate(Base.get(s, "maximise", true), s["weights"])

aggregate(maximise::Bool, weights::Vector{Float}) =
  AggregateFitnessScheme(maximise, weights)

uses(s::AggregateFitnessScheme) = AggregateFitness
maximise(s::AggregateFitnessScheme) = s.maximise
assign(s::AggregateFitnessScheme, o::Vector{Float}) =
  AggregateFitness(sum(o .* weights), o)


text:
function jac2cheb(c_jac::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPlan{FORWARD,T}) where T<:AbstractFloat
    M,N,nM₀,αN,K = getconstants(plan)
    i₁,i₂,j₁,j₂ = getindices(plan)
    p₁,p₂,rp,c₁,c₂,um,vm,cfs,θ,tempcos,tempsin,tempcosβsinα,tempmindices,cnαβ,cnmαβ = getplan(plan)

    v_cheb = zero(c_jac)

    k=0
    for j = 1:K
        k = j
        copyto!(cnmαβ,cnαβ)

        j₁[k] < nM₀ && (k-=1; break)

        @inbounds for i=i₁[k]:i₁[k+1]-1 v_cheb[i] += backward_recurrence(c_jac,j₂[k],θ[i],tempcos[i],tempsin[i],rp) end
        @inbounds for i=i₂[k+1]+1:i₂[k] v_cheb[i] += backward_recurrence(c_jac,j₂[k],θ[i],tempcos[i],tempsin[i],rp) end

        for m=0:M-1
            # Initialize diagonal N-scaling multiplied by local coefficients and zero out excess
            init_c₁c₂!(c₁,c₂.parent,cnmαβ,c_jac,j₁[k],j₂[k])

            # Apply planned DCT-I and DST-I in-place
            applyTN!(c₁,p₁);applyUN!(c₂,p₂)

            # Compute u_m(θ) and v_m(θ)
            compute_umvm!(um,vm,cfs,α,β,tempcos,tempsin,tempcosβsinα,m,θ,i₁[k+1]:i₂[k+1])

            # Multiply point-wise by u_m(θ) and v_m(θ) for valid indices
            @inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += um[i]*c₁[i] + vm[i]*c₂.parent[i] end

            # Update C_{n,m}^{α,β} by recurrence in m
            @inbounds for i=1:N+1 cnmαβ[i] /= 2(2i+α+β+m) end
        end
    end

    # Finish off recurrence
    @inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += backward_recurrence(c_jac,j₂[k+1],θ[i],tempcos[i],tempsin[i],rp) end

    # perform IDCT-I
    applyTNinv!(v_cheb,p₁)
end


text:
# Julia wrapper for header: /usr/include/scip/scip_dialog.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeDialog(scip, dialog, dialogcopy, dialogexec, dialogdesc, dialogfree, name, desc, issubmenu, dialogdata)
    ccall((:SCIPincludeDialog, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_DIALOG}}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring, UInt32, Ptr{SCIP_DIALOGDATA}), scip, dialog, dialogcopy, dialogexec, dialogdesc, dialogfree, name, desc, issubmenu, dialogdata)
end

function SCIPexistsDialog(scip, dialog)
    ccall((:SCIPexistsDialog, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}), scip, dialog)
end

function SCIPcaptureDialog(scip, dialog)
    ccall((:SCIPcaptureDialog, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}), scip, dialog)
end

function SCIPreleaseDialog(scip, dialog)
    ccall((:SCIPreleaseDialog, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_DIALOG}}), scip, dialog)
end

function SCIPsetRootDialog(scip, dialog)
    ccall((:SCIPsetRootDialog, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}), scip, dialog)
end

function SCIPgetRootDialog(scip)
    ccall((:SCIPgetRootDialog, libscip), Ptr{SCIP_DIALOG}, (Ptr{SCIP_},), scip)
end

function SCIPaddDialogEntry(scip, dialog, subdialog)
    ccall((:SCIPaddDialogEntry, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_DIALOG}, Ptr{SCIP_DIALOG}), scip, dialog, subdialog)
end

function SCIPaddDialogInputLine(scip, inputline)
    ccall((:SCIPaddDialogInputLine, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, inputline)
end

function SCIPaddDialogHistoryLine(scip, inputline)
    ccall((:SCIPaddDialogHistoryLine, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, inputline)
end

function SCIPstartInteraction(scip)
    ccall((:SCIPstartInteraction, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end


text:
# To do: Add a structure a make a surchage of 'plot' of Makie

"""
    iplsr(Xtrain, Ytrain, X, Y; nint = 5, score = rmsep, nlv)
Interval PLSR (iPLS) (Nørgaard et al. 2000)
* `Xtrain` : Training X-data.
* `Ytrain` : Training Y-data.
* `X` : Validation X-data.
* `Y` : Validation Y-data.
* `nint` : Nb. intervals.
* `score` : Function computing the prediction score (= error rate; e.g. msep).
* `nlv` : Nb. latent variables (LVs) in the PLSR models.

The range `1:p` (where `p` is the nb. of columns of `Xtrain`) is segmented
to `nint` column-intervals of (when possible) equal size. Then, the validation
score is computed for each of the `nint` PLSR models and compared to the
one of the overal `Xtrain`.

## References

- Nørgaard, L., Saudland, A., Wagner, J., Nielsen, J.P., Munck, L.,
    Engelsen, S.B., 2000. Interval Partial Least-Squares Regression (iPLS):
    A Comparative Chemometric Study with an Example from Near-Infrared
    Spectroscopy. Appl Spectrosc 54, 413–419. https://doi.org/10.1366/0003702001949500
"""
function iplsr(Xtrain, Ytrain, X, Y;
        nint = 5, score = rmsep, nlv)
    Y = ensure_mat(Y)
    p = size(Xtrain, 2)
    q = size(Y, 2)
    z = collect(round.(range(1, p + 1, length = nint + 1)))
    int = Int64.([z[1:nint] z[2:(nint + 1)] .- 1])
    pred = similar(Y)
    res = list(nint, Matrix{Float64})
    @inbounds for i = 1:nint
        u = int[i, 1]:int[i, 2]
        fm = plskern(vcol(Xtrain, u), Ytrain; nlv = nlv)
        pred .= predict(fm, vcol(X, u)).pred
        res[i] = score(pred, Y)
    end
    res = reduce(vcat, res)
    fm = plskern(Xtrain, Ytrain; nlv = nlv)
    pred = predict(fm, X).pred
    res0 = score(pred, Y)
    dat = DataFrame(int, [:lo, :up])
    namy = map(string, repeat(["y"], q), 1:q)
    res = DataFrame(res, Symbol.(namy))
    res = hcat(dat, res)
    res0 = DataFrame(res0, Symbol.(namy))
    (res = res, res0 = res0, int = int)
end


text:
using ValidatedNumerics

using Test

@testset "Load and briefly test packages" begin

    @testset "IntervalArithmetic" begin
        @test 1..2 == Interval(1, 2)
    end

    @testset "IntervalRootFinding" begin
        roots = IntervalRootFinding.roots(x->x^2 - 2, -10..10)
        @test length(roots) == 2
    end

    @testset "IntervalConstraintProgramming" begin
        S = @constraint x^2 + y^2 <= 1
        X = (-∞..∞) × (-∞..∞)
        paving = pave(S, X, 1.0)

        @test length(paving.inner) == 3
        @test length(paving.boundary) == 7
    end

    @testset "IntervalContractors" begin
        X = 1..4
        Y = 5..10

        Z = X + Y

        Z = Z ∩ (13..15)

        Z, X, Y = plus_rev(Z, X, Y)

        @test X == 3..4
        @test Y == 9..10
    end

    @testset "IntervalOptimisation" begin
        f(x) = x^4 - 3x^3 + 2x

        globalmin, minimisers = minimise(f, -1e10..1e10, tol=1e-5)

        @test globalmin ⊆ (-4.15.. -4.14)
    end

    @testset "TaylorModels" begin
        x0 = Interval(0.0)
        ii0 = Interval(-0.5, 0.5)

        tpol = exp( Taylor1(2) )
        @test TaylorModels.bound_taylor1( tpol, ii0) == tpol(ii0.lo) .. tpol(ii0.hi)
    end
end


text:
module PhysicsBox

using SpecialFunctions

include("./Integrators.jl")

include("./RandomDistributions.jl")

include("./PhysicsFunctions.jl")

end # module


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

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

# ╔═╡ ae20362a-4c60-11eb-2e5a-5b3904d7ef11
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add("InformationDistances")
	Pkg.add("CodecXz")
	Pkg.add("CodecBzip2")
	Pkg.add("StringDistances")
	Pkg.add("BioServices")
	Pkg.add("EzXML")
	Pkg.add("DataFrames")
	Pkg.add("PlutoUI")
	Pkg.add("Plots")
	Pkg.add("PyPlot")
	Pkg.add(name="PhyloPlots", rev="master")
	Pkg.add(name="PhyloNetworks", version="0.12.0")
end

# ╔═╡ 38a62858-4c67-11eb-312a-cb83457c1b8a
begin
	using InformationDistances
	using CodecXz: XzCompressor
	using CodecBzip2: Bzip2Compressor

	# Additional string distances for reference
	import StringDistances

	# Loading Mitochondrial DNA
	import BioServices
	import EzXML
	import DataFrames

	# Pluto UI elements
	import PlutoUI

	# Plotting heatmaps and phylogentic trees
	import Plots
	import PhyloNetworks
	import PhyloPlots
end

# ╔═╡ a361f525-31f3-445d-9ad1-5f4344f8e2dd
PlutoUI.TableOfContents()

# ╔═╡ 8fb6f778-63d3-11eb-37bd-81c244553b6c
md"""
# Mitochondrial Genome Phylogency

This notebook tries to replicate a result from [1] where the autors used the [Normalized Compression Distance](https://en.wikipedia.org/wiki/Normalized_compression_distance) to calculate pairwaise distaces between the mitochondrial DNA of various animals and then uses these distances to create a [phylogenetic tree](https://en.wikipedia.org/wiki/Phylogenetic_tree).

"""

# ╔═╡ dfaaa16b-04ee-4a28-bb62-5f7c93fc4c4c
Markdown.MD(Markdown.Admonition("note", "Note", [
md"""The results here might be slightly different from the paper as we do not use the exact same string compressors, and the algorithm used to reconstruct the phylogenetic tree might also be slightly different.
"""]))

# ╔═╡ a08db6ca-63dc-11eb-2d6f-99961e06ef3c
md"""
## Setup
We load a bunch of packages, this might take a while...
"""

# ╔═╡ dca8b740-63dc-11eb-2bd6-2753b38c90f8
md"""
We start by loading the mitochondiral DNA of 20 different animals from [nuccore](https://www.ncbi.nlm.nih.gov/nucleotide/). Using the `efetch` function from [BioServices.jl](https://github.com/BioJulia/BioServices.jl) this is a very easy thing to do, one has just to know the correct sequence id for each animal.
"""

# ╔═╡ 2b43945a-53f6-11eb-3555-7509d07e5fe4
mitochondrial_dna_sources = DataFrames.DataFrame([
	(name = "rat", sequence_id = "NC_001665.2"),
	(name = "house mouse", sequence_id = "NC_005089.1"),
	(name = "gray seal", sequence_id = "NC_001602.1"),
	(name = "harbor seal", sequence_id = "NC_001325.1"),
	(name = "cat", sequence_id = "NC_001700.1"),
	(name = "white rhino", sequence_id = "NC_001808.1"),
	(name = "horse", sequence_id = "NC_001640.1"),
	(name = "finback whale", sequence_id = "NC_001321.1"),
	(name = "blue whale", sequence_id = "NC_001601.1"),
	(name = "cow", sequence_id = "NC_006853.1"),
	(name = "gibbon", sequence_id = "NC_002082.1"),
	(name = "gorilla", sequence_id = "NC_001645.1"),
	(name = "human", sequence_id = "NC_012920.1"),
	(name = "chimpanzee", sequence_id = "NC_001643.1"),
	(name = "pigmy chimpanzee", sequence_id = "NC_001644.1"),
	(name = "orangutan", sequence_id = "NC_001646.1"),
	(name = "Sumatran orangutan", sequence_id = "NC_002083.1"),
	(name = "opossum", sequence_id = "NC_001610.1"),
	(name = "wallaroo", sequence_id = "NC_001794.1"),
	(name = "platypus", sequence_id = "NC_000891.1"),
])

# ╔═╡ 188e2e0e-5404-11eb-04b8-936f9971bd95
animals = map(eachrow(mitochondrial_dna_sources)) do row
	res = EzXML.parsexml(BioServices.EUtils.efetch(
			db = "nuccore", id=row.sequence_id, retmode="xml").body)
	latin_name = EzXML.nodecontent(findfirst("//GBSeq_organism", res))
	dna = EzXML.nodecontent(findfirst("//GBSeq_sequence", res))
	return (name=row.name, latin_name=latin_name, DNA=dna)
end |> DataFrames.DataFrame

# ╔═╡ 3d59590e-63de-11eb-1275-d3545baf1a6f
md"""
## Calculating pairwise distances

The Normalized Compressor Distance uses a loseless string compressor - changing the compressor should have an influence on the quality of the result.

For comparison we also use the [Levenshtein Distance](https://en.wikipedia.org/wiki/Levenshtein_distance) which is provided by [StringDistances.jl](https://github.com/matthieugomez/StringDistances.jl). Calculating this istance might need a few minutes.
"""

# ╔═╡ 9aca90d2-45dc-4424-a6ae-769c62e80bf4
Markdown.MD(Markdown.Admonition("tip", "Tip", [
md"""When viewing the static html version of this notebook it might not be possible to choose another distance - try running the interactive notebook instead.
"""]))

# ╔═╡ 564ac642-4c95-11eb-1329-a15bb5b6bef0
@bind distance_key PlutoUI.Radio([
		"default" => "XzCompressor(level = 9) (default)",
		"XzCompressor_level1" => "XzCompressor(level = 1)",
		"Bzip2Compressor" => "Bzip2Compressor",
		"LibDeflateCompressor" => "LibDeflateCompressor",
		"LevenshteinDistance" => "Levenshtein Distance"
], default="default")

# ╔═╡ 5f9ccc3a-4c96-11eb-3857-65954810d179
distance = Dict(
	"default" =>
		NormalizedCompressionDistance(),
	"XzCompressor_level1" =>
		NormalizedCompressionDistance(CodecCompressor{XzCompressor}(; level=1)),
	"Bzip2Compressor" =>
		NormalizedCompressionDistance(CodecCompressor{Bzip2Compressor}(;workfactor=250)),
	"LibDeflateCompressor" =>
		NormalizedCompressionDistance(LibDeflateCompressor()),
	"LevenshteinDistance" => StringDistances.Levenshtein()
)[distance_key]

# ╔═╡ f4c5057a-63de-11eb-0745-293580109763
md"""
We can now compare the distance between some animals. One would expect that rats and mice are much closer than either of them to seals.
"""

# ╔═╡ 9e85ae72-4c97-11eb-2aa4-dfed6f457274
(rat, mouse, seal) = (animals[1, :], animals[2, :], animals[3, :])

# ╔═╡ 49d74218-5416-11eb-0c31-9f1f42fe0498
distance(rat.DNA, mouse.DNA)

# ╔═╡ cbb4f58c-4c6e-11eb-3a82-a38fdbeb1836
distance(rat.DNA, seal.DNA)

# ╔═╡ 26226824-63df-11eb-3650-bdf6ab44a134
md"""
We can now calulate all pairwise distances and visualize the matrix. As close animal groups are already next to each other in the input list, we can immediately spot some clusters.
"""

# ╔═╡ 5894c63e-4c70-11eb-0978-71432cc5cda1
distance_matrix = Float64[distance(x, y) for x in animals[!, :DNA], y in animals[!, :DNA]]

# ╔═╡ 41a44156-4c76-11eb-3675-11ca0568656b
begin
	Plots.pyplot()
	animal_names = replace.(animals[!, :name], ' ' => "\\ ")
	Plots.heatmap(
		animal_names,
		animal_names,
		distance_matrix;
		xtickfontrotation=90,
		yflip=true,
		yticks=:all,
		xticks=:all,
		right_margin=2.5(Plots.mm),
	)
end

# ╔═╡ 51fd3fea-63d6-11eb-2263-196018bcbff6
md"""
## Creating a phylogenetic tree

Finally, we can pass the distance matrix to the `nj!` function from [PhyloNetworks.jl](https://github.com/crsl4/PhyloNetworks.jl) to create a phylogenetic tree and then plot this tree using [PhyloPlots.jl](https://github.com/cecileane/PhyloPlots.jl).
"""

# ╔═╡ 8b6e1a10-4c70-11eb-07d9-113e4f2866a0
begin
	tree = PhyloNetworks.nj!(copy(distance_matrix), collect(animals[!, :name]))
	PhyloPlots.plot(tree)
end

# ╔═╡ ae2812f8-63d2-11eb-2e31-a7e555cb1e24
md"""
## References
[1]: [Li, Ming, Xin Chen, Xin Li, Bin Ma, and Paul MB Vitányi. "The similarity metric." IEEE transactions on Information Theory 50, no. 12 (2004): 3250-3264.](www.google.ch)
"""

# ╔═╡ Cell order:
# ╟─a361f525-31f3-445d-9ad1-5f4344f8e2dd
# ╟─8fb6f778-63d3-11eb-37bd-81c244553b6c
# ╟─dfaaa16b-04ee-4a28-bb62-5f7c93fc4c4c
# ╟─a08db6ca-63dc-11eb-2d6f-99961e06ef3c
# ╠═ae20362a-4c60-11eb-2e5a-5b3904d7ef11
# ╠═38a62858-4c67-11eb-312a-cb83457c1b8a
# ╟─dca8b740-63dc-11eb-2bd6-2753b38c90f8
# ╠═2b43945a-53f6-11eb-3555-7509d07e5fe4
# ╠═188e2e0e-5404-11eb-04b8-936f9971bd95
# ╟─3d59590e-63de-11eb-1275-d3545baf1a6f
# ╟─9aca90d2-45dc-4424-a6ae-769c62e80bf4
# ╠═564ac642-4c95-11eb-1329-a15bb5b6bef0
# ╠═5f9ccc3a-4c96-11eb-3857-65954810d179
# ╟─f4c5057a-63de-11eb-0745-293580109763
# ╠═9e85ae72-4c97-11eb-2aa4-dfed6f457274
# ╠═49d74218-5416-11eb-0c31-9f1f42fe0498
# ╠═cbb4f58c-4c6e-11eb-3a82-a38fdbeb1836
# ╟─26226824-63df-11eb-3650-bdf6ab44a134
# ╠═5894c63e-4c70-11eb-0978-71432cc5cda1
# ╠═41a44156-4c76-11eb-3675-11ca0568656b
# ╟─51fd3fea-63d6-11eb-2263-196018bcbff6
# ╠═8b6e1a10-4c70-11eb-07d9-113e4f2866a0
# ╟─ae2812f8-63d2-11eb-2e31-a7e555cb1e24


text:
using OmegaGrav
using CSV, DataFrames
using Dierckx
using Plots, LaTeXStrings
# %% Compute the matter power spectrum using CLASS
include("compute_pk_class.jl")
# %% Compute `Ωgrav = Ωm * W / 2` and Ωtherm at redshifts of the data points
d = CSV.read("data/d16_Omega_th_data.csv", DataFrame)
nred = size(d)[1]
redshift = zeros(nred + 1)
Ωghalo = zeros(nred + 1)
Ωth = zeros(nred + 1)
Ωthu = zeros(nred + 1)
Ωthl = zeros(nred + 1)
ΩthB1 = zeros(nred + 1)
for ired = 1:nred+1
   if ired == 1
      z = 0 # Add z=0 which is not included in the data file
   else
      z = d.z[ired-1]
   end
   redshift[ired] = z
   # Define functions to return linear baryon+CDM power spectrum
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pkcb_class(kovh) = cosmo.pk_cb_lin(kovh * h0, z) * h0^3
   # Spline interpolate in log(k)
   lnk = log(1e-4):0.05:log(100)
   pkcb = Spline1D(lnk, pkcb_class.(exp.(lnk)))
   # %% Compute Ωgrav from Halos, excluding the neutrino contribution
   Ωghalo[ired] = ograv_halo(x -> pkcb(log(x)), z, Ωm, Ωcb)
   # %% Compute Ωtherm from Halos, excluding the neutrino contribution
   Ωth[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb)
   # %% Compute Ωtherm for upper and lower 68% confidence level in B
   Ωthl[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.315)
   Ωthu[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.221)
   # %% Compute Ωtherm for no mass bias case (B=1) for comparison
   ΩthB1[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1)
   # %% Compute Ωtherm for Komatsu-Seljak profile
   # ΩthB1[ired] = otherm_ks(x -> pkcb(log(x)), z, Ωm, 0.157, Ωcb)
end

#%% Plot results and save to figure3.pdf
fb = params["omega_b"] / (params["omega_cdm"] + params["omega_b"])
ii = findall(x -> x < 1, d.z)
p = scatter(
   d.z[ii],
   -1.5 * d.Omega_th[ii] ./ Ωghalo[ii.+1] / fb,
   yerror = (
      -1.5 * (d.Omega_th .- d.Omega_th_low) ./ Ωghalo[2:end] / fb,
      -1.5 * (d.Omega_th_up .- d.Omega_th) ./ Ωghalo[2:end] / fb,
   ),
   ms = 5,
   c = 1,
   ylims = [0, 1.5],
   xlab = "Redshift, z",
   ylab = L"\Omega_{th}/\Omega_W^{ref}",
   lab = "Data",
   legend = :topright,
   legendfontsize = 12,
   labelfontsize = 15,
)
ii = findall(x -> x > 1, d.z)
u = zeros(length(ii))
v = -0.1 * ones(length(ii))
p = quiver!(
   d.z[ii],
   -1.5 * d.Omega_th_up[ii] ./ Ωghalo[ii.+1] / fb,
   quiver = (u, v),
   c = 1,
)
p = plot!(
   redshift,
   -1.5 * Ωth ./ Ωghalo / fb,
   ribbon = (
      -1.5 * (Ωth - Ωthl) ./ Ωghalo / fb,
      -1.5 * (Ωthu - Ωth) ./ Ωghalo / fb,
   ),
   lab = "Halo Model",
   lw = 5,
   c = 1,
)
p = plot!(
   redshift,
   -1.5 * ΩthB1 ./ Ωghalo / fb,
   ls = :dashdot,
   lab = "Halo Model (B=1)",
   lw = 5,
)
p = hline!([1], c = :grey, lab="", lw = 2)
savefig("figure3.pdf")
display(p)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()


text:
abstract type TimeEvolution end

const UNDEF_SAMPLEPETIME = -1 # For handling promotion of Matrix to LTISystem

struct Discrete{T} <: TimeEvolution
    Ts::T
    function Discrete{T}(Ts::T) where T
        if Ts <= 0 && Ts != UNDEF_SAMPLEPETIME
            throw(ErrorException("Creating a continuous time system by setting sample time to 0 is no longer supported."))
        end
        new{T}(Ts)
    end
end
Discrete{T}(x) where T = Discrete{T}(T(x))

struct Continuous <: TimeEvolution end
# Basic Constructors
Discrete(x::T) where T<:Number = Discrete{T}(x)
Continuous(x::Continuous) = x
# Simple converseion
Discrete{T}(x::Discrete) where T = Discrete{T}(x.Ts)

undef_sampletime(::Type{Discrete{T}}) where T = Discrete{T}(UNDEF_SAMPLEPETIME)
undef_sampletime(::Type{Continuous}) where T = Continuous()

# Promotion
# Fallback to give useful error
Base.promote_rule(::Type{<:Discrete}, ::Type{Continuous}) = throw(ErrorException("Sampling time mismatch"))
Base.promote_rule(::Type{Discrete{T1}}, ::Type{Discrete{T2}}) where {T1,T2}= Discrete{promote_type(T1,T2)}

Base.convert(::Type{Discrete{T1}}, x::Discrete{T2}) where {T1,T2} = Discrete{T1}(x.Ts)

# Promoting two or more systems systems should promote sample times
common_timeevol(x::TimeEvolution) = x
common_timeevol(x::TimeEvolution, y::TimeEvolution) = throw(ErrorException("Sampling time mismatch"))
common_timeevol(x::TimeEvolution, y::TimeEvolution, z...) = common_timeevol(common_timeevol(x, y), z...)
common_timeevol(a::Base.Generator) = reduce(common_timeevol, a)

function common_timeevol(x::Discrete{T1}, y::Discrete{T2}) where {T1,T2}
    if x != y && x.Ts != UNDEF_SAMPLEPETIME && y.Ts != UNDEF_SAMPLEPETIME
         throw(ErrorException("Sampling time mismatch"))
    end

    if x.Ts == UNDEF_SAMPLEPETIME
        return Discrete{promote_type(T1,T2)}(y)
    else
        return Discrete{promote_type(T1,T2)}(x)
    end
end

common_timeevol(x::Continuous, ys::Continuous...) = Continuous()

# Check equality
==(x::TimeEvolution, y::TimeEvolution) = false
==(x::Discrete, y::Discrete) = (x.Ts == y.Ts)
==(::Continuous, ::Continuous) = true

isapprox(x::TimeEvolution, y::TimeEvolution, args...; kwargs...) = false
isapprox(x::Discrete, y::Discrete, args...; kwargs...) = isapprox(x.Ts, y.Ts, args...; kwargs...)
isapprox(::Continuous, ::Continuous, args...; kwargs...) = true


text:
using Logging

@testset "return type" begin
    @jyro function m()::Int
        x ~ Normal()
        return 2
    end

    @test typeof(m) == MinijyroModel
    m_return_types = Base.return_types(m.model_fn)
    @test length(m_return_types) == 1
    @test m_return_types[1] == Dict{Any,Any}

    t = m()
    @test haskey(t, :_return)
    @test t[:_return] == 2

    @jyro function m2()
        x ~ Normal()
    end

    m_return_types = Base.return_types(m2.model_fn)
    @test length(m_return_types) == 1
    @test m_return_types[1] == Dict{Any,Any}

    t = m2()
    @test haskey(t, :_return)
    @test t[:_return] == nothing
end

@testset "multiple return statements" begin
    @jyro function m(switch::Bool)
        x ~ Normal()
        if switch
            return 1
        end

        y ~ Normal()
        return 2
    end

    trace!(m)
    t = m(true)
    @test t[:_return] == 1
    @test haskey(t[:msgs], :x)
    @test !haskey(t[:msgs], :y)

    t = m(false)
    @test t[:_return] == 2
    @test haskey(t[:msgs], :x)
    @test haskey(t[:msgs], :y)
end


text:
@doc raw"""
    lentz_thompson(b₀, a, b)

Lentz–Thompson algorithm for the forward evaluation of continued
fractions:

```math
f_n = b_0 + \frac{a_1}{b_1+}\frac{a_2}{b_2+}...\frac{a_{n-1}}{b_{n-1}+}\frac{a_n}{b_n}
```

As described in

- Barnett, A. R. (1996). The Calculation of Spherical Bessel and
  Coulomb Functions. In (Eds.), Computational Atomic Physics
  (pp. 181–202). Springer Berlin Heidelberg.

"""
function lentz_thompson(b₀::T, a, b;
                        max_iter=20000,
                        ϵ = eps(T), tol=100eps(T),
                        verbosity=0) where T
    val_or_ϵ(x) = abs(x) < ϵ ? ϵ : x

    f = val_or_ϵ(b₀)
    C = f
    D = zero(T)
    δ = zero(T)

    fmt = if verbosity>1
        printfmtln("{1:>4s} {2:>10s} {3:>10s} {4:>10s} {5:>10s} {6:>10s} {7:>10s} {8:>10s}",
                   "n", "a", "b", "C", "D", "|Δ-1|", "f", "δf")
        FormatExpr("{1:>4d} {2:10.3e} {3:10.3e} {4:10.3e} {5:10.3e} {6:10.3e} {7:10.3e} {8:10.3e}")
    else
        nothing
    end

    s = true

    for n = 1:max_iter
        aₙ = a(n)
        bₙ = b(n)
        C = val_or_ϵ(bₙ + aₙ/C)
        D = inv(val_or_ϵ(bₙ + aₙ*D))
        if isreal(D) && D < 0
            # See section 3.2 of
            #
            # - Barnett, A., Feng, D., Steed, J., & Goldfarb, L. (1974). Coulomb
            #   wave functions for all real $\eta$ and ρ. Computer Physics
            #   Communications, 8(5),
            #   377–395. http://dx.doi.org/10.1016/0010-4655(74)90013-7
            s = !s
        end
        Δ = C*D
        fn = f*Δ
        δ = abs(Δ-1)
        verbosity>1 && printfmtln(fmt, n, aₙ, bₙ, C, D, δ, fn, fn-f)
        f = fn
        δ < tol && return f,n,δ,s,true
    end

    verbosity > 0 && @warn "Lentz–Thompson did not converge in $(max_iter) iterations, |Δ-1| = $(δ)"
    f,max_iter,δ,s,false
end


text:
using Pkg
Pkg.activate(".")

using Revise
using ExponentialFamilies
using BenchmarkTools
using Distributions
@btime logpdf(normal(0,1), 0.0)

const x = randn(10000);

logpdf(iid(normal(0,1)),x)

dist = zeros(10000)
expf = zeros(10000)

for j=1:10000
    @views dist[j] = @elapsed sum(Distributions.logpdf.(Normal(),x[1:j]))
    @views expf[j] = @elapsed logpdf(iid(normal(0,1)),x[1:j])
end

using Plots

js = 1:10000
plot(js, dist ./ js * 1e9, label="Using Distributions.jl", legend=:bottomleft)
plot!(js, expf ./ js * 1e9, label="Using ExponentialFamilies.jl")
ylims!(0,10) #quantile(dist ./ js * 1e9, 0.99))
xlabel!("Number of Observations")
ylabel!("Nanoseconds per Observation")
# using Cthulhu


text:
# ---
# title: 1000. Minimum Cost to Merge Stones
# id: problem1000
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-cost-to-merge-stones/description/>
# hidden: true
# ---
#
# There are `N` piles of stones arranged in a row.  The `i`-th pile has
# `stones[i]` stones.
#
# A _move_ consists of merging **exactly  `K` consecutive** piles into one pile,
# and the cost of this move is equal to the total number of stones in these `K`
# piles.
#
# Find the minimum cost to merge all piles of stones into one pile.  If it is
# impossible, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: stones = [3,2,4,1], K = 2
#     Output: 20
#     Explanation:
#     We start with [3, 2, 4, 1].
#     We merge [3, 2] for a cost of 5, and we are left with [5, 4, 1].
#     We merge [4, 1] for a cost of 5, and we are left with [5, 5].
#     We merge [5, 5] for a cost of 10, and we are left with [10].
#     The total cost was 20, and this is the minimum possible.
#
#
# **Example 2:**
#
#
#
#     Input: stones = [3,2,4,1], K = 3
#     Output: -1
#     Explanation: After any merge operation, there are 2 piles left, and we can't merge anymore.  So the task is impossible.
#
#
# **Example 3:**
#
#
#
#     Input: stones = [3,5,1,2,6], K = 3
#     Output: 25
#     Explanation:
#     We start with [3, 5, 1, 2, 6].
#     We merge [5, 1, 2] for a cost of 8, and we are left with [3, 8, 6].
#     We merge [3, 8, 6] for a cost of 17, and we are left with [17].
#     The total cost was 25, and this is the minimum possible.
#
#
#
#
# **Note:**
#
#   * `1 <= stones.length <= 30`
#   * `2 <= K <= 30`
#   * `1 <= stones[i] <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
using ParametricMachinesDemos, Flux

poly(x, y) = (2x-1)^2 + 2y + x * y - 3
poly((x, y)) = poly(x, y)
rg = 0:0.01:1
N = length(rg)
flat = hcat(repeat(rg, inner=N), repeat(rg, outer=N))
truth = map(poly, eachrow(flat))

## Let us generate a `6 x 6` trainig grid.

N_train = 6
rg_train = range(0, 1, length=N_train)
X_train = hcat(repeat(rg_train, inner=N_train), repeat(rg_train, outer=N_train)) |> permutedims
Y_train = map(poly, eachcol(X_train)) |> permutedims

dimensions = [2, 4, 4, 4]

machine = DenseMachine(dimensions, sigmoid)

model = Flux.Chain(machine, Dense(sum(dimensions), 1)) |> f64

function loss(X_train, Y_train)
    Flux.Losses.mse(model(X_train), Y_train)
end

opt = ADAM(0.1)
ps = Flux.params(model)

# check that learning happens correctly

N = 14
rg = range(0, 1, length=N)
X = hcat(repeat(rg, inner=N), repeat(rg, outer=N)) |> permutedims
Y = map(poly, eachcol(X)) |> permutedims

for i in 1:10000
    gs = gradient(ps) do
        loss(X_train, Y_train)
    end
    Flux.Optimise.update!(opt, ps, gs)
    if i % 500 == 0
        @show loss(X_train, Y_train)
        @show loss(X, Y)
    end
end


text:
using BioSeq
type Hsp
	hsp_num::Int
	bitScore::Float64
	evalue::Float64
	queryFrom::Int
	queryTo::Int
	queryFrame::Int
	hitFrame::Int
	identity::Int
	positive::Int
	gaps::Int
	alignLen::Int

	qseq::Array{AminoAcid,1}
	qseq_str::AbstractString

	hseq::Array{AminoAcid,1}
	midline::Array{AminoAcid,1}

end

type Hit
	hit_num::Int
	id::AbstractString
	def::AbstractString
	accession::AbstractString
	len::Int

	hsps::Array{Hsp,1}
end

function fastarepresentation(hit::Hit)
	fasta = ""

	for hsp in hit.hsps
		fasta = (hit.id, hsp.qseq_str)
	end
	return fasta
end

function check_threshold(threshold::Float64, hit::Hit)
	# check single hsp, for now
	for hsp in hit.hsps
		if hsp.evalue < threshold
			return true
		else
			return false
		end
	end
end


text:
CoordinateSystemCardType = CardType{:coordinate_system}()
const CoordinateSystemCard = NastranCard{CoordinateSystemCardType}

type CORD1R <: CoordinateSystemCard
    id::Int64
    gridA::Int64
    gridB::Int64
    gridC::Int64
end
function convert(::Type{CORD1R},card)
    CORD1R(card[2]::Int64,
           card[3]::Int64,
           card[4]::Int64,
           card[5]::Int64)
end

type CORD2R <: CoordinateSystemCard
    id::Int64
    csys_id::Int64
    gridAx::Float64
    gridAy::Float64
    gridAz::Float64
    gridBx::Float64
    gridBy::Float64
    gridBz::Float64
    gridCx::Float64
    gridCy::Float64
    gridCz::Float64
end
function convert(::Type{CORD2R},card)
    CORD2R(card[2]::Int64,
           card[3] == "" ? 0 : card[3]::Int64,
           card[4]::Float64,
           card[5]::Float64,
           card[6]::Float64,
           card[7]::Float64,
           card[8]::Float64,
           card[9]::Float64,
           card[10]::Float64,
           card[11]::Float64,
           card[12]::Float64)
end


text:
# OptiGas: Optimize Gas Network using Benders Algorithm
# Version: 2.0
# Author: Edward J. Xu, edxu96@outlook.com
# Date: March 20th, 2019
push!(LOAD_PATH, "$(homedir())/Desktop/OptiGas, DTU42136")
cd("$(homedir())/Desktop/OptiGas, DTU42136")
# using BendersMilp_EDXU
using Benders
using LinearAlgebra
# -------------------------------------------------------------------
# 1. Parameters
vec_nameNodes = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J"]
numNodes = length(vec_nameNodes)
vec_xNode = [1    3    8    2    4    5    6    7    9    9]
vec_yNode = [1    2    1    7    9    5    3    7    9    4]
mat_arcTwoNodes = [0    1    0    0    0    0    0    0    0    0;
                   0    0    0    0    0    1    0    0    0    0;
                   0    0    0    0    0    0    0    0    0    1;
                   1    0    0    0    1    0    0    0    0    0;
                   0    0    0    0    0    1    0    0    0    0;
                   1    0    0    1    0    0    0    1    0    0;
                   0    1    1    0    0    0    0    0    0    0;
                   0    0    0    0    0    0    0    0    1    0;
                   0    0    0    0    1    0    0    0    0    0;
                   0    0    0    0    0    1    1    0    1    0]
vec_netEject =  [0    0   0    0    0   90    0    0   0   80 ]
vec_netInject = [10   30  20   10   40   0   25   20  15    0 ]
mat_distance = zeros(Float64, numNodes, numNodes)
for m = 1: numNodes
    for n = 1: numNodes
        mat_distance[m, n] = floor(sqrt( (vec_xNode[m] - vec_xNode[n]) * (vec_xNode[m] - vec_xNode[n]) +
                                   (vec_yNode[m] - vec_yNode[n]) * (vec_yNode[m] - vec_yNode[n]) ) )
    end
end
mat_fixedCost = zeros(Float64, numNodes, numNodes)
for m = 1: numNodes
    for n = 1: numNodes
        mat_fixedCost[m, n] = 10 * mat_distance[m, n]
    end
end
# ----------------------------------------------------------------------------------------------------------------------
# Question 4
# Fomulation of Matrix A
mat_a1_1 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_a1_1[mn, mn] = -1
end
# mat_a2_1 = zeros(numNodes^2, numNodes^2)
mat_a3_1 = zeros(numNodes, numNodes^2)
seq_zeroTenNinety = collect(0:10:90)
for m = 1: numNodes
    mat_a3_1[m, (10 * m - 9): (10 * m)] = repeat([-1], 10)  # sent out
    mat_a3_1[m, (seq_zeroTenNinety + repeat([m],10))] = repeat([1], 10)  # sent in
    mat_a3_1[:, (11 * m - 10)] = repeat([0], 10)  # self-sending doesn't count
end
mat_a_1 = vcat(mat_a1_1, mat_a3_1)
# Fomulation of Matrix B
mat_b1_1 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_b1_1[mn, mn] = 170
end
# mat_b2_1 = Diagnal(repeat([1], numNodes^2))
mat_b3_1 = zeros(numNodes, numNodes^2)
mat_b_1 = vcat(mat_b1_1, mat_b3_1)
# Fomulation of Vector b
vec_b1_1 = zeros(numNodes^2)
# vec_b2_1 = zeros(numNodes^2)
# for m = 1: numNodes
#     for n = 1: numNodes
#         vec_b2_1[10 * (m - 1) + n] = mat_arcTwoNodes[m, n]
#     end
# end
vec_b3_1 = zeros(numNodes)
for n = 1: numNodes
    vec_b3_1[n] = vec_netInject[n] - vec_netEject[n]
end
vec_b_1 = vcat(vec_b1_1, vec_b3_1)
vec_b_1 = hcat(vec_b_1)
#
vec_max_y_1 = repeat([1], numNodes^2)
vec_max_y_1 = hcat(vec_max_y_1)
vec_min_y_1 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_min_y_1[10 * (m - 1) + n] = mat_arcTwoNodes[m, n]
    end
end
#
vec_c_1 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_c_1[10 * (m - 1) + n] = mat_distance[m, n]
    end
end
vec_c_1 = hcat(vec_c_1)
#
vec_f_1 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_f_1[10 * (m - 1) + n] = mat_fixedCost[m, n]
    end
end
vec_f_1 = hcat(vec_f_1)
#
Benders.milp(n_x = numNodes^2,
             n_y = numNodes^2,
             vec_min_y = vec_min_y_1,
             vec_max_y = vec_max_y_1,
             vec_c = vec_c_1,
             vec_f = vec_f_1,
             vec_b = vec_b_1,
             mat_a = mat_a_1,
             mat_b = mat_b_1,
             epsilon = 0.0001,
             timesIterationMax = 1000)
# 1280 - sum(mat_arcTwoNodes .* mat_fixedCost)
# ----------------------------------------------------------------------------------------------------------------------
# Question 5
# Fomulation of Matrix A
mat_a1_2 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_a1_2[mn, mn] = -1
end
mat_a3_2 = zeros(numNodes, numNodes^2)
seq_zeroTenNinety = collect(0:10:90)
for m = 1: numNodes
    mat_a3_2[m, (10 * m - 9): (10 * m)] = repeat([-1], 10)  # sent out
    mat_a3_2[m, (seq_zeroTenNinety + repeat([m],10))] = repeat([1], 10)  # sent in
    mat_a3_2[:, (11 * m - 10)] = repeat([0], 10)  # self-sending doesn't count
end
mat_a_2 = vcat(mat_a1_2, mat_a3_2)
# Fomulation of Matrix B
mat_b1_2 = zeros(numNodes^2, numNodes^2)
for mn = 1: numNodes^2
    mat_b1_2[mn, mn] = 170
end
mat_b3_2 = zeros(numNodes, numNodes^2)
mat_b_2 = vcat(mat_b1_2, mat_b3_2)
# Fomulation of Vector b
vec_b1_2 = zeros(numNodes^2)
vec_b3_2 = zeros(numNodes)
for n = 1: numNodes
    vec_b3_2[n] = vec_netInject[n] - vec_netEject[n]
end
vec_b_2 = vcat(vec_b1_2, vec_b3_2)
vec_b_2 = hcat(vec_b_2)
#
vec_max_y_2 = repeat([1], numNodes^2)
vec_max_y_2 = hcat(vec_max_y_2)
#
vec_c_2 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_c_2[10 * (m - 1) + n] = mat_distance[m, n]
    end
end
vec_c_2 = hcat(vec_c_2)
#
vec_f_2 = zeros(numNodes^2)
for m = 1: numNodes
    for n = 1: numNodes
        vec_f_2[10 * (m - 1) + n] = mat_fixedCost[m, n]
    end
end
vec_f_2 = hcat(vec_f_2)
#
Benders.milp(n_x = numNodes^2,
            n_y = numNodes^2,
            vec_min_y = repeat([0], numNodes^2),
            vec_max_y = vec_max_y_2,
            vec_c = vec_c_2,
            vec_f = vec_f_2,
            vec_b = vec_b_2,
            mat_a = mat_a_2,
            mat_b = mat_b_2,
            epsilon = 0.0001,
            timesIterationMax = 1000)
# obj - sum(mat_arcTwoNodes .* mat_fixedCost)
# ----------------------------------------------------------------------------------------------------------------------
# using LightGraphs
# using GraphPlot
# using Compose
# using Cairo
# using Fontconfig
# mat_result_y = [0 1 0 0 0 0 0 0 0 0;
#                 0 0 0 0 0 1 0 0 0 0;
#                 0 0 0 0 0 0 0 0 0 1;
#                 1 0 0 0 1 0 0 0 0 0;
#                 0 0 0 0 0 1 0 0 0 0;
#                 1 1 0 1 1 0 0 1 0 0;
#                 0 1 1 0 0 0 0 0 0 0;
#                 0 0 0 0 0 0 0 0 1 0;
#                 0 0 0 0 1 0 0 0 0 0;
#                 0 0 1 0 0 1 1 1 1 0]
# graph_initial = DiGraph(mat_arcTwoNodes)
# plot_graph_initial = gplot(graph_initial, nodelabel = vec_nameNodes)
# graph_result = graph_initial
# mat_arcTwoNodes_new = mat_result_y - mat_arcTwoNodes
# for m = 1: numNodes
#     for n = 1: numNodes
#         if mat_arcTwoNodes_new[m,n] == 1
#             add_edge!(graph_result, m, n)
#         end
#     end
# end
# gplot(graph_result, nodelabel = vec_nameNodes)
# draw(PNG("graph_initial.png", 16cm, 16cm), gplot(graph_initial))
# draw(PNG("graph_result.png", 16cm, 16cm), gplot(result))

# Final Corrected
vec_x= [0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        10 10 0 10 40 0 0 20 0 0
        0 20 20 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 65 0 15 0]
vec_y= [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 1 0 1 0]

vec_y: [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 1 0 0 1 1 1 1 0]
result_q: 600.0

vec_x= [0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        10 30 0 10 40 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 15 0
        0 0 0 0 0 0 0 0 0 0
        0 0 20 0 0 0 25 35 0 0]

# vec_result_x = reshape(transpose(vec_x), numNodes^2)
vec_result_node = zeros(numNodes)
for m = 1: numNodes
    vec_result_node[m] = sum(vec_x[(10 * m - 9): (10 * m)]) - sum(vec_x[(seq_zeroTenNinety + repeat([m],10))])
end
vec_y= [0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 1
        1 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0
        1 1 0 1 1 0 0 1 0 0
        0 1 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0
        0 0 0 0 1 0 0 0 0 0
        0 0 1 0 0 1 1 1 1 0]

max_stars_repo_name: eitt/BendersDecomposition


text:
module DataReader
export Data
  mutable struct Data
    N_vertex::Int64
    N_edges::Int64
    initial_node::Int64
    end_node::Int64
    max_possible_flow::Array{Int64}
    start_nodes::Array{Int64}
  end

  N_vertex(d::Data) = d.N_vertex
  N_edges(d::Data) = d.N_edges
  initial_node(d::Data) = d.initial_node
  end_node(d::Data) = d.end_node
  max_possible_flow(d::Data) = d.max_possible_flow
  start_nodes(d::Data) = d.start_nodes

  d1 = Data(0,0,0,0,[],[]) # initializing all with 0

  function readData(filepath)
    open(filepath) do f
      d1.N_vertex = parse(Int64,readline(f))
      d1.N_edges = parse(Int64,readline(f))
      d1.initial_node = parse(Int64,readline(f))
      d1.end_node = parse(Int64,readline(f))

      d1.start_nodes = zeros(Int64, d1.N_vertex)

      global capacity = zeros(Int64,(d1.N_vertex, d1.N_vertex))

      vertices1 = Any[]
      vertices2 = Any[]

      for n in 1:d1.N_edges
        line = readline(f)
        vertex1 = parse(Int64,split(line, " ")[1])
        vertex2 = parse(Int64,split(line, " ")[2])
        cap = parse(Int64,split(line, " ")[3])

        capacity[vertex1, vertex2] = cap

        if ! (vertex1 in vertices1)
          push!(vertices1, vertex1)
        end

        if ! (vertex2 in vertices2)
          push!(vertices2, vertex2)
        end

      end

      for element in vertices2
        if element in vertices1
          filter!(e->e!=element,vertices1) #get all the initial_nodes
        end
      end

      d1.start_nodes = vertices1

      d1.max_possible_flow = zeros(Int64,length(d1.start_nodes))

      for i in d1.start_nodes
        for j in 1:d1.N_vertex
          d1.max_possible_flow[i] += capacity[i,j]
        end
        capacity[i,d1.end_node] = d1.max_possible_flow[i] #creates the new edges
      end
    end

    return capacity, d1
  end
end

max_stars_repo_name: ItamarRocha/Operations-Research

max_issues_repo_name: ItamarRocha/Operations-Research

max_forks_repo_name: ItamarRocha/Operations-Research


text:
module TestCompositionsBase
import Compat
using CompositionsBase
using CompositionsBase: @var_str
using Test

@testset "binary" begin
    @test tuple ∘ inv ===
          compose(tuple, inv) ===
          var"⨟"(inv, tuple) ===
          opcompose(inv, tuple)
end

@testset "unary" begin
    @test ∘(tuple) === compose(tuple) === var"⨟"(tuple) === opcompose(tuple) === tuple
end

if VERSION >= v"1.5.0-DEV.302"  # for ⨟
    include("test_julia15.jl")
end
if VERSION >= v"1.6"
    @testset "decompose, deopcompose" begin
        @test decompose(sin) === (sin,)
        @test decompose(sin∘cos) === (sin,cos)
        @test decompose(sin∘cos∘ tan) === (sin,cos,tan)
        @test decompose((sin∘cos)∘ tan) === (sin,cos,tan)
        @test decompose(sin∘(cos∘ tan)) === (sin,cos,tan)
        @test decompose((sqrt∘sin)∘(cos∘ tan)) === (sqrt,sin,cos,tan)

        @test deopcompose((sqrt∘sin)∘(cos∘ tan)) === (tan,cos,sin,sqrt)
        @inferred deopcompose((sqrt∘sin)∘(cos∘ tan))
        @inferred decompose((sqrt∘sin)∘(cos∘ tan))
    end
end

end  # module


text:
using SparseArrays, LinearAlgebra, Luxor, Colors, ProgressMeter

function NonLinSolver(Fm::Matrix{Float64}, Ns, L::Float64, E::Float64, nu::Float64, DT::Float64)
    L1, L2 = [L, 0.], [0., L]
    U, Fext, Fint = zeros(Float64, Ns[6], 1), zeros(Float64, Ns[6], 1), zeros(Float64, Ns[6], 1)
    gps = gaussPoints2D(2)
    function Assembly(U::Matrix{Float64})
        # Update the elements and assemble the system matrices.
        Fint, Krow, Kcol, K = zeros(Float64, Ns[6], 1), Int64[], Int64[], zeros(Float64, Ns[6], Ns[6])
        for i = 1:Ns[5]
            Ce = vec(DOFs[i, :, :]')
            elFint, elK = Quad4(U[Ce, :], i, elCoords[i, :, :], E, nu)
            Fint[Ce] += elFint
            K[Ce, Ce] += elK
        end
        K = sparse(K)
        Kims = K[2*Ns[1]+1:end, 2*Ns[1]+1:end]
        Rims = (Fint-Fext)[2*Ns[1]+1:end, :]
        return Fint, Kims, Rims
    end
    function generateEq(Kims::SparseMatrixCSC, Rims::Matrix{Float64})
        Kii = Kims[1:2*Ns[2], 1:2*Ns[2]]
        Kim = Kims[1:2*Ns[2], 2*Ns[2]+1:2*Ns[8]]
        Kis = Kims[1:2*Ns[2], 2*Ns[8]+1:2*Ns[10]]

        Kmi = Kims[2*Ns[2]+1:2*Ns[8], 1:2*Ns[2]]
        Kmm = Kims[2*Ns[2]+1:2*Ns[8], 2*Ns[2]+1:2*Ns[8]]
        Kms = Kims[2*Ns[2]+1:2*Ns[8], 2*Ns[8]+1:2*Ns[10]]

        Ksi = Kims[2*Ns[8]+1:2*Ns[10], 1:2*Ns[2]]
        Ksm = Kims[2*Ns[8]+1:2*Ns[10], 2*Ns[2]+1:2*Ns[8]]
        Kss = Kims[2*Ns[8]+1:2*Ns[10], 2*Ns[8]+1:2*Ns[10]]

        Ri = Rims[1:2*Ns[2], :]
        Rm = Rims[2*Ns[2]+1:2*Ns[8], :]
        Rs = Rims[2*Ns[8]+1:2*Ns[10], :]

        Kupper = hcat(Kii, Kim + Kis)
        Klower = hcat(Kmi + Ksi, Kmm + Ksm + Kms + Kss)

        K = vcat(Kupper, Klower)
        R = vcat(Ri, Rm + Rs)
        RNorm = norm(R)

        return K, R, RNorm
    end

    function Iterate(TOL::Float64, DT::Float64, Draw::Bool)
        t, tFinal = (0.0, 1.0)
        A(t) = (t/tFinal)*(Fm-I)[1:2, 1:2]

        p1 = Progress(Int64(tFinal/DT))
        while round(t, digits = 2) < tFinal
            t += DT
            Ua, Ub = A(t)*L1, A(t)*L2
            Ucorner = vcat(zeros(2, 1), Ua, Ua + Ub, Ub)
            Uinternal = U[2*Ns[1]+1:2*Ns[7], :]
            Umaster = U[2*Ns[7]+1:2*Ns[9], :]
            Uslave = Umaster + D * Ucorner
            U = vcat(Ucorner, Uinternal, Umaster, Uslave)
            Fint, Kims, Rims = Assembly(U)
            RNorm = Inf
            #@printf("-----t: %.2f-----\n",t)
            #p2 = ProgressThresh(TOL, "Minimizing: ")
            #update!(p2, RNorm)
            while RNorm > TOL
                K, R, RNorm = generateEq(Kims, Rims)
                dU = K \ -R
                Uinternal += dU[1:2*Ns[2], :]    # dUInternal
                Umaster += dU[2*Ns[2]+1:end, :]  # dUmaster
                Uslave = Umaster + D * Ucorner
                U = vcat(Ucorner, Uinternal, Umaster, Uslave)
                Fint, Kims, Rims = Assembly(U)
                #@printf("  | RNorm: %.3f \n", RNorm)
                #update!(p2, RNorm)
            end
            next!(p1)
        end
        if Draw
            deformed = elCoords + U[DOFs]
            maxUx, maxUy = max(maximum(elCoords[:,:,1]), maximum(deformed[:,:,1])), max(maximum(elCoords[:,:,2]), maximum(deformed[:,:,2]))
            scl = 200
            Drawing(Int(round(maxUx*scl*1.4)), Int(round(maxUy*scl*1.4)))
            background("white")
            origin(maxUx*scl*0.2, maxUy*scl*1.2)
            for i = 1:Ns[5]
                undefPoints = [Point(elCoords[i, j, 1]*scl, -elCoords[i, j, 2]*scl) for j in 1:4]
                defPoints = [Point(deformed[i, j, 1]*scl, -deformed[i, j, 2]*scl) for j in 1:4]
                sethue("black")
                poly(undefPoints, :stroke, close = true)
                sethue("blue")
                poly(defPoints, :stroke, close = true)
            end
            finish();
        end
        return U
    end
    return Iterate(1e-6, DT, true)
end

max_stars_repo_name: likemaestro/Nonlinear-Homogenization

max_issues_repo_name: likemaestro/Nonlinear-Homogenization

max_forks_repo_name: likemaestro/Nonlinear-Homogenization


text:
# exports
export create, Disk, show

mutable struct Disk <: Shape
    radius::Float64
    thickness::Float64
end

function Base.show(io::IO, x::Disk)
    println(io, "Disk")
    println(io, "Radius: $(x.radius)")
    println(io, "Thickness: $(x.thickness)")
end

function create(c::Disk; resolution=nothing, rand_=0.0, type::Int64=1)
    if isa(resolution, Nothing)
        resolution = min(c.radius/10, c.thickness/3)
    end
    radius=c.radius
    thickness=c.thickness
    e_size=resolution
    nm_radial_ele = round((radius-e_size/2)/e_size)
    er_size = (radius-e_size/2) / nm_radial_ele
    nm_thickness = round(thickness/e_size)
    et_size = (thickness/nm_thickness)
    mesh = Vector{Float64}[]
    vol = Float64[]
    total_vol = 0
    for i in 1:nm_thickness
        z = (i - 0.5)*et_size
        total_vol += et_size*pi*e_size^2
        push!(mesh, [0, 0.0, z + rand_*randn()*et_size])
        push!(vol, et_size*pi*e_size^2)
        for j in 1:nm_radial_ele
            r_in = e_size/2 + (j-1) * er_size
            r_out = e_size/2 + (j) * er_size

            nm_cir_ele = round(( pi * (r_in+r_out) ) / e_size )
            ec_size = ( pi * (r_in+r_out) ) / nm_cir_ele

            for k in 1:nm_cir_ele

                theta = 2*pi / nm_cir_ele
                e_area = pi* (r_out^2 - r_in^2) / nm_cir_ele
                r_g = (2 * pi * (r_out^3 - r_in^3) * sin(theta)/ 3 / theta ) / e_area / nm_cir_ele

                x = r_g * cos((k-1)*theta) + rand_*randn()*e_size
                y = r_g * sin((k-1)*theta) + rand_*randn()*e_size

                push!(vol, e_area * et_size)
                push!(mesh, [x, y, z + rand_*randn()*et_size])
                total_vol += e_area * et_size
            end
        end
    end
    mesh = hcat(mesh...)
    return mesh, zeros(size(mesh)), copy(mesh), vol, type*ones(Int64, length(vol))
end

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh


text:
formatters = Dict{ ASCIIStr, Function }()

cfmt( fmt::ASCIIStr, x ) = m_eval(Expr(:call, generate_formatter( fmt ), x))

function checkfmt(fmt)
    test = Base.Printf.parse( fmt )
    (length( test ) == 1 && typeof( test[1] ) <: Tuple) ||
        error( "Only one AND undecorated format string is allowed")
end

function generate_formatter( fmt::ASCIIStr )
    global formatters

    haskey( formatters, fmt ) && return formatters[fmt]

    if !occursin("'", fmt)
        checkfmt(fmt)
        return (formatters[ fmt ] = @eval(x->@sprintf( $fmt, x )))
    end

    conversion = fmt[end]
    conversion in "sduifF" ||
        error( string("thousand separator not defined for ", conversion, " conversion") )

    fmtactual = replace( fmt, "'" => ""; count=1 )
    checkfmt( fmtactual )
    conversion in "sfF" ||
        return (formatters[ fmt ] = @eval(x->checkcommas(@sprintf( $fmtactual, x ))))

    formatters[ fmt ] =
        if endswith( fmtactual, 's')
            @eval((x::Real)->((eltype(x) <: Rational)
                              ? addcommasrat(@sprintf( $fmtactual, x ))
                              : addcommasreal(@sprintf( $fmtactual, x ))))
        else
            @eval((x::Real)->addcommasreal(@sprintf( $fmtactual, x )))
        end
end

function addcommasreal(s)
    dpos = Compat.findfirst( isequal('.'), s )
    dpos !== nothing && return string(addcommas( s[1:dpos-1] ), s[ dpos:end ])
    # find the rightmost digit
    for i in length( s ):-1:1
        isdigit( s[i] ) && return string(addcommas( s[1:i] ), s[i+1:end])
    end
    s
end

function addcommasrat(s)
    # commas are added to only the numerator
    spos = Compat.findfirst( isequal('/'), s )
    string(addcommas( s[1:spos-1] ), s[spos:end])
end

function checkcommas(s)
    for i in length( s ):-1:1
        if isdigit( s[i] )
            s = string(addcommas( s[1:i] ), s[i+1:end])
            break
        end
    end
    s
end

function addcommas( s::ASCIIStr )
    len = length(s)
    t = ""
    for i in 1:3:len
        subs = s[max(1,len-i-1):len-i+1]
        if i == 1
            t = subs
        elseif match( r"[0-9]", subs ) != nothing
            t = string(subs, ',', t)
        else
            t = string(subs, t)
        end
    end
    t
end

function generate_format_string(;
                                width::Int=-1,
                                precision::Int= -1,
                                leftjustified::Bool=false,
                                zeropadding::Bool=false,
                                commas::Bool=false,
                                signed::Bool=false,
                                positivespace::Bool=false,
                                alternative::Bool=false,
                                conversion::ASCIIStr="f" #aAdecEfFiosxX
                                )
    s = ['%'%UInt8]
    commas &&
        push!(s, '\'')
    alternative && in( conversion[1], "aAeEfFoxX" ) &&
        push!(s, '#')
    zeropadding && !leftjustified && width != -1 &&
        push!(s, '0')
    if signed
        push!(s, '+')
    elseif positivespace
        push!(s, ' ')
    end
    if width != -1
        leftjustified && push!(s, '-')
        append!(s, _codeunits(string( width )))
    end
    precision != -1 &&
        append!(s, _codeunits(string( '.', precision )))
    String(append!(s, _codeunits(conversion)))
end

function format( x::T;
                 width::Int=-1,
                 precision::Int= -1,
                 leftjustified::Bool=false,
                 zeropadding::Bool=false, # when right-justified, use 0 instead of space to fill
                 commas::Bool=false,
                 signed::Bool=false, # +/- prefix
                 positivespace::Bool=false,
                 stripzeros::Bool=(precision== -1),
                 parens::Bool=false, # use (1.00) instead of -1.00. Used in finance
                 alternative::Bool=false, # usually for hex
                 mixedfraction::Bool=false,
                 mixedfractionsep::UTF8Str="_",
                 fractionsep::UTF8Str="/", # num / den
                 fractionwidth::Int = 0,
                 tryden::Int = 0, # if 2 or higher,
                                  # try to use this denominator, without losing precision
                 suffix::UTF8Str="", # useful for units/%
                 autoscale::Symbol=:none, # :metric, :binary or :finance
                 conversion::ASCIIStr=""
                 ) where {T<:Real}
    checkwidth = commas
    if conversion == ""
        if T <: AbstractFloat || T <: Rational && precision != -1
            actualconv = "f"
        elseif T <: Unsigned
            actualconv = "x"
        elseif T <: Integer
            actualconv = "d"
        else
            conversion = "s"
            actualconv = "s"
        end
    else
        actualconv = conversion
    end
    signed && commas && error( "You cannot use signed (+/-) AND commas at the same time")

    T <: Rational && conversion == "s" && (stripzeros = false)
    if ( T <: AbstractFloat && actualconv == "f" || T <: Integer ) && autoscale != :none
        actualconv = "f"
        if autoscale == :metric
            scales = [
                (1e24, "Y" ),
                (1e21, "Z" ),
                (1e18, "E" ),
                (1e15, "P" ),
                (1e12, "T" ),
                (1e9,  "G"),
                (1e6,  "M"),
                (1e3,  "k") ]
            if abs(x) > 1
                for (mag, sym) in scales
                    if abs(x) >= mag
                        x /= mag
                        suffix = string(sym, suffix)
                        break
                    end
                end
            elseif T <: AbstractFloat
                smallscales = [
                    ( 1e-12, "p" ),
                    ( 1e-9,  "n" ),
                    ( 1e-6,  "μ" ),
                    ( 1e-3,  "m" ) ]
                for (mag,sym) in smallscales
                    if abs(x) < mag*10
                        x /= mag
                        suffix = string(sym, suffix)
                        break
                    end
                end
            end
        else
            if autoscale == :binary
                scales = [
                    (1024.0 ^8,  "Yi" ),
                    (1024.0 ^7,  "Zi" ),
                    (1024.0 ^6,  "Ei" ),
                    (1024.0 ^5,  "Pi" ),
                    (1024.0 ^4,  "Ti" ),
                    (1024.0 ^3,  "Gi"),
                    (1024.0 ^2,  "Mi"),
                    (1024.0,     "Ki")
                ]
            else # :finance
                scales = [
                    (1e12, "t" ),
                    (1e9,  "b"),
                    (1e6,  "m"),
                    (1e3,  "k") ]
            end
            for (mag, sym) in scales
                if abs(x) >= mag
                    x /= mag
                    suffix = string(sym, suffix)
                    break
                end
            end
        end
    end

    nonneg = x >= 0
    fractional = 0
    if T <: Rational && mixedfraction
        actualconv = "d"
        actualx = trunc( Int, x )
        fractional = abs(x) - abs(actualx)
    else
        actualx = (parens && !in( actualconv[1], "xX" )) ? abs(x) : x
    end
    s = cfmt( generate_format_string( width=width,
                                      precision=precision,
                                      leftjustified=leftjustified,
                                      zeropadding=zeropadding,
                                      commas=commas,
                                      signed=signed,
                                      positivespace=positivespace,
                                      alternative=alternative,
                                      conversion=actualconv
                                      ),
              actualx)

    if T <:Rational && conversion == "s"
        if mixedfraction && fractional != 0
            num = fractional.num
            den = fractional.den
            if tryden >= 2 && mod( tryden, den ) == 0
                num *= div(tryden,den)
                den = tryden
            end
            fs = string( num, fractionsep, den)
            length(fs) < fractionwidth &&
                (fs = string(repeat( "0", fractionwidth - length(fs) ), fs))
            s = (actualx != 0
                 ? string(rstrip(s), mixedfractionsep, fs)
                 : (nonneg ? fs : string('-', fs)))
            checkwidth = true
        elseif !mixedfraction
            s = replace( s, "//" => fractionsep )
            checkwidth = true
        end
    elseif stripzeros && in( actualconv[1], "fFeEs" )
        dpos = Compat.findfirst( isequal('.'), s )
        dpos === nothing && (dpos = length(s))
        if actualconv[1] in "eEs"
            epos = Compat.findfirst(isequal(actualconv[1] == 'E' ? 'E' : 'e'), s)
            rpos = (epos === nothing) ? length( s ) : (epos-1)
        else
            rpos = length(s)
        end
        # rpos at this point is the rightmost possible char to start
        # stripping
        stripfrom = rpos+1
        for i = rpos:-1:dpos+1
            if s[i] == '0'
                stripfrom = i
            elseif s[i] ==' '
                continue
            else
                break
            end
        end
        if stripfrom <= rpos
            # everything after decimal is 0, so strip the decimal too
            s = string(s[1:stripfrom-1-(stripfrom == dpos+1)], s[rpos+1:end])
            checkwidth = true
        end
    end

    s = string(s, suffix)

    if parens && !in( actualconv[1], "xX" )
        # if zero or positive, we still need 1 white space on the right
        s = nonneg ? string(' ', strip(s), ' ') : string('(', strip(s), ')')
        checkwidth = true
    end

    if checkwidth && width != -1
        if (len = length(s) - width) > 0
            s = replace( s, " " => ""; count=len )
            if (len = length(s) - width) > 0
                endswith(s, " ") && (s = reverse(replace(reverse(s), " " => ""; count=len)))
                (len = length(s) - width) > 0 && (s = replace( s, "," => ""; count=len ))
            end
        elseif len < 0
            # Todo: should use lpad or rpad here, can be more efficient
            s = leftjustified ? string(s, repeat( " ", -len )) : string(repeat( " ", -len), s)
        end
    end

    s
end