text:
# Simplification.jl

# [[file:~/Documents/Julia/scrap.org::*Simplification.jl][Simplification.jl:1]]
#----------------------------------------------------------------------------
# Some utility functions
#----------------------------------------------------------------------------

function Base.vcat(ex::T, arr::Array) where {T<:AbstractSymExpr}
    T(ex.op, vcat(ex.args, arr))
end

function trim(ex::T, index::Int) where {T<:AbstractSymExpr}
    T(ex.op, ex.args[1:end .!= index])
end

function find_first_duplicates(a::Array)
    for i in 1:(length(a)-1)
        out = [i]
        for j in (i+1):length(a)
            if a[i] == a[j]
                push!(out, j)
            end
            if (j == length(a)) && (length(out) > 1)
                return out
            end
        end
    end
    []
end

apply_f_to_pairs(f, x) = x
function apply_f_to_pairs(f, ex::T) where {T<:AbstractSymExpr}
    if length(ex.args) == 2
        return f(ex)
    elseif length(ex.args) > 2
        for i in 1:(length(ex.args)-1), j in (i+1):length(ex.args)
            old_arg = SymExpr(ex.op, [ex.args[i], ex.args[j]])
            newarg = f(old_arg)
            if newarg != old_arg
                newarglist = ex.args[1:end .!= i][1:end .!= j-1] #This should select everything but the i and j components since j>i
                return SymExpr(ex.op, vcat([newarg], newarglist))
            end
        end
    end
    ex
end

#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# The main loop for recursively going through a symbolic expression and simplifying it
#----------------------------------------------------------------------------

simplify(x) = x
function simplify(ex::AbstractSymExpr)
    out1 = simplify_pass(ex)
    out2 = simplify_pass(out1)
    counter = 1
    while (out2 isa SymExpr) && (out1 isa SymExpr) && (out2 != out1)
        out1 = simplify_pass(out2)
        out2 = simplify_pass(out1)
        counter += 1
        if counter > 1000
            @warn "possible infinite loop in simplification rules. Breaking"
            return out2
        end
    end
    out2
end

simplify_pass(x) = x
function simplify_pass(ex::AbstractSymExpr)
    postwalk(apply_algebraic_rules, ex)
end

apply_algebraic_rules(x) = x
function apply_algebraic_rules(expr::AbstractSymExpr)
    @> expr begin
        stripiden
        add_rules
        mult_rules
        exp_rules
        distribute_negation
        denest
        collect_numeric_terms
        collect_identical
        eval_numeric
    end
end

#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# The simplification functions
#----------------------------------------------------------------------------

#----------------------------------------------------------------------------
# Rules for Addition
#----------------------------------------------------------------------------
add_rules(x) = x
function add_rules(ex::T) where {T<:AbstractSymExpr}
    @> ex begin
        distribute_negation
        remove_addative_identity
        ex -> apply_f_to_pairs(factor_addition, ex)
    end
end

remove_addative_identity(x) = x
function remove_addative_identity(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == :+) && (0 in ex.args)
        return +(ex.args[ex.args .!= 0]...)
    end
    ex
end

distribute_negation(x) = x
function distribute_negation(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == Sym(:-)) && (length(ex.args) == 1) && (ex.args[1] isa AbstractSymExpr) && (ex.args[1].op == Sym(:+))
        T(:+, [-(arg) for arg in ex.args[1].args])
    else
        ex
    end
end

factor_addition(ex) = ex
function factor_addition(ex::T) where {T<:AbstractSymExpr}
    @match ex begin
        #SymExpr(:+, [x, SymExpr(:*, [a, x])]) => (a+1)*x
        SymExpr(:+, [SymExpr(:*, [a,b]), SymExpr(:*, [c,d])]), if b==d end => (a+c)*b
        SymExpr(:+, [SymExpr(:*, [a,b]), SymExpr(:*, [c,d])]), if b==c end => (a+d)*b
        SymExpr(:+, [SymExpr(:*, [a,b]), SymExpr(:*, [c,d])]), if a==c end => (b+d)*a
        SymExpr(:+, [SymExpr(:*, [a,b]), SymExpr(:*, [c,d])]), if a==d end => (b+c)*a
        SymExpr(:+, [a, SymExpr(:*, [b, c])]), if a==b end                 => (c+1)*a
        SymExpr(:+, [a, SymExpr(:*, [b, c])]), if a==c end                 => (b+1)*a
        _                                                                  => ex
    end
end

#----------------------------------------------------------------------------
# Rules for Multiplication
#----------------------------------------------------------------------------
mult_rules(x) = x
function mult_rules(ex::T) where {T<:AbstractSymExpr}
    @> ex begin
        mult_by_zero
        mult_by_one
        x -> apply_f_to_pairs(factor_powers, x)
    end
end

mult_by_zero(ex) = ex
function mult_by_zero(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == Sym(:*)) && (0 in ex.args)
        0
    else
        ex
    end
end

mult_by_one(ex) = ex
function mult_by_one(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == Sym(:*)) && (1 in ex.args)
        newargs = ex.args[ex.args .!= 1]
        if newargs == []
            1
        else
            *(newargs...)
        end
    else
        ex
    end
end

factor_powers(x) = x
function factor_powers(ex::T) where {T<:AbstractSymExpr}
    @match ex begin
        SymExpr(:*, [SymExpr(:^, [a,b]), SymExpr(:^, [c,d])]), if a==c end => a^(b+d)
        SymExpr(:*, [a, SymExpr(:^, [b,c])]),                  if a==b end => a^(c+1)
        SymExpr(:*, [SymExpr(:^, [a,b]), c]),                  if a==c end => a^(b+1)
        SymExpr(:*, [a, b]),                                   if a==b end => a^2
        _                                                                  => ex
    end
end

#----------------------------------------------------------------------------
# Rules for Exponentiation
#----------------------------------------------------------------------------
exp_rules(x) = x
function exp_rules(ex::T) where {T<:AbstractSymExpr}
    @> ex begin
        exp_one
        exp_zero
        exp_exp
    end
end

exp_one(x) = x
function exp_one(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == Sym(:^)) && (ex.args[2] == 1)
        return ex.args[1]
    end
    ex
end

exp_zero(x) = x
function exp_zero(ex::T) where {T<:AbstractSymExpr}
    if (ex.op == Sym(:^)) && (ex.args[2] == 0)
        return 1
    end
    ex
end

exp_exp(x) = x
function exp_exp(ex::T) where {T<:AbstractSymExpr}
     @match ex begin
         T(:^, [T(:^, [a,b]), c]) => a^(b*c)
         _                        => ex
     end
end

#----------------------------------------------------------------------------
# More Misc. Rules
#----------------------------------------------------------------------------
eval_numeric(x) = x
function eval_numeric(ex::T) where {T<:AbstractSymExpr}
    if ([!(arg isa Symbolic) for arg in ex.args] == trues(length(ex.args))) && (try eval(Symbol(ex.op)) catch e; false end isa Function)
        eval(Expr(ex))
    else
        ex
    end
end

remove_identity_operations(x) = x
function remove_identity_operations(expr::T) where {T<:AbstractSymExpr}
    if (expr.op == Sym(:^)) && (expr.args[2] == 1)
        return expr.args[1]
    elseif expr.op == :+
        if (0 in expr.args)
            return +(expr.args[expr.args .!= 0]...)
        end

    end
    expr
end

denest(x) = x
function denest(ex::AbstractSymExpr)
    (ex.op in [Sym(:+), Sym(:*)]) ? denest_op(ex, ex.op) : ex
end

function denest_op(ex, op)
    for arg in ex.args
        if (arg isa AbstractSymExpr) && (arg.op == op)
            ii = findall(ex.args .== arg)[1]
            return typeof(arg)(op, vcat(arg.args, ex.args[1:end .!= ii]))
        end
    end
    ex
end

collect_numeric_terms(x) = x
function collect_numeric_terms(ex::T) where {T<:SymExpr}
    if (ex.op in [:+, :*]) && (length(ex.args) > 2)
        is_numeric = isa.(ex.args, Union{Real, Complex})
        numbers = ex.args[is_numeric]
        if length(numbers) > 1
            new_arg = eval(ex.op.name)(numbers...)
            T(ex.op, [new_arg, ex.args[.~(is_numeric)]...])
        else
            ex
        end
    else
        ex
    end
end

collect_identical(x) = x
function collect_identical(ex::T) where {T<:AbstractSymExpr}
    if (ex.op in [:+, :*]) && (length(ex.args) > 1)
        dup_inds = find_first_duplicates(ex.args)
        is_dup = [(i in dup_inds) for i in 1:length(ex.args)]
        identical_terms = ex.args[is_dup]
        if length(identical_terms) > 1
            if ex.op == :+
                new_arg = length(identical_terms)*identical_terms[1]
            elseif ex.op == :*
                new_arg = identical_terms[1]^length(identical_terms)
            end
            eval(ex.op.name)(new_arg, ex.args[.~(is_dup)]...)
        else
            ex
        end
    else
        ex
    end
end
# Simplification.jl:1 ends here


text:
# Code for solving assuming that all other cars have a single identical behavior model

type SingleBehaviorSolver <: Solver
    inner_solver::Solver
    behavior::BehaviorModel
end

abstract type Policy{MLState} end
struct SingleBehaviorPolicy{Policy,BehaviorModel} #<: Policy{MLState}
    inner_policy::Policy
    behavior::BehaviorModel
end

#
#type SingleBehaviorPolicy <: Policy{MLState}
#    inner_policy::Policy{MLState}
#    behavior::BehaviorModel
#end

set_rng!(solver::SingleBehaviorSolver, rng::AbstractRNG) = set_rng!(solver.inner_solver, rng)

function solve(solver::SingleBehaviorSolver, mdp::NoCrashProblem)
    single_behavior_mdp = deepcopy(mdp)
    single_behavior_mdp.dmodel.behaviors = DiscreteBehaviorSet(BehaviorModel[solver.behavior], Weights([1.0]))

    inner_policy = solve(solver.inner_solver, single_behavior_mdp)
    return SingleBehaviorPolicy(inner_policy, solver.behavior)
end

function action(p::SingleBehaviorPolicy, s::MLState, a::MLAction=MLAction())
    as = actions(p.inner_policy.mdp, s, actions(p.inner_policy.mdp))
    if length(as) == 1
        return first(as)
    end
    return action(p.inner_policy, single_behavior_state(s, p.behavior))
end

function single_behavior_state(s::MLState, behavior)
    new_cars = Array(CarState, length(s.cars))
    for (i,c) in enumerate(s.cars)
        new_cars[i] = CarState(c.x, c.y, c.vel, c.lane_change, behavior, c.id)
    end
    return MLState(s, new_cars)
end

ext: jl

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
module GraphNeuralNetworks

using Statistics: mean
using LinearAlgebra, Random
using SparseArrays
import KrylovKit
using Base: tail
using CUDA
using Flux
using Flux: glorot_uniform, leakyrelu, GRUCell, @functor, batch
using MacroTools: @forward
import LearnBase
using LearnBase: getobs
using NNlib, NNlibCUDA
using ChainRulesCore
import LightGraphs
using LightGraphs: AbstractGraph, outneighbors, inneighbors, adjacency_matrix, degree

export
    # gnngraph
    GNNGraph,
    edge_index,
    node_features, edge_features, graph_features,
    adjacency_list, normalized_laplacian, scaled_laplacian,
    add_self_loops, remove_self_loops,
    getgraph,

    # from LightGraphs
    adjacency_matrix,
    # from SparseArrays
    sprand, sparse,

    # msgpass
    update_node, update_edge, compute_message, propagate,

    # layers/basic
    GNNLayer,
    GNNChain,

    # layers/conv
    ChebConv,
    EdgeConv,
    GATConv,
    GatedGraphConv,
    GCNConv,
    GINConv,
    GraphConv,
    NNConv,
    ResGatedGraphConv,
    SAGEConv,

    # layers/pool
    GlobalPool,
    TopKPool,
    topk_index

include("gnngraph.jl")
include("graph_conversions.jl")
include("utils.jl")
include("msgpass.jl")
include("layers/basic.jl")
include("layers/conv.jl")
include("layers/pool.jl")
include("deprecations.jl")

end


text:
function run_config(config, init=Dict())
  e = Interpreter(config, "test_salt", init)
  get_params(e)
end

function run_config_single(config)
  x_config = Dict("op" => "set", "var" => "x", "value" => config)
  getindex(run_config(x_config), "x")
end

expected_val = "x_val"
c = Dict("op" => "set", "value" => expected_val, "var" => "x")
d = run_config(c)
@test getindex(d.data, "x") == expected_val

expected_x_val = "x_val"
expected_y_val = "y_val"
c = Dict("op" => "seq", "seq" =>
      [Dict("op" => "set", "value" => expected_x_val, "var" => "x"),
       Dict("op" => "set", "value" => expected_y_val, "var" => "y")
     ])
d = run_config(c)
@test getindex(d.data, "x") == expected_x_val
@test getindex(d.data, "y") == expected_y_val

arr = [4, 5, "a"]
d = run_config_single(Dict("op" => "array", "values" => arr))
@test d == arr

x = run_config_single(Dict("op" => "coalesce", "values" => [Nullable()]))
@test isequal(x, Nullable())

x = run_config_single(Dict("op" => "coalesce", "values" => [Nullable(), 42, Nullable()]))
@test x == 42

x = run_config_single(Dict("op" => "or", "values" => [0, 0, 0]))
@test x == false

x = run_config_single(Dict("op" => "or", "values" => [0, 0, 1]))
@test x == true

x = run_config_single(Dict("op" => "and", "values" => [1, 1, 0]))
@test x == false

x = run_config_single(Dict("op" => "and", "values" => [1, 1, 1]))
@test x == true

map_val = Dict("a" => 2, "b" => "c", "d" => false)
map_op = run_config_single(Dict("op" => "map", "a" => 2, "b" => "c", "d" => false))
@test map_op == map_val

empty_map = Dict()
map_op2 = run_config_single(Dict("op" => "map"))
@test empty_map == map_op2

arr = [33, 7, 18, 21, -3]

sum_test = run_config_single(Dict("op" => "sum", "values" => arr))
@test sum_test == 76

product_test = run_config_single(Dict("op" => "product", "values" => arr))
@test product_test == -261954

min_test = run_config_single(Dict("op" => "min", "values" => arr))
@test min_test == -3

max_test = run_config_single(Dict("op" => "max", "values" => arr))
@test max_test == 33

core_dict = Dict("op" => "equals", "left" => 1, "right" => 2)
eq = run_config_single(core_dict)
@test eq == false

eq = run_config_single(Dict("op" => "equals", "left" => 2, "right" => 2))
@test eq == true

core_dict["op"] = ">"
gt = run_config_single(core_dict)
@test gt == false

core_dict["op"] = ">="
gte = run_config_single(core_dict)
@test gte == false

core_dict["op"] = "<"
lt = run_config_single(core_dict)
@test lt == true

core_dict["op"] = "<="
lte = run_config_single(core_dict)
@test lte == true

mod = run_config_single(Dict("op" => "%", "left" => 11, "right" => 3))
@test mod == (11 % 3)

div = run_config_single(Dict("op" => "/", "left" => 3, "right" => 4))
@test div == 0.75

round = run_config_single(Dict("op" => "round", "value" => 3.1))
@test round == 3.0

arr = [0, 1, 2, 3]
length_test = run_config_single(Dict("op" => "length", "value" => arr))
@test length(arr) == length_test

length_test = run_config_single(Dict("op" => "length", "value" => []))
@test 0 == 0

not_test = run_config_single(Dict("op" => "not", "value" => 0))
@test not_test == true

not_test = run_config_single(Dict("op" => "not", "value" => false))
@test not_test == true

not_test = run_config_single(Dict("op" => "not", "value" => 1))
@test not_test == false

not_test = run_config_single(Dict("op" => "not", "value" => true))
@test not_test == false

array_literal = [10, 20, 30]
dict_literal = Dict("a" => 42, "b" => 43)
x = run_config_single(Dict("op" => "index", "index" => 0, "base" => array_literal))
@test x == 10

x = run_config_single(Dict("op" => "index", "index" => 2, "base" => array_literal))
@test x == 30

x = run_config_single(Dict("op" => "index", "index" => "a", "base" => dict_literal))
@test x == 42

x = run_config_single(Dict("op" => "index", "index" => 6, "base" => array_literal))
@test isequal(x, Nullable())

x = run_config_single(Dict("op" => "index", "index" => "c", "base" => dict_literal))
@test isequal(x, Nullable())

x = run_config_single(Dict("op" => "index", "index" => 2, "base" => Dict("op" => "array", "values" => array_literal)))
@test x == 30

function return_runner(val)
  config = Dict("op" => "seq", "seq" =>
              [Dict("op" => "set", "var" => "x", "value" => 2),
               Dict("op" => "return", "value" => val),
               Dict("op" => "set", "var" => "y", "value" => 4)]
           )
  Interpreter(config, "test_salt")
end

i = return_runner(true)
@test getindex(get_params(i), "x") == 2
@test i.in_experiment == true

i = return_runner(42)
@test getindex(get_params(i), "x") == 2
@test i.in_experiment == true

i = return_runner(false)
@test getindex(get_params(i), "x") == 2
@test i.in_experiment == false

i = return_runner(0)
@test getindex(get_params(i), "x") == 2
@test i.in_experiment == false


text:
"""
    MatrixWorkspace

This is a data structure for the efficient repeated solution of a square or
overdetermined linear system `Ax=b`.
"""
struct MatrixWorkspace{M<:AbstractMatrix{ComplexF64}} <: AbstractMatrix{ComplexF64}
    A::M # Matrix
    d::Vector{Float64} # Inverse of scaling factors
    factorized::Base.RefValue{Bool}
    lu::LA.LU{ComplexF64,M} # LU Factorization of D * J
    qr::LA.QR{ComplexF64,Matrix{ComplexF64}}
    row_scaling::Vector{Float64}
    scaled::Base.RefValue{Bool}
    # mixed precision iterative refinement
    x̄::Vector{ComplexDF64}
    r::Vector{ComplexF64}
    r̄::Vector{ComplexDF64}
    δx::Vector{ComplexF64}
    inf_norm_est_work::Vector{ComplexF64}
    inf_norm_est_rwork::Vector{Float64}
end

MatrixWorkspace(m::Integer, n::Integer; kwargs...) =
    MatrixWorkspace(zeros(ComplexF64, m, n); kwargs...)
function MatrixWorkspace(Â::AbstractMatrix; optimize_data_structure = true)
    m, n = size(Â)
    m ≥ n || throw(ArgumentError("Expected system with more rows than columns."))

    A = Matrix{ComplexF64}(Â)
    d = ones(m)
    factorized = Ref(false)
    qr = LA.qrfactUnblocked!(copy(A))
    # experiments show that for m > 25 the data layout as a
    # struct array is beneficial
    if m > 25 && optimize_data_structure
        A = StructArrays.StructArray(A)
    end
    row_scaling = ones(m)
    scaled = Ref(false)

    lu = LinearAlgebra.LU{eltype(A),typeof(A)}(copy(A), zeros(Int, m), 0)
    r = zeros(ComplexF64, m)
    r̄ = zeros(ComplexDF64, m)
    x̄ = zeros(ComplexDF64, n)
    δx = zeros(ComplexF64, n)
    inf_norm_est_work = Vector{ComplexF64}(undef, n)
    inf_norm_est_rwork = Vector{Float64}(undef, n)

    MatrixWorkspace(
        A,
        d,
        factorized,
        lu,
        qr,
        row_scaling,
        scaled,
        x̄,
        r,
        r̄,
        δx,
        inf_norm_est_work,
        inf_norm_est_rwork,
    )
end

Base.size(MW::MatrixWorkspace) = size(MW.A)

import Base: @propagate_inbounds
@propagate_inbounds Base.getindex(MW::MatrixWorkspace, i::Integer) = getindex(MW.A, i)
@propagate_inbounds Base.getindex(MW::MatrixWorkspace, i::Integer, j::Integer) =
    getindex(MW.A, i, j)
@propagate_inbounds Base.setindex!(MW::MatrixWorkspace, x, i::Integer) =
    setindex!(MW.A, x, i)
@propagate_inbounds Base.setindex!(MW::MatrixWorkspace, x, i::Integer, j::Integer) =
    setindex!(MW.A, x, i, j)
@propagate_inbounds Base.copyto!(MW::MatrixWorkspace, A::AbstractArray) = copyto!(MW.A, A)

matrix(M::MatrixWorkspace) = M.A
matrix(M::AbstractMatrix) = M

"""
    updated!(MW::MatrixWorkspace)

Indicate that the matrix `MW` got updated.
"""
function updated!(MW::MatrixWorkspace)
    MW.factorized[] = false
    MW.scaled[] = false
    m, n = size(MW)
    if m == n
        @inbounds copyto!(MW.lu.factors, MW.A)
    else
        @inbounds copyto!(MW.qr.factors, MW.A)
    end
    MW
end
updated!(M::AbstractMatrix) = M

"""
    update!(MW::MatrixWorkspace, A::Matrix)

Update the matrix in `MW` with `A`.
"""
@inline function update!(MW::AbstractMatrix, A::Matrix)
    @boundscheck (size(MW.A) == size(A) || throw(ArgumentError("Matrix of invalid size.")))
    @inbounds copyto!(MW.A, A)
    updated!(MW)
end

##################
# Factorizations #
##################

###########################
# CUSTOM LU Factorization #
###########################
# This is an adoption of LA.generic_lufact!
# with 3 changes:
# 1) For choosing the pivot we use abs2 instead of abs
# 2) Instead of using the robust complex division we
#    just use the naive division. This shouldn't
#    lead to problems since numerical diffuculties only
#    arise for very large or small exponents
# 3) We fold lu! and ldiv! into one routine
#    this has the effect that we do not need to allocate
#    the pivot vector anymore and also avoid the allocations
#    coming from the LU wrapper
function lu!(
    A::AbstractMatrix{T},
    ipiv::Union{Vector{I},Nothing} = nothing,
) where {T,I<:Integer}
    m, n = size(A)
    minmn = min(m, n)
    # LU Factorization
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if T <: Complex
                amax = abs2(A[k, k])
            else
                amax = abs(A[k, k])
            end
            for i = k+1:m
                if T <: Complex
                    absi = abs2(A[i, k])
                else
                    absi = abs(A[i, k])
                end
                kp, amax = ifelse(absi > amax, (i, absi), (kp, amax))
            end
            if !(ipiv isa Nothing)
                ipiv[k] = kp
            end
            if !iszero(amax)
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k, i]
                        A[k, i] = A[kp, i]
                        A[kp, i] = tmp
                    end
                end
                # Scale first column
                Akk = A[k, k]
                for i = k+1:m
                    # we assume the matrix is decently scaled
                    # so that the naive division algorithm doesn't overflow
                    @fastmath A[i, k] = A[i, k] / Akk
                end
            end
            # Update the rest
            for j = k+1:n
                A_kj = A[k, j]
                for i = k+1:m
                    A[i, j] -= A[i, k] * A_kj
                end
            end
        end
    end
    A
end

####################
# QR Factorization #
####################
# This is the generic base implementation with the small change that we introduce
# a @fastmath for the division step since results in an approx 50% speedup

# Elementary reflection similar to LAPACK. The reflector is not Hermitian but
# ensures that tridiagonalization of Hermitian matrices become real. See lawn72
@inline function reflector!(x::AbstractVector)
    n = length(x)
    @inbounds begin
        ξ1 = x[1]
        normu = abs2(ξ1)
        for i = 2:n
            normu += abs2(x[i])
        end
        if iszero(normu)
            return zero(ξ1 / normu)
        end
        normu = sqrt(normu)
        ν = copysign(normu, real(ξ1))
        ξ1 += ν
        x[1] = -ν
        for i = 2:n
            @fastmath x[i] = x[i] / ξ1
        end
    end
    ξ1 / ν
end

# apply reflector from left
@inline function reflectorApply!(x::AbstractVector, τ::Number, A::StridedMatrix)
    m, n = size(A)
    @inbounds begin
        for j = 1:n
            # dot
            vAj = A[1, j]
            for i = 2:m
                vAj += x[i]' * A[i, j]
            end

            vAj = conj(τ) * vAj

            # ger
            A[1, j] -= vAj
            for i = 2:m
                A[i, j] -= x[i] * vAj
            end
        end
    end
    return A
end

function qr!(qr::LA.QR{T}) where {T}
    A = qr.factors
    τ = qr.τ
    m, n = size(A)
    for k = 1:min(m - 1 + !(T <: Real), n)
        x = view(A, k:m, k)
        τk = reflector!(x)
        τ[k] = τk
        reflectorApply!(x, τk, view(A, k:m, k+1:n))
    end
    qr
end

function factorize!(WS::MatrixWorkspace)
    m, n = size(WS)
    if m == n
        lu!(WS.lu.factors, WS.lu.ipiv)
    else
        qr!(WS.qr)
    end
    WS.factorized[] = true
    WS
end

##########
## ldiv ##
##########
@inline _ipiv!(A::LA.LU, b::AbstractVector) = apply_ipiv!(b, 1:length(A.ipiv), A.ipiv)
@inline _inverse_ipiv!(A::LA.LU, b::StridedVecOrMat) =
    apply_ipiv!(b, length(A.ipiv):-1:1, A.ipiv)
@inline function apply_ipiv!(b::AbstractVector, range::OrdinalRange, ipiv)
    @inbounds for i in range
        if i != ipiv[i]
            _swap_rows!(b, i, ipiv[i])
        end
    end
    b
end
@propagate_inbounds function _swap_rows!(B::StridedVector, i::Integer, j::Integer)
    B[i], B[j] = B[j], B[i]
    B
end

@inline function ldiv_upper!(A::AbstractMatrix, b::AbstractVector, x::AbstractVector = b)
    n = size(A, 2)
    @inbounds for j = n:-1:1
        # singular_exception && iszero(A[j, j]) && throw(LA.SingularException(j))
        xj = x[j] = (@fastmath A[j, j] \ b[j])
        for i = 1:(j-1)
            b[i] -= A[i, j] * xj
        end
    end
    b
end
@inline function ldiv_unit_lower!(
    A::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector = b,
)
    n = size(A, 2)
    @inbounds for j = 1:n
        xj = x[j] = b[j]
        for i = j+1:n
            b[i] -= A[i, j] * xj
        end
    end
    x
end

function lu_ldiv!(x, LU::LA.LU, b::AbstractVector)
    x === b || copyto!(x, b)
    _ipiv!(LU, x)
    ldiv_unit_lower!(LU.factors, x)
    ldiv_upper!(LU.factors, x)
    x
end

function ldiv_adj_unit_lower!(A::AbstractMatrix, b::AbstractVector, x::AbstractVector = b)
    n = size(A, 1)
    @inbounds for j = n:-1:1
        z = b[j]
        for i = n:-1:j+1
            z -= conj(A[i, j]) * x[i]
        end
        x[j] = z
    end
    x
end

function ldiv_adj_upper!(
    A::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector = b;
    singular_exception::Bool = false,
)
    n = size(A, 1)
    @inbounds for j = 1:n
        z = b[j]
        for i = 1:j-1
            z -= conj(A[i, j]) * x[i]
        end
        iszero(A[j, j]) && singular_exception && throw(SingularException(j))
        x[j] = @fastmath conj(A[j, j]) \ z
    end
    x
end

function lu_ldiv_adj!(x, LU::LA.LU, b::AbstractVector; check::Bool = false)
    x === b || copyto!(x, b)
    ldiv_adj_upper!(LU.factors, x; singular_exception = check)
    ldiv_adj_unit_lower!(LU.factors, x)
    _inverse_ipiv!(LU, x)
    x
end

# QR
# adaption of base implementation
function lmul_Q_adj!(A::LA.QR, b::AbstractVector)
    mA, nA = size(A.factors)
    mB = length(b)
    Afactors = A.factors
    @inbounds begin
        for k = 1:min(mA, nA)
            vBj = b[k]
            for i = k+1:mB
                vBj += conj(Afactors[i, k]) * b[i]
            end
            vBj = conj(A.τ[k]) * vBj
            b[k] -= vBj
            for i = k+1:mB
                b[i] -= Afactors[i, k] * vBj
            end
        end
    end
    b
end
function qr_ldiv!(x, QR::LA.QR, b::AbstractVector) where {T}
    # overwrites b
    # assumes QR is a tall matrix
    lmul_Q_adj!(QR, b)
    @inbounds for i = 1:length(x)
        x[i] = b[i]
    end
    ldiv_upper!(QR.factors, x)
    return x
end

function LA.ldiv!(x::AbstractVector, WS::MatrixWorkspace, b::AbstractVector)
    m, n = size(WS)
    if (m, n) == (1, 1)
        x[1] = b[1] / WS[1, 1]
        return x
    end
    WS.factorized[] || factorize!(WS)
    if m == n
        if WS.scaled[]
            x .= WS.row_scaling .* b
            lu_ldiv!(x, WS.lu, x)
        else
            lu_ldiv!(x, WS.lu, b)
        end
    else
        WS.r .= b
        qr_ldiv!(x, WS.qr, WS.r)
    end
    x
end

function LA.inv!(B::AbstractMatrix, M::MatrixWorkspace)
    M.factorized[] || factorize!(M)
    B .= M.lu.factors
    LA.inv!(M.lu)
    @inbounds for i in eachindex(B)
        M.lu.factors[i], B[i] = B[i], M.lu.factors[i]
    end
    return B
end

"""
    skeel_row_scaling!(W::MatrixWorkspace, c)
    skeel_row_scaling!(d, A, c)

Compute optimal scaling factors `d` for the matrix `A` following Skeel [^S79]
if `c` is approximately of the order of the solution of the linear system
of interest.
The scaling factors are rounded to powers of the base radix 2.

[^S79]: Skeel, Robert D. "Scaling for numerical stability in Gaussian elimination." Journal of the ACM (JACM) 26.3 (1979): 494-526.
"""
function skeel_row_scaling!(
    d::AbstractVector{<:Real},
    A::AbstractMatrix{<:Complex},
    c::AbstractVector{<:Real},
)
    n = length(c)
    @inbounds d .= zero(eltype(d))
    @inbounds for j = 1:n
        cj = c[j]
        for i = 1:n
            d[i] += fast_abs(A[i, j]) * cj
        end
    end
    @inbounds for i = 1:n
        d[i] = exp2(-last(frexp(d[i])))
    end

    d
end

function skeel_row_scaling!(W::MatrixWorkspace, c::AbstractVector{<:Real})
    skeel_row_scaling!(W.row_scaling, W.A, c)
    W
end

function row_scaling!(
    d::AbstractVector{<:Real},
    WS::MatrixWorkspace,
    c::AbstractVector{<:Real},
)
    m, n = size(WS)
    if m == n
        skeel_row_scaling!(d, WS.A, c)
    else
        d .= 1.0
    end
    d
end

"""
    apply_row_scaling!(W::MatrixWorkspace)

Apply the computed row scaling.
"""
function apply_row_scaling!(W::MatrixWorkspace)
    A, d = W.lu.factors, W.row_scaling
    m, n = size(A)
    @inbounds for j = 1:n, i = 1:m
        A[i, j] = A[i, j] * d[i]
    end
    W.scaled[] = true
    W
end

## Iterative Refinement
"""
    residual!(r, A, x, b)

Compute the residual `Ax-b` and store it in `r`
"""
function residual!(r::AbstractVector, A::AbstractMatrix, x::AbstractVector, b)
    m, n = size(A)
    @boundscheck m == length(b) == length(r) && n == length(x)
    r .= 0
    @inbounds for j = 1:n
        x_j = x[j]
        for i = 1:m
            r[i] += A[i, j] * x_j
        end
    end
    @inbounds for i = 1:m
        r[i] -= b[i]
    end
    r
end

"""
    mixed_precision_iterative_refinement!(x, A, b, norm = nothing)

Perform one step of mixed precision iterative refinement.
If `norm` is an `AbstractNorm` then the normwise relative error before
the refinement step is returned otherwise `x` is returned.
"""
function mixed_precision_iterative_refinement!(
    x::AbstractVector,
    M::MatrixWorkspace,
    b::AbstractVector,
    norm::Union{AbstractNorm,Nothing} = nothing,
)
    M.x̄ .= x
    residual!(M.r̄, M.A, M.x̄, b)
    M.r .= M.r̄
    LA.ldiv!(M.δx, M, M.r)
    x .-= M.δx
    if norm isa Nothing
        x
    else
        norm(M.δx) / norm(x)
    end
end

"""
    fixed_precision_iterative_refinement!(x, A, b, norm = nothing)

Perform one step of mixed precision iterative refinement.
If `norm` is an `AbstractNorm` then the normwise relative error before
the refinement step is returned otherwise `x` is returned.
"""
function fixed_precision_iterative_refinement!(
    x::AbstractVector,
    M::MatrixWorkspace,
    b::AbstractVector,
    norm::Union{AbstractNorm,Nothing} = nothing,
)
    residual!(M.r, M.A, x, b)
    LA.ldiv!(M.δx, M, M.r)
    x .-= M.δx
    if norm isa Nothing
        x
    else
        norm(M.δx) / norm(x)
    end
end

# CONDITION NUMBERS

"""
    inverse_inf_norm_est(WS::MatrixWorkspace,
                         d_l::Union{Nothing,Vector{<:Real}}=nothing
                         d_r::Union{Nothing,Vector{<:Real}}=nothing)

Estimate of the infinity norm of `diag(d_r)⁻¹A⁻¹diag(d_l)⁻¹` where `d_l` and `d_r` are
optional positive vectors.
If `d_l` or `d_r` is `nothing` the all one vector is used.
This uses the 1-norm lapack condition estimator described by Highahm in [^H88].

[^H88]: Higham, Nicholas J. "FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation." ACM Transactions on Mathematical Software (TOMS) 14.4 (1988): 381-396.
"""
function inverse_inf_norm_est(
    WS::MatrixWorkspace,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    WS.factorized[] || factorize!(WS)
    work, rwork = WS.inf_norm_est_work, WS.inf_norm_est_rwork
    if !WS.scaled[]
        inverse_inf_norm_est(WS.lu, d_l, d_r, nothing, work, rwork)
    else
        inverse_inf_norm_est(WS.lu, d_l, d_r, WS.row_scaling, work, rwork)
    end
end
function inverse_inf_norm_est(
    lu::LA.LU,
    d_l::Union{Nothing,Vector{<:Real}},
    d_r::Union{Nothing,Vector{<:Real}},
    row_scaling::Union{Nothing,Vector{<:Real}},
    work::Vector{<:Complex},
    rwork::Vector{<:Real},
)
    z = ξ = y = work
    x = rwork

    n = size(lu.factors, 1)
    x .= inv(n)
    if d_r !== nothing
        x ./= d_r
    end
    lu_ldiv_adj!(y, lu, x)
    if d_l !== nothing
        y ./= d_l
    end
    if row_scaling !== nothing
        y .*= row_scaling
    end
    γ = sum(fast_abs, y)
    ξ ./= fast_abs.(y)
    if d_l !== nothing
        ξ ./= d_l
    end
    if row_scaling !== nothing
        ξ ./= row_scaling
    end
    lu_ldiv!(z, lu, ξ)

    if d_r !== nothing
        x .= real.(z) ./ d_r
    else
        x .= real.(z)
    end
    k = 2
    while true
        j, max_xᵢ = 1, abs(x[1])
        for i = 2:n
            abs_xᵢ = abs(x[i])
            if abs_xᵢ > max_xᵢ
                j, max_xᵢ = i, abs_xᵢ
            end
        end
        x .= zero(eltype(x))
        x[j] = one(eltype(x))
        if d_r !== nothing
            x ./= d_r
        end
        lu_ldiv_adj!(y, lu, x)
        if d_l !== nothing
            y ./= d_l
        end
        if row_scaling !== nothing
            y .*= row_scaling
        end
        γ̄ = γ
        γ = sum(fast_abs, y)
        if γ ≤ γ̄
            γ = γ̄
            break
        end
        ξ .= y ./ fast_abs.(y)
        if d_l !== nothing
            ξ ./= d_l
        end
        if row_scaling !== nothing
            ξ ./= row_scaling
        end
        lu_ldiv!(z, lu, ξ)
        if d_r !== nothing
            x .= real.(z) ./ d_r
        else
            x .= real.(z)
        end
        k += 1
        if x[j] == LA.norm(x, Inf) || k > 2
            break
        end
    end
    nanmin(γ, Inf)
end

function inf_norm(
    WS::MatrixWorkspace,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    norm = -Inf
    A = WS.A
    n, m = size(A)
    @inbounds for i = 1:n
        normᵢ = 0.0
        for j = 1:m
            if d_r isa Nothing
                normᵢ += fast_abs(A[i, j])
            else
                normᵢ += fast_abs(A[i, j]) * d_r[j]
            end
        end
        if !(d_l isa Nothing)
            normᵢ *= d_l[i]
        end
        norm = @fastmath max(norm, normᵢ)
    end
    norm
end

function max_min_row(
    WS::MatrixWorkspace,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    max_row = -Inf
    min_row = Inf
    A = WS.A
    n, m = size(A)
    @inbounds for i = 1:n
        normᵢ = 0.0
        for j = 1:m
            if d_r isa Nothing
                normᵢ += fast_abs(A[i, j])
            else
                normᵢ += fast_abs(A[i, j]) * d_r[j]
            end
        end
        if !(d_l isa Nothing)
            normᵢ *= d_l[i]
        end
        max_row = @fastmath max(max_row, normᵢ)
        min_row = @fastmath min(min_row, normᵢ)
    end
    max_row, min_row
end

"""
    cond(A::MatrixWorkspace,
         d_l::Union{Nothing,Vector{<:Real}} = nothing,
         d_r::Union{Nothing,Vector{<:Real}} = nothing)

Compute the condition number w.r.t. the infinity-norm of `diag(d_l) * A * diag(d_r)`.
If `d_l` or `d_r` is `nothing` the all one vector is used.
If `size(A) == (1,1)` then just the norm of the inverse is returned.
"""
function LA.cond(
    WS::MatrixWorkspace,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    m, n = size(WS)
    if m == n == 1
        if isa(d_l, Nothing) && isa(d_r, Nothing)
            inv(abs(WS.A[1, 1]))
        elseif isa(d_l, Nothing)
            inv(abs(WS.A[1, 1]) * d_r[1])
        elseif isa(d_r, Nothing)
            inv(d_l[1] * abs(WS.A[1, 1]))
        else
            inv(d_l[1] * abs(WS.A[1, 1]) * d_r[1])
        end
    elseif m > n
        WS.factorized[] || factorize!(WS)
        rmax, rmin = -Inf, Inf
        for i = 1:n
            rᵢ = fast_abs(WS.qr.factors[i, i]) * d_r[i]
            rmax = max(rmax, rᵢ)
            rmin = min(rmin, rᵢ)
        end
        rmax / rmin
    else
        inverse_inf_norm_est(WS, d_l, d_r) * inf_norm(WS, d_l, d_r)
    end
end

function egcond(
    WS::MatrixWorkspace,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    m, n = size(WS)
    if m == n == 1
        a = abs(WS.A[1, 1])
        if d_l !== nothing
            a *= d_l[1]
        end
        if d_r !== nothing
            a *= d_r[1]
        end
        max(a, 1.0) / a
    elseif m > n
        WS.factorized[] || factorize!(WS)
        rmax, rmin = -Inf, Inf
        for i = 1:n
            rᵢ = abs(WS.qr.factors[i, i]) * d_r[i]
            rmax = max(rmax, rᵢ)
            rmin = min(rmin, rᵢ)
        end
        rmax / rmin
    else
        max_row, min_row = max_min_row(WS, d_l, d_r)
        inverse_inf_norm_est(WS, d_l, d_r) * max(max_row, 1 / min_row)
    end
end

#####################
## Jacobian ##
#####################
struct Jacobian{M}
    workspace::MatrixWorkspace{M}
    # stats
    factorizations::Base.RefValue{Int}
    ldivs::Base.RefValue{Int}
end
Jacobian(A::AbstractMatrix) = Jacobian(MatrixWorkspace(A), Ref(0), Ref(0))

updated!(J::Jacobian) = (updated!(J.workspace); J)

workspace(J::Jacobian) = J.workspace
matrix(J::Jacobian) = matrix(workspace(J))

function Base.show(io::IO, J::Jacobian{T}) where {T}
    println(io, "Jacobian{$T}:")
    println(io, " • # factorizations → ", J.factorizations[])
    println(io, " • # ldivs → ", J.ldivs[])
end

"""
    ldiv!(x̂::AbstractVector, J::Jacobian, b::AbstractVector)

solve the linear system `matrix(J)x̂=b`.
"""
function LA.ldiv!(
    x̂::AbstractVector,
    J::Jacobian,
    b::AbstractVector,
    norm = nothing;
    row_scaling::Bool = true,
)
    # stats update
    J.factorizations[] += !workspace(J).factorized[]
    J.ldivs[] += 1
    LA.ldiv!(x̂, workspace(J), b)
    x̂
end
function LA.ldiv!(
    x̂::AbstractVector,
    J::Jacobian,
    b::AbstractVector,
    norm::WeightedNorm;
    row_scaling::Bool = true,
)
    if !workspace(J).factorized[] && row_scaling
        skeel_row_scaling!(workspace(J), weights(norm))
        apply_row_scaling!(workspace(J))
    end
    # stats update
    J.factorizations[] += !workspace(J).factorized[]
    J.ldivs[] += 1
    LA.ldiv!(x̂, workspace(J), b)
    x̂
end

function iterative_refinement!(
    x::AbstractVector,
    J::Jacobian,
    b::AbstractVector,
    norm::AbstractNorm;
    max_iters::Int = 3,
    tol::Float64 = sqrt(eps()),
)
    J.ldivs[] += 1
    δ = mixed_precision_iterative_refinement!(x, workspace(J), b, norm)
    for i = 2:max_iters
        J.ldivs[] += 1
        δ′ = mixed_precision_iterative_refinement!(x, workspace(J), b, norm)
        if δ′ < tol
            return (accuracy = δ′, diverged = false)
        elseif δ′ > 0.5 * δ
            return (accuracy = δ′, diverged = true)
        end
        δ = δ′
    end
    (accuracy = δ, diverged = false)
end

"""
    init!(J::Jacobian)

(Re-)initialize the `Jacobian`.
"""
function init!(J::Jacobian; keep_stats::Bool = false)
    J.factorizations[] = 0
    J.ldivs[] = 0
    J
end

function LA.cond(
    J::Jacobian,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    LA.cond(workspace(J), d_l, d_r)
end

function egcond(
    J::Jacobian,
    d_l::Union{Nothing,Vector{<:Real}} = nothing,
    d_r::Union{Nothing,Vector{<:Real}} = nothing,
)
    egcond(workspace(J), d_l, d_r)
end


text:
@testset "Configuration" begin
    @testset "constructors" begin
        c = Configuration("test")
        p = [StringParameter("value",  "test1"),
             IntegerParameter(1, 5, 4, "test2")]
        c = Configuration(p, "test")
        @test c["test1"] == p[1]
        @test c["test2"] == p[2]
        p = Dict{String, Parameter}()
        p["test1"] = StringParameter("value",  "test")
        p["test2"] = IntegerParameter(1, 5, 4, "test")
        c = Configuration(p, "test")
        @test c["test1"] == p["test1"]
        @test c["test2"] == p["test2"]
        p = [StringParameter("valuea",  "test1"),
             StringParameter("valueb", "test2"),
             IntegerParameter(1, 44, 3, "test3")]
        c = Configuration(p, "test")
        @test c["test1"] == p[1]
        @test c["test2"] == p[2]
    end
    @testset "perturb!" begin
        p = [StringParameter("valuea",  "test1"),
             StringParameter("valueb", "test2"),
             IntegerParameter(1, 44, 3, "test3")]
        t = Configuration(p, "test1")
        l = [IntegerParameter(1, 38, 3, "i0"),
             IntegerParameter(4, 66, 55, "i1"),
             StringParameter("value", "i2"),
             FloatParameter(2.33, 85.33, 22.2, "i3"),
             t]
        c = Configuration(l, "test")
        r = perturb!(c)
        r = perturb!(c)
        r = perturb!(c)
        r = perturb!(c)
        @test typeof(r["test1"]) <: Configuration
        @test r["test1"]["test1"].value == "valuea"
        @test r["test1"]["test2"].value == "valueb"
        @test r["i2"].value             == "value"
        @test typeof(r["test1"]["test3"]) <: NumberParameter
        @test typeof(r["i0"])             <: NumberParameter
        @test typeof(r["i1"])             <: NumberParameter
        @test typeof(r["i3"])             <: NumberParameter
        d = Dict{String, Any}()
        d["i0"]                    = 3
        d["i1"]                    = 4
        d["i3"]                    = 2.3
        d["test1"]                 = Dict{String, Any}()
        d["test1"]["test3"]        = 2
        @test (typeof(perturb!(c, d)) <: Configuration)
        r                          = perturb!(c, d)
        r                          = perturb!(c, d)
        r                          = perturb!(c, d)
        r                          = perturb!(c, d)
        @test r["i0"].value       == c["i0"].value
        @test r["i1"].value       == c["i1"].value
        @test r["i3"].value       == c["i3"].value
        @test r["test1"]["test3"] == c["test1"]["test3"]
        d["test1"]["test1"]        = 23
        @test_throws MethodError perturb!(c, d)
    end
    @testset "neighbor!" begin
        l = [IntegerParameter(1, 38, 3, "i0"),
             IntegerParameter(4, 66, 55, "i1"),
             StringParameter("value", "i2"),
             FloatParameter(2.33, 85.33, 22.2, "i3")]
        d        = Dict{String, Any}()
        c        = Configuration(l, "test")
        println(c)
        d["i0"]  = 5
        v        = c["i0"].value
        neighbor!(c, d)
        @test v != c["i0"].value
        v        = c["i0"].value
        neighbor!(c, d, distance = 9)
        d        = Dict{String, Any}()
        d["i2"]  = 5
        @test_throws MethodError neighbor!(c, d)
        @test_throws MethodError neighbor!(c, d, 2)
        l        = [IntegerParameter(1, 38, 3, "i0")]
        c        = Configuration(l, "test2")
        v        = c["i0"].value
        neighbor!(c)
        @test v != c["i0"].value
    end
    @testset "update! and convert" begin
        c                    = Configuration("test")
        c["test"]            = IntegerParameter(1, 2323, 2, "i0")
        a                    = convert(Array{Parameter}, c)
        @test a[1]          == c["test"]
        push!(a, StringParameter("name", "i2"))
        update!(c, a)
        @test c["i2"]       != a[1]
        @test c["i2"].value == a[2].value
    end
end


text:
struct Unsubscribe <: HasId
    header::UInt8
    id::UInt16
    topics::Array{String}
end
Unsubscribe(topics::Array{String}) = Unsubscribe(convert(UInt8, UNSUBSCRIBE) | 0x02, 0x0000, topics)
Unsubscribe(packet::Unsubscribe, id::UInt16) = Unsubscribe(packet.header, id, packet.topics)

function write(s::IO, packet::Unsubscribe)
    mqtt_write(s, packet.id)
    for topic in packet.topics
        mqtt_write(s, topic)
    end
end

Base.show(io::IO, x::Unsubscribe) = print(io, "unsubscribe[id: ", x.id, ", topics: ", join(x.topics, ", "), "]")

struct Unsuback <: Ack
    header::UInt8
    id::UInt16
end
Unsuback(id) = Unsuback(UNSUBACK, id)
Base.show(io::IO, x::Unsuback) = print(io, "unsuback[id: ", x.id, "]")

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
# Combined 2D plots across multiple frames, multi-processing version.
#
# To run on a single node,
# julia -p $ncores demo_2d_mp_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed, ParallelDataTransfer, Glob
@everywhere using Vlasiator, PyPlot, Printf, LaTeXStrings

@everywhere function init_figure(cmaps, norms, ticks, pArgs1, extent)
   fig, axs = plt.subplots(2, 3; num=myid(),
      figsize=(13.5, 9.5), sharex=true, sharey=true, constrained_layout=true)

   for ax in axs
      ax.set_aspect("equal")
      # Set border line widths
      for loc in ("left", "bottom", "right", "top")
         edge = get(ax.spines, loc, nothing)
         edge.set_linewidth(2.0)
      end
      ax.xaxis.set_tick_params(width=2.0, length=3)
      ax.yaxis.set_tick_params(width=2.0, length=3)
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   end

   axs[2,1].set_xlabel(pArgs1.strx; fontsize)
   axs[2,2].set_xlabel(pArgs1.strx; fontsize)
   axs[2,3].set_xlabel(pArgs1.strx; fontsize)
   axs[1,1].set_ylabel(pArgs1.stry; fontsize)
   axs[2,1].set_ylabel(pArgs1.stry; fontsize)

   axs[1,1].set_title("Density"; fontsize)
   axs[1,2].set_title("Velocity x"; fontsize)
   axs[1,3].set_title("Velocity z"; fontsize)
   axs[2,1].set_title("Thermal pressure"; fontsize)
   axs[2,2].set_title("Magnetic field z"; fontsize)
   axs[2,3].set_title("Electric field"; fontsize)

   x1, x2 = Vlasiator.get_axis(pArgs1)

   range1 = searchsortedfirst(x1, extent[1]):searchsortedlast(x1, extent[2])
   range2 = searchsortedfirst(x2, extent[3]):searchsortedlast(x2, extent[4])

   fakedata = zeros(Float32, length(range2), length(range1))

   c1 = axs[1,1].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[1], cmap=cmaps[1])
   c2 = axs[1,2].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[2], cmap=cmaps[2])
   c3 = axs[1,3].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[3], cmap=cmaps[2])
   c4 = axs[2,1].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[4], cmap=cmaps[1])
   c5 = axs[2,2].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[5], cmap=cmaps[2])
   c6 = axs[2,3].pcolormesh(x1[range1], x2[range2], fakedata; norm=norms[6], cmap=cmaps[1])

   format = matplotlib.ticker.FormatStrFormatter("%.1f")
   cb1 = colorbar(c1; ax=axs[1,1], ticks=ticks[1], format)
   cb1.ax.set_ylabel("[amu/cc]"; fontsize)
   cb1.outline.set_linewidth(1.0)

   cb2 = colorbar(c2; ax=axs[1,2], ticks=ticks[2])
   cb2.ax.set_ylabel("[km/s]"; fontsize)
   cb2.outline.set_linewidth(1.0)

   cb3 = colorbar(c3; ax=axs[1,3], ticks=ticks[3])
   cb3.ax.set_ylabel("[km/s]"; fontsize)
   cb3.outline.set_linewidth(1.0)

   cb4 = colorbar(c4; ax=axs[2,1], ticks=ticks[4], format)
   cb4.ax.set_ylabel("[nPa]"; fontsize)
   cb4.outline.set_linewidth(1.0)

   cb5 = colorbar(c5; ax=axs[2,2], ticks=ticks[5], format, extend="both")
   cb5.ax.set_ylabel("[nT]"; fontsize)
   cb5.outline.set_linewidth(1.0)

   cb6 = colorbar(c6; ax=axs[2,3], ticks=ticks[6], extend="max")
   cb6.ax.set_ylabel("[mV/m]"; fontsize)
   cb6.outline.set_linewidth(1.0)

   cs = (c1, c2, c3, c4, c5, c6)

   return fig, cs, range1, range2
end

@everywhere function update_plot!(fig, cs, range1, range2, outdir, file)
   isfile(outdir*file[end-8:end-5]*".png") && return

   println("file = $file")
   meta = load(file)

   data = Vlasiator.prep2d(meta, "proton/vg_rho", :mag)'
   cs[1].set_array(data[range2,range1] ./ 1f6)

   data = Vlasiator.prep2d(meta, "proton/vg_v", :x)'
   cs[2].set_array(data[range2,range1] ./ 1f3)

   data = Vlasiator.prep2d(meta, "proton/vg_v", :z)'
   cs[3].set_array(data[range2,range1] ./ 1f3)

   data = Vlasiator.prep2d(meta, "vg_pressure", :mag)'
   cs[4].set_array(data[range2,range1] .* 1f9)

   data = Vlasiator.prep2d(meta, "vg_b_vol", :z)'
   cs[5].set_array(data[range2,range1] .* 1f9)

   data = Vlasiator.prep2d(meta, "vg_e_vol", :mag)'
   cs[6].set_array(data[range2,range1] .* 1f3)

   str_title = @sprintf "Density pulse run, t= %4.1fs" meta.time
   fig.suptitle(str_title; fontsize="xx-large")

   savefig(outdir*file[end-8:end-5]*".png", bbox_inches="tight")
end

function make_jobs(files)
   for f in files
       put!(jobs, f)
   end
end

@everywhere function do_work(jobs, status, outdir, cmaps, norms, ticks, pArgs1, extent)
   fig, cs, range1, range2 = init_figure(cmaps, norms, ticks, pArgs1, extent)
   while true
      file = take!(jobs)
      update_plot!(fig, cs, range1, range2, outdir, file)
      put!(status, true)
   end
   close(fig)
end

############################################################################################
files = glob("bulk*.vlsv", ".")

nfile = length(files)
# Set output directory
outdir = "contour/"

# Set colormaps for continuous and divergent data
cmaps = (matplotlib.cm.turbo, matplotlib.cm.RdBu_r)
axisunit = EARTH
extent = [0., 20., -20., 20.] # [RE], default full domain: [-Inf32, Inf32, -Inf32, Inf32]

# Plotting range for each variable
ρmin, ρmax   = 0.1, 15.0     # [amu/cc]
vxmin, vxmax = -650.0, 650.0 # [km/s]
vzmin, vzmax = -500.0, 500.0 # [km/s]
pmin, pmax   = 0.0, 3.6      # [nPa]
bmin, bmax   = -60.0, 60.    # [nT]
emin, emax   = 0.0, 20.      # [mV/m]

meta = load(files[1])

pArgs1 = Vlasiator.set_args(meta, "proton/vg_rho", axisunit; normal=:none)

norm1, ticks1 = Vlasiator.set_colorbar(Linear, ρmin, ρmax)
norm2, ticks2 = Vlasiator.set_colorbar(Linear, vxmin, vxmax)
norm3, ticks3 = Vlasiator.set_colorbar(Linear, vzmin, vzmax)
norm4, ticks4 = Vlasiator.set_colorbar(Linear, pmin, pmax)
norm5, ticks5 = Vlasiator.set_colorbar(Linear, bmin, bmax)
norm6, ticks6 = Vlasiator.set_colorbar(Linear, emin, emax)

norms = (norm1, norm2, norm3, norm4, norm5, norm6)
ticks = (ticks1, ticks2, ticks3, ticks4, ticks5, ticks6)

const jobs   = RemoteChannel(()->Channel{String}(nfile))
const status = RemoteChannel(()->Channel{Bool}(nworkers()))

@broadcast begin
   const fontsize = "x-large"
end

println("Total number of files: $nfile")
println("Running with $(nworkers()) workers...")

@async make_jobs(files) # Feed the jobs channel with all files to process.

@sync for p in workers()
   @async remote_do(do_work, p, jobs, status, outdir, cmaps, norms, ticks, pArgs1, extent)
end

let n = nfile
   t = @elapsed while n > 0 # wait for all jobs to complete
      take!(status)
      n -= 1
   end
   println("Finished in $(round(t, digits=2))s.")
end


text:
using Distributed
using Distributions
using Genie, Genie.Router, Genie.Renderer, Genie.Requests
import Genie.Router: @params

form = """

Approximation of Pi

 Enter an integer in the form, I will calculate an approximation of Pi.

<form action="/" method="POST" enctype="multipart/form-data">
  <input type="text" name="nb_draws" value="" placeholder="Enter number of random draws" />
  <input type="submit" value="Enter" />
</form>
"""

route("/") do
  html(form)
end

route("/", method = POST) do

	function inside_circle(x::Float64, y::Float64)
	    output = 0
	    if x^2 + y^2 <= 1
		output = 1
	    end
	    return output
	end

	function pi_serial(nbPoints::Int64 = 128 * 1000; d=Uniform(-1.0,1.0))

	   #draw NbPoints from within the square centered in 0
	   #with side length equal to 2
	   xDraws = rand(d, nbPoints)
	   yDraws = rand(d, nbPoints)
	   sumInCircle = 0

	   for (xValue, yValue) in zip(xDraws, yDraws)
		sumInCircle+=inside_circle(xValue, yValue)
	   end

	   return 4*sumInCircle/nbPoints

	end

	pi_approx = pi_serial(parse(Int64, @params(:nb_draws)));

	"""

Approximation of Pi

Approximation for Pi $(pi_approx) with $(parse(Int64, @params(:nb_draws))) random draws
 True value is $(pi)"""

end

up()

max_stars_repo_name: JulienPascal/MyFirstDockerApp

max_issues_repo_name: JulienPascal/MyFirstDockerApp

max_forks_repo_name: JulienPascal/MyFirstDockerApp


text:
using DataFrames, FixedEffectModels, Random

# Very simple setup
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
@time reg(df, @formula(y ~ x1 + x2))
# 0.582029 seconds (852 allocations: 535.311 MiB, 18.28% gc time)
@time reg(df, @formula(y ~ x1 + x2),  Vcov.cluster(:id2))
# 0.809652 seconds (1.55 k allocations: 649.821 MiB, 14.40% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)))
# 1.655732 seconds (1.21 k allocations: 734.353 MiB, 16.88% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)), Vcov.cluster(:id1))
#  2.113248 seconds (499.66 k allocations: 811.751 MiB, 15.08% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
# 3.553678 seconds (1.36 k allocations: 1005.101 MiB, 10.55% gc time))

# More complicated setup
N = 800000 # number of observations
M = 40000 # number of workers
O = 5000 # number of firms
id1 = rand(1:M, N)
id2 = [rand(max(1, div(x, 8)-10):min(O, div(x, 8)+10)) for x in id1]
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
# 15.738214 seconds (929.97 k allocations: 147.737 MiB, 0.95% gc time)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ a9bbd3ea-fc34-11ea-282a-43f4fc159d39
using Pkg, DrWatson

# ╔═╡ a9e16b46-fc34-11ea-36de-ef647e4e4f6f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanQuap
	using StatisticalRethinking
end

# ╔═╡ feebbe4e-fc33-11ea-2ee7-d14f977c6497
md"## Clip-05-01-02s.jl"

# ╔═╡ a9e54766-fc34-11ea-1246-93a25f748d1d
md"### snippet 5.1"

# ╔═╡ 10471c94-fc55-11ea-0fec-656bd563d513
md"##### D (Divorce rate), A (MediumAgeMarriage) and M (Marriage rate) are all standardized."

# ╔═╡ a9f32456-fc34-11ea-3af8-c7e0de78a85d
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
	scale!(df, [:Divorce, :MedianAgeMarriage, :Marriage])
end;

# ╔═╡ 4d0ca900-fc53-11ea-0f6c-e7d09d775387
md"##### The model m5.1s represents a regression of Divorce on MedianAgeMarriage and is defined as:"

# ╔═╡ 25c70d5a-fc54-11ea-3910-a9276dc7b696
md"
```
model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         // Priors
  bA ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bA * A;
  D ~ normal(mu , sigma);     // Likelihood
}
```
"

# ╔═╡ b8ae47cc-01af-11eb-3cbd-a749b6e3d327
md"##### The Stan language model."

# ╔═╡ b5431248-01af-11eb-03ac-ed6a17e17e39
stan5_1_priors = "
	data {
	 int N; // Sample size
	 vector[N] D; // Outcome
	 vector[N] A; // Predictor
	}

	parameters {
	 real a; // Intercept
	 real bA; // Slope (regression coefficients)
	 real < lower = 0 > sigma;    // Error SD
	}

	model {
	  vector[N] mu;               // mu is a vector
	  a ~ normal(0, 0.2);         // Priors
	  bA ~ normal(0, 0.5);
	  sigma ~ exponential(1);
	  mu = a + bA * A;
	  D ~ normal(mu , sigma);   // Likelihood
	}
";

# ╔═╡ db59a030-423e-11eb-3eb2-5d1f7d3dd0f5
begin
	m5_1s_priors = SampleModel("m5.1.priors", stan5_1_priors)
	rc5_1s_priors = stan_sample(m5_1s_priors;
		data = Dict("N" => 0, "D" => [], "A" => []))
	success(rc5_1s_priors) &&
		(part5_1s_priors = read_samples(m5_1s_priors, :particles))
end

# ╔═╡ 55ef3d46-423f-11eb-1d83-55276a40b702
if success(rc5_1s_priors)
	priors5_1s_df = read_samples(m5_1s_priors, :dataframe)
	xi = -3.0:0.1:3.0
	plot(xlab="Medium age marriage (scaled)", ylab="Divorce rate (scaled)",
		title="Showing 50 regression lines")
	for i in 1:50
		local yi = mean(priors5_1s_df[i, :a]) .+ priors5_1s_df[i, :bA] .* xi
		plot!(xi, yi, color=:lightgrey, leg=false)
	end
	scatter!(df[:, :MedianAgeMarriage_s], df[!, :Divorce_s], color=:darkblue)
end

# ╔═╡ fad134cc-423f-11eb-3398-b7d0e22aad9c
stan5_1 = "
	data {
	 int N; // Sample size
	 vector[N] D; // Outcome
	 vector[N] A; // Predictor
	}

	parameters {
	 real a; // Intercept
	 real bA; // Slope (regression coefficients)
	 real < lower = 0 > sigma;    // Error SD
	}

	model {
	  vector[N] mu;               // mu is a vector
	  a ~ normal(0, 0.2);         // Priors
	  bA ~ normal(0, 0.5);
	  sigma ~ exponential(1);
	  mu = a + bA * A;
	  D ~ normal(mu , sigma);   // Likelihood
	}
";

# ╔═╡ aa241f5e-fc34-11ea-21d5-1b03dfafd34f
md"##### Result rethinking:"

# ╔═╡ aa2a005e-fc34-11ea-09ab-290a8d6b6a9a
rethinking = "
	   mean   sd  5.5% 94.5%
a      0.00 0.10 -0.16  0.16
bA    -0.57 0.11 -0.74 -0.39
sigma  0.79 0.08  0.66  0.91
";

# ╔═╡ cbd5fe10-01af-11eb-055b-bd009d0dca55
begin
	data = (N = size(df, 1), D = df.Divorce_s, A = df.MedianAgeMarriage_s,
		M = df.Marriage_s)
	init = (a=1.0, bA=1.0, bM=1.0, sigma=10.0)
	q5_1s, m5_1s, om5_1s = stan_quap("m5.1s", stan5_1; data, init)
	if !isnothing(q5_1s)
		quap5_1s_df = sample(q5_1s)
		PRECIS(quap5_1s_df)
	end
end

# ╔═╡ 9547a2a0-01c4-11eb-169c-17d6dbace4d9
md"##### Compare below figure with the corresponding figure in clip-05-01-05s.jl."

# ╔═╡ 07cfcaec-01c4-11eb-38b7-31dca6cafefb
if !isnothing(m5_1s)
	post5_1s_df = read_samples(m5_1s, :dataframe)
	plot(xlab="Medium age marriage (scaled)", ylab="Divorce rate (scaled)",
		title="Showing 50 regression lines")
	for i in 1:50
		local yi = mean(post5_1s_df[i, :a]) .+ post5_1s_df[i, :bA] .* xi
		plot!(xi, yi, color=:lightgrey, leg=false)
	end
	scatter!(df[:, :MedianAgeMarriage_s], df[!, :Divorce_s], color=:darkblue)
end

# ╔═╡ aa064588-fc34-11ea-2153-59bc6c7f18a8
md"### snippet 5.2"

# ╔═╡ aa36edc8-fc34-11ea-3701-434917b6f7a3
if !isnothing(q5_1s)

	# Plot regression line D on A

	title1 = "Divorce rate vs. median age at marriage" *
		"\nshowing predicted and quantile range"
	fig1 = plotbounds(
		df, :MedianAgeMarriage, :Divorce,
		quap5_1s_df, [:a, :bA, :sigma];
		title=title1,
		colors=[:lightblue, :darkgrey]
	)
end

# ╔═╡ 683093dc-fc54-11ea-3be9-fdad0a8812f6
md"##### The model m5.2s represents a regression of Divorce on Marriage and is defined as:"

# ╔═╡ 77cb9670-fc54-11ea-38fe-3b8c48e2ce09
md"
```
model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         // Priors
  bM ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bM * M;
  D ~ normal(mu , sigma);     // Likelihood
}
```
"

# ╔═╡ 3c3888ec-fc55-11ea-3542-c7df3ecc74b4
md"##### Both D (Divorce rate) and A (Marriage rate) are standardized."

# ╔═╡ 864e52fc-fc39-11ea-2c30-43e6b888bf60
stan5_2 = "
data {
  int N;
  vector[N] D;
  vector[N] M;
}
parameters {
  real a;
  real bM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + bM * M;
  a ~ normal( 0 , 0.2 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ normal( mu , sigma );
}
";

# ╔═╡ 4a4b1284-80f7-11eb-01a4-452b4e2597e7
begin
	q5_2s, m5_2s, om = stan_quap("m5.2s", stan5_2; data)
	if !isnothing(q5_2s)
		quap5_2s_df = sample(q5_2s)
	end
	PRECIS(quap5_2s_df)
end

# ╔═╡ c206cdc6-fc55-11ea-2835-71c3067485f8
if !isnothing(m5_2s)
	post5_2s_df = read_samples(m5_2s, :dataframe)
	PRECIS(post5_2s_df)
end

# ╔═╡ 279006d8-fc56-11ea-2792-d1c2678a1e08
if !isnothing(q5_2s)

	# Plot regression line D on M

	title2 = "Divorce rate vs. marriage rate" * "\nshowing predicted and hpdi range"
	fig2 = plotbounds(
		df, :Marriage, :Divorce,
		quap5_2s_df, [:a, :bM, :sigma];
		title=title2,
		colors=[:lightblue, :darkgrey]
	)

end

# ╔═╡ 41bc0716-fc56-11ea-1cfe-3db82349a2d2
	plot(fig2, fig1, layout=(1,2), title="")

# ╔═╡ aa37b78a-fc34-11ea-11e5-7d1ef7bdf603
md"## End of clip-05-01-02s.jl"

# ╔═╡ Cell order:
# ╟─feebbe4e-fc33-11ea-2ee7-d14f977c6497
# ╠═a9bbd3ea-fc34-11ea-282a-43f4fc159d39
# ╠═a9e16b46-fc34-11ea-36de-ef647e4e4f6f
# ╟─a9e54766-fc34-11ea-1246-93a25f748d1d
# ╟─10471c94-fc55-11ea-0fec-656bd563d513
# ╠═a9f32456-fc34-11ea-3af8-c7e0de78a85d
# ╟─4d0ca900-fc53-11ea-0f6c-e7d09d775387
# ╟─25c70d5a-fc54-11ea-3910-a9276dc7b696
# ╟─b8ae47cc-01af-11eb-3cbd-a749b6e3d327
# ╠═b5431248-01af-11eb-03ac-ed6a17e17e39
# ╠═db59a030-423e-11eb-3eb2-5d1f7d3dd0f5
# ╠═55ef3d46-423f-11eb-1d83-55276a40b702
# ╠═fad134cc-423f-11eb-3398-b7d0e22aad9c
# ╟─aa241f5e-fc34-11ea-21d5-1b03dfafd34f
# ╠═aa2a005e-fc34-11ea-09ab-290a8d6b6a9a
# ╠═cbd5fe10-01af-11eb-055b-bd009d0dca55
# ╟─9547a2a0-01c4-11eb-169c-17d6dbace4d9
# ╠═07cfcaec-01c4-11eb-38b7-31dca6cafefb
# ╟─aa064588-fc34-11ea-2153-59bc6c7f18a8
# ╠═aa36edc8-fc34-11ea-3701-434917b6f7a3
# ╟─683093dc-fc54-11ea-3be9-fdad0a8812f6
# ╟─77cb9670-fc54-11ea-38fe-3b8c48e2ce09
# ╟─3c3888ec-fc55-11ea-3542-c7df3ecc74b4
# ╠═864e52fc-fc39-11ea-2c30-43e6b888bf60
# ╠═4a4b1284-80f7-11eb-01a4-452b4e2597e7
# ╠═c206cdc6-fc55-11ea-2835-71c3067485f8
# ╠═279006d8-fc56-11ea-2792-d1c2678a1e08
# ╠═41bc0716-fc56-11ea-1cfe-3db82349a2d2
# ╟─aa37b78a-fc34-11ea-11e5-7d1ef7bdf603


text:
# CR3BP Indirect with STM parameters
mutable struct CR3BPIndirectWithSTMParams <: AbstractCR3BPIndirectParams
    # Spacecraft parameters
    sp::SimpleSpacecraft

    # CR3BP parameters
    crp::CR3BPParams

    # Continuation variable
    ϵ::Float64

    # Thrust type
    utype::Int64

    # Allocate matricies
    m1::Matrix{Float64}
    m2::Matrix{Float64}

    function CR3BPIndirectWithSTMParams(sp::Spacecraft, crp::CR3BPParams, ϵ)
        new(sp, crp, ϵ, 0, zeros(14,14), zeros(14,14))
    end
end

function cr3bpEomIndirectWithSTM!(du::AbstractVector, u::AbstractVector,
                                  p::CR3BPIndirectWithSTMParams, t,
                                  homotopyFlag::MEMF)
    @inbounds begin
        # Evaluate state/co-state dynamics
        GVec = cr3bpEomIndirect!(du, u, p, t, homotopyFlag)

        # Dynamics Jacobian
        jac = p.m1; jac .= 0.0
        cr3bpEOMJac!(jac, u, p, GVec, homotopyFlag)

        # Compute STM derivs
        Φ   = reshape(view(u,  15:210), (14,14))
        dΦ  = reshape(view(du, 15:210), (14,14))
        Octavian.matmul!(dΦ, jac, Φ, 1.0, 0.0)
    end
    return nothing
end

function cr3bpEOMJac!(jac, u, p, GVec, homotopyFlag::MEMF)

    # Get requirements
    TU  = p.crp.TU
    LU  = p.crp.LU
    MU  = p.crp.MU
    μ   = p.crp.μ
    c   = p.sp.c

    # Scale requirements
    tMaxSc = p.sp.tMax * TU * TU / (MU*LU*1000.0)
    cSc = c*TU / (LU*1000.0)

    # Compute thrust direction
    λv = sqrt(u[11]*u[11] + u[12]*u[12] + u[13]*u[13])

    # Compute throttline
    S = computeS(u, λv, cSc)
    γ = computeU(S, p.utype, p.ϵ)

    # Compute r1 and r2 with inverse powers
    xpmu    = u[1] + μ
    xpmum1  = u[1] + μ - 1
    r1      = sqrt(xpmu*xpmu + u[2]*u[2] + u[3]*u[3])
    r2      = sqrt(xpmum1*xpmum1 + u[2]*u[2] + u[3]*u[3])
    invr1   = inv(r1)
    invr2   = inv(r2)
    invr15  = invr1*invr1*invr1*invr1*invr1
    invr25  = invr2*invr2*invr2*invr2*invr2
    invr17  = invr15*invr1*invr1
    invr27  = invr25*invr2*invr2

    # Partials of G matrix wrt r
    dG11dr          = @SVector [(1 - μ)*(μ + u[1])*(9*invr15 - 15*(u[1] + μ)^2*invr17) +
                    μ*(μ + u[1] - 1)*(9*invr25 - 15*(u[1] + μ - 1)^2*invr27),
                    (1 - μ)*u[2]*(3*invr15 - 15*(u[1] + μ)^2*invr17) +
                    μ*u[2]*(3*invr25 - 15*(u[1] + μ - 1)^2*invr27),
                    (1 - μ)*u[3]*(3*invr15 - 15*(u[1] + μ)^2*invr17) +
                    μ*u[3]*(3*invr25 - 15*(u[1] + μ - 1)^2*invr27)];

    dG12dr          = @SVector [u[2]*(1 - μ)*(3*invr15 - 15*(u[1] + μ)^2*invr17) +
                    u[2]*μ*(3*invr25 - 15*(u[1] + μ - 1)^2*invr27),
                    (1 - μ)*(u[1] + μ)*(3*invr15 - 15*u[2]^2*invr17) +
                    μ*(u[1] + μ - 1)*(3*invr25 - 15*u[2]^2*invr27),
                    -15*(1 - μ)*(u[1] + μ)*u[2]*u[3]*invr17 -
                        15*μ*(u[1] + μ - 1)*u[2]*u[3]*invr27];

    dG13dr          = @SVector [u[3]*(1 - μ)*(3*invr15 - 15*(u[1] + μ)^2*invr17) +
                    u[3]*μ*(3*invr25 - 15*(u[1] + μ - 1)^2*invr27),
                    -15*(1 - μ)*(u[1] + μ)*u[2]*u[3]*invr17 -
                    15*μ*(u[1] + μ - 1)*u[2]*u[3]*invr27,
                    (1 - μ)*(u[1] + μ)*(3*invr15 - 15*u[3]^2*invr17) +
                    μ*(u[1] + μ - 1)*(3*invr25 - 15*u[3]^2*invr27)];

    dG22dr          = @SVector [(1 - μ)*(u[1] + μ)*(3*invr15 - 15*u[2]^2*invr17) +
                    μ*(u[1] + μ - 1)*(3*invr25 - 15*u[2]^2*invr27),
                    (1 - μ)*u[2]*(9*invr15 - 15*u[2]^2*invr17) +
                    μ*u[2]*(9*invr25 - 15*u[2]^2*invr27),
                    (1 - μ)*u[3]*(3*invr15 - 15*u[2]^2*invr17) +
                    μ*u[3]*(3*invr25 - 15*u[2]^2*invr27)];

    dG23dr          = @SVector [-15*(1 - μ)*(u[1] + μ)*u[2]*u[3]*invr17 -
                    15*μ*(u[1] + μ - 1)*u[2]*u[3]*invr27,
                    (1 - μ)*u[3]*(3*invr15 - 15*u[2]^2*invr17) +
                    μ*u[3]*(3*invr25 - 15*u[2]^2*invr27),
                    (1 - μ)*u[2]*(3*invr15 - 15*u[3]^2*invr17) +
                    μ*u[2]*(3*invr25 - 15*u[3]^2*invr27)];

    dG33dr          = @SVector [(1 - μ)*(u[1] + μ)*(3*invr15 - 15*u[3]^2*invr17) +
                    μ*(u[1] + μ - 1)*(3*invr25 - 15*u[3]^2*invr27),
                    (1 - μ)*u[2]*(3*invr15 - 15*u[3]^2*invr17) +
                    μ*u[2]*(3*invr25 - 15*u[3]^2*invr27),
                    (1 - μ)*u[3]*(9*invr15 - 15*u[3]^2*invr17) +
                    μ*u[3]*(9*invr25 - 15*u[3]^2*invr27)];

    # Fill Jacobian
    jac[8,   1] = -dG11dr[1]*u[11] - dG12dr[1]*u[12] - dG13dr[1]*u[13]
    jac[8,   2] = -dG11dr[2]*u[11] - dG12dr[2]*u[12] - dG13dr[2]*u[13]
    jac[8,   3] = -dG11dr[3]*u[11] - dG12dr[3]*u[12] - dG13dr[3]*u[13]
    jac[9,   1] = -dG12dr[1]*u[11] - dG22dr[1]*u[12] - dG23dr[1]*u[13]
    jac[9,   2] = -dG12dr[2]*u[11] - dG22dr[2]*u[12] - dG23dr[2]*u[13]
    jac[9,   3] = -dG12dr[3]*u[11] - dG22dr[3]*u[12] - dG23dr[3]*u[13]
    jac[10,  1] = -dG13dr[1]*u[11] - dG23dr[1]*u[12] - dG33dr[1]*u[13]
    jac[10,  2] = -dG13dr[2]*u[11] - dG23dr[2]*u[12] - dG33dr[2]*u[13]
    jac[10,  3] = -dG13dr[3]*u[11] - dG23dr[3]*u[12] - dG33dr[3]*u[13]
    jac[1,   4] = 1.0
    jac[2,   5] = 1.0
    jac[3,   6] = 1.0
    jac[4,   1] = GVec[1]
    jac[4,   2] = GVec[4]
    jac[4,   3] = GVec[5]
    jac[5,   1] = GVec[4]
    jac[5,   2] = GVec[2]
    jac[5,   3] = GVec[6]
    jac[6,   1] = GVec[5]
    jac[6,   2] = GVec[6]
    jac[6,   3] = GVec[3]
    jac[4,   5] = 2.0
    jac[5,   4] = -2.0
    jac[8,  11] = -GVec[1]
    jac[8,  12] = -GVec[4]
    jac[8,  13] = -GVec[5]
    jac[9,  11] = -GVec[4]
    jac[9,  12] = -GVec[2]
    jac[9,  13] = -GVec[6]
    jac[10, 11] = -GVec[5]
    jac[10, 12] = -GVec[6]
    jac[10, 13] = -GVec[3]
    jac[11,  8] = -1.0
    jac[12,  9] = -1.0
    jac[13, 10] = -1.0
    jac[11, 12] = 2.0
    jac[12, 11] = -2.0

    # Fill remaining depending on utype
    if p.utype != 1
        compute_G1!(jac, u, λv, tMaxSc, γ)
        compute_G2!(jac, u, λv, tMaxSc, γ)
        compute_G7!(jac, u, λv, tMaxSc, γ)
        compute_G8!(jac, u, λv, tMaxSc, γ)
    else
        compute_G1!(jac, u, λv, tMaxSc, cSc, p.ϵ, γ)
        compute_G2!(jac, u, λv, tMaxSc, cSc, p.ϵ, γ)
        compute_G3!(jac, u, λv, tMaxSc, p.ϵ)
        compute_G4!(jac, u, λv, tMaxSc, p.ϵ)
        compute_G5!(jac, u, λv, tMaxSc, p.ϵ)
        compute_G6!(jac, p.ϵ, tMaxSc, cSc)
        compute_G7!(jac, u, λv, tMaxSc, cSc, p.ϵ, γ)
        compute_G8!(jac, u, λv, tMaxSc, cSc, p.ϵ, γ)
        compute_G9!(jac, u, λv, tMaxSc, p.ϵ)
    end
    return nothing
end

function compute_G1!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, cSc, ϵ, γ)
    @inbounds begin
        temp = γ*tMaxSc / (λv*u[7]*u[7]) + cSc*tMaxSc / (2.0*ϵ*u[7]^3)
        jac[4, 7] = temp*u[11]
        jac[5, 7] = temp*u[12]
        jac[6, 7] = temp*u[13]
    end
    return nothing
end

function compute_G1!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, γ)
    @inbounds begin
        temp    = γ*tMaxSc / (λv*u[7]*u[7])
        jac[4, 7] = temp*u[11]
        jac[5, 7] = temp*u[12]
        jac[6, 7] = temp*u[13]
    end
    return nothing
end

function compute_G2!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, γ)
    @inbounds begin
        λvInv    = 1.0 / λv
        λvInv3   = λvInv^3
        temp1    = -γ*tMaxSc/u[7]
        jac[4, 11] = temp1*(λvInv - u[11]*u[11]*λvInv3)
        jac[5, 12] = temp1*(λvInv - u[12]*u[12]*λvInv3)
        jac[6, 13] = temp1*(λvInv - u[13]*u[13]*λvInv3)
        jac[4, 12] = -temp1*u[11]*u[12]*λvInv3
        jac[4, 13] = -temp1*u[11]*u[13]*λvInv3
        jac[5, 11] = jac[4, 12]
        jac[5, 13] = -temp1*u[12]*u[13]*λvInv3
        jac[6, 11] = jac[4, 13]
        jac[6, 12] = jac[5, 13]
    end
    return nothing
end

function compute_G2!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, cSc, ϵ, γ)
    @inbounds begin
        λvInv   = 1.0 / λv
        λvInv2  = λvInv*λvInv
        λvInv3  = λvInv2*λvInv
        temp1   = -γ*tMaxSc / u[7]
        temp2   = -cSc*tMaxSc / (2.0*ϵ*u[7]^2)
        jac[4, 11] = temp1*(λvInv - u[11]*u[11]*λvInv3) + temp2*u[11]*u[11]*λvInv2
        jac[5, 12] = temp1*(λvInv - u[12]*u[12]*λvInv3) + temp2*u[12]*u[12]*λvInv2
        jac[6, 13] = temp1*(λvInv - u[13]*u[13]*λvInv3) + temp2*u[13]*u[13]*λvInv2
        jac[4, 12] = -temp1*u[11]*u[12]*λvInv3 + temp2*u[11]*u[12]*λvInv2
        jac[4, 13] = -temp1*u[11]*u[13]*λvInv3 + temp2*u[11]*u[13]*λvInv2
        jac[5, 11] = jac[4, 12]
        jac[5, 13] = -temp1*u[12]*u[13]*λvInv3 + temp2*u[12]*u[13]*λvInv2
        jac[6, 11] = jac[4, 13]
        jac[6, 12] = jac[5, 13]
    end
    return nothing
end

function compute_G3!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, ϵ)
    @inbounds begin
        temp        = -tMaxSc / (2.0*ϵ*λv*u[7])
        jac[4, 14]    = temp*u[11]
        jac[5, 14]    = temp*u[12]
        jac[6, 14]    = temp*u[13]
    end
    return nothing
end

function compute_G4!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, ϵ)
    @inbounds jac[7,7] = λv*tMaxSc / (2.0*ϵ*u[7]^2)
    return nothing
end

function compute_G5!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, ϵ)
    @inbounds begin
        temp = -tMaxSc / (2.0*ϵ*λv*u[7])
        jac[7, 11] = temp*u[11]
        jac[7, 12] = temp*u[12]
        jac[7, 13] = temp*u[13]
    end
    return nothing
end

function compute_G6!(jac::AbstractMatrix, ϵ, tMaxSc, cSc)
    @inbounds jac[7,14] = -tMaxSc / (2.0*ϵ*cSc)
    return nothing
end

function compute_G7!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, γ)
    @inbounds jac[14, 7] = 2.0*λv*γ*tMaxSc / (u[7]^3)
    return nothing
end

function compute_G7!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, cSc, ϵ, γ)
    @inbounds jac[14,7] = 2.0*λv*γ*tMaxSc / (u[7]^3) + λv^2*cSc*tMaxSc / (2.0*ϵ*u[7]^4)
    return nothing
end

function compute_G8!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, γ)
    @inbounds begin
        temp = -γ*tMaxSc / (λv*u[7]^2)
        jac[14, 11] = temp*u[11]
        jac[14, 12] = temp*u[12]
        jac[14, 13] = temp*u[13]
    end
    return nothing
end

function compute_G8!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, cSc, ϵ, γ)
    @inbounds begin
        temp = -γ*tMaxSc / (λv*u[7]^2) - cSc*tMaxSc / (2.0*ϵ*u[7]^3)
        jac[14, 11] = temp*u[11]
        jac[14, 12] = temp*u[12]
        jac[14, 13] = temp*u[13]
    end
    return nothing
end

function compute_G9!(jac::AbstractMatrix, u::AbstractVector, λv, tMaxSc, ϵ)
    @inbounds jac[14,14] = -λv*tMaxSc / (2.0*ϵ*u[7]^2)
    return nothing
end


text:
function [F]=FcnDiff(U,CG,Param)
nz=Param.Grid.nz;
RhoPos=Param.RhoPos;
uPos=Param.uPos;
vPos=Param.vPos;
wPos=Param.wPos;
ThPos=Param.ThPos;
W=[zeros(size(U,1),1)...
   ,U(:,:,wPos)];

if Param.HyperD>0
  if strcmp(Param.Thermo,'Energy')
    [FU,FW,FTh]=DiffusionVec(U(:,:,uPos:vPos),W...
      ,U(:,:,ThPos)+Param.Pres...
      ,U(:,:,RhoPos),CG,Param);
  else
    [FU,FW,FTh]=DiffusionVec(U(:,:,uPos:vPos),W,U(:,:,ThPos)...
      ,U(:,:,RhoPos),CG,Param);
  end
end

F=zeros(size(U));
F(:,:,uPos:vPos)=FU;
F(:,:,wPos)=FW(:,2:nz+1);
F(:,:,ThPos)=FTh;
end


text:
@testset "Roundtrip" begin

  i, f, s, b = 10, 12.3, "a", false

  cl = Jib.ComboLeg(ratio=1, exchange="ex")
  cv = Jib.ConditionVolume("o", true, 1, 2, "ex")

  o = Jib.Requests.Encoder(IOBuffer())

  o(i, f, s, b, Jib.Requests.splat(cl)..., cv)

  m = split(String(take!(o.buf)), '\0')

  it = Iterators.Stateful(Jib.Reader.Decoder.Field.(m))

  j::Int,
  g::Float64,
  z::String,
  l::Bool = it

  cc = Jib.Reader.Decoder.slurp(Jib.ComboLeg, it)

  cw = Jib.Reader.Decoder.slurp(Jib.condition_map[Jib.Reader.Decoder.slurp(Jib.ConditionType, it)], it)

  @test i == j
  @test f == g
  @test s == z
  @test b == l
  @test cl == cc
  @test cv == cw

  popfirst!(it)
  @test isempty(it)

end


text:
function TTest()
    return TTest(())
end

function homoscedastic_t(obj::TTest, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(obj, "homoscedasticT", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function homoscedastic_t(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "homoscedasticT", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function homoscedastic_t_test(obj::TTest, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(obj, "homoscedasticTTest", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function homoscedastic_t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "homoscedasticTTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function homoscedastic_t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "homoscedasticTTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function paired_t(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "pairedT", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function paired_t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "pairedTTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function paired_t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "pairedTTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function t(obj::TTest, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(obj, "t", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function t(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "t", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function t(obj::TTest, arg0::jdouble, arg1::StatisticalSummary)
    return jcall(obj, "t", jdouble, (jdouble, StatisticalSummary), arg0, arg1)
end

function t(obj::TTest, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "t", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end

function t_test(obj::TTest, arg0::StatisticalSummary, arg1::StatisticalSummary)
    return jcall(obj, "tTest", jdouble, (StatisticalSummary, StatisticalSummary), arg0, arg1)
end

function t_test(obj::TTest, arg0::StatisticalSummary, arg1::StatisticalSummary, arg2::jdouble)
    return jcall(obj, "tTest", jboolean, (StatisticalSummary, StatisticalSummary, jdouble), arg0, arg1, arg2)
end

function t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "tTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function t_test(obj::TTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "tTest", jboolean, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function t_test(obj::TTest, arg0::jdouble, arg1::StatisticalSummary)
    return jcall(obj, "tTest", jdouble, (jdouble, StatisticalSummary), arg0, arg1)
end

function t_test(obj::TTest, arg0::jdouble, arg1::StatisticalSummary, arg2::jdouble)
    return jcall(obj, "tTest", jboolean, (jdouble, StatisticalSummary, jdouble), arg0, arg1, arg2)
end

function t_test(obj::TTest, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "tTest", jdouble, (jdouble, Vector{jdouble}), arg0, arg1)
end

function t_test(obj::TTest, arg0::jdouble, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "tTest", jboolean, (jdouble, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end


text:
#InspectDR: Functionality for user-overwritable styles (ex: used in stylesheets)
#-------------------------------------------------------------------------------

#=
Also contains mechanism to overwrite certain properties of a structure.
=#

#==Main types
===============================================================================#
abstract type AbstractStyle end #Type that can be made into a "StyleType"

#Mechanism to overwrite only certain properties from a function call:
struct KeepPrevType end #Keep previous value
const KEEP_PREV = KeepPrevType()

struct UseDefaultType end #Use style default
const USE_DFLT = UseDefaultType()

#Defines a "style" object (to implement user-overwritable styles):
mutable struct StyleType{T<:AbstractStyle}
	values::T #Active values
	defaults::T #Defaults for this style
	overwrite::Set{Symbol} #Masks out which elements user wishes to overwrite
end
StyleType(defaults::T) where T<:AbstractStyle =
	StyleType(deepcopy(defaults), defaults, Set{Symbol}([]))
StyleType(values::T, defaults::T) where T<:AbstractStyle =
	StyleType(values, defaults, Set{Symbol}([]))

#==Main functions
===============================================================================#
overwriteprop!(s, key::Symbol, value::UseDefaultType) = nothing #Don't change value
overwriteprop!(s, key::Symbol, value::KeepPrevType) = nothing #Don't change value
function overwriteprop!(s, key::Symbol, value)
	value = convert(typeof(getfield(s, key)), value) #Need proper type
	setfield!(s, key, value)
end

#Want to use default, so ensure key not in overwrite "set"
function overwritestyleattrib(overwrite::Set{Symbol}, key::Symbol, value::UseDefaultType)
	if key in overwrite
		pop!(overwrite, key)
	end
	return
end
#Make sure marked as overwritten:
overwritestyleattrib(overwrite::Set{Symbol}, key::Symbol, value) = push!(overwrite, key)
#Mark given keys as overwritten:
markoverwritten!(s::StyleType, keylist::Set{Symbol}) = union!(s.overwrite, keylist)

function overwrite!(s::StyleType, key::Symbol, value)
	overwriteprop!(s.values, key, value)
	overwritestyleattrib(s.overwrite, key, value)
end

#Synchronize values with defaults if they are out of sync.
function refresh!(s::StyleType{T}) where T
	proplist = fieldnames(T)
	for key in proplist
		if !(key in s.overwrite)
			value = deepcopy(getfield(s.defaults, key))
			setfield!(s.values, key, value)
		end
	end
	return
end

#Changes the default to the given value, and refreshes .value fields, if requested
#TODO: rename setdefaults!()??
function setstyle!(s::StyleType{T}, value::T; refresh::Bool=true) where T<:AbstractStyle
	s.defaults = value
	if refresh
		refresh!(s)
	end
	return
end

#==Mutators
===============================================================================#
Base.setindex!(s::StyleType, value, key::Symbol) = overwrite!(s, key, value)
Base.getindex(s::StyleType, key::Symbol) = getfield(s.values, key)

#Last line