text:
# by lazarusA
using GLMakie
using GeometryBasics, LinearAlgebra, StatsBase
using Makie: get_dim, surface_normals
ϵ = 1e-7
# thanks to @jkrumbiegel for the lift
let
    function getMesh(x,y,z)
        positions = vec(map(CartesianIndices(z)) do i
        GeometryBasics.Point{3, Float32}(
            get_dim(x, i, 1, size(z)),
            get_dim(y, i, 2, size(z)),
            z[i])
        end)
        faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
        normals = surface_normals(x, y, z)
        vertices = GeometryBasics.meta(positions; normals=normals)
        meshObj = GeometryBasics.Mesh(vertices, faces)
        meshObj
    end
    function Π(v, n)
        v - dot(v, n) / norm(n)^2 * n
    end
    function multiply(q₁, q₂)
        z₁ = Complex(q₁[1], q₁[2])
        z₂ = Complex(q₁[3], q₁[4])
        A = [z₁ -conj(z₂); z₂ conj(z₁)]
        w₁ = Complex(q₂[1], q₂[2])
        w₂ = Complex(q₂[3], q₂[4])
        B = [w₁ -conj(w₂); w₂ conj(w₁)]
        C = A * B
        x₁, x₂ = C[1,1], C[2,1]
        [real(x₁); imag(x₁); real(x₂); imag(x₂)]
    end
    d(Φ, p) = begin
        rows = length(p)
        D = Array{Float64,1}(undef, rows)
        E = Array{Float64,2}(I, rows, rows) .* ϵ
        P = repeat(p, 1, rows)
        P′ = P + E
        for i in 1:rows
            D[i] = (Φ(P′[:, i]) - Φ(P[:, i])) ./ ϵ
        end
        D
    end
    directional(Φ, p, x) = begin
        rows = length(x)
        D = Array{Float64,1}(undef, rows)
        E = Array{Float64,2}(I, rows, rows) .* ϵ
        P = repeat(x, 1, rows)
        P′ = P + E
        for i in 1:rows
            D[i] = (Φ(p, P′[:, i]) - Φ(p, P[:, i])) ./ ϵ
        end
        D
    end
    getz₀(longitude, latitude, phase) = ℯ^(im * (longitude + phase)) * √((1 - sin(latitude)) / 2)
    getz₁(longitude, latitude, phase) = ℯ^(im * phase) * √((1 + sin(latitude)) / 2)
    TS³(z₀, z₁) =  begin
        # initialization
        X₀ = Complex(rand() + im * rand())
        X₁ = Complex(rand() + im * rand())
        getloss(z, x) = begin
            z₀ = Complex(z[1] + im * z[2])
            z₁ = Complex(z[3] + im * z[4])
            X₀ = Complex(x[1] + im * x[2])
            X₁ = Complex(x[3] + im * x[4])
            abs(conj(z₀) * X₀ + conj(z₁) * X₁)
        end
        x = [real(X₀); imag(X₀); real(X₁); imag(X₁)]
        z = [real(z₀); imag(z₀); real(z₁); imag(z₁)]
        loss = getloss(z, x)
        threshold = 1e-4
        η = 1e-2
        iterations = 150
        i = 0
        while loss > threshold && i < iterations
            g = directional(getloss, z, x)
            X₀ = Complex(x[1] - η * g[1] + im * (x[2] - η * g[2]))
            X₁ = Complex(x[3] - η * g[3] + im * (x[4] - η * g[4]))
            x = [real(X₀); imag(X₀); real(X₁); imag(X₁)]
            loss = getloss(z, x)
            i += 1
        end
        X₀, X₁
    end
    connection1form(z₀, z₁) = begin
        X₀, X₁ = TS³(z₀, z₁)
        1/2 * (conj(z₀) * X₀ - z₀ * conj(X₀) + conj(z₁) * X₁ - z₁ * conj(X₁))
    end
    tangent_tail(x, f) = [x[1]; x[2]; f(x)]
    tangent_head(x, f, ξ) = begin
        dΦ = d(f, x)
        ξΦ = dot(dΦ, ξ)
        dξΦ = dΦ .* ξ
        [dξΦ[1]; dξΦ[2]; (f([x[1] + dξΦ[1] * ϵ; x[2] + dξΦ[2] * ϵ]) - f(x)) ./ ϵ] .* ξΦ
    end
    f₁(i) = begin
        √(1 - (i[1])^2 - (i[2])^2)
    end
    f₂(i) = -f₁(i)
    getpoints(start, stop, segments1, segments2, point, radius) = begin
        getp(point) = begin
            ψ, ϕ, θ = point
            z₀ = getz₀(ϕ, θ, ψ)
            z₁ = getz₁(ϕ, θ, ψ)
            q₁, q₂ = real(z₀), imag(z₀)
            q₃, q₄ = real(z₁), imag(z₁)
            q = [q₁; q₂; q₃; q₄]
            p = q[1:3] ./ (1 - q₄)
            magnitude = norm(p)
            normalize(p) * tanh(magnitude)
        end
        getlengths(points) = begin
            lengths = []
            segments = length(points)
            for i in 1:segments-1
                p₁ = getp(points[i])
                p₂ = getp(points[i + 1])
                distance = norm(p₁ - p₂)
                push!(lengths, distance)
            end
            lengths
        end
        getcircumference(points) = begin
            circumference = 0
            segments = length(points)
            for i in 1:segments-1
                p₁ = getp(points[i])
                p₂ = getp(points[i + 1])
                distance = norm(p₁ - p₂)
                circumference += distance
            end
            circumference
        end
        getdistances(points, base) = begin
            distances = []
            segments1, segments2 = size(points)
            for i in 1:segments1
                for j in 1:segments2
                    p₁ = getp(points[i, j])
                    p₂ = getp([points[i, j][1]; base])
                    distance = norm(p₁ - p₂)
                    push!(distances, distance)
                end
            end
            distances
        end
        inputs = Array{Array{Float64,1},2}(undef, segments1, segments2)
        #ϕ, θ = point
        ϕ, θ = [point[1], -π/2]
        linearspace1 = collect(range(start, stop = stop, length = segments1))
        linearspace2 = collect(range(0, stop = 2π, length = segments2))
        for i in 1:segments1
            for j in 1:segments2
                ψ = linearspace1[i]
                ϕ = point[1] + cos(linearspace2[j]) * radius
                θ = point[2] + sin(linearspace2[j]) * radius
                inputs[i, j] = [ψ; ϕ; θ]
            end
        end
        getvariance(X) = begin
            N = length(X)
            μ = sum(X) / N
            σ² = sum((X .- μ).^2) / N
            σ²
        end
        getvariance1(X) = begin
            N = length(X)
            μ = sum(X) / N
            σ² = sum(exp.((X .- μ).^2) .- 1) / N
            σ²
        end
        converttovector(points) = begin
            vector = []
            segments1, segments2 = size(points)
            for i in 1:segments1
                for j in 1:segments2
                    vector = [vector; points[i, j]]
                end
            end
            vector
        end
        converttomatrix(vector, shape) = begin
            matrix = Array{Array{Float64,1},2}(undef, shape...)
            index = 1
            for i in 1:shape[1]
                for j in 1:shape[2]
                    matrix[i, j] = vector[index: index + 2]
                    index += 3
                end
            end
            matrix
        end
        getloss(points, base) = begin
            segments1, segments2 = size(points)
            loss = 0
            # get the length of segments along a fiber
            lengths = [getlengths(points[:, i]) for i in 1:segments2]
            loss += sum(getvariance.(lengths)) / segments2
            # get the length of segments along the sectional cut of a fiber
            lengths = [getlengths(points[i, :]) for i in 1:segments1]
            loss += sum(getvariance.(lengths)) / segments1
            # get the circumference of the sectional cut of a fiber
            circumferences = [getcircumference(points[i, :]) for i in 1:segments1]
            loss += getvariance(circumferences) / segments1
            # get the the distance between the central point of a sectional cut of a fiber and the surrounding star-shaped points
            distances = getdistances(points, base)
            loss += sum(getvariance.(distances)) / (segments1 * segments2)
            error1 = sum((points[begin, :][begin] .- start).^2) / segments2
            error2 = sum((points[end, :][begin] .- stop).^2) / segments2
            loss + error1 + error2
        end
        getloss1(v₁, v₂) = getloss(converttomatrix(v₁ + v₂, size(inputs)), point)
        v₂ = converttovector(inputs)
        v₁ = zeros(size(v₂))
        loss = getloss1(v₁, v₂)
        threshold = 1e-3
        η = 1e-3
        iterations = 300
        α = 0.9
        i = 1
        while loss > threshold && i < iterations
            heat = rand(size(v₂)...)
            v₂ = v₂ + (0.001 .* heat)
            g = directional(getloss1, v₁, v₂)
            new = v₂ - (η .* g)
            v₂ = α .* v₂ + (1 - α) .* new
            loss = getloss1(v₁, v₂)
            i += 1
        end
        outputs = converttomatrix(v₂, size(inputs))
        [getp(outputs[i, j]) for i in 1:segments1, j in 1:segments2]
    end

    segments = 30
    segments2 = 30
    radius = 0.05
    # Grids of polar and azimuthal angles
    θ = LinRange(0, π, segments)
    ϕ = LinRange(0, 2π, segments)
    xx = [sin(θ)*sin(ϕ) for θ in θ, ϕ in ϕ]
    yy = [sin(θ)*cos(ϕ) for θ in θ, ϕ in ϕ]
    zz = [cos(θ)        for θ in θ, ϕ in ϕ]

    l = Node(4)
    m = Node(1)

    lon = Node(1.0)
    lat = Node(1.5)
    ang = Node(1.0)
    vector = Node([1.0; 0.0; 0.0])

    #w₀ = @lift(ℯ^(im * $ang) * sqrt((1 + sin($lat)) / 2))
    #w₁ = @lift(ℯ^(im * ($lon + $ang)) * sqrt((1 - sin($lat)) / 2))
    #w = @lift($w₁ / $w₀)
    w = lift(lon, lat) do lon, lat
        a₁ = cos(lat) * cos(lon)
        a₂ = cos(lat) * sin(lon)
        a₃ = sin(lat)
        a₁ / (1 - a₃) + im * a₂ / (1 - a₃)
    end
    fiber = lift(lon, lat, ang) do lon, lat, ang
        #points = [[lon; lat] + [cos(i); sin(i)] .* radius for i in range(0, stop = 2π, length = segments2)]
        start = 0
        stop = ang
        #ψ = range(0, stop = ang, length = segments)
        #actions = [gets1action(start, stop, segments, point) for point in points]
        xyz(i, j) = begin
            point = points[j]
            action = actions[j][i]
            v₀ = getz₀(point[1], point[2], action)
            v₁ = getz₁(point[1], point[2], action)
            q₁, q₂ = real(v₀), imag(v₀)
            q₃, q₄ = real(v₁), imag(v₁)
            q = [q₁; q₂; q₃; q₄]
            p = q[1:3] ./ (1 - q₄)
            magnitude = norm(p)
            normalize(p) * tanh(magnitude)
        end
        #a = [xyz(i, j) for i in 1:segments, j in 1:length(points)]
        point = [lon; lat]
        a = getpoints(start, stop, segments, segments2, point, radius)
        a₁ = map(i -> i[1], a)
        a₂ = map(i -> i[2], a)
        a₃ = map(i -> i[3], a)
        a₁, a₂, a₃
    end
    fiber₁ = @lift($fiber[1])
    fiber₂ = @lift($fiber[2])
    fiber₃ = @lift($fiber[3])

    fiber2 = lift(lon, lat, ang) do lon, lat, ang
        #points = [[lon; lat] + [cos(i); sin(i)] .* radius for i in range(0, stop = 2π, length = segments2)]
        start = ang
        stop = 2π
        #ψ = range(start, stop = stop, length = segments)
        #actions = [gets1action(start, stop, segments, point) for point in points]
        xyz(i, j) = begin
            point = points[j]
            action = actions[j][i]
            v₀ = getz₀(point[1], point[2], action)
            v₁ = getz₁(point[1], point[2], action)
            q₁, q₂ = real(v₀), imag(v₀)
            q₃, q₄ = real(v₁), imag(v₁)
            q = [q₁; q₂; q₃; q₄]
            p = q[1:3] ./ (1 - q₄)
            magnitude = norm(p)
            normalize(p) * tanh(magnitude)
        end
        #a = [xyz(i, j) for i in 1:segments, j in 1:length(points)]
        point = [lon; lat]
        a = getpoints(start, stop, segments, segments2, point, radius)
        a₁ = map(i -> i[1], a)
        a₂ = map(i -> i[2], a)
        a₃ = map(i -> i[3], a)
        a₁, a₂, a₃
    end
    fiber2₁ = @lift($fiber2[1])
    fiber2₂ = @lift($fiber2[2])
    fiber2₃ = @lift($fiber2[3])

    tail = lift(w) do w
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        tangent_tail([u; v], f)
    end
    xconst = lift(w, tail) do w, tail
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [1; 0]))
    end
    yconst = lift(w, tail) do w, tail
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [0; 1]))
    end
    head = lift(w, ang, tail, xconst, yconst) do w, ang, tail, x, y
        # parallel transport
        #v = to_value(vector)
        #n = normalize(cross(x, y))
        #if -1 < dot(x, y) < 1
        #    projected = Π(v, n)
        #else
        #    projected = v
        #end
        #projected = normalize(projected)
        #vector[] = projected
        #projected
        u, v = real(w), imag(w)
        denominator = u^2 + v^2 + 1
        u = 2u / denominator
        v = 2v / denominator
        f = abs(w) ≤ 1 ? f₂ : f₁
        normalize(tangent_head([u; v], f, [cos(ang); sin(ang)]))

        #p₁ = [lon; lat]
        #p₂ = p₁ + [0.1; 0.1]
        #v₀ = ℯ^(im * (p₁[1] + ψ)) * sqrt((1 - sin(p₁[2])) / 2)
        #v₁ = ℯ^(im * ψ) * sqrt((1 + sin(p₁[2])) / 2)

    end
    connectiontail = lift(lon, lat, ang) do lon, lat, ang
        point = [lon; lat]
        v₀ = getz₀(point[1], point[2], ang)
        v₁ = getz₁(point[1], point[2], ang)
        q₁, q₂ = real(v₀), imag(v₀)
        q₃, q₄ = real(v₁), imag(v₁)
        q = [q₁; q₂; q₃; q₄]
        p = q[1:3] ./ (1 - q₄)
        magnitude = norm(p)
        normalize(p) * tanh(magnitude)
    end
    connectionhead = lift(lon, lat, ang) do lon, lat, ang
        point = [lon; lat]
        v₀ = getz₀(point[1], point[2], ang)
        v₁ = getz₁(point[1], point[2], ang)
        oneform = connection1form(v₀, v₁)
        v₀, v₁ = oneform * v₀, oneform * v₁
        q₁, q₂ = real(v₀), imag(v₀)
        q₃, q₄ = real(v₁), imag(v₁)
        q = [q₁; q₂; q₃; q₄]
        p = q[1:3] ./ (1 - q₄)
        magnitude = norm(p)
        h = normalize(p) * tanh(magnitude)
        t = to_value(connectiontail)
        h - t
    end

    tails₁ = @lift([$tail[1], $tail[1], $tail[1]])
    tails₂ = @lift([$tail[2], $tail[2], $tail[2]])
    tails₃ = @lift([$tail[3], $tail[3], $tail[3]])
    heads₁ = @lift([$head[1], $xconst[1], $yconst[1]])
    heads₂ = @lift([$head[2], $xconst[2], $yconst[2]])
    heads₃ = @lift([$head[3], $xconst[3], $yconst[3]])

    conn_tail₁ = @lift([$connectiontail[1]])
    conn_tail₂ = @lift([$connectiontail[2]])
    conn_tail₃ = @lift([$connectiontail[3]])
    conn_head₁ = @lift([$connectionhead[1]])
    conn_head₂ = @lift([$connectionhead[2]])
    conn_head₃ = @lift([$connectionhead[3]])

    x₁ = @lift(xx .* 0.1 .+ $tail[1])
    y₁ = @lift(yy .* 0.1 .+ $tail[2])
    z₁ = @lift(zz .* 0.1 .+ $tail[3])

    ambient =  Vec3f0(0.75, 0.75, 0.75)
    with_theme(theme_black()) do
        fig = Figure(resolution = (1400, 800))

        ax1 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        ax2 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        pltobj1 = mesh!(ax1, getMesh(xx, yy, zz), color = :blue, ambient = ambient)
        pltobj3 = mesh!(ax1, @lift(getMesh($x₁, $y₁, $z₁)), color = :gold, ambient = ambient)
        pltobj4 = arrows!(ax1, tails₁, tails₂, tails₃, heads₁, heads₂, heads₃, color = [:red, :silver, :gray], lengthscale = 0.3f0)
        pltobj5 = mesh!(ax2, @lift(getMesh($fiber₁, $fiber₂, $fiber₃)), color = :red, ambient = ambient)
        pltobj6 = mesh!(ax2, @lift(getMesh($fiber2₁, $fiber2₂, $fiber2₃)), color = :blue, ambient = ambient)
        pltobj7 = arrows!(ax2, conn_tail₁, conn_tail₂, conn_tail₃, conn_head₁, conn_head₂, conn_head₃, color = [:gold], lengthscale = 1.0f0)
        fig[1,1] = ax1
        fig[1,2] = ax2
        round3(n) = round(n, sigdigits = 3)
        degree(n) = round3(180(n / π))
        fig[0,1:2] = Label(fig, @lift("Hopf fiber: longitude = $(degree($lon)), latitude = $(degree($lat)), S¹action = $(degree($ang))."), textsize = 20)
        sl3 = Slider(fig[end+1, 1:2], range = range(-π, stop = π, length = 100))
        sl4 = Slider(fig[end+1, 1:2], range = range(-π/2, stop = 0.99 * π/2, length = 100))
        sl5 = Slider(fig[end+1, 1:2], range = range(0, stop = 2π, length = 100))
        set_close_to!(sl3, to_value(lon))
        set_close_to!(sl4, to_value(lat))
        set_close_to!(sl5, to_value(ang))
        connect!(lon, sl3.value)
        connect!(lat, sl4.value)
        connect!(ang, sl5.value)
        display(fig)
        framerate = 30
        totaltime = 30
        timestamps = range(0, totaltime, step=1/framerate)
#=
        record(fig, joinpath("gallery", "connection1form_a.mp4"), timestamps;
            framerate = framerate) do t
            step = t / totaltime
            longitude = π/2
            latitude = π/4
            action = step * 2π
            #set_close_to!(sl3, longitude)
            #set_close_to!(sl4, latitude)
            set_close_to!(sl5, action)
        end =#
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module AsapMapper

using Mapper2
using IterTools
using JSON
using Logging
using DataStructures

import Base: parse
# Set up directory paths
const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)

#set_logging(level) = configure_logging(AsapMapper, min_level=level)
set_logging(level) = nothing

export  place_and_route,
        testmap,
        asap_place,
        parallel_map_and_save,
        place!,
        route!,
        set_logging,
        # Taskgraph constructors
        load_taskgraph,
        # Architecture constructors
        asap4,
        asap3,
        asap2,
        Rectangular,
        Hexagonal,
        # Architecture types
        KC, Asap2,
        # Project Manger interface
        PMConstructor,
        SimConstructor,
        build_map

# Helper Functions
function oneofin(a,b)
    for i in a
        i in b && return true
    end
    return false
end

const push_to_dict = Mapper2.Helper.push_to_dict

################################################################################
# Attributes to determine what tasks may be mapped to which components in
# the architecture.
################################################################################
abstract type MapConstructor end

################################################################################
# Custom Architecture used by this Framework
################################################################################

include("Mapper2_Interface.jl")
include("Helper.jl")
include("Metadata.jl")

# Architectures
include("cad_models/cad_models.jl")
#include("experimental_models/experimental_models.jl")
#include("models/models.jl")

include("PM_Interface/PM_Interface.jl")
include("Simulator_Interface.jl")

# Customize placement/routing plus architectures.
include("PNR.jl")

#include("IP_Router/Router.jl")
include("Plots/MappingPlots.jl")

################################################################################
# Generic Place and Route function.
################################################################################

function place_and_route(profile_path, dump_path)
    # swoop(profile_path)
    # Initialize an uncompressed taskgraph constructor
    c = PMConstructor(profile_path)
    m = build_map(c)

    # Run place-and-route
    if typeof(m.options[:existing_map]) <: Nothing
        m = asap_pnr(m)
    end
    # Dump mapping to given dump path
    dump_map(m, dump_path)
    # Return the mapping object for possible later plotting
    return m
end

function parallel_map_and_save(input_file, output_dir, num_mappings)
    # Make output directory if it does not exist.
    if !isdir(output_dir)
        mkdir(output_dir)
    end

    # Launch a parallel mapping
    c = PMConstructor(input_file)
    maps = pmap((i) -> (asap_pnr ∘ build_map)(c), 1:num_mappings)

    for (i,m) in enumerate(maps)
        savepath = joinpath(output_dir, "map_$(i).jls")
        Mapper2.MapperCore.save(m, savepath)
    end
end

report_routing_stats(m::Map) = Mapper2.MapperCore.report_routing_stats(m)

function approximate_energy(m::Map)
    taskgraph = m.taskgraph
    total = 0
    for (index, edge) in enumerate(getedges(taskgraph))
        # Get the number of reads
        measurements = edge.metadata["measurements_dict"]
        reads = get(measurements, "num_reads", nothing)
        isnothing(reads) && continue

        # Get the total number of links used for this edge
        links = Mapper2.MapperCore.count_global_links(getpath(m, index))
        total += reads * links
    end
    return total
end

end # module

ext: jl


text:
using MusicTransformer
using NoteSequences
using Test

include("musicencoders.jl")


text:
"""
    module PermPlain

Permuation types and operations.

Types representing permutations: `PList`, `PCycles`, `PMatrix`, `PDict`

Functions for converting among reprentations: `permlist`, `permcycles`, `permmatrix`, `permsparse`.

Functions: `isperm`, `permsgn`, `canoncycles`, `cycletype`, `cyclelengths`, `numcycles`, `permlength`,
      `permorder`, `permpower`, `permpower2`, `permcommute`, `permcompose`, `permcompose!`,
     `permkron`, `ltpermlist`, `support`, `supportsize`, `greatestmoved`, `leastmoved`, `preimage`, `same`, `fixed`
"""
module PermPlain

using DataStructures: counter
import LinearAlgebra
using SparseArrays: sparse
using Random
import Base: isperm, ==

include("collect.jl")

export permlist, permcycles, permsparse, permmatrix, canoncycles, permpower, permpower2,
    permsgn, cycletype, support, supportsize, greatestmoved, leastmoved, fixed
export cyclelengths, numcycles, permlength, permorder, getdata, cyc_pow_perm, permcompose, permcompose!,
    permcommute, permkron, ltpermlist, preimage
export PCycles, PList, PMatrix, PDict

randperm(::Type{T}, n::Integer) where T = collect(T, randperm(n))

eye(::Type{T}, n::Integer) where T =  Matrix{T}(LinearAlgebra.I, n, n)
speye(::Type{T}, n) where T = sparse(LinearAlgebra.I * one(T), n, n)

PCyclesVec{T} = AbstractVector{<:AbstractVector{T}} where {T <: Real}

"""
    PCycles{T}

Cycle representation of a permutation.
"""
mutable struct PCycles{T}
    data::PCyclesVec{T}
    len::Int
end

"""
    permlength(p)

The number of elements in the permutation.
"""
permlength(c::PCycles) = c.len
getdata(c::PCycles) = c.data
==(c1::PCycles, c2::PCycles) = c1.data == c2.data && c1.len == c2.len

"""
    PList{T} = AbstractVector{T} where {T <: Real}

List representation of a permtuation.
"""
PList{T} = AbstractVector{T} where {T <: Real}

permlength(p::PList) = length(p)

"""
    PMatrix{T} = AbstractMatrix{T} where {T <: Real}

Permutation matrix.
"""
PMatrix{T} = AbstractMatrix{T} where {T <: Real}
permlength(m::PMatrix) = size(m, 1)

PDictData{T} = Dict{T, T} where T

"""
    PDict{T}

Sparse representation of a permtuation as a `Dict`.
"""
mutable struct PDict{T}
    data::PDictData{T}
    len::Int
end

permlength(c::PDict) = c.len
getdata(c::PDict) = c.data
==(c1::PDict, c2::PDict) = c1.data == c2.data && c1.len == c2.len

# Compute cyclic decomposition.
# Builds a cycle list in the canonical order.
# See pari for possible efficiency improvement.
"""
    permcycles(p::AbstractVector{T}) where {T <: Real}

Return a cycle representation of the permuation `p`.

### Example
```jldoctest
julia> p = [10, 8, 5, 1, 9, 4, 6, 3, 2, 7];

julia> permcycles(p)
2-element Array{Array{Int64,1},1}:
 [1, 10, 7, 6, 4]
 [2, 8, 3, 5, 9]
```
"""
function permcycles(p::PList{T}) where T
    n = length(p)
    visited = falses(n)
    cycles = Vector{Vector{T}}(undef, 0)
@inbounds for k in convert(T,1):convert(T,n)
        if ! visited[k]
            knext = k
            cycle = Array{T}(undef, 0)
            while ! visited[knext]
                push!(cycle,knext)
                visited[knext] = true
                knext = p[knext]
            end
            length(cycle) > 1 && push!(cycles,cycle)
        end
    end
    return PCycles(cycles, n)
end

"""
    permcycles(m::AbstractMatrix)

Return a cycle representation of the permutation `m`.

### Example
```jldoctest
julia> p = [6, 1, 5, 3, 4, 2];

julia> m = permmatrix(p)
6×6 Array{Int64,2}:
 0  0  0  0  0  1
 1  0  0  0  0  0
 0  0  0  0  1  0
 0  0  1  0  0  0
 0  0  0  1  0  0
 0  1  0  0  0  0

julia> permcycles(m)
2-element Array{Array{Int64,1},1}:
 [1, 6, 2]
 [3, 5, 4]
```
"""
permcycles(m::PMatrix) = permcycles(permlist(m))

"""
    permcycles(sp::Dict{T, T})

Return a cycle representation of the permutation `sp`.

### Example
```jldoctest
julia> p = [6, 1, 5, 3, 4, 2];

julia> sp = permsparse(p)
(Dict(4=>3,2=>1,3=>5,5=>4,6=>2,1=>6), 6)

julia> p = [10, 8, 5, 1, 9, 4, 6, 3, 2, 7];

julia> sp, maxn = permsparse(p)
(Dict(7=>6,9=>2,4=>1,10=>7,2=>8,3=>5,8=>3,5=>9,6=>4,1=>10…), 10)

julia> permcycles(sp)
2-element Array{Array{Int64,1},1}:
 [7, 6, 4, 1, 10]
 [9, 2, 8, 3, 5]
```
"""
permcycles(sp::PDict) = sparsetocycles(sp)

function sparsetocycles(dsp::PDict{T}) where T
    cycs = Array{Vector{T}}(undef, 0)
    permlength(dsp) == 0 && return PCycles(cycs, permlength(dsp))
    sp = getdata(dsp)
    ks = collect(keys(sp))
    n = length(ks)
    seen = Dict{T, Bool}()
    for k in ks seen[k] = false end
@inbounds   k = ks[1]
    nseen = 0
@inbounds while nseen <= n
        foundunseen = false
        for i in ks  # could be more efficient starting search after last found
            if seen[i] == false
                k = i
                foundunseen = true
                break
            end
        end
        foundunseen == false && break
        kcyclestart = k
        cyc = Array{T}(undef, 0)
        while true
            push!(cyc,k)
            if seen[k] == true
                error("Algorithm error: double setting seen k=$k, nseen=$nseen, kcyclestart=$kcyclestart")
            end
            seen[k] = true
            k = sp[k]
            nseen = nseen + 1
            if k == kcyclestart
                break
            end
        end
        push!(cycs,cyc)
    end
    return PCycles(cycs, permlength(dsp))
end

"""
    permlist(m::PMatrix)

Return a list representation of the permutation `m`.

### Example
```jldoctest
julia> m = [0 0 0 0 0 1; 1 0 0 0 0 0; 0 0 0 0 1 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 1 0 0 0 0]
6×6 Array{Int64,2}:
 0  0  0  0  0  1
 1  0  0  0  0  0
 0  0  0  0  1  0
 0  0  1  0  0  0
 0  0  0  1  0  0
 0  1  0  0  0  0

julia> print(permlist(m))
[6, 1, 5, 3, 4, 2]
```
"""
permlist(m::PMatrix) = mattoperm(m)
mattoperm(m::PMatrix{T}) where T = mattoperm!(Vector{T}(undef, size(m, 1)), m)

function mattoperm!(p, m::PMatrix{T}) where T
    n = size(m)[1]
    maxk = zero(T)
@inbounds  for i in 1:n
        for j in 1:n
            if m[j,i] != 1
                continue
            end
            p[j] = i
        end
    end
    p
end

"""
    permlist(cycs::PCycles, pmax::Real = 0)

Return the list representation of the permuatation `cycs`, where `pmax` is
the length of the permutation.

### Examples
```jldoctest
julia> pcyc = [[4, 3, 5], [2, 1, 6]]
2-element Array{Array{Int64,1},1}:
 [4, 3, 5]
 [2, 1, 6]

julia> print(permlist(pcyc))
[6, 1, 5, 3, 4, 2]
```
"""
permlist(cycs::PCycles) = _permlist(getdata(cycs), permlength(cycs))

_permlist(cycsvec::PCyclesVec, pmax::Real = 0) = cycstoperm(cycsvec, pmax)

function cycstoperm(cycsvec::PCyclesVec{T}, pmax::Integer = 0) where T
    isempty(cycsvec) && return collect(one(T):convert(T,pmax))
    cmaxes = [maximum(c) for c in cycsvec]
    cmax = maximum(cmaxes)  # must be a faster way
#    perm = [one(T): (pmax > cmax ? convert(T,pmax) : convert(T,cmax))]
    perm = collect(one(T): (pmax > cmax ? convert(T,pmax) : convert(T,cmax)))
@inbounds  for c in cycsvec
        for i in 2:length(c)
            perm[c[i-1]] = c[i]
        end
        perm[c[end]] = c[1]
    end
    perm[cmax+1:pmax] = collect(cmax+1:pmax)
    return perm
end

"""
    permlist(sp::PDict)

Return the list representation of the permuatation `cycs`.

### Examples
```jldoctest
julia> sp = Dict(4=>3,2=>1,3=>5,5=>4,6=>2,1=>6);

julia> print(permlist(sp))
[6, 1, 5, 3, 4, 2]
```
"""
permlist(sp::PDict) = sparsetolist(sp)

function sparsetolist(sp::PDict{T}) where T
    p = collect(one(T):convert(T, permlength(sp)))
    data = getdata(sp)
@inbounds for (i,v) in data
        p[i] = v
    end
    return p
end

## permsparse ##
"""
    permsparse(m::PMatrix)

Return a sparse (`PDict`) representation of the permutation `m`.
"""
permsparse(m::PMatrix) = mattosparse(m)

function mattosparse(m::PMatrix{T}) where T
    p = Dict{T, T}()
    return PDict{T}(mattoperm!(p, m), permlength(m))
end

permsparse(p::PList) = listtosparse(p)

function listtosparse(p::PList{T}) where T
    data = Dict{T,T}()
    maxk = zero(T)
    length(p) == 0 && return PDict{T}(data, maxk)
@inbounds for i in p
        pv = p[i]
        if pv != i
            data[i] = pv
            pv > maxk ? maxk = pv : nothing
        end
    end
    return PDict(data, permlength(p))
end

permsparse(cycs::PCycles) = cycstosparse(cycs)

function cycstosparse(cycs::PCycles{T}) where T
    data = Dict{T,T}()
    maxk = zero(T)
@inbounds for c in getdata(cycs)
        pv = c[1]
        data[c[end]] = pv
        c[end] > maxk ? maxk = c[end] : nothing
        for i in 1:length(c)-1
            pv = c[i]
            data[pv] = c[i+1]
            pv > maxk ? maxk = pv : nothing
        end
    end
    return PDict(data, permlength(cycs))
end

## permmatrix ##

"""
    permmatrix(p, sparse::Bool = false)

Convert `p` to a matrix representation.
"""
permmatrix(p::PList, sparse::Bool = false) = permtomat(p, sparse)
# Convert PLIST to PMAT
function permtomat(p::PList{T}, sparse::Bool = false) where T <: Real
    n::T = length(p)
    A = sparse ? speye(T, n) : eye(T, n)
    return A[p,:]
end

function permmatrix(sp::PDict{T}) where {T}
    n = convert(T, permlength(sp))
    ot = one(T)
    z = zero(T)
    m = eye(Int, n)
    data = getdata(sp)
@inbounds for (i, v) in data
        m[i,i] = z
        m[i,v] = ot
    end
    return m
end

# Convert cyclic decomposition to canonical form
# used by gap, Mma, and Arndt thesis.
# Note that Arndt uses a different order internally to store the cycles as a single list.
"""
    canoncycles(cycs::PCycles)

Convert cyclic decomposition `cycs` to a canonical form.
"""
function canoncycles(cycs::PCycles{T}) where T
    ocycsdata = Array{Array{T,1}}(undef, 0)
    for cyc in getdata(cycs)
        push!(ocycsdata, circshift(cyc, -argmin(cyc) + 1))
    end
    sort!(ocycsdata, by=(x)->x[1])
    return PCycles(ocycsdata, permlength(cycs))
end

## cyclelengths. Find cyclic decomposition, but only save cycle lengths ##

"""
    cyclelengths(p)

Return a `Vector` of the lengths of the cycles in the permutation `p`.

`cyclelengths(p)` is equivalent to, but more efficient than `lenth(permcycles(p))`.
"""
cyclelengths(c::PCycles) = [length(x) for x in getdata(c)]
cyclelengths(m::PMatrix) = cyclelengths(permlist(m))
function cyclelengths(p::PList{T}) where T
    n = length(p)
    visited = falses(n)
    lengths = Array{Int}(undef, 0)
@inbounds for k in one(T):convert(T,n)
        if ! visited[k]
            knext = k
            len = 0
            while ! visited[knext]
                len += 1
                visited[knext] = true
                knext = p[knext]
            end
            len > 1 && push!(lengths,len)
        end
    end
    return lengths
end

# Gives cyclelengths in canonical order.
function cyclelengths(spd::PDict{T}) where T
    sp = getdata(spd)
    cyclens = Array{Int}(undef, 0)
    isempty(sp) && return cyclens
    ks = sort(collect(keys(sp)))
    n = length(ks)
    seen = Dict{T, Bool}()
    for k in ks seen[k] = false end
    k = ks[1]
    nseen = 0
    while nseen <= n
        didsee = false
        for i in ks
            if seen[i] == false
                k = i
                didsee = true
                break
            end
        end
        didsee == false && break
        k1 = k
        nincyc = 0
        while true
            nincyc += 1
            if seen[k] == true
                error("Algorithm error: double setting seen k=$k, nseen=$nseen, k1=$k1")
            end
            seen[k] = true
            k = sp[k]
            nseen = nseen + 1
            if k == k1
                break
            end
        end
        push!(cyclens,nincyc)  # inefficient
    end
    return cyclens
end

"""
    numcycles(p)

Return the number of cycles in the permutation `p`.

`numcycles(p)` is in general more efficient than
`length(cyclelengths(p))`.
"""
numcycles(p) = length(cyclelengths(p))
numcycles(c::PCycles) = length(getdata(c))
numcycles(m::PMatrix) = numcycles(permlist(m))
function numcycles(p::PList{T}) where T
    n = length(p)
    visited = falses(n)
    Ncycles = 0
@inbounds for k in one(T):convert(T, n)
        if ! visited[k]
            knext = k
            len = 0
            while ! visited[knext]
                len += 1
                visited[knext] = true
                knext = p[knext]
            end
            len > 1 ? Ncycles += 1 : nothing
        end
    end
    return Ncycles
end

# Wikipedia says some authors also require no fixed points.
# Here, we allow fixed points
"""
    iscyclic(p)

Return `true` if permutation `p` is cyclic.
"""
iscyclic(p) = numcycles(p) == 1

## cycle type, sign ##

"""
    cycletype(p)

The cycle type of permutation `p`.
"""
cycletype(p) = counter(cyclelengths(p))

permsgn_from_lengths(lens) = (-1)^(length(lens)+sum(lens))

"""
    permsgn(p)

The sign of permutation `p`.
"""
permsgn(p) = permsgn_from_lengths(cyclelengths(p))

function permorder_from_lengths(clengths)
    result = 1
    for c in clengths
        result = lcm(result, c)
    end
    return result
end

"""
    permorder(p)

Return the order of permutation `p`.
"""
permorder(p) = permorder_from_lengths(cyclelengths(p))

# distance between two PLISTs
# could use a macro for this and ==.
# is there a penalty for using swap macro on p and q instead of two branches ?

"""
    permdistance(p::PList, q::PList)

Compute the distance between permutations `p` and `q`.
"""
function permdistance(p::PList{T}, q::PList{T}) where T
    lp = length(p)
    lq = length(q)
    count = 0
    if lp < lq
        for i in 1:lp
            p[i] != q[i] ? count += 1 : nothing
        end
        for i in lp+1:lq
            q[i] != i ? count += 1 : nothing
        end
    else  # could factor code with refs, prbly not worth the trouble
        for i in 1:lq
            p[i] != q[i] ? count += 1 : nothing
        end
        for i in lq+1:lp
            p[i] != i ? count += 1 : nothing
        end
    end
    return count
end

# composition (multiplication) of two PLISTs

"""
    permcompose(p::PList, q::PList)

Compute the composition of permutations `p` and `q`.
"""
function permcompose(p::PList{T}, q::PList) where {T}
    lq = length(q)
    lp = length(p)
    lq == lp && return p[q]  # Much faster to use the Base method
    lr = lq < lp ? lp : lq
    r = Array{T}(undef, lr)
    @inbounds  if lq <= lp
        r[1:lq] = p[q]
        r[lq+1:lp] = p[lq+1:lp]
    else
      @inbounds  for k in 1:lq
            i = q[k]
            r[k] = (i < lp ? p[i] : i)
        end
    end
    return r
end

function permcompose(q::PDict{T}, p::PDict{V}) where {T, V}
    dout = Dict{T,T}()
    z = zero(T)
    maxk = z
    seen = Dict{T,Bool}()
    for k in keys(q) seen[k] = false end  # wasteful!
    for (k,pofk) in p
        qofpofk = get(q,pofk,z)
        seen[pofk] = true
        k == qofpofk ? continue : nothing  # ignore 1-cycles
        dout[k] = (qofpofk == z ? pofk : qofpofk)
          qofpofk > maxk ? maxk = qofpofk : nothing
    end
    for j in keys(q)
        seen[j] ? nothing : dout[j] = q[j]
    end
    return PDict(dout, max(permlength(q), permlength(p)))
end

permcompose(p::PMatrix, q::PMatrix) = p * q

## permapply ##

# Vector ?
function permapply(q::PDict{T}, a::AbstractArray) where T
    aout = copy(a)
    len = length(aout)
    for (k,v) in getdata(q)
        if k <= len && v <= len
            aout[k] = a[v]
        end
    end
    return aout
end

function permapply(q::PList, a::AbstractArray)
    aout = copy(a)
    lenq = length(q)
    lena = length(a)
    len =  lenq < lena ? lenq : lena
    for k in 1:len
        v = q[k]
        v <= len  && (aout[k] = a[v])
    end
    return aout
end

permapply(q::Union{PDict{T}, PList{T}}, a::String) where T = String(permapply(q, a.data))

# power of PLIST. output is PLIST
# The method permpower below is usually preferred because it does less allocation
function permpower2(p::PList{T}, n::Integer) where T
    n == 0 && return [one(T):convert(T,length(p))]
    n == 1 && return copy(p) # for consistency, don't return ref
    n < 0  && return permpower2(invperm(p),-n)
    q = permpower2(p, Int(floor(n/2)))
    q = q[q]
    return iseven(n) ? q : p[q]
end

function permpower!(p::PList{T},
                    pret::PList{T},
                    ptmp::PList{T},
                    n::Integer) where T
    onep = one(T)
    lenp = convert(T,length(p))
    n == 0 && (for i in onep:lenp pret[i] = i end; return )
    n < 0  && (permpower!(invperm(p), pret, ptmp, -n); return )
    n == 1 && (copy!(pret,p); return)
    permpower!(p, ptmp, pret,Int(floor(n/2)))
    if iseven(n)
        for i in onep:lenp pret[i] = ptmp[ptmp[i]] end
    else
        for i in onep:lenp pret[i] = p[ptmp[ptmp[i]]] end
    end
end

# This does less allocation (in general) than permpower2, and is faster in benchmarks
function permpower(p::PList{T}, n::Integer) where T
    n == 0 && return [one(T):convert(T,length(p))]
    n == 1 && return copy(p) # for consistency, don't return ref
    pret = similar(p)
    ptmp = similar(p)
    permpower!(p,pret,ptmp,n)
    return pret
end

# Compute power of permutation. Both input and output are PCYC
# Translated from pari perm.c
# Careful of degeneracy, and empty array may be returned.
function permpower(pcyc::PCycles{T}, exp::Integer) where {T}
    cyc = getdata(pcyc)
    r = 1
    for j in 1:length(cyc)
        r += gcd(length(cyc[j])-1,exp)
    end
    c = Vector{Vector{T}}(undef, 0)
    for j in 1:length(cyc)
        v = cyc[j]
        n = length(v)
        e = mod(exp,n)
        g = gcd(n,e)
        m = div(n,g)
        if m == 1 continue end
        for i in 1:g
            p = Array{Int}(undef, 0)
            l = i
            for k in 1:m
                push!(p,v[l])
                l += e
                if l > n  l -= n  end
            end
            push!(c,p)
        end
    end
    return PCycles(c, permlength(pcyc))
end

permpower(q::PDict{T}, exp::Integer) where T = cycstosparse(permpower(sparsetocycles(q),exp))

# power of PCYC. output is PLIST
# see pari perm.c
function cyc_pow_perm(pcyc::PCycles{T}, exp::Integer) where {T}
    n = 0
    cyc = getdata(pcyc)
    # cmaxes = [maximum(c) for c in cyc]
    # n = maximum(cmaxes)
    n = permlength(pcyc)
# This routine assumes cycles of length 1 are included. Still don't understand it.
#    for j = 1:length(cyc)
#        n += length(cyc[j])+1
#    end
    p = collect(T, 1:n) # wasteful
    for j = 1:length(cyc)
        v = cyc[j]
        n = length(v)
        e = mod(exp,n)
        l = e
        for k in 1:n
            nv = v[l+1]
            ind = v[k]
            if (ind>length(p))
                error("ind is $ind length p is $(length(p))")
            end
            p[ind] = nv
            l += 1
            l == n ? l = 0 : nothing
        end
    end
    return p
end

"""
    permcommute(p, q)

Return `true` if permutations `p` and `q` commuted.
"""
function permcommute(p::PList, q::PList)
    length(q) < length(p) ? (p,q) = (q,p) : nothing
    for i in length(p)
        q[p[i]] == p[q[i]] || return false
    end
    for i in length(p)+1:length(q)
        q[i] == i || return false
    end
    return true
end

permcommute(p::PMatrix, q::PMatrix) = p * q == q * p

"""
    isperm(m::AbstractMatrix)

Return `true` if `m` is a permutation matrix.
"""
function isperm(m::PMatrix{T}) where T
    sx,sy = size(m)
    sx == sy || return false
    z = zero(T)
    o = one(T)
    seen = falses(sx)
    for i in 1:sx
        for j in 1:sx
            val = m[j,i]
            if val == o
                seen[j] == true && return false
                seen[j] = true
            elseif val != z
                return false
            end
        end
    end
    true
end

# is cyclic decomposition a permutation
function isperm(pcycs::PCycles)
    cycs = getdata(pcycs)
    isempty(cycs) && return true
    seen = counter(eltype(cycs[1])) # inefficient
    for c in cycs
        for i in c
            push!(seen,i) == 1 || return false
        end
    end
    return true
end

# is sparse representation a permutation
function isperm(dsp::PDict)
    sp = getdata(dsp)
    sort(collect(keys(sp))) == sort(collect(values(sp)))  # inefficient
end

# TODO: Could be isone ? No,..
function isid(p::PList)
    for i in 1:length(p)
        i == p[i] || return false
    end
    return true
end

function permlistisequal(p::PList{T}, q::PList) where T
    lp = length(p)
    lq = length(q)
    s = one(T)
    if lp < lq
        for i in s:lp
            p[i] == q[i] || return false
        end
        for i in lp+s:lq
            q[i] == i || return false
        end
    else  # could factor code with refs, prbly not worth the trouble
        for i in s:lq
            p[i] == q[i] || return false
        end
        for i in lq+s:lp
            p[i] == i || return false
        end
    end
    return true
end

# is p < q with lexicographical ordering ?
for (f,ret) in ((:ltpermlist, :false), (:lepermlist, :true))
    @eval begin
function ($f)(p::PList{T}, q::PList{V}) where {T, V}
    o = one(T)
    lp = length(p)
    lq = length(q)
    if ( lp < lq )
        for k in o:lp
            if q[k] != p[k]
                return p[k] < q[k] ? true : false
            end
        end
        for k in lp+o:lq
            q[k] != k && return false
        end
    else
        for k in o:lq
            if q[k] != p[k]
                return p[k] < q[k] ? true : false
            end
        end
        for k in lq+o:lp
            p[k] != k && return false
        end
    end
    return $ret
end
end
end

"""
    preimage(p, k::Int)

The preimage of `k` under permutation `p`.
"""
function preimage(p::PList, k::Int)
    k > length(p) && return k
    for i in 1:length(p)
        if p[i] == k
            return i
        end
    end
    return k  #  make preimage consistent with image
#    error("Can't find inverse image of $k.")
end

function preimage(p::PDict, i::Int)
    for (k,v) in getdata(p)
        v == i && return k
    end
    return i  #  make preimage consistent with image
#    error("Can't find inverse image of $k.")
end

# List of points mapped to same point by p and q
#function same(pin::PermList, qin::PermList)
"""
    same(p::PList, q::PList)

A list of points mapped to same point by p and q
"""
function same(p::PList, q::PList)
    lp = length(p)
    lq = length(q)
    d = Array{eltype(p)}(undef, 0)
    if lp < lq
        for i in 1:lp
            p[i] == q[i] ? push!(d,p[i]) : nothing
        end
        for i in lp+1:lq
            q[i] == i ? push!(d,q[i]) : nothing
        end
    else  # could factor code with refs, prbly not worth the trouble
        for i in 1:lq
            p[i] == q[i] ? push!(d,p[i]) : nothing
        end
        for i in lq+1:lp
            p[i] != i ? push!(d,p[i]) : nothing
        end
    end
    return d
end

# The return type depends on value of input. How to get around this ?
# There is no integer Inf.
# agrees with gap (except definition, use of Inf is different)
"""
    leastmoved(p::PList)

The least element in permutation `p` that is moved.
"""
function leastmoved(p::PList)
    lp = length(p)
    lm = lp+1
    for i in 1:lp
        k = p[i]
        k == i ? nothing :
           k < lm ? lm = k : nothing
    end
    return lm > lp ? Inf : lm
end

# agrees with gap
"""
    greatestmoved(p::PList)

The greatest element in permutation `p` that is moved.
"""
function greatestmoved(p::PList)
    lp = length(p)
    gm = 0
    for i in 1:lp
        k = p[i]
        k == i ? nothing :
           k > gm ? gm = k : nothing
    end
    return gm
end

"""
    supportsize(p::PList)

The number of elements in the support of permutation `p`.
"""
function supportsize(p::PList)
    lp = length(p)
    count = 0
    for i in 1:lp
        k = p[i]
        k != i ? count += 1 : nothing
    end
    count
end

"""
    support(p::PList)

The support of permutation `p`. That is, the elements
that are moved under the permutation.
"""
function support(p::PList)
    lp = length(p)
    mov = Array{eltype(p)}(undef, 0)
    for i in 1:lp
        k = p[i]
        k != i ? push!(mov,i) : nothing
    end
    return mov
end

"""
    fixed(p::PList)

The elements that are not moved under permutation `p`.

`fixed` is the complement of `support`.
"""
function fixed(p::PList)
    lp = length(p)
    fixedel = Array{eltype(p)}(undef, 0)
    for i in 1:lp
        k = p[i]
        k == i ? push!(fixedel,i) : nothing
    end
    return fixedel
end

## Output ##

function cycleprint(io::IO, v::Vector)
    len = length(v)
    print(io,"(")
    for i in 1:len
        print(io,v[i])
        i != len ? print(" ") : nothing
    end
    print(io,")")
end

function permarrprint(io::IO, v::Array)
    if isempty(v)
        print(io, "()")
    else
        pd = ndigits(maximum(v))
        len = length(v)
        print(io,"( ")
        for k in v
            print(io,rpad(k,pd,' ')," ")
        end
        print(io,")")
    end
end

cycleprint(v::Array) = cycleprint(STDOUT,v)
permarrprint(v::Array) = permarrprint(STDOUT,v)

# Copied from Base string.jl
function pprint_to_string(printfunc::Function, xs...)
    s = IOBuffer(Array{Uint8}(isa(xs[1],String) ? endof(xs[1]) : 0), true, true)
    truncate(s,0)
    for x in xs
        printfunc(s, x)
    end
    takebuf_string(s)
end

permarrstring(xs...) = pprint_to_string(permarrprint, xs...)
cyclestring(xs...) = pprint_to_string(cycleprint, x...)

include("matrixops.jl")
include("pivot.jl")

end # module Permplain


text:
@generated function struct_as_namedtuple(st)
  A = ( Expr(:(=), n, :(st.$n)) for n in fieldnames(st))
  Expr(:tuple, A...)
end

remaker_of(prob::T) where {T} = parameterless_type(T){isinplace(prob)}

# Define `remaker_of` for the types that does not (make sense to)
# implement `isinplace` trait:
for T in [
    NoiseProblem,
    SplitFunction,  # TODO: use isinplace path for type-stability
    TwoPointBVPFunction,
    ]
  @eval remaker_of(::$T) = $T
end

"""
    remake(thing; <keyword arguments>)

Re-construct `thing` with new field values specified by the keyword
arguments.
"""
function remake(thing; kwargs...)
  T = remaker_of(thing)
  T(; struct_as_namedtuple(thing)...,kwargs...)
end

function remake(thing::AbstractJumpProblem; kwargs...)
  parameterless_type(thing)(remake(thing.prob;kwargs...))
end


text:
#
# Client/Server common declarations. Is used in Client, Server
# and parent code.
# TODO: describe this module as a glue between Client and Server modules
#
# @author DeadbraiN
#
module Connection
  import Config.@if_debug
  import Helper
  import Event
  import FastApi

  export REQUEST_BIT
  export RESPONSE_BIT

  export EVENT_BEFORE_RESPONSE
  export EVENT_AFTER_REQUEST

  export Answer
  export Command
  export CMD_NO_FUNC

  export request
  export fastRequest
  export answer
  export fastAnswer
  export notEmpty
  #
  # First bit in fast requests, which shows if current request
  # is a request or response
  #
  const REQUEST_BIT       = UInt8(0b10000000)
  const RESPONSE_BIT      = UInt8(0b01111111)
  #
  # Name of the event, which is fired if answer from client's
  # request is obtained.
  #
  const EVENT_AFTER_REQUEST  = "after-response"
  #
  # Name of the event, which is fired if client sent us a command. If
  # this event fires, then specified command should be runned here - on
  # server side.
  #
  const EVENT_BEFORE_RESPONSE = "before-response"
  #
  # Stub if Command.fn is not set
  #
  const CMD_NO_FUNC       = 0

  #
  # This command is used for sending some data from client
  # to server. It contains from function (or command id) and
  # arguments. An answer will back in Answer type.
  #
  type Command
    fn::Int
    args::Array{Any, 1}
  end
  #
  # This typeis used for answering from Server to Clients.
  # data property contains Server response value.
  #
  type Answer
    id::Int
    data::Any
  end
  #
  # Makes request to client. This method is not blocking. It returns
  # just after the call. Answer will be obtained in run() method
  # async loop.
  # @param sock Client socket returned by accept() function
  # @param fn Remote function id
  # @param args Custom fn arguments
  # @return true - request was sent, false wasn't
  #
  function request(sock::Base.TCPSocket, fn::Int, args...)
    if !Helper.isopen(sock) return false end
    try
      serialize(sock, Command(fn, collect(args)))
    catch e
      Helper.warn("Connection.request(): $e")
      @if_debug showerror(STDOUT, e, catch_backtrace())
      close(sock)
      return false
    end

    true
  end
  #
  # Fast version of request function. It uses native read()/write()
  # function for data transfer. This function is non blocking one.
  # TODO: describe dataIndex and data
  # @param sock Client socket returned by accept() function
  # @param dataIndex Index in api array
  # @param args Data to transfer in request
  #
  function fastRequest(sock::Base.TCPSocket, dataIndex::UInt8, args...)
    if !Helper.isopen(sock) return false end
    try
      #
      # This index will be used on client side for correct
      # reading of (data) parameters. First bit is used as
      # request bit.
      #
      _fastWrite(sock, dataIndex | REQUEST_BIT, args)
    catch e
      Helper.warn("Connection.fastRequest(): $e")
      @if_debug showerror(STDOUT, e, catch_backtrace())
      close(sock)
      return false
    end

    true
  end
  #
  # Reads one command from socket and fires an event to main
  # code. After that, writes an answer into the socket
  # back.
  # @param sock Socket
  # @param obs Observer for firing an event to "parent" code
  # @param excFn Callback function for socket exceptions
  #
  function answer(sock::Base.TCPSocket, obs::Event.Observer, excFn::Function)
    local data::Any
    local typ::DataType

    try
      #
      # Right now, only two types of responses are supported:
      # answers after client request(Answer) and
      # server requests (Command).
      #
      data = deserialize(sock)
      typ  = typeof(data)
      if typ === Answer
        Event.fire(obs, EVENT_AFTER_REQUEST, data)
      elseif typ === Command
        local ans::Answer = Answer(CMD_NO_FUNC, nothing)
        Event.fire(obs, EVENT_BEFORE_RESPONSE, sock, data, ans)
        if notEmpty(ans) serialize(sock, ans) end
      end
    catch e
      @if_debug showerror(STDOUT, e, catch_backtrace())
      return excFn(sock, e)
    end

    true
  end
  #
  # This is fast version of server answer, because it uses native
  # read()/write() functions. It listens socket and reads header
  # byte, which consists of 8 bits. First bit is request/response
  # bit. Last 7 - data index value. Request/response bit means# that
  # this is request or response from client. For example: let's
  # imagin we have a request from server to client with a response
  # as an answer. In this case first bit in request from server to
  # client will be REQUEST_BIT, but back response will contain RESPONSE_BIT.
  # The same for client to server request. Input request will cantain
  # REQUEST_BIT, but response (from server to client) - RESPONSE_BIT.
  # @param sock Client's socket
  # @param obs Observer for firing an event to "parent" code
  # @param excFn Callback function for socket exceptions
  #
  function fastAnswer(sock::Base.TCPSocket, obs::Event.Observer, excFn::Function)
    local dataIndex::UInt8
    local isRequest::Bool

    try
      #
      # Right now, only two types of responses are supported:
      # answers after server request(REQUEST_BIT) and client
      # requests (RESPONSE_BIT).
      #
      isRequest = (dataIndex = read(sock, UInt8)) & REQUEST_BIT > 0
      dataIndex = dataIndex & RESPONSE_BIT
      if isRequest
        local ans::Answer = Answer(CMD_NO_FUNC, nothing)
        Event.fire(obs, EVENT_BEFORE_RESPONSE, sock, _fastRead(sock, dataIndex), ans)
        if notEmpty(ans) _fastWrite(sock, UInt8(ans.id) & RESPONSE_BIT, ans.data) end
      else
        Event.fire(obs, EVENT_AFTER_REQUEST, _fastRead(sock, dataIndex))
      end
    catch e
      @if_debug showerror(STDOUT, e, catch_backtrace())
      return excFn(sock, e)
    end

    true
  end
  #
  # Checks if answer is empty - wasn't modified in parent code.
  # @return {Bool} empty status
  #
  function notEmpty(ans::Answer)
    ans.id !== CMD_NO_FUNC || ans.data !== nothing
  end
  #
  # Reads data types from socket and returns them in array
  # @param sock socket
  # @param dataIndex Index of data in FastApi array
  # @return {Any} read data
  #
  function _fastRead(sock::Base.TCPSocket, dataIndex::UInt8)
    local types::Array{DataType, 1} = FastApi.getTypes()
    local data::Array{Any, 1} = Any[]
    local t::DataType

    if dataIndex < UInt8(1) || dataIndex > UInt(length(types)) error("Protocol error: invalid dataIndex header") end
    t = types[dataIndex]
    while t !== Void
      push!(data, t === String ? readuntil(sock, '\x0')[1:end-1] : read(sock, t))
      t = types[(dataIndex += 1)]
    end

    data
  end
  #
  # Writes data in a socket in a fast way using native read()
  # and write() functions
  # TODO: arguments
  #
  function _fastWrite(sock::Base.TCPSocket, dataIndex::UInt8, data::Any)
    local vec::Vector{UInt8} = UInt8[dataIndex]

    for d::Any in data
      append!(vec, Helper.toBytes(d))
      #
      # Every string should be ended with '\x0' - zero character
      #
      if typeof(d) === String append!(vec, '\x0') end
    end
    write(sock, vec)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
    NoTransform()

Identity transform.

Specifies that no transformation should be applied.
"""
struct NoTransform <: AbstractTransform end

"""
    NoTransform!()

In-place version of [`NoTransform`](@ref).
"""
struct NoTransform! <: AbstractTransform end

const AnyNoTransform = Union{NoTransform, NoTransform!}

is_inplace(::NoTransform) = false
is_inplace(::NoTransform!) = true

binv(::T) where {T <: AnyNoTransform} = T()
length_output(::AnyNoTransform, length_in::Integer) = length_in
eltype_output(::AnyNoTransform, ::Type{T}) where T = T
eltype_input(::AnyNoTransform, ::Type) = Nothing
expand_dims(::T, ::Val{N}) where {T <: AnyNoTransform, N} =
    N == 0 ? () : (T(), expand_dims(T(), Val(N - 1))...)
scale_factor(::AnyNoTransform, A, dims) = 1

plan(::NoTransform, A, dims; kwargs...) = IdentityPlan()
plan(::NoTransform!, A, dims; kwargs...) = IdentityPlan!()

"""
    IdentityPlan

Type of plan associated to [`NoTransform`](@ref).
"""
struct IdentityPlan <: AbstractCustomPlan end

LinearAlgebra.mul!(y, ::IdentityPlan, x) = (y === x) ? y : copy!(y, x)
LinearAlgebra.ldiv!(y, ::IdentityPlan, x) = mul!(y, IdentityPlan(), x)
Base.:*(::IdentityPlan, x) = copy(x)
Base.:\(::IdentityPlan, x) = copy(x)

"""
    IdentityPlan!

Type of plan associated to [`NoTransform!`](@ref).
"""
struct IdentityPlan! end

function LinearAlgebra.mul!(y, ::IdentityPlan!, x)
    if x !== y
        throw(ArgumentError("in-place IdentityPlan applied to out-of-place data"))
    end
    y
end
LinearAlgebra.ldiv!(y, ::IdentityPlan!, x) = mul!(y, IdentityPlan!(), x)
Base.:*(::IdentityPlan!, x) = x
Base.:\(::IdentityPlan!, x) = x


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function calz_unred(wave::T, flux::T, ebv::T, r_v::T) where {T<:AbstractFloat}
    x  = 10000 / wave # Wavelength in inverse microns
    if 6300 <= wave <= 22000
        klam = 2.659*(-1.857 + 1.040*x) + r_v
    elseif 912 <= wave < 6300
        klam = 2.659*(@evalpoly(x, -2.156, 1.509, -0.198, 0.011)) + r_v
    else
        return flux
    end
    return flux * exp10(0.4 * klam * ebv)
end

"""
    calz_unred(wave, flux, ebv[, r_v]) -> deredden_wave

### Purpose ###

Deredden a galaxy spectrum using the Calzetti et al. (2000) recipe.

### Explanation ###

Calzetti et al.  (2000, ApJ 533, 682;
http://adsabs.harvard.edu/abs/2000ApJ...533..682C) developed a recipe for
dereddening the spectra of galaxies where massive stars dominate the radiation
output, valid between \$0.12\$ to \$2.2\$ microns.  (`calz_unred` extrapolates
between \$0.12\$ and \$0.0912\$ microns.)

### Arguments ###

* `wave`: wavelength (Angstroms)
* `flux`: calibrated flux.
* `ebv`: color excess E(B-V).  If a negative `ebv` is supplied, then
  fluxes will be reddened rather than deredenned.  Note that the supplied color
  excess should be that derived for the stellar continuum, EBV(stars), which is
  related to the reddening derived from the gas, EBV(gas), via the Balmer
  decrement by EBV(stars) = 0.44*EBV(gas).
* `r_v` (optional): ratio of total to selective extinction, default is
  4.05.  Calzetti et al. (2000) estimate ``r_v = 4.05 \\pm 0.80`` from optical-IR
  observations of 4 starbursts.

### Output ###

Unreddened flux, same units as `flux`.  Flux values will be left unchanged outside valid
domain (\$0.0912\$ - \$2.2\$ microns).

### Example ###

Estimate how a flat galaxy spectrum (in wavelength) between \$1200 Å\$ and
\$3200 Å\$ is altered by a reddening of E(B-V) = 0.1.

```julia
wave = collect(1200:50:3150);
flux = ones(wave);
flux_new = calz_unred.(wave, flux, -0.1);
```

Using a plotting tool you can visualize the unreddend flux.  For example, with
[PyPlot.jl](https://github.com/JuliaPy/PyPlot.jl)

```julia
using PyPlot
plot(wave, flux_new)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
calz_unred(wave::Real, flux::Real, ebv::Real, r_v::Real=4.05) =
    calz_unred(promote(float(wave), float(flux), float(ebv), float(r_v))...)


text:
using Images
using CUDA

function get_array_of_image(filename, resized = true)
    if resized
        image_path = "../../resources/drawables/" * filename
    else
        image_path = "../../resources/drawables/size_original/" * filename
    end
    img = load(image_path)
    return channelview(RGB.(img))
end

function includes_point(
    click_coordinate_x,
    click_coordinate_y,
    absolute_coordinate_x,
    absolute_coordinate_y,
    width,
    height,
)
    return (
        absolute_coordinate_x <= click_coordinate_x <= (absolute_coordinate_x + width) &&
        absolute_coordinate_y <= click_coordinate_y <= (absolute_coordinate_y + height)
    )
end

#written for cpu computation, probably needs to be changed for GPU
#function __blit_single_channel_inplace(dest,src,x,y)
#    dest[x:x+size(src,1)-1,y:y+size(src,2)-1] = src
#end

#number threads needs to be bigger than number of rows of src
#TODO make work for src images with more rows than threadnumber
function __blit_single_channel_inplace(threadID, dest, src, channel, x, y)
    if threadID <= size(src, 3)
        for i = 1:size(src, 2)
            @inbounds dest[channel, y-1+i, x-1+threadID] = src[channel,i, threadID]
        end
    end
    return
end
function __blit_single_channel_inplace(threadID,blockID, dest, src, channel, x, y)
    if ndims(src) == 3
        if threadID <= size(src, 3)
            for i = 1:size(src, 2)
                @inbounds dest[channel, y-1+i, x-1+threadID,blockID] = src[channel,i, threadID]
            end
        end
    else
        if threadID <= size(src, 3)
            for i = 1:size(src, 2)
                @inbounds dest[channel, y-1+i, x-1+threadID,blockID] = src[channel,i, threadID,blockID]
            end
        end
    end
    return
end

#function kernel_blit_image_inplace(threadID, dest, src, x, y)

#    __blit_single_channel_inplace(threadID, dest, src, 1, x, y)
#    __blit_single_channel_inplace(threadID, dest, src, 2, x, y)
#    __blit_single_channel_inplace(threadID, dest, src, 3, x, y)
#    return
#end
function kernel_blit_image_inplace(threadID,blockID, dest, src, x, y)

    __blit_single_channel_inplace(threadID, blockID, dest, src, 1, x, y)
    __blit_single_channel_inplace(threadID, blockID, dest, src, 2, x, y)
    __blit_single_channel_inplace(threadID, blockID, dest, src, 3, x, y)
    return
end

####
function get_blank_image_as_array(color_red,color_green,_color_blue,width,height)

    result_array = @cuStaticSharedMem(Float32,(height,width,3))
    function fill_array(array,color,channel)
        for i in 1:height
            for j in 1:width
                array[channel,i ,j] = color
            end
        end
    end
    fill_array(result_array,color_red,1)
    fill_array(result_array,color_green,2)
    fill_array(result_array,_color_blue,3)
end

function grayscale(dest,src)
    if threadIdx().x <= size(src,3)
        for i in 1:size(src,2)
            @inbounds dest[i,threadIdx().x,blockIdx().x] = src[1,i,threadIdx().x,blockIdx().x] * (0.21f0) + convert(Float32,src[2,i,threadIdx().x,blockIdx().x]) *(0.71f0) + convert(Float32,src[3,i,threadIdx().x,blockIdx().x]) * (0.07f0)
            if blockIdx().x == 1
            #@cushow(convert(Float32,src[1,i,threadIdx().x,blockIdx().x]))
            end
        end
    end
end

function fill_input(dest,src)
    if threadIdx().x <= size(src,3)
        for i in 1:size(src,1)
            for j in 1:size(src,2)
                @inbounds dest[i,j,threadIdx().x] = convert(Float32,src[i,j,threadIdx().x,blockIdx().x])
            end
        end
    end
end


text:
module NanStats

# package code goes here
include("nan_stats_def.jl")

export nanmean, nanvar, nanskew, nankurt

end # module


text:
"""
    NexDualStateButton(nexSerial, name; pid=pid, cid=cid)

A Nextion Button UI control.

On Nextion side, it have the following properties:
- bco0
- bco1
- pco
- font
- xcen
- ycen
- val
- txt
- isbr
"""
struct NexDualStateButton <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    booleanvalued::IBooleanValued
    colourable::IColourable
    touchable::ITouchable

    function NexDualStateButton(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), IBooleanValued(nid), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexDualStateButton, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexDualStateButton, property::Symbol)
    _getcommonproperty(obj, property)
end


text: struct income :> Capital end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using JSON

export deserialize

# Deserialize a JSON string to a given type.
#
# Example:
# ```
# using Base.Test
#
# type Foo
#    bar::UTF8String
# end
#
# @test DandelionSlack.deserialize(Foo, """{"bar": "baz"}""") == Foo("baz")
# ```
#
# Nullable fields are considered optional, but all other fields are mandatory.
#
DandelionSlack.deserialize{T}(::Type{T}, json::AbstractString) = DandelionSlack.deserialize(T, JSON.parse(json))
DandelionSlack.deserialize{T}(::Type{T}, json::Dict) =
    T([deserialize_field(T, f, json) for f in fieldnames(T)]...)
function DandelionSlack.deserialize{T}(::Type{Array{T,1}}, a::Array{Any,1})
    result = Array{T, 1}()
    for e in a
        element_value = DandelionSlack.deserialize(T, e)
        push!(result, element_value)
    end
    result
end
DandelionSlack.deserialize{T<:Integer}(::Type{T}, i::Integer) = T(i)
DandelionSlack.deserialize{T<:Integer}(::Type{T}, i::AbstractString) = parse(T, i)
DandelionSlack.deserialize{T<:AbstractFloat}(::Type{T}, i::AbstractFloat) = T(i)
DandelionSlack.deserialize{T<:AbstractFloat}(::Type{T}, i::AbstractString) = parse(T, i)
DandelionSlack.deserialize{T<:AbstractString}(::Type{T}, i::AbstractString) = T(i)
DandelionSlack.deserialize{T}(::Type{Nullable{T}}, json::AbstractString) = DandelionSlack.deserialize(T, json)
DandelionSlack.deserialize{T}(::Type{Nullable{T}}, json::Dict) = DandelionSlack.deserialize(T, json)
DandelionSlack.deserialize{T}(::Type{Nullable{T}}, ::Void) = Nullable{T}()
DandelionSlack.deserialize{T}(::Type{Nullable{T}}, x) = DandelionSlack.deserialize(T, x)

# Deserialize a given field (name or symbol) from a Dict, into a given type.
deserialize_field{Tf}(::Type{Tf}, field::AbstractString, json::Dict) = DandelionSlack.deserialize(Tf, json[field])
deserialize_field{Tf<:Nullable}(::Type{Tf}, field::AbstractString, json::Dict) =
    haskey(json, field) ? DandelionSlack.deserialize(Tf, json[field]) : Tf()
deserialize_field{Tf}(ta::Type{Array{Tf}}, field::AbstractString, json::Dict) =
    DandelionSlack.deserialize(ta, json[field])

deserialize_field{T}(::Type{T}, field::Symbol, json::Dict) =
    deserialize_field(fieldtype(T,field), string(field), json)

flatten{T<:Array}(A::Array{T,1}) = vcat(flatten(A)...)
flatten{T}(A::Array{T,1}) = vcat(A...)
flatten(A) = A