text:
@testset ExtendedTestSet "euc" begin
    @testset ExtendedTestSet "euc" begin
        A = CartesianIndex(1, 3, 1)
        B = CartesianIndex(1, 3, 1)
        @test euc(A, B) == 0
    end

    @testset ExtendedTestSet "euc" begin
        A = CartesianIndex(1, 3, 1)
        C = CartesianIndex(1, 3, 2)
        @test euc(A, C) == 1.0
    end
end

@testset ExtendedTestSet "find_edges" begin
    @testset ExtendedTestSet "find_edges" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y2 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y2; dims=3)

        C = CartesianIndex{3}[
            CartesianIndex(1, 3, 1),
            CartesianIndex(2, 3, 1),
            CartesianIndex(3, 3, 1),
            CartesianIndex(4, 3, 1),
            CartesianIndex(1, 3, 2),
            CartesianIndex(2, 3, 2),
            CartesianIndex(3, 3, 2),
            CartesianIndex(4, 3, 2),
        ]
        @test find_edges(y) == C
    end
    @testset ExtendedTestSet "find_edges" begin
        x1 = [1 1 1 0; 1 1 1 0]
        x2 = [1 1 1 0; 1 1 1 0]
        x = cat(x1, x2; dims=3)

        C = CartesianIndex{3}[
            CartesianIndex(1, 3, 1)
            CartesianIndex(2, 3, 1)
            CartesianIndex(1, 3, 2)
            CartesianIndex(2, 3, 2)
        ]
        @test find_edges(x) == C
    end
end

@testset ExtendedTestSet "detect_edges_3D" begin
    @testset ExtendedTestSet "detect_edges_3D" begin
        x1 = [1 1 1 0; 1 1 1 0]
        x2 = [1 1 1 0; 1 1 1 0]
        x3 = [1 1 1 0; 1 1 1 0]
        x = cat(x1, x2, x3; dims=3)

        e1 = [0 0 1 0; 0 0 1 0]
        e2 = [0 0 1 0; 0 0 1 0]
        e3 = [0 0 1 0; 0 0 1 0]
        edge = cat(e1, e2, e3; dims=3)
        @test detect_edges_3D(x, Canny()) ≈ edge
    end
end


text:
function compute_residual(L, x, ∂, b, d::Domain)
    ϵ = new_rhs(d)
    mul!(ϵ, L, x)
    @. ϵ = ϵ + ∂ - b
    return ϵ
end

ext: jl


text:
module Newick

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos

import TranscodingStreams: TranscodingStreams, TranscodingStream

using AbstractTrees

mutable struct Record
    name::Union{Missing, String}
    distance::Union{Missing, Float64}
    parent::Union{Missing, Record}
    children::Union{Missing, Vector{Record}}
end

"Root constructor"
function Record()
    return Record(missing, missing, missing, missing)
end

"Child node constructor"
function Record(parent::Record)
    return Record(missing, missing, parent, missing)
end

function prenatal!(record)
    record.children = Vector{Record}()
    return record
end

function empty!(record::Record)
    record.name = missing
    record.distance = missing
    empty!(record.children)
    record.children = missing
    return record
end

function name(record::Record)
    return record.name
end

function name!(record::Record, name)
    record.name = name
    return record
end

function distance(record::Record)
    return record.distance
end

function distance!(record::Record, distance)
    record.distance = distance
    return record
end

function haschildren(record::Record)
    return !ismissing(record.children)
end

const machine = (function ()

    re = Automa.RegExp

    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    opt = Automa.RegExp.opt
    alt = Automa.RegExp.alt

    name = re"[0-9A-Za-z_-]+?"
    name.actions[:enter] = [:mark]
    name.actions[:exit] = [:name]

    distance = let
        len = cat(re"[0-9]+", opt(re"\.[0-9]*"))
        len.actions[:enter] = [:mark]

        opt(cat(re":", len))
    end
    distance.actions[:exit] = [:distance]

    nodefinish = re","
    nodefinish.actions[:enter] = [:nodefinish]

    cladestart = re"\("
    cladestart.actions[:enter] = [:cladestart]

    cladefinish = re"\)"
    cladefinish.actions[:enter] = [:nodefinish, :cladefinish]

    finish = re";"
    finish.actions[:enter] = [:finish]

    node = cat(
        opt(cladestart),
        opt(name),
        opt(distance),
        alt(nodefinish, cladestart, cladefinish, finish)
    )

    newick = rep(node)

    Automa.compile(newick)
end)()

const actions = Dict(
    :mark => :(@mark),

    :name => quote
        @debug "name $p"
        r = @markpos():p - 1
        if length(r) > 0
            str = String(data[r])
            name!(record, str)
        end
    end,

    :distance => quote
        @debug "distance $p"
        r = @markpos():p - 1
        if length(r) > 0
            str = String(data[r])
            parsed = parse(Float64, str)
            distance!(record, parsed)
        end
    end,

    :nodefinish => quote
        @debug "nodefinish $p"
        push!(record.parent.children, record)
        record = Record(record.parent)
    end,

    :cladestart => quote
        @debug "cladestart $p"
        prenatal!(record)
        record = Record(record) # Setup first child.
    end,

    :cladefinish => quote
        @debug "cladefinish $p"
        record = record.parent # Move to parent.
    end,

    :finish => quote
        @debug "finish $p"
    end
)

initcode = quote
    root = record
end

returncode = :(return cs, root)

context = Automa.CodeGenContext(generator = :goto, checkbounds = false, loopunroll = 8)

Automa.Stream.generate_reader(
    :read!,
    machine,
    arguments = (:(record::Record),),
    actions = actions,
    context = context,
    initcode = initcode,
    returncode = returncode
) |> eval

function Base.show(io::IO, node::Record)
    return print(io, (name=node.name, distance=node.distance, children=node.children))
end

function AbstractTrees.printnode(io::IO, node::Record)
    return print(io, (name=node.name, distance=node.distance))
end

function AbstractTrees.children(r::Record)
    if !haschildren(r)
        return ()
    end
    return r.children
end

end # module


text:
using ImagePhaseCongruency, Images, PyPlot

img = circsine(offset = pi/4, ampexponent = -1.5);

imshow(img)

fimg = imfilter(img, KernelFactors.gaussian((2, 2)))

(grad_y, grad_x) = imgradients(fimg, KernelFactors.ando3)

imshow(grad_y)

orient = orientation(grad_x, grad_y)
mag = magnitude(grad_x, grad_y)

thinned, subpix = thin_edges_nonmaxsup(mag, orient, radius=1.35, theta=pi/180)

imshow(thinned)

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
export spg_lasso

"""
SPG_LASSO  Solve the LASSO problem\n

   SPG_LASSO is designed to solve the LASSO problem\n

   (LASSO)  minimize  ||AX - B||_2  subject to  ||X||_1 <= tau,\n

   where A is an M-by-N matrix, B is an M-vector, and TAU is a
   nonnegative scalar.  In all cases below, A can be an explicit M-by-N
   matrix or matrix-like object for which the operations  A*x  and  A'*y
   are defined (i.e., matrix-vector multiplication with A and its
   adjoint.)\n

   Also, A can be a function handle that points to a function with the
   signature\n

   v = A(w,mode)   which returns:\n
                                  v = A *w  if mode == 1;\n
                                  v = A'*w  if mode == 2. \n

   X = SPG_LASSO(A,B,TAU) solves the LASSO problem.\n

   X = SPG_LASSO(A,B,TAU,OPTIONS) specifies options that are set using
   SPGSETPARMS.\n

   [X,R,G,INFO] = SPG_LASSO(A,B,TAU,OPTIONS) additionally returns the
   residual R = B - A*X, the objective gradient G = A'*R, and an INFO
   structure.  (See SPGL1 for a description of this last output argument.)\n

   See also spgl1, spgSetParms, spg_bp, spg_bpdn.\n

   Copyright 2008, Ewout van den Berg and Michael P. Friedlander\n
   http://www.cs.ubc.ca/labs/scl/spgl1\n
   Id: spg_lasso.m 1074 2008-08-19 05:24:28Z ewout78 \n
"""
function spg_lasso(A ,b; tau::AbstractFloat=NaN ,options::spgOptions=spgOptions())

    x, r, g, info = spgl1(A, b, tau = tau, options = options)

end


text:
export
ruleSPGaussianMeanPrecisionOutNPP,
ruleSPGaussianMeanPrecisionMPNP,
ruleSPGaussianMeanPrecisionOutNGP,
ruleSPGaussianMeanPrecisionMGNP,
ruleVBGaussianMeanPrecisionM,
ruleVBGaussianMeanPrecisionW,
ruleVBGaussianMeanPrecisionOut,
ruleSVBGaussianMeanPrecisionOutVGD,
ruleSVBGaussianMeanPrecisionW,
ruleSVBGaussianMeanPrecisionMGVD,
ruleMGaussianMeanPrecisionGGD,
ruleMGaussianMeanPrecisionGGN

ruleSPGaussianMeanPrecisionOutNPP(  msg_out::Nothing,
                                    msg_mean::Message{PointMass, V},
                                    msg_prec::Message{PointMass}) where V<:VariateType =
    Message(V, GaussianMeanPrecision, m=deepcopy(msg_mean.dist.params[:m]), w=deepcopy(msg_prec.dist.params[:m]))

ruleSPGaussianMeanPrecisionMPNP(msg_out::Message{PointMass}, msg_mean::Nothing, msg_prec::Message{PointMass}) =
    ruleSPGaussianMeanPrecisionOutNPP(msg_mean, msg_out, msg_prec)

function ruleSPGaussianMeanPrecisionOutNGP( msg_out::Nothing,
                                            msg_mean::Message{F, V},
                                            msg_prec::Message{PointMass}) where {F<:Gaussian, V<:VariateType}

    d_mean = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_mean.dist)

    Message(V, GaussianMeanVariance, m=d_mean.params[:m], v=d_mean.params[:v] + cholinv(msg_prec.dist.params[:m]))
end

ruleSPGaussianMeanPrecisionMGNP(msg_out::Message{F}, msg_mean::Nothing, msg_prec::Message{PointMass}) where F<:Gaussian =
    ruleSPGaussianMeanPrecisionOutNGP(msg_mean, msg_out, msg_prec)

ruleVBGaussianMeanPrecisionM(   dist_out::ProbabilityDistribution{V},
                                dist_mean::Any,
                                dist_prec::ProbabilityDistribution) where V<:VariateType =
    Message(V, GaussianMeanPrecision, m=unsafeMean(dist_out), w=unsafeMean(dist_prec))

function ruleVBGaussianMeanPrecisionW(  dist_out::ProbabilityDistribution{Univariate},
                                        dist_mean::ProbabilityDistribution{Univariate},
                                        dist_prec::Any)

    (m_mean, v_mean) = unsafeMeanCov(dist_mean)
    (m_out, v_out) = unsafeMeanCov(dist_out)

    Message(Univariate, Gamma, a=1.5, b=0.5*(v_mean + v_out + (m_mean - m_out)^2))
end

function ruleVBGaussianMeanPrecisionW(  dist_out::ProbabilityDistribution{Multivariate},
                                        dist_mean::ProbabilityDistribution{Multivariate},
                                        dist_prec::Any)

    (m_mean, v_mean) = unsafeMeanCov(dist_mean)
    (m_out, v_out) = unsafeMeanCov(dist_out)

    Message(MatrixVariate, Wishart, v=cholinv( v_mean + v_out + (m_mean - m_out)*(m_mean - m_out)' ), nu=dims(dist_out) + 2.0)
end

ruleVBGaussianMeanPrecisionOut( dist_out::Any,
                                dist_mean::ProbabilityDistribution{V},
                                dist_prec::ProbabilityDistribution) where V<:VariateType =
    Message(V, GaussianMeanPrecision, m=unsafeMean(dist_mean), w=unsafeMean(dist_prec))

ruleSVBGaussianMeanPrecisionOutVGD(dist_out::Any,
                                   msg_mean::Message{F, V},
                                   dist_prec::ProbabilityDistribution) where{F<:Gaussian, V<:VariateType} =
    Message(V, GaussianMeanVariance, m=unsafeMean(msg_mean.dist), v=unsafeCov(msg_mean.dist) + cholinv(unsafeMean(dist_prec)))

function ruleSVBGaussianMeanPrecisionW(
    dist_out_mean::ProbabilityDistribution{Multivariate, F},
    dist_prec::Any) where F<:Gaussian

    joint_dims = dims(dist_out_mean)
    d_out_mean = convert(ProbabilityDistribution{Multivariate, GaussianMeanVariance}, dist_out_mean)
    (m, V) = unsafeMeanCov(d_out_mean)
    if joint_dims == 2
        return Message(Univariate, Gamma, a=1.5, b=0.5*(V[1,1] - V[1,2] - V[2,1] + V[2,2] + (m[1] - m[2])^2))
    else
        d = Int64(joint_dims/2)
        return Message(MatrixVariate, Wishart, v=cholinv( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end, 1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ), nu=d + 2.0)
    end
end

function ruleSVBGaussianMeanPrecisionMGVD(  msg_out::Message{F, V},
                                            dist_mean::Any,
                                            dist_prec::ProbabilityDistribution) where {F<:Gaussian, V<:VariateType}

    d_out = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_out.dist)

    Message(V, GaussianMeanVariance, m=d_out.params[:m], v=d_out.params[:v] + cholinv(unsafeMean(dist_prec)))
end

function ruleMGaussianMeanPrecisionGGD(
    msg_out::Message{<:Gaussian, V},
    msg_mean::Message{<:Gaussian, V},
    dist_prec::ProbabilityDistribution) where V<:VariateType

    d_mean = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_mean.dist)
    d_out = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_out.dist)

    xi_y = d_out.params[:xi]
    W_y = d_out.params[:w]
    xi_m = d_mean.params[:xi]
    W_m = d_mean.params[:w]
    W_bar = unsafeMean(dist_prec)

    return ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[xi_y; xi_m], w=[W_y+W_bar -W_bar; -W_bar W_m+W_bar])
end

function ruleMGaussianMeanPrecisionGGN(
    msg_out::Message{<:Gaussian, V},
    msg_mean::Message{<:Gaussian, V},
    msg_prec::Message{PointMass}) where V<:VariateType

    d_mean = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_mean.dist)
    d_out = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_out.dist)

    xi_y = d_out.params[:xi]
    W_y = d_out.params[:w]
    xi_m = d_mean.params[:xi]
    W_m = d_mean.params[:w]
    W_bar = msg_prec.dist.params[:m]

    return ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[xi_y; xi_m], w=[W_y+W_bar -W_bar; -W_bar W_m+W_bar])
end


text:
_throw_table_error() = throw(ArgumentError("Please specify the column that contains the targets explicitly, or provide a target-extraction-function as first parameter. see parameter 'f' in ?targets."))

# required data container interface
LearnBase.nobs(dt::AbstractDataFrame) = DataFrames.nrow(dt)
LearnBase.getobs(dt::AbstractDataFrame, idx) = dt[idx,:]

LearnBase.nobs(dt::DataFrameRow) = 1  # it is a observation
function LearnBase.getobs(dt::DataFrameRow, idx)
    idx == 1:1 || throw(ArgumentError(
         "Attempting to read multiple rows ($idx) with a single row"))

    return dt
end

# custom data subset in form of SubDataFrame
LearnBase.datasubset(dt::AbstractDataFrame, idx, ::ObsDim.Undefined) =
    @view dt[idx, :]

# throw error if no target extraction function is supplied
LearnBase.gettarget(::typeof(identity), dt::AbstractDataFrame) =
    _throw_table_error()
LearnBase.gettarget(::typeof(identity), dt::DataFrameRow) =
    _throw_table_error()

# convenience syntax to allow column name
LearnBase.gettarget(col::Symbol, dt::AbstractDataFrame) = dt[1, col]
LearnBase.gettarget(col::Symbol, dt::DataFrameRow) = dt[col]
LearnBase.gettarget(fun, dt::AbstractDataFrame) = fun(dt)

# avoid copy when target extraction function is supplied
MLDataPattern.getobs_targetfun(dt::AbstractDataFrame) = dt

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

#= Use Requires.jl once fixed for 0.6

import DataTables: DataTables, AbstractDataTable

# required data container interface
LearnBase.nobs(dt::AbstractDataTable) = DataTables.nrow(dt)
LearnBase.getobs(dt::AbstractDataTable, idx) = dt[idx,:]

# custom data subset in form of SubDataFrame
LearnBase.datasubset(dt::AbstractDataTable, idx, ::ObsDim.Undefined) =
    view(dt, idx)

# throw error if no target extraction function is supplied
LearnBase.gettarget(dt::AbstractDataTable) = _throw_table_error()

# convenience syntax to allow column name
LearnBase.gettarget(col::Symbol, dt::AbstractDataTable) = dt[1, col]

# avoid copy when target extraction function is supplied
MLDataPattern.getobs_targetfun(dt::AbstractDataTable) = dt

=#


text:
module TestDoctest

import SyncBarriers
using Documenter: doctest
using Test

function test()
    doctest(SyncBarriers)
end

end  # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Ghostscript_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Ghostscript")
JLLWrappers.@generate_main_file("Ghostscript", UUID("61579ee1-b43e-5ca0-a5da-69d92c66a64b"))
end  # module Ghostscript_jll


text:
#!/usr/bin/env julia

using Pkg3.TOML
using Base: LibGit2

function sha1map(pkgs::Dict{String,Package})
    f = joinpath(@__DIR__, "sha1map.toml")
    d = ispath(f) ? TOML.parsefile(f) : Dict()
    for (pkg, p) in pkgs
        isempty(p.versions) && continue
        uuid = string(p.uuid)
        haskey(d, uuid) || (d[uuid] = Dict())
        updated = false
        repo = nothing
        for (ver, v) in p.versions
            haskey(d[uuid], v.sha1) && continue
            git_commit_hash = LibGit2.GitHash(v.sha1)
            if repo == nothing
                repo_path = joinpath(homedir(), ".julia", "upstream", uuid)
                repo = ispath(repo_path) ? LibGit2.GitRepo(repo_path) : begin
                    updated = true
                    info("Cloning [$uuid] $pkg")
                    LibGit2.clone(p.url, repo_path, isbare=true)
                end
            end
            if !updated
                try LibGit2.GitObject(repo, git_commit_hash)
                catch err
                    err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
                    info("Updating $pkg from $(p.url)")
                    LibGit2.fetch(repo, remoteurl=p.url, refspecs=["+refs/*:refs/remotes/cache/*"])
                end
            end
            git_commit = try LibGit2.GitObject(repo, git_commit_hash)
            catch err
                err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
                error("$pkg: git object $(v.sha1) could not be found")
            end
            git_commit isa LibGit2.GitCommit || git_commit isa LibGit2.GitTag ||
                error("$pkg: git object $(v.sha1) not a commit – $(typeof(git_commit))")
            git_tree = LibGit2.peel(LibGit2.GitTree, git_commit)
            @assert git_tree isa LibGit2.GitTree
            git_tree_hash = string(LibGit2.GitHash(git_tree))
            d[uuid][v.sha1] = git_tree_hash
        end
    end
    open(f, "w") do io
        TOML.print(io, d, sorted=true)
    end
    return d
end


text:
function dfs(maze::Maze, cur_pos::Tuple{Int, Int})::Array{Tuple{Int, Int}, 1}
	if cur_pos == maze.exit
		set_mazecolor(maze, cur_pos, 0xffff0000)
		return [cur_pos]
	end
	set_mazecolor(maze, cur_pos, 0xffafafaf)
	dirs = get_dirs(maze, cur_pos)

	ret = [cur_pos]
	if dirs & 1 != 0
		chk_pos = (cur_pos[1] + 1, cur_pos[2])
		temp = dfs(maze, chk_pos)
		if length(temp) > 0
			set_mazecolor(maze, cur_pos, 0xffff9400)
			append!(ret, temp)
			return ret
		end
	end
	if dirs & 2 != 0
		chk_pos = (cur_pos[1], cur_pos[2] + 1)
		temp = dfs(maze, chk_pos)
		if length(temp) > 0
			set_mazecolor(maze, cur_pos, 0xffff9400)
			append!(ret, temp)
			return ret
		end
	end
	if dirs & 4 != 0
		chk_pos = (cur_pos[1] - 1, cur_pos[2])
		temp = dfs(maze, chk_pos)
		if length(temp) > 0
			set_mazecolor(maze, cur_pos, 0xffff9400)
			append!(ret, temp)
			return ret
		end
	end
	if dirs & 8 != 0
		chk_pos = (cur_pos[1], cur_pos[2] - 1)
		temp = dfs(maze, chk_pos)
		if length(temp) > 0
			set_mazecolor(maze, cur_pos, 0xffff9400)
			append!(ret, temp)
			return ret
		end
	end
	return []
end

function dfs(maze::Maze)::Array{Tuple{Int, Int}, 1}
	sx = maze.entr[1]
	sy = maze.entr[2]

	maze.data[sy, sx] = 0xff00ff00

	if sx == 1
		sx += 1
	elseif sx == maze.w
		sx -= 1
	elseif sy == 1
		sy += 1
	else
		sy -= 1
	end

	return dfs(maze, (sx, sy))
end


text:
def TEXTWITHATTRIBUTES (TEXTALIGNMENT='centre', TEXTANGLE=0, TEXTWIDTH=1.0, TEXTHEIGHT=2.0, TEXTSPACING=0.25):
    PRINT(TEXTALIGNMENT)
    ALIGN = IF([ K(TEXTALIGNMENT == 'centre'),
				 COMP([ APPLY, CONS([ COMP([ T(1), RAISE(DIV), CONS([ SIZE(1), K(-2) ]) ]), ID ]) ]),
				 IF([ K(TEXTALIGNMENT == 'right'),
					  COMP([ APPLY, CONS([ COMP([ T(1), RAISE(DIFF), SIZE(1) ]), ID ]) ]),
					  ID ]) ])
    HANDLE = CONS([
				COMP([ AA(S([1,2])([TEXTWIDTH/FONTWIDTH, TEXTHEIGHT/FONTHEIGHT])), charpols ]),
				K(T(1)(TEXTSPACING + TEXTWIDTH)) ])
    fn = COMP([ R([1,2])(TEXTANGLE), ALIGN, STRUCT, CAT, DISTR, HANDLE, CHARSEQ ])
    return fn

	function textWithAttributes(textalignment="centre", textangle=0, textwidth=1.0, textheight=2.0, textspacing=0.25)
		id = x->x
		align = if textalignment == 'centre'
					COMP([ APPLY, CONS([ COMP([ T(1), RAISE(DIV), CONS([ SIZE(1), K(-2) ]) ]), ID ]) ]),
				elseif textalignment == 'right'
					COMP([ APPLY, CONS([ COMP([ T(1), RAISE(DIFF), SIZE(1) ]), id ]) ]),
				else
					id
				end
		handle = CONS([
					COMP([ AA( Lar.s(GL.textwidth/GL.fontwidth, GL.textheight/GL.fontheight) ), GL.charpols ]),
					GL.k(Lar.t(textwidth+textspacing,0)) ])

		fn = COMP([ R([1,2])(TEXTANGLE), ALIGN, STRUCT, CAT, DISTR, HANDLE, GL.charseq ])
		return fn
	end

		function textWithAttributes0(strand)
			mat = Lar.s(GL.textwidth/GL.fontwidth,GL.textheight/GL.fontheight)
			Lar.apply(Lar.r(textangle),
			COMP([ align(textalignment),
			   Lar.struct2lar,
			   Lar.Struct,
			   GL.cat,
			   GL.distr,
			   GL.cons([ GL.a2a(mat) ∘ GL.charpols,
					GL.k(Lar.t(textwidth+textspacing,0)) ]),
			   GL.charseq ]))(strand)
		end
		return textWithAttributes0
	end

	if __name__=="__main__":
	    VIEW(TEXT('PLASM'))
	    VIEW(TEXTWITHATTRIBUTES('centre',PI/4,0.1,0.4,0.05)('PLASM'))

"""
	COMP(Funs)

Compose an Array of functions.
```
julia> COMP([sin,cos,tan])(pi/4)
0.5143952585235492

julia> COMP([sin,cos,tan])(pi/4) == (sin∘cos)(1)
true
```
"""
function COMP(Funs)
	function ∘(f,g)
		function h(x)
			return f(g(x))
		end
		return h
	end
	return reduce(∘,Funs)
end

"""
	lift(f::Function)::Function
```
```
"""
function lift(f::Function)::Function
	return funs-> COMP([f, CONS(funs)])
end

"""
	raise(f::Function)::Function
```
```
"""
function raise(f::Function)::Function
	function raise0(args)
		return (issequence(isfun) ? lift(f) : f)(args)
	end
	return raise0
end

"""
	julia> isfun(x::Any)::Bool
# Example
```
julia> isfun(sin)
true

julia> isfun(x->2x)
true

julia> isfun("a string")
false
```
"""
isfun(x::Any)::Bool = isa(x, Function)

"""
```
julia> args = [sin,cos,tan]
3-element Array{Function,1}:
 sin
 cos
 tan

julia> isSequenceOf(isfun)(args)
true
```
"""
function isSequenceOf(predicate)
	function isSequenceof0(args)
		return (&)(map(predicate,args)...)
	end
	return isSequenceof0
end

"""
	ispol(x::Any)::Bool
```
julia> ispol(Lar.cuboid([1,1,1]))
true
```
"""
ispol(x::Any)::Bool = isa(x,Lar.LAR)

"""
	isnum(x::Any)::Bool
```
julia> isnum("234")
false

julia> isnum(234)
true
```
"""
isnum(x::Any)::Bool = isa(x,Number)

"""
	isstring(x::Any)::Bool
```
julia> isstring("234")
true

julia> isstring(234)
false
```
"""
isstring(x::Any)::Bool = isa(x,String)

function diff(args):
	if isinstance(args,list) and ISPOL(args[0]):
		return DIFFERENCE(args)
	elseif isnum(args):
		return -1 * args
	elseif isinstance(args,list) and  ISNUM(args[0]):
		return reduce(lambda x,y: x - y, args)
	elseif isinstance(args,list) and  isinstance(args[0],list):
		#matrix difference
		if isinstance(args[0][0],list):
			return AA(VECTDIFF)(zip(*args))
		# vector diff
		else:
			return VECTDIFF(args)
		end
	end

	raise Exception("\'-\' function has been applied to %s!" % repr(args))

end

if __name__ == "__main__":
	assert(DIFF(2)==-2 and DIFF([1,2,3])==-4 and DIFF([[1,2,3],[1,2,3]])==[0,0,0])


text:
analyze!(g::EGraph, an::Type{<:AbstractAnalysis}, id::EClassId) =
    analyze!(g, an, reachable(g, id))

function analyze!(g::EGraph, an::Type{<:AbstractAnalysis})
    analyze!(g, an, collect(keys(g.classes)))
end

"""

**WARNING**. This function is unstable.
An [`EGraph`](@ref) can only contain one analysis of type `an`.
"""
function analyze!(g::EGraph, an::Type{<:AbstractAnalysis}, ids::Vector{EClassId})
    push!(g.analyses, an)
    ids = sort(ids)
    # @assert isempty(g.dirty)

    did_something = true
    while did_something
        did_something = false

        for id ∈ ids
            eclass = geteclass(g, id)
            id = eclass.id
            pass = mapreduce(x -> make(an, g, x), (x, y) -> join(an, x, y), eclass)
            # pass = make_pass(G, analysis, find(G,id))

            # if pass !== missing
            if !isequal(pass, getdata(eclass, an, missing))
                setdata!(eclass, an, pass)
                did_something = true
                push!(g.dirty, id)
            end
        end
    end

    for id ∈ ids
        eclass = geteclass(g, id)
        id = eclass.id
        if !hasdata(eclass, an)
            # display(g.classes[id]); println()
            # display(analysis.data); println()
            error("failed to compute analysis for eclass ", id)
        end
    end

    # rebuild!(g)

    # display(analysis.data); println()

    return true
end


text:
module Stripeline

include("instrumentdb.jl")
include("scanning.jl")
include("timesplit.jl")
include("mapmaker.jl")
include("tod_splitter.jl")
include("DBforsimulation.jl")
include("noisegeneration.jl")
include("adc.jl")

end


text:
module AtomEnergyLevels

import LinearAlgebra: Symmetric, eigen, diagind, I
import Printf: @sprintf

# rectangle rule for integration is accurate enough
@inline ∫(h, f) = h*sum(f)

include("sincdif.jl")
include("conf_parse.jl")
include("periodic_table.jl")
include("xclib.jl")
include("xc_func.jl")
include("rshreq.jl")
include("TF.jl")
include("lda.jl")

export laplacian, radial_shr_eq, TF, lda
export LDA_X, LDA_C_CHACHIYO, LDA_C_VWN
export XC_Type, XC_Functional, xc_func
export xc_type_lda, xc_type_gga, xc_type_mgga, xc_type_hybrid
export conf_enc, Nₑ, @c_str
export shells, atom
export laplacian
export radial_shr_eq

end # module


text:
module Layout

import Parameters: @with_kw

import Term:
    int,
    get_lr_widths,
    textlen,
    cint,
    fint,
    rtrim_str,
    ltrim_str,
    calc_nrows_ncols,
    do_by_line
import ..Renderables: RenderablesUnion, Renderable, AbstractRenderable, RenderableText
import ..Consoles: console_width, console_height
import ..Boxes: get_lrow, get_rrow
import ..Style: apply_style
import ..Measures: Measure, height, width
import ..Segments: Segment
using ..Boxes

export Padding, vstack, hstack, pad, pad!, vertical_pad, vertical_pad!
export Spacer, vLine, hLine, PlaceHolder
export leftalign!, center!, rightalign!
export leftalign, center, rightalign
export lvstack, cvstack, rvstack
export grid

# ---------------------------------------------------------------------------- #
#                                    PADDING                                   #
# ---------------------------------------------------------------------------- #
"""
    Padding

Stores information about ammount of padding.
"""
@with_kw struct Padding
    left::Int
    right::Int
    top::Int
    bottom::Int
end

"""
    pad(text::AbstractString, target_width::Int, method::Symbol)::String

Pad a string to width: `target_width` by adding empty spaces strings " ".
Where the spaces are added depends on the justification `method` ∈ (`:left`, `:center`, `:right`).

#### Examples

``` julia
julia> pad("ciao", 10, :left)
"ciao      "

julia> pad("ciao", 10, :center)
"   ciao   "

julia> pad("ciao", 10, :right)
"      ciao"
```
"""
function pad(text::AbstractString, target_width::Int, method::Symbol; bg = nothing)::String
    if occursin('\n', text)
        return do_by_line(ln -> pad(ln, target_width, method), text)
    end

    # get total padding size
    lw = width(text)
    lw ≥ target_width && return text

    npads = target_width - lw
    if method == :left
        p = isnothing(bg) ? " "^npads : "{$bg}" * " "^npads * "{/$bg}"
        return text * p
    elseif method == :right
        p = isnothing(bg) ? " "^npads : "{$bg}" * " "^npads * "{/$bg}"
        return p * text
    else
        nl, nr = get_lr_widths(npads)
        l = isnothing(bg) ? " "^nl : "{$bg}" * " "^nl * "{/$bg}"
        r = isnothing(bg) ? " "^nr : "{$bg}" * " "^nr * "{/$bg}"
        return l * text * r
    end
end

"""
---
    pad(text::AbstractString, left::Int = 0, right::Int = 0)

Pad a string by a fixed ammount to the left and to the right.
"""
function pad(text::AbstractString, left::Int = 0, right::Int = 0; bg = nothing)
    if isnothing(bg)
        return " "^max(0, left) * text * " "^max(0, right)
    else
        l = "{$bg}" * " "^max(0, left) * "{/$bg}"
        r = "{$bg}" * " "^max(0, right) * "{/$bg}"
        return l * text * r
    end
end
"""
---
    pad(segments::AbstractVector{Segment}, left::Int = 0, right::Int = 0)

Pad a renderable's segments to the left and the right.
"""
function pad(segments::AbstractVector{Segment}, left::Int = 0, right::Int = 0)
    return map((s) -> Segment(pad(s.text, left, right)), segments)
end

"""
---
    pad(ren::AbstractRenderable, left::Int, right::Int)

Pad an `AbstractRenderable` by padding each of its segments.
"""
function pad(ren::AbstractRenderable, left::Int, right::Int)
    segments = pad(ren.segments, left, right)
    measure = Measure(segments)
    return Renderable(segments, measure)
end

"""
---
    pad(ren::AbstractRenderable; width::Int)

Pad a renderable to achieve a target width.

!!! note
    The padding is added to the left and to the right as needed to achieve the target width.
    The resulting renderable object will be **center** in the target space.

#### Example
```
julia> pad(RenderableText("ciao"); width=10)
    ciao
```
"""
function pad(ren::AbstractRenderable; width::Int)
    npads = width - ren.measure.w
    nl, nr = get_lr_widths(npads)
    return pad(ren, nl, nr)
end

"""
    pad!(ren::AbstractRenderable, left::Int, right::Int)

In place version for padding a renderable.
"""
function pad!(ren::AbstractRenderable, left::Int, right::Int)
    ren.segments = pad(ren.segments, left, right)
    return ren.measure = Measure(ren.segments)
end

"""
    pad!(ren::AbstractRenderable; width::Int)

In place version for padding a renderable to achieve a given width.
"""
function pad!(ren::AbstractRenderable; width::Int)
    npads = width - ren.measure.w
    nl, nr = get_lr_widths(npads)
    return pad!(ren, nl, nr)
end

# ------------------------------- vertical pad ------------------------------- #

"""
vertical_pad(text, target_height::Int, method::Symbol)::String

Vertically pad a string to height: `target_height` by adding empty strings above/below " ".
Where the spaces are added depends on the justification `method` ∈ (:top, :center, :bottom).
"""
function vertical_pad(text::AbstractString, target_height::Int, method::Symbol)::String
    # get total padding size
    h = height(text)
    h ≥ target_height && return text

    space = " "^(width(text))
    npads = target_height - h
    if method == :bottom
        return vertical_pad(text, npads, 0)
    elseif method == :top
        return vertical_pad(text, 0, npads)
    else
        above, below = get_lr_widths(npads)
        return vertical_pad(text, above, below)
    end
end

"""
    vertical_pad(text::AbstractString, above::Int = 0, below::Int = 0)

Vertical pad a string by a fixed ammount to above and below.
"""
function vertical_pad(text::AbstractString, above::Int = 0, below::Int = 0)
    space = " "^(width(text))
    return string(vstack(repeat([space], above)..., text, repeat([space], below)...))
end

"""
    vertical_pad(segments::AbstractVector{Segment}, above::Int = 0, below::Int = 0)

Pad a renderable's segments to the above and the below.
"""
function vertical_pad(segments::AbstractVector{Segment}, above::Int = 0, below::Int = 0)
    space = " "^(segments[1].measure.w)
    above::Vector{Segment} = repeat([Segment(space)], above)
    below::Vector{Segment} = repeat([Segment(space)], below)
    return [above..., segments..., below...]
end

"""
    vertical_pad(ren::AbstractRenderable, above::Int, below::Int)

Pad a renderable, vertically.
"""
function vertical_pad(ren::AbstractRenderable, above::Int, below::Int)
    segments = vertical_pad(ren.segments, above, below)
    measure = Measure(segments)
    return Renderable(segments, measure)
end

"""
vertical_pad(ren::AbstractRenderable; height::Int)

Vertical pad a renderable to achieve a target height.
"""
function vertical_pad(ren::AbstractRenderable; height::Int)
    npads = height - ren.measure.h
    nl, nr = get_lr_widths(npads)
    return vertical_pad(ren, nl, nr)
end

"""
    verti0cal_pad!(ren::AbstractRenderable, above::Int, below::Int)

In place version for vertically padding a renderable.
"""
function vertical_pad!(ren::AbstractRenderable, above::Int, below::Int)
    ren.segments = vertical_pad(ren.segments, above, below)
    return ren.measure = Measure(ren.segments)
end

"""
    pad!(ren::AbstractRenderable; width::Int)

In place version for padding a renderable to achieve a given width.
"""
function vertical_pad!(ren::AbstractRenderable; height::Int)
    npads = height - ren.measure.h
    nl, nr = get_lr_widths(npads)
    return vertical_pad!(ren, nl, nr)
end

# ---------------------------------------------------------------------------- #
#                                    JUSTIFY                                   #
# ---------------------------------------------------------------------------- #
"""
    leftalign(renderables::RenderablesUnion...)

Pad two (or more) renderables so that they have the same width and they
are left-aligned.

NOTE: the renderables returned  have different type and potentially different size
    from the inputs

# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
p1, p2 = leftalign(p1, p2)
print(p1/p2)

╭───────────────────────╮
╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
```
"""
function leftalign(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)
    renderables = map(r -> pad(r, 0, width - r.measure.w), renderables)
    return renderables
end

"""
    leftalign!(renderables::RenderablesUnion...)

In place version of leftalign.

# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
leftalign!(p1, p2)
print(p1/p2)

╭───────────────────────╮
╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
"""
function leftalign!(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)
    for ren in renderables
        pad!(ren, 0, width - ren.measure.w)
    end
end

"""
    center(renderables::RenderablesUnion... )

Pad two (or more) renderables so that they have the same width and they
are centered.

NOTE: the renderables returned  have different type and potentially different size
    from the inputs

# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
p1, p2 = center(p1, p2)
print(p1/p2)

             ╭───────────────────────╮
             ╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
```
"""
function center(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)
    renderables = map(r -> pad(r; width = width), renderables)
    return renderables
end

"""
    center!(renderables::RenderablesUnion... )

In place version of `center`.

# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
center!(p1, p2)
print(p1/p2)

             ╭───────────────────────╮
             ╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
```
"""
function center!(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)
    for ren in renderables
        pad!(ren; width = width)
    end
end

"""
    rightalign(renderables::RenderablesUnion... )

Pad two (or more) renderables so that they have the same width and they
are right aligned.

NOTE: the renderables returned  have different type and potentially different size
    from the inputs

# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
p1, p2 = rightalign(p1, p2)
print(p1/p2)

                         ╭───────────────────────╮
                         ╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
```
"""
function rightalign(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)
    renderables = map(r -> pad(r, width - r.measure.w, 0), renderables)
    return renderables
end

"""
    rightalign!(renderables::RenderablesUnion... )

In place version of `rightalign`.
# Examples
```julia
p1 = Panel(; width=25)
p2 = Panel(; width=50)
rightalign!(p1, p2)
print(p1/p2)

                         ╭───────────────────────╮
                         ╰───────────────────────╯
╭────────────────────────────────────────────────╮
╰────────────────────────────────────────────────╯
```
"""
function rightalign!(renderables::RenderablesUnion...)
    length(renderables) < 2 && return renderables
    renderables = Renderable.(renderables)
    width = max(map(r -> r.measure.w, renderables)...)

    for ren in renderables
        pad!(ren, width - ren.measure.w, 0)
    end
end

# ---------------------------------------------------------------------------- #
#                                   STACKING                                   #
# ---------------------------------------------------------------------------- #

vstack(s1::String, s2::String; pad::Int = 0) = s1 * "\n"^(pad + 1) * s2

"""
    vstack(renderables...)

Vertically stack a variable number of renderables to give a new renderable
"""
function vstack(renderables...; pad::Int = 0)
    renderables = leftalign(renderables...)

    segments::Vector{Segment} = []
    if pad > 0
        renderables = foldl((a, b) -> a / ("\n"^pad) / b, renderables)
        segments = renderables.segments
    else
        segments = vcat(getfield.(renderables, :segments)...)
    end
    measure = Measure(segments)
    return Renderable(segments, measure)
end

vstack(renderables::Union{Vector,Tuple}; kwargs...) = vstack(renderables...; kwargs...)

"""
    hstack(r1::RenderablesUnion, r2::RenderablesUnion)

Horizontally stack two renderables to give a new renderable.
"""
function hstack(r1::RenderablesUnion, r2::RenderablesUnion; pad::Int = 0)
    r1 = r1 isa AbstractRenderable ? r1 : Renderable(r1)
    r2 = r2 isa AbstractRenderable ? r2 : Renderable(r2)

    # get dimensions of final renderable
    h1 = r1.measure.h
    h2 = r2.measure.h
    Δh = abs(h2 - h1)

    # make sure both renderables have the same number of segments
    if h1 > h2
        s1 = r1.segments
        s2 = vcat(r2.segments, [Segment(" "^(r2.measure.w)) for i in 1:(Δh)])
    elseif h1 < h2
        s1 = vcat(r1.segments, [Segment(" "^(r1.measure.w)) for i in 1:(Δh)])
        s2 = r2.segments
    else
        s1, s2, = r1.segments, r2.segments
    end

    # combine segments
    segments = [Segment(s1.text * " "^pad * s2.text) for (s1, s2) in zip(s1, s2)]

    return Renderable(segments, Measure(segments))
end

"""
    hstack(renderables...)

Horizonatlly stack a variable number of renderables.
"""
function hstack(renderables...; pad::Int = 0)
    renderable = Renderable()

    for (i, ren) in enumerate(renderables)
        _pad = i == 1 ? 0 : pad
        renderable = hstack(renderable, ren; pad = _pad)
    end
    return renderable
end

hstack(renderables::Union{Vector,Tuple}; kwargs...) = hstack(renderables...; kwargs...)

# --------------------------------- operators -------------------------------- #

"""
    Operators for more coincise layout of renderables
"""

Base.:/(r1::RenderablesUnion, r2::RenderablesUnion) = vstack(r1, r2)
Base.:/(rr::Tuple{RenderablesUnion,RenderablesUnion}) = vstack(rr...)

Base.:*(r1::AbstractRenderable, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractString, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractRenderable, r2::AbstractString) = hstack(r1, r2)

# --------------------------- convenience functions -------------------------- #

"""
    lvstack(renderables::RenderablesUnion...)

Left align renderables and then vertically stack.
"""
function lvstack(renderables::RenderablesUnion...; kwargs...)::Renderable
    renderables = leftalign(renderables...)
    return vstack(renderables...; kwargs...)
end

"""
    lvstack(renderables::RenderablesUnion...)

Center align renderables and then vertically stack.
"""
function cvstack(renderables::RenderablesUnion...; kwargs...)::Renderable
    renderables = center(renderables...)
    return vstack(renderables...; kwargs...)
end

"""
    lvstack(renderables::RenderablesUnion...)

Right align renderables and then vertically stack.
"""
function rvstack(renderables::RenderablesUnion...; kwargs...)::Renderable
    renderables = rightalign(renderables...)
    return vstack(renderables...; kwargs...)
end

rvstack(renderables::Union{Tuple,Vector}; kwargs...) = rvstack(renderables...; kwargs...)
cvstack(renderables::Union{Tuple,Vector}; kwargs...) = cvstack(renderables...; kwargs...)
lvstack(renderables::Union{Tuple,Vector}; kwargs...) = lvstack(renderables...; kwargs...)

# ---------------------------------------------------------------------------- #
#                                LINES & SPACER                                #
# ---------------------------------------------------------------------------- #
abstract type AbstractLayoutElement <: AbstractRenderable end

# ---------------------------------- spacer ---------------------------------- #
"""
        Spacer

A box of empty text with given width and height.
"""
mutable struct Spacer <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

function Spacer(width::Int, height::Int; char::Char = ' ')
    line = char^width
    seg = Segment(line)
    segments = repeat([seg], height)
    return Spacer(segments, Measure(seg.measure.w, height))
end

function Spacer(width::Number, height::Number; char::Char = ' ')
    return Spacer(int(width), int(height); char = char)
end

# ----------------------------------- vline ---------------------------------- #
"""
    vLine

A multi-line renderable with each line made of a | to create a vertical line.
"""
mutable struct vLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

"""
    vLine(height::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a `vLine` given a height and, optionally, style information.
"""
function vLine(
    height::Int;
    style::String = "default",
    box::Symbol = :ROUNDED,
    char::Union{Char,Nothing} = nothing,
)
    char = isnothing(char) ? getfield(Boxes, box).head.left : char
    line = apply_style("{" * style * "}" * char * "{/" * style * "}\e[0m")
    segments = repeat([Segment(line)], height)
    return vLine(segments, Measure(1, height))
end

"""
    vLine(ren::AbstractRenderable; kwargs...)

Construct a vLine with the same height as a renderable.
"""
vLine(ren::AbstractRenderable; kwargs...) = vLine(ren.measure.h; kwargs...)

"""
    vLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Create a `vLine` as tall as the `stdout` console.
"""
function vLine(; style::String = "default", box::Symbol = :ROUNDED)
    return vLine(console_height(stdout); style = style, box = box)
end

# ----------------------------------- hLine ---------------------------------- #

"""
    hLine

A 1-line renderable made of repeated character from a Box.
"""
mutable struct hLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

"""
    hLine(width::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a styled `hLine` of given width.
"""
function hLine(width::Int; style::String = "default", box::Symbol = :ROUNDED)
    char = eval(box).row.mid
    segments = [Segment(char^width * "\e[0m", style)]
    return hLine(segments, Measure(width, 1))
end

"""
    hLine(width::Number, text::String; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Creates an hLine object with texte centered horizontally.
"""
function hLine(
    width::Number,
    text::String;
    style::String = "default",
    box::Symbol = :ROUNDED,
)
    box = eval(box)
    text = apply_style(text) * "\e[0m"
    tl, tr = get_lr_widths(textlen(text))
    lw, rw = get_lr_widths(width)

    line =
        get_lrow(box, lw - tl, :top; with_left = false) *
        " " *
        text *
        " " *
        "{$style}" *
        get_rrow(box, rw - tr, :top; with_right = false) *
        "\e[0m"

    return hLine([Segment(line, style)], Measure(width, 1))
end

"""
    hLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout`.
"""
hLine(; style::String = "default", box::Symbol = :ROUNDED) =
    hLine(console_width(stdout); style = style, box = box)

"""
    hLine(text::AbstractString; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout` with centered text.
"""
hLine(text::AbstractString; style::String = "default", box::Symbol = :ROUNDED) =
    hLine(console_width(stdout), text; style = style, box = box)

"""
    hLine(ren::AbstractRenderable; kwargs)

Construct an hLine with same width as a renderable.
"""
hLine(ren::AbstractRenderable; kwargs...) = hLine(ren.measure.w; kwargs...)

# -------------------------------- PlaceHolder ------------------------------- #

"""
    mutable struct PlaceHolder <: AbstractLayoutElement
        segments::Vector{Segment}
        measure::Measure
    end
A `renderable` used as a place holder when creating layouts (e.g. with `grid`).

# Examples
```julia

println(PlaceHolder(25, 10))

╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ (25 × 10) ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲
```
"""
mutable struct PlaceHolder <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

function place_holder_line(w, i)
    if w > 1
        if w % 2 == 0
            left, right = w ÷ 2, w ÷ 2
            line = i == 0 ? "╲ "^(left) : " ╲"^(right)
        else
            left = cint(w / 2)
            line = i == 0 ? "╲ "^(left) : " ╲"^(left)
            line = ltrim_str(line, w)
        end

    else
        line = i == 0 ? "╲" : " "
    end
    return line
end

"""
    PlaceHolder(
        w::Int,
        h::Int;
        style::String = "dim",
        text::Union{Nothing,String} = nothing,
    )

Create a `PlaceHolder` with additional style information.
"""
function PlaceHolder(
    w::Int,
    h::Int;
    style::String = "dim",
    text::Union{Nothing,String} = nothing,
)
    # create lines of renderable
    b1 = place_holder_line(w, 0)
    b2 = place_holder_line(w, 1)
    lines::Vector{Segment} = map(i -> Segment(i % 2 != 0 ? b1 : b2, style), 1:h)

    # insert renderable size at center
    text = isnothing(text) ? "($w × $h)" : text
    text = "  " * text * "  "
    text = width(text) % 2 == 0 ? " " * text : text
    l = width(text)

    if l < (w / 2) && w > 13
        f = w < 30 ? 2.5 : 3
        _w = cint(ncodeunits(lines[1].text) / f)
        _l = cint(l / 2)

        original = lines[cint(h / 2)].text
        lines[cint(h / 2)] = Segment(
            ltrim_str(original, _w - _l) *
            "{default bold white}" *
            text *
            "{/default bold white}" *
            rtrim_str(original, _w + _l),
        )
    end

    return PlaceHolder(lines, Measure(lines))
end

PlaceHolder(ren::AbstractRenderable) = PlaceHolder(ren.measure.w, ren.measure.h)

# ---------------------------------------------------------------------------- #
#                                     GRID                                     #
# ---------------------------------------------------------------------------- #

"""
    grid(
        rens::Union{Nothing,AbstractVector{<:AbstractRenderable}} = nothing;
        placeholder::Union{Nothing,AbstractRenderable} = nothing,
        aspect::Union{Number,NTuple} = (16, 9),
        layout::Union{Nothing,Tuple} = nothing,
        pad::Union{Tuple,Integer} = 0,
    )

Construct a square grid from a `AbsractVector` of `AbstractRenderable`.
"""
function grid(
    rens::Union{Nothing,AbstractVector{<:AbstractRenderable}} = nothing;
    placeholder::Union{Nothing,AbstractRenderable} = nothing,
    aspect::Union{Number,NTuple} = (16, 9),
    layout::Union{Nothing,Tuple} = nothing,
    pad::Union{Tuple,Integer} = 0,
)
    nrows, ncols = isnothing(layout) ? calc_nrows_ncols(length(rens), aspect) : layout
    if isnothing(rens)
        isnothing(layout) &&
            throw(ArgumentError("`layout` must be given as `Tuple` of `Integer`s"))
        rens = fill(PlaceHolder(aspect...), prod(layout))
    else
        if isnothing(nrows)
            nrows, r = divrem(length(rens), ncols)
            r == 0 || (nrows += 1)
        elseif isnothing(ncols)
            ncols, r = divrem(length(rens), nrows)
            r == 0 || (ncols += 1)
        end
        fill_in = isnothing(placeholder) ? PlaceHolder(first(rens)) : placeholder
        rens = vcat(rens, repeat([fill_in], nrows * ncols - length(rens)))
    end
    return grid(reshape(rens, nrows, ncols); pad = pad)
end

"""
    grid(rens::AbstractMatrix{<:AbstractRenderable}; pad::Union{Nothing,Integer} = 0))

Construct a grid from a `AbstractMatrix` of `AbstractRenderable`.
"""
function grid(rens::AbstractMatrix{<:AbstractRenderable}; pad::Union{Tuple,Integer} = 0)
    hpad, vpad = if pad isa Integer
        (pad, pad)
    else
        pad
    end
    rows = collect(
        foldl((a, b) -> a * ' '^hpad * b, col[2:end]; init = first(col)) for
        col in eachrow(rens)
    )
    if vpad > 0
        vspace = vpad > 1 ? vstack(repeat([" "], vpad)...) : " "
        cat = (a, b) -> a / vspace / b
    else
        cat = (a, b) -> a / b
    end
    return foldl(cat, rows[2:end]; init = first(rows))
end

end


text:
#=
	author: Jonathan Payne
=#

using Roots

function rk4step(fdash, t, x, h)
    # Calculates the the nth step in the classical fourth order Runga-Kutta
    # procedure for solving a differential equation of the form:
    #   x'(t) = fdash(x, t),
    # Inputs:
    #   * fdash: user defined function
    #   * t: previous value of independent variable
    #   * x: previous value of dependent variable
    #   * h: increment
    # Outputs:
    #   * xnew: next step
    #   * tnew: new time
    k1 = h*fdash(x, t)
    k2 = h*fdash(x + k1/2, t + h/2)
    k3 = h*fdash(x + k2/2, t + h/2)
    k4 = h*fdash(x + k3, t + h)
    xnew = x + (k1 + 2*k2 + 2k3 + k4) / 6
    tnew = t + h
    return xnew, tnew
end

function rk4uniform(fdash, x0, tspan, N=100)
    # Approximate the solution of the initial value problem using the
    # classical fourth order Runge-Kutta method with uniform step sizes:
    #   * x'(t) = fdash(x, t),
    #   * x(t0) = x0
    # Inputs:
    #   * fdash: user defined function
    #   * t0: initial value of the independent variable
    #   * x0: initial value of the dependent variable(s) (column vector)
    #   * tf: final value of the independent variable
    #   * N: number of universally sized time steps to be taken to advance the
    #        the solution from t = t0 to t = tf
    # Outputs:
    #   * wi: vector / matrix containing the values of the approximate solution
    #         to the ode
    #   * ti: number of uniformly sized time steps to be taken to advance t
    t0, tf = tspan
    neqn = length(x0)
    ti = linspace(t0, tf, N+1)
    wi = zeros(neqn, N+1)
    wi[1:neqn, 1] = x0
    h = (tf - t0) / N
    for i = 1:N
        x0, t0 = rk4step(fdash, t0, x0, h)
        wi[1:neqn, i+1] = x0
    end
    return wi, ti
end

# ------------------------------
# Usage example 1:
# ------------------------------
t0, x0, tf, N = 0, [0,-4], 10, 100
A = [[1 2], [3 2]]
fdash(x, t) = A*x
xtrue(t) = [ 8/5*exp(-t) - 8/5*exp(4*t), -8/5*exp(-t) - 12/5*exp(4*t) ]
tspan = [t0, tf]
wi2, ti2 = rk4uniform(fdash, x0, tspan, N)

using ODE
t0, x0, tf, N = 0, [0,-4], 10, 100
A = [[1 2], [3 2]]
fdash(t, x) = A*x
xtrue(t) = [ 8/5*exp(-t) - 8/5*exp(4*t), -8/5*exp(-t) - 12/5*exp(4*t) ]
tspan = [t0, tf]
span = linspace(t0, tf, N)
ode23s(fdash, x0, span)

# ------------------------------
# odelay (from the pycse package)
# ------------------------------
# Solve an ODE with events.
#
# Inputs:
# func        = is callable, with signature func(Y, x)
# y0          = are the initial conditions
# xspan       = what you want to integrate over
# events      = is a list of callable functions with signature event(Y, x).
# 			    These functions return zero when an event has happened.
# [value, isterminal, direction] = event(Y, x)
# value       = the value of the event function. When value = 0, an event
# 			    is triggered
# isterminal  = True if the integration is to terminate at a zero of
# 			    this event function, otherwise, False.
# direction   = 0 if all zeros are to be located (the default),
#			  = +1 if only zeros where the event function is increasing, and
#             = -1 if only zeros where the event function is decreasing.
# TOLERANCE   = what is used to identify when an event has occurred.
# fsolve_args = a dictionary of options for fsolve
# kwargs  	  = additional options you want to send to odeint.
#
# Returns [x, y, te, ye, ie]
# x           = the independent variable array
# y           = the solution
# te          = an array of independent variable values where events occurred
# ye          = an array of the solution at the points where events occurred
# ie          = an array of indices indicating which event function occurred.

function odelayold(func, y0, xspan,
				 events, TOLERANCE = 1e-6)#, fsolve_args=[], args...)
	#if fsolve_args == None
    #    fsolve_args = []
    #end
    N = length(xspan)
    x0 = xspan[1] # initial point

    sol = [float64(y0)]
    TE, YE, IE = [], [], [] # to store where events occur

    # initial value of events
    e = zeros(length(events), N)
    for (i, event) in enumerate(events)
        e[i, 1], isterminal, direction = event(y0, x0)
    end

    # now we step through the integration
	for (i, x1) in enumerate(xspan[1:end-1])
		println("i", i)
        x2 = xspan[i + 1]
        f1 = sol[i]

        f2, wz = rk4uniform(func, f1', [x1, x2], 1)#, args...)

        #X[i+1] = x2
        #println(X[i+1])
        X = push!(X, x2)
        sol = push!(sol, f2[:,end][1])

        # check event functions. At each step we compute the event
        # functions, and check if they have changed sign since the
        # last step. If they changed sign, it implies a zero was
        # crossed.
        for (j, event) in enumerate(events)
     		e[j, i + 1], isterminal, direction = event(sol[i + 1], X[i + 1])

            if ((e[j, i + 1] * e[j, i] < 0) || abs(e[j, i + 1]) < TOLERANCE || abs(e[j, i]) < TOLERANCE)
                xLt = X[end]       # Last point
                fLt = sol[end]
                eLt = e[j, i+1]

	            function objective(x)
	            	# evaluate ode from xLT to x
	             	txspan = [xLt, x]
	             	tempsol, wyz = rk4uniform(func, fLt, txspan, length(txspan))#, args...)
	             	tsol = tempsol[:, end]
	             	val, isterminal, direction = event(tsol, x)
	             	return val[1]
	            end
	            xZ = fzero(objective, xLt)#, **fsolve_args)  # this should be the
                                              # value of x that makes
                                              # the event zero
                # now evaluate solution at this point, so we can
                # record the function values here.
                txspan = [xLt, xZ]
                tempsol, wyz = rk4uniform(func, fLt, txspan)#, args...)
                fZ = tempsol[:,end]

                vZ, isterminal, direction = event(fZ, xZ)

                COLLECTEVENT = false
                if direction == 0
                    COLLECTEVENT = true
                elseif (e[j, i + 1] > e[j, i] ) && direction == 1
                    COLLECTEVENT = true
                elseif (e[j, i + 1] < e[j, i] ) && direction == -1
                    COLLECTEVENT = true
                end

                if COLLECTEVENT
                    push!(TE, xZ)
                    push!(YE, fZ)
                    push!(IE, j)

                    if isterminal
                        X[end] = xZ[1]
                        sol[end] = fZ[1]
                        return X, sol, TE, YE, IE
                    end
                end

            end
        end
    end
    # at the end, return what we have
    return X, sol, TE, YE, IE

end

function odelay(func, y0, xspan,
                 events, TOLERANCE = 1e-6)#, fsolve_args=[], args...)
    #if fsolve_args == None
    #    fsolve_args = []
    #end
    N = length(xspan)
    x0 = xspan[1] # initial point

    X = [x0]
    Nf = length(y0)
    sol = zeros(N, Nf)
    sol[1,:] = float64(y0)
    TE, YE, IE = [], [], [] # to store where events occur

    # initial value of events
    e = zeros(length(events), N)
    for (i, event) in enumerate(events)
        e[i, 1], isterminal, direction = event(y0, x0)
    end

    # now we step through the integration
    for (i, x1) in enumerate(xspan[1:end-1])
        x2 = xspan[i + 1]
        f1 = sol[i, :]

        f2, wz = rk4uniform(func, f1', [x1, x2], 1)#, args...)

        #X[i+1] = x2
        #println(X[i+1])
        X = push!(X, x2)
        sol[i+1,:] = f2[:, end]

        # check event functions. At each step we compute the event
        # functions, and check if they have changed sign since the
        # last step. If they changed sign, it implies a zero was
        # crossed.
        for (j, event) in enumerate(events)
            e[j, i + 1], isterminal, direction = event(sol[i + 1], X[i + 1])

            if ((e[j, i + 1] * e[j, i] < 0) || abs(e[j, i + 1]) < TOLERANCE || abs(e[j, i]) < TOLERANCE)
                xLt = X[end]       # Last point
                fLt = sol[i+1,:]
                eLt = e[j, i+1]

                function objective(x)
                    # evaluate ode from xLT to x
                    txspan = [xLt, x]
                    tempsol, wyz = rk4uniform(func, fLt', txspan, length(txspan))#, args...)
                    tsol = tempsol[:, end]
                    val, isterminal, direction = event(tsol, x)
                    return val[1]
                end
                xZ = fzero(objective, xLt)#, **fsolve_args)  # this should be the
                                              # value of x that makes
                                              # the event zero
                # now evaluate solution at this point, so we can
                # record the function values here.
                txspan = [xLt, xZ]
                tempsol, wyz = rk4uniform(func, fLt', txspan, length(txspan))#, args...)
                fZ = tempsol[:,end]

                vZ, isterminal, direction = event(fZ, xZ)

                COLLECTEVENT = false
                if direction == 0
                    COLLECTEVENT = true
                elseif (e[j, i + 1] > e[j, i] ) && direction == 1
                    COLLECTEVENT = true
                elseif (e[j, i + 1] < e[j, i] ) && direction == -1
                    COLLECTEVENT = true
                end

                if COLLECTEVENT
                    push!(TE, xZ)
                    push!(YE, fZ)
                    push!(IE, j)

                    if isterminal
                        X[end] = xZ[1]
                        sol[i+1,:] = fZ[1]
                        return X, sol[1:i+1,:], TE, YE, IE
                    end
                end

            end
        end
    end
    # at the end, return what we have
    return X, sol, TE, YE, IE

end

# ------------------------------
# Example usage : 1 dimension
# ------------------------------

# myode(f, x) = 3*x^2 + 12*x -4

# function event1(f, x)
#     # An event is when f = 0 and event is decreasing'
#     isterminal = true
#     direction = -1
#     return f, isterminal, direction
# end

# function event2(f, x)
#     # An event is when f = 0 and increasing'
#     isterminal = false
#     direction = 1
#     return f, isterminal, direction
# end

# f0 = -120

# xspan = linspace(-8, 4, 50)
# events=[event1, event2]
# func = myode; y0 = f0;
# X, F, TE, YE, IE = odelay(myode, f0, xspan, events)

# ------------------------------
# Example usage : 2 dimensions
# ------------------------------

# function myode(z, t)
#     # PROJ: ODE for projectile motion with linear air resistance
#     # z = [z0, z1, z2, z3]
#     g = 9.81 # Specify approximate graviatational constant
#     b = 0.28 # Representative value
#     return [z[3], z[4], -g-b*z[3], -b*z[4]]
# end

# function linevent(z, t)
#     # LINEVENT: Contains the event we are looking for
#     # In this event, z(1) = 0 (hitting the ground)
#     lookfor = z[1] # Sets this to 0
#     isterminal = true # Stop when event is located
#     direction = -1 # specify downward direction
#     return lookfor, isterminal, direction
# end

# f0 = [0.18, 0, 5, 10]
# xspan = linspace(0, 10, 50)
# events = [linevent]
# func = myode
# y0 = f0
# Xm, Fm, TEm, YEm, IEm = odelay(myode, f0, xspan, events)

max_stars_repo_name: jepayne/NYU

max_issues_repo_name: jepayne/NYU

max_forks_repo_name: jepayne/NYU


text:
module Config

"""
    Configuration

Interface used by the `CentralNode` and the `Node` to work with any
configuration.
"""
abstract type Configuration end

# --------------------------------
# Include
# --------------------------------

include("base_config.jl")
include("vanilla_config.jl")
include("quantized_config.jl")
include("qdiff_config.jl")
include("gd_config.jl")

# --------------------------------
# Export
# --------------------------------

export Configuration, BaseConfig, VanillaConfig, QuantizedConfig,
    QDiffConfig, QDiffStaticConfig, GDConfig, GDStaticConfig

end # module


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

include("positivity_zhang_shu.jl")

end # @muladd


text:
@testset "GrpAbFinGen" begin
  @testset "Type stuff" begin
    @test elem_type(GrpAbFinGen) == GrpAbFinGenElem
    @test parent_type(GrpAbFinGenElem) == GrpAbFinGen
  end

  @testset "Constructor" begin
    M1 = matrix(FlintZZ, 2, 3, [1, 2, 3, 4, 5, 6])
    G = @inferred AbelianGroup(M1)
    @test isa(G, GrpAbFinGen)
    @test G.rels == M1

    M = FlintZZ[1 2 3; 4 5 6] # fmpz_mat
    G = @inferred AbelianGroup(M)
    @test isa(G, GrpAbFinGen)
    @test G.rels == M1

    M = fmpz[1 2 3; 4 5 6]
    G = @inferred AbelianGroup(M)
    @test isa(G, GrpAbFinGen)
    @test G.rels == M1

    M = [1 2 3; 4 5 6]
    G = @inferred AbelianGroup(M)
    @test isa(G, GrpAbFinGen)
    @test G.rels == M1

    M = fmpz[1, 2, 3, 4, 5, 6]
    G = @inferred AbelianGroup(M)
    @test isa(G, GrpAbFinGen)
    @test G.rels == matrix(FlintZZ, 1, 6, M)

    M = [1, 2, 3, 4, 5, 6]
    G = @inferred AbelianGroup(M)
    @test isa(G, GrpAbFinGen)
    @test G.rels == matrix(FlintZZ, 1, 6, M)

    M = [3, 0]
    G = @inferred DiagonalGroup(M)
    @test isa(G, GrpAbFinGen)

    M = fmpz[3, 0]
    G = @inferred DiagonalGroup(M)
    @test isa(G, GrpAbFinGen)

    M = matrix(FlintZZ, 1, 2, [3, 0])
    G = @inferred DiagonalGroup(M)
    @test isa(G, GrpAbFinGen)

    N = [3, 5]
    G = @inferred DiagonalGroup(N)
    @test isa(G, GrpAbFinGen)
    @test G.rels == matrix(FlintZZ, 2, 2, [3, 0, 0, 5])

    N = fmpz[3, 5]
    G = @inferred DiagonalGroup(N)
    @test isa(G, GrpAbFinGen)
    @test G.rels == matrix(FlintZZ, 2, 2, [3, 0, 0, 5])

    N = matrix(FlintZZ, 1, 2, [3, 5])
    G = @inferred DiagonalGroup(N)
    @test isa(G, GrpAbFinGen)
    @test G.rels == matrix(FlintZZ, 2, 2, [3, 0, 0, 5])

    @test_throws ErrorException DiagonalGroup(FlintZZ[1 2; 3 4])
    @test issnf(M)
  end

  @testset "String I/O" begin
    N = matrix(FlintZZ, 1, 2, [3, 5])
    G = DiagonalGroup(N)
    @test isa(string(G), String)
  end

  @testset "Field access" begin
    S = DiagonalGroup([3, 0])
    @test @inferred issnf(S)
    @test @inferred ngens(S) == 2
    @test @inferred nrels(S) == 2
    @test @inferred rels(S) == matrix(FlintZZ, 2, 2, [3, 0, 0, 0])

    G = DiagonalGroup([3, 5])
    @test @inferred !issnf(G)
    @test @inferred ngens(G) == 2
    @test @inferred nrels(G) == 2
    @test @inferred rels(G) == matrix(FlintZZ, 2, 2, [3, 0, 0, 5])
  end

  @testset "Hermite normal form" begin
    M   = FlintZZ[1 2 3; 4 5 6]
    HNF = FlintZZ[1 2 3; 0 3 6]
    G = AbelianGroup(M)
    Hecke.assure_has_hnf(G)
    @test G.hnf == HNF
  end

  @testset "Smith normal form" begin
    M = FlintZZ[16 17 2 ; 19 23 8 ; 16 17 2]
    G = AbelianGroup(M)
    S, mS = @inferred snf(G)
    @test issnf(S)
    @test S.snf == fmpz[45, 0]
    @test codomain(mS) == G
    @test domain(mS) == S
  end

  @testset "Finiteness" begin
    G = DiagonalGroup([3, 15])
    @test issnf(G)
    @test @inferred isfinite(G)
    @test @inferred !isinfinite(G)

    G = DiagonalGroup([3, 5])
    @test @inferred isfinite(G)
    @test @inferred !isinfinite(G)

    G = DiagonalGroup([3, 15, 0])
    @test issnf(G)
    @test @inferred !isfinite(G)
    @test @inferred isinfinite(G)

    G = DiagonalGroup([3, 5, 0])
    @test @inferred !isfinite(G)
    @test @inferred isinfinite(G)
  end

  @testset "Rank" begin
    G = DiagonalGroup([3, 15])
    @test @inferred rank(G) == 0

    G = DiagonalGroup([3, 5])
    @test @inferred rank(G) == 0

    G = DiagonalGroup([3, 15, 0])
    @test @inferred rank(G) == 1

    G = DiagonalGroup([3, 5, 0])
    @test @inferred rank(G) == 1
  end

  @testset "Order" begin
    G = DiagonalGroup([3, 5])
    @test @inferred order(G) == 15
    G = DiagonalGroup([3, 15])
    @test @inferred order(G) == 45
    G = DiagonalGroup([3, 5, 0])
    @test_throws ErrorException order(G)
  end

  @testset "Exponent" begin
    G = DiagonalGroup([3, 5])
    @test @inferred exponent(G) == 15
    G = DiagonalGroup([3, 15])
    @test @inferred exponent(G) == 15
  end

  @testset "Trivial" begin
    G = DiagonalGroup([1])
    @test @inferred istrivial(G)
    G = DiagonalGroup([1, 1, 1])
    @test @inferred istrivial(G)
    G = DiagonalGroup([3, 3])
    @test @inferred !istrivial(G)
    G = DiagonalGroup([3, 5])
    @test @inferred !istrivial(G)
  end

  @testset "Isomorphism" begin
    b = @inferred isisomorphic(DiagonalGroup(Int[]), DiagonalGroup(Int[]))
    @test b

    G = DiagonalGroup([2, 3, 5])
    H = DiagonalGroup([30])
    @test @inferred isisomorphic(G, H)
  end

  @testset "Direct product" begin
    G = DiagonalGroup([5, 3])
    H = DiagonalGroup([4])
    K = direct_product(G, H)[1]
    @test isisomorphic(K, DiagonalGroup([60]))
  end

  @testset "Torsion" begin
    G = DiagonalGroup([5, 4])
    @test @inferred istorsion(G)
    H, mH = torsion_subgroup(G)
    @test order(H) == 20

    G = DiagonalGroup([5, 0, 4, 0])
    @test @inferred !istorsion(G)
    H, mH = torsion_subgroup(G)
    @test isisomorphic(H, DiagonalGroup([5, 4]))
  end

  @testset "Subgroup" begin
    @test_throws ErrorException sub(GrpAbFinGenElem[])

    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0])
    g1 = G[1]
    g2 = G[2]
    g3 = G[3]
    S, S_map = @inferred sub([g1, g2, g3])
    @test isisomorphic(G, S)

    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0])
    S, mS = snf(G)
    s1 = S[1]
    s2 = S[2]
    s3 = S[3]
    H, mH = @inferred sub(S, [s1, s2, s3])
    @test isisomorphic(H, G)

    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0])
    g1 = G[1]
    H, mH = @inferred sub(G, [g1])
    @test isisomorphic(H, DiagonalGroup([3]))

    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0])
    S, mS = snf(G)
    s1 = S[1]
    H, mH = @inferred sub(S, [s1])
    @test isisomorphic(H, DiagonalGroup([3]))

    # G contains empty relation
    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0 ; 0 0 30 ; 0 0 0])
    g1 = G[3]
    S, mS = @inferred sub(G, [g1])
    @test isisomorphic(S, DiagonalGroup([30]))

    # n*G

    G = DiagonalGroup([6, 6, 12, 5])
    H, mH = @inferred sub(G, 2)
    @test isisomorphic(H, DiagonalGroup([3, 3, 6, 5]))

    H, mH = @inferred sub(G, fmpz(2))
    @test isisomorphic(H, DiagonalGroup([3, 3, 6, 5]))
  end

  @testset "Quotient" begin
    G = AbelianGroup(FlintZZ[3 0 0 ; 0 15 0])

    Q, mQ = @inferred quo(G, GrpAbFinGenElem[])
    @test isisomorphic(Q, G)

    g2 = G[2]
    Q, mQ = @inferred quo(G, [g2])
    @test isisomorphic(Q, DiagonalGroup([3, 0]))

    S = DiagonalGroup([3, 15, 0])
    @test issnf(S)
    g2 = S[2]
    Q, mQ = @inferred quo(S, [g2])
    @test isisomorphic(Q, DiagonalGroup([3, 0]))

    G = DiagonalGroup([6, 6, 12, 5, 0])
    H, mH = @inferred quo(G, 2)
    @test isisomorphic(H, DiagonalGroup([2, 2, 2, 2]))

    H, mH = @inferred quo(G, fmpz(2))
    @test isisomorphic(H, DiagonalGroup([2, 2, 2, 2]))
  end

  @testset "Cyclic" begin
    G = DiagonalGroup([3, 5])
    @test @inferred iscyclic(G)

    G = DiagonalGroup([3, 15])
    @test @inferred !iscyclic(G)
  end

  @testset "p-Sylow subgroup" begin
    G = DiagonalGroup([1, 3, 9, 5, 15, 20, 7])
    P, mP = psylow_subgroup(G, 3)
    @test order(P) == 3^valuation(order(G), 3)
    P, mP = psylow_subgroup(G, 5)
    @test order(P) == 5^valuation(order(G), 5)
    P, mP = psylow_subgroup(G, 11)
    @test order(P) == 11^valuation(order(G), 11)
  end

  @testset "All abelian groups" begin
    l = collect(abelian_groups(1))
    @test length(l) == 1
    l = collect(abelian_groups(2))
    @test length(l) == 1
    l = collect(abelian_groups(2^6 * 3^3 * 5^2 * 7))
    @test length(l) == 66
  end

  @testset "HomAlg" begin
    G = AbelianGroup([3 1; 0 3])
    S, mS = snf(G)
    H, mH = hom(G,G)
    T, p = tensor_product(G, G)
    D = direct_product(G, G, task = :none)
    for i=1:5
      Th = hom(T, T, map(mH, [rand(H) for x = 1:2])) #induced map in tensor product
      Dh = hom(D, D, map(mH, rand(H, (2,2)))) #induced map on direct prod
    end
    C = free_resolution(G)
    D = free_resolution(S)
    @test isexact(C)
    hom(D, C, mS)
    E = hom(D, T)
    @test 9 == prod(order(x) for x = homology(E))
    @test !isexact(E)
    E = hom(T, C)
    @test !isexact(E)
    E = tensor_product(C, T)
    @test !isexact(E)

    A = AbelianGroup([3 1; 0 3])
    B = AbelianGroup([9 2 1; 0 12 1; 0 0 25])
    C = DiagonalGroup([3, 4, 0])
    @test isisomorphic(hom(tensor_product(A, B, task = :none), C)[1],
                       hom(A, hom(B, C)[1])[1])
  end

end


text:
using Base: Float64
const π = 3.1415926535
const μ = 39.4234021
const a = 0.38709893
const α = 1.09778201e-8
const e = 0.20563069

rad_to_arcsec(radians::Float64)::Float64 = radians*180/π*60*60

struct Point
	x::Float64
	y::Float64
	Point(x,y) = new(x,y)
	Point() = new(0.0,0.0)
end

norm(p::Point)::Float64 = hypot(p.x, p.y)
angle(p::Point)::Float64 = mod(atan(p.y, p.x),2π)
dot(p1::Point,p2::Point)::Float64 = p1.x*p2.x + p1.y*p2.y

import Base: +,*,/
(p1::Point + p2::Point)::Point = Point(p1.x+p2.x, p1.y+p2.y)
(p::Point * scalar)::Point = Point(p.x*scalar,p.y*scalar)
(scalar * p::Point)::Point = p*scalar
(p::Point / scalar)::Point = Point(p.x/scalar,p.y/scalar)

function acceleration(r_vec::Point)::Point
	r::Float64 = norm(r_vec)
	β::Float64 = -μ/r^3 * (1+α/r^2)
	return r_vec * β
end

function verlet(tmax::Float64,dt_max::Float64,dt_min::Float64)::Nothing
	r_past = Point(a*(1+e),0.)
	v_past = Point(0.,sqrt((μ/a)*(1-e)/(1+e)))
	a_past = acceleration(r_past)

	r_now = Point()
	v_now = Point()
	a_now = Point()

	event_past = dot(r_past,v_past)
	event_now = 0.0

	dt = dt_max
	t = 0.0
	while t < tmax
		r_now = r_past + v_past*dt + a_past*dt^2/2
		a_now = acceleration(r_now)
		v_now = v_past + dt*(a_now+a_past)/2
		event_now = dot(r_now,v_now)

		if event_past<0 && event_now>0
			if dt<dt_min
				println(t," ",rad2deg(angle((r_now+r_past)/2)-π)*3600)
				dt = dt_max
			else
				# rewind time
				t -= dt
				r_now = r_past
				v_now = v_past
				a_now = a_past
				event_now = event_past
				# adjust dt
				dt /= 2
			end
		end

		r_past = r_now
		v_past = v_now
		a_past = a_now
		event_past = event_now
		t += dt
	end
	return nothing
end

# compile (?)
verlet(1e-1,1e-2,1e-3)

println()

verlet(10.0,1e-6,1e-12)


text:
@static if Sys.iswindows()
	ts = Int(floor(time()))
	if Sys.ARCH == :x86_64
		# 64-bit version is not available in an end-user package, so we download the SDK
		srcUrl = "http://downloads.sourceforge.net/project/openil/DevIL%20Windows%20SDK/1.7.8/DevIL-SDK-x64-1.7.8.zip?r=&ts=$ts&use_mirror=auto_select"
		fileName = "DevIL-SDK-x64-1.7.8.zip"
	elseif Sys.ARCH == :x86
		srcUrl = "http://downloads.sourceforge.net/project/openil/DevIL%20Win32/1.7.8/DevIL-EndUser-x86-1.7.8.zip?r=&ts=$ts&use_mirror=auto_select"
		fileName = "DevIL-EndUser-x86-1.7.8.zip"
	else
		error("DevIL: Unsupported Windows architecture")
	end

	dstDir = string(Sys.ARCH)
	if (!isdir(dstDir))
		mkdir(dstDir)
	end
	dstFile = dstDir * "\\" * fileName
	download(srcUrl, dstFile)
	run(`"$(Sys.BINDIR)\\7z.exe" e "$dstFile" *.dll -o"$dstDir" -y`)
	rm(dstFile)
end

@static if Sys.islinux()
    run(`sudo apt-get update`)
    run(`sudo apt-get install libdevil1c2`)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

#Period types
value(x::Period) = x.value

# The default constructors for Periods work well in almost all cases
# P(x) = new((convert(Int64,x))
# The following definitions are for Period-specific safety
for period in (:Year, :Month, :Week, :Day, :Hour, :Minute, :Second, :Millisecond)
    period_str = string(period)
    accessor_str = lowercase(period_str)
    # Convenience method for show()
    @eval _units(x::$period) = " " * $accessor_str * (abs(value(x)) == 1 ? "" : "s")
    # periodisless
    @eval periodisless(x::$period,y::$period) = value(x) < value(y)
    # AbstractString parsing (mainly for IO code)
    @eval $period(x::AbstractString) = $period(Base.parse(Int64,x))
    # Period accessors
    typ_str = period in (:Hour, :Minute, :Second, :Millisecond) ? "DateTime" : "TimeType"
    description = typ_str == "TimeType" ? "`Date` or `DateTime`" : "`$typ_str`"
    reference = period == :Week ? " For details see [`$accessor_str(::$typ_str)`](:func:`$accessor_str`)." : ""
    @eval begin
        @doc """
            $($period_str)(dt::$($typ_str)) -> $($period_str)

        The $($accessor_str) part of a $($description) as a `$($period_str)`.$($reference)
        """ ->
        $period(dt::$(Symbol(typ_str))) = $period($(Symbol(accessor_str))(dt))

        @doc """
            $($period_str)(v)

        Construct a `$($period_str)` object with the given `v` value. Input must be
        losslessly convertible to an `Int64`.
        """ $period(v)
    end
end
# Now we're safe to define Period-Number conversions
# Anything an Int64 can convert to, a Period can convert to
Base.convert{T<:Number}(::Type{T},x::Period) = convert(T,value(x))
Base.convert{T<:Period}(::Type{T},x::Real) = T(x)

#Print/show/traits
Base.string{P<:Period}(x::P) = string(value(x),_units(x))
Base.show(io::IO,x::Period) = print(io,string(x))
Base.zero{P<:Period}(::Union{Type{P},P}) = P(0)
Base.one{P<:Period}(::Union{Type{P},P}) = P(1)
Base.typemin{P<:Period}(::Type{P}) = P(typemin(Int64))
Base.typemax{P<:Period}(::Type{P}) = P(typemax(Int64))

# Default values (as used by TimeTypes)
"""
    default(p::Period) -> Period

Returns a sensible "default" value for the input Period by returning `one(p)` for Year,
Month, and Day, and `zero(p)` for Hour, Minute, Second, and Millisecond.
"""
function default end

default{T<:DatePeriod}(p::Union{T,Type{T}}) = one(p)
default{T<:TimePeriod}(p::Union{T,Type{T}}) = zero(p)

(-){P<:Period}(x::P) = P(-value(x))
Base.isless{P<:Period}(x::P,y::P) = isless(value(x),value(y))
=={P<:Period}(x::P,y::P) = value(x) == value(y)

# Period Arithmetic, grouped by dimensionality:
import Base: div, fld, mod, rem, gcd, lcm, +, -, *, /, %, .+, .-, .*, .%
for op in (:+,:-,:lcm,:gcd)
    @eval ($op){P<:Period}(x::P,y::P) = P(($op)(value(x),value(y)))
end

for op in (:/,:div,:fld)
    @eval begin
        ($op){P<:Period}(x::P,y::P) = ($op)(value(x),value(y))
        ($op){P<:Period}(x::P,y::Real) = P(($op)(value(x),Int64(y)))
    end
end

for op in (:rem,:mod)
    @eval begin
        ($op){P<:Period}(x::P,y::P) = P(($op)(value(x),value(y)))
        ($op){P<:Period}(x::P,y::Real) = P(($op)(value(x),Int64(y)))
    end
end

/{P<:Period}(X::StridedArray{P}, y::P) = X ./ y
%{P<:Period}(X::StridedArray{P}, y::P) = X .% y
*{P<:Period}(x::P,y::Real) = P(value(x) * Int64(y))
*(y::Real,x::Period) = x * y
.*{P<:Period}(y::Real, X::StridedArray{P}) = X .* y
for (op,Ty,Tz) in ((:.*,Real,:P),
                   (:./,:P,Float64), (:./,Real,:P),
                   (:div,:P,Int64), (:div,Integer,:P),
                   (:.%,:P,:P),
                   (:mod,:P,:P))
    sop = string(op)
    op_ = sop[1] == '.' ? Symbol(sop[2:end]) : op
    @eval begin
        function ($op){P<:Period}(X::StridedArray{P},y::$Ty)
            Z = similar(X, $Tz)
            for (Idst, Isrc) in zip(eachindex(Z), eachindex(X))
                @inbounds Z[Idst] = ($op_)(X[Isrc],y)
            end
            return Z
        end
    end
end

# intfuncs
Base.gcdx{T<:Period}(a::T,b::T) = ((g,x,y)=gcdx(value(a),value(b)); return T(g),x,y)
Base.abs{T<:Period}(a::T) = T(abs(value(a)))

periodisless(::Period,::Year)        = true
periodisless(::Period,::Month)       = true
periodisless(::Year,::Month)         = false
periodisless(::Period,::Week)        = true
periodisless(::Year,::Week)          = false
periodisless(::Month,::Week)         = false
periodisless(::Period,::Day)         = true
periodisless(::Year,::Day)           = false
periodisless(::Month,::Day)          = false
periodisless(::Week,::Day)           = false
periodisless(::Period,::Hour)        = false
periodisless(::Minute,::Hour)        = true
periodisless(::Second,::Hour)        = true
periodisless(::Millisecond,::Hour)   = true
periodisless(::Period,::Minute)      = false
periodisless(::Second,::Minute)      = true
periodisless(::Millisecond,::Minute) = true
periodisless(::Period,::Second)      = false
periodisless(::Millisecond,::Second) = true
periodisless(::Period,::Millisecond) = false

# return (next coarser period, conversion factor):
coarserperiod{P<:Period}(::Type{P}) = (P,1)
coarserperiod(::Type{Millisecond}) = (Second,1000)
coarserperiod(::Type{Second}) = (Minute,60)
coarserperiod(::Type{Minute}) = (Hour,60)
coarserperiod(::Type{Hour}) = (Day,24)
coarserperiod(::Type{Day}) = (Week,7)
coarserperiod(::Type{Month}) = (Year,12)

# Stores multiple periods in greatest to least order by type, not values,
# canonicalized to eliminate zero periods, merge equal period types,
# and convert more-precise periods to less-precise periods when possible
"""
    CompoundPeriod

A `CompoundPeriod` is useful for expressing time periods that are not a fixed multiple of
smaller periods. For example, \"a year and a  day\" is not a fixed number of days, but can
be expressed using a `CompoundPeriod`. In fact, a `CompoundPeriod` is automatically
generated by addition of different period types, e.g. `Year(1) + Day(1)` produces a
`CompoundPeriod` result.
"""
type CompoundPeriod <: AbstractTime
    periods::Array{Period,1}
    function CompoundPeriod(p::Vector{Period})
        n = length(p)
        if n > 1
            sort!(p, rev=true, lt=periodisless)
            # canonicalize p by merging equal period types and removing zeros
            i = j = 1
            while j <= n
                k = j+1
                while k <= n
                    if typeof(p[j]) == typeof(p[k])
                        p[j] += p[k]
                        k += 1
                    else
                        break
                    end
                end
                if p[j] != zero(p[j])
                    p[i] = p[j]
                    i += 1
                end
                j = k
            end
            n = i - 1 # new length
        elseif n == 1 && value(p[1]) == 0
            p = Period[]
            n = 0
        end
        # canonicalize Periods by pushing "overflow" into a coarser period.
        if n > 0
            pc = sizehint!(Period[], n)
            P = typeof(p[n])
            v = value(p[n])
            i = n - 1
            while true
                Pc, f = coarserperiod(P)
                if i > 0 && typeof(p[i]) == P
                    v += value(p[i])
                    i -= 1
                end
                v0 = f == 1 ? v : rem(v, f)
                v0 != 0 && push!(pc, P(v0))
                if v != v0
                    P = Pc
                    v = div(v - v0, f)
                elseif i > 0
                    P = typeof(p[i])
                    v = value(p[i])
                    i -= 1
                else
                    break
                end
            end
            p = reverse!(pc)
            n = length(p)
        else
            return new(resize!(p, n))
        end
        # reduce the amount of mixed positive/negative Periods.
        if n > 0
            pc = sizehint!(Period[], n)
            i = n
            while i > 0
                j = i

                # Determine sign of the largest period in this group which
                # can be converted into via coarserperiod.
                last = Union{}
                current = typeof(p[i])
                while i > 0 && current != last
                    if typeof(p[i]) == current
                        i -= 1
                    end
                    last, current = current, coarserperiod(current)[1]
                end
                s = sign(value(p[i + 1]))

                # Adjust all the periods in the group based upon the
                # largest period sign.
                P = typeof(p[j])
                v = 0
                while j > i
                    Pc, f = coarserperiod(P)
                    if j > 0 && typeof(p[j]) == P
                        v += value(p[j])
                        j -= 1
                    end
                    v0 = f == 1 ? v : mod(v, f * s)
                    v0 != 0 && push!(pc, P(v0))
                    if v != v0
                        P = Pc
                        v = div(v - v0, f)
                    elseif j > 0
                        P = typeof(p[j])
                        v = 0
                    else
                        break
                    end
                end
            end
            p = reverse!(pc)
        end
        return new(p)
    end
end

"""
    CompoundPeriod(periods) -> CompoundPeriod

Construct a `CompoundPeriod` from a `Vector` of `Period`s. The constructor will
automatically simplify the periods into a canonical form according to the following rules:

* All `Period`s of the same type will be added together
* Any `Period` large enough be partially representable by a coarser `Period` will be broken
  into multiple `Period`s (eg. `Hour(30)` becomes `Day(1) + Hour(6)`)
* `Period`s with opposite signs will be combined when possible
  (eg. `Hour(1) - Day(1)` becomes `-Hour(23)`)

Due to the canonicalization, `CompoundPeriod` is also useful for converting time periods
into more human-comprehensible forms.

# Examples
```julia
julia> Dates.CompoundPeriod([Dates.Hour(12), Dates.Hour(13)])
1 day, 1 hour

julia> Dates.CompoundPeriod([Dates.Hour(-1), Dates.Minute(1)])
-59 minutes

julia> Dates.CompoundPeriod([Dates.Month(1), Dates.Week(-2)])
1 month, -2 weeks

julia> Dates.CompoundPeriod(Dates.Minute(50000)))
4 weeks, 6 days, 17 hours, 20 minutes
```
"""
CompoundPeriod{P<:Period}(p::Vector{P}) = CompoundPeriod(Array{Period}(p))

Base.convert(::Type{CompoundPeriod}, x::Period) = CompoundPeriod(Period[x])
function Base.string(x::CompoundPeriod)
    if isempty(x.periods)
        return "empty period"
    else
        s = ""
        for p in x.periods
            s *= ", " * string(p)
        end
        return s[3:end]
    end
end
Base.show(io::IO,x::CompoundPeriod) = print(io,string(x))

# E.g. Year(1) + Day(1)
(+)(x::Period,y::Period) = CompoundPeriod(Period[x,y])
(+)(x::CompoundPeriod,y::Period) = CompoundPeriod(vcat(x.periods,y))
(+)(y::Period,x::CompoundPeriod) = x + y
(+)(x::CompoundPeriod,y::CompoundPeriod) = CompoundPeriod(vcat(x.periods,y.periods))
# E.g. Year(1) - Month(1)
(-)(x::Period,y::Period) = CompoundPeriod(Period[x,-y])
(-)(x::CompoundPeriod,y::Period) = CompoundPeriod(vcat(x.periods,-y))
(-)(x::CompoundPeriod) = CompoundPeriod(-x.periods)
(-)(y::Union{Period,CompoundPeriod},x::CompoundPeriod) = (-x) + y

GeneralPeriod = Union{Period,CompoundPeriod}
(+)(x::GeneralPeriod) = x
(+){P<:GeneralPeriod}(x::StridedArray{P}) = x

for op in (:.+, :.-)
    op_ = Symbol(string(op)[2:end])
    @eval begin
        function ($op){P<:GeneralPeriod}(X::StridedArray{P},y::GeneralPeriod)
            Z = similar(X, CompoundPeriod)
            for (Idst, Isrc) in zip(eachindex(Z), eachindex(X))
                @inbounds Z[Idst] = ($op_)(X[Isrc],y)
            end
            return Z
        end
        ($op){P<:GeneralPeriod}(x::GeneralPeriod,Y::StridedArray{P}) = ($op)(Y,x) |> ($op_)
        ($op_){P<:GeneralPeriod}(x::GeneralPeriod,Y::StridedArray{P}) = ($op)(Y,x) |> ($op_)
        ($op_){P<:GeneralPeriod}(Y::StridedArray{P},x::GeneralPeriod) = ($op)(Y,x)
        ($op_){P<:GeneralPeriod, Q<:GeneralPeriod}(X::StridedArray{P}, Y::StridedArray{Q}) =
            reshape(CompoundPeriod[($op_)(x,y) for (x,y) in zip(X, Y)], promote_shape(size(X),size(Y)))
    end
end

(==)(x::CompoundPeriod, y::Period) = x == CompoundPeriod(y)
(==)(y::Period, x::CompoundPeriod) = x == y
(==)(x::CompoundPeriod, y::CompoundPeriod) = x.periods == y.periods

# Capture TimeType+-Period methods
(+)(a::TimeType,b::Period,c::Period) = (+)(a,b+c)
(-)(a::TimeType,b::Period,c::Period) = (-)(a,b-c)
(+)(a::TimeType,b::Period,c::Period,d::Period...) = (+)((+)(a,b+c),d...)
(-)(a::TimeType,b::Period,c::Period,d::Period...) = (-)((-)(a,b-c),d...)

function (+)(x::TimeType,y::CompoundPeriod)
    for p in y.periods
        x += p
    end
    return x
end
(+)(x::CompoundPeriod,y::TimeType) = y + x

function (-)(x::TimeType,y::CompoundPeriod)
    for p in y.periods
        x -= p
    end
    return x
end
(-)(x::CompoundPeriod,y::TimeType) = y - x

# Fixed-value Periods (periods corresponding to a well-defined time interval,
# as opposed to variable calendar intervals like Year).
typealias FixedPeriod Union{Week,Day,Hour,Minute,Second,Millisecond}

# like div but throw an error if remainder is nonzero
function divexact(x,y)
    q,r = divrem(x, y)
    r == 0 || throw(InexactError())
    return q
end

# FixedPeriod conversions and promotion rules
const fixedperiod_conversions = [(Week,7),(Day,24),(Hour,60),(Minute,60),(Second,1000),(Millisecond,1)]
for i = 1:length(fixedperiod_conversions)
    (T,n) = fixedperiod_conversions[i]
    N = 1
    for j = i-1:-1:1 # less-precise periods
        (Tc,nc) = fixedperiod_conversions[j]
        N *= nc
        vmax = typemax(Int64) ÷ N
        vmin = typemin(Int64) ÷ N
        @eval function Base.convert(::Type{$T}, x::$Tc)
            $vmin ≤ value(x) ≤ $vmax || throw(InexactError())
            return $T(value(x)*$N)
        end
    end
    N = n
    for j = i+1:length(fixedperiod_conversions) # more-precise periods
        (Tc,nc) = fixedperiod_conversions[j]
        @eval Base.convert(::Type{$T}, x::$Tc) = $T(divexact(value(x), $N))
        @eval Base.promote_rule(::Type{$T},::Type{$Tc}) = $Tc
        N *= nc
    end
end
# have to declare thusly so that diagonal dispatch above takes precedence:
(==){T<:FixedPeriod,S<:FixedPeriod}(x::T,y::S) = (==)(promote(x,y)...)
Base.isless{T<:FixedPeriod,S<:FixedPeriod}(x::T,y::S) = isless(promote(x,y)...)

# other periods with fixed conversions but which aren't fixed time periods
typealias OtherPeriod Union{Month,Year}
let vmax = typemax(Int64) ÷ 12, vmin = typemin(Int64) ÷ 12
    @eval function Base.convert(::Type{Month}, x::Year)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError())
        Month(value(x)*12)
    end
end
Base.convert(::Type{Year}, x::Month) = Year(divexact(value(x),12))
Base.promote_rule(::Type{Year}, ::Type{Month}) = Month
(==){T<:OtherPeriod,S<:OtherPeriod}(x::T,y::S) = (==)(promote(x,y)...)
Base.isless{T<:OtherPeriod,S<:OtherPeriod}(x::T,y::S) = isless(promote(x,y)...)

# truncating conversions to milliseconds and days:
toms(c::Millisecond) = value(c)
toms(c::Second)      = 1000*value(c)
toms(c::Minute)      = 60000*value(c)
toms(c::Hour)        = 3600000*value(c)
toms(c::Day)         = 86400000*value(c)
toms(c::Week)        = 604800000*value(c)
toms(c::Month)       = 86400000.0*30.436875*value(c)
toms(c::Year)        = 86400000.0*365.2425*value(c)
toms(c::CompoundPeriod) = isempty(c.periods)?0.0 : Float64(sum(toms,c.periods))
days(c::Millisecond) = div(value(c),86400000)
days(c::Second)      = div(value(c),86400)
days(c::Minute)      = div(value(c),1440)
days(c::Hour)        = div(value(c),24)
days(c::Day)         = value(c)
days(c::Week)        = 7*value(c)
days(c::Year)        = 365.2425*value(c)
days(c::Month)       = 30.436875*value(c)
days(c::CompoundPeriod) = isempty(c.periods)?0.0 : Float64(sum(days,c.periods))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Libdl

oldwdir = pwd()

@show M4_VERSION = "1.4.17"
@show YASM_VERSION = "1.3.0"
@show MPIR_VERSION = "3.0.0"
@show MPFR_VERSION = "4.0.0"
@show FLINT_VERSION = "adf1583c6bd92a454f3f92a18adf9063d14637a0"
@show ARB_VERSION = "4ddd6771d94e602462e1fc470ec9d55b248c04d7"

pkgdir = dirname(dirname(@__FILE__))
wdir = joinpath(pkgdir, "deps")
vdir = joinpath(pkgdir, "local")

println("\n\t$(vdir)\n")

if Sys.isapple() && !("CC" in keys(ENV))
   ENV["CC"] = "clang"
   ENV["CXX"] = "clang++"
end

if !ispath(vdir)

    mkdir(vdir)

    if !ispath(joinpath(vdir, "lib"))
        mkdir(joinpath(vdir, "lib"))
    end
else
    println("Deleting old $vdir")
    rm(vdir, force=true, recursive=true)
    mkdir(vdir)
    mkdir(joinpath(vdir, "lib"))
end

LDFLAGS = "-Wl,-rpath,$vdir/lib -Wl,-rpath,\$\$ORIGIN/../share/julia/site/v$(VERSION.major).$(VERSION.minor)/Nemo/local/lib"
DLCFLAGS = "-fPIC -fno-common"

cd(wdir)

function download_dll(url_string, location_string)
   try
      run(`curl -o $(location_string) -L $(url_string)`)
   catch
      download(url_string, location_string)
   end
end

#install libpthreads

if Sys.iswindows()
   println("Downloading libpthread ... ")
   if Int == Int32
      download_dll("http://nemocas.org/binaries/w32-libwinpthread-1.dll", joinpath(vdir, "lib", "libwinpthread-1.dll"))
   else
      download_dll("http://nemocas.org/binaries/w64-libwinpthread-1.dll", joinpath(vdir, "lib", "libwinpthread-1.dll"))
   end
   println("DONE")
end

cd(wdir)

# install M4

if !Sys.iswindows()
   try
      run(`m4 --version`)
   catch
      println("Building m4 ... ")
      M4_FILE = "m4-" * M4_VERSION * ".tar.bz2"
      download("http://ftp.gnu.org/gnu/m4/$M4_FILE", joinpath(wdir, "$M4_FILE"))
      run(`tar -xvf $M4_FILE`)
      run(`rm $M4_FILE`)
      cd(joinpath("$wdir", "m4-$M4_VERSION"))
      run(`./configure --prefix=$vdir`)
      run(`make`)
      run(`make install`)
      println("DONE")
   end
end

cd(wdir)

# install yasm

if !Sys.iswindows()
   if !ispath(joinpath(wdir, "yasm-$YASM_VERSION"))
      println("Building yasm ... ")
      YASM_FILE = "yasm-" * YASM_VERSION * ".tar.gz"
      download("http://www.tortall.net/projects/yasm/releases/$YASM_FILE", YASM_FILE)
      run(`tar -xvf $YASM_FILE`)
      run(`rm $YASM_FILE`)
      cd(joinpath("$wdir","yasm-$YASM_VERSION"))
      run(`./configure`)
      run(`make`)
      println("DONE")
   end
end

cd(wdir)

# install GMP/MPIR

MPIR_FILE = "mpir-" * MPIR_VERSION * ".tar.bz2"

if !ispath(joinpath(wdir, "mpir-$MPIR_VERSION"))
   println("Downloading MPIR sources ... ")
   download("http://mpir.org/$MPIR_FILE", joinpath(wdir, MPIR_FILE))
   println("DONE")
end

if Sys.iswindows()
   println("Downloading MPIR ... ")
   if Int == Int32
      download_dll("http://nemocas.org/binaries/w32-libgmp-16.dll", joinpath(vdir, "lib", "libgmp-16.dll"))
   else
      download_dll("http://nemocas.org/binaries/w64-libgmp-16.dll", joinpath(vdir, "lib", "libgmp-16.dll"))
   end
   println("DONE")
else
   println("Building MPIR ... ")
   if isfile(joinpath(wdir, MPIR_FILE))
      run(`tar -xvf $MPIR_FILE`)
      run(`rm $MPIR_FILE`)
   end
   cd("$wdir/mpir-$MPIR_VERSION")
   try
      run(`m4 --version`)
      run(`./configure --with-yasm=$wdir/yasm-$YASM_VERSION/yasm --prefix=$vdir --enable-gmpcompat --disable-static --enable-shared`)
   catch
      run(`./configure --with-yasm=$wdir/yasm-$YASM_VERSION/yasm --prefix=$vdir M4=$vdir/bin/m4 --enable-gmpcompat --disable-static --enable-shared`)
   end
   run(`make -j4`)
   run(`make install`)
   cd(wdir)
   run(`rm -rf bin`)
   println("DONE")
end

cd(wdir)

# install MPFR

MPFR_FILE = "mpfr-" * MPFR_VERSION * ".tar.bz2"

if !ispath(joinpath(wdir, "mpfr-$MPFR_VERSION"))
   println("Downloading MPFR sources ... ")
   download("http://ftp.gnu.org/gnu/mpfr/$MPFR_FILE", joinpath(wdir, MPFR_FILE))

   println("DONE")
end

if Sys.iswindows()
   println("Downloading MPFR ... ")
   if Int == Int32
      download_dll("http://nemocas.org/binaries/w32-libmpfr-4.dll", joinpath(vdir, "lib", "libmpfr-4.dll"))
   else
      download_dll("http://nemocas.org/binaries/w64-libmpfr-4.dll", joinpath(vdir, "lib", "libmpfr-4.dll"))
   end
   println("DONE")
else
   println("Building MPFR ... ")
   if isfile(joinpath(wdir, MPFR_FILE))
      run(`tar -xvf $MPFR_FILE`)
      run(`rm $MPFR_FILE`)
   end
   cd("$wdir/mpfr-$MPFR_VERSION")
   withenv("LD_LIBRARY_PATH"=>"$vdir/lib", "LDFLAGS"=>LDFLAGS) do
      run(`./configure --prefix=$vdir --with-gmp=$vdir --disable-static --enable-shared`)
      run(`make -j4`)
      run(`make install`)
   end
   cd(wdir)
   println("DONE")
end

cd(wdir)

# install FLINT
if !Sys.iswindows()
  try
    println("Cloning flint2 ... ")
    run(`git clone https://github.com/wbhart/flint2.git`)
    cd(joinpath("$wdir", "flint2"))
    run(`git checkout $FLINT_VERSION`)
    cd(wdir)
  catch
    if ispath(joinpath("$wdir", "flint2"))
       cd(joinpath("$wdir", "flint2"))
       run(`git fetch`)
       run(`git checkout $FLINT_VERSION`)
       cd(wdir)
    end
  end
  println("DONE")
end

if Sys.iswindows()
   println("Downloading flint ... ")
   if Int == Int32
      download_dll("http://nemocas.org/binaries/w32-libflint.dll", joinpath(vdir, "lib", "libflint.dll"))
   else
      download_dll("http://nemocas.org/binaries/w64-libflint.dll", joinpath(vdir, "lib", "libflint.dll"))
   end
   try
      run(`ln -sf $vdir\\lib\\libflint.dll $vdir\\lib\\libflint-13.dll`)
   catch
      cp(joinpath(vdir, "lib", "libflint.dll"), joinpath(vdir, "lib", "libflint-13.dll"), force=true)
   end
   println("DONE")
else
   println("Building flint ... ")
   cd(joinpath("$wdir", "flint2"))
   withenv("LD_LIBRARY_PATH"=>"$vdir/lib", "LDFLAGS"=>LDFLAGS) do
      run(`./configure --prefix=$vdir --disable-static --enable-shared --with-mpir=$vdir --with-mpfr=$vdir`)
      run(`make -j4`)
      run(`make install`)
   end
   println("DONE")
end

cd(wdir)

# INSTALL ARB

if !Sys.iswindows()
  println("Cloning arb ... ")
  try
    run(`git clone https://github.com/fredrik-johansson/arb.git`)
    cd(joinpath("$wdir", "arb"))
    run(`git checkout $ARB_VERSION`)
    cd(wdir)
  catch
    if ispath(joinpath("$wdir", "arb"))
      cd(joinpath("$wdir", "arb"))
      run(`git fetch`)
      run(`git checkout $ARB_VERSION`)
      cd(wdir)
    end
  end
  println("DONE")
end

cd(wdir)

if Sys.iswindows()
   println("Downloading arb ... ")
   if Int == Int32
      download_dll("http://nemocas.org/binaries/w32-libarb.dll", joinpath(vdir, "lib", "libarb.dll"))
   else
      download_dll("http://nemocas.org/binaries/w64-libarb.dll", joinpath(vdir, "lib", "libarb.dll"))
   end
   println("DONE")
else
   println("Building arb ... ")
   cd(joinpath("$wdir", "arb"))
   withenv("LD_LIBRARY_PATH"=>"$vdir/lib", "LDFLAGS"=>LDFLAGS) do
      run(`./configure --prefix=$vdir --disable-static --enable-shared --with-mpir=$vdir --with-mpfr=$vdir --with-flint=$vdir`)
      run(`make -j4`)
      run(`make install`)
   end
   println("DONE")
end

cd(wdir)

push!(Libdl.DL_LOAD_PATH, joinpath(vdir, "lib"))

cd(oldwdir)


text:
using Test
import MechanicalUnits: @import_expand, ∙, Quantity
@import_expand(m, s, N) # Takes long if defined late)
using MechGluecode
import Plots
using Plots: plot
import MechGluecode

@testset "Plots mixed quantity" begin
    s1x = [1.0,2,3]
    s2x = [1.0,2,3.5]
    s1y = [1.0,2,4]
    s2y = [1.0,2,2]
    mxqm = hcat(s1x∙m, s2x∙m²)
    myqm = hcat(s1y∙s, s2y∙s²)
    plot(mxqm, myqm, seriestype = [:path :scatter], label = ["path" "scatter"], xguide = "x", yguide = "y")
end

@testset "Plots functions with mixed units" begin
    f_q_q(t::Quantity) = sin(0.3∙2π∙t / s) ∙ 9.81N
    s_1_down  = range(20, 0, length = 100)
    s_1_up = range(0, 20, length = 100)
    s_2_down = range(8, 1, length = 100)
    s_2_up =  range(1, 8, length = 100)
    plot( xlims =(-5s, 5s),
        [f_q_q  x-> 0.5f_q_q(1.5x)∙m],
        ribbon = ([s_1_down∙N  s_2_down∙N∙m], [s_1_up∙N   s_2_up∙N∙m]),
        yguide = "Load", xguide = "Time")
end


text:
module CircuitModelZoo

using LazyArtifacts

export zoo_file, zoo_clt_file, zoo_cnf_file, zoo_dnf_file, zoo_jlc_file, zoo_nnf_file, zoo_sdd_file, zoo_vtree_file, zoo_jpc_file, zoo_psdd_file, zoo_spn_file

const ZOO_VERSION = "/Circuit-Model-Zoo-0.1.6"

zoo_file(folder, name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/$folder/$name"

zoo_clt_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/clts/$name"

zoo_cnf_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/cnfs/$name"

zoo_dnf_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/dnfs/$name"

zoo_sdd_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/sdds/$name"

zoo_vtree_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/vtrees/$name"

zoo_jpc_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/jpcs/$name"

zoo_psdd_file(name) =
    artifact"circuit_model_zoo" * ZOO_VERSION * "/psdds/$name"

end # module