text:
#=
Train the model using gradient descent.

Training parameters
num_epochs - number of epochs
α - learning rate
μ - momentum
TODO - decay - function of epochs
TODO - early stopping

Store and display parameters
last_train_every - selected epoch intervals for last training batch
full_train_every - selected epoch intervals for full training set
val_every - selected epoch intervals for validation
display - choose whether to display progress
store - choose whether to store progress
=#

################################################################################
#
# Main Training Function
#
################################################################################

function train(net::NeuralNet, num_epochs::Integer, α::Real, μ::Real;
               nesterov::Bool = true, shuffle::Bool=false,
               last_train_every::Integer = 1,
               full_train_every::Integer = num_epochs,
               val_every::Integer = num_epochs)

    # Some counts
    num_obs = size(net.data.X_train)[end]
    num_batches = div(num_obs, net.batch_size)

    # Floating point data type
    F = net.dtype

    # Learning parameters - convert type
    α = F(α)
    μ = F(μ)

    # Nesterov parameters
    a = μ^2
    b = (μ + F(1)) * α

    # Choose momentum gradient descent or Nesterov accelerated gradient descent
    nesterov ? (update! = nag_update!) : (update! = mgd_update!)

    # Choose coffee break function - either regression or classification problem
    coffee_break = choose_coffee_break(net)

    # Permutation of observations - to be shuffled if selected
    perm = collect(1:num_obs)

    # Loop over epochs
    for epoch in 1:num_epochs

        # Shuffle permutation (if selected) and initialise batch unit range
        shuffle && shuffle!(perm)
        batch_range = 1:net.batch_size

        # Loop over mini batches
        for batch in 1:num_batches

            # Update inputs & target with next batch and increment range
            batch_ids = perm[batch_range]
            net.layers[1].A[:] = viewbatch(net.data.X_train, batch_ids)
            net.layers[end].y[:] = viewbatch(net.data.y_train, batch_ids)
            batch_range += net.batch_size

            # Forward propogation through each layer
            fwd_prop!(net)

            # Back propogation through each layer to calculate gradients
            bwd_prop!(net)

            # Update parameters using the gradients
            for ℓ in net.layers[2:end]
                update!(ℓ, μ, α, a, b)
            end

        end

        # Display and store progress updates
        # Last training loss
        last_batch_cost(net, epoch, last_train_every)
        # Coffee break - full training and validation loss
        coffee_break(net, epoch, num_epochs, full_train_every, val_every)

    end

    # Completed training - show final results
    show_completed(net)

    return
end

################################################################################
#
# Momentum Gradient Descent Update
#
################################################################################

function mgd_update!{F<:AbstractFloat}(W::Array{F}, W_vel::Array{F},
                                       ∇W::Array{F}, μ::F, α::F)
    for i in eachindex(W)
        @inbounds W_vel[i] = μ * W_vel[i] - α * ∇W[i]
        @inbounds W[i] += W_vel[i]
    end
end

function mgd_update!{F<:AbstractFloat}(ℓ::NetLayer, μ::F, α::F, a::F, b::F)
    mgd_update!(ℓ.b, ℓ.b_vel, ℓ.∇b, μ, α)
    mgd_update!(ℓ.W, ℓ.W_vel, ℓ.∇W, μ, α)
end

function mgd_update!{F<:AbstractFloat}(ℓ::WordEmbeddingLayer, μ::F, α::F, a::F, b::F)
    mgd_update!(ℓ.W, ℓ.W_vel, ℓ.∇W, μ, α)
end

function mgd_update!{F<:AbstractFloat}(ℓ::PoolLayer, μ::F, α::F, a::F, b::F)
end

################################################################################
#
# Nesterov's Accelerated Gradient Descent Update
#
################################################################################

function nag_update!{F<:AbstractFloat}(W::Array{F}, W_vel::Array{F},
                                       ∇W::Array{F}, μ::F, α::F, a::F, b::F)
    for i in eachindex(W)
        @inbounds W[i] += a * W_vel[i] - b * ∇W[i]
        @inbounds W_vel[i] = μ * W_vel[i] - α * ∇W[i]
    end
end

function nag_update!{F<:AbstractFloat}(ℓ::NetLayer, μ::F, α::F, a::F, b::F)
    nag_update!(ℓ.b, ℓ.b_vel, ℓ.∇b, μ, α, a, b)
    nag_update!(ℓ.W, ℓ.W_vel, ℓ.∇W, μ, α, a, b)
end

function nag_update!{F<:AbstractFloat}(ℓ::WordEmbeddingLayer, μ::F, α::F, a::F, b::F)
    nag_update!(ℓ.W, ℓ.W_vel, ℓ.∇W, μ, α, a, b)
end

function nag_update!{F<:AbstractFloat}(ℓ::PoolLayer, μ::F, α::F, a::F, b::F)
end

################################################################################
#
# Display and store progress
#
################################################################################

########################################
# Last training batch
########################################

# Time
t() = Dates.format(now(), "HH:MM:SS")

# Last batch training error - cost with regularisation
function last_batch_cost(net, epoch, last_train_every)

    if mod(epoch, last_train_every) == 0
        # Cost
        L = net.layers[end]
        J = L.cost(L.A, L.y) / net.batch_size

        # Add regularisation cost
        J += net.reg_cost(net)

        # Add sparsity cost
        J += sparse_cost(net)

        # Display
        print(t())
        @printf(" : Epoch %d, last batch training error (with regⁿ) - %.3f\n",
                epoch, J)

        # Store
        push!(net.last_train_loss, J)
    end

end

########################################
# Coffee Breaks
########################################

# Choos classification or regression
function choose_coffee_break(net)
    L = net.layers[end]
    typeof(L) <: LinOutputLayer ? coffee_break_r : coffee_break_c
end

# Coffee break for classification task - shows accuracy
function coffee_break_c(net, epoch, num_epochs, full_train_every, val_every)

    # Don't run coffee break after last epoch
    epoch == num_epochs && return

    # Both training and validation sets
    if mod(epoch, full_train_every) == 0 && mod(epoch, val_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of training set
        J, acc = loss_and_accuracy(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f", J)
        @printf("  |  Training accuracy - %.1f\n", acc)
        push!(net.full_train_loss, J)

        # Cost and display of validation set
        J, acc = val_loss_and_accuracy(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f", J)
        @printf("  |  Validation accuracy - %.1f\n\n", acc)
        push!(net.val_loss, J)

    # Training set only
    elseif mod(epoch, full_train_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of training set
        J, acc = loss_and_accuracy(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f", J)
        @printf("  |  Training accuracy - %.1f\n\n", acc)
        push!(net.full_train_loss, J)

    # Validation set only
    elseif mod(epoch, val_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of validation set
        J, acc = val_loss_and_accuracy(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f", J)
        @printf("  |  Validation accuracy - %.1f\n\n", acc)
        push!(net.val_loss, J)

    end
end

# Coffee break for regression task - only loss (no accuracy)
function coffee_break_r(net, epoch, full_train_every, val_every)

    # Don't run coffee break after last epoch
    epoch == num_epochs && return

    # Both training and validation sets
    if mod(epoch, full_train_every) == 0 && mod(epoch, val_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of training set
        J, acc = loss(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f\n", J)

        # Cost and display of validation set
        J, acc = val_loss_and_accuracy(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f\n\n", J)

    # Training set only
    elseif mod(epoch, full_train_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of training set
        J, acc = loss_and_accuracy(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f\n\n", J)

    # Validation set only
    elseif mod(epoch, val_every) == 0
        @printf("\nCoffee break:\n")

        # Cost and display of validation set
        J, acc = val_loss_and_accuracy(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f\n\n", J)

    end
end

################################################################################
#
# Completed training - results
#
################################################################################

function show_completed(net)
    if typeof(net.layers[end]) <: LinOutputLayer
        show_completed_regression(net)
    else
        show_completed_classification(net)
    end
end

function show_completed_classification(net)

    @printf "\nCompleted Training:\n"

    if typeof(net.data) <: TData

        # Cost and display of training set
        J, acc = loss_and_accuracy(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f", J)
        @printf("  |  Training accuracy - %.1f\n", acc)
        push!(net.full_train_loss, J)

    else
        # Cost and display of training set
        J, acc = loss_and_accuracy(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f", J)
        @printf("  |  Training accuracy - %.1f\n", acc)
        push!(net.full_train_loss, J)

        # Cost and display of validation set
        J, acc = val_loss_and_accuracy(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f", J)
        @printf("  |  Validation accuracy - %.1f\n", acc)
        push!(net.val_loss, J)
    end
end

function show_completed_regression(net)

    @printf "\nCompleted Training:\n"

    if typeof(net.data) <: TData

        # Cost and display of training set
        J = loss(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f\n", J)
        push!(net.full_train_loss, J)

    else
        # Cost and display of training set
        J = loss(net, net.data.X_train, net.data.y_train)
        @printf("Training error (with regⁿ) - %.3f\n", J)
        push!(net.full_train_loss, J)

        # Cost and display of validation set (without regularisation)
        J = val_loss(net, net.data.X_val, net.data.y_val)
        @printf("Validation error (without regⁿ) - %.3f\n", J)
        push!(net.val_loss, J)
    end
end


text:
include("attributes.jl")
include("mechanical-problems.jl")
include("chemical-problems.jl")
include("power-electronics-problems.jl")
include("power-system-problems.jl")
include("process-design-and-synthesis.jl")

# Mechanical problems
export pressure_vessel,
vibrating_platform,
two_bar_Truss_design_problems,
weldan_beam_design,
disc_brake_design,
speed_reducer_design,
gear_train_design,
car_side_impact_design,
four_bar_plane_truss,
two_bar_plane_truss,
water_tesource_management,
simply_supported_I_beam_design,
gear_box_design,
multiple_disk_clutch_brake_design,
spring_design,
cantilever_beam_design,
bulk_carriers_design,
front_rail_design,
multi_product_batch_plant,
hydro_static_thrust_bearing_design,
crash_energy_management_for_high_speed_train

# chemical problems
export haverlys_pooling_problem,
reactor_network_design,
heat_exchanger_network_design

# Process Design and Synthesis
export process_synthesis,
process_sythesis_and_design,
process_flow_sheeting_problem,
two_reactor_problem,
process_synthesis_problem

# Power electronics Problem
export SOPWM_for_3_level_invereters,
SOPWM_for_5_level_Inverters,
SOPWM_for_7_level_inverters,
SOPWM_for_9_level_inverters,
SOPWM_for_11_level_inverters,
SOPWM_for_13_level_inverters

export sizing_of_single_phase_distributed_generation_1,
 sizing_of_single_phase_distributed_generation_2,
 sizing_of_single_phase_distributed_generation_3,
 sizing_of_single_phase_distributed_generation_4,
 active_and_reactive_power_loss,
 optimal_power_flow_1,
 optimal_power_flow_2,
 optimal_power_flow_3,
 optimal_power_flow_4,
 optimal_power_flow_5,
 optimal_power_flow_6,
 minimization_of_active_power_loss_in_islanded_microgrids_1,
 minimization_of_active_power_loss_in_islanded_microgrids_2,
 minimization_of_active_power_loss_in_islanded_microgrids_3,
 power_distribution_system_planning

export get_RW_MOP_problem

ext: jl


text:
@testset "TigerProblemEnv" begin

    rng = StableRNG(123)
    obs_prob = 0.85
    env = TigerProblemEnv(; rng = rng, obs_prob = obs_prob)

    RLBase.test_interfaces!(env)

    rewards = []
    for _ in 1:50_000
        @test state(env) in state_space(env)
        r = 0.0
        while !is_terminated(env)
            env(rand(rng, action_space(env)))
            r += reward(env)
        end
        push!(rewards, r)
        reset!(env)
    end

    # R̂ = 1/3 (-1 + R̂) + 2 * 1/3 (-100/2 + 10 / 2)
    expectation = -45.5
    @test isapprox(mean(rewards), expectation; atol = 0.5)
end


text:
# SAM Record
# ==========

type SAMRecord
    name::Bio.StringFields.StringField
    flag::UInt16
    refname::Bio.StringFields.StringField
    pos::Int64
    mapq::Int8
    cigar::Bio.StringFields.StringField
    next_refname::Bio.StringFields.StringField
    next_pos::Int64
    tlen::Int32
    seq::DNASequence
    qual::Vector{UInt8}
    optional_fields::Dict{String,Any}
end

function SAMRecord()
    return SAMRecord("", 0x0000, "*", 0, 0, "*", "*", 0, 0, "", UInt8[], Dict())
end

function Base.show(io::IO, rec::SAMRecord)
    println(summary(rec), ':')
    println(io, "reference name: ", refname(rec))
    println(io, "next reference name: ", nextrefname(rec))
    println(io, "position: ", position(rec))
    println(io, "next position: ", nextposition(rec))
    println(io, "mapping quality: ", mappingquality(rec))
    println(io, "flag: ", flag(rec))
    println(io, "template length: ", templatelength(rec))
    println(io, "sequence name: ", seqname(rec))
    println(io, "CIGAR string: ", cigar(rec))
    println(io, "sequence: ", sequence(rec))
    println(io, "base qualities: ", qualities(rec))
      print(io, "optional fields: ", rec.optional_fields)
end

function Base.copy(rec::SAMRecord)
    return deepcopy(rec)
end

function ismapped(r::SAMRecord)
    return r.pos != 0
end

function seqname(r::SAMRecord)
    return r.name
end

function flag(r::SAMRecord)
    return r.flag
end

function refname(r::SAMRecord)
    return r.refname
end

function nextrefname(r::SAMRecord)
    return r.next_refname
end

function Base.position(r::SAMRecord)
    return r.pos
end

function nextposition(r::SAMRecord)
    return r.next_pos
end

function mappingquality(r::SAMRecord)
    return r.mapq
end

function templatelength(r::SAMRecord)
    return r.tlen
end

function cigar(r::SAMRecord)
    return r.cigar
end

function sequence(r::SAMRecord)
    return r.seq
end

function qualities(r::SAMRecord)
    return r.qual
end

function Base.getindex(r::SAMRecord, field::AbstractString)
    return r.optional_fields[field]
end


text:
type MALPDeployment <: DeploymentModel
    m::JuMP.Model
    x::Vector{JuMP.Variable}
end
deployment(m::MALPDeployment) = Int[round(Int,x) for x in JuMP.getValue(m.x)]

function MALPDeployment(p::DeploymentProblem,
                        q::Float64; # busy fraction;
                        α::Float64 = 0.99, # reliability level
                        tol = params.δ,
                        solver = GurobiSolver(OutputFlag=0))
    demand = vec(mean(p.demand[p.train,:],1))
    @assert length(demand) == p.nregions

    I = 1:p.nlocations
    J = 1:p.nregions
    b = ceil(Int, log(1-α)/log(q))

    m = JuMP.Model(solver=solver)
    JuMP.@variable(m, x[1:p.nlocations] >= 0, Int)
    JuMP.@variable(m, z[1:p.nregions, 1:b], Bin)

    JuMP.@objective(m, Max, sum(demand[j]*z[j,b], j in J))

    JuMP.@constraint(m, sum(x[i] for i in I) <= p.nambulances)
    for j in J # coverage over all regions
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= 1)
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= sum(z[j,k] for k in 1:b))
        for k in 2:b
            JuMP.@constraint(m, z[j,k] <= z[j,k-1])
        end
    end

    MALPDeployment(m, x)
end

solve(model::MALPDeployment) = JuMP.solve(model.m)
solve(model::MALPDeployment, p::DeploymentProblem) = JuMP.solve(model.m)


text:
using StaticArrays, NLSolvers, LinearAlgebra
OPT_PROBS = Dict()

#### Himmelblau

OPT_PROBS["himmelblau"] = Dict()
OPT_PROBS["himmelblau"]["array"] = Dict()

function himmelblau!(x)
    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
    return fx
end
function himmelblau_batched_f!(X)
    F = map(himmelblau!, X)
    return F
end
function himmelblau_batched_f!(F, X)
    map!(himmelblau!, F, X)
    return F
end
function himmelblau_g!(∇f, x)
    ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
    44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
    4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    ∇f
end

function himmelblau_fg!(∇f, x)
    ∇f = himmelblau_g!(∇f, x)
    fx = himmelblau!(x)
    return fx, ∇f
end
function himmelblau_fgh!(∇f, ∇²f, x)
    ∇²f = himmelblau_h!(∇²f, x)
    fx, ∇f = himmelblau_fg!(∇f, x)
    return fx, ∇f, ∇²f
end
function himmelblau_h!(∇²f, x)
    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
    ∇²f[2, 1] = ∇²f[1, 2]
    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
    return ∇²f
end
function himmelblau_hv!(hv, x, v)
    hv[1] = (12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0)*v[1] + (4.0 * x[1] + 4.0 * x[2])*v[2]
    hv[2] =  (4.0 * x[1] + 4.0 * x[2])*v[1] + (2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0)*v[2]
    return hv
end

OPT_PROBS["himmelblau"]["array"]["x0"] = [3.0, 1.0]
OPT_PROBS["himmelblau"]["array"]["mutating"] = ScalarObjective(himmelblau!, himmelblau_g!, himmelblau_fg!, himmelblau_fgh!, himmelblau_h!, himmelblau_hv!, himmelblau_batched_f!, nothing)

OPT_PROBS["himmelblau"]["staticarray"] = Dict()

function himmelblau_f(x)
    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
    return fx
end

function himmelblau_fg(∇f, x)
    ∇f = himmelblau_g(∇f, x)
    fx = himmelblau_f(x)
    return fx, ∇f
end

function himmelblau_g(∇f, x)
    ∇f1 = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
    44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    ∇f2 = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
    4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    ∇f = @SVector([∇f1, ∇f2])
    return ∇f
end

function himmelblau_fgh(∇f, ∇²f, x)
    ∇²f11 = 12.0 * x[1]^2 + 4.0 * x[2] - 42.0
    ∇²f12 = 4.0 * x[1] + 4.0 * x[2]
    ∇²f21 = 4.0 * x[1] + 4.0 * x[2]
    ∇²f22 = 12.0 * x[2]^2 + 4.0 * x[1] - 26.0
    ∇²f = @SMatrix([∇²f11 ∇²f12; ∇²f21 ∇²f22])
    fx, ∇f = himmelblau_fg(∇f, x)

    return fx, ∇f, ∇²f
end

sl = @SVector([3.0,1.0])
state0newton = (@SVector([2.0,2.0]), I+sl*sl')
OPT_PROBS["himmelblau"]["staticarray"]["x0"] = sl
OPT_PROBS["himmelblau"]["staticarray"]["state0"] = state0newton
OPT_PROBS["himmelblau"]["staticarray"]["static"] = ScalarObjective(himmelblau_f, himmelblau_g, himmelblau_fg, himmelblau_fgh, nothing, nothing, nothing, nothing)

### Exponential

exponential!(x) = exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
function exponential_g!(g, x)
    g[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    g[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
    return g
end
function exponential_h!(H, x)
    H[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    H[1, 2] = 0.0
    H[2, 1] = 0.0
    H[2, 2] = 2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
    return H
end
function exponential_hv!(Hv, x, v)
    Hv[1, 1] = (2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9))*v[1]
    Hv[2, 2] = (2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19))*v[2]
    return Hv
end
function exponential_fg!(g, x)
    fx = exponential!(x)
    g = exponential_g!(g, x)
    return fx, g
end
function exponential_fgh!(g, H, x)
    fx, g = exponential_fg!(g, x)
    H = exponential_h!(H, x)
    return fx, g, H
end

OPT_PROBS["exponential"] = Dict()
OPT_PROBS["exponential"]["array"] = Dict()
# Byttet om på x og H
OPT_PROBS["exponential"]["array"]["x0"] = [0.0, 0.0]
OPT_PROBS["exponential"]["array"]["mutating"] = ScalarObjective(exponential!, exponential_g!, exponential_fg!, exponential_fgh!, exponential_h!, exponential_hv!, nothing, nothing)

OPT_PROBS["laplacian"] = Dict()
OPT_PROBS["laplacian"]["array"] = Dict()
# Byttet om på x og H
OPT_PROBS["laplacian"]["array"]["x0(n)"] = n->zeros(n)

plap(U; n=length(U)) = (n-1) * sum((0.1 .+ diff(U).^2).^2) - sum(U) / (n-1)
plap1(U; n=length(U), dU = diff(U), dW = 4 .* (0.1 .+ dU.^2) .* dU) =
(n - 1) .* ([0.0; dW] .- [dW; 0.0]) .- ones(n) / (n-1)
precond(x::Vector) = precond(length(x))
precond(n::Number) = spdiagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1)) * (n+1)
_f(x) = plap([0;x;0])
function _fg(∇f, x)
    fx = _f(x)
    copyto!(∇f, (plap1([0;x;0]))[2:end-1])
    fx, ∇f
end

OPT_PROBS["laplacian"]["array"]["mutating"] = OptimizationProblem(ScalarObjective(_f, nothing, _fg, nothing, nothing, nothing, nothing, nothing))

function theta(x)
    if x[1] > 0
        return atan(x[2] / x[1]) / (2.0 * pi)
    else
        return (pi + atan(x[2] / x[1])) / (2.0 * pi)
    end
end

function fletcher_powell_f(x)
    theta_x = theta(x)
    fx = 100.0 * ((x[3] - 10.0 * theta_x)^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
    return fx
end
function fletcher_powell_g(∇f, x)
    T = eltype(x)
    theta_x = theta(x)

    if x[1]^2 + x[2]^2 == 0
        dtdx1 = T(0)
        dtdx2 = T(0)
    else
        dtdx1 = - x[2] / ( T(2) * pi * ( x[1]^2 + x[2]^2 ) )
        dtdx2 =   x[1] / ( T(2) * pi * ( x[1]^2 + x[2]^2 ) )
    end
    ∇f1 = -2000.0*(x[3]-10.0*theta_x)*dtdx1 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 )
    ∇f2 = -2000.0*(x[3]-10.0*theta_x)*dtdx2 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 )
    ∇f3 =  200.0*(x[3]-10.0*theta_x) + 2.0*x[3];
    ∇f = @SVector[∇f1, ∇f2, ∇f3]

    return ∇f
end
function fletcher_powell_fg(∇f, x)
    T = eltype(x)
    theta_x = theta(x)

    if x[1]^2 + x[2]^2 == 0
        dtdx1 = T(0)
        dtdx2 = T(0)
    else
        dtdx1 = - x[2] / ( T(2) * pi * ( x[1]^2 + x[2]^2 ) )
        dtdx2 =   x[1] / ( T(2) * pi * ( x[1]^2 + x[2]^2 ) )
    end
    ∇f1 = -2000.0*(x[3]-10.0*theta_x)*dtdx1 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 )
    ∇f2 = -2000.0*(x[3]-10.0*theta_x)*dtdx2 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 )
    ∇f3 =  200.0*(x[3]-10.0*theta_x) + 2.0*x[3];
    ∇f = @SVector[∇f1, ∇f2, ∇f3]

    fx = 100.0 * ((x[3] - 10.0 * theta_x)^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

    return fx, ∇f
end
const sv3 = @SVector[0.1,0.0,0.0]

const fp_static = ScalarObjective(fletcher_powell_f, fletcher_powell_g, fletcher_powell_fg, nothing, nothing, nothing, nothing, nothing)

const state0 = (@SVector[-0.5, 0.0, 0.0], I+sv3*sv3')
OPT_PROBS["fletcher_powell"] = Dict()
OPT_PROBS["fletcher_powell"]["staticarray"] = Dict()
# Byttet om på x og H
OPT_PROBS["fletcher_powell"]["staticarray"]["x0"] = state0
OPT_PROBS["fletcher_powell"]["staticarray"]["static"] = OptimizationProblem(fp_static; inplace=false)


text:
_MsgJointLikelihood(;relatives::MsgRelativeType=MsgRelativeType(),priors::MsgPriorType=MsgPriorType() ) = _MsgJointLikelihood(relatives, priors)

function Base.show(io::IO, x::_MsgJointLikelihood)
  println(io, )
  printstyled(io, " _MsgJointLikelihood:\n", color=:blue)
  print(io, "  .relatives: ")
  for tp in x.relatives
    print(io, tp.variables, "::", typeof(tp.likelihood).name)
    print(io, "; ")
  end
  println(io)
  print(io, "  .priors: ")
  for k in keys(x.priors)
    print(io, k, ", ")
  end
  println(io)
end

Base.show(io::IO, ::MIME"text/plain", x::_MsgJointLikelihood) = show(io, x)

LikelihoodMessage(; sender::NamedTuple{(:id,:step),Tuple{Int,Int}}=(;id=0, step=0),
                    status::CliqStatus=NULL,
                    beliefDict::Dict=Dict{Symbol, TreeBelief}(),
                    variableOrder::Vector{Symbol}=Symbol[],
                    cliqueLikelihood=nothing,
                    msgType::T=NonparametricMessage(),
                    hasPriors::Bool=true,
                    childSolvDims::Dict{Int, Float64}=Dict{Int, Float64}(),
                    jointmsg::_MsgJointLikelihood=_MsgJointLikelihood(),
                  ) where {T <: MessageType} =
        LikelihoodMessage{T}(sender, status, beliefDict, variableOrder, cliqueLikelihood, msgType, hasPriors, childSolvDims, jointmsg)
#

function Base.show(io::IO, msg::LikelihoodMessage)
  t = typeof(msg)
  fields = fieldnames(t)
  nf = nfields(msg)
  println(io, "LikelihoodMessage:")
  for f in fields
    printstyled(io, f,": ", color=:blue)
    show(io, getproperty(msg, f))
    println(io)
  end

end

Base.show(io::IO, ::MIME"text/plain", msg::LikelihoodMessage) = show(io, msg)

function compare( l1::LikelihoodMessage,
                  l2::LikelihoodMessage;
                  skip::Vector{Symbol}=Symbol[] )
  #
  TP = true
  TP = TP && l1.status == l2.status
  TP = TP && l1.variableOrder == l2.variableOrder
  TP = TP && l1.msgType == l2.msgType
  TP = TP && l1.cliqueLikelihood |> typeof == l2.cliqueLikelihood |> typeof
  for (k,v) in l1.belief
    TP = TP && haskey(l2.belief, k)
    TP = TP && compare(v, l2.belief[k])
  end
  return TP
end

# overload
==(l1::LikelihoodMessage,l2::LikelihoodMessage) = compare(l1,l2)

function BayesTreeNodeData(;status::CliqStatus=NULL,
                            frontalIDs=Symbol[],
                            separatorIDs=Symbol[],
                            inmsgIDs=Symbol[],
                            potIDs=Symbol[],
                            potentials=Symbol[],
                            partialpotential=Bool[],
                            dwnPotentials=Symbol[],
                            dwnPartialPotential=Bool[],
                            cliqAssocMat=Array{Bool}(undef, 0,0),
                            cliqMsgMat=Array{Bool}(undef, 0,0),
                            directvarIDs=Int[],
                            directFrtlMsgIDs=Int[],
                            msgskipIDs=Int[],
                            itervarIDs=Int[],
                            directPriorMsgIDs=Int[],
                            debug=nothing,
                            debugDwn=nothing,
                            allmarginalized=false,
                            initialized=:NULL,
                            upsolved=false,
                            downsolved=false,
                            isCliqReused=false,
                            messages = MessageBuffer()
                          )
  btnd = BayesTreeNodeData(status,
                        frontalIDs,
                        separatorIDs,
                        inmsgIDs,
                        potIDs,
                        potentials,
                        partialpotential,
                        dwnPotentials,
                        dwnPartialPotential,
                        cliqAssocMat,
                        cliqMsgMat,
                        directvarIDs,
                        directFrtlMsgIDs,
                        msgskipIDs,
                        itervarIDs,
                        directPriorMsgIDs,
                        debug,
                        debugDwn,
                        allmarginalized,
                        initialized,
                        upsolved,
                        downsolved,
                        isCliqReused,
                        messages  )
  #
  return btnd
end
#

function compare( c1::BayesTreeNodeData,
                  c2::BayesTreeNodeData;
                  skip::Vector{Symbol}=[] )
  #
  TP = true

  TP = TP && c1.frontalIDs == c2.frontalIDs
  TP = TP && c1.separatorIDs == c2.separatorIDs
  TP = TP && c1.inmsgIDs == c2.inmsgIDs
  TP = TP && c1.potIDs == c2.potIDs
  TP = TP && c1.potentials == c2.potentials
  TP = TP && c1.partialpotential == c2.partialpotential
  TP = TP && c1.dwnPotentials == c2.dwnPotentials
  TP = TP && c1.dwnPartialPotential == c2.dwnPartialPotential
  TP = TP && c1.cliqAssocMat == c2.cliqAssocMat
  TP = TP && c1.cliqMsgMat == c2.cliqMsgMat
  TP = TP && c1.directvarIDs == c2.directvarIDs
  TP = TP && c1.directFrtlMsgIDs == c2.directFrtlMsgIDs
  TP = TP && c1.msgskipIDs == c2.msgskipIDs
  TP = TP && c1.itervarIDs == c2.itervarIDs
  TP = TP && c1.directPriorMsgIDs == c2.directPriorMsgIDs
  TP = TP && c1.debug == c2.debug
  TP = TP && c1.debugDwn == c2.debugDwn
  TP = TP && c1.allmarginalized == c2.allmarginalized
  TP = TP && c1.initialized == c2.initialized
  TP = TP && c1.upsolved == c2.upsolved
  TP = TP && c1.downsolved == c2.downsolved
  TP = TP && c1.isCliqReused == c2.isCliqReused

  return TP
end

function convert(::Type{PackedBayesTreeNodeData}, btnd::BayesTreeNodeData)
  return PackedBayesTreeNodeData(
    btnd.frontalIDs,
    btnd.separatorIDs,
    btnd.inmsgIDs,
    btnd.potIDs,
    btnd.potentials,
    btnd.partialpotential,
    btnd.dwnPotentials,
    btnd.dwnPartialPotential,
    btnd.cliqAssocMat,
    btnd.cliqMsgMat,
    btnd.directvarIDs,
    btnd.directFrtlMsgIDs,
    btnd.msgskipIDs,
    btnd.itervarIDs,
    btnd.directPriorMsgIDs  )
end

function convert(::Type{BayesTreeNodeData}, pbtnd::PackedBayesTreeNodeData)
  btnd = BayesTreeNodeData()
    btnd.frontalIDs = pbtnd.frontalIDs
    btnd.separatorIDs = pbtnd.separatorIDs
    btnd.inmsgIDs = pbtnd.inmsgIDs
    btnd.potIDs = pbtnd.potIDs
    btnd.potentials = pbtnd.potentials
    btnd.partialpotential = pbtnd.partialpotential
    btnd.dwnPotentials = pbtnd.dwnPotentials
    btnd.dwnPartialPotential = pbtnd.dwnPartialPotential
    btnd.cliqAssocMat = pbtnd.cliqAssocMat
    btnd.cliqMsgMat = pbtnd.cliqMsgMat
    btnd.directvarIDs = pbtnd.directvarIDs
    btnd.directFrtlMsgIDs = pbtnd.directFrtlMsgIDs
    btnd.msgskipIDs = pbtnd.msgskipIDs
    btnd.itervarIDs = pbtnd.itervarIDs
    btnd.directPriorMsgIDs = pbtnd.directPriorMsgIDs
  return btnd
end

##==============================================================================
## Cliques
## TreeClique
##==============================================================================

Base.getindex(cId::CliqueId) = cId.value
Base.show(io::IO, ::MIME"text/plain", x::CliqueId) = print(io, x.value)
Base.show(io::IO, x::CliqueId) = print(io, x.value)

getId(c::TreeClique) = c.id

function Base.getproperty(x::TreeClique,f::Symbol)
  if f == :index
    Base.depwarn("`TreeCliqe` field `index` is deprecated, use `id`", :getproperty)
    f = :id
  end
  getfield(x,f)
end

function Base.setproperty!(x::TreeClique, f::Symbol, val)
  if f == :index
    Base.depwarn("`TreeCliqe` field `index` is deprecated, use `id`", :setproperty!)
    f = :id
  end
  return setfield!(x, f, convert(fieldtype(typeof(x), f), val))
end

TreeClique(i::Int) = TreeClique(CliqueId(i), BayesTreeNodeData(), Dict{String,Any}())
TreeClique(id::CliqueId) = TreeClique(id, BayesTreeNodeData(), Dict{String,Any}())

DFG.getLabel(cliq::TreeClique) = cliq.attributes["label"]
function setLabel!(cliq::TreeClique, lbl::String)
  cliq.attributes["label"] = lbl
  lbl
end


text:
using CoordinateTransformations, Rotations, Combinatorics

_getproperty(x::Track, ::Val{:homing}) = x.coords[1:x.tp]
_getproperty(x::Track, ::Val{:searching}) = x.coords[x.tp:end]
_getproperty(x::Track, ::Val{:searchcenter}) = mean(x.searching)
_getproperty(x::Track, ::Val{:turningpoint}) = x.coords[x.tp]
_getproperty(x::Track, ::Val{T}) where {T} = getfield(x, T)
Base.getproperty(x::Track, sym::Symbol) = _getproperty(x, Val(sym))

_getproperty(x::Common, ::Val{:homing}) = x.track.homing
_getproperty(x::Common, ::Val{:searching}) = x.track.searching
_getproperty(x::Common, ::Val{:searchcenter}) = x.track.searchcenter
_getproperty(x::Common, ::Val{:turningpoint}) = x.track.turningpoint
_getproperty(x::Common, ::Val{:pellet1}) = first(x.pellet.xy)
_getproperty(x::Common, ::Val{T}) where {T} = getfield(x, T)
Base.getproperty(x::Common, sym::Symbol) = _getproperty(x, Val(sym))

function get_rotation(nest, feeder)
    v = feeder - nest
    α = -atan(v[2], v[1]) - π/2
    rot = LinearMap(Angle2d(α))
    rot, α
end

(p::LinearMap{Angle2d{Float64}})(x::Missing) = missing

# LinearMap{Angle2d{Float64}}(x::Missing) = missing
function rotate!(x::Common)
    rot, α = get_rotation(x.nest, x.feeder)
    x.nest = rot(x.nest)
    x.feeder = rot(x.feeder)
    x.track.coords .= rot.(x.track.coords)
    x.track.direction .-= α
    x.pellet.xy .= rot.(x.pellet.xy)
    x.originalnest = rot(x.originalnest)
end

function center2!(x::Common, c)
    x.nest -= c
    x.feeder -= c
    for i in eachindex(x.track.coords)
        x.track.coords[i] -= c
    end
    for i in eachindex(x.pellet)
        x.pellet[i].xy -= c
    end
    if !ismissing(x.originalnest)
        x.originalnest -= c
    end
end

function plottrack(template_tex, r)

    p = @pgf Axis(
                  {smooth, legend_style =
                   {
                    at = Coordinate(1.15, 0.5),
                    anchor = "west",
                    legend_columns = 1,
                   },
                   xmajorgrids,
                   ymajorgrids,
                   "ytick distance=25",
                   "xtick distance=25",
                   axis_equal,
                   scale_only_axis,
                   xlabel = "X (cm)",
                   ylabel = "Y (cm)",
                  },
                  PlotInc({only_marks}, Coordinates([Tuple(r.feeder)])),
                  LegendEntry("feeder"),
                  PlotInc({no_markers}, Coordinates(Tuple.(r.homing))),
                  LegendEntry("homing"),
                  PlotInc({only_marks}, Coordinates([Tuple(r.turningpoint)])),
                  LegendEntry("turning point"),
                  PlotInc({no_markers}, Coordinates(Tuple.(r.searching))),
                  LegendEntry("searching"),
                  PlotInc({only_marks}, Coordinates([Tuple(r.searchcenter)])),
                  LegendEntry("search center")
                 );

    if !isempty(r.pellet)
        @pgf push!(p, PlotInc({only_marks}, Coordinates(Tuple.(r.pellet.xy))), LegendEntry("pellet"))
    end
    @pgf push!(p, PlotInc({only_marks}, Coordinates([Tuple(r.nest)])), LegendEntry("nest"));
    if !ismissing(r.originalnest)
        @pgf push!(p, PlotInc({only_marks}, Coordinates([Tuple(r.originalnest)])), LegendEntry("original nest"))
    end
    name = tempname()*".tikz"
    pgfsave(name, p)
    replace(template_tex, "{{track}}" => name)
end

distance(x::Common, from::Symbol, to::Symbol) = round(Int, norm(getproperty(x, from) .- getproperty(x, to)))

getthings(x::Common{Missing}) = isempty(x.pellet) ? (:feeder, :turningpoint, :searchcenter, :nest) : (:feeder, :turningpoint, :searchcenter, :nest, :pellet1)
getthings(x::Common{T}) where {T} = isempty(x.pellet) ? (:feeder, :turningpoint, :searchcenter, :nest, :originalnest) : (:feeder, :turningpoint, :searchcenter, :nest, :originalnest, :pellet1)

function getdistancetable(x::Common)
    things = getthings(x)
    rows = Vector{Vector{String}}()
    for (from, to) in combinations(things, 2)
        push!(rows, string.([from, to, distance(x, from, to)]))
    end
    distancestable = latex_tabular(String, Tabular("llr"),
                                   [Rule(:top),
                                    ["From", "To", "cm"],
                                    Rule(:mid),
                                    rows...,
                                    Rule(:bottom)])
end

function plotdirection(template_tex, r)
    t = time(r)
    tl = range(0, stop = t[end], length = 5length(t))
    d = rad2deg.(r.track.direction .- π)
    itp = interpolate((t, ), d, Gridded(Linear()))
    dl = itp.(tl)
    tpl = findfirst(x -> x ≥ t[r.track.tp], tl)
    p = @pgf PolarAxis(
                       {legend_style =
                        {
                         at = Coordinate(1.15, 0.5),
                         anchor = "west",
                         legend_columns = 1,
                        },
                        # axis_equal,
                        scale_only_axis,
                        xlabel = "Direction (°)",
                        ylabel = "Time (min.)",
                       },
                       PlotInc({no_markers}, Coordinates(dl[1:tpl], tl[1:tpl])),
                       LegendEntry("homing"),
                       PlotInc({only_marks}, Coordinates([dl[tpl]], [tl[tpl]])),
                       LegendEntry("turning point"),
                       PlotInc({no_markers}, Coordinates(dl[tpl:end], tl[tpl:end])),
                       LegendEntry("searching")
                      )
    name = tempname()*".tikz"
    pgfsave(name, p)
    replace(template_tex, "{{direction}}" => name)
end

function plotpathlength(template_tex, r)
    t = time(r)
    pl = cumsum(r.track.distance)/100
    p = @pgf Axis(
                  {legend_style =
                   {
                    at = Coordinate(1.15, 0.5),
                    anchor = "west",
                    legend_columns = 1,
                   },
                   # axis_equal,
                   scale_only_axis,
                   xlabel = "Time (min.)",
                   ylabel = "Path length (m)",
                  },
                  PlotInc({no_markers}, Coordinates(t[1:r.track.tp], pl[1:r.track.tp])),
                  LegendEntry("homing"),
                  PlotInc({only_marks}, Coordinates([t[r.track.tp]], [pl[r.track.tp]])),
                  LegendEntry("turning point"),
                  PlotInc({no_markers}, Coordinates(t[r.track.tp:end], pl[r.track.tp:end])),
                  LegendEntry("searching")
                 )
    name = tempname()*".tikz"
    pgfsave(name, p)
    replace(template_tex, "{{pathlength}}" => name)
end

function plotspeed(template_tex, r)

    t = time(r)
    s = r.track.distance/r.track.Δt
    p = @pgf Axis(
                  {
                   legend_style =
                   {
                    at = Coordinate(1.15, 0.5),
                    anchor = "west",
                    legend_columns = 1,
                   },
                   # axis_equal,
                   scale_only_axis,
                   xlabel = "Time (min.)",
                   ylabel = "Speed (cm/s)",
                  },
                  PlotInc({no_markers}, Coordinates(t[1:r.track.tp], s[1:r.track.tp])),
                  LegendEntry("homing"),
                  PlotInc({only_marks}, Coordinates([t[r.track.tp]], [s[r.track.tp]])),
                  LegendEntry("turning point"),
                  PlotInc({no_markers}, Coordinates(t[r.track.tp:end], s[r.track.tp:end])),
                  LegendEntry("searching")
                 )

    name = tempname()*".tikz"
    pgfsave(name, p)
    template_tex = replace(template_tex, "{{speed}}" => name)

    p2 = @pgf Axis(
                  {colormap_name = "viridis", smooth, legend_style =
                   {
                    at = Coordinate(1.15, 0.5),
                    anchor = "west",
                    legend_columns = 1,
                   },
                   xmajorgrids,
                   ymajorgrids,
                   "ytick distance=25",
                   "xtick distance=25",
                   axis_equal,
                   scale_only_axis,
                   xlabel = "X (cm)",
                   ylabel = "Y (cm)",
                   "colorbar horizontal",
                   colorbar_style={
                                   xticklabel_pos="right", xlabel = "Speed (cm/s)", at=Coordinate(0,1.1), anchor="north west"
                                  },
                  },
                  Plot(
                       {
                        mesh,
                        # shader="faceted interp",
                        point_meta = raw"\thisrow{s}",
                        line_width = 3,
                       },
                       Table(x = first.(r.track.coords), y = last.(r.track.coords), s = s)
                      ),
                  LegendEntry("track"),
                  PlotInc({only_marks}, Coordinates([Tuple(r.feeder)])),
                  LegendEntry("feeder"),
                  PlotInc({only_marks}, Coordinates([Tuple(r.turningpoint)])),
                  LegendEntry("turning point"),
                  PlotInc({only_marks}, Coordinates([Tuple(r.searchcenter)])),
                  LegendEntry("search center")
                 );
    if !isempty(r.pellet)
        @pgf push!(p2, PlotInc({only_marks}, Coordinates(Tuple.(r.pellet.xy))), LegendEntry("pellet"))
    end
    @pgf push!(p2, PlotInc({only_marks}, Coordinates([Tuple(r.nest)])), LegendEntry("nest"));
    if !ismissing(r.originalnest)
        @pgf push!(p2, PlotInc({only_marks}, Coordinates([Tuple(r.originalnest)])), LegendEntry("original nest"))
    end

    name = tempname()*".tikz"
    pgfsave(name, p2)
    replace(template_tex, "{{trackspeed}}" => name)

end

function plotangularchange(template_tex, r)
    t = time(r)
    tl = range(0, stop = t[end-1], length = 1000)
    d = diff(rad2deg.(r.track.direction))
    itp = interpolate((t[1:end-1], ), d, Gridded(Linear()))
    dl = itp.(tl)
    tpl = findfirst(x -> x ≥ t[r.track.tp], tl)
    if isnothing(tpl)
        tpl = length(tl)
    end
    p = @pgf PolarAxis(
                       {legend_style =
                        {
                         at = Coordinate(1.15, 0.5),
                         anchor = "west",
                         legend_columns = 1,
                        },
                        # axis_equal,
                        scale_only_axis,
                        xlabel = "Direction change (°)",
                        ylabel = "Time (min.)",
                        xmin = -90,
                        xmax = 90,
                       },
                       PlotInc({no_markers}, Coordinates(dl[1:tpl], tl[1:tpl])),
                       LegendEntry("homing"),
                       PlotInc({only_marks}, Coordinates([dl[tpl]], [tl[tpl]])),
                       LegendEntry("turning point"),
                       PlotInc({no_markers}, Coordinates(dl[tpl:end], tl[tpl:end])),
                       LegendEntry("searching")
                      )
    name = tempname()*".tikz"
    pgfsave(name, p)
    replace(template_tex, "{{angularchange}}" => name)
end

# time(x::Common, s::Symbol) = range(0, step = x.track.Δt/60, length = length(getfield(x.track, s)))
# time(x::Common, length::Int) = range(0, step = x.track.Δt/60, length = length)

time(x::Common) = range(0, step = x.track.Δt/60, length = length(x.track.coords))

function plotresult(r::Common, targetfile)

    template_tex = read("/home/yakir/dungProject/DungAnalyse/src/results.tex", String)

    template_tex = plottrack(template_tex, r)
    template_tex = plotdirection(template_tex, r)
    template_tex = plotspeed(template_tex, r)
    template_tex = plotangularchange(template_tex, r)
    template_tex = plotpathlength(template_tex, r)

    template_tex = replace(template_tex,  "{{distancestable}}" => getdistancetable(r))

    base = tempname()
    texfile = "$base.tex"
    write(texfile, template_tex)
    latex_success, log, latexcmd = PGFPlotsX.run_latex_once(texfile, latexengine(), vcat(PGFPlotsX.DEFAULT_FLAGS, PGFPlotsX.CUSTOM_FLAGS))
    if latex_success
        mv("$base.pdf", targetfile, force=true)
    else
        error("failed to compile latex file $base.tex")
    end
end

# r = last(first(trackdata)).runs[3].data
# rotate!(r)
# center2!(r, r.nest)
# plotresult(r, "a.pdf")


text:
### utilities.jl

#
# Strings
#
const rchar = "rml"
const rmchar = "rmean"
const rqchar = "rquant"
const wchar = "waic"
const routchar = "rout"
const runchar = "run"
const chichar = "chi"
const chioutchar = "chiout"
const rlastchar = "rlast"
const onchar = "LCON"
const offchar = "LCOFF"
const fishchar = "FISH"
const txtstr = ".txt"

"""
eig_decompose(M)

Take matrix M and return values and vectors
"""
function eig_decompose(M)
    Meiv = eigen(M)
    return Meiv.values, Meiv.vectors
end

"""
nonzero_rows(T)

Returns an array of row indices that have at least one nonzero element for matrix T
"""
function nonzero_rows(T)
    n = size(T)[2]
    nonzeros = Array{Int}(undef,0)
    for  a in 1:size(T)[1]
        if T[a,:] != zeros(n)
            push!(nonzeros,a)
        end
    end
    return nonzeros
end

"""
normalize_histogram(hist)

Returns normalized histogram hist
"""
normalize_histogram(hist) = hist/sum(hist)

"""
var_update(vartuple, newValue)

vartuple = count, mean, M2
iteration step for online algorithm to compute
mean and variance
returns count, mean, M2
"""
function var_update(vartuple, newValue)
    count, mean, M2 = vartuple
    count += 1
    delta = newValue - mean
    mean .+= delta / count
    delta2 = newValue - mean
    M2 .+= delta .* delta2
    return count, mean, M2
end
"""
cov_update(covtuple,newValue1,newValue2)

update covtuple for online algorithm to compute
    covariance
"""
function cov_update(covtuple,newValue1,newValue2)
    count, meanx, meany, C = covtuple
    count += 1
    dx = newx - meanx
    meanx .+= dx / count
    meany .+= (newy - meany) / count
    C .+= dx * (newy - meany)
    return count, meanx, meany, C
end
"""
finalize(vartuple)

Retrieve the mean, variance and sample variance from an aggregate
"""
function finalize(vartuple)
    count, mean, M2 = vartuple
    if count < 2
        return mean, M2
    else
        return mean, M2 ./ (count - 1)
    end
end

"""
truncatehistogram(x::Array{Int,1},yield::Float64,nhistmax::Int)

only keep columns of histogram x that accounts for yield fraction of total sum
"""
function truncate_histogram(x::Array,yield::Float64,nhistmax::Int)
    counts = sum(x)
    if counts > 0
        nhist = 1
        ratio = 0.
        # Only keep yield of all data to limit support of histogram
        while ratio < yield && nhist <= nhistmax
            ratio = sum(x[1:nhist])/counts
            nhist += 1
        end
        return x[1:min(nhist,length(x))]
    else
        return 0
    end
end

"""
combine_histogram(x1::Array,x2::Array)

add histograms x1 and x2

"""
function combine_histogram(x1::Array,x2::Array)
    n = min(length(x1),length(x2))
    x1[1:n] + x2[1:n]
end

"""
pooled_mean(means,counts)

compute weighted average of means given count totals
"""
pooled_mean(means,counts) = counts'*means/sum(counts)

"""
pooled_variance(vars,counts)

computed weighted average of variances given counts
"""
pooled_variance(vars,counts) = (counts .- 1)'*vars/sum(counts .- 1)

"""
pooled_std(std,counts)

compute weighted average of standard deviations given counts
"""
pooled_std(std,counts) = sqrt(pooled_variance(std.^2,counts))

"""
var_ratio(mua,mub,vara,varb,cov)

Compute ratio of variances of two variables given means, variances, and covariance of variables
"""
var_ratio(mua,mub,vara,varb,cov) = mua^2/mub^2*(vara/mua^2 - 2*cov/(mua*mub) + varb/mub^2)

"""
decimal(x::Vector)

convert digits of base 3 number to base 10
"""
function decimal(x::Vector,base::Int=3)
    nr = length(x)
    x' * base .^ collect(0:nr-1)
end
"""
findbase(l,n,nr)

Find the number of G states for transition rate matrix of size l
"""
function findbase(l,n,nr)
    if l == (n+1)*3^nr
        return 3
    elseif l == (n+1)*2^nr
        return 2
    else
        throw("wrong length")
    end
end

"""
KL(x,y)

Kullback-Leibler distance
"""
KL(x,y) = sum(x .* (log.(max.(x,eps(Float64))) - log.(max.(y,eps(Float64)))))

"""
trim(h::Array,nh::Array)

"""
function trim(h::Array,nh::Array)
    for i in eachindex(h)
        h[i] = h[i][1:nh[i]]
    end
    return h
end

"""
LogNormal_array(param,cv)
Prior distribution arrays
"""
LogNormal_array(param,cv) = distributionarray(param,cv,LogNormal)

"""
Gamma_array(param,cv)

"""
Gamma_array(param,cv) = distributionarray(param,cv,Gamma)

"""
function LogNormalBeta_array(param,cv,ind)

"""
function LogNormalBeta_array(param,cv,ind)
    if ind == 1
        d = [setBeta(param[ind],cv[ind])]
    else
        barind = 1:ind-1
        d = LogNormalarray(param[barind],cv[barind])
        push!(d,setBeta(param[ind],cv[ind]))
    end
    return d
end

"""
distribution_array(param,cv,dist)

fill an array with dist(param,cv)
"""
function distributionBeta_array(param::Vector,cv::Vector,ind::Int,dist=LogNormal)
    if ind == 1
        d = [Beta_meancv(param[ind],cv[ind])]
    else
        barind = 1:ind-1
        d = distribution_array(param[barind],cv[barind],dist)
        push!(d,Beta_meancv(param[ind],cv[ind]))
    end
    return d
end

"""
distribution_array(param::Vector,cv,dist=LogNormal)

"""
function distribution_array(param::Vector,cv,dist=LogNormal)
    d = []
    for i in eachindex(param)
        if dist == LogNormal
            push!(d,LogNormal_meancv(param[i],cv[i]))
        elseif dist == Gamma
            push!(d,Gamma_meancv(param[i],cv[i]))
        else
            push!(d,dist(param[i],cv[i]))
        end

    end
    return d
end

"""
Gamma_meancv(mean,cv)
LogNormal_meancv(mean,cv)
Beta_meancv(mean,cv)

Reparameterized distributions to mean and coefficient of variation arguments
"""
Gamma_meancv(mean,cv) = Gamma(1/cv^2,mean * cv^2)
LogNormal_meancv(mean,cv) = LogNormal(mulognormal(mean,cv), sigmalognormal(cv))
function Beta_meancv(m,cv)
    cv2 = cv^2
    fac = (1-m)/cv2/m - 1
    if fac <= 0.
        fac = 2/m
    end
    alpha = m*fac
    beta = (1-m)*fac
    Beta(alpha,beta)
end
"""
sigmalognormal(cv)
mulognormal(mean,cv)

compute LogNormal mu and sigma parameters
given desired mean and coefficient of variation
"""
sigmalognormal(cv) = sqrt.(log.(1 .+ cv .^ 2))

mulognormal(mean,cv) = log.(mean) - .5*log.(1 .+ cv .^ 2)

"""
Parameter information for GR models
"""
function num_rates(G,R)
    n = G - 1
    2*n + 2*R + 2
end
function Grange(G)
    n = G - 1
    1 : 2*n
end
function initiation(G)
    n = G - 1
    2*n + 1
end
function Rrange(G,R)
    n = G - 1
    2*n + 2 : 2*n + R + 1
end
function Srange(G,R)
    n = G - 1
    2*n + R + 2 : 2*n + 2*R + 1
end

"""
mean_histogram(x)

"""
mean_histogram(x::Vector{Float64}) = collect(0:length(x)-1)' * x/sum(x)

function mean_histogram(x::Vector{Array})
    y = Vector{Float64}(undef,length(x))
    for i in eachindex(x)
        y[i] = mean_histogram(x[i])
    end
    return y
end

"""
    Difference_Zscore(x1,x2,sig1,sig2) = (x1 .- x2)./ sqrt(sig1.^2 + sig2.^2)

"""
Difference_Zscore(x1,x2,sig1,sig2) = (x1 - x2) / sqrt(sig1^2 + sig2^2)

"""
m2_histogram(x)

"""
m2_histogram(x) = (collect(0:length(x)-1).^2)' * x/sum(x)

"""
var_histogram(x)

"""
function var_histogram(x)
    m2_histogram(x) - mean_histogram(x)^2
end

function moment_histogram(x,n)
    y = collect(1:length(x)) .- 1
    v = (y .- mean_histogram(x)).^n
    v' * x/sum(x)
end

function factorial_moment(h::Vector,n)
    m = 0
    for i in n:length(h)
        a = i-1
        for j in 1:n-1
            a *= i-j-1
        end
        m += h[i]*a
    end
    return m / sum(h)
end

function moment_param_estimates(h)
    e1 = factorial_moment(h,1)
    e2 = factorial_moment(h,2)
    e3 = factorial_moment(h,3)
    r1 = e1
    r2 = e2/e1
    r3 = e3/e2
    println(r1,", ",r2,", ",r3)

    lambda = 2*r1*(r3-r2)/(r1*r2-2*r1*r3+r2*r3)
    mu = 2*(r2-r1)*(r1-r3)*(r3-r2)/(r1*r2-2*r1*r3+r2*r3)/(r1-2*r2+r3)
    nu = (-r1*r2 + 2*r1*r3-r2*r3)/(r1-2*r2+r3)

    return lambda,mu,nu
end
"""
tstat_2sample(x1,x2)

Compute t statistics of histograms x1 and x2
"""
tstat_2sample(x1,x2) = (mean_histogram(x1) - mean_histogram(x2))/(sqrt(var_histogram(x1)/length(x1) + var_histogram(x2)/length(x2)))

"""
log_2sample(x1,x2)

compute the log of the ratio of means of histograms x1 and x2
"""
log_2sample(x1,x2) = log(mean_histogram(x1)) - log(mean_histogram(x2))

"""
delta_2sample(x1,x2)

compute the difference of means of x1 and x2
"""
delta_2sample(x1,x2) = mean_histogram(x1)-mean_histogram(x2)

"""
delta_2frac(x1,x2)

compute E(x1) - E(x2) of means normalize dot mean of x1
"""
delta_2frac(x1,x2) = delta_2sample(x1,x2)/mean_histogram(x1)

"""
mediansmooth(xin,window)

median smooth histogram xin over a window
"""
function mediansmooth(xin,window)
    x = copy(xin)
    halfwin = div(window-1,2)
    for i in 1:length(x)-window+1
        x[i+halfwin] = median(x[i:i+window-1])
    end
    x
end

"""
residenceprob_G(file,G,header)
residenceprob_G(r,n)

Residence probability of G states
given by exact steady state solution
of master equation

"""
function residenceprob_G(file::String,G,header=false)
    r = get_all_rates(file,header)
    m = size(r)[1]
    p = Array{Any,2}(undef,m,G+1)
    n = G-1
    for i in 1:m
        p[i,1] = r[i,1]
        p[i,2:end] = residenceprob_G(r[i,2:2*n+1],n)
    end
    return p
end

function residenceprob_G(r::Vector,n::Int)
    Gss = Array{Float64,2}(undef,1,n+1)
    Gss[1,1] = 1.
    for k in 1:n
        Gss[1,k+1] = Gss[1,k]*r[2*k-1]/r[2*k]
    end
    Gss ./= sum(Gss)
end

"""
splicesiteusage()

splice site usage probability is the probabilty of ejection times
the probability of not ejection prior to that point
"""
splicesiteusage(model::GRSMmodel) = splicesiteusage(model.rates,model.G-1,model.R)
function splicesiteusage(r::Vector,n::Int,nr::Int)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
	ssf = zeros(nr)
	survival = 1
	for j in 1:nr
		ssf[j] = eta[j]/(nu[j+1]+eta[j])*survival
		survival *= nu[j+1]/(nu[j+1]+eta[j])
	end
	return ssf
end

"""
burstoccupancy(n,nr,r)
Burst size distribution of GRS  model
for total pre-mRNA occupancy and
unspliced (visible) occupancy
obtained by marginalizing over conditional steady state distribution
"""
burstoccupancy(model::GRSMmodel) = burstsize(model.G-1,model.R,model.rates)

function burstoccupancy(n::Int,nr::Int,r::Vector)
    T =  mat_GSR_T(r,n,nr)
    pss = normalized_nullspace(T)
	Rss = zeros(nr)
	Rssvisible = zeros(nr)
	ssf = zeros(nr)
	asum = 0
	for w=1:nr
		for i = 1:n+1, z = 1:3^nr
			zdigits = digits(z-1,base=3,pad=nr)
			a = i + (n+1)*(z-1)
			if sum(zdigits .== 2) == w
				Rssvisible[w] += pss[a]
			end
			if sum(zdigits .> 0) == w
				Rss[w] += pss[a]
			end
		end
	end
	Rss ./= sum(Rss), Rssvisible ./= sum(Rssvisible)
end

model2_mean(ron,roff,eject,decay,nalleles) = 2*model2_mean(ron,roff,eject,decay)

model2_mean(ron,roff,eject,decay) = ron/(ron+roff)*eject/decay

model2_variance(ron,roff,eject,decay,nalleles) = 2*model2_variance(ron,roff,eject,decay)

model2_variance(ron,roff,eject,decay) = ron/(ron+roff)*eject/decay + ron*roff/(ron+roff)^2 * eject^2/decay/(ron + roff + decay)

"""
test_steadystatemodel(model,nhist)

Compare chemical master solution to Gillespie simulation for steadystate mRNA distribution

"""
function test_steadystatemodel(model::AbstractGMmodel,nhist)
    G = model.G
    r = model.rates
    g1 = steady_state(r[1:2*G],G-1,nhist,model.nalleles)
    g2 = simulatorGM(r[1:2*G],G-1,nhist,model.nalleles)
    return g1,g2
end
#
# function fit_rna_test(root)
#     gene = "CENPL"
#     cell = "HCT116"
#     fish = false
#     data = data_rna(gene,"MOCK","data/HCT116_testdata",fish,"scRNA_test",root)
#     model = model_rna(gene,cell,2,fish,.01,[1,2,3],(),"scRNA_test","scRNA_test",1,".",data,.05,1.0)
#     options = MHOptions(10000,2000,0,120.,1.,100.)
#     fit,stats,waic = run_mh(data,model,options,1);
#     return stats.meanparam, fit.llml
# end
#

# function online_covariance(data1, data2)
#     meanx = meany = C = n = 0
#     for x in data1, y in data2
#         n += 1
#         dx = x - meanx
#         meanx += dx / n
#         meany += (y - meany) / n
#         C += dx * (y - meany)
#         population_covar = C / n
#         # Bessel's correction for sample variance
#         sample_covar = C / (n - 1)
#     end
# end


text:
using CSV, Primes, DataFrames

kappa=1.4;
Q=[1 2 3];

for a=1:2000
    for b=a+1:100000
        c=a+b;
        X=factor(Vector,a*b*c);
        L=length(X);
        rad=X[1];
        if [gcd(a,b) gcd(a,c) gcd(b,c)]==[1 1 1]
            for i=2:L
                if X[i]==X[i-1]
                    rad=rad*1;
                else
                    rad=rad*X[i];
                end
            end
            if (c >= rad^kappa)
                q=[a b c];
                Q=[Q;q];
            end
        end
    end
end

R=DataFrame(Q,:auto);
CSV.write("abc14.csv", R);


text:
"""
    BernoulliLikelihood

Bernoulli likelihood is to be used if we assume that the
uncertainity associated with the data follows a Bernoulli distribution.

```math
    p(y|f) = Bernoulli(y | f)
```
On calling, this would return a Bernoulli distribution with `f` probability of `true`.
"""
struct BernoulliLikelihood end

(l::BernoulliLikelihood)(f::Real) = Bernoulli(logistic(f))

(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(Bernoulli.(logistic.(fs)))


text:
# Feature Selection
# =================

"""
A list of features.

Fields
------

- `names`: feature names.
- `selected`: a bit vector of selected features.

`names` and `selected` must have the same number of elements.
"""
struct Features
    names::Vector{String}
    selected::BitVector

    function Features(names, selected)
        if length(names) != length(selected)
            throw(ArgumentError("mismatching length"))
        end
        return new(names, selected)
    end
end

"""
    Features(names)

Create a features object from feature names `names`.

All features are not selected by default. To select and unselect some features,
use `addfeatures!` and `dropfeatures!`, respectively.

`selectfeatures` algorithmically selects important features from a count matrix.
"""
Features(names::AbstractVector{String}) = Features(names, falses(length(names)))

Base.copy(features::Features) = Features(copy(features.names), copy(features.selected))

function Base.show(io::IO, features::Features)
    print(io, summary(features), "(<#features=$(nfeatures(features)),#selected=$(nselected(features))>)")
end

nfeatures(features::Features) = length(features.names)
nselected(features::Features) = sum(features.selected)
selectedfeatures(features::Features) = features.names[features.selected]

"""
    addfeatures(features, names) -> Features

Create a new features object by adding `names` to `features`.

See also `dropfeatures` and `addfeatures!`.
"""
addfeatures(features::Features, featurenames::AbstractVector{String}) = addfeatures!(copy(features), featurenames)

"""
    addfeatures!(features, names) -> Features

Add `names` to `features` in place.

See also `dropfeatures!` and `addfeatures`.
"""
function addfeatures!(features::Features, featurenames::AbstractVector{String})
    for name in featurenames
        i = findfirst(isequal(name), features.names)
        if i == nothing
            throw(ArgumentError("not found feature '$(name)'"))
        end
        features.selected[i] = true
    end
    return features
end

"""
    dropfeatures(features, names) -> Features

Create a new features object by dropping `names` from `features`.

See also `addfeatures` and `dropfeatures!`.
"""
dropfeatures(features::Features, featurenames::AbstractVector{String}) = dropfeatures!(copy(features), featurenames)

"""
    dropfeatures!(features, names) -> Features

Drop `names` from `features` in place.

See also `addfeatures!` and `dropfeatures`.
"""
function dropfeatures!(features::Features, featurenames::AbstractVector{String})
    for name in featurenames
        i = findfirst(isequal(name), features.names)
        if i == nothing
            throw(ArgumentError("not found feature '$(name)'"))
        end
        features.selected[i] = false
    end
    return features
end

"""
    selectfeatures(counts, featurenames; n_min_features=cld(size(counts, 1), 10)) -> Features

Select features from `counts`.

Arguments
---------

- `counts`: transcriptome expression matrix (features x cells).
- `featurenames`: vector of feature names.
- `n_min_features`: number of minimum features [default: 10% of all features].
"""
function selectfeatures(counts::AbstractMatrix,
                        featurenames::AbstractVector{String};
                        n_min_features::Real=cld(size(counts, 1), 10))
    M = size(counts, 1)
    if M != length(featurenames)
        throw(ArgumentError("mismatching featurenames size"))
    elseif !(1 ≤ n_min_features ≤ M)
        throw(ArgumentError("invalid n_min_features"))
    end
    mincounts, maxcounts = feature_extrema(counts)
    lowerbound = sort(maxcounts, rev=true)[n_min_features]
    return Features(featurenames, (maxcounts .≥ lowerbound) .& (maxcounts .> mincounts))
end

function feature_extrema(counts::AbstractMatrix)
    return vec(minimum(counts, dims=2)), vec(maximum(counts, dims=2))
end

# https://github.com/JuliaLang/julia/issues/27836
function feature_extrema(counts::SparseMatrixCSC{T}) where T <: Real
    m, n = size(counts)
    mincounts = fill(typemax(T), m)
    maxcounts = fill(typemin(T), m)
    for col in 1:n
        for j in counts.colptr[col]:counts.colptr[col+1]-1
            i = counts.rowval[j]
            val = counts.nzval[j]
            mincounts[i] = min(val, mincounts[i])
            maxcounts[i] = max(val, maxcounts[i])
        end
    end
    return mincounts, maxcounts
end


text:
function _start_tracking_file(filename::AbstractString)
    _realpath::String = convert(String, realpath(filename))::String
    return _start_tracking_file(_realpath)
end

function _start_tracking_file(filename::String)::Nothing
    _realpath::String = convert(String, realpath(filename))::String
    original_block = _read_julia_file(_realpath)
    always_assert( length(original_block.args) >= 1,
                  "length(original_block.args) >= 1")
    always_assert( original_block.args[1] isa LineNumberNode,
                  "original_block.args[1] isa LineNumberNode")
    instrumented_block = instrument(_realpath, original_block)
    backup_filename = string(_realpath, ".backup")
    rm(backup_filename; force = true, recursive = true)
    @debug("Moving file from src to dst",
           src = _realpath,
           dst = backup_filename)
    mv(_realpath, backup_filename; force = true)
    _write_julia_file(_realpath, instrumented_block)
    return nothing
end


text:
include("cg_sem.jl")
using StaticArrays: SVector

wave_speed(ρ, p, γ) = sqrt(γ * p ./ ρ)

Base.@kwdef struct Params{T}
    T_virt_surf::T = 280
    T_min_ref::T = 230
    MSLP::T = 10^5
    grav::T = 98 // 100
    R_d::T = 287
    γ::T = 7 // 5
end

# Taken from Daniel's code
function decaying_temperature_profile(params::Params, z)
    # Scale height for surface temperature
    H_sfc = params.R_d * params.T_virt_surf / params.grav
    H_t = H_sfc

    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = params.T_virt_surf - params.T_min_ref
    Tv = params.T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / params.T_virt_surf
    p = -H_t * (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
    p /= H_sfc * (1 - ΔTv′^2)
    p = params.MSLP * exp(p)
    ρ = p / (params.R_d * Tv)
    return (p = p, ρ = ρ)
end

function element_split_form!(∂q, q, O, H, z)
  Nq, K, nflds = size(q)

  Q = O.Q
  @inbounds for e in 1:K
    for i in 1:Nq, j in 1:Nq
      # get the two states
      qi = @view q[i, e, :]
      qj = @view q[j, e, :]
      zi = z[i, e]
      zj = z[j, e]
      # Evaluate the flux difference form:
      #   ((H ∘ Qᵀ) - (Q ∘ H)) 𝟙
      ∂q[i, e, :] += (Q[i, j] - Q[j, i]) * H(qi, qj, zi, zj)
    end
  end
end

pressure(ρ, w, ρe, z, g, γ) = (γ - 1) * (ρe - ρ * w^2 / 2 - ρ * g * z)

function KennedyGruber(q_i, q_j, z_i, z_j, grav, pressure)
  @inbounds (ρ_i, w_i, ρe_i) = q_i[1:3]
  @inbounds (ρ_j, w_j, ρe_j) = q_j[1:3]
  e_i, e_j = ρe_i / ρ_i, ρe_j / ρ_j
  p_i = pressure(ρ_i, w_i, ρe_i, z_i, grav)
  p_j = pressure(ρ_j, w_j, ρe_j, z_j, grav)
  ϕ_i = z_i * grav
  ϕ_j = z_j * grav
  (fρ, fρw, fρe) = KennedyGruber(ρ_i, w_i, e_i, p_i, ϕ_i, ρ_j, w_j, e_j, p_j, ϕ_j)
  return SVector(fρ, fρw / ρ_i, fρe)
end

function KennedyGruber(ρ_i, w_i, e_i, p_i, ϕ_i, ρ_j, w_j, e_j, p_j, ϕ_j)
  ρ̄ = (ρ_i + ρ_j) / 2
  w̄ = (w_i + w_j) / 2
  Δw = (w_j - w_i) / 2
  ē = (e_i + e_j) / 2
  Δϕ = (ϕ_j - ϕ_i) / 2

  p̄ = (p_i + p_j) / 2

  fρ  = ρ̄ * w̄
  fρw = ρ̄ * w̄ * Δw + p̄ + ρ̄ * Δϕ
  fρe = (ρ̄ * ē + p̄) * w̄
  return (fρ, fρw, fρe)
end

function tendency!(
    ∂q,
    q,
    t,
    H,
    mesh,
    elem,
    grav,
    p_,
    bc = (nothing, nothing),
    forcing = nothing,
)
    fill!(∂q, 0)

    # Evaluate the volume terms
    element_split_form!(∂q, q, elem, H, mesh.zdg)

    # XXX: Hack in the bcs
    ρ = q[1, 1, 1]
    w = q[1, 1, 2]
    ρe = q[1, 1, 3]
    p = p_(ρ, w, ρe, mesh.zdg[1], grav)
    ∂q[1, 1, 2] -= p / ρ
    ρ = q[end, end, 1]
    w = q[end, end, 2]
    ρe = q[end, end, 3]
    p = p_(ρ, w, ρe, mesh.zdg[end], grav)
    ∂q[end, end, 2] += p / ρ

    # Add in MMS formcing
    if !isnothing(forcing)
        error("not implemented yet")
    end

    # Use upwind boundary treatment
    if !isnothing(bc[1])
        error("not implemented yet")
    end

    if !isnothing(bc[2])
        error("not implemented yet")
    end

    # DG -> CG projection and then scatter for storage
    dg2cg_scatter(∂q, mesh)
end


text:
println("---- Benchmarking ----")
include("probability_distribution_calculations.jl")


text:
to_int(s) = parse(Int, s)

str2tuple(s) = (x -> x + 1).(to_int.(split(s, ',')))

parse_line(line) = str2tuple.(split(line, " -> "))

function line_to_points((x1, y1), (x2, y2))
    if x1 == x2
        [(y, x1) for y = min(y1, y2):max(y1, y2)]
    elseif y1 == y2
        [(y1, x) for x = min(x1, x2):max(x1, x2)]
    else
        Δx = x2 > x1 ? 1 : -1
        Δy = y2 > y1 ? 1 : -1
        [(y1 + Δy * n, x1 + Δx * n) for n = 0:abs(x1 - x2)]
    end
end

function solve2(filename)
    board = zeros(1000, 1000)
    for line in readlines(filename)
        points = line_to_points(parse_line(line)...)
        for point in points
            board[point...] += 1
        end
    end
    length(filter(>=(2), board))
end

solve2("5.txt")


text:
"""
Generates various SearchLight files.
"""
module Generator

import Revise
import Unicode, Logging
using SearchLight, SearchLight.FileTemplates, SearchLight.Inflector, SearchLight.Configuration, SearchLight.Migration

"""
    newmodel(cmd_args::Dict{String,Any}) :: Nothing

Generates a new SearchLight model file and persists it to the resources folder.
"""
function newmodel(cmd_args::Dict{String,Any}; pluralize::Bool = true) :: Nothing
  resource_name = uppercasefirst(cmd_args["model:new"])

  SearchLight.Inflector.is_singular(resource_name) &&
    (resource_name = SearchLight.Inflector.to_plural(resource_name))

  resource_path = setup_resource_path(resource_name)
  mfn = model_file_name(resource_name)
  write_resource_file(resource_path, mfn, resource_name, :model, pluralize = pluralize) &&
    @info "New model created at $(abspath(joinpath(resource_path, mfn)))"

  nothing
end
function newmodel(resource_name::Union{String,Symbol}) :: Nothing
  newresource(resource_name)
end

"""
    newresource(resource_name::Union{String,Symbol}) :: Nothing

Generates all the files associated with a new resource and persists them to the resources folder.
"""
function newresource(resource_name::Union{String,Symbol}; pluralize::Bool = true) :: Nothing
  resource_name = string(resource_name)

  sf = SearchLight.Inflector.tosingular(resource_name)

  model_name = (sf === nothing ? resource_name : sf) |> uppercasefirst
  newmodel(Dict{String,Any}("model:new" => model_name), pluralize = pluralize)
  new_table_migration(Dict{String,Any}("migration:new" => resource_name), pluralize = pluralize)

  resource_name = uppercasefirst(resource_name)
  pluralize && SearchLight.Inflector.is_singular(resource_name) &&
    (resource_name = SearchLight.Inflector.to_plural(resource_name))

  resource_path = setup_resource_path(resource_name)
  for (resource_file, resource_type) in [(validator_file_name(resource_name), :validator)]
    write_resource_file(resource_path, resource_file, resource_name, resource_type, pluralize = pluralize) &&
      @info "New $resource_type created at $(abspath(joinpath(resource_path, resource_file)))"
  end

  isdir(SearchLight.TEST_PATH_UNIT) || mkpath(SearchLight.TEST_PATH_UNIT)
  test_file = resource_name * SearchLight.TEST_FILE_IDENTIFIER |> lowercase
  write_resource_file(SearchLight.TEST_PATH_UNIT, test_file, resource_name, :test, pluralize = pluralize) &&
    @info "New unit test created at $(abspath(joinpath(SearchLight.TEST_PATH_UNIT, test_file)))"

  SearchLight.load_resources()

  nothing
end

"""
"""
function new_table_migration(cmd_args::Dict{String,Any}; pluralize::Bool = true) :: Nothing
  resource_name = uppercasefirst(cmd_args["migration:new"])

  SearchLight.Inflector.is_singular(resource_name) && pluralize &&
    (resource_name = SearchLight.Inflector.to_plural(resource_name))

  migration_name = "create_table_" * lowercase(resource_name)
  SearchLight.Migration.new_table(migration_name, lowercase(resource_name))

  nothing
end
function new_table_migration(migration_name::String) :: Nothing
  new_table_migration(Dict{String,Any}("migration:new" => migration_name))
end

function new_migration(cmd_args::Dict{String,Any}) :: Nothing
  migration_name = replace(uppercasefirst(cmd_args["migration:new"]) |> lowercase, " "=>"_")

  Migration.new(migration_name)

  nothing
end
function new_migration(migration_name::String) :: Nothing
  new_migration(Dict{String,Any}("migration:new" => migration_name))
end

"""
    setup_resource_path(resource_name::String) :: String

Computes and creates the directories structure needed to persist a new resource.
"""
function setup_resource_path(resource_name::String) :: String
  resources_dir = SearchLight.RESOURCES_PATH
  resource_path = joinpath(resources_dir, lowercase(resource_name))

  if ! isdir(resource_path)
    mkpath(resource_path)
    push!(LOAD_PATH, resource_path)
  end

  resource_path
end

"""
    write_resource_file(resource_path::String, file_name::String, resource_name::String) :: Bool

Generates all resouce files and persists them to disk.
"""
function write_resource_file(resource_path::String, file_name::String, resource_name::String, resource_type::Symbol; pluralize::Bool = true) :: Bool
  resource_name = (pluralize ? SearchLight.Inflector.tosingular(resource_name) : resource_name) |> SearchLight.Inflector.from_underscores

  try
    if resource_type == :model
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, SearchLight.FileTemplates.newmodel(resource_name, pluralize = pluralize))
      end
    end
  catch ex
    @error ex
  end

  try
    if resource_type == :validator
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, SearchLight.FileTemplates.newvalidator(resource_name, pluralize = pluralize))
      end
    end
  catch ex
    @error ex
  end

  try
    if resource_type == :test
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        uname = SearchLight.Inflector.from_underscores(resource_name)
        uname = pluralize ? SearchLight.Inflector.to_plural(uname) : uname

        write(f, SearchLight.FileTemplates.newtest(uname, resource_name, pluralize = pluralize))
      end
    end
  catch ex
    @error ex
  end

  true
end

"""
"""
function new_db_config(app_name::String = "App", adapter::Symbol = :mysql; create_folder::Bool = true, folder_name = lowercase(app_name)) :: Nothing
  if create_folder
    mkdir(folder_name)
    cd(folder_name)
  end

  isdir(SearchLight.CONFIG_PATH) || mkpath(SearchLight.CONFIG_PATH)
  isdir(SearchLight.APP_PATH) || mkpath(SearchLight.APP_PATH)
  isdir(SearchLight.config.db_migrations_folder) || mkpath(SearchLight.config.db_migrations_folder)
  if ! ispath(SearchLight.LOG_PATH)
    mkpath(SearchLight.LOG_PATH)
  end

  open(joinpath(SearchLight.DB_PATH, SearchLight.SEARCHLIGHT_DB_CONFIG_FILE_NAME), "w") do f
    write(f, SearchLight.FileTemplates.new_db_config(adapter))
  end

  open(app_name * ".jl", "w") do f
    write(f, SearchLight.FileTemplates.new_app_loader(app_name))
  end

  open(SearchLight.SEARCHLIGHT_BOOTSTRAP_FILE_NAME, "w") do f
    write(f, SearchLight.FileTemplates.new_app_bootstrap(app_name))
  end

  open(SearchLight.SEARCHLIGHT_INFO_FILE_NAME, "w") do f
    write(f, SearchLight.FileTemplates.new_app_info(app_name))
  end

  @info "New app ready at $(pwd())"

  nothing
end

const new_app = new_db_config

"""
"""
function resource_does_not_exist(resource_path::String, file_name::String) :: Bool
  if isfile(joinpath(resource_path, file_name))
    @warn "File already exists, $(joinpath(resource_path, file_name)) - skipping"

    return false
  end

  true
end

"""
"""
function model_file_name(resource_name::Union{String,Symbol})
  string(resource_name) * ".jl"
end

"""
"""
function validator_file_name(resource_name::Union{String,Symbol})
  string(resource_name) * SearchLight.SEARCHLIGHT_VALIDATOR_FILE_POSTFIX
end

"""
"""
function create_migrations_table()
  SearchLight.create_migrations_table()
end

"""
    init() :: Bool

Sets up the DB tables used by SearchLight.
"""
function init() :: Bool
  SearchLight.create_migrations_table(SearchLight.config.db_migrations_table_name)
end

end

const Generators = Generator


text:
"""
Load a DFTK-compatible lattice object from a supported python object (e.g. pymatgen or ASE)
"""
function load_lattice(T, pyobj::PyObject)
    mg = pyimport_e("pymatgen")
    ase = pyimport_e("ase")

    if !ispynull(mg)
        mg_supported = (mg.Structure, mg.Lattice)
        if any(pyisinstance(pyobj, s) for s in mg_supported)
            return load_lattice_pymatgen(T, pyobj)
        end
    end

    if !ispynull(ase)
        ase_supported = (ase.Atoms, ase.cell.Cell)
        if any(pyisinstance(pyobj, s) for s in ase_supported)
            return load_lattice_ase(T, pyobj)
        end
    end

    error("load_lattice not implemented for python type $pyobj")
end

function load_atoms(T, pyobj::PyObject)
    mg = pyimport_e("pymatgen")
    ase = pyimport_e("ase")

    if !ispynull(mg) && pyisinstance(pyobj, mg.Structure)
        return load_atoms_pymatgen(T, pyobj)
    end
    if !ispynull(ase) && pyisinstance(pyobj, ase.Atoms)
        return load_atoms_ase(T, pyobj)
    end

    error("load_atoms not implemented for python type $pyobj")
end


text:
mutable struct BifDiagNode{Tγ, Tc}
	# current level of recursion
	level::Int64

	# code for finding this node in the tree, this is the index of the bifurcation point
	# from which γ branches off
	code::Int64

	# branch associated to the current node
	γ::Tγ

	# childs of current node. These are the different branches off the bifurcation points
	# in γ
	child::Tc
end

@inline level(node::BifDiagNode) = node.level
hasbranch(tree::BifDiagNode) = ~isnothing(tree.γ)
from(tree::BifDiagNode) = from(tree.γ)
add!(tree::BifDiagNode, γ::AbstractBranchResult, level::Int, code::Int) = push!(tree.child, BifDiagNode(level, code, γ, BifDiagNode[]))
add!(tree::BifDiagNode, γ::Vector{ <: AbstractBranchResult}, level::Int, code::Int) = map(x -> add!(tree, x, level, code), γ)
add!(tree::BifDiagNode, γ::Nothing, level::Int, code::Int) = nothing
getContResult(br::ContResult) = br
getContResult(br::Branch) = br.γ
Base.show(io::IO, tree::BifDiagNode) = (println(io, "Bifurcation diagram. Root branch (level $(tree.level)) has $(length(tree.child)) children and is such that:"); show(io, tree.γ))

# total size of the tree
_size(tree::BifDiagNode) = length(tree.child) > 0 ? 1 + mapreduce(size, +, tree.child) : 1

"""
$(SIGNATURES)

Return the size of the bifurcation diagram. The argument `code` is the same as in `getBranch`.
"""
Base.size(tree::BifDiagNode, code = ()) = _size(getBranch(tree, code))

"""
$(SIGNATURES)

Return the part of the tree (bifurcation diagram) by recursively descending down the tree using the `Int` valued tuple `code`. For example `getBranch(tree, (1,2,3,))` returns `tree.child[1].child[2].child[3]`.
"""
function getBranch(tree::BifDiagNode, code)
	isempty(code) && return tree
	return getBranch(tree.child[code[1]], code[2:end])
end

"""
$(SIGNATURES)

Return the part of the tree corresponding to the indbif-th bifurcation point on the root branch.
"""
function getBranchesFromBP(tree::BifDiagNode, indbif::Int)
	# parameter value at the bp
	p = tree.γ.specialpoint[indbif].param
	return BifDiagNode[br for br in tree.child if from(br).p == p]
end

"""
$(SIGNATURES)

Compute the bifurcation diagram associated with the problem `F(x, p) = 0` recursively.

# Arguments
- `F, dF, d2F, d3F` functional and its derivatives
- `x0` initial guess
- `par0` parameter values at `x0`
- `lens` lens to select the parameter axis
- `level` maximum branching (or recursion) level for computing the bifurcation diagram
- `options = (x, p, level) -> contparams` this function allows to change the [`continuation`](@ref) options depending on the branching `level`. The argument `x, p` denotes the current solution to `F(x, p)=0`.
- `kwargs` optional arguments. Look at [`bifurcationdiagram!`](@ref) for more details.

# Simplified call:

We also provide the call

`bifurcationdiagram(F, dF, d2F, d3F, br::ContResult, level::Int, options; kwargs...)`

where `br` is a branch computed after a call to [`continuation`](@ref) from which we want to compute the bifurcating branches recursively.
"""
function bifurcationdiagram(F, dF, d2F, d3F, x0, par0, lens::Lens, level::Int, options; kwargs...)
	γ, u = continuation(F, dF, x0, par0, lens, options(x0, par0, 1); kwargs...)
	bifurcationdiagram(F, dF, d2F, d3F, γ, level, options; code = "0", kwargs...)
end

# TODO, BifDiagNode[] makes it type unstable it seems
function bifurcationdiagram(F, dF, d2F, d3F, br::AbstractBranchResult, maxlevel::Int, options; kwargs...)
	printstyled(color = :magenta, "#"^50 * "\n---> Automatic computation of bifurcation diagram\n\n")
	bifurcationdiagram!(F, dF, d2F, d3F, BifDiagNode(1, 0, br, BifDiagNode[]), maxlevel, options; code = "0", kwargs...)
end

"""
$(SIGNATURES)

Similar to [`bifurcationdiagram`](@ref) but you pass a previously computed `node` from which you want to further compute the bifurcated branches. It is usually used with `node = getBranch(diagram, code)` from a previously computed bifurcation `diagram`.

# Arguments
- `node::BifDiagNode` a node in the bifurcation diagram
- `maxlevel = 1` required maximal level of recursion.
- `options = (x, p, level) -> contparams` this function allows to change the [`continuation`](@ref) options depending on the branching `level`. The argument `x, p` denotes the current solution to `F(x, p)=0`.

# Optional arguments
- `code = "0"`
- `usedeflation = false`
- `halfbranch = false` for Pitchfork/Transcritical bifurcations, compute only half of the branch. Can be useful when there are symmetries.
- `kwargs` optional arguments as for [`continuation`](@ref) but also for the different versions listed in [Continuation](https://rveltz.github.io/BifurcationKit.jl/dev/library/#Continuation-1).
"""
function bifurcationdiagram!(F, dF, d2F, d3F,
		node::BifDiagNode,
		maxlevel::Int,
		options;
		code = "0",
		usedeflation = false,
		halfbranch = false,
		kwargs...)
	if node.level >= maxlevel || isnothing(node.γ); return node; end

	# current level of recursion
	level = node.level

	# convenient function for branching
	function letsbranch(_id, _pt, _level; _dsfactor = 1, _ampfactor = 1)
		plotfunc = get(kwargs, :plotSolution, (x, p; kws...) -> plot!(x; kws...))
		optscont = options(_pt.x, _pt.param, _level + 1)
		@set! optscont.ds *= _dsfactor

		continuation(F, dF, d2F, d3F, getContResult(node.γ), _id, optscont;
			nev = optscont.nev, kwargs...,
			ampfactor = _ampfactor,
			usedeflation = usedeflation,
			plotSolution = (x, p; kws...) -> (plotfunc(x, p; ylabel = code*"-$_id", xlabel = "level = $(_level+1), dim = $(kernelDim(_pt))", label="", kws...);plot!(node.γ; subplot = 1, legend=:topleft, putspecialptlegend = false, markersize = 2)))
	end

	for (id, pt) in enumerate(node.γ.specialpoint)
		# we put this condition in case the specialpoint at step = 0 corresponds to the one we are branching from. If we remove this, we keep computing the same branch (possibly).
		if pt.step > 1 && pt.type in (:bp, :nd)
			try
				println("─"^80*"\n--> New branch, level = $(level+1), dim(Kernel) = ", kernelDim(pt), ", code = $code, from bp #",id," at p = ", pt.param, ", type = ", type(pt))
				γ, = letsbranch(id, pt, level)
				add!(node, γ, level+1, id)
				 ~isnothing(γ) && printstyled(color = :green, "----> From ", type(from(γ)), "\n")

				# in the case of a Transcritical bifurcation, we compute the other branch
				if ~isnothing(γ) && ~(γ isa Vector)
					if ~halfbranch && from(γ) isa Transcritical
						γ, = letsbranch(id, pt, level; _dsfactor = -1)
						add!(node, γ, level+1, id)
					end
					if ~halfbranch && from(γ) isa Pitchfork
						γ, = letsbranch(id, pt, level; _ampfactor = -1)
						add!(node, γ, level+1, id)
					end
				end

			catch ex
				# @error ex
			# 	return node
			end
		end
	end
	for (ii, _node) in enumerate(node.child)
		bifurcationdiagram!(F, dF, d2F, d3F, _node, maxlevel, options; code = code*"-$ii", kwargs...)
	end
	return node
end


text:
function ADC(VT,ALT)
   R0 = 2.377E-3
   TFAC = 1.0 - 0.703E-5*ALT
   T = 519.0*TFAC
   if ALT > 35000.0
      T= 390.0
   end
   RHO = R0*(TFAC^4.14)
   dum = 1.4*1716.3*T
   MACH = VT/sqrt(dum)
   QBAR = 0.5*RHO*VT*VT
   PS = 1715.0*RHO*T
return MACH, QBAR
end


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_translated)" --startup-file=no %s 2>&1 | FileCheck %s

f() = return
emit(f)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.f)}, svec(), f() in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/translate/return_nothing.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.f)}"(%arg0: !jlir<"typeof(Main.f)">) -> !jlir.Nothing attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     "jlir.return"(%0) : (!jlir.Nothing) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
"""
Defines [`Document`](@ref) and its supporting types

- [`Page`](@ref)
- [`User`](@ref)
- [`Internal`](@ref)
- [`Globals`](@ref)

"""
module Documents

import ..Documenter:
    Anchors,
    Formats,
    Utilities

using Compat, DocStringExtensions

import Compat: String

# Pages.
# ------

"""
[`Page`](@ref)-local values such as current module that are shared between nodes in a page.
"""
type Globals
    mod  :: Module
    meta :: Dict{Symbol, Any}
end
Globals() = Globals(current_module(), Dict())

"""
Represents a single markdown file.
"""
immutable Page
    source   :: Compat.String
    build    :: Compat.String
    """
    Ordered list of raw toplevel markdown nodes from the parsed page contents. This vector
    should be considered immutable.
    """
    elements :: Vector
    """
    Each element in `.elements` maps to an "expanded" element. This may be itself if the
    element does not need expanding or some other object, such as a `DocsNode` in the case
    of `@docs` code blocks.
    """
    mapping  :: ObjectIdDict
    globals  :: Globals
end
function Page(source::AbstractString, build::AbstractString)
    elements = Base.Markdown.parse(readstring(source)).content
    Page(source, build, elements, ObjectIdDict(), Globals())
end

# Document Nodes.
# ---------------

## IndexNode.

immutable IndexNode
    pages    :: Vector{String} # Which pages to include in the index? Set by user.
    modules  :: Vector{Module} # Which modules to include? Set by user.
    order    :: Vector{Symbol} # What order should docs be listed in? Set by user.
    build    :: String         # Path to the file where this index will appear.
    source   :: String         # Path to the file where this index was written.
    elements :: Vector         # (object, doc, page, mod, cat)-tuple for constructing links.

    function IndexNode(;
            # TODO: Fix difference between uppercase and lowercase naming of keys.
            #       Perhaps deprecate the uppercase versions? Same with `ContentsNode`.
            Pages   = [],
            Modules = [],
            Order   = [:module, :constant, :type, :function, :macro],
            build   = error("missing value for `build` in `IndexNode`."),
            source  = error("missing value for `source` in `IndexNode`."),
            others...
        )
        new(Pages, Modules, Order, build, source, [])
    end
end

## ContentsNode.

immutable ContentsNode
    pages    :: Vector{String} # Which pages should be included in contents? Set by user.
    depth    :: Int            # Down to which level should headers be displayed? Set by user.
    build    :: String         # Same as for `IndexNode`s.
    source   :: String         # Same as for `IndexNode`s.
    elements :: Vector         # (order, page, anchor)-tuple for constructing links.

    function ContentsNode(;
            Pages  = [],
            Depth  = 2,
            build  = error("missing value for `build` in `ContentsNode`."),
            source = error("missing value for `source` in `ContentsNode`."),
            others...
        )
        new(Pages, Depth, build, source, [])
    end
end

## Other nodes

immutable MetaNode
    dict :: Dict{Symbol, Any}
end

immutable MethodNode
    method  :: Method
    visible :: Bool
end

immutable DocsNode
    docstr  :: Any
    anchor  :: Anchors.Anchor
    object  :: Utilities.Object
    page    :: Documents.Page
end

immutable DocsNodes
    nodes :: Vector{DocsNode}
end

immutable EvalNode
    code   :: Base.Markdown.Code
    result :: Any
end

immutable RawHTML
    code::String
end

immutable RawNode
    name::Symbol
    text::Compat.String
end

# Navigation
# ----------------------

"""
Element in the navigation tree of a document, containing navigation references
to other page, reference to the [`Page`](@ref) object etc.
"""
type NavNode
    """
    `null` if the `NavNode` is a non-page node of the navigation tree, otherwise
    the string should be a valid key in `doc.internal.pages`
    """
    page           :: Nullable{Compat.String}
    """
    If not `null`, specifies the text that should be displayed in navigation
    links etc. instead of the automatically determined text.
    """
    title_override :: Nullable{Compat.String}
    parent         :: Nullable{NavNode}
    children       :: Vector{NavNode}
    visible        :: Bool
    prev           :: Nullable{NavNode}
    next           :: Nullable{NavNode}
end
NavNode(page, title_override, parent) = NavNode(page, title_override, parent, [], true, nothing, nothing)

"""
Constructs a list of the ancestors of the `navnode` (inclding the `navnode` itself),
ordered so that the root of the navigation tree is the first and `navnode` itself
is the last item.
"""
navpath(navnode::NavNode) = isnull(navnode.parent) ? [navnode] :
    push!(navpath(get(navnode.parent)), navnode)

# Inner Document Fields.
# ----------------------

"""
User-specified values used to control the generation process.
"""
immutable User
    root    :: Compat.String  # An absolute path to the root directory of the document.
    source  :: Compat.String  # Parent directory is `.root`. Where files are read from.
    build   :: Compat.String  # Parent directory is also `.root`. Where files are written to.
    format  :: Vector{Symbol} # What format to render the final document with?
    clean   :: Bool           # Empty the `build` directory before starting a new build?
    doctest :: Bool           # Run doctests?
    linkcheck::Bool           # Check external links..
    linkcheck_ignore::Vector{Union{String,Regex}}  # ..and then ignore (some of) them.
    checkdocs::Symbol         # Check objects missing from `@docs` blocks. `:none`, `:exports`, or `:all`.
    strict::Bool              # Throw an exception when any warnings are encountered.
    modules :: Set{Module}    # Which modules to check for missing docs?
    pages   :: Vector{Any}    # Ordering of document pages specified by the user.
    assets  :: Vector{Compat.String}
    repo    :: Compat.String  # Template for URL to source code repo
    sitename:: Compat.String
    authors :: Compat.String
    analytics::Compat.String
    version :: Compat.String # version string used in the version selector by default
    html_prettyurls :: Bool # Use pretty URLs in the HTML build?
end

"""
Private state used to control the generation process.
"""
immutable Internal
    assets  :: Compat.String             # Path where asset files will be copied to.
    remote  :: Compat.String             # The remote repo on github where this package is hosted.
    pages   :: Dict{Compat.String, Page} # Markdown files only.
    navtree :: Vector{NavNode}           # A vector of top-level navigation items.
    navlist :: Vector{NavNode}           # An ordered list of `NavNode`s that point to actual pages
    headers :: Anchors.AnchorMap         # See `modules/Anchors.jl`. Tracks `Markdown.Header` objects.
    docs    :: Anchors.AnchorMap         # See `modules/Anchors.jl`. Tracks `@docs` docstrings.
    bindings:: ObjectIdDict              # Tracks insertion order of object per-binding.
    objects :: ObjectIdDict              # Tracks which `Utilities.Objects` are included in the `Document`.
    contentsnodes :: Vector{ContentsNode}
    indexnodes    :: Vector{IndexNode}
    locallinks :: Dict{Base.Markdown.Link, Compat.String}
    errors::Set{Symbol}
end

# Document.
# ---------

"""
Represents an entire document.
"""
immutable Document
    user     :: User     # Set by the user via `makedocs`.
    internal :: Internal # Computed values.
end

function Document(;
        root     :: AbstractString   = Utilities.currentdir(),
        source   :: AbstractString   = "src",
        build    :: AbstractString   = "build",
        format   :: Any              = :markdown,
        clean    :: Bool             = true,
        doctest  :: Bool             = true,
        linkcheck:: Bool             = false,
        linkcheck_ignore :: Vector   = [],
        checkdocs::Symbol            = :all,
        strict::Bool                 = false,
        modules  :: Utilities.ModVec = Module[],
        pages    :: Vector           = Any[],
        assets   :: Vector           = Compat.String[],
        repo     :: AbstractString   = "",
        sitename :: AbstractString   = "",
        authors  :: AbstractString   = "",
        analytics :: AbstractString = "",
        version :: AbstractString = "",
        html_prettyurls :: Bool = false,
        others...
    )
    Utilities.check_kwargs(others)

    local fmt = Formats.fmt(format)
    @assert !isempty(fmt) "No formats provided."

    if version == "git-commit"
        version = "git:$(Utilities.get_commit_short(root))"
    end

    user = User(
        root,
        source,
        build,
        fmt,
        clean,
        doctest,
        linkcheck,
        linkcheck_ignore,
        checkdocs,
        strict,
        Utilities.submodules(modules),
        pages,
        assets,
        repo,
        sitename,
        authors,
        analytics,
        version,
        html_prettyurls,
    )
    internal = Internal(
        Utilities.assetsdir(),
        Utilities.getremote(root),
        Dict{Compat.String, Page}(),
        [],
        [],
        Anchors.AnchorMap(),
        Anchors.AnchorMap(),
        ObjectIdDict(),
        ObjectIdDict(),
        [],
        [],
        Dict{Base.Markdown.Link, Compat.String}(),
        Set{Symbol}(),
    )
    Document(user, internal)
end

## Methods

function addpage!(doc::Document, src::AbstractString, dst::AbstractString)
    page = Page(src, dst)
    # page's identifier is the path relative to the `doc.user.source` directory
    name = normpath(relpath(src, doc.user.source))
    doc.internal.pages[name] = page
end

"""
$(SIGNATURES)

Populates the `ContentsNode`s and `IndexNode`s of the `document` with links.

This can only be done after all the blocks have been expanded (and nodes constructed),
because the items have to exist before we can gather the links to those items.
"""
function populate!(document::Document)
    for node in document.internal.contentsnodes
        populate!(node, document)
    end
    for node in document.internal.indexnodes
        populate!(node, document)
    end
end

function populate!(index::IndexNode, document::Document)
    # Filtering valid index links.
    for (object, doc) in document.internal.objects
        page = relpath(doc.page.build, dirname(index.build))
        mod  = object.binding.mod
        # Include *all* signatures, whether they are `Union{}` or not.
        cat  = Symbol(lowercase(Utilities.doccat(object.binding, Union{})))
        if _isvalid(page, index.pages) && _isvalid(mod, index.modules) && _isvalid(cat, index.order)
            push!(index.elements, (object, doc, page, mod, cat))
        end
    end
    # Sorting index links.
    pagesmap   = precedence(index.pages)
    modulesmap = precedence(index.modules)
    ordermap   = precedence(index.order)
    comparison = function(a, b)
        (x = _compare(pagesmap,   3, a, b)) == 0 || return x < 0 # page
        (x = _compare(modulesmap, 4, a, b)) == 0 || return x < 0 # module
        (x = _compare(ordermap,   5, a, b)) == 0 || return x < 0 # category
        string(a[1].binding) < string(b[1].binding)              # object name
    end
    sort!(index.elements, lt = comparison)
    return index
end

function populate!(contents::ContentsNode, document::Document)
    # Filtering valid contents links.
    for (id, filedict) in document.internal.headers.map
        for (file, anchors) in filedict
            for anchor in anchors
                page = relpath(anchor.file, dirname(contents.build))
                if _isvalid(page, contents.pages) && Utilities.header_level(anchor.object) ≤ contents.depth
                    push!(contents.elements, (anchor.order, page, anchor))
                end
            end
        end
    end
    # Sorting contents links.
    pagesmap   = precedence(contents.pages)
    comparison = function(a, b)
        (x = _compare(pagesmap, 2, a, b)) == 0 || return x < 0 # page
        a[1] < b[1]                                            # anchor order
    end
    sort!(contents.elements, lt = comparison)
    return contents
end

## Utilities.

function buildnode(T::Type, block, doc, page)
    mod  = get(page.globals.meta, :CurrentModule, current_module())
    dict = Dict{Symbol, Any}(:source => page.source, :build => page.build)
    for (ex, str) in Utilities.parseblock(block.code, doc, page)
        if Utilities.isassign(ex)
            cd(dirname(page.source)) do
                dict[ex.args[1]] = eval(mod, ex.args[2])
            end
        end
    end
    T(; dict...)
end

function _compare(col, ind, a, b)
    x, y = a[ind], b[ind]
    haskey(col, x) && haskey(col, y) ? _compare(col[x], col[y]) : 0
end
_compare(a, b)  = a < b ? -1 : a == b ? 0 : 1
_isvalid(x, xs) = isempty(xs) || x in xs
precedence(vec) = Dict(zip(vec, 1:length(vec)))

end


text:
module DistSlicingExample
    using PicoQuant
    using MPI

    MPI.Init()
    comm = MPI.COMM_WORLD
    number_partitions = MPI.Comm_size(comm)
    partition = MPI.Comm_rank(comm) + 1

    MPI.Barrier(comm)

    n = 4
    circ = create_simple_preparation_circuit(n, 2).compose(create_qft_circuit(n))

    if partition == 1
        # contract the regular way
        tn = convert_qiskit_circ_to_network(circ, decompose=true, transpile=true)
        add_input!(tn, "0"^n)
        wf = full_wavefunction_contraction!(tn, "vector")
    end

    tn = convert_qiskit_circ_to_network(circ, decompose=true, transpile=true)
    add_input!(tn, "0"^n)
    bond_labels, bond_values = partition_network_on_virtual_bonds(tn,
                                                                  number_partitions,
                                                                  partition)
    slice_tensor_network(tn, bond_labels, bond_values)
    wf_part = full_wavefunction_contraction!(tn, "vector")

    MPI.Reduce!(wf_part, MPI.SUM, 0, MPI.COMM_WORLD)
    if partition == 1
        if wf ≈ wf_part
            println("Results match with overlap: ", sum(wf .* conj(wf)))
        else
            println("Result mismatch with overlap: ", sum(wf .* conj(wf)))
        end
    end
end