text:
using DataFrames
using CSV
using Calculus
using Plots
using Statistics

plotly()
@doc "
    traffic_assignment(;networkName, tol=0.001, maxIters=20)
"
function traffic_assignment(;networkName, tol=0.001, maxIters=20)
    networkFile = "Network\\$networkName\\network.csv"
    demandFile = "Network\\$networkName\\demand.csv"

    𝒩 = Int64[]                            # Nodes
    𝒜 = Array{Int64,1}[]                   # Arcs as adjacency list
    v₀ = Array{Float64,1}[]                 # Link volume capcity
    t₀ = Array{Float64,1}[]                 # Link free flow travel time
    α₀ = Array{Float64,1}[]                 # BPR parameters
    β₀ = Array{Float64,1}[]                 # BPR parameters

    W = Array{Int64,1}[]                    # OD pairs
    q = Float64[]                           # Demand between OD pairs
    P = Array{Array{Int64,1},1}[]           # Paths between OD pairs

    cⁱ = Float64[]                  # Stores total cost from each iteration
    xⁱₐ = Array{Array{Float64,1}}[] # Stores link flows from each iteration
    eⁱ = Float64[]                  # Stores convergent from each iteration

    # Enumerates all paths from pivot node i using nodes from a vector of open nodes 𝕏
    function enumerate(i, 𝕏)
        paths = Array{Int64,1}[]          # Vector of paths
        𝕏 = filter(x -> (x ≠ i), 𝕏)

        # Dead ends: if there is no outgoing arc from node i or if all outgoing nodes from i are closed
        if isempty(𝒜[i]) push!(paths, [i]) end
        if sum([if 𝒜[i][j] in 𝕏 1 else 0 end for j in 1:length(𝒜[i])]) == 0 push!(paths, [i]) end

        # Dveveloping paths recursively with subpaths
        for j in 1:length(𝒜[i])
            if 𝒜[i][j] in 𝕏
                for subpath in enumerate(𝒜[i][j], 𝕏)
                    push!(paths, [i])
                    append!(paths[end], subpath)
                end
            end
        end
        return paths
    end

    # Network build
    function build()
        println()
        # File import
        network = CSV.File(networkFile, types=[Int64, Int64, Float64, Float64, Float64, Float64, Float64])
        df = DataFrame(network)
        head = df[1]::Array{Int64,1}
        tail = df[2]::Array{Int64,1}
        linkcapacity = df[3]::Array{Float64,1}
        linklength = df[4]::Array{Float64,1}
        linkfft = df[5]::Array{Float64,1}
        alpha = df[6]::Array{Float64,1}
        beta = df[7]::Array{Float64,1}

        n = max(maximum(head), maximum(tail))
        for i in 1:n
            append!(𝒩, i)
            push!(𝒜, [])
            push!(v₀, [])
            push!(t₀, [])
            push!(α₀, [])
            push!(β₀, [])
        end

        for i in 1:length(head)
            append!(𝒜[head[i]], tail[i])
            append!(v₀[head[i]], linkcapacity[i])
            append!(t₀[head[i]], linkfft[i])
            append!(α₀[head[i]], alpha[i])
            append!(β₀[head[i]], beta[i])
        end

        demand = CSV.File(demandFile, types=[Int64, Int64, Float64])
        df′ = DataFrame(demand)
        origin = df′[1]::Array{Int64,1}
        destination = df′[2]::Array{Int64,1}
        demand = df′[3]::Array{Float64,1}

        for i in 1:nrow(df′)
            push!(W, [origin[i], destination[i]])
            append!(q, demand[i])
            push!(P, [])
        end

        for r in unique(origin)
            S = [W[k][2] for k in findall(x -> (x == r), origin)]
            for s in S
                index = findfirst(x -> (x == [r,s]), W)
                paths = Array{Int64,1}[]
                for path in enumerate(r, copy(𝒩))
                    if s in path
                        i = findfirst(x -> (x == s), path)
                        push!(paths, path[1:i])
                    end
                end
                append!(P[index], paths)
            end
        end
    end

    # Arc cost
    function 𝐶ₐ(i, j, x)
        α = α₀[i][j]
        β = β₀[i][j]
        t = t₀[i][j]
        v = v₀[i][j]
        return t * (1 + α * (x/v) ^ β)
    end

    # Path cost for path p given arc flows x and arc cost function cₐ
    function 𝐶ₚ(p, x, cₐ)
        cost = 0.0
        for k in 2:length(p)
            i = p[k-1]
            j = findfirst(x -> (x == p[k]), 𝒜[i])::Int64
            cost += cₐ(i, j, x[i][j])
        end
        return cost
    end

    # Djikstra's shortest path between r and s
    function sp(𝒞, r, s)
        p = Int64[]
        ℒ = [[if i==r r else -1 end, if i==r 0.0 else Inf end] for i in 𝒩]    # predecessor and cost
        𝕏 = copy(𝒩)   # Set of open nodes
        i = r
        deleteat!(𝕏, i)
        while !isempty(𝕏)
            for j in 1:length(𝒜[i])
                k = 𝒜[i][j]
                C = ℒ[i][2] + 𝒞[i][j]
                if C < ℒ[k][2] && k in 𝕏
                    ℒ[k][1] = i
                    ℒ[k][2] = C
                end
            end
            index = argmin([ℒ[i][2] for i in 𝕏])
            i = 𝕏[index]
            deleteat!(𝕏, index)
        end

        i = s
        append!(p, i)
        while i ≠ r
            i = Int(ℒ[i][1])
            append!(p, i)
        end
        reverse!(p)
        return p
    end

    # General Equilibration Algorithm
    function gea(;cₐ, ε, kₘₐₓ)
        cₚ = 𝐶ₚ
        ε′ = 0.000001
        k′ₘₐₓ = 20

        xᵏₐ = Array{Float64,1}[[0.0 for j in 𝒜[i]] for i in 𝒩]
        fᵏₚ = Array{Float64,1}[[0.0 for _ in P[n]] for n in 1:length(W)]

        # Intialization
        cᵒₐ = [[cₐ(i, j, 0.0) for j in 1:length(𝒜[i])] for i in 𝒩]
        for n in 1:length(W)
            r, s = W[n]
            qᵣₛ = q[n]
            pᵣₛ = sp(cᵒₐ, r, s)
            m = findfirst(x -> (x == pᵣₛ), P[n])
            fᵏₚ[n][m] = qᵣₛ
            for l in 2:length(pᵣₛ)
                i = pᵣₛ[l-1]
                j = findfirst(x -> (x == pᵣₛ[l]), 𝒜[i])
                xᵏₐ[i][j] += fᵏₚ[n][m]
            end
        end

        # Iteration
        allConverged = false
        k = 1
        while !allConverged
            xᵏ⁺¹ₐ = deepcopy(xᵏₐ)
            fᵏ⁺¹ₚ = deepcopy(fᵏₚ)

            for n in 1:length(W)
                odConverged=false
                k′ = 1
                while !odConverged
                    c = [cₚ(P[n][m], xᵏ⁺¹ₐ, cₐ) for m in 1:length(P[n])]
                    indices = findall(x -> (x > 0), fᵏ⁺¹ₚ[n])

                    # Active path with highest costs
                    m₁ = indices[argmax([c[m] for m in indices])]
                    # Path with least cost
                    m₂ = argmin(c)

                    X = [[P[n][m₁][l-1], P[n][m₁][l]] for l in 2:length(P[n][m₁])]
                    Y = [[P[n][m₂][l-1], P[n][m₂][l]] for l in 2:length(P[n][m₂])]
                    Z = intersect(X, Y)
                    num = c[m₁] - c[m₂]
                    den = Float64(length(X) + length(Y) - 2*length(Z))
                    if num == den == 0 δ = 0.0
                    else δ = min(num/den, fᵏ⁺¹ₚ[n][m₁])
                    end

                    # OD convergence test
                    if δ < ε′ || k′ == k′ₘₐₓ odConverged = true end

                    # Updation
                    fᵏ⁺¹ₚ[n][m₁] -= δ
                    fᵏ⁺¹ₚ[n][m₂] += δ
                    for l in 2:length(P[n][m₁])
                        i = P[n][m₁][l-1]
                        j = findfirst(x -> (x == P[n][m₁][l]), 𝒜[i])
                        xᵏ⁺¹ₐ[i][j] -= δ
                    end
                    for l in 2:length(P[n][m₂])
                        i = P[n][m₂][l-1]
                        j = findfirst(x -> (x == P[n][m₂][l]), 𝒜[i])
                        xᵏ⁺¹ₐ[i][j] += δ
                    end

                    k′ += 1
                end
            end

            # Overall convergence test
            num = 0.0
            den = 0.0
            for i in 𝒩
                for j in 1:length(𝒜[i])
                    num += (xᵏ⁺¹ₐ[i][j] - xᵏₐ[i][j])^2
                    den += xᵏₐ[i][j]
                end
            end
            Δ = sqrt(num)/den
            if Δ < ε || k == kₘₐₓ allConverged = true end

            # Updation
            xᵏₐ = deepcopy(xᵏ⁺¹ₐ)
            fᵏₚ = deepcopy(fᵏ⁺¹ₚ)

            k += 1
        end
        return xᵏₐ
    end

    # Variational Inequality: Projection method
    function vi(;cₐ=𝐶ₐ, ε=tol, kₘₐₓ=maxIters)
        xᵏₐ = Array{Float64,1}[[0.0 for j in 𝒜[i]] for i in 𝒩]

        # Intialization
        cₒ = [[cₐ(i, j, 0.0) for j in 1:length(𝒜[i])] for i in 𝒩]
        for n in 1:length(W)
            r, s = W[n]
            qᵣₛ = q[n]
            pᵣₛ = sp(cₒ, r, s)
            for l in 2:length(pᵣₛ)
                i = pᵣₛ[l-1]
                j = findfirst(x -> (x == pᵣₛ[l]), 𝒜[i])
                xᵏₐ[i][j] += qᵣₛ
            end
        end

        f′(i,j,x) = derivative(x -> cₐ(i,j,x))(x)
        L = sum([(f′(i, j, xᵏₐ[i][j]))^2 for i in 𝒩 for j in 1:length(𝒜[i])])^0.5
        γ = 1/L
        println("gamma: $γ")

        println("\niter: 0")
        println("total flow: ", sum(sum.(xᵏₐ)))
        println("total cost: ", sum([xᵏₐ[i][j] * cₐ(i, j, xᵏₐ[i][j]) for i in 𝒩 for j in 1:length(𝒜[i])]))
        println("convergent: NA")

        append!(cⁱ, sum([xᵏₐ[i][j] * cₐ(i, j, xᵏₐ[i][j]) for i in 𝒩 for j in 1:length(𝒜[i])]))
        push!(xⁱₐ, xᵏₐ)

        # Iteration
        converged = false
        k = 1
        while !converged
            xᵏ⁺¹ₐ = Array{Float64,1}[[0.0 for j in 𝒜[i]] for i in 𝒩]

            # Auxilary problem
            c(i,j,x) = x - (xᵏₐ[i][j] - γ*cₐ(i,j,xᵏₐ[i][j]))
            xᵏ⁺¹ₐ = gea(cₐ=c, ε=0.000001, kₘₐₓ=20)

            # Convergence test
            Δ = sqrt(sum(vcat((xᵏ⁺¹ₐ - xᵏₐ)...) .^2))/sum(vcat(xᵏₐ...))
            if Δ < ε || k == kₘₐₓ converged = true end

            # Updating link flows and cost
            xᵏₐ = deepcopy(xᵏ⁺¹ₐ)

            println("\niter: $k")
            println("total flow: ", sum(sum.(xᵏₐ)))
            println("total cost: ", sum([xᵏₐ[i][j] * cₐ(i, j, xᵏₐ[i][j]) for i in 𝒩 for j in 1:length(𝒜[i])]))
            println("convergent: $Δ")

            append!(cⁱ, sum([xᵏₐ[i][j] * cₐ(i, j, xᵏₐ[i][j]) for i in 𝒩 for j in 1:length(𝒜[i])]))
            push!(xⁱₐ, xᵏₐ)
            append!(eⁱ, Δ)

            k += 1
        end
        return xᵏₐ
    end

    build()
    println(vi())
    return cⁱ, xⁱₐ, eⁱ
end
traffic_assignment(networkName="Braess-Example", tol=0.001, maxIters=20)


text:
function rho_divergence(a,b)
    b .- a .+ a .* log.( a./b )
end

"""
wasserstein_distance_greenkhorn(p, q, C, ε, L)
p: right marginal N vector
q: left marginal N vector
C: Cost NxN matrix
ε: entropic regularization parameter
L: number of Sinkhorn iterations

returns a tupel of entropic Wasserstein distance W and scaling potentials u,v as N vectors

TODO: number of grid points different in marginals
"""
function wasserstein_distance_greenkhorn(p, q, C, ε, L)

    @assert length(p) == size(C,1)
    @assert length(q) == size(C,2)

    a = ones(length(p)) # a = exp(u)
    b = ones(length(q)) # b = exp(v)

    K = get_gibbs_matrix(C, ε)

    local r_Gamma = K * b # initial row sum
    local c_Gamma = K * a # initial column sum

    for l in 1:L
        local m_i, I = findmax( rho_divergence(p, r_Gamma) )
        local m_j, J = findmax( rho_divergence(q, c_Gamma) )
        if m_i > m_j
            a_old = a[I]
            a[I] *= p[I] / r_Gamma[I] # update Ith entry of left scaling vector
            r_Gamma[I] *= a[I] / a_old # update row sum
            # c_Gamma .+= (a[I] - a_old) .* K[I,:] .* b # update column sum - O(N)
            c_Gamma .= b .* K * a # O(N^2)
        else
            b_old = b[J]
            b[J] *= q[J] / c_Gamma[J] # update Jth entry of right scaling vector
            c_Gamma[J] *= b[J] / b_old # update column sum
            # r_Gamma .+= a .* K[:,J] .* (b[J] - b_old) # update row sum
            r_Gamma .= a .* K * b
        end
    end

    return get_distance_from_scaling(a, b, C, ε), a, b
end

"""
wasserstein_distance_greenkhorn_logstable(p, q, C, ε, L)
p: right marginal N vector
q: left marginal N vector
C: Cost NxN matrix
ε: entropic regularization parameter
L: number of Sinkhorn iterations

returns a tupel of entropic Wasserstein distance W and scaling potentials u,v as N vectors

TODO: number of grid points different in marginals
"""
function wasserstein_distance_greenkhorn_logstable(p, q, C, ε, L)

    @assert length(p) == size(C,1)
    @assert length(q) == size(C,2)

    u = zeros(length(p)) # a = exp(u)
    v = zeros(length(q)) # b = exp(v)
    log_p = log.(p)
    log_q = log.(q)

    K = get_gibbs_matrix(C, ε)

    local r_Gamma = K * exp.(v) # initial row sum
    local c_Gamma = K * exp.(u) # initial column sum

    for l in 1:L
        local m_i, I = findmax( rho_divergence(p, r_Gamma) )
        local m_j, J = findmax( rho_divergence(q, c_Gamma) )
        if m_i > m_j
            u_old = u[I]
            u[I] += log_p[I] - log(r_Gamma[I]) # update Ith entry of left scaling potential
            r_Gamma[I] *= exp(u[I] - u_old) # update row sum
            @. c_Gamma += exp(u[I] - C[I,:]/ε + v) - exp(u_old - C[I,:]/ε + v) # update column sum
        else
            v_old = v[J]
            v[J] += log_q[J] - log(c_Gamma[J]) # update Jth entry of right scaling vector
            c_Gamma[J] *= exp(v[J] - v_old) # update column sum
            @. r_Gamma += exp(u - C[:,J]/ε + v[J]) - exp(u - C[:,J]/ε + v_old) # update row sum
        end
    end

    return get_distance_from_scaling_logstable(u, v, C, ε), u, v
end

ext: jl

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
#=
In the card game poker, a hand consists of five cards and are ranked, from lowest to highest, in the following way:

High Card: Highest value card.
One Pair: Two cards of the same value.
Two Pairs: Two different pairs.
Three of a Kind: Three cards of the same value.
Straight: All cards are consecutive values.
Flush: All cards of the same suit.
Full House: Three of a kind and a pair.
Four of a Kind: Four cards of the same value.
Straight Flush: All cards are consecutive values of same suit.
Royal Flush: Ten, Jack, Queen, King, Ace, in same suit.
The cards are valued in the order:
2, 3, 4, 5, 6, 7, 8, 9, 10, Jack, Queen, King, Ace.

If two players have the same ranked hands then the rank made up of the highest value wins; for example, a pair of eights beats a pair of fives (see example 1 below). But if two ranks tie, for example, both players have a pair of queens, then highest cards in each hand are compared (see example 4 below); if the highest cards tie then the next highest cards are compared, and so on.

Consider the following five hands dealt to two players:

Hand	 	Player 1	 	Player 2	 	Winner
1	 	5H 5C 6S 7S KD
Pair of Fives
 	2C 3S 8S 8D TD
Pair of Eights
 	Player 2
2	 	5D 8C 9S JS AC
Highest card Ace
 	2C 5C 7D 8S QH
Highest card Queen
 	Player 1
3	 	2D 9C AS AH AC
Three Aces
 	3D 6D 7D TD QD
Flush with Diamonds
 	Player 2
4	 	4D 6S 9H QH QC
Pair of Queens
Highest card Nine
 	3D 6D 7H QD QS
Pair of Queens
Highest card Seven
 	Player 1
5	 	2H 2D 4C 4D 4S
Full House
With Three Fours
 	3C 3D 3S 9S 9D
Full House
with Three Threes
 	Player 1
The file, poker.txt, contains one-thousand random hands dealt to two players. Each line of the file contains ten cards (separated by a single space): the first five are Player 1's cards and the last five are Player 2's cards. You can assume that all hands are valid (no invalid characters or repeated cards), each player's hand is in no specific order, and in each hand there is a clear winner.

How many hands does Player 1 win?
=#

using DelimitedFiles
import Base.<

struct PokerHand
  cards::Dict{Int, Int}
  suits::Dict{Int, Int}
  function PokerHand(cardArray, suitArray)
    cards, suits = Dict{Int, Int}(), Dict{Int, Int}()
    for i in zip(cardArray, suitArray)
      cardvalue, suitvalue = get(cards, i[1], 0), get(suits, i[2], 0)
      cards[i[1]], suits[i[2]] = cardvalue + 1, suitvalue + 1
    end
    new(cards, suits)
  end
end

function handRank(hand::PokerHand)
  if isFlush(hand) && collect(values(hand.cards)) |> length == 5 && sort(collect(keys(hand.cards))) == sort(collect(10:14))
    return 10 # Royal Flush
  end
  if isFlush(hand) && isStraight(hand)
    return 9 # Straight Flush
  end
  if filter(x->x==4, collect(values(hand.cards))) |> length == 1
    return 8 # Four of a Kind
  end
  if filter(x->x==3, collect(values(hand.cards))) |> length == 1 && filter(x->x==2, collect(values(hand.cards))) |> length == 1
    return 7 # Full House
  end
  if isFlush(hand)
    return 6 # Flush
  end
  if isStraight(hand)
    return 5 # Straight
  end
  if filter(x->x==3, collect(values(hand.cards))) |> length == 1
    return 4 # Three of A Kind
  end
  if filter(x->x==2, collect(values(hand.cards))) |> length == 2
    return 3 # Two Pair
  end
  if filter(x->x==2, collect(values(hand.cards))) |> length == 1
    return 2 # One Pair
  end
  1 # High Card
end

function isFlush(hand::PokerHand)
  length(keys(hand.suits)) == 1
end

function isStraight(hand::PokerHand)
  k = reverse(sort(collect(keys(hand.cards))))
  if k |> length != 5 return false end
  return k[1] - k[5] == 4
end

function histogram(hand::PokerHand)
  sortedWithinRanks = Dict{Int, Array{Int, 1}}()
  for (key, value) in hand.cards
    existing = get(sortedWithinRanks, value, Int[])
    if existing |> length == 0
      sortedWithinRanks[value] = [key]
    else
      sortedWithinRanks[value] = push!(existing, key)
    end
  end
  for (key, value) in sortedWithinRanks
    sortedWithinRanks[key] = reverse(sort(value))
  end
  result, hist = Int[], Int[]
  for key in reverse(sort(collect(keys(sortedWithinRanks))))
    for rank in sortedWithinRanks[key] |> sort |> reverse
      push!(hist, key)
      append!(result, fill(rank, key))
    end
  end
  (hist, result)
end

function <(x::PokerHand, y::PokerHand)
  rank1, rank2, rank = handRank(x), handRank(y), handRank(x)
  if rank1 < rank2
    return true
  elseif rank1 != rank2
    return false
  end
  xRanks, yRanks = histogram(x)[2], histogram(y)[2]
  for i in zip(xRanks, yRanks)
    if i[1] < i[2] return true end
    if i[1] > i[2] return false end
  end
  false
end

function cardvalue(n)
  values = ["2", "3", "4", "5", "6", "7", "8", "9", "T", "J", "Q", "K", "A"]
  findfirst(x -> x == n, values)
end

function suitvalue(n)
  values = ["H", "C", "D", "S"]
  findfirst(x -> x == n, values)
end

function calc()
  count = 0
  hands = readdlm(dirname(@__FILE__()) * "/../Resources/p54.txt", ' ', AbstractString)
  for i in 1:1000
      suits, cards = Int[], Int[]
      for j in 1:10
        card = hands[i,j]
        cardstring = chop(card)
        suitstring = card[length(card):length(card)]
        push!(cards, cardvalue(cardstring))
        push!(suits, suitvalue(suitstring))
      end
      hand1 = PokerHand(cards[1:5], suits[1:5])
      hand2 = PokerHand(cards[6:10], suits[6:10])
      if (hand1 < hand2) == false
        count = count + 1
      end
  end
  count
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
export j, pol, ∠, phasor, phasorsine, angulardimension, phasordimension

"""
`j = 1im` equals the imaginary unit
"""
const j=1im

"""
# Function call

`pol(r, phi)`

# Description

Creates a complex quantity with length `r` and angle `phi`

# Variables

`r` Length of complex quantity; r may be utilized including a unit generated by
Unitful

`phi` Angle of complex quantity; if module Unitful is utilized, the angle may be
specified in degrees, by using unit `°`

# Examples

```julia
julia> using Unitful, Unitful.DefaultSymbols, ElectricalEngineering
julia> U1 = pol(2V,pi)
-2 + 0im V
julia> U2 = pol(sqrt(2)*1V,45°)
1 + 1im V
```
"""
function pol(r, phi)
  return r*cos(phi) + 1im*r*sin(phi)
end

"""
# Function call

`∠(phi)`

# Description

Creates a complex quantity with length 1 and angle `phi`

# Variables

`phi` Angle of complex quantity; if module Unitful is utilized, the angle may be
specified in degrees, by using unit `°`

# Examples

```julia
julia> using Unitful, Unitful.DefaultSymbols, ElectricalEngineering
julia> U1 = 2V*∠(90°)
-0.0 + 2.0im V
```
"""
function ∠(phi)
  return cos(phi) + 1im*sin(phi)
end

@doc raw"""
# Function call

```
phasor(c;origin=0.0+0.0im, ref=1, par=0,
    labelrsep=0.5, labeltsep=0.1, label="", ha="center", va="center",
    labelrelrot=false, labelrelangle=0,
    color="black", backgroundcolor="none", linesstyle="-", linewidth=1,
    width=0.2, headlength=10, headwidth=5)
```

# Description

This function draws a phasor from a starting point `origin` and end point
`origin`+`c`. The phasor consists of a shaft and an arrow head.

Each phasor `c` is plotted as a relative quantity, i.e., `c/ref` is actually
shown the plot figure. This concept of plotting a per unit phasor is used to
be able to plot phasor with different quantities, e.g., voltage and current
phasors. It is important that the variables `c`, `origin` and `ref` have the
same units (defined through Unitful).

# Variables

`c` Complex phasor, drawn relative relative to `origin`

`origin` Complex number representing the origin of the phasor; this variable
needs to have the same unit as `c`

`ref` Reference length of scaling; this is required as in a phasor diagram
voltages and currents may be included; in order to account for the different
voltage and current scales, one (constant)  `ref` is used for voltage
phasors and another (constant) `ref` is used for current phasors; this
variable needs to have the same unit as `c`

`par` In order to be able to plot parallel phasors, par is used to specify
the per unit tangential shift (offset) of a phasor, with respect to `ref`;
so typically `ref` will be selected to be around 0.05 to 0.1;
default value = 0 (no shift of phasor)

`labelrsep` Radial per unit location of label (in direction of the phasor):
`labelrsep = 0` represents the shaft of the phasor and `labelrsep = 1` represents
the arrow hear of the phasor; default value = 0.5, i.e., the radial center
of the phasor

`labeltsep` Tangential per unit location of label: `labeltsep = 0` means that the
label is plotted onto the phasor; `labeltsep = 0.1` plots the label on top of
the phasor applying a displacement of 10% with respect to `ref`;
`labeltsep = -0.2` Plots the label below the
phasor applying a displacement of 20% with respect to `ref`; default value
= 0.1

`ha` Horizontal alignment of label; actually represents the tangential
alignment of the label; default value = "center"

`va` Vertical alignment of label; actually represents the radial
alignment of the label; default value = "center"

`labelrelrot` Relative rotation of label; if `labelrelrot == false` (default
value) then the label is not rotated relative to the orientation of the phasor;
otherwise the label is rotated relative to the phasor by the angle
`labelrelangle`

`labelrelangle` Relative angle of label  with respect to phasor orientation;
this angle is only applied, it `labelrelrot == true`; this angle the indicates
the relative orientation of the label with respect to the orientation of the
phasor; default value = 0

`color` Color of the phasor; i.e., shaft and arrow head color; default
value = "black"

`backgroundcolor` Background color of the label; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`linestyle` Line style of the phasor; default value = "-"

`linewidth` Line width of the phasor; default value = 1

`width` Line width of the shaft line; default value = 0.2

`headlength` Length of arrow head; default value = 10

`headwidth` Width of arrow head; default value = 5

# Example

Copy and paste code:

```julia
using Unitful, Unitful.DefaultSymbols, PyPlot, ElectricalEngineering
figure(figsize=(3.3, 2.5))
rc("text", usetex=true); rc("font", family="serif")
V1 = 100V + j*0V # Voltage
Z1 = 30Ω + j*40Ω # Impedance
I1 = V1/Z1 # Current
Vr = real(Z1)*I1
Vx = V1 - Vr
refV = abs(V1); refI=abs(I1)*0.8
phasor(V1, label=L"$\underline{V}_1$", labeltsep=0.1, ref=refV,
    labelrelrot=true)
phasor(Vr, label=L"$\underline{V}_r$", labeltsep=-0.25, ref=refV,
    labelrelrot=true)
phasor(Vx, origin=Vr, label=L"$\underline{V}_x$", labeltsep=-0.2, ref=refV,
    labelrelrot=true)
phasor(I1, label=L"$\underline{I}_1$", labeltsep=-0.2, labelrsep=0.7, ref=refI,
    labelrelrot=true, linestyle="--", par=-0.05)
phi1=angle(I1)
phi2=angle(V1)
angulardimension(0.3,phi1,phi2,arrowstyle1=".",arrowstyle2="-|>",ha="left",
    label=L"$\varphi_1$", labelrsep=0.05)
axis("square"); xlim(-1,1); ylim(-1,1)
removeaxes(); # Remove axis
# save3fig("phasordiagram",crop=true);
```
"""
function phasor(c;
    origin = (0.0+0.0im).*c./ustrip(c),
    ref = abs(c./ustrip(c)),
    par = 0.0,
    labelrsep = 0.5,
    labeltsep = 0.1,
    label = "",
    ha = "center",
    va = "center",
    labelrelrot = false,
    labelrelangle = 0.0,
    color = "black",
    backgroundcolor = "none",
    linestyle = "-",
    linewidth = 1,
    width = 0.2,
    headlength = 10.0,
    headwidth = 5.0)

    # Check if units of c, origin and ref are compatible
    # Starting point (origin) of phase
    xorigin = 0.0
    yorigin = 0.0
    xend = 0.0
    yend = 0.0
    try
        xorigin = uconvert(Unitful.NoUnits, real(origin)./ref)
        yorigin = uconvert(Unitful.NoUnits, imag(origin)./ref)
        # End point of phasor
        xend = uconvert(Unitful.NoUnits, real(origin+c)./ref)
        yend = uconvert(Unitful.NoUnits, imag(origin+c)./ref)
    catch err
        error("module ElectricalEngineering: function phasor: Dimension mismatch of arguments `c`, `origin` and `ref`\n    The arguments `c`, `origin` and `ref` must have the same dimension (koherent SI unit)")
    end

    # Draw phasor only if length of c is greater than zero, see
    # https://github.com/christiankral/ElectricalEngineering.jl/issues/1
    if upstrip(abs(c))>0
        # Length of phasor
        dr = sqrt((xend-xorigin)^2 + (yend-yorigin)^2)
        # Real part of phasor
        drx = (xend - xorigin)/dr # = real(c)./ref
        # Imag part of phasor
        dry = (yend - yorigin)/dr # = imag(c)./ref
        # Angle of phasor
        absangle = atan(dry, drx)
        # Orientation tangential to phasor (lagging by 90°)
        # Real part of tangential component with respect to length
        dtx = +dry
        # Imag part of tangential component with respect to length
        dty = -drx
        # Real part of parallel shift of phasor
        dpx = -par*dtx
        # Imag part of parallel shift of phasor
        dpy = -par*dty
        # Origin of head
        xoriginHead = xorigin + dr*drx*0.999
        yoriginHead = yorigin + dr*dry*0.999
        # Draw arrow shaft and head
        # https://matplotlib.org/api/_as_gen/matplotlib.pyplot.annotate.html?highlight=annotate#matplotlib.pyplot.annotate
        # Draw shaft separately: otherwise, the arrow contour will be drawn as in
        # https://matplotlib.org/users/annotations.html#basic-annotation
        # so that the back and forth paths overlap and the line style does not
        # appear correctly; replace
        plot([xorigin+dpx,xend+dpx], [yorigin+dpy,yend+dpy],
            color=color, linestyle=linestyle, linewidth=linewidth, clip_on=false)
        # Code based on plot replaces the previous implementation inspired by:
        # https://stackoverflow.com/questions/51746400/linestyle-in-plot-and-annotate-are-not-equal-in-matplotlib
        # Previous (obsolete) implementation:
        # annotate("", xy=(xend+dpx,yend+dpy),
        #     xytext=(xorigin+dpx,yorigin+dpy), xycoords="data",
        #     arrowprops=Dict("arrowstyle"=>"-",
        #         "linestyle"=>linestyle,
        #         "linewidth"=>linewidth,
        #         "color"=>color,
        #         "facecolor"=>color),
        #     annotation_clip=false)

        # Draw arrow head without line style; this is a workaround explained in
        # https://stackoverflow.com/questions/47180328/pyplot-dotted-line-with-fancyarrowpatch/47205418#47205418
        annotate("", xy=(xend+dpx, yend+dpy),
            xytext=(xoriginHead+dpx, yoriginHead+dpy), xycoords="data",
            arrowprops=Dict("edgecolor"=>color, "facecolor"=>color,
                "width"=>width, "linestyle"=>"-",
                "headlength"=>headlength,
                "headwidth"=>headwidth),
            annotation_clip=false)

        # Plot label
        if labelrelrot == false
            # Without relative rotation of label
            text(xorigin + dr*drx*labelrsep - dtx*labeltsep + dpx,
                yorigin + dr*dry*labelrsep - dty*labeltsep + dpy,
                label, ha=ha, va=va, rotation=labelrelangle*180/pi,
                backgroundcolor=backgroundcolor)
        else
            # Applying relative rotation of label
            text(xorigin + dr*drx*labelrsep - dtx*labeltsep + dpx,
                yorigin + dr*dry*labelrsep - dty*labeltsep + dpy,
                label, ha=ha, va=va, rotation=(absangle+labelrelangle)*180/pi,
                backgroundcolor=backgroundcolor)
        end
    end
end

@doc raw"""
# Function call

```
phasorsine(mag = 1, phi = 0; add = false, figsize = (6.6,2.5),
    xlabel = L"$\omega t$\,($^\circ $)", ylabel = "", maglabel = "",
    phasorlabel = maglabel, labeltsep = 0.1, labelrsep = 0.5,
    labelrelrot = true, labelrelangle = 0,
    color = "black", linewidth = 1, linestyle = "-",
    colorDash="gray", left=0.20, right=0.80, bottom=0.20, top=0.80,
    showsine=true, showdashline=true)
```

# Description

This function draws a phasor with magnitude between 0 and 1 on the left subplot
of a figure and a sine diagram on the right subplot of the figure. Such graph
is used in electrical engineering to explain the relationship between phasors
and time domain wave forms.

# Variables

`mag` Magnitude of displayed phasor; shall be between 0 and 1;
default value = 1

`phi` Phase angle of the phasor; default value = 0

`add` When calling this function the first time, `add` shall be equal to
`false`, which is the default value. In this case a new figure with the
dimensions specified in `figsize` is created. In order to add a second phasor
and sine diagram to an existing figure, `add` has to be set `true`

`figsize` Size of the new figure, if `add` is equal to `false`; default value
= (6.6,2.5)

`xlabel` Label of x-axis of right subplot; default value ="ωt(°)" in LaTeX
notation

`ylabel` label of y-axis of right subplot; default value = "";  if more than one
phasors and sine diagram shall be drawn (`add = true`), this label is displayed
only once; therefore, one has to create the label of all phasors when function
`phasorsine` is called the first time for creating a figure

`maglabel` Label of positive and negative magnitude `mag` on the right subplot;
dafault value = "";

`phasorlabel` Label of phasor of left subplot; default value = `maglabel`

`labelrsep` Radial per unit location of label (in direction of the phasor):
`labelrsep = 0` represents the shaft of the phasor and `labelrsep = 1`
represents the arrow hear of the phasor; default value = 0.5, i.e., the radial
center of the phasor

`labeltsep` Tangential per unit location of label: `labeltsep = 0` means that
the label is plotted onto the phasor; `labeltsep = 0.1` plots the label on top
of the phasor applying a displacement of 10% with respect to `ref`; `labeltsep =
-0.2` Plots the label below the phasor applying a displacement of 20% with
respect to `ref`; default value = 0.1

`labelrelrot` Relative rotation of label; if `labelrelrot == false` (default value)
then the label is not rotated relative to the orientation of the phasor;
otherwise the label is rotated relative to the phasor by the angle
`labelrelangle`

`labelrelangle` Relative angle of label with respect to phasor
orientation; this angle is only applied, it `labelrelrot == true`; this angle the
indicates the relative orientation of the label with respect to the
orientation of the phasor; default value = 0

`color` Color of the phasor; i.e., shaft and arrow head color; default
value = "black"

`backgroundcolor` Background color of all labels; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`linewidth` Line width of the phasor; default value = 1

`linestyle` Line style of the phasor; default value = "-"

`colorDash` Color of the dashed circle (left subplot) and the horizontal dashed
lines between the left and right subplot; default value = colorBlack4

`left` Left side of the figure; default value = 0.2

`right` Right side of the figure; default value = 0.85

`bottom` Bottom side of the figure; default value = 0.2

`top` Top side of the figure; default value = 0.85

`showsine` If `true`, the sinewave is shown; default value = `true`

`showdashline` If `true`, the dashed lines are shown; default value = `true`

# Example

Copy and paste code:

```julia
using Unitful, Unitful.DefaultSymbols, PyPlot, ElectricalEngineering
rc("text", usetex=true); rc("font", family="serif")
phasorsine(1, 45°, ylabel=L"$u,i$", maglabel=L"$\hat{U}$", labelrsep=0.3,
    color="gray", linestyle="--")
phasorsine(0.55, 0, add=true, maglabel=L"$\hat{I}$")
# save3fig("phasorsine",crop=true);
```
"""
function phasorsine(mag = 1,
    phi = 0;
    add = false,
    figsize = (6.6,2.5),
    xlabel = L"$\omega t$\,($^\circ $)",
    ylabel = "",
    maglabel = "",
    phasorlabel = maglabel,
    color = "black",
    backgroundcolor = "none",
    linewidth = 1,
    linestyle = "-",
    labeltsep = 0.1,
    labelrsep = 0.5,
    labelrelrot = true,
    labelrelangle = 0,
    colorDash="gray",
    left=0.20,
    right=0.80,
    bottom=0.20,
    top=0.80,
    showsine = true,
    showdashline = true)
    # https://matplotlib.org/tutorials/text/annotations.html#plotting-guide-annotation
    # https://matplotlib.org/users/annotations.html
    # https://stackoverflow.com/questions/17543359/drawing-lines-between-two-plots-in-matplotlib

    # Create figure
    if !add
        # Create new figure
        fig = figure(figsize=figsize)
        subplots_adjust(left=left, right=right, bottom=bottom, top=top)
    end
    # Create left subplot
    subplot(121)
    # Angle vector to draw circle
    psi = collect(0:pi/500:2*pi)
    # Coordinates of circle
    x = mag*cos.(psi)
    y = mag*sin.(psi)
    # Plot circle
    plot(x, y, color=colorDash, linewidth=1, linestyle=":",
        dash_capstyle="round")
    # Plot phasor
    phasor(pol(mag,phi), ref=1,
        label=phasorlabel, labelrsep=labelrsep, labeltsep=labeltsep,
        labelrelrot=labelrelrot, labelrelangle=labelrelangle,
        color=color, backgroundcolor=backgroundcolor,
        linestyle=linestyle, linewidth=linewidth)
    axis("square")
    ax1 = gca()
    if !add
        xlim(-1.1,1.1)
        ylim(-1.1,1.1)
        ax1.spines["left"].set_visible(false)
        ax1.spines["right"].set_visible(false)
        ax1.spines["bottom"].set_visible(false)
        ax1.spines["top"].set_visible(false)
        # Remove ticks
        xticks([])
        yticks([])
    end

    # Create right subplot
    subplot(122)
    # Plot sine if selected by showsine = true
    if showsine
        yphi = mag*sin.(psi.+phi)
        plot(psi*180/pi, yphi,
            color=color, linewidth=linewidth, linestyle=linestyle)
    end
    ax2 = gca()
    # First phasor and sine wave diagram must be drawn with add = false
    if !add
        # Scale and tick x-axis
        xlim(0, 360)
        ylim(-1.1, 1.1)
        xticks(collect(90:90:360),
            backgroundcolor=backgroundcolor)
        if maglabel != ""
            yticks([-mag, 0, mag],[L"$-$"*maglabel, L"$0$", maglabel],
                backgroundcolor=backgroundcolor)
        else
            yticks([-mag, 0, mag],["","",""])
        end
        # Create arrows and labels of axes
        arrowaxes(xlabel=xlabel, ylabel=ylabel)
    else
        # If additional phasor and sine wave diagram are added, further
        # y-axis ticks have to beee added; for this purpose, the existing (old)
        # ticks are stored
        yticks_old = ax2.get_yticks()
        # Store existing (old) labels
        ytickslabel_old = ax2.get_yticklabels()
        # Extend old ticks and labels by addition ticks and labels
        if maglabel != ""
            yticks(cat(yticks_old, [-mag,mag],dims=1),
                cat(ytickslabel_old, [L"$-$"*maglabel,maglabel],dims=1),
                    backgroundcolor=backgroundcolor)
        end
    end
    # Change vertical scaling of left subplot to match
    ax1.set_ylim(ax2.get_ylim())

    # Plot dashed lines, if showdashline = true
    if showdashline
        # Dotted line from phasor arrow to begin of sine wave, split in two
        # parts, to avoid overlay effects if multiple dashed lines are drawn
        con = matplotlib.patches."ConnectionPatch"(
            xyB=(360, mag*sin(phi)), xyA=(0, mag*sin(phi)),
            coordsA="data", coordsB="data",
            axesA=ax2, axesB=ax2, color=colorDash,
            linewidth=lineWidth4, linestyle=":", clip_on=false)
        ax2.add_artist(con)
        con = matplotlib.patches."ConnectionPatch"(
            xyA=(0, mag*sin(phi)), xyB=(mag*cos(phi), mag*sin(phi)),
            coordsA="data", coordsB="data",
            axesA=ax2, axesB=ax1, color=colorDash,
            linewidth=lineWidth4, linestyle=":", clip_on=false)
        ax2.add_artist(con)
        # Dotted line of y-axis of right diagram to maximum of sine wave
        con = matplotlib.patches."ConnectionPatch"(
            xyB=(mod(90-phi*180/pi,360), mag), xyA=(0, mag),
            coordsA="data", coordsB="data",
            axesA=ax2, axesB=ax2, color=colorDash,
            linewidth=lineWidth4, linestyle=":", clip_on=false)
        ax2.add_artist(con)
        # Dotted line of y-axis of right diagram to minimum of sine wave
        con = matplotlib.patches."ConnectionPatch"(
            xyB=(mod(270-phi*180/pi, 360),-mag), xyA=(0, -mag),
            coordsA="data", coordsB="data",
            axesA=ax2, axesB=ax2, color=colorDash,
            linewidth=lineWidth4, linestyle=":", clip_on=false)
            ax2.add_artist(con)
    end
end

@doc raw"""
# Function call

```
angulardimension(r = 1, phi1 = 0, phi2 = pi/2; origin = 0.0im,
    label= "", labelphisep = 0.5, labelrsep = 0.1,
    labelrelrot = false, labelrelangle = 0, ha = "center", va = "center",
    color="black", backgroundcolor="none",
    arrowstyle1 = ".", arrowstyle2 = "-|>", dot90 = false,
    linewidth = 0.6, linestyle = "-", width = 0.2,
    headlength = 5, headwidth = 2.5)
```

# Description

This function draws an arrowed arc, intended to label the angle between phasors.
The arc is drawn from angle `phi1` (begin) to `phi2` (end). The begin and end
arrow shapes can be set. The arc can be labeled and optionally a dot maker can
be used to indicate right angles (90°).

# Variables

`r` Radius if the arc with not unit; shall be between 0 and 1; default value = 1

`phi1` Phase angle of the begin of the arc; default value = 0

`phi2` Phase angle of the end of the arc; default value = pi/2

`origin` Complex quantity, indicating the origin of the arc with no unit;
default value = 0.0 + 0.0im

`label` Label of the angle; default value =""

`labelphisep` Angular separation of the label with respect to the arc; if
`labelphisep == 0`, the label is located at angle `phi1` and if `labelphisep ==
1`, the label is located at angle `phi2`; default value = 0.5, right in the
middle between `phi1` and `phi2`

labelrsep` Radial per unit location of label (in direction of the phasor):
`labelrsep = 0` locates the label right on the arc. A positive value locates the
label outside the arc, a negative value locates the label inside the arc;
default value = 0.1

`labelrelrot` Relative rotation of label; if `labelrelrot == false` (default
value) then the label is not rotated relative to the center of the arc;
otherwise the label is rotated relative to the  angle `labelrelangle`

`labelrelangle` Relative angle of label with respect to center of the arc; this
angle is only applied, it `labelrelrot == true`; this angle indicates the
relative orientation of the label with respect to the center of the arc;
default value = 0

`ha` Horizontal alignment of label; actually represents the radial alignment of
the label; default value = "center"

`va` Vertical alignment of label; actually represents the tangential alignment
of the label; default value = "center"

`color` Color of the arc; default value = "black"

`backgroundcolor` Background color of the label; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`arrowstyle1` Arrow style of the begin of the arc; default value = "."; valid
strings are:

- `.` dot marker
- `<|-` arrow

`arrowstyle2` Arrow style of the end of the arc; default value = "-|>"; valid
strings are:

- `.` dot marker
- `-|>` arrow

`headlength` Length of arrow head; default value = 5

`headwidth` Width of arrow head; default value = 2.5

# Example

Copy and paste code:

```julia
using Unitful, Unitful.DefaultSymbols, PyPlot, ElectricalEngineering
figure(figsize=(3.3, 2.5))
rc("text", usetex=true); rc("font", family="serif")
V1 = 100V + j*0V # Voltage
Z1 = 30Ω + j*40Ω # Impedance
I1 = V1/Z1 # Current
Vr = real(Z1)*I1
Vx = V1 - Vr
refV = abs(V1); refI=abs(I1)*0.8
phasor(V1, label=L"$\underline{V}_1$", labeltsep=0.1, ref=refV,
    labelrelrot=true)
phasor(Vr, label=L"$\underline{V}_r$", labeltsep=-0.25, ref=refV,
    labelrelrot=true)
phasor(Vx, origin=Vr, label=L"$\underline{V}_x$", labeltsep=-0.2, ref=refV,
    labelrelrot=true)
phasor(I1, label=L"$\underline{I}_1$", labeltsep=-0.2, labelrsep=0.7, ref=refI,
    labelrelrot=true, linestyle="--", par=-0.05)
phi1=angle(I1)
phi2=angle(V1)
angulardimension(0.3,phi1,phi2,arrowstyle1=".",arrowstyle2="-|>",ha="left",
    label=L"$\varphi_1$", labelrsep=0.05)
axis("square"); xlim(-1,1); ylim(-1,1)
removeaxes(); # Remove axis
# save3fig("phasordiagram",crop=true);
```
"""
function angulardimension(r = 1,
    phi1 = 0,
    phi2 = pi/2;
    origin = 0.0im,
    label = "",
    labelphisep = 0.5,
    labelrsep = 0.1,
    labelrelrot = false,
    labelrelangle = 0,
    ha = "center",
    va = "center",
    color="black",
    backgroundcolor="none",
    arrowstyle1 = ".",
    arrowstyle2 = "-|>",
    dot90 = false,
    linewidth = 0.6,
    linestyle = "-",
    width = 0.2,
    headlength = 5,
    headwidth = 2.5)

    # Determine sign of the difference of angles in order to revert oriention
    # angle, if required
    sig = sign(phi2-phi1)
    # Number of segments based on every 2 degrees
    segs = Int(round((phi2-phi1)*180/pi/    2)*sig)
    # Create vector of angles to draw arc
    phi = collect(range(phi1,stop=phi2,length=segs))
    # Arc coordinates
    x = upstrip(r)*cos.(phi) + fill(real(upstrip(origin)),length(phi))
    y = upstrip(r)*sin.(phi) + fill(imag(upstrip(origin)),length(phi))
    # Draw arc
    plot(x,y,color=color,linewidth=linewidth,linestyle=linestyle)

    # Draw arrows at begin and end of arc, depending on the specified arrow
    # styles

    # Begin of arc
    if arrowstyle1 == "<|-"
        # Arrow head
        annotate("", xy=(x[1], y[1]),
            xytext=(x[1]-sig*0.01*r*cos(phi1-pi/2*0.98),
                y[1]-sig*0.01*r*sin(phi1-pi/2*0.98)),
            xycoords="data",
            arrowprops=Dict("edgecolor"=>color, "facecolor"=>color,
                "width"=>width, "linestyle"=>"-",
                "headlength"=>headlength,
                "headwidth"=>headwidth),
            annotation_clip=false)
    elseif arrowstyle1 == "."
        # Dot marker
        plot(x[1],y[1],marker=".",color=color, clip_on=false)
    end

    # End of arc
    if arrowstyle2 == "-|>"
        # Arrow head
        annotate("", xy=(x[end], y[end]),
            xytext=(x[end]-sig*0.01*r*cos(phi2+pi/2*0.98),
                y[end]-sig*0.01*r*sin(phi2+pi/2*0.98)),
            xycoords="data",
            arrowprops=Dict("edgecolor"=>color, "facecolor"=>color,
                "width"=>width, "linestyle"=>"-",
                "headlength"=>headlength,
                "headwidth"=>headwidth),
            annotation_clip=false)
    elseif arrowstyle2 == "."
        # Marker dot
        plot(x[end],y[end],marker="." ,color=color, clip_on=false)
    end

    # Plot label
    # Angle in the middle of the arc
    phim = (phi1+phi2)/2.0
    # Angle difference
    dphi = phi2-phi1
    # Radial position of text
    rlabel = r + labelrsep
    # Angular position of text
    philabel = phi1 + dphi*labelphisep
    # Rectangular position of text
    if labelrelrot == false
        # Without relative rotation of label
        text(rlabel*cos(philabel) + real(origin),
            rlabel*sin(philabel) + imag(origin),
            label, ha=ha, va=va, rotation=labelrelangle*180/pi,
            backgroundcolor=backgroundcolor)
    else
        # Applying relative rotation of label
        text(rlabel*cos(philabel) + real(origin),
            rlabel*sin(philabel) + imag(origin),
            label, ha=ha, va=va, rotation=(phim+labelrelangle)*180/pi,
            backgroundcolor=backgroundcolor)
    end

    # Optionally create a dot in the center of the arc to indicate 90°
    if dot90
        plot(r/2*cos(phim) + real(origin),
            r/2*sin(phim) + imag(origin),
            marker=".", color=color, clip_on=false)
    end
end

@doc raw"""
# Function call

```
phasordimension(c; origin = 0im, ref = 1;
    label = "", labeltsep = 0.1, labelrsep=0.5, labelrelrot=false,
    labelrelangle=0, ha = "center", va = "center",
    color="black", backgroundcolor = "none",
    arrowstyle1 = "<|-", arrowstyle2 = "-|>", linewidth=0.6, linestyle = "-",
    width=0.2, headlength=5, headwidth=2.5,
    par=0, paroverhang = 0.02, parcolor = "black",
    parlinewidth=0.6, parlinestyle = "-")

```

# Description

This function draws length dimension with label based on rectangular x- and
y-coordinates. The length dimension arrow can be shifted parallel to the
specified  coordinates by means the single parameter `par`. Additionally,
auxiliary lines  from the specified coordinates to the length dimension are
drawn.

# Variables

`c` Complex phasor to be dimensioned

`origin` Complex number representing the origin of the phasor dimension; this
variable needs to have the same unit as `c`

`ref` Reference length of scaling; this is required as in a phasor diagram
voltages and currents may be included; in order to account for the different
voltage and current scales, one (constant)  `ref` is used for voltage
phasors and another (constant) `ref` is used for current phasors; this
variable needs to have the same unit as `c`

`label` Label of the angle; default value =""

`labelphisep` Angular separation of the label with respect to the arc; if
`labelphisep == 0`, the label is located at angle `phi1` and if `labelphisep ==
1`, the label is located at angle `phi2`; default value = 0.5, right in the
middle between `phi1` and `phi2`

labelrsep` Radial per unit location of label (in direction of the phasor):
`labelrsep = 0` locates the label right on the arc. A positive value locates the
label outside the arc, a negative value locates the label inside the arc;
default value = 0.1

`labelrelrot` Relative rotation of label; if `labelrelrot == false` (default
value) then the label is not rotated relative to the center of the arc;
otherwise the label is rotated relative to the  angle `labelrelangle`

`labelrelangle` Relative angle of label with respect to center of the arc; this
angle is only applied, it `labelrelrot == true`; this angle indicates the
relative orientation of the label with respect to the center of the arc;
default value = 0

`ha` Horizontal alignment of label; actually represents the radial alignment of
the label; default value = "center"

`va` Vertical alignment of label; actually represents the tangential alignment
of the label; default value = "center"

`color` Color of the arc; default value = "black"

`backgroundcolor` Background color of the label; if labelrsep is equal to 0, the
background color "white" can be used; default value = "none"

`arrowstyle1` Arrow style of the begin of the arc; default value = "."; valid
strings are:

- `.` dot marker
- `<|-` arrow

`arrowstyle2` Arrow style of the end of the arc; default value = "-|>"; valid
strings are:

- `.` dot marker
- `-|>` arrow

`headlength` Length of arrow head; default value = 5

`headwidth` Width of arrow head; default value = 2.5

`par` In order to be able to draw the length dimension parallel to the specified
coordiantes `x1`, `y1`, `x2`, `y2`, `par` is used to specify the per unit
tangential shift (offset) of the length dimension; default value = 0 (no shift)

`paroverhang` The auxiliary lines, drawn in case of `par != 0`, show the absolute
overhang `paroverhang`; default value = 0.02

`parcolor` Color of the auxiliary lines; default value = "black"

`parlinewidth` Line width of the auxiliary line; default value = 0.6

`parlinestyle` Line style of the auxiliary line; default value = "-"

# Examples

Copy and paste the following code:

```julia
using Unitful, Unitful.DefaultSymbols, PyPlot, ElectricalEngineering
figure(figsize=(3.3, 2.5))
rc("text", usetex=true); rc("font", family="serif")
Z1 = pol(1,30°)
phasor(Z1, label=L"$\underline{Z}$", labeltsep = 0.05, labelrelrot=true)
phasordimension(real(Z1), label=L"$R$", arrowstyle1="",
    linestyle="-", linewidth=1, headwidth=5, headlength=10,
    labeltsep=0, color="gray", backgroundcolor="white")
    phasordimension(j*imag(Z1), origin=real(Z1), label=L"j$\cdot X$",
    arrowstyle1="", linestyle="-", linewidth=1, headwidth=5, headlength=10,
    labeltsep=0, color="gray", backgroundcolor="white")
arrowaxes(xmin = real(Z1)+0.1, xlabel="Re", ylabel=L"j$\cdot$Im")
xlim([-0.5,1]);ylim([-0.5,1]); axis("square"); removeaxes();
# save3fig("phasordimension", crop=true)
```
"""
function phasordimension(c;
    origin = (0.0+0.0im).*c./ustrip(c),
    ref = abs(c./ustrip(c)),
    label = "",
    labeltsep = 0.1,
    labelrsep = 0.5,
    labelrelrot = false,
    labelrelangle = 0,
    ha = "center",
    va = "center",
    color="black",
    backgroundcolor = "none",
    arrowstyle1 = "<|-",
    arrowstyle2 = "-|>",
    linewidth = 0.6,
    linestyle = "-",
    width = 0.2,
    headlength = 5,
    headwidth = 2.5,
    par = 0,
    paroverhang = 0.02,
    parcolor = "black",
    parlinewidth = 0.6,
    parlinestyle = "-")

    lengthdimension(real(origin/ref), imag(origin/ref),
        real((c+origin)/ref), imag((c+origin)/ref),
        label=label, labeltsep = labeltsep, labelrsep=labelrsep,
        labelrelrot=labelrelrot, labelrelangle = labelrelangle, ha=ha, va=va,
        color=color, backgroundcolor = backgroundcolor,
        arrowstyle1=arrowstyle1, arrowstyle2=arrowstyle2,
        linewidth = linewidth, linestyle=linestyle, width=width,
        headlength=headlength, headwidth=headwidth,
        par=par, paroverhang=paroverhang, parcolor = parcolor,
        parlinewidth=parlinewidth, parlinestyle=parlinestyle)
end


text:
@testset "Element stiffness matrix" begin

    # 1x1 linear quadrilateral cell
    # Isoparameteric
    # 2nd order quadrature rule

    E = 1
    nu = 0.3
    k = [1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8, -1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8]
    Ke = E/(1-nu^2)*[k[1] k[2] k[3] k[4] k[5] k[6] k[7] k[8];
                     k[2] k[1] k[8] k[7] k[6] k[5] k[4] k[3];
                     k[3] k[8] k[1] k[6] k[7] k[4] k[5] k[2];
                     k[4] k[7] k[6] k[1] k[8] k[3] k[2] k[5];
                     k[5] k[6] k[7] k[8] k[1] k[2] k[3] k[4];
                     k[6] k[5] k[4] k[3] k[2] k[1] k[8] k[7];
                     k[7] k[4] k[5] k[2] k[3] k[8] k[1] k[6];
                     k[8] k[3] k[2] k[5] k[4] k[7] k[6] k[1]]

    problem = HalfMBB((2,2), (1.,1.), 1., 0.3, 1.);

    #@test ElementFEAInfo(problem).Kes[1] ≈ Ke
end


text:
# Parent file to run all scripts which may generate
# some output that you want to display on the website.
# this can be used as a tester to check that all the code
# on your website runs properly.

dir = @__DIR__

"""
    genplain(s)

Small helper function to run some code and redirect the output (stdout) to a file.
"""
function genplain(s::String)
    open(joinpath(dir, "output", "$(splitext(s)[1]).out"), "w") do outf
        redirect_stdout(outf) do
            include(joinpath(dir, s))
        end
    end
end
genplain("geoFieldCountries.jl")
genplain("geoCoastlinesField.jl")
genplain("geoCoastlines.jl")
genplain("geoCoastlinesStatesUS.jl")

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Clang, Clang.cindex
# We'll replace some of the wrap functions, and we need to call the following:
import Clang.wrap_c: repr_jl, rep_type, rep_args, name_safe

context=wrap_c.init(output_file="gen_libcufft.jl",
                    common_file="gen_libcufft_h.jl",
                    header_library=x->"libcufft",
                    clang_includes=["/usr/include"],
                    header_wrapped=(x,y)->(contains(x,"cuda")))

# context.options = wrap_c.InternalOptions(true)  # wrap structs, too
path = "/usr/local/cuda-5.0/include"
headers = [joinpath(path,"cufft.h")]

# Customize the wrap function for functions. This was copied
# from Clang/src/wrap_c.jl, with the following customizations:
#   - error-check each function that returns a cudaError_t
#   - omit types from function prototypes
function wrap_c.wrap(buf::IO, funcdecl::FunctionDecl, libname::ASCIIString)
    function print_args(buf::IO, cursors, types)
        i = 1
#        for (c,t) in zip(cursors,types)
#            print(buf, name_safe(c), "::", t)
#            (i < length(cursors)) && print(buf, ", ")
#            i += 1
#        end
        for c in cursors
            print(buf, name_safe(c))
            (i < length(cursors)) && print(buf, ", ")
            i += 1
        end
    end

    cu_spelling = spelling(funcdecl)

    funcname = spelling(funcdecl)
    arg_types = cindex.function_args(funcdecl)
    args = [x for x in search(funcdecl, ParmDecl)]
    arg_list = tuple( [repr_jl(x) for x in arg_types]... )
    ret_type = repr_jl(return_type(funcdecl))

    print(buf, "function ")
    print(buf, spelling(funcdecl))
    print(buf, "(")
    print_args(buf, args, [myrepr_jl(x) for x in arg_types])
    println(buf, ")")
    print(buf, "  ")
    ret_type == "cufftResult" && print(buf, "checkerror(")
    print(buf, "ccall( (:", funcname, ", ", libname, "), ")
    print(buf, rep_type(ret_type))
    print(buf, ", ")
    print(buf, rep_args(arg_list), ", ")
    for (i,arg) in enumerate(args)
        print(buf, name_safe(arg))
        (i < length(args)) && print(buf, ", ")
    end
    ret_type == "cufftResult" && print(buf, ")")
    println(buf, ")")
    println(buf, "end")
end

function myrepr_jl(x)
    str = repr_jl(x)
    return (str == "Ptr{Cint}") ? "Array{Cint}" : str
end

wrap_c.wrap_c_headers(context, headers)


text: datadir = "Astrobee"


text:
module Queries
  using Catlab: @present
  using Catlab.Programs.RelationalPrograms
  using Catlab.Programs.RelationalPrograms: TheoryTypedRelationDiagram, TheoryTypedNamedRelationDiagram, TypedNamedRelationDiagram, parse_relation_diagram
  using Catlab.Graphics
  using Catlab.WiringDiagrams
  using Catlab.CategoricalAlgebra.CSets
  using ..DB

  export TheoryQuery, Query, @query, to_sql, draw_query, to_prepared_sql, infer!

  const SQLOperators = Dict(:<    => ("<", [:first, :second]),
                            :>    => (">", [:first, :second]),
                            :(==) => ("=", [:first, :second]),
                            :<=   => ("<=", [:first, :second]),
                            :>=   => (">=", [:first, :second]),
                            :(!=) => ("<>", [:first, :second]),
                           )

  @present TheoryQuery <: TheoryTypedNamedRelationDiagram begin
    Comparison::Ob
    comp_port1::Hom(Comparison, Port)
    comp_port2::Hom(Comparison, Port)
    # comp_port1⋅port_type == comp_port2⋅port_type
    # port_type == junction⋅junction_type
    # subpart(q, :port_type) == subpart(q, subpart(q, :junction), :junction_type)
  end

  @acset_type Query(TheoryQuery,
                    index=[:box, :junction, :outer_junction, :port_name],
                    unique_index=[:variable])

  NullableSym = Union{Symbol, Nothing}
  Query() = Query{NullableSym, NullableSym, NullableSym}()

  function infer!(wd, rels::Array{Tuple{Array{Symbol,1}, Array{Symbol,1}},1}; max_iter=2*length(rels))
    # Perform multiple steps to fill in chains of relations
    for i in 1:max_iter
      inf_step = [infer_indiv!(wd, rel) for rel in rels]
      if all([!s[2] for s in inf_step])
        if all([s[1] for s in inf_step])
          # Successfully defined all values to be defined
          return true
        end
        # Did not define all values, but stopped changing
        throw(ErrorException("Not all relations were able to be defined. Original system insufficiently defined"))
      end
    end
    return false
  end

  function infer_indiv!(wd, rel::Tuple{Array{Symbol}, Array{Symbol}})
    a = subpart(wd, rel[1][1])

    for i in 2:length(rel[1])
      if i == length(rel[1])
        a = view(subpart(wd, rel[1][i]), a)
      else
        a = subpart(wd, a, rel[1][i])
      end
    end

    b = subpart(wd, rel[2][1])
    for i in 2:length(rel[2])
      if i == length(rel[2])
        b = view(subpart(wd, rel[2][i]), b)
      else
        b = subpart(wd, b, rel[2][i])
      end
    end
    changed = false
    length(a) == length(b) || throw(DimensionMismatch("The arguments of $(rel[1])==$(rel[2]) have different lengths, $(length(a)) and $(length(b))"))
    is_defined = map(1:length(a)) do i
      if a[i] == nothing
        if b[i] == nothing
          return false
        else
          a[i] = b[i]
          changed = true
        end
      else
        if b[i] == nothing
          b[i] = a[i]
          changed = true
        else
          a[i] == b[i] || throw(DimensionMismatch("The arguments of $(rel[1])==$(rel[2]) have inconsistent values at index $i, $(a[i]) and $(b[i])"))
        end
      end
      true
    end
    return all(is_defined), changed
  end

  function Query(schema, wd)
    port_names = get_fields(schema)
    q = Query()
    copy_parts!(q, wd)

    # Add types to each of the junctions by iterating through each box
    for b in 1:nparts(q, :Box)
      # Set junction and outer_port types (these will be inferred from schema types)
      q[:outer_port_type] .= nothing
      q[:outer_port_name] .= nothing
      q[:junction_type] .= nothing
      name = q[b, :name]

      ports = incident(wd, b, :box)
      # add comparison references for later type-inference
      if name in keys(SQLOperators)
        add_part!(q, :Comparison, comp_port1=ports[1], comp_port2=ports[2])
        set_subparts!(q, ports, port_name=SQLOperators[name][2][1:2], port_type=[nothing, nothing])
      else
        fields = [port_names[name][i][1] for i in 1:length(ports)]
        types = [Symbol(port_names[name][i][2]) for i in 1:length(ports)]
        set_subparts!(q, ports, port_name=fields, port_type=types)
      end
    end
    infer!(q, [([:port_type],[:junction, :junction_type]),
               ([:outer_junction, :junction_type],[:outer_port_type]),
               ([:comp_port1,:port_type],[:comp_port2,:port_type]),
               ([:outer_port_name],[:outer_junction, :variable])]);
    q
  end

  macro query(schema, exprs...)
    Expr(:call, GlobalRef(Queries, :parse_query_statement),
                esc(schema), exprs)
  end

  function parse_query_statement(schema, exprs)
    wd = parse_relation_diagram((expr for expr in exprs)...)
    Query(schema, wd)
  end

  function port_indices(wd)
    box_sizes = zeros(Int, nparts(wd, :Box))
    map(subparts(wd, [:box])) do b
      box_sizes[b[1]] += 1
      return box_sizes[b[1]]
    end
  end

  # Generates a new array where all strings are unique. The new strings are based
  # off of the strings provided in the original array.
  uniquify(a::Array{String,1}) = begin
    a_n = Array{String,1}()
    # Fill a_n with unique values
    for i in 1:length(a)
      cur_val = a[i]

      # Iterate the digit until the value is unique
      while cur_val in a_n
        if occursin(r"_\d*$", cur_val)
          cur_val = replace(cur_val, r"_\d*$" => (c) -> "_"*string(parse(Int, c[2:end])+1))
        else
          cur_val = cur_val*"_0"
        end
      end
      append!(a_n, [cur_val])
    end
    return a_n
  end

  function to_sql(q::Query)

    box_names = subpart(q, :name)
    table_inds = findall(x -> !(x in keys(SQLOperators)), box_names)
    op_inds = findall(x -> x in keys(SQLOperators), box_names)

    outer_juncs = subpart(q, :outer_junction)
    port_info = subparts(q, [:box, :junction, :port_name])
    junctions = zeros(Int, nparts(q, :Junction))
    variables = subpart(q, :variable)
    prepared_junctions = findall(x -> string(x)[1] == '_', variables)

    junctions[prepared_junctions] = nparts(q, :Port) .+ (1:length(prepared_junctions))

    field_name(port) = begin
      if port > length(port_info)
        "\$$(port - length(port_info))"
      else
        p = port_info[port]
        "t$(p[1]).$(p[3])"
      end
    end

    # Construct the conditions by iterating through each port
    conditions = Array{String,1}()
    for (i, p) in enumerate(port_info)
      if p[1] in table_inds
        prev_pind = junctions[p[2]]
        if prev_pind == 0
          junctions[p[2]] = i
        else
          push!(conditions, "$(field_name(i))=$(field_name(prev_pind))")
        end
      end
    end

    # Construct the operator relations
    append!(conditions, map(op_inds) do i
              p_ind = junctions[subpart(q, incident(q,i,:box), :junction)]
              op = SQLOperators[box_names[i]][1]
              "$(field_name(p_ind[1]))$(op)$(field_name(p_ind[2]))"
            end)

    # Construct the aliases for each table
    alias = map(table_inds) do i "$(box_names[i]) AS t$i" end

    # Construct the field output selection statement
    selector_alias = uniquify(map(i -> "$(variables[i])",outer_juncs))
    selectors = map(enumerate(outer_juncs)) do (i,j)
      "$(field_name(junctions[j])) AS $(selector_alias[i])"
    end

    query = "SELECT $(join(selectors, ", "))\nFROM $(join(alias, ", "))"

    if length(conditions) > 0
      query = string(query, "\nWHERE $(join(conditions, " AND "))")
    end

    query
  end

  function to_prepared_sql(q::Query, uid::String)
    junc_types = map(x -> typeToSQL(x), subpart(q, :junction_type))
    prepared_junctions = findall(x -> string(x)[1] == '_', subpart(q, :variable))

    return "PREPARE \"$uid\" ($(join(junc_types[prepared_junctions], ","))) AS\n$(to_sql(q))",
           length(prepared_junctions)
  end

  function draw_query(q; kw...)
    uwd = TypedNamedRelationDiagram{NullableSym, NullableSym, NullableSym}()
    copy_parts!(uwd, q)
    to_graphviz(uwd; box_labels=:name, junction_labels=:variable, kw...)
  end

  # Replication of CSet functionality
  # TODO: Find best way to copy objects and attributes between CSets of
  #       different types
  function subparts(acs::ACSet, names::Array{Symbol,1})
    collect(zip([subpart(acs, name) for name in names]...))
  end
end


text:
function legendre_next(x::Float64, j::Int64, p1::Float64, p0::Float64; tol::Float64=1e-4)
    if x < -1.0 - tol || x > 1.0 + tol
        throw(DomainError("x lies outside [- 1 - $tol, 1 + $tol]"))
    elseif tol < 0.0
        throw(ArgumentError("tol must be non-negative"))
    end
    if j == 0 || j == 1
        return j == 0 ? 1.0 : x
    end
    x = min(x, 1.0)
    x = max(x, -1.0)
    return (2 * j - 1) * x * p1 / j - (j - 1) * p0 / j
end

"""
    Index(dim::Vector{Int64}, deg::Vector{Int64})

Multivariate monomial index.

The monomial `x[1] * x[3]^2`` can be encoded using `dim = [1, 3]`, `deg = [1, 2]`
"""
struct Index
    dim::Vector{Int64}
    deg::Vector{Int64}

    function Index(dim::Vector{Int64}, deg::Vector{Int64})
        if length(unique(dim)) != length(dim)
            throw(ArgumentError("dim contains repeated values"))
        elseif !all(dim .> 0)
            throw(ArgumentError("dim contains values < 1"))
        elseif !all(deg .>= 0)
            throw(ArgumentError("deg contains values < 0"))
        elseif length(dim) != length(deg)
            throw(DimensionMismatch("dim and deg have mismatched length"))
        else
            return new(dim, deg)
        end
    end
end

udeg(index) = length(index.deg) == 0 ? 0 : maximum(index.deg)

function _fill_tpbasis!(
    d::Int64, J::Int64, basis::Vector{Index}, dim::Vector{Int64}, deg::Vector{Int64}
)
    # Reached the bottom of the recursion. Append the basis.
    if d == 0 || J == 0
        push!(basis, Index(dim, deg))
        return nothing
    end

    # Append the current basis function with every remaining j.
    for j in 0:J
        dim1, deg1 = deepcopy(dim), deepcopy(deg)
        if j > 0
            push!(dim1, d)
            push!(deg1, j)
        end
        _fill_tpbasis!(d - 1, J - j, basis, dim1, deg1)
    end
end

"""
    tpbasis(d::Int64, J::Int64)

Construct the `d`-dimensional truncated tensor-product basis.

All index terms have a degree ≤ `J`.

See also [`Index`](@ref)
"""
function tpbasis(d::Int64, J::Int64)
    if d < 1
        throw(ArgumentError("recieved d < 1"))
    elseif J < 0
        throw(ArgumentError("recieved negative J"))
    end
    basis = Vector{Index}(undef, 0)
    dim, deg = Int64[], Int64[]
    _fill_tpbasis!(d, J, basis, dim, deg)
    return basis
end

function check_limits(limits::Matrix{Float64})
    if size(limits, 2) != 2
        throw(ArgumentError("limits should have exactly 2 columns"))
    elseif !all(limits[:, 1] .<= limits[:, 2])
        throw(ArgumentError("limits[:, 1] must be <= limits[:, 2]"))
    end
end

"""
    expand(X::AbstractMatrix, basis::Vector{Index}, limits::Matrix{Float64};
           J::Union{Nothing,Int64}=nothing)

Expand the columns of X into a rescaled legendre polynomial basis.

# Arguments
- `X::AbstractMatrix`: Matrix with observations stored as columns.
- `basis::Vector{Index}`: Vector of monomial indices.
- `limits::Matrix{Float64}`: Matrix with two columns defining the lower/upper limits of the space.
- `J::Union{Nothing, Int64}=nothing`: The maximum degree of the basis. Inferred if not specified.
"""
function expand(
    X::AbstractMatrix,
    basis::Vector{Index},
    limits::Matrix{Float64};
    J::Union{Nothing,Int64}=nothing,
)
    check_limits(limits)
    if size(limits, 1) != size(X, 1)
        throw(ArgumentError("X and limits don't match in their first dimension"))
    elseif !isnothing(J) && J < 0
        throw(ArgumentError("J should be >= 0"))
    end

    J = isnothing(J) ? maximum(udeg.(basis)) : J # Infer J if not supplied
    d, n = size(X)
    nbas = length(basis)
    Y = ones(nbas, n)
    U = ones(J + 1, d)
    for i in 1:n
        # Univariate expansions
        for l in 1:d
            p1 = p0 = 1.0
            x = (2 * X[l, i] - limits[l, 1] - limits[l, 2]) / (limits[l, 2] - limits[l, 1])
            for j in 0:J
                U[j + 1, l] = legendre_next(x, j, p1, p0)
                p0 = p1
                p1 = U[j + 1, l]
            end
        end

        # Mulitplicative combinations
        for b in 1:nbas
            for l in 1:length(basis[b].dim)
                ud = basis[b].dim[l]
                uj = basis[b].deg[l]
                Y[b, i] *= U[uj + 1, ud]
            end
        end
    end
    return Y
end

"""
    BayesPM(basis::Vector{Index}, limits::Matrix{Float64}; <keyword arguments>)

Construct a Bayesian linear model on polynomial features.

# Arguments
- `basis::Vector{Index}`: Vector of monomial indices.
- `limits::Matrix{Float64}`: Matrix with two columns defining the lower/upper limits of the space.
- `λ::Float64=1.0`: Prior covariance scale factor.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e=3`: Inverse-gamma prior scale hyperparameter.
"""
mutable struct BayesPM <: AbstractBayesianLM
    J::Int64
    basis::Vector{Index}
    limits::Matrix{Float64}
    lm::BayesLM

    function BayesPM(
        basis::Vector{Index},
        limits::Matrix{Float64};
        λ::Float64=1.0,
        shape0::Float64=1e-3,
        scale0::Float64=1e-3,
    )
        check_limits(limits)
        lm = BayesLM(length(basis); λ=λ, shape0=shape0, scale0=scale0)
        J = maximum(udeg.(basis))
        return new(J, basis, limits, lm)
    end
end

function fit!(pm::BayesPM, X::AbstractMatrix, y::AbstractMatrix)
    check_regression_data(X, y)
    Z = expand(X, pm.basis, pm.limits; J=pm.J)
    return fit!(pm.lm, Z, y)
end

function (pm::BayesPM)(X::AbstractMatrix)
    Z = expand(X, pm.basis, pm.limits; J=pm.J)
    return pm.lm(Z)
end

function shape_scale(pm::BayesPM)
    return pm.lm.shape, pm.lm.scale
end


text:
export PendingSubject, PendingSubjectFactory

import Base: show, similar

"""
    PendingSubject(::Type{D}) where D
    PendingSubject(::Type{D}, factory::F) where { D, F <: AbstractFactory }
    PendingSubject(::Type{D}, subject::S) where { D, S }

A variant of Subject that emits its last value on completion.
Reemits last value on further subscriptions and then completes.

See also: [`PendingSubjectFactory`](@ref), [`Subject`](@ref), [`SubjectFactory`](@ref)
"""
function PendingSubject end

"""
    PendingSubjectFactory(factory::F) where { F <: AbstractSubjectFactory }
    PendingSubjectFactory(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

A variant of SubjectFactory that creates an instance of PendingSubject.

See also: [`SubjectFactory`](@ref), [`AbstractSubjectFactory`](@ref), [`PendingSubject`](@ref), [`Subject`](@ref)
"""
function PendingSubjectFactory end

##

mutable struct PendingSubjectProps{D}
    last :: Union{Nothing, D}

    PendingSubjectProps{D}() where D = new(nothing)
end

struct PendingSubjectInstance{D, S} <: AbstractSubject{D}
    subject :: S
    props   :: PendingSubjectProps{D}
end

Base.show(io::IO, ::PendingSubjectInstance{D, S}) where { D, S } = print(io, "PendingSubject($D, $S)")

Base.similar(subject::PendingSubjectInstance{D, S}) where { D, S } = PendingSubject(D, similar(subject.subject))

function PendingSubject(::Type{D}) where D
    return PendingSubject(D, SubjectFactory(AsapScheduler()))
end

function PendingSubject(::Type{D}, factory::F) where { D, F <: AbstractSubjectFactory }
    return PendingSubject(D, create_subject(D, factory))
end

function PendingSubject(::Type{D}, subject::S) where { D, S }
    return as_pending_subject(D, as_subject(S), subject)
end

as_pending_subject(::Type{D},  ::InvalidSubjectTrait,    subject)    where D          = throw(InvalidSubjectTraitUsageError(subject))
as_pending_subject(::Type{D1}, ::ValidSubjectTrait{D2},  subject)    where { D1, D2 } = throw(InconsistentSubjectDataTypesError{D1, D2}(subject))
as_pending_subject(::Type{D},  ::ValidSubjectTrait{D},   subject::S) where { D, S }   = PendingSubjectInstance{D, S}(subject, PendingSubjectProps{D}())

##

getlast(subject::PendingSubjectInstance)         = subject.props.last
setlast!(subject::PendingSubjectInstance, value) = subject.props.last = value

iscompleted(subject::PendingSubjectInstance) = iscompleted(subject.subject)
isfailed(subject::PendingSubjectInstance)    = isfailed(subject.subject)

##

function on_next!(subject::PendingSubjectInstance{D}, data::D) where D
    if !iscompleted(subject) && !isfailed(subject)
        setlast!(subject, data)
    end
end

function on_error!(subject::PendingSubjectInstance, err)
    if !iscompleted(subject) && !isfailed(subject)
        error!(subject.subject, err)
    end
end

function on_complete!(subject::PendingSubjectInstance)
    if !iscompleted(subject) && !isfailed(subject)
        last = getlast(subject)
        if last !== nothing
            next!(subject.subject, last)
        end
        complete!(subject.subject)
    end
end

##

function on_subscribe!(subject::PendingSubjectInstance, actor)
    if iscompleted(subject)
        last = getlast(subject)
        if last !== nothing
            next!(actor, last)
        end
        complete!(actor)
        return voidTeardown
    else
        return subscribe!(subject.subject, actor)
    end
end

##

struct PendingSubjectFactoryInstance{ F <: AbstractSubjectFactory } <: AbstractSubjectFactory
    factory :: F
end

Base.show(io::IO, subject::PendingSubjectFactoryInstance{F}) where F = print(io, "PendingSubjectFactoryInstance($F)")

create_subject(::Type{L}, factory::PendingSubjectFactoryInstance) where L = PendingSubject(L, factory.factory)

function PendingSubjectFactory(factory::F) where { F <: AbstractSubjectFactory }
    return PendingSubjectFactoryInstance(factory)
end

function PendingSubjectFactory(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }
    return PendingSubjectFactoryInstance(SubjectFactory{H}(scheduler))
end


text: const Caching = String


text:
# PyPlot initialization — the hardest part is finding a working backend.
using VersionParsing

###########################################################################

# global PyObject constants that get initialized at runtime.  We
# initialize them here (rather than via "global foo = ..." in __init__)
# so that their type is known at compile-time.

const matplotlib = PyNULL()
const plt = PyNULL()
const Gcf = PyNULL()
const orig_draw = PyNULL()
const orig_gcf = PyNULL()
const orig_figure = PyNULL()
const orig_show = PyNULL()

###########################################################################
# file formats supported by Agg backend, from MIME types
const aggformats = Dict("application/eps" => "eps",
                        "image/eps" => "eps",
                        "application/pdf" => "pdf",
                        "image/png" => "png",
                        "application/postscript" => "ps",
                        "image/svg+xml" => "svg")

# In 0.6, TextDisplay can show e.g. image/svg+xml as text (#281).
# Any "real" graphical display should support PNG, I hope.
isdisplayok() = displayable(MIME("image/png"))

###########################################################################
# We allow the user to turn on or off the Python gui interactively via
# pygui(true/false).  This is done by loading pyplot with a GUI backend
# if possible, then switching to a Julia-display backend (if available)

# like get(dict, key, default), but treats a value of "nothing" as a missing key
function getnone(dict, key, default)
    ret = get(dict, key, default)
    return ret === nothing ? default : ret
end

# return (backend,gui) tuple
function find_backend(matplotlib::PyObject)
    gui2matplotlib = Dict(:wx=>"WXAgg",:gtk=>"GTKAgg",:gtk3=>"GTK3Agg",
                          :qt_pyqt4=>"Qt4Agg", :qt_pyqt5=>"Qt5Agg",
                          :qt_pyside=>"Qt4Agg", :qt4=>"Qt4Agg",
                          :qt5=>"Qt5Agg", :qt=>"Qt4Agg",:tk=>"TkAgg")
    conda = PyCall.conda || !isempty(PyCall.anaconda_conda())
    if Sys.islinux()
        guis = [:tk, :gtk3, :gtk, :qt5, :qt4, :wx]
    elseif Sys.isapple()
        guis = [:qt5, :qt4, :wx, :gtk, :gtk3, :tk] # issue #410
    else
        guis = [:tk, :qt5, :qt4, :wx, :gtk, :gtk3]
    end
    options = [(g,gui2matplotlib[g]) for g in guis]

    matplotlib2gui = Dict("wx"=>:wx, "wxagg"=>:wx,
                          "gtkagg"=>:gtk, "gtk"=>:gtk,"gtkcairo"=>:gtk,
                          "gtk3agg"=>:gtk3, "gtk3"=>:gtk3,"gtk3cairo"=>:gtk3,
                          "qt5agg"=>:qt5, "qt4agg"=>:qt4, "tkagg"=>:tk,
                          "agg"=>:none,"ps"=>:none,"pdf"=>:none,
                          "svg"=>:none,"cairo"=>:none,"gdk"=>:none,
                          "module://gr.matplotlib.backend_gr"=>:gr)

    qt2gui = Dict("pyqt5"=>:qt_pyqt5, "pyqt4"=>:qt_pyqt4, "pyside"=>:qt_pyside)

    rcParams = PyDict(matplotlib."rcParams")
    default = lowercase(get(ENV, "MPLBACKEND",
                            getnone(rcParams, "backend", "none")))
    if haskey(matplotlib2gui,default)
        defaultgui = matplotlib2gui[default]

        # if the user explicitly requested a particular GUI,
        # it makes sense to ensure that the relevant Conda
        # package is installed (if we are using Conda).
        if conda
            if defaultgui == :qt || defaultgui == :qt4
                # default to pyqt rather than pyside, as below:
                defaultgui = qt2gui[lowercase(getnone(rcParams,"backend.qt4", "pyqt4"))]
                if defaultgui == :qt_pyside
                    pyimport_conda("PySide", "pyside")
                else
                    try
                        pyimport_conda("PyQt5", "pyqt")
                    catch
                        pyimport("PyQt4")
                    end
                end
            elseif defaultgui == :qt5
                pyimport_conda("PyQt5", "pyqt")
            elseif defaultgui == :wx
                pyimport_conda("wx", "wxpython")
            end
        end

        insert!(options, 1, (defaultgui,default))
    end

    try
        # We will get an exception when we import pyplot below (on
        # Unix) if an X server is not available, even though
        # pygui_works and matplotlib.use(backend) succeed, at
        # which point it will be too late to switch backends.  So,
        # throw exception (drop to catch block below) if DISPLAY
        # is not set.  [Might be more reliable to test
        # success(`xdpyinfo`), but only if xdpyinfo is installed.]
        if options[1][1] != :none && Sys.isunix() && !Sys.isapple()
            ENV["DISPLAY"]
        end

        if PyCall.gui == :default
            # try to ensure that GUI both exists and has a matplotlib backend
            for (g,b) in options
                if Sys.isapple() && g == :tk && default != "tkagg" && conda
                    # tkagg backend crashes badly on MacOS
                    g, b = try
                        @info "Installing pyqt package to avoid buggy tkagg backend."
                        pyimport_conda("PyQt5", "pyqt")
                        :qt5, gui2matplotlib[:qt5]
                    catch
                        :none, "agg"
                    end
                end
                if g == :none # Matplotlib is configured to be non-interactive
                    pygui(:default)
                    matplotlib."use"(b)
                    matplotlib."interactive"(false)
                    return (b, g)
                elseif g == :gr
                    return (b, g)
                elseif PyCall.pygui_works(g)
                    # must call matplotlib.use *before* loading backends module
                    matplotlib."use"(b)
                    if g == :qt || g == :qt4
                        if haskey(rcParams,"backend.qt4")
                            g = qt2gui[lowercase(rcParams."backend.qt4")]
                        elseif !PyCall.pyexists("PyQt5") && !PyCall.pyexists("PyQt4")
                            # both Matplotlib and PyCall default to PyQt4
                            # if it is available, but we need to tell
                            # Matplotlib to use PySide otherwise.
                            rcParams."backend.qt4" = "PySide"
                        end
                    end
                    if pyexists("matplotlib.backends.backend_" * lowercase(b))
                        isjulia_display[1] || pygui_start(g)
                        matplotlib."interactive"(!isjulia_display[1] && Base.isinteractive())
                        return (b, g)
                    end
                end
            end
            error("no gui found") # go to catch clause below
        else # the user specified a desired backend via pygui(gui)
            gui = pygui()
            matplotlib."use"(gui2matplotlib[gui])
            if (gui==:qt && !PyCall.pyexists("PyQt5") && !PyCall.pyexists("PyQt4")) || gui==:qt_pyside
                rcParams."backend.qt4" = "PySide"
            end
            isjulia_display[1] || pygui_start(gui)
            matplotlib."interactive"(!isjulia_display[1] && Base.isinteractive())
            return (gui2matplotlib[gui], gui)
        end
    catch e
        if !isjulia_display[1]
            @warn("No working GUI backend found for matplotlib")
            isjulia_display[1] = true
        end
        pygui(:default)
        matplotlib."use"("Agg") # GUI not available
        matplotlib."interactive"(false)
        return ("Agg", :none)
    end
end

# initialization -- anything that depends on Python has to go here,
# so that it occurs at runtime (while the rest of PyPlot can be precompiled).
function __init__()
    global isjulia_display = Bool[isdisplayok()]
    copy!(matplotlib, pyimport_conda("matplotlib", "matplotlib"))
    mvers = matplotlib.__version__
    global version = try
        vparse(mvers)
    catch
        v"0.0.0" # fallback
    end

    backend_gui = find_backend(matplotlib)
    # workaround JuliaLang/julia#8925
    global backend = backend_gui[1]
    global gui = backend_gui[2]
    if Sys.isapple() && gui == :tk
        @warn "PyPlot is using tkagg backend, which is known to cause crashes on MacOS (#410); use the MPLBACKEND environment variable to request a different backend."
    end

    copy!(plt, pyimport("matplotlib.pyplot")) # raw Python module

    pytype_mapping(plt."Figure", Figure)

    copy!(Gcf, pyimport("matplotlib._pylab_helpers")."Gcf")
    copy!(orig_gcf, plt."gcf")
    copy!(orig_figure, plt."figure")
    plt."gcf" = gcf
    plt."figure" = figure

    if isdefined(Main, :IJulia) && Main.IJulia.inited
        Main.IJulia.push_preexecute_hook(force_new_fig)
        Main.IJulia.push_postexecute_hook(display_figs)
        Main.IJulia.push_posterror_hook(close_figs)
    end

    if isjulia_display[1] && gui != :gr && backend != "Agg"
        plt."switch_backend"("Agg")
        plt."ioff"()
    end

    init_colormaps()
end

function pygui(b::Bool)
    if !b != isjulia_display[1]
        if backend != "Agg"
            plt."switch_backend"(b ? backend : "Agg")
            if b
                pygui_start(gui) # make sure event loop is started
                Base.isinteractive() && plt."ion"()
            else
                plt."ioff"()
            end
        elseif b
            error("No working GUI backend found for matplotlib.")
        end
        isjulia_display[1] = !b
    end
    return b
end


text:
# dcontract, dot, tdot, otimes, cross
"""
```julia
dcontract(::SecondOrderTensor, ::SecondOrderTensor)
dcontract(::SecondOrderTensor, ::FourthOrderTensor)
dcontract(::FourthOrderTensor, ::SecondOrderTensor)
dcontract(::FourthOrderTensor, ::FourthOrderTensor)
```
Computes the double contraction between two tensors.
The symbol `⊡`, written `\\boxdot`, is overloaded for double contraction.
The reason `:` is not used is because it does not have the same precedence as multiplication.

**Example:**

```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> B = rand(SymmetricTensor{2, 2});

julia> dcontract(A,B)
1.9732018397544984

julia> A ⊡ B
1.9732018397544984
```
"""
@inline dcontract{dim}(S1::Tensor{2, dim}, S2::Tensor{2, dim}) = tovector(S1) ⋅ tovector(S2)

@inline function dcontract{dim}(S1::Tensor{4, dim}, S2::Tensor{4, dim})
    Tensor{4, dim}(tomatrix(S1) * tomatrix(S2))
end

@inline function dcontract{dim}(S1::Tensor{4, dim}, S2::Tensor{2, dim})
    Tensor{2, dim}(tomatrix(S1) * tovector(S2))
end

@inline function dcontract{dim}(S1::Tensor{2, dim}, S2::Tensor{4, dim})
    Tensor{2, dim}(tomatrix(S2)' * tovector(S1))
end

const ⊡ = dcontract

# Specialized methods for symmetric tensors
@generated function dcontract{dim}(S1::SymmetricTensor{2, dim}, S2::SymmetricTensor{2, dim})
    idx2(i,j) = compute_index(SymmetricTensor{2, dim}, i, j)
    ex1, ex2 = Expr[], Expr[]
    for j in 1:dim, i in j:dim
        if i == j
            push!(ex1, :(get_data(S1)[$(idx2(i, j))]))
            push!(ex2, :(get_data(S2)[$(idx2(i, j))]))
        else
            push!(ex1, :(2 * get_data(S1)[$(idx2(i, j))]))
            push!(ex2, :(    get_data(S2)[$(idx2(i, j))]))
        end
    end
    exp = make_muladd_exp(ex1, ex2)
    return quote
        $(Expr(:meta, :inline))
        $exp
    end
end

@generated function dcontract{dim}(S1::Tensor{2, dim}, S2::SymmetricTensor{4, dim})
    idx4(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    idx2(i,j) = compute_index(Tensor{2, dim}, i, j)
    exps = Expr(:tuple)
    for l in 1:dim, k in l:dim
        ex1, ex2 = Expr[], Expr[]
        for j in 1:dim, i in 1:dim
            push!(ex1, :(data2[$(idx2(i, j))]))
            push!(ex2, :(data4[$(idx4(i, j, k, l))]))
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    quote
        $(Expr(:meta, :inline))
        data2 = get_data(S1)
        data4 = get_data(S2)
        @inbounds r = $exps
        SymmetricTensor{2, dim}(r)
    end
end

@generated function dcontract{dim}(S1::SymmetricTensor{4, dim}, S2::Tensor{2, dim})
    idx4(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    idx2(k,l) = compute_index(Tensor{2, dim}, k, l)
    exps = Expr(:tuple)
    for j in 1:dim, i in j:dim
        ex1, ex2 = Expr[], Expr[]
        for l in 1:dim, k in 1:dim
            push!(ex1, :(data4[$(idx4(i, j, k, l))]))
            push!(ex2, :(data2[$(idx2(k, l))]))
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    quote
        $(Expr(:meta, :inline))
        data2 = get_data(S2)
        data4 = get_data(S1)
        @inbounds r = $exps
        SymmetricTensor{2, dim}(r)
    end
end

@generated function dcontract{dim}(S1::SymmetricTensor{2, dim}, S2::SymmetricTensor{4, dim})
    idx4(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    idx2(i,j) = compute_index(SymmetricTensor{2, dim}, i, j)
    exps = Expr(:tuple)
    for l in 1:dim, k in l:dim
        ex1, ex2 = Expr[], Expr[]
        for j in 1:dim, i in j:dim
            if i == j
                push!(ex1, :(data2[$(idx2(i, j))]))
                push!(ex2, :(data4[$(idx4(i, j, k, l))]))
            else
                push!(ex1, :(2 * data2[$(idx2(i, j))]))
                push!(ex2, :(    data4[$(idx4(i, j, k, l))]))
            end
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    quote
        $(Expr(:meta, :inline))
        data2 = get_data(S1)
        data4 = get_data(S2)
        @inbounds r = $exps
        SymmetricTensor{2, dim}(r)
    end
end

@generated function dcontract{dim}(S1::SymmetricTensor{4, dim}, S2::SymmetricTensor{2, dim})
    idx4(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    idx2(k,l) = compute_index(SymmetricTensor{2, dim}, k, l)
    exps = Expr(:tuple)
    for j in 1:dim, i in j:dim
        ex1, ex2 = Expr[], Expr[]
        for l in 1:dim, k in l:dim
            if k == l
                push!(ex1, :(data4[$(idx4(i, j, k, l))]))
                push!(ex2, :(data2[$(idx2(k, l))]))
            else
                push!(ex1, :(2 * data4[$(idx4(i, j, k, l))]))
                push!(ex2, :(    data2[$(idx2(k, l))]))
            end
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    quote
        $(Expr(:meta, :inline))
        data2 = get_data(S2)
        data4 = get_data(S1)
        @inbounds r = $exps
        SymmetricTensor{2, dim}(r)
    end
end

@generated function dcontract{dim}(S1::SymmetricTensor{4, dim}, S2::SymmetricTensor{4, dim})
    idx4(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    exps = Expr(:tuple)
    for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
        ex1, ex2 = Expr[], Expr[]
        for n in 1:dim, m in n:dim
            if m == n
                push!(ex1, :(data1[$(idx4(i, j, m, n))]))
                push!(ex2, :(data2[$(idx4(m, n, k, l))]))
            else
                push!(ex1, :(2 * data1[$(idx4(i, j, m, n))]))
                push!(ex2, :(    data2[$(idx4(m, n, k, l))]))
            end
        end
        push!(exps.args, make_muladd_exp(ex1, ex2))
    end
    quote
        $(Expr(:meta, :inline))
        data2 = get_data(S2)
        data1 = get_data(S1)
        @inbounds r = $exps
        SymmetricTensor{4, dim}(r)
    end
end

# Promotion
@inline dcontract{dim}(S1::Tensor{2, dim}, S2::SymmetricTensor{2, dim}) = dcontract(promote(S1, S2)...)
@inline dcontract{dim}(S1::Tensor{4, dim}, S2::SymmetricTensor{2, dim}) = dcontract(S1, convert(Tensor, S2))

@inline dcontract{dim}(S1::SymmetricTensor{2, dim}, S2::Tensor{2, dim}) = dcontract(promote(S1, S2)...)
@inline dcontract{dim}(S1::SymmetricTensor{2, dim}, S2::Tensor{4, dim}) = dcontract(convert(Tensor, S1), S2)

@inline dcontract{dim}(S1::Tensor{4, dim}, S2::SymmetricTensor{4, dim}) = dcontract(promote(S1, S2)...)
@inline dcontract{dim}(S1::SymmetricTensor{4, dim}, S2::Tensor{4, dim}) = dcontract(promote(S1, S2)...)

"""
```julia
otimes(::Vec, ::Vec)
otimes(::SecondOrderTensor, ::SecondOrderTensor)
```
Computes the open product between two tensors.
The symbol `⊗`, written `\\otimes`, is overloaded for tensor products.

**Example:**

```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> B = rand(SymmetricTensor{2, 2});

julia> A ⊗ B
2×2×2×2 ContMechTensors.SymmetricTensor{4,2,Float64,9}:
[:, :, 1, 1] =
 0.271839  0.352792
 0.352792  0.260518

[:, :, 2, 1] =
 0.469146  0.608857
 0.608857  0.449607

[:, :, 1, 2] =
 0.469146  0.608857
 0.608857  0.449607

[:, :, 2, 2] =
 0.504668  0.654957
 0.654957  0.48365
```
"""
@inline function otimes{dim}(S1::Tensor{2, dim}, S2::Tensor{2, dim})
    Tensor{4, dim}(tovector(S1) * tovector(S2)')
end

@inline function otimes{dim}(v1::Vec{dim}, v2::Vec{dim})
    Tensor{2, dim}(tovector(v1) * tovector(v2)')
end

const ⊗ = otimes

# Specialized methods for symmetric tensors
@inline function otimes{dim}(S1::SymmetricTensor{2, dim}, S2::SymmetricTensor{2, dim})
    SymmetricTensor{4, dim}(tovector(S1) * tovector(S2)')
end

# Promotion
@inline otimes{dim}(S1::SymmetricTensor{2, dim}, S2::Tensor{2, dim}) = otimes(promote(S1, S2)...)
@inline otimes{dim}(S1::Tensor{2, dim}, S2::SymmetricTensor{2, dim}) = otimes(promote(S1, S2)...)

"""
```julia
dot(::Vec, ::Vec)
dot(::Vec, ::SecondOrderTensor)
dot(::SecondOrderTensor, ::Vec)
dot(::SecondOrderTensor, ::SecondOrderTensor)
```
Computes the dot product (single contraction) between two tensors.
The symbol `⋅`, written `\\cdot`, is overloaded for single contraction.

**Example:**

```jldoctest
julia> A = rand(Tensor{2, 2})
2×2 ContMechTensors.Tensor{2,2,Float64,4}:
 0.590845  0.566237
 0.766797  0.460085

julia> B = rand(Tensor{1, 2})
2-element ContMechTensors.Tensor{1,2,Float64,2}:
 0.794026
 0.854147

julia> dot(A, B)
2-element ContMechTensors.Tensor{1,2,Float64,2}:
 0.952796
 1.00184

julia> A ⋅ B
2-element ContMechTensors.Tensor{1,2,Float64,2}:
 0.952796
 1.00184
```
"""
@inline Base.dot{dim}(v1::Vec{dim}, v2::Vec{dim}) = tovector(v1) ⋅ tovector(v2)

@inline function Base.dot{dim}(S1::Tensor{2, dim}, v2::Vec{dim})
    return Vec{dim}(tomatrix(S1) * tovector(v2))
end

@inline function Base.dot{dim}(v1::Vec{dim}, S2::Tensor{2, dim})
    return Vec{dim}(tomatrix(S2)' * tovector(v1))
end

@inline function Base.dot{dim}(S1::Tensor{2, dim}, S2::Tensor{2, dim})
    return Tensor{2, dim}(tomatrix(S1) * tomatrix(S2))
end

@inline function Base.dot{dim}(S1::SymmetricTensor{2, dim}, S2::SymmetricTensor{2, dim})
    S1_t = convert(Tensor{2, dim}, S1)
    S2_t = convert(Tensor{2, dim}, S2)
    return Tensor{2, dim}(tomatrix(S1_t) * tomatrix(S2_t))
end

@inline Base.dot{dim}(S1::SymmetricTensor{2, dim}, v2::Vec{dim}) = dot(convert(Tensor{2, dim}, S1), v2)

@inline Base.dot{dim}(v2::Vec{dim}, S1::SymmetricTensor{2, dim}) = dot(S1, v2)

# Promotion
Base.dot{dim}(S1::Tensor{2, dim}, S2::SymmetricTensor{2, dim}) = dot(promote(S1, S2)...)
Base.dot{dim}(S1::SymmetricTensor{2, dim}, S2::Tensor{2, dim}) = dot(promote(S1, S2)...)

"""
```julia
tdot(::SecondOrderTensor)
```
Computes the transpose-dot of a second order tensor with itself.
Returns a `SymmetricTensor`.

**Example:**

```jldoctest
julia> A = rand(Tensor{2,3})
3×3 ContMechTensors.Tensor{2,3,Float64,9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> tdot(A)
3×3 ContMechTensors.SymmetricTensor{2,3,Float64,6}:
 1.2577   1.36435   0.48726
 1.36435  1.57172   0.540229
 0.48726  0.540229  0.190334
```
"""
@generated function tdot{dim}(S1::Tensor{2, dim})
    idx(i,j) = compute_index(Tensor{2, dim}, i, j)
    ex = Expr(:tuple)
    for j in 1:dim, i in j:dim
        ex1, ex2 = Expr[], Expr[]
        for k in 1:dim
            push!(ex1, :(get_data(S1)[$(idx(k,i))]))
            push!(ex2, :(get_data(S1)[$(idx(k,j))]))
        end
        push!(ex.args, make_muladd_exp(ex1, ex2))
    end
    return quote
        $(Expr(:meta, :inline))
        @inbounds r = $ex
        SymmetricTensor{2, dim}(r)
    end
end

@inline tdot{dim}(S1::SymmetricTensor{2,dim}) = tdot(convert(Tensor{2,dim}, S1))

"""
```julia
dott(::SecondOrderTensor)
```
Computes the dot-transpose of a second order tensor with itself.
Returns a `SymmetricTensor`.

**Example:**

```jldoctest
julia> A = rand(Tensor{2,3})
3×3 ContMechTensors.Tensor{2,3,Float64,9}:
 0.590845  0.460085  0.200586
 0.766797  0.794026  0.298614
 0.566237  0.854147  0.246837

julia> dott(A)
3×3 ContMechTensors.SymmetricTensor{2,3,Float64,6}:
 0.601011  0.878275  0.777051
 0.878275  1.30763   1.18611
 0.777051  1.18611   1.11112
```
"""
@inline dott(S::SecondOrderTensor) = tdot(transpose(S))

"""
```julia
cross(::Vec, ::Vec)
```
Computes the cross product between two `Vec` vectors, returns a `Vec{3}`. For dimensions 1 and 2 the `Vec`'s
are expanded to 3D first. The infix operator `×` (written `\\times`) can also be used.

**Example:**

```jldoctest
julia> a = rand(Vec{3})
3-element ContMechTensors.Tensor{1,3,Float64,3}:
 0.590845
 0.766797
 0.566237

julia> b = rand(Vec{3})
3-element ContMechTensors.Tensor{1,3,Float64,3}:
 0.460085
 0.794026
 0.854147

julia> a × b
3-element ContMechTensors.Tensor{1,3,Float64,3}:
  0.20535
 -0.24415
  0.116354
```
"""
function Base.cross{T}(u::Vec{3, T}, v::Vec{3, T})
    @inbounds w = Vec{3}((u[2]*v[3] - u[3]*v[2], u[3]*v[1] - u[1]*v[3], u[1]*v[2] - u[2]*v[1]))
    return w
end
function Base.cross{T}(u::Vec{2, T}, v::Vec{2, T})
    @inbounds w = Vec{3}((zero(T), zero(T), u[1]*v[2] - u[2]*v[1]))
    return w
end
function Base.cross{T}(::Vec{1, T}, ::Vec{1, T})
    return zero(Vec{3,T})
end


text:
Base.⋅(d,z) = d[1]*z[1]+d[2]*z[2]

function g3neumann(x,y)
  z = k*abs(y-x)
  if z < 1/16
      C,z2 = 2log(2)-2γ,z*z
      z4=z2*z2
      z6=z2*z4
      ret = k^2*complex( ( (C+2)/8 - (C+3)*z2/64 + (3C+11)*z4/4608 - (6C+25)*z6/442386 -(1/4-z2/32+z4/768-z6/36864)*log(k) )/π, 1/8-z2/64+3z4/4608-6z6/442368)
  else
      ret = im*k/4*hankelh1(1,k*abs(y-x))./abs(y-x) - g1(x,y)./abs(y-x).^2/π -g2(x,y).*logabs(y-x)/π
  end
  ret
end
@vectorize_2arg Number g3neumann

function makegif(x,y,u,L;plotfunction=Main.Plots.PyPlot.contourf,seconds=1,cmap="seismic",vert=1)
    tm=string(time_ns())
    dr = pwd()*"/"*tm*"mov"
    mkdir(dr)

    umax = maxabs(u)
    fps = 24
    MLen = seconds*fps
    for k=1:MLen
        t = 2π/ω*(k-1)/fps
        Main.Plots.PyPlot.clf()
        Main.Plots.PyPlot.axes(aspect="equal")
        plot!(dom;color=:black,legend=false)
        plotfunction(x,y,real(u*exp(-im*ω*t)),L;vmin=-umax*vert,vmax=umax*vert,cmap="seismic")
        xlabel!("\$x\$");ylabel!("\$y\$")
        Main.Plots.PyPlot.savefig(dr * "/" * lpad(k,max(4,ceil(Int,log10(MLen))),0) * ".png";dpi=150,bbox_inches="tight")
    end
    # Requires: brew install imagemagick
    run(`convert -delay 6 -loop 0 $dr/*.png $dr/scattering.gif`)
    run(`open $dr/scattering.gif`)
end

#=
# This is for the landing Helmholtz gif. Solve time on the order of ~70 seconds.
N = 10
r = 5e-2
cr = exp(im*2π*[-0.5:N-1.5]/N)
crl = (1-2im*r)cr
crr = (1+2im*r)cr
dom = ∪(Interval,crl[1:2:end],crr[1:2:end]) ∪ ∪(Circle,cr[2:2:end],ones(length(cr[2:2:end]))r)∪Circle(0.,0.5)
function ui(x,y)
    c = 2exp(im*2π*(0.)/N)
    val = hankelh1(0,k*abs(complex(x,y)-c))
    for i=2:N
        c = 2exp(im*2π*(i-1.)/N)
        val = val + hankelh1(0,k*abs(complex(x,y)-c))
    end
    im/4*val
end
=#

#=
# This is the random screen in Figure 6.3 of Slevinsky & Olver 2015. Solve time on the order of ~70 seconds.
ccr = [-3.0,-2.4710248798864565,-1.7779535080542614,-0.999257770563108,-0.9160576190726175,-0.5056650643725802,0.7258681480228484,1.2291671942613505,1.3417993440008456,1.485081132919861,1.7601585357456848,2.9542404467603642,3.0]
dom = ∪(Interval,(ccr+(3-ccr[end-1])/2)[1:2:end-1],(ccr+(3-ccr[end-1])/2)[2:2:end])
=#