text:
# Enhanced RVQ, aka Stacked Quantizers: https://arxiv.org/abs/1411.2173
export train_ervq, quantize_ervq, experiment_ervq

"""
    quantize_ervq(X, C, V=false) -> B, singletons

Given data and full-dimensional codebooks, quantize.
This methods is identical to [`quantize_rvq`](@ref)

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook.
- `V::Bool`: Whether to print progress

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` codes that approximate `X`
- `singletons::Vector{Matrix{T}}`: `m` matrices with unused codebook entries
"""
function quantize_ervq(
  X::Matrix{T},         # d-by-n. Data to encode
  C::Vector{Matrix{T}}, # codebooks
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # The quantization method is the same as RVQ
  quantize_rvq(X, C, V)
end

"""
    train_ervq(X, B, C, m, h, niter, V=false) -> C, B, error

Train an enhanced residual quantizer / stacked quantizer.
This method is typically initialized by [Residual vector quantization (RVQ)](@ref)

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `B::Matrix{T2}`: `m`-by-`n` matrix with pre-trained codes
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrices. Each matrix is a pretrained codebook of size approximately `d`-by-`h`.
- `m::Integer`: Number of codebooks
- `h::Integer`: Number of entries in each codebook (typically 256)
- `niter::Integer`: Number of iterations to use
- `V::Bool`: Whether to print progress

`T <: AbstractFloat` and `T2 <: Integer`

# Returns
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook of size approximately `d`-by-`h`.
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `error::T`: The quantization error after training
"""
function train_ervq(
  X::Matrix{T},  # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of iterations for training
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  B = convert(Matrix{Int}, B)

  # Then we do the fine-tuning part of https://arxiv.org/abs/1411.2173
  # TODO make sure that C is full-dimensional
  error = qerror(X, B, C)
  if V print("Error after init is $error \n"); end

  for i = 1:niter
    if V print("=== Iteration $i / $niter ===\n"); start_time = time_ns(); end

    # Dummy singletons
    singletons = Vector{Matrix{T}}(undef, m)

    Xr = copy(X)
    Xd = X .- reconstruct(B[2:end,:], C[2:end])

    for j=1:m
      if V print("Updating codebook $j... "); end

      if j == m
        Xd = Xr .- reconstruct(B[j-1,:], C[j-1])
      elseif j > 1
        Xd = Xr .- reconstruct(vcat(B[j-1,:]', B[j+1:end,:]), [C[j-1],C[j+1:end]...])
      end

      # Update the codebook C[j]
      weights = nothing  # use unweighted version of update_centers!
      to_update = zeros(Bool, h)
      to_update[B[j,:]] .= true # In. Whether a codebook entry needs update
      cweights = zeros(Float64, h)   # Out. Cluster weights. We do not use this.
      Clustering.update_centers!(Xd, weights, B[j,:], to_update, C[j], cweights)

      # Check if some centres are unasigned
      if sum(to_update) < h && j == 1
        # If this happens with the first codebook, we have no precomputed
        # singletons, so we have to compute them ourselves.
        # Empirically this is an edge case but it has happened to me so.

        unused = .!to_update # In. Unused centres are centres that were not updated.
        costs = sum((X .- Xd).^2, 1) # In. The cost of assigning the ith vector
        Clustering.repick_unused_centers(Xd, costs, C[j], unused)

      elseif sum(to_update) < h
        # In other cases we already have the precomputed singletons
        ii = 1
        for idx in findall(.!to_update)
          C[j][:,idx] = singletons[2][:, ii]
          ii = ii + 1
        end
      end
      if V print("done.\n"); end

      # Update the residual
      if j > 1
        Xr .-= reconstruct(B[j-1,:], C[j-1])
      end

      if V print("Updating codes... "); end
      B[j:end,:], singletons = quantize_ervq(Xr, C[j:end], V)
      if V print("done. "); end
      error = qerror(X, B, C)
      print("Qerror is $error.\n")

    end # End loop over codebooks

    if V
      error = qerror(X, B, C)
      time_spent = (time_ns() - start_time) / 1e9
      print("Iteration $i / $niter done in $time_spent seconds. Qerror is $error.\n");
    end

  end

  error = qerror(X, B, C)
  C, B, qerror(X, B, C)
end

function train_ervq(
  X::Matrix{T},  # d-by-n. Data to learn codebooks from
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  V::Bool=false) where T <: AbstractFloat# whether to print progress

  # Initialize with RVQ
  C, B, _ = train_rvq(X, m, h, niter, V)
  train_ervq(X, B, C, m, h, niter, V)
end

function experiment_ervq(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # === ERVQ train ===
  d, _ = size(Xt)
  @time C, B, train_error = Rayuela.train_ervq(Xt, B, C, m, h, niter, V)
  norms_B, norms_C = get_norms_codebook(B, C)

  # === Encode the base set ===
  B_base, _ = Rayuela.quantize_ervq(Xb, C, V)
  base_error = qerror(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # Compute and quantize the database norms
  B_base_norms, _ = quantize_norms( B_base, C, norms_C )
  db_norms        = vec( norms_C[ B_base_norms ] )

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, Matrix{Float32}(1.0I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, B_base, recall
end

function experiment_ervq_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # === ERVQ train ===
  d, _ = size(Xt)
  C, B, train_error = Rayuela.train_ervq(Xt, B, C, m, h, niter, V)
  norms_B, norms_C = get_norms_codebook(B, C)
  db_norms     = vec( norms_C[ norms_B ] )

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B, Xq, C, db_norms, Matrix{Float32}(1.0I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, train_error, recall
end

function experiment_ervq(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  C, B, _ = Rayuela.train_rvq(Xt, m, h, niter, V)
  experiment_ervq(Xt, B, C, Xb, Xq, gt, m, h, niter, knn, V)
end

function experiment_ervq_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  C, B, _ = Rayuela.train_rvq(Xt, m, h, niter, V)
  experiment_ervq_query_base(Xt, B, C, Xq, gt, m, h, niter, knn, V)
end


text:
#=
Example 2:
Fit a linear model first, then add a further component and fit a quadratic model.

This example is very similar to #1, but it shows how to check for
goodness of fit, and how to add a new component to obtain a better fit.
=#
using Random, DataFitting

# "True" parameter values:
p0 = 1.
p1 = 2.
p2 = 3.

# Evaluate the model with "True" parameters
x = 0.:0.1:10;
y = collect(p0 .+ p1 .* x .+ p2 .* x.^2);

# Add some noise to simulate a measurement process
noise = 1.
y .+= noise .* randn(MersenneTwister(0), length(x));

# Wrap empirical measure and uncertainties in a `Measure` object
data = Measures(y, noise)

# Prepare a model with just two `ScalarParam` components, add a model
# domain and an expression involving components:
model = Model()                           # Start with empty Model
addcomp!(model, :p0 => ScalarParam(p0))   # add 1st component
addcomp!(model, :p1 => ScalarParam(p1))   # add 2nd component
add_dom!(model, x)                        # add domain
addexpr!(model, :(p0 .+ p1 .* domain[1])) # add model expression

# Fit the model to empirical data
result = fit!(model, data)

# Plot data and best fit model
using Gnuplot
@gp    model[].domain data.val data.unc "w yerr t 'Data'" :-
@gp :- model[].domain model[].expr1 "w line t 'Model'"

# Print the reduced chi-squared
println("Red. χ^2 = ", result.cost / result.dof)

# Clearly this is not a good fit.  Let's add a 3rd `ScalarParam`
# component to obtain a quadratic model
addcomp!(model, :p2 => ScalarParam(p2))   # add 3rd component
addexpr!(model, :(expr1 .+ p2 .* domain[1].^2)) # add model expression
setflag!(model, :expr1, false)

result = fit!(model, data)
@gp    model[].domain data.val data.unc "w yerr t 'Data'" :-
@gp :- model[].domain model[].expr2 "w line t 'Model'"

# Compare "true" parameter values with best fit ones:
using Printf
@printf "%f ∼ %f  (accuracy: %5.3f%%)" p0 result.p0.par.val 100. * (result.p0.par.val-p0) / p0
@printf "%f ∼ %f  (accuracy: %5.3f%%)" p1 result.p1.par.val 100. * (result.p1.par.val-p1) / p1
@printf "%f ∼ %f  (accuracy: %5.3f%%)" p2 result.p2.par.val 100. * (result.p2.par.val-p2) / p2

model = Model(:p0 => ScalarParam(p0), :p1 => ScalarParam(p1))
add_dom!(model, x)                        # add domain
addexpr!(model, :(p0 .+ p1 .* domain[1])) # add model expression
addcomp!(model, :p2 => ScalarParam(p2))   # add 3rd component
replaceexpr!(model, :expr1, :(p0 .+ p1 .* domain[1] .+ p2 .* domain[1].^2))
result = fit!(model, data)
@gp    model[].domain data.val data.unc "w yerr t 'Data'" :-
@gp :- model[].domain model[].expr1 "w line t 'Model'"

ext: jl


text:
using LinearAlgebra, SparseArrays
using Plots

mkpath("$(@__DIR__)/plots")


text:
function drawline!(data::IBuffer, p1::Vector{Float64}, p2::Vector{Float64},
                   color::Vector{Float64})
    all(isfinite.(p1)) && all(isfinite.(p2)) || return    # refuse to draw NaN and Inf

    x = p1[1]
    y = p1[2]
    z = p1[3]
    dx = p2[1] - p1[1]
    dy = p2[2] - p1[2]
    dz = p2[3] - p1[3]

    if abs(dy) < abs(dx)
        if p1[1] > p2[1]
            drawline!(data, p2, p1, color)
        else    # drawline: low
            mx = p2[1]
            yi = 1
            if dy < 0
                yi = -1
                dy = -dy
            end

            d = 2dy - dx
            dist = abs(mx - x)
            dz = dist > 0 ? dz / dist : zero(dz)
            for i in x:mx
                plot!(data, RInt[z, y, i], color)
                d += 2dy
                z += dz

                if d > 0
                    y += yi
                    d -= 2dx
                end
            end
        end
    else
        if p1[2] > p2[2]
            drawline!(data, p2, p1, color)
        else    # drawline: high
            my = p2[2]
            xi = 1
            if dx < 0
                xi = -1
                dx = -dx
            end

            d = 2dx -dy
            dist = abs(my - y)
            dz = dist > 0 ? dz / dist : zero(dz)

            for i in y:my
                plot!(data, RInt[z, i, x], color)
                d += 2dx
                z += dz

                if d > 0
                    x += xi
                    d -= 2dy
                end
            end
        end
    end
end

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
module matrix

    using LinearAlgebra

    # simple function for building a
    # strickly increasing matrix,
    # which is useful for testing.
    function get_inc_mat(T, n::Integer)::Matrix
        mat = Matrix{T}(undef, n, n)

        for row in 1:n
            for col in 1:n
                mat[row, col] = (row - 1)*n + col
            end
        end

        return mat
    end

    #= Rotate Matrix.

    Given an NxN image represented by an NxN matrix,
    where each pixel in the image is 4 bytes, write a method
    to rotate the image by 90°.

    Can you do it in place?
    =#

    export rotate_matrix1

    # long way, basically a brute force solution
    # with O(n²) run-time and O(N) memory.
    # Which, is bad.
    function rotate_matrix1(mat::Array)::Array

        length(mat) <= 1 && return mat

        width, height = size(mat)
        result = Matrix(undef, width, height)

        for row in 1:height
            for col in 1:width
                result[col, height - row + 1] = mat[row, col]
            end
        end

        return result
    end

    export rotate_matrix2

    # Rotate 90+ = transpose, rev
    # in place sln
    function rotate_matrix2(mat::Array)::Array
        length(mat) <= 1 && return mat

        return reverse(transpose(Matrix(mat)), dims=2)
    end

    #= Zero Matrix.

    Write an algorithm such that if an element in an NxM matrix is
    0, its entire row and column are set to 0.
    =#

    export zero_matrix1

    # O(N) solution
    function zero_matrix1(mat::Matrix)::Array

        function zero_row(mat::Matrix, row::Integer)

            height, width = size(mat)

            for index in 1:width
                mat[row, index] = 0
            end
        end

        function zero_column(mat::Matrix, col::Integer)

            height, width = size(mat)

            for index in 1:height
                mat[index, col] = 0
            end
        end

        length(mat) <= 1 && return mat

        height, width = size(mat)

        zero_rows = Set{Integer}()
        zero_cols = Set{Integer}()

        for row in 1:height
            for col in 1:width
                if mat[row, col] == 0
                    zero_rows = union(zero_rows, row)
                    zero_cols = union(zero_cols, col)
                end
            end
        end

        for row in zero_rows
            zero_row(mat, row)
        end

        for col in zero_cols
            zero_column(mat, col)
        end

        return mat
    end

end # end module

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
# types
abstract AbstractCell

immutable PreCell <: AbstractCell
  data::AbstractArray
  width::Int
  height::Int
end

immutable Vertical <: AbstractCell
  data::String
  width::Int
  height::Int
  Vertical(height) = new("|", 1, height)
end

immutable Dash <: AbstractCell
  data::String
  repeat::Int
  Dash(dash::String, n::Int) = new(dash, n)
end

immutable Connector <: AbstractCell
  data::String
  Connector() = new("+")
end

type Cell <: AbstractCell
  data::AbstractArray
  width::Int
  height::Int
end

type Margin
  leftside::Int
  rightside::Int
end

type Mill
  board::AbstractArray
  option::Dict
  Mill(board, option::Dict) = new(board, option)
end

typealias Linear{T<:Union{AbstractCell}} AbstractVector{T}
typealias Horizontal{T<:Union{Dash,Connector}} AbstractVector{T}
typealias PlateVector{T<:Union{Linear,Horizontal}} AbstractVector{T}


text:
module ApproxFunIntrospect
    using ApproxFun, TikzGraphs, LightGraphs

    importall ApproxFun

    # ApproxFun general import
    import ApproxFun: BandedMatrix,order,
                      linesum,complexlength,
                      real, eps, isapproxinteger, ∞, FiniteRange

    # Operator import
    import ApproxFun:    bandinds,SpaceOperator, ConversionWrapper, DerivativeWrapper,
                      rangespace, domainspace,
                      promotedomainspace,  CalculusOperator, interlace, Multiplication,
                       choosedomainspace, IntegralWrapper, LeftIntegralWrapper, RightIntegralWrapper,
                        Dirichlet, Neumann, Laplacian, ConstantTimesOperator, Conversion,
                        dirichlet, neumann, Derivative, ConcreteMultiplication, ConcreteConversion, ConcreteLaplacian,
                        ConcreteDerivative, TimesOperator, MultiplicationWrapper, TridiagonalOperator,
                        InterlaceOperator, ConcreteIntegral, ConcreteLeftIntegral, ZeroOperator, PermutationOperator,
                        ConcreteRightIntegral, LeftDerivativeWrapper, RightDerivativeWrapper

    # Spaces import
    import ApproxFun: PolynomialSpace, ConstantSpace, NoSpace,
                        SumSpace,PiecewiseSpace, ArraySpace,RealBasis,ComplexBasis,AnyBasis,
                        UnsetSpace, canonicalspace, domain, evaluate,
                        AnyDomain, plan_transform,plan_itransform,
                        transform,itransform,transform!,itransform!,
                        isambiguous, fromcanonical, tocanonical, checkpoints, ∂, spacescompatible,
                       mappoint, UnivariateSpace, setdomain, Space, points, space, conversion_rule, maxspace_rule,
                       coefficients, RealUnivariateSpace

    # Multivariate import
    import ApproxFun: BivariateDomain,DirectSumSpace,TupleSpace, AbstractProductSpace,
                        BivariateFun,  ProductFun, LowRankFun, lap, columnspace,
                        blockbandinds, subblockbandinds, fromtensor, totensor, isbandedblockbanded,
                        Tensorizer, tensorizer, block, blockstart, blockstop, blocklengths,
                        domaintensorizer, rangetensorizer, blockrange

    # Jacobi import
    import ApproxFun: jacobip, JacobiSD

    # Singularities
    import ApproxFun: WeightSpace, weight

## General

export introspect

function treeadd_edges!(str,ops,node,M,labels)
    labels[node]=str

    nd=node+1
    for op in ops
        LightGraphs.add_edge!(M,node,nd)
        add_edges!(op,nd,M,labels)
        nd+=treecount(op)
    end
end

function introspect(A::Union{Operator,Space})
    m=treecount(A)

    M=LightGraphs.DiGraph(m)
    labels=Array(String,m)

    add_edges!(A,1,M,labels)

    TikzGraphs.plot(M,labels)
end

## Spaces

texname(FS)=string(typeof(FS))

texname(::Line)="\$(-\\infty,\\infty)\$"
texname(d::Interval)="\$["*string(d.a)*","*string(d.b)*"]\$"

texname(C::ConstantSpace) = "Const."
texname(CS::CosSpace)="Cos"
texname(CS::SinSpace)="Sin"
texname(CS::Chebyshev)="\$T\$"  #*string(first(domain(CS)))*","*string(last(domain(CS)))*"]"
texname(CS::Ultraspherical)="\$U\^{"*string(order(CS))*"}\$"
texname(CS::Taylor)="\$H\^+\$"
texname(CS::Hardy{false})="\$H\^-\$"
texname(::SumSpace)="\$\\oplus\$"
texname(::PiecewiseSpace)="\$\\bigcup\$"
texname{T}(A::ArraySpace{T,1})="["*string(length(A))*"]"
texname{T}(A::ArraySpace{T,2})="["*string(size(A,1))*"\$\\times\$"*string(size(A,2))*"]"
texname(J::JacobiWeight)="\$(1+x)\^{"*string(J.α)*"}(1-x)\^{"*string(J.β)*"}\$"
texname(J::Jacobi) = "\$P^{($(J.a),$(J.b))}\$"

treecount(S::Union{DirectSumSpace,PiecewiseSpace})=1+mapreduce(treecount,+,S.spaces)
treecount(S::Union{ArraySpace,JacobiWeight})=1+treecount(S.space)
treecount(::Space)=1

add_edges!(FS::Space,nd,M,labels)=(labels[nd]=string(nd)*":"*texname(FS))

for (OP) in (:DirectSumSpace,:PiecewiseSpace)
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),A.spaces,nd,M,labels)
end

for (OP) in (:ArraySpace,:JacobiWeight)
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),[A.space],nd,M,labels)
end

## Operators

treecount(::Operator) = 1

treecount(::Derivative) = 1
treecount(::ConstantOperator) = 1
treecount(::Conversion) = 1
treecount(::ConcreteMultiplication) = 1

treecount(M::Union{MultiplicationWrapper,DerivativeWrapper,IntegralWrapper,LeftIntegralWrapper,RightIntegralWrapper,
                   ConversionWrapper,SpaceOperator,ConstantTimesOperator,
                   LeftDerivativeWrapper,RightDerivativeWrapper})=1+treecount(M.op)
treecount(A::Union{PlusOperator,TimesOperator,InterlaceOperator,KroneckerOperator})=1+mapreduce(treecount,+,A.ops)

domainrangestr(A)=texname(domainspace(A))*"\$\\rightarrow\$"*texname(rangespace(A))

texname(A::Operator)=string(typeof(A))
texname(D::ConcreteDerivative)=(D.order==1?"\$D":"\$D\^"*string(D.order))*"\$:"*domainrangestr(D)
texname(A::ConstantOperator)=string(A.λ)*"I"
texname(C::ConcreteConversion)="C:"*domainrangestr(C)
texname(A::ConcreteMultiplication)=space(A.f)==domainspace(A)==rangespace(A)?"M["*texname(space(A.f))*"]":"M["*texname(space(A.f))*"]:"domainrangestr(A)

texname(D::ConcreteIntegral)=(D.order==1?"\$Q":"\$Q\^"*string(D.order))*"\$:"*domainrangestr(D)
texname(D::ConcreteLeftIntegral)=(D.order==1?"\$Q_{$(first(domain(D)))}":"\$Q_{$(first(domain(D)))}\^\{"*string(D.order))*"\}\$:"*domainrangestr(D)

texname(D::DerivativeWrapper)=(D.order==1?"\$(D":"\$(D\^"*string(D.order))*")\$"
texname(D::LeftDerivativeWrapper)=(D.order==1?"\$(D_{$(first(domain(D)))}":"\$(D_{$(first(domain(D)))}\^\{"*string(D.order))*"\})\$"
texname(D::RightDerivativeWrapper)=(D.order==1?"\$(D_{$(last(domain(D)))}":"\$(D_{$(last(domain(D)))}\^\{"*string(D.order))*"\})\$"
texname(Q::LeftIntegralWrapper)=(Q.order==1?"\$(Q_{$(first(domain(Q)))}":"\$(Q_{$(first(domain(Q)))}\^\{"*string(Q.order))*"\})\$"
texname(Q::RightIntegralWrapper)=(Q.order==1?"\$(Q_{$(last(domain(Q)))}":"\$(Q_{$(last(domain(Q)))}\^\{"*string(Q.order))*"\})\$"

texname(A::SpaceOperator)="("*domainrangestr(A)*")"
texname(A::ConstantTimesOperator)=string(A.c)
texname(A::ZeroOperator) = "0"
texname(P::PermutationOperator) = "Perm"

add_edges!(A::Operator,nd,M,labels)=(labels[nd]=string(nd)*":"*texname(A))

for (WRAP,STR) in ((:MultiplicationWrapper,:"(M)"),(:ConversionWrapper,"(C)"))
    @eval add_edges!(A::$WRAP,nd,M,labels)=treeadd_edges!(string(nd)*":"*$STR,[A.op],nd,M,labels)
end

add_edges!(A::Union{SpaceOperator,DerivativeWrapper,
                    LeftDerivativeWrapper,RightDerivativeWrapper,
                    LeftIntegralWrapper,RightIntegralWrapper},nd,M,labels) =
    treeadd_edges!(string(nd)*":"*texname(A),[A.op],nd,M,labels)

for (OP,STR) in ((:PlusOperator,:"+"),(:TimesOperator,:"*"),(:InterlaceOperator,:"Interlace"),
                 (:KroneckerOperator,:"\$\\otimes\$"))
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*$STR,A.ops,nd,M,labels)
end

add_edges!(A::ConstantTimesOperator,nd,M,labels)=treeadd_edges!(string(nd)*":"*string(A.c),[A.op],nd,M,labels)

end #Module


text:
using compathelper_integration_test_repo
using Test

@testset "compathelper_integration_test_repo.jl" begin
    @test 1 == 1
end


text:
function resize_filllast!(v::Vector, n::Integer)
    vlen = length(v)
    (n <= vlen) && return resize!(v, n)
    resize!(v, n)
    (@view v[(vlen+1):end]) .= Ref(v[vlen])
    v
end

function simulate_all!(s0::Simulation; reps=30, max_idle=500)::NamedTuple{(:stepstats, :endstats),Tuple{DataFrame,DataFrame}}
    statefile = tempname()
    println(statefile)
    serialize(statefile, s0)
    res_infected_count = Vector{Vector{Int}}(undef, reps)
    res_m = Vector{Vector{Float64}}(undef, reps)
    worker_ids = Vector{Int}(undef, reps)
    #ss = [deepcopy(s0) for i in 1:Threads.nthreads()]
    #Threads.@threads for rep in 1:reps
    rrs = @distributed (append!) for rep in 1:reps
        local s = deserialize(statefile)
        Random.seed!(rep);
        init!(s)
        simulate_steps!(s; max_idle = max_idle)
        r_inf_count = copy(s.step_infected_agents_count)
        r_m = copy(s.step_total_avg_m_driven)
        [(rep=rep, r_inf_count=r_inf_count, r_m=r_m, worker_id=myid())]
    end
    for rr in rrs
        res_infected_count[rr.rep] = rr.r_inf_count
        res_m[rr.rep] = rr.r_m
        worker_ids[rr.rep] = rr.worker_id
    end
    dfend = DataFrame(rep_no=Int[],steps=Int[],total_m_driven=Float64[],
        median_step=Int[],
        left_tail_sum=Float64[],
        right_tail_sum=Float64[],
        n_agents=Int[], discretize_m=Float64[] ,jump=Bool[],
        g_nodes=Int[], g_edges=Int[],worker_id=Int[])
    half_agents = ceil(Int, s0.p.n_agents/2)
    for rep in 1:reps
        med = findfirst(>=(half_agents), res_infected_count[rep])
        steps_no = length(res_infected_count[rep])
        push!(dfend,(rep, steps_no, res_m[rep][end],
            med,
            sum(res_infected_count[rep][1:med])/s0.p.n_agents,
            steps_no-med+1-sum(res_infected_count[rep][med:steps_no])/s0.p.n_agents,
            s0.p.n_agents, s0.p.discretize_m, s0.p.jump_infect,
            nv(s0.g), ne(s0.g), worker_ids[rep] ))
    end

    maxlen = maximum(length.(res_infected_count))
    resize_filllast!.(res_infected_count,maxlen)
    resize_filllast!.(res_m,maxlen)
    df = DataFrame(step=Int[],total_avg_m_driven=Float64[],
        reps=Int[],infected_avg=Float64[],
        infected_std=Float64[],infected_lo=Float64[],infected_hi=Float64[],
        infected_q05=Float64[],infected_q95=Float64[],
        infected_median=Float64[],
        n_agents=Int[], discretize_m=Float64[] ,jump=Bool[],
        g_nodes=Int[], g_edges=Int[])

    for step in 1:maxlen
        data_for_step = getindex.(res_infected_count,step)
        data_for_meters = getindex.(res_m,step)
        stats = OneSampleTTest(data_for_step)
        conf_interv = confint(stats)
        meanv = mean(data_for_step)
        push!(df,(step,mean(data_for_meters),reps,meanv,std(data_for_step; mean=meanv),
            max(0,conf_interv[1]), min(conf_interv[2],s0.p.n_agents),
            quantile(data_for_step,0.05), quantile(data_for_step,0.95),
            median(data_for_step),
            s0.p.n_agents, s0.p.discretize_m, s0.p.jump_infect,
            nv(s0.g), ne(s0.g)))
    end
    (stepstats=df, endstats=dfend)
end


text:
using WAV
using MFCC
using Plots

function load_data()
    return wavread("0_george_9.wav")
end

function main()
    y, fs = wavread("0_george_9.wav")
    wavplay(y, fs)
    m = mfcc(y, fs)
    heatmap(1:size(m[1],1), 1:size(m[1],2), m[1])
end

main()

max_stars_repo_name: seantfish/chatbot

max_issues_repo_name: seantfish/chatbot

max_forks_repo_name: seantfish/chatbot


text:
struct VonMisesSampler <: Sampleable{Univariate,Continuous}
    μ::Float64
    κ::Float64
    r::Float64

    function VonMisesSampler(μ::Float64, κ::Float64)
        τ = 1.0 + sqrt(1.0 + 4 * abs2(κ))
        ρ = (τ - sqrt(2.0 * τ)) / (2.0 * κ)
        new(μ, κ, (1.0 + abs2(ρ)) / (2.0 * ρ))
    end
end

# algorithm from
#     DJ Best & NI Fisher (1979). Efficient Simulation of the von Mises
#     Distribution. Journal of the Royal Statistical Society. Series C
#     (Applied Statistics), 28(2), 152-157.
function rand(rng::AbstractRNG, s::VonMisesSampler)
    f = 0.0
    local x::Float64
    if s.κ > 700.0
        x = s.μ + randn(rng) / sqrt(s.κ)
    else
        while true
            t, u = 0.0, 0.0
            while true
                d = abs2(rand(rng) - 0.5)
                e = abs2(rand(rng) - 0.5)
                if d + e <= 0.25
                    t = d / e
                    u = 4 * (d + e)
                    break
                end
            end
            z = (1.0 - t) / (1.0 + t)
            f = (1.0 + s.r * z) / (s.r + z)
            c = s.κ * (s.r - f)
            if c * (2.0 - c) > u || log(c / u) + 1 >= c
                break
            end
        end
        acf = acos(f)
        x = s.μ + (rand(rng, Bool) ? acf : -acf)
    end
    return x
end


text:
using Yfjd
using Test

@testset "Yfjd.jl" begin
    # Write your tests here.
end


text:
"""
    A mixed continuous-discrete distribution
"""
struct MixedDistribution{DT<:Dst.ContinuousUnivariateDistribution} <: Dst.ContinuousUnivariateDistribution
    mass_probs::Vector{Float64}
    mass_points::Vector{Float64}
    cont_weights::Vector{Float64}
    cont_dists::Vector{DT}
    function MixedDistribution(mass_probs::Vector{Float64}, mass_points::Vector{Float64}, cont_weights::Vector{Float64}, cont_dists::Vector{DT}) where {DT <: Dst.ContinuousUnivariateDistribution}
        length(mass_probs) == length(mass_points)  || ArgumentError("Mass length inconsistency")
        length(cont_weights) == length(cont_dists) || ArgumentError("Continuous length inconsistency")
        # sort mass points
        pairs = sort(zip(mass_points, mass_probs) |> collect, by = t -> t[1])
        mass_points = [t[1] for t in pairs]
        mass_probs =  [t[2] for t in pairs]
        s = sum(mass_probs) + sum(cont_weights)
        s ≈ one(eltype(mass_probs)) && return new{DT}(mass_probs, mass_points, cont_weights, cont_dists)
        @warn "Prior probabilities not normalized"
        new{DT}(mass_probs./s, mass_points, cont_weights./s, cont_dists)
    end
end

function Base.maximum(md::MixedDistribution)
    max_mass = isempty(md.mass_points) ? -Inf64 : maximum(md.mass_points)
    max_cont = isempty(md.cont_dists)  ? -Inf64 : maximum.(md.cont_dists) |> maximum
    return max(max_mass, max_cont)
end

function Base.minimum(md::MixedDistribution)
    min_mass = isempty(md.mass_points) ? +Inf64 : minimum(md.mass_points)
    min_cont = isempty(md.cont_dists)  ? +Inf64 : minimum.(md.cont_dists) |> minimum
    return min(min_mass, min_cont)
end

function Dst.ncategories(md::MixedDistribution)
    return length(md.mass_probs)
end

function Statistics.mean(md::MixedDistribution)
	μ = sum(md.mass_points.*md.mass_probs)
	for idx in eachindex(md.cont_dists)
		μ += mean(md.cont_dists[idx]) * md.cont_weights[idx]
	end
	return μ
end

"""
	Variance of mixed continuous discrete computed as a mixture:
	Var[X] = E[Var[X|θ]] + Var[E[X|θ]]
"""
function Statistics.var(md::MixedDistribution)
	disc_weight = sum(md.mass_probs)
    μ_discrete = sum(md.mass_points.*md.mass_probs)/disc_weight
    σ2d = sum((md.mass_points.-μ_discrete).^2 .* md.mass_probs)
    μ = mean(md)
    tot_var = isempty(md.mass_probs) ? 0.0 : disc_weight * (σ2d + μ_discrete^2)
    for idx in eachindex(md.cont_dists)
		tot_var += md.cont_weights[idx] * (var(md.cont_dists[idx]) + mean(md.cont_dists[idx])^2)
    end
    return tot_var - μ^2
end

function Dst.insupport(md::MixedDistribution, x::Real)
    return minimum(md) <= x && x <= maximum(md)
end

function Dst.pdf(md::MixedDistribution, x::Real)
    any(md.mass_points .≈ x) && DomainError("PDF undefined at a mass point")
end

function Dst.cdf(md::MixedDistribution, x::Real)
    cdf_point = 0.0
    for idx in eachindex(md.mass_probs)
        if md.mass_points[idx] <= x
            cdf_point += md.mass_probs[idx]
        end
    end
    for idx in eachindex(md.cont_dists)
        cdf_point += md.cont_weights[idx] * Dst.cdf(md.cont_dists[idx], x)
    end
    return cdf_point
end

"""
	Number of components of the mixed distribution:
	number of discrete points + continuous distributions
"""
Dst.ncomponents(md::MixedDistribution) = Dst.ncategories(md) + length(md.cont_dists)

"""
	Returns a vector of prior probabilities [discrete priors, continuous priors]
"""
Dst.probs(md::MixedDistribution) = vcat(md.mass_probs, md.cont_weights)

Base.rand(rng::AbstractRNG, md::MixedDistribution) = ordered_rand(md, rand(rng))

Base.rand(md::MixedDistribution) = rand(Random.GLOBAL_RNG, md)

"""
Integration of function including Dirac peaks
∫x⋅f(x)dx
Uses the trapezoid method with given step dx
"""
function partial_expectancy(md::MixedDistribution; xlow = minimum(md), xhigh = maximum(md), dx = min((xhigh - xlow) * 0.05, 0.01))
    xhigh > xlow || throw(ArgumentError("xlow = $xlow greater than xhigh = $xhigh"))
    disc_idx = 1
    while disc_idx <= length(md.mass_points) && md.mass_points[disc_idx] < xlow
        disc_idx += 1
    end
    still_mass = disc_idx <= length(md.mass_points) && md.mass_points[disc_idx] <= xhigh
    x = xlow
    integrator = 0.0
    while x <= xhigh
        δ = min(dx, xhigh - x)
        # check if mass point in chunk
        if still_mass && md.mass_points[disc_idx] <= x + δ
            integrator += md.mass_points[disc_idx] * md.mass_probs[disc_idx]
            disc_idx += 1
            if disc_idx > length(md.mass_points) || md.mass_probs[disc_idx] > xhigh
                still_mass = false
            end
        end
        # add continuous chunk
        for (d,p) in zip(md.cont_dists, md.cont_weights)
            if Dst.insupport(d, x) && Dst.insupport(d, x + δ)
                integrator += p * 0.5 * δ * (x + δ*0.5) * (Dst.pdf(d,x) + Dst.pdf(d,x+δ))
            end
        end
        x += dx
    end
    return integrator
end

"""
Computing integral of function including Dirac peaks in intervals between eval_points
"""
function partial_expectancy(md::MixedDistribution, eval_points::V; dx = min((xhigh - xlow) * 0.05, 0.01)) where {V<:AbstractVector{<:Real}}
    intervals = begin
        sorted_eval = sort(eval_points)
        collect(zip(sorted_eval[1:end-1], sorted_eval[2:end]))
    end
    ys = [partial_expectancy(md, xlow = xl, xhigh = xh, dx = dx) for (xl,xh) in intervals]
    (intervals, ys)
end

"""
Exception used when Cumulated Distribution Function does not integrate to 1 in a computation
"""
struct CDFException{D<:Dst.Distribution} <: Exception
	d::D
	CDFException(d::D) where {D<:Dst.Distribution} = new{D}(d)
end

"""
	Ordered_rand returns the quantity corresponding to
	an already-generated random number
"""
function ordered_rand(md::MixedDistribution, r)
	cs_discrete = cumsum(md.mass_probs)
	for idx in eachindex(cs_discrete)
		if r <= cs_discrete[idx]
			return md.mass_points[idx]
		end
	end
	prob_offset = isempty(cs_discrete) ? 0.0 : last(cs_discrete)
	cs_cont = vcat(cs_discrete, cumsum(md.cont_weights) .+ prob_offset)
	for idx in length(cs_discrete)+1:length(cs_cont)
		if r <= cs_cont[idx]
			cont_idx = idx-length(cs_discrete)
			last_prob = idx == 1 ? 0.0 : cs_cont[idx-1]
			scaled_point = (r - last_prob) / (cs_cont[idx] - last_prob)
			return Dst.quantile(md.cont_dists[cont_idx], scaled_point)
		end
	end
	# all points reached beforehand
	throw(CDFException(md))
end


text:
using Revise
using ADCME
using PyPlot
using AdFem
using Test

mesh = Mesh(2, 2, 0.5)
@test mesh.nodes≈[ 0.0  0.0
                0.5  0.0
                1.0  0.0
                0.0  0.5
                0.5  0.5
                1.0  0.5
                0.0  1.0
                0.5  1.0
                1.0  1.0]

@test get_ngauss(mesh)≈24

nodes = gauss_nodes(mesh)
close("all")
visualize_mesh(mesh)
scatter(nodes[:,1], nodes[:,2])
savefig("g.png")

@test get_area(mesh)≈ones(8)*0.25/2

@PDATA
Mesh("twoholes.stl")

@testset "fem_to_gauss_points" begin
    mesh = Mesh(10, 10, 1.0)
    x, y = mesh.nodes[:,1], mesh.nodes[:,2]
    u = @. x^2 + y
    out = fem_to_gauss_points(constant(u), mesh)
    out2 = fem_to_gauss_points(u, mesh)
    sess = Session(); init(sess)
    out1 = run(sess, out)
    @test out1≈out2
    # close("all")
    # visualize_scalar_on_gauss_points(out1, mesh)
    # savefig("test2.png")

    # close("all")
    # visualize_scalar_on_fem_points(u, mesh)
    # savefig("test.png")
end

@testset "visualize_scalar_on_fvm_points" begin
    mesh = Mesh(10, 10, 1.0)
    u = rand(size(mesh.elems, 1))
    visualize_scalar_on_fvm_points(u, mesh)
    # close("all")
    # visualize_scalar_on_gauss_points(out1, mesh)
    # savefig("test2.png")

    # close("all")
    # visualize_scalar_on_fem_points(u, mesh)
    # savefig("test.png")
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

@testset "Universal Generating Operator" begin

    @testset "Series-Parallel Network" begin
        # Example from: Multi-State System Reliability Analysis and Optimization
        # for Engineers and Industrial Managers (p. 11/209)
        val_sol = [0u"m^3/hr", 1500u"m^3/hr", 1800u"m^3/hr", 2000u"m^3/hr", 3500u"m^3/hr"]
        prb_sol = [0.01207138, 0.10333981359, 0.07371573369, 0.05304777100, 0.75782529972]

        global λ¹, μ¹ = 7.0u"1/yr", 100.0u"1/yr"
        global λ², μ² = 10.0u"1/yr", 80.0u"1/yr"
        global λ³ᵃ, λ³ᵇ = 7.0u"1/yr", 10.0u"1/yr"
        global μ³ᵃ, μ³ᵇ = 120.0u"1/yr", 110.0u"1/yr"

        ntwᶠᵗˢ = include(joinpath(_MSS.BASE_DIR,"test/networks/flow_transmission_system.jl"))
        solve!(ntwᶠᵗˢ)
        val_ntw = ntwᶠᵗˢ.usr[1][:ugf].val
        prb_ntw = ntwᶠᵗˢ.usr[1][:ugf].prb

        @test isapprox.(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))
    end

    @testset "Series-Parallel Network with Measurements Uncertainty" begin
        # Example from: Multi-State System Reliability Analysis and Optimization
        # for Engineers and Industrial Managers (p. 11/209)
        val_sol = [0u"m^3/hr", 1500u"m^3/hr", 1800u"m^3/hr", 2000u"m^3/hr", 3500u"m^3/hr"]
        prb_sol = [0.01207138, 0.10333981359, 0.07371573369, 0.05304777100, 0.75782529972]

        global λ¹, μ¹ = (7.0±1.0)u"1/yr", 100.0u"1/yr"
        global λ², μ² = 10.0u"1/yr", (80.0±10.0)u"1/yr"
        global λ³ᵃ, λ³ᵇ = 7.0u"1/yr", (10.0±2.0)u"1/yr"
        global μ³ᵃ, μ³ᵇ = (120.0±5.0)u"1/yr", 110.0u"1/yr"

        ntwᶠᵗˢ = include(joinpath(_MSS.BASE_DIR,"test/networks/flow_transmission_system.jl"))
        solve!(ntwᶠᵗˢ)
        val_ntw = ntwᶠᵗˢ.usr[1][:ugf].val
        prb_ntw = ntwᶠᵗˢ.usr[1][:ugf].prb

        @test isapprox.(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))

        @test isapprox(_MSM.uncertainty(sum(prb_ntw)), 0.0, atol=1e-6)
    end

    @testset "Bridge Network" begin
        val_sol = [0.0u"m^3/hr", 1.0u"m^3/hr", 2.0u"m^3/hr"]
        prb_sol = [0.02152, 0.32238, 0.6561]

        ntwᵇⁿ = include(joinpath(_MSS.BASE_DIR,"test/networks/bridge_network.jl"))
        solve!(ntwᵇⁿ)
        val_ntw = ntwᵇⁿ.usr[1][:ugf].val
        prb_ntw = ntwᵇⁿ.usr[1][:ugf].prb

        @test isapprox.(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))
    end

    @testset "Double Bridge Network" begin
        val_sol = [0.0u"m^3/hr", 1.0u"m^3/hr", 2.0u"m^3/hr"]
        prb_sol = [0.0102772, 0.1827198, 0.807003]

        ntwᵈᵇⁿ = include(joinpath(_MSS.BASE_DIR,"test/networks/double_bridge_network.jl"))
        solve!(ntwᵈᵇⁿ)
        val_ntw = ntwᵈᵇⁿ.usr[1][:ugf].val
        prb_ntw = ntwᵈᵇⁿ.usr[1][:ugf].prb

        @test isapprox.(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))
    end

    @testset "Disconnected Sources and Users" begin
        # Example from: Industrial Energy System Availability Management (p. 92)
        val_sol = [0.0u"MW", 1.0u"MW"]
        prb_sol = [0.1, 0.9]

        ntwᵈˢᶠ = include(joinpath(_MSS.BASE_DIR,"test/networks/double_single_feeder.jl"))
        solve!(ntwᵈˢᶠ)
        val_usr_1, prb_usr_1 = ntwᵈˢᶠ.usr[1][:ugf].val, ntwᵈˢᶠ.usr[1][:ugf].prb
        val_usr_2, prb_usr_2 = ntwᵈˢᶠ.usr[2][:ugf].val, ntwᵈˢᶠ.usr[2][:ugf].prb

        @test isapprox.(sum(prb_usr_1), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_usr_1, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_usr_1, rtol=1e-6))

        @test isapprox.(sum(prb_usr_2), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_usr_2, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_usr_2, rtol=1e-6))
    end

    @testset "Source Dependence" begin
        # Example from: Analytical Model for Availability Assessment of Large-
        # Scale Offshore Wind Farms including their Collector System (p. 5)
        val_sol = [0.0u"MW", 2.0u"MW", 4.0u"MW", 6.0u"MW", 8.0u"MW"]
        prb_sol = [0.307,0.0126,0.11907,0.10206,0.45927]

        ntwʷᶠ = include(joinpath(_MSS.BASE_DIR,"test/networks/wind_farm.jl"))
        solve!(ntwʷᶠ)
        val_ntw = ntwʷᶠ.usr[1][:ugf].val
        prb_ntw = ntwʷᶠ.usr[1][:ugf].prb

        @test isapprox.(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))
    end

    @testset "Evaluation Dependence" begin
        ugfᵉ¹ = UGF(:power, [0.0u"MW",1.0u"MW",0.0u"MW",1.0u"MW"],
                            [0.1,0.2,0.3,0.4], rdc=false)
        ugfᵉ² = UGF(:power, [0.0u"MW",0.0u"MW",1.0u"MW",1.0u"MW"],
                            [0.1,0.2,0.3,0.4], rdc=false)

        # parallel src system: str = ugfᵉ¹.val + ugfᵉ².val
        # 1) val = 0.0u"MW" → prb = 0.1
        # 2) val = 1.0u"MW" → prb = 0.2 + 0.3 = 0.5
        # 2) val = 2.0u"MW" → prb = 0.4
        val_sol = [0.0u"MW", 1.0u"MW", 2.0u"MW"]
        prb_sol = [0.1, 0.5, 0.4]

        ntwˣ = Network()
        add_sources!(ntwˣ, node = [1,2],
                           ugf = [ugfᵉ¹,ugfᵉ²],
                           eval_dep = true)
        add_user!(ntwˣ, node = 3)
        add_components!(ntwˣ, edge = [(1,3),(2,3)])
        solve!(ntwˣ)
        val_ntw = ntwˣ.usr[1][:ugf].val
        prb_ntw = ntwˣ.usr[1][:ugf].prb

        @test isapprox(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))

        # series cmp system: str = min(ugfᵉ¹.val,ugfᵉ².val)
        # 1) val = 0.0u"MW" → prb = 0.1 + 0.2 + 0.3 = 0.6
        # 2) val = 1.0u"MW" → prb = 0.4
        val_sol = [0.0u"MW", 1.0u"MW"]
        prb_sol = [0.6, 0.4]

        ntwˢ = Network()
        add_source!(ntwˢ, node = 1)
        add_user!(ntwˢ, node = 3)
        add_components!(ntwˢ, edge = [(1,2),(2,3)],
                              ugf = [ugfᵉ¹,ugfᵉ²],
                              eval_dep = true)
        solve!(ntwˢ)
        val_ntw = ntwˢ.usr[1][:ugf].val
        prb_ntw = ntwˢ.usr[1][:ugf].prb

        @test isapprox(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))

        # parallel cmp system: str = ugfᵉ¹.val + ugfᵉ².val
        # 1) val = 0.0u"MW" → prb = 0.1
        # 2) val = 1.0u"MW" → prb = 0.2 + 0.3 = 0.5
        # 2) val = 2.0u"MW" → prb = 0.4
        val_sol = [0.0u"MW", 1.0u"MW", 2.0u"MW"]
        prb_sol = [0.1, 0.5, 0.4]

        ntwᵖ = Network()
        add_source!(ntwᵖ, node = 1)
        add_user!(ntwᵖ, node = 2)
        add_components!(ntwᵖ, edge = [(1,2),(1,2)],
                              ugf = [ugfᵉ¹,ugfᵉ²],
                              eval_dep = true)
        solve!(ntwᵖ)
        val_ntw = ntwᵖ.usr[1][:ugf].val
        prb_ntw = ntwᵖ.usr[1][:ugf].prb

        @test isapprox(sum(prb_ntw), 1.0, rtol=1e-6)
        @test all(isapprox.(val_sol, val_ntw, rtol=1e-6))
        @test all(isapprox.(prb_sol, prb_ntw, rtol=1e-6))
    end

end


text:
abstract type AbstractOutput{T} end
"""

Preallocates and holds arrays for positions, velocities, and Jacobian at every integrator step
"""
# Should add option to choose out intervals, checkpointing, etc.
struct CartesianOutput{T<:AbstractFloat} <: AbstractOutput{T}
    states::Vector{State{T}}
    nstep::Int64
    filename::String
    file::Bool

    function CartesianOutput(::Type{T},nbody::Int64,nstep::Int64;filename="data.jld2",file=false) where T <: AbstractFloat
        states = Array{State,1}(undef,nstep)
        return new{T}(states,nstep,filename,file)
    end
end

# Allow user to not have to specify type. Defaults to Float64
CartesianOutput(nbody::T,nstep::T;filename="data.jld2") where T<:Int64 = CartesianOutput(Float64,nbody,nstep,filename=filename)

"""

Runs integrator like (*insert doc reference here*) and output positions, velocities, and Jacobian to a JLD2 file.
"""
function (intr::Integrator)(s::State{T},o::CartesianOutput{T}) where T<:AbstractFloat
    t0 = s.t[1] # Initial time
    time = intr.tmax
    nsteps = o.nstep

    # Integrate in proper direction
    h = intr.h * check_step(t0,time)
    tmax = t0 + (h * nsteps)

    # Preallocate struct of arrays for derivatives (and pair)
    d = Derivatives(T,s.n)

    for i in 1:nsteps
        # Save State from current step
        o.states[i] = deepcopy(s)
        # Take integration step and advance time
        intr.scheme(s,d,h)
        s.t[1] = t0 +  (h * i)
    end

    # Output to jld2 file if desired
    if o.file; save(o.filename,Dict("states" => o.states)); end
    return
end

# Includes
const outputs = ["elements"]
for i in outputs; include("$(i).jl"); end


text:
#Functions to plot matrix structure of optigraphs
rectangle(w, h, x, y) = Plots.Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

"""
    PlasmoPlots.matrix_plot(graph::OptiGraph;node_labels = false,labelsize = 24,subgraph_colors = false,node_colors = false,markersize = 1)

Plot a matrix visualization of the optigraph: `graph`. The following keyword arguments can be provided to customize the matrix visual.

* `node_labels = false`: whether to label nodes using the corresponding optinode label attribute.
* `labelsize`: the size for each node label.  Only active if `node_labels = true`.
* `subgraph_colors = false`: whether to color nodes according to their subgraph.
* `node_colors = false`: whether to color nodes.  Only active if `subgraph_colors = false`.
* `markersize = 1`: Size of the linking constraints in the matrix representation.
"""
function matrix_plot(graph::OptiGraph;node_labels = false,labelsize = 24,subgraph_colors = false,node_colors = false,markersize = 1)

    n_graphs = length(graph.subgraphs)
    if subgraph_colors
        cols = Colors.distinguishable_colors(n_graphs + 1)
        if cols[1] == Colors.colorant"black"
            cols[1] = Colors.colorant"grey"
        end
        colors = cols[2:end]
    else
        colors = [Colors.colorant"grey" for _= 1:n_graphs]
    end

    if node_colors
        cols = Colors.distinguishable_colors(length(all_nodes(graph)) + 1)
        if cols[1] == Colors.parse(Colorant,"black")
            cols[1] = Colors.parse(Colorant,"grey")
        end
        node_cols = cols[2:end]
    end

    #Plot limits
    n_vars_total = num_all_variables(graph)
    n_cons_total = num_all_constraints(graph)
    n_linkcons_total = num_all_linkconstraints(graph)

    n_all_cons_total = n_cons_total + n_linkcons_total #n_link_edges_total

    if n_all_cons_total >= 5
        yticks = Int64.(round.(collect(range(0,stop = n_all_cons_total,length = 5))))
    else
        yticks = Int64.(round.(collect(range(0,stop = n_all_cons_total,length = n_all_cons_total + 1))))
    end

    #Setup plot dimensions
    plt = Plots.plot(;xlims = [0,n_vars_total],ylims = [0,n_all_cons_total],legend = false,framestyle = :box,xlabel = "Node Variables",ylabel = "Constraints",size = (800,800),
    guidefontsize = 24,tickfontsize = 18,grid = false,yticks = yticks)

    #plot top level nodes, then start going down subgraphs
    n_link_constraints = num_link_constraints(graph)  #local links
    col = 0
    node_indices = Dict()
    node_col_ranges = Dict()
    for (i,node) in enumerate(all_nodes(graph))
        node_indices[node] = i
        node_col_ranges[node] = [col,col + num_variables(node)]
        col += num_variables(node)
    end

    row = n_all_cons_total  - n_link_constraints #- height_initial
    #draw node blocks for this graph
    for (i,node) in enumerate(getnodes(graph))
        height = num_constraints(node)
        row -= height
        #row_start,row_end = node_row_ranges[node]
        row_start = row
        col_start,col_end = node_col_ranges[node]
        width = col_end - col_start

        row_end = row - height
        rec = rectangle(width,height,col_start,row_start)

        if !(node_colors)
            Plots.plot!(plt,rec,opacity = 1.0,color = :grey)
        else
            Plots.plot!(plt,rec,opacity = 1.0,color = node_cols[i])
        end
        if node_labels
            Plots.annotate!(plt,(col_start + width + col_start)/2,(row + height + row)/2,Plots.text(node.label,labelsize))
        end
    end

    #plot link constraints for highest level using rectangles
    row = n_all_cons_total
    #recs = []

    link_rows = []
    link_cols = []
    for link in getlinkconstraints(graph)
        #row -= 1

        linkcon = constraint_object(link)
        vars = keys(linkcon.func.terms)
        for var in vars
            node = getnode(var)

            col_start,col_end = node_col_ranges[node]
            col_start = col_start + var.index.value - 1 + 0.5

            #these are just points.
            #rec = rectangle(1,1,col_start,row)
            push!(link_rows,row - 0.5)
            push!(link_cols,col_start)
            # Plots.plot!(plt,rec,opacity = 1.0,color = :blue);
        end
        row -= 1
    end
    Plots.scatter!(plt,link_cols,link_rows,markersize = markersize,markercolor = :blue,markershape = :rect);

    if length(graph.optinodes) > 0
        row -= 1
    end

    _plot_subgraphs!(graph,plt,node_col_ranges,row,node_labels = node_labels,labelsize = labelsize,colors = colors,markersize = markersize)
    return plt
end

function _plot_subgraphs!(graph::OptiGraph,plt,node_col_ranges,row_start_graph;node_labels = false,labelsize = 24,colors = nothing,markersize = 1)
    if colors == nothing
        colors = [Colors.parse(Colorant,"grey") for _= 1:length(graph.subgraphs)]
    end

    row_start_graph = row_start_graph
    for (i,subgraph) in enumerate(getsubgraphs(graph))

        link_rows = []
        link_cols = []
        row = row_start_graph#

        for link in getlinkconstraints(subgraph)
            #row -= 1
            #nodes = getnodes(link)
            linkcon = constraint_object(link)
            vars = keys(linkcon.func.terms)
            for var in vars
                node = getnode(var)
                col_start,col_end = node_col_ranges[node]
                col_start = col_start + var.index.value - 1 + 0.5
                # rec = rectangle(1,1,col_start,row)
                # Plots.plot!(plt,rec,opacity = 1.0,color = :blue)
                push!(link_rows,row - 0.5)
                push!(link_cols,col_start)
            end
            row -= 1
        end
        Plots.scatter!(plt,link_cols,link_rows,markersize = markersize,markercolor = :blue,markershape = :rect);

        if !(isempty(subgraph.optinodes))
            subgraph_col_start = node_col_ranges[subgraph.optinodes[1]][1]
        else
            subgraph_col_start = 0
        end

        #draw node blocks for this graph
        for node in getnodes(subgraph)
            height = num_constraints(node)
            row -= height
            row_start = row
            col_start,col_end = node_col_ranges[node]
            width = col_end - col_start

            rec = rectangle(width,height,col_start,row_start)
            Plots.plot!(plt,rec,opacity = 1.0,color = colors[i])
            if node_labels
                Plots.annotate!(plt,(col_start + width + col_start)/2,(row + height + row)/2,Plots.text(node.label,labelsize))
            end

        end

        _plot_subgraphs!(subgraph,plt,node_col_ranges,row,node_labels = node_labels,labelsize = labelsize)

        num_cons = num_all_constraints(subgraph) + num_all_linkconstraints(subgraph)
        num_vars = num_all_variables(subgraph)
        row_start_graph -= num_cons
        subgraph_row_start = row_start_graph

        rec = rectangle(num_vars,num_cons,subgraph_col_start,subgraph_row_start)
        Plots.plot!(plt,rec,opacity = 0.1,color = colors[i])

    end
end

#Overlap spy
function matrix_plot(graph::OptiGraph,subgraphs::Vector{OptiGraph};node_labels = false,labelsize = 24,subgraph_colors = true)

    n_graphs = length(subgraphs)
    if subgraph_colors
        cols = Colors.distinguishable_colors(n_graphs + 1)
        if cols[1] == Colors.parse(Colorant,"black")
            cols[1] = Colors.parse(Colorant,"grey")
        end
        colors = cols[2:end]
    else
        colors = [Colors.parse(Colorant,"grey") for _= 1:n_graphs]
    end

    #Plot limits
    n_vars_total = sum(num_variables.(subgraphs)) #+ sum(num_variables.(getnodes(graph))) #master
    n_cons_total = sum(num_all_constraints.(subgraphs)) #+ sum(num_constraints.(getnodes(graph))) #+ num_linkconstraints(graph)
    n_linkcons_total = sum(num_all_linkconstraints.(subgraphs)) #+ num_all_linkconstraints(graph)

    n_all_cons_total = n_cons_total + n_linkcons_total

    if n_all_cons_total >= 5
        yticks = Int64.(round.(collect(range(0,stop = n_all_cons_total,length = 5))))
    else
        yticks = Int64.(round.(collect(range(0,stop = n_all_cons_total,length = n_all_cons_total + 1))))
    end

    #Setup plot dimensions
    plt = Plots.plot(;xlims = [0,n_vars_total],ylims = [0,n_all_cons_total],legend = false,framestyle = :box,xlabel = "Node Variables",ylabel = "Constraints",size = (800,800),
    guidefontsize = 24,tickfontsize = 18,grid = false,yticks = yticks)

    row_start_graph = n_all_cons_total - 1
    col_start_graph = 1
    for i = 1:length(subgraphs)
        subgraph = subgraphs[i]
        #column data for subgraph
        node_indices = Dict()
        node_col_ranges = Dict()

        col = col_start_graph
        for (i,node) in enumerate(all_nodes(subgraph))
            node_indices[node] = i
            node_col_ranges[node] = [col,col + num_variables(node)]
            col += num_variables(node)
        end

        #Now just plot columns of overlap nodes
        nodes = all_nodes(subgraph)
        overlap_nodes = Dict()
        for j = 1:length(subgraphs)
            if j != i
                other_subgraph = subgraphs[j]
                other_nodes = all_nodes(other_subgraph)
                overlap = intersect(nodes,other_nodes)
                overlap_nodes[j] = overlap
            end
        end
                #plot local column overlap
        link_rows = []
        link_cols = []
        row = row_start_graph
        for link in getlinkconstraints(subgraph)
            row -= 1
            linkcon = constraint_object(link)
            vars = keys(linkcon.func.terms)
            for var in vars
                node = getnode(var)
                col_start,col_end = node_col_ranges[node]
                col_start = col_start + var.index.value - 1
                push!(link_rows,row)
                push!(link_cols,col_start)
            end
        end
        Plots.scatter!(plt,link_cols,link_rows,markersize = 1,markercolor = :blue,markershape = :rect);

        #draw node blocks for this graph
        for node in getnodes(subgraph)
            height = num_constraints(node)
            row -= height
            row_start = row
            col_start,col_end = node_col_ranges[node]
            width = col_end - col_start

            rec = rectangle(width,height,col_start,row_start)
            Plots.plot!(plt,rec,opacity = 1.0,color = colors[i])
            if node_labels
                Plots.annotate!(plt,(col_start + width + col_start)/2,(row + height + row)/2,Plots.text(node.label,labelsize))
            end
        end

        num_cons = num_all_constraints(subgraph) + num_all_linkconstraints(subgraph)
        num_vars = num_all_variables(subgraph)

        subgraph_plt_start = row
        rec = rectangle(num_vars,num_cons,col_start_graph,subgraph_plt_start)
        Plots.plot!(plt,rec,opacity = 0.1,color = colors[i])

        #overlap rectanges
        for (j,overlap) in overlap_nodes
            for node in overlap
                col_start,col_end = node_col_ranges[node]
                rec = rectangle(num_variables(node),num_cons,col_start,subgraph_plt_start)
                Plots.plot!(plt,rec,opacity = 0.1,color = colors[j])
            end
        end

        col_start_graph = col
        row_start_graph = row
    end

    return plt
end

Plots.spy(graph::OptiGraph;kwargs...) = matrix_plot(graph;kwargs...)
Plots.spy(graph::OptiGraph,subgraphs::Vector{OptiGraph};kwargs...) = matrix_plot(graph,subgraphs;kwargs...)


text:
mutable struct _Linear <: DiffableFunction
    grad_field :: GradField
end

function forward(::_Linear, x, W, b)
    return x * W .+ b
end

function backward(f::_Linear, gy)
    x, W, b = f.grad_field.inputs
    gx = matmul(gy, transpose(W))
    gW = matmul(transpose(x), gy)
    gb = sum_to(gy, size(b.values))
    return gx, gW, gb
end

linear(x, W, b) = _Linear(GradField())(x, W, b)


text:
#-----------------------------------------------------------------------------------------------------
# Projected Basis
#-----------------------------------------------------------------------------------------------------
"""
    ProjectedBasis

Basis for subspace that is spanned only by product states.

Properties:
-----------
- dgt : Vector{Int}
- I   : Vector{Int}, list of indicies.
- B   : Int, levels on each site.
"""
struct ProjectedBasis <: AbstractBasis
    dgt::Vector{Int}
    I::Vector{Int}
    B::Int
end
eltype(::ProjectedBasis) = Int
#-----------------------------------------------------------------------------------------------------
function change!(b::ProjectedBasis, i::Integer)::Int
    change!(b.dgt, b.I[i], base=b.B)
    1
end
#-----------------------------------------------------------------------------------------------------
function index(b::ProjectedBasis)::Tuple{Int, Int}
    ind = binary_search(b.I, index(b.dgt, base=b.B))
    1, ind
end
#-----------------------------------------------------------------------------------------------------
size(b::ProjectedBasis, i::Integer) = (i == 2 || i == 1) ? length(b.I) : 1
size(b::ProjectedBasis) = (l=length(b.I); (l,l))

#-----------------------------------------------------------------------------------------------------
# Construct basis
#-----------------------------------------------------------------------------------------------------
export projectedbasis
function projectedbasis(f, L::Integer; base::Integer=2)
    dgt = zeros(Int, L)
    I = []
    for i = 1:base^L
        change!(dgt, i, base=base)
        if f(dgt)
            append!(I, i)
        end
    end
    ProjectedBasis(dgt, I, base)
end

#-----------------------------------------------------------------------------------------------------
# Helper functions
#-----------------------------------------------------------------------------------------------------
"""
    binary_search(list::AbstractVector{<:Integer}, i<:Integer)

Return the position of i in a sorted list using binary search algorithm.
"""
function binary_search(list::AbstractVector{<:Integer}, i::Integer)
    l, r = 1, length(list)
    c = (l+r) ÷ 2
    while true
        t = list[c]
        if i < t
            r = c - 1
        elseif i > t
            l = c + 1
        else
            break
        end
        if l > r
            return 0
        else
            c = (l+r) ÷ 2
        end
    end
    c
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CustomResourceColumnDefinition specifies a column for server side printing.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition(;
        description=nothing,
        format=nothing,
        jsonPath=nothing,
        name=nothing,
        priority=nothing,
        type=nothing,
    )

    - description::String : description is a human readable description of this column.
    - format::String : format is an optional OpenAPI type definition for this column. The 'name' format is applied to the primary identifier column to assist in clients identifying column is the resource name. See https://github.com/OAI/OpenAPI-Specification/blob/master/versions/2.0.md#data-types for details.
    - jsonPath::String : jsonPath is a simple JSON path (i.e. with array notation) which is evaluated against each custom resource to produce the value for this column.
    - name::String : name is a human readable name for the column.
    - priority::Int32 : priority is an integer defining the relative importance of this column compared to others. Lower numbers are considered higher priority. Columns that may be omitted in limited space scenarios should be given a priority greater than 0.
    - type::String : type is an OpenAPI type definition for this column. See https://github.com/OAI/OpenAPI-Specification/blob/master/versions/2.0.md#data-types for details.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition <: SwaggerModel
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    format::Any # spec type: Union{ Nothing, String } # spec name: format
    jsonPath::Any # spec type: Union{ Nothing, String } # spec name: jsonPath
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    priority::Any # spec type: Union{ Nothing, Int32 } # spec name: priority
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition(;description=nothing, format=nothing, jsonPath=nothing, name=nothing, priority=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("format"), format)
        setfield!(o, Symbol("format"), format)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("jsonPath"), jsonPath)
        setfield!(o, Symbol("jsonPath"), jsonPath)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition = Dict{Symbol,Symbol}(Symbol("description")=>Symbol("description"), Symbol("format")=>Symbol("format"), Symbol("jsonPath")=>Symbol("jsonPath"), Symbol("name")=>Symbol("name"), Symbol("priority")=>Symbol("priority"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition = Dict{Symbol,String}(Symbol("description")=>"String", Symbol("format")=>"String", Symbol("jsonPath")=>"String", Symbol("name")=>"String", Symbol("priority")=>"Int32", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition)
    (getproperty(o, Symbol("jsonPath")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition }, name::Symbol, val)
end


text:
"""This collection of methods and definitions allow for easy location determination and comparison between differnt points on any sphere defined by their latitude and longitude."""
module EarthSphere
abstract Coordinate
#A coordinate describes a point in any coordinate system.

abstract LatLong_Coord <: Coordinate
#This is a coordinate in a standard latitude-longitude coordinate system

type Degree
	#Abstraction for dealing with sexagesimal units of subtending the circle
	degree::Float64
	minute::Float64
	second::Float64
end
export Degree

type Sphere_LatLong_Coord <: LatLong_Coord
	#Deals with a spherical latlong coordinate, similar to that used on good old planet earth
	latitude::Float64
	longitude::Float64
	#Example bearing: SW, NE
	bearing::String
	radius::Float64 # 6371000 Meters for earth
end
function Earth_LatLong_Coord(latitude::Degree,longitude::Degree,bearing::String)
	latnumber::Float64 = degree_to_decimal_degree(latitude)
	longnumber::Float64 = degree_to_decimal_degree(longitude)
	Sphere_LatLong_Coord(latnumber,longnumber,bearing,6371)
end
export Earth_LatLong_Coord

Earth_LatLong_Coord(latitude::Float64,longitude::Float64,bearing::Float64) = Sphere_LatLong_Coord(latitude,longitude,bearing,6371)
function Earth_LatLong_Coord(latitude::Float64,longitude::Float64)
	if sign(latitude) == 1
		inclination = "N"
	elseif sign(latitude) == -1
		inclination = "S"
	else
		inclination = ""
	end
	if sign(longitude) == 1
		azimuth = "E"
	elseif sign(longitude) == -1
		azimuth = "W"
	else
		azimuth = ""
	end
	my_sphere = Sphere_LatLong_Coord(abs(latitude),abs(longitude),inclination*azimuth,6371)
	return my_sphere
end

type Spherical_Coord <: Coordinate
	#Mathematical formulation for a 3D spherical coordinate system
	radius::Float64
	theta::Float64
	phi::Float64
end
export Spherical_Coord

function sign_lat_long(latitudeangle::Float64,longitudeangle::Float64,bearing::String)
	"""Returns a signed version of latitude and longitude based on the bearing given,
	uses a coordinate system where east is positive, west is negative, north is positive
	south is negative."""
	"""This will break if you feed it a bearing like "ES" or "WN", but such a case would be
	trivial to build up, build against, and is also not idiomatic"""
	"""UPDATE: This will no longer break if you do that thanks to a cheap and poorly executed hack. Please feel free to make this clean"""
	if length(bearing) == 1
		if bearing[1] == 'S'
			latitudeangle = (-1)*latitudeangle
		elseif bearing[1] == 'W'
			longitudeangle = (-1)*longitudeangle
		else
			return (latitudeangle,longitudeangle)
		end
	end
	if bearing[1] == 'E' | bearing[1] == 'W'
		bearing[1] = temp
		bearing[1] = bearing[2]
		bearing[2] = bearing[1]
	end

	if bearing[1] == 'S'
		latitudeangle = (-1)*latitudeangle
	end
	if bearing[2] == 'W'
		longitudeangle = (-1)*longitudeangle
	end
	return (latitudeangle,longitudeangle)
end

function sign_lat_long(latlong::LatLong_Coord)
	"Signs a latlong coord using the general sign_lat_long"
	(newlat::Float64,newlong::Float64) = sign_lat_long(latlong.latitude,latlong.longitude,latlong.bearing)
end

function degree_to_decimal_degree(degree::Degree)
	"Takes a degree in sextagesimal and returns a decimal degree"
	return degree.degree + degree.minute/60 + degree.second/(60^2)
end

function convert(Spherical_Coord,latlong::LatLong_Coord)
	(signedlat,signedlong) = sign_lat_long(latlong.latitude,latlong.longitude,latlong.bearing)
	phi::Float64 = to_radians(signedlat)
	theta::Float64 =  to_radians(signedlong)
	return Spherical_Coord(latlong.radius,theta,phi)
end
export convert
function tuple(s::Sphere_LatLong_Coord)
	#Takes a Spherical Coordinate and returns the signed latitude and longitude
	first = s.longitude
	second = s.latitude
	if ismatch(r"S",s.bearing)
		first = (-1)*first
	end
	if ismatch(r"W",s.bearing)
		second = (-1)*second
	end
	return (first,second)
end

function to_radians(degree::Real)
	"Takes a degree, returns radians"
	_sign = sign(degree)
	reduced_degree::Float64 = mod(abs(degree),360)
	radians::Float64 = reduced_degree*(pi/180)
	return radians*_sign
end

function get_distance(c1::Spherical_Coord,c2::Spherical_Coord)
	"Gets the distance between two different points on a sphere."
	if c1.radius != c2.radius
		error("These two points are not located on the same sphere")
	end
	#We're using the Vincenty formula, because it appears to be very computationally and algorithmically reasonable.
	#=central_angle = spherical_law_of_cosines(c1,c2)=#
	central_angle::Float64 = vincenty(c1,c2)
	#=central_angle = haversine(c1,c2)=#
	#=central_angle = alternative_haversine(c1,c2)=#
	distance::Float64 = c1.radius*central_angle

	return distance
end
function get_distance(c1::Sphere_LatLong_Coord,c2::Sphere_LatLong_Coord)
	c1s = convert(Spherical_Coord,c1)
	c2s = convert(Spherical_Coord,c2)
	return get_distance(c1s,c2s)::Float64
end

export get_distance

function vincenty(c1::Spherical_Coord,c2::Spherical_Coord)
	#An implementation of the Vincenty formula: http://en.wikipedia.org/wiki/Great-circle_distance
	longdiff::Float64 = c1.theta-c2.theta
	#We could take absolute value of long (the difference of longitudes) but we only ever take the cosine of it, so we can just leave it (cosine is even, sign doesn't matter)
	vin_numerator::Float64 = ( (cos(c2.phi)*sin(longdiff))^2 + ( cos(c1.phi)*sin(c2.phi) - sin(c1.phi)*cos(c2.phi)*cos(longdiff))^2 )^(1/2)
	vin_denominator::Float64 = sin(c1.phi)*sin(c2.phi)+cos(c1.phi)*cos(c2.phi)*cos(longdiff)
	central_angle::Float64 = atan2(vin_numerator,vin_denominator)
	return central_angle
end
vincenty(c1,c2) = vincenty(convert(Spherical_Coord,c1),convert(Spherical_Coord,c2))
spherical_law_of_cosines(c1,c2) = spherical_law_of_cosines(convert(Spherical_Coord,c1),convert(Spherical_Coord,c2))
function spherical_law_of_cosines(c1::Spherical_Coord,c2::Spherical_Coord)
	#Implements Spherical Law of cosines, neither vincenty nor this are getting me close to google maps
	#UPDATE!: It's okay, I just accidentally had latitude for longitude and longitude for latitude, not really a big deal except it broke everything
	longdiff::Float64 = abs(c1.theta-c2.theta)
	central_angle::Float64 = acos(sin(c1.phi)*sin(c2.phi)+cos(c1.phi)*cos(c2.phi)*cos(longdiff))
	return central_angle
end
function haversine(c1::Spherical_Coord,c2::Spherical_Coord)
	longdiff::Float64 = c1.theta-c2.theta
	latdiff::Float64 = c1.phi-c2.phi
	a::Float64 = (sin(latdiff/2))^2 + cos(c1.phi)*cos(c2.phi)*(sin(longdiff/2))^2
	central_angle::Float64 = 2*atan2(sqrt(a),sqrt(1-a))
	return central_angle
end
function alternative_haversine(c1::Spherical_Coord,c2::Spherical_Coord)
	#This is another formulation of the haversine formula from wikipedia

	latdiff::Float64 = abs(c1.phi-c2.phi)
	longdiff::Float64 = abs(c1.theta-c2.theta)
	argument::Float64 = sqrt(sin((latdiff)/2)*sin((latdiff)/2) - cos(c1.phi)*cos(c2.phi)*sin((longdiff)/2)*sin((longdiff)/2))
	central_angle::Float64 = 2*asin(argument)
	return central_angle
end

alternative_haversine(c1,c2) =alternative_haversine(convert(Spherical_Coord,c1),convert(Spherical_Coord,c2))
haversine(c1,c2) = haversine(convert(Spherical_Coord,c1),convert(Spherical_Coord,c2))

end #End module


text:
# Mixture models
distset = [
    MixtureModel([Uniform(), Normal(), LogNormal()]),
    MixtureModel([Uniform(1, 10), Gamma()]),
    MixtureModel([Uniform(1, 100), Normal(1, 1000)]),
    MixtureModel([DiscreteUniform(1, 10), DiscreteUniform(1, 10)])
]

for dist in distset
    println(dist)
    μ = mean(dist)
    σ = std(dist)
    # No convenience call.
    E = expectation(dist)
    @test E(x -> x) ≈ μ
    @test E(x -> x^2) - μ^2 ≈ σ^2
    # Convenience call.
    @test expectation(x -> x, dist) ≈ μ
    @test expectation(x -> x^2, dist) - μ^2 ≈ σ^2
    # Stress tests.
    # Many nodes.
    E2 = expectation(dist, n = 100)
    @test E2(x -> x) ≈ μ
    @test E2(x -> x^2) - μ^2 ≈ σ^2
end

# Addition of expectations
e1 = expectation(Normal())
e2 = expectation(Uniform())
e3 = expectation(Gamma())

@test ((e1 + e2) + e3)(identity) == ((e1 + e2) + e3)(identity) == (e1 + e2 + e3)(identity)
@test (e1 + e2 + e3)(x -> x^2) ≈ 3.333333333333342

# Linear operator behavior.

# Error handling.


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[9], 0.0)
JuMP.set_lower_bound(x[8], 0.0)
JuMP.set_lower_bound(x[7], 0.0)
JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_lower_bound(x[11], 0.0)
JuMP.set_lower_bound(x[10], 0.0)
JuMP.set_lower_bound(x[12], 0.0)
JuMP.set_lower_bound(x[1], 10.0)
JuMP.set_upper_bound(x[1], 350.0)
JuMP.set_lower_bound(x[2], 10.0)
JuMP.set_upper_bound(x[2], 350.0)
JuMP.set_lower_bound(x[3], 10.0)
JuMP.set_upper_bound(x[3], 200.0)
JuMP.set_lower_bound(x[4], 10.0)
JuMP.set_upper_bound(x[4], 200.0)
JuMP.set_upper_bound(x[5], 10.0)
JuMP.set_upper_bound(x[6], 10.0)
JuMP.set_upper_bound(x[7], 10.0)
JuMP.set_upper_bound(x[8], 10.0)
JuMP.set_upper_bound(x[9], 10.0)
JuMP.set_upper_bound(x[10], 10.0)
JuMP.set_upper_bound(x[11], 10.0)
JuMP.set_upper_bound(x[12], 10.0)
JuMP.set_lower_bound(x[13], 150.0)
JuMP.set_upper_bound(x[13], 310.0)
JuMP.set_lower_bound(x[14], 150.0)
JuMP.set_upper_bound(x[14], 310.0)
JuMP.set_lower_bound(x[15], 150.0)
JuMP.set_upper_bound(x[15], 310.0)
JuMP.set_lower_bound(x[16], 150.0)
JuMP.set_upper_bound(x[16], 310.0)

# ----- Constraints ----- #
@constraint(m, e1, x[5]+x[9] == 10.0)
@constraint(m, e2, x[5]-x[6]+x[11] == 0.0)
@constraint(m, e3, x[7]+x[9]-x[10] == 0.0)
@constraint(m, e4, -x[6]+x[7]+x[8] == 0.0)
@constraint(m, e5, -x[10]+x[11]+x[12] == 0.0)
@NLconstraint(m, e6, x[16]*x[11]-x[13]*x[6]+150*x[5] == 0.0)
@NLconstraint(m, e7, x[15]*x[7]-x[14]*x[10]+150*x[9] == 0.0)
@NLconstraint(m, e8, x[6]*x[15]-x[6]*x[13] == 1000.0)
@NLconstraint(m, e9, x[10]*x[16]-x[10]*x[14] == 600.0)
@constraint(m, e10, x[1]+x[15] == 500.0)
@constraint(m, e11, x[2]+x[13] == 250.0)
@constraint(m, e12, x[3]+x[16] == 350.0)
@constraint(m, e13, x[4]+x[14] == 200.0)

# ----- Objective ----- #
@NLobjective(m, Min, (1300*(1000/(0.0333333333333333*x[1]*x[2]+0.166666666666667*x[1]+0.166666666666667*x[2]))^0.6+1300*(600/(0.0333333333333333*x[3]*x[4]+0.166666666666667*x[3]+0.166666666666667*x[4]))^0.6))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
objval = objective_value(m)
println("run time: $run_time")


text:
# NOTE setup too for `cthulhu_ast`
include("sandbox.jl")
using .CthulhuTestSandbox
Revise.track(CthulhuTestSandbox, normpath(@__DIR__, "sandbox.jl"))

@testset "printer test" begin
    _, src, infos, mi, rt, slottypes = process(testf_revise);
    tf = (true, false)

    @testset "codeview: $codeview" for codeview in Cthulhu.CODEVIEWS
        if !@isdefined(Revise)
            codeview == Cthulhu.cthulhu_ast && continue
        end
        @testset "optimize: $optimize" for optimize in tf
            @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
                params = Cthulhu.current_params()
                config = Cthulhu.CONFIG

                io = IOBuffer()
                codeview(io, mi, optimize, debuginfo, params, config)
                @test !isempty(String(take!(io))) # just check it works
            end
        end
    end

    @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
        @testset "iswarn: $iswarn" for iswarn in tf
            @testset "hide_type_stable: $hide_type_stable" for hide_type_stable in tf
                @testset "inline_cost: $inline_cost" for inline_cost in tf
                    io = IOBuffer()
                    Cthulhu.cthulhu_typed(io, debuginfo,
                        src, rt, mi;
                        iswarn, hide_type_stable, inline_cost)
                    @test !isempty(String(take!(io))) # just check it works
                end
            end
        end
    end
end

foo = Core.eval(Module(), Meta.parseall(raw"""
function foo(x, y)
    z = x + y
    while z < 4
        z += 1
    end
    u = (x -> x + z)(x)
    v = Ref{Union{Int, Missing}}(x) + y
    return z * 7
end
"""; filename="foobar.jl"))

function strip_base_linenums(s)
    s = replace(s, r"((boot|int|refvalue|refpointer)\.jl:)\d+" => s"\1")
    return s
end

function filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost)
    return "test_output/foo-$(optimize ? :opt : :unopt)-$(debuginfo)-$(iswarn ? :warn : :nowarn)\
-$(hide_type_stable ? :hide_type_stable : :show_type_stable)-$(inline_cost ? :inline_cost : :nocost)"
end

#=
# to generate test cases
function generate_test_cases(foo)
    outputs = Dict()
    tf = (true, false)
    for optimize in tf
        _, src, infos, mi, rt, slottypes = process(foo, (Int, Int); optimize);
        for (debuginfo, iswarn, hide_type_stable, inline_cost) in Iterators.product(
            instances(Cthulhu.DInfo.DebugInfo), tf, tf, tf,
        )
            !optimize && debuginfo === Cthulhu.DInfo.compact && continue
            !optimize && inline_cost && continue

            s = sprint(; context=:color=>true) do io
                Cthulhu.cthulhu_typed(io, debuginfo,
                                      src, rt, mi;
                                      iswarn, hide_type_stable, inline_cost)
            end
            s = strip_base_linenums(s)
            write(filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost), s)

            in(s, values(outputs)) && (@show((optimize, debuginfo, iswarn, hide_type_stable, inline_cost)); @show(findfirst(==(s), pairs(outputs))))
            @assert !in(s, values(outputs))
            outputs[(optimize, debuginfo, iswarn, hide_type_stable, inline_cost)] = s
        end
    end
end

generate_test_cases(foo)
=#

using DeepDiffs

@testset "IRShow tests" begin
    tf = (true, false)

    @testset "optimize: $optimize" for optimize in tf
        _, src, infos, mi, rt, slottypes = process(foo, (Int, Int); optimize);

        # if coverage is enabled, src may contain `code_coverage_effect` expressions
        if optimize && Base.JLOptions().code_coverage != 0
            for bb in src.cfg.blocks
                for idx in bb.stmts
                    if src.stmts.inst[idx] == Expr(:code_coverage_effect)
                        src.stmts.inst[idx] = nothing
                    end
                end
            end
            src = Cthulhu.dce!(src, mi)
        end

        @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
            @testset "iswarn: $iswarn" for iswarn in tf
                @testset "hide_type_stable: $hide_type_stable" for hide_type_stable in tf
                    @testset "inline_cost: $inline_cost" for inline_cost in tf
                        !optimize && debuginfo === Cthulhu.DInfo.compact && continue
                        !optimize && inline_cost && continue

                        s = sprint(; context=:color=>true) do io
                            Cthulhu.cthulhu_typed(io, debuginfo,
                                                  src, rt, mi;
                                                  iswarn, hide_type_stable, inline_cost)
                        end
                        s = strip_base_linenums(s)

                        ground_truth = read(filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost), String)
                        if Sys.iswindows()
                            ground_truth = replace(ground_truth, "\r\n" => "\n")
                        end
                        @test s == ground_truth
                        s != ground_truth && println(deepdiff(s, ground_truth))
                    end
                end
            end
        end
    end
end

@testset "hide type-stable statements" begin
    let # optimize code
        _, src, infos, mi, rt, slottypes = @eval Module() begin
            const globalvar = Ref(42)
            $process() do
                a = sin(globalvar[])
                b = sin(undefvar)
                return (a, b)
            end
        end
        function prints(; kwargs...)
            io = IOBuffer()
            Cthulhu.cthulhu_typed(io, :none, src, rt, mi; kwargs...)
            return String(take!(io))
        end

        let # by default, should print every statement
            s = prints()
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let # should omit type stable statements
            s = prints(; hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
    end

    let # unoptimize code
        _, src, infos, mi, rt, slottypes = @eval Module() begin
            const globalvar = Ref(42)
            $process(; optimize=false) do
                a = sin(globalvar[])
                b = sin(undefvar)
                return (a, b)
            end
        end
        function prints(; kwargs...)
            io = IOBuffer()
            Cthulhu.cthulhu_typed(io, :none, src, rt, mi; kwargs...)
            return String(take!(io))
        end

        let # by default, should print every statement
            s = prints()
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let # should omit type stable statements
            s = prints(; hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end

        # should work for warn mode
        let
            s = prints(; iswarn=true)
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let
            s = prints(; iswarn=true, hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
    end
end


text:
using QuadGK

"""
    diagnose_buoyancy_flux(b, z, time_period, depth)

Returns the buoyancy flux (m²/s³) required to extract all the buoyancy from a 1D buoyancy profile `b(z)`
down to some z-coordinate `depth` (< 0) over a `time_period` τ (s) by integrating the buoyancy profile
and converting the integrated buoyancy content to an equivalent buoyancy flux:

          1  ⌠ z_surface
    Qb = --- |           [b(z) - b(depth)] dz
          τ  ⌡ depth
"""
function diagnose_buoyancy_flux(b, z, time_period, depth)

    Δz = z[2] - z[1]
    z_surface = last(z)
    @assert first(z) < depth < z_surface

    ℑb = InterpolatedProfile(b, z, Δz, z_surface)
    b₀ = ℑb(depth)

    integrand(z) = ℑb(z) - b₀
    B, err = quadgk(integrand, depth, z_surface)

    return B / time_period
end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 1.5.
#

# this utility function fills up the mantissa
function sum_(bits)
    # start at 0
    output = fl(0)
    # add 1, 2, 4, \dots
    for i in range(0, length=bits)
        output += fl(2^i)
    end
    return output
end

# assert the program’s arguments
if length(ARGS) == 0
    println("error: provide the floating point precision (16, 32 or 64 bits)")
    exit(1)
end

fl_precision_string = ARGS[1]

if fl_precision_string == "16"
    fl = Float16
    # fill up the mantissa
    mach_max = sum_(11)
elseif fl_precision_string == "32"
    fl = Float32
    # fill up the mantissa
    mach_max = sum_(24)
elseif fl_precision_string == "64"
    fl = Float64
    # fill up the mantissa
    mach_max = sum_(53)
else
    println("error: provide one of [16, 32, 64]")
    exit(1)
end
# assert end

# iteratively find the max value
for i in range(1, length=2500)
    # if the next value would be considered ∞ break the loop
    if isinf(fl(mach_max * 2))
        break
    end
    # otherwise continue
    global mach_max = fl(mach_max * 2)
end

# print out the results
o1, o2 = mach_max, floatmax(fl)

println("iterative method = ", o1, "; floatmax(fl) = ", o2, "; equal = ", (o1 == o2))

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
const code_object_to_vk_prefix = Dict(
    :struct => "Vk",
    :function => "vk",
    :constant => "VK",
)
const vulkan_convention_rev = Dict(v => k for (k, v) in vulkan_convention)

vk_prefix(::Type{T}) where {T <: NamingConvention} = code_object_to_vk_prefix[vulkan_convention_rev[T]]

"""Prefix a function with the Vulkan convention prefix vk, Vk or VK_ depending on the code object.
"""
function prefix_vk(name)
   convention = detect_convention(name)
   prefix(convention(name), vk_prefix(convention)).value
end
prefix_vk(name::T) where {T <: NamingConvention} = prefix(name, vk_prefix(typeof(name)))

# remove_vk_prefix(name) = replace(name, Regex("^" * prefix_vk(name)) => "")
remove_vk_prefix(name) = replace(name, r"^(?:vk|Vk|VK_)" => "")


text:
module Encode
using ..Meta, ..Patterns, ..Dispatch
import ..Nodes
import ..Patterns.resultof, ..Nodes.encode

export code_dispatch, seq_dispatch!

# ---- Result: Node type for instantiated results -----------------------------

type Result{T} <: Node{T}
    node::Node{T}
    name::Union{Symbol,Void}
    nrefs::Int
    ex

    Result(node::Node{T}, name) = new(node, name, 1, nothing)
end
Result{T}(node::Node{T}, name) = Result{T}(node, name)
Result(node::Node)             = Result(node, nothing)

resultof(node::Result) = (@assert node.ex != nothing; node.ex)

# ---- sequence!: Create evaluation order, instantiate nodes into Result's ----

typealias ResultsDict Dict{Node,Node}
type Sequence
    intent::Intension
    results::ResultsDict
    namer::Function
    seq::Vector{Node}

    Sequence(i::Intension, r::ResultsDict, namer) = new(i, r, namer, Node[])
end

wrap(s::Sequence, node::Guard, orig_node::Node) = node
function wrap(s::Sequence, node::Node, orig_node::Node)
    Result(node, s.namer === nothing ? nothing : s.namer(orig_node))
end

function sequence!(s::Sequence, node::Node)
    if haskey(s.results, node)
        newnode = s.results[node]
        if isa(newnode, Result); newnode.nrefs += 1; end
        return
    end

    for dep in depsof(s.intent, node);  sequence!(s, dep)  end
    newnode = s.results[node] = wrap(s, subs(s.results, node), node)
    push!(s.seq, newnode)
end

preguard!(results::ResultsDict, g::Guard) = (results[g] = Guard(always))
function provide!(results::ResultsDict, node::Node, ex)
    res = Result(node)
    res.ex = ex
    results[node] = res
end

# ---- sequence decision tree -------------------------------------------------

function seq_dispatch!(d::DNode, methods, hullT::Tuple)
    # todo: extract methods from d instead
    results = ResultsDict()
    for pred in predsof(intension(hullT)); preguard!(results, Guard(pred)); end

    argsyms = []
    for k=1:length(hullT)
        node, name = Nodes.refnode(Nodes.argnode, k), nothing
        for method in methods
            rb = method.sig.rev_bindings
            if haskey(rb, node)
                name = symbol(string(rb[node], '_', method.id))
                break
            end
        end
        if name === nothing;  name = symbol("arg$k");  end

        push!(argsyms, name)
        provide!(results, node, name)
    end

    seq_dispatch!(results, d)
    argsyms
end

seq_dispatch!(results::ResultsDict, d::DNode) = nothing
function seq_dispatch!(results::ResultsDict, m::MethodCall)
    s = Sequence(domainof(m.m), results, make_namer([m.m]))
    for node in values(m.m.sig.bindings);  sequence!(s, node)  end
    m.bind_seq = s.seq
    m.bindings = Node[results[node] for node in m.m.bindings]
end
function seq_dispatch!(results::ResultsDict, d::Decision)
    results_fail = copy(results)
    s = Sequence(d.domain, results, make_namer(d.methods))
    for p in predsof(d.domain);  sequence!(s, Guard(p))  end

    k=findfirst(node->isa(node,Guard), s.seq)
    d.pre = s.seq[1:(k-1)]
    d.seq = s.seq[k:end]

    seq_dispatch!(results, d.pass)
    seq_dispatch!(results_fail, d.fail)
end

# ---- decision tree -> dispatch code -----------------------------------------

blockwrap(ex)       = ex
blockwrap(exprs...) = Expr(:block, exprs...)

function code_dispatch(::NoMethodNode)
    :( error("No matching pattern method found") )
end
function code_dispatch(m::MethodCall)
    prebind = encoded(m.bind_seq)
    args = Any[resultof(node) for node in m.bindings]
    blockwrap(prebind...,
         Expr(:call, quot(m.m.body), args...))
end
function code_dispatch(d::Decision)
    pre = encoded(d.pre)
    pred = code_predicate(d.seq)
    pass = code_dispatch(d.pass)
    fail = code_dispatch(d.fail)
    code = Expr(:if, pred, pass, fail)
    blockwrap(pre..., code)
end

# ---- generate code from (instantiated) node sequence ------------------------

function code_predicate(seq::Vector{Node})
    code = []
    pred = nothing
    for node in seq
        if isa(node, Guard)
            factor = isempty(code) ? node.pred.ex : Expr(:block, code..., node.pred.ex)
            pred = pred === nothing ? factor : (:($pred && $factor))
            code = []
        else
            encode!(code, node)
        end
    end
    pred
end

function encoded(seq::Vector{Node})
    code = []
    for node in seq;  encode!(code, node);  end
    code
end

encode!(code::Vector, ::Guard) = error("Undefined!")
function encode!(code::Vector, node::Result)
    if node.ex != nothing;  return  end
    ex = encode(node.node)

    if isa(ex, Symbol) || node.nrefs == 1
        node.ex = ex
    else
        var = node.name === nothing ? gensym("t") : node.name
        node.ex = var
        push!(code, :( $var = $ex ))
    end
end

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to convert anomalies related to the orbit.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export M_to_E, M_to_f
export E_to_f, E_to_M
export f_to_E, f_to_M

################################################################################
#                              From Mean Anomaly
################################################################################

"""
    M_to_E(e::T1, M::T2; kwargs...) where {T1, T2}

Compute the eccentric anomaly (0,2π) \\[rad] given the eccentricity `e` and the
mean anomaly `M` [rad].

This function uses the Newton-Raphson algorithm to solve the Kepler's equation.

# Keywords

- `tol::Union{Nothing, Number}`: Tolerance to accept the solution from
    Newton-Raphson algorithm. If `tol` is `nothing`, then it will be
    `eps(T)`, where `T` is a floating-point type obtained from the promotion of
    `T1` and `T2`. (**Default** = `nothing`)
- `max_iterations::Number`: Maximum number of iterations allowed for the
    Newton-Raphson algorithm. If it is lower than 1, then it is set to 10.
    (**Default** = 10)
"""
@inline function M_to_E(
    e::T1,
    M::T2;
    max_iterations::Integer = 10,
    tol::Union{Nothing, Number} = nothing
) where {T1, T2}
    T = float(promote_type(T1, T2))

    # Compute the eccentric anomaly using the Newton-Raphson method.
    # ==============================================================

    # Make sure that M is in the interval [0,2π].
    M = mod(M, T(2π))

    # Initial guess.
    #
    # See [1, p. 75].
    E = (M > π) ? M - e : M + e

    sin_E, cos_E = sincos(E)

    # Check the tolerance.
    if tol === nothing
        δ = eps(T)
    else
        δ = T(tol)
    end

    # Check the maximum number of iterations.
    if max_iterations < 1
        max_iterations = 10
    end

    # Newton-Raphson iterations.
    for i in 1:max_iterations
        abs(E - e*sin_E - M) ≤ δ && break
        E = E - (E - e * sin_E - M) / (1 - e * cos_E)
        sin_E, cos_E = sincos(E)
    end

    # Return the eccentric anomaly in the interval [0, 2π].
    return mod(E, T(2π))
end

"""
    M_to_f(e::T1, M::T2; kwargs...) where {T1, T2}

Compute the true anomaly (0,2π) \\[rad] given the eccentricity `e` and the mean
anomaly `M` [rad].

This function uses the Newton-Raphson algorithm to solve the Kepler's equation.

# Keywords

- `tol::Union{Nothing, Number}`: Tolerance to accept the solution from
    Newton-Raphson algorithm. If `tol` is `nothing`, then it will be
    `eps(T)`, where `T` is a floating-point type obtained from the promotion of
    `T1` and `T2`. (**Default** = `nothing`)
- `max_iterations::Number`: Maximum number of iterations allowed for the
    Newton-Raphson algorithm. If it is lower than 1, then it is set to 10.
    (**Default** = 10)
"""
@inline function M_to_f(e::T1, M::T2; kwargs...) where {T1, T2}
    # Compute the eccentric anomaly.
    E = M_to_E(e, M; kwargs...)

    # Compute the true anomaly in the interval [0,2π].
    return E_to_f(e, E)
end

################################################################################
#                            From Eccentric Anomaly
################################################################################

"""
    E_to_f(e::T1, E::T2) where {T1, T2}

Compute the true anomaly (0,2π) \\[rad] given the eccentricity `e` and the
eccentric anomaly `E` [rad].
"""
@inline function E_to_f(e::T1, E::T2) where {T1, T2}
    T = float(promote_type(T1, T2))

    sin_Eo2, cos_Eo2 = sincos(E / 2)

    # Compute the true anomaly in the interval [0, 2*π].
    return mod(2atan(sqrt(1 + e) * sin_Eo2, sqrt(1 - e) * cos_Eo2), T(2π))
end

"""
    E_to_M(e::T1, E::T2) where {T1, T2}

Compute the mean anomaly (0,2π) \\[rad] given the eccentricity `e` and the
eccentric anomaly `E` [rad].
"""
@inline function E_to_M(e::T1, E::T2) where {T1, T2}
    T = float(promote_type(T1, T2))
    return mod(E - e * sin(E), T(2π))
end

################################################################################
#                              From True Anomaly
################################################################################

"""
    f_to_E(e::T1, f::T2) where {T1, T2}

Compute the eccentric anomaly (0,2π) \\[rad] given the eccentricity `e` and
the true anomaly `f` [rad].
"""
@inline function f_to_E(e::T1, f::T2) where {T1, T2}
    T = float(promote_type(T1, T2))
    sin_fo2, cos_fo2 = sincos(f / 2)
    return mod(2atan(sqrt(1 - e) * sin_fo2, sqrt(1 + e) * cos_fo2), T(2π) )
end

"""
    f_to_M(e::T1, f::T2) where {T1, T2}

Compute the mean anomaly (0,2π) \\[rad] given the eccentricity `e` and the
true anomaly `f` [rad].
"""
@inline function f_to_M(e::T1, f::T2) where {T1, T2}
    # Compute the eccentric anomaly.
    E = f_to_E(e, f)

    # Compute the true anomaly in the interval [0, 2π].
    return E_to_M(e, E)
end