text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

function MOI.features_available(d::NLPEvaluator)
    # Check if we are missing any hessians for user-defined multivariate
    # operators, in which case we need to disable :Hess and :HessVec.
    d.disable_2ndorder = any(
        op -> op.∇²f === nothing,
        d.data.operators.registered_multivariate_operators,
    )
    if d.disable_2ndorder
        return [:Grad, :Jac, :JacVec]
    end
    return [:Grad, :Jac, :JacVec, :Hess, :HessVec]
end

function MOI.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
    # Check that we support the features requested by the user.
    available_features = MOI.features_available(d)
    for feature in requested_features
        if !(feature in available_features)
            error("Unsupported feature $feature")
        end
    end
    moi_index_to_consecutive_index =
        Dict(x => i for (i, x) in enumerate(d.ordered_variables))
    N = length(moi_index_to_consecutive_index)
    #
    largest_user_input_dimension = 1
    for op in d.data.operators.registered_multivariate_operators
        largest_user_input_dimension = max(largest_user_input_dimension, op.N)
    end
    d.objective = nothing
    d.user_output_buffer = zeros(largest_user_input_dimension)
    d.jac_storage = zeros(max(N, largest_user_input_dimension))
    d.constraints = _FunctionStorage[]
    d.last_x = fill(NaN, N)
    d.want_hess = :Hess in requested_features
    want_hess_storage = (:HessVec in requested_features) || d.want_hess
    coloring_storage = Coloring.IndexedSet(N)
    max_expr_length = 0
    #
    main_expressions = [c.expression.nodes for (_, c) in d.data.constraints]
    if d.data.objective !== nothing
        pushfirst!(main_expressions, d.data.objective.nodes)
    end
    d.subexpression_order, individual_order = _order_subexpressions(
        main_expressions,
        [expr.nodes for expr in d.data.expressions],
    )
    num_subexpressions = length(d.data.expressions)
    d.subexpression_linearity = Vector{Linearity}(undef, num_subexpressions)
    subexpression_variables = Vector{Vector{Int}}(undef, num_subexpressions)
    subexpression_edgelist =
        Vector{Set{Tuple{Int,Int}}}(undef, num_subexpressions)
    d.subexpressions = Vector{_SubexpressionStorage}(undef, num_subexpressions)
    d.subexpression_forward_values = zeros(num_subexpressions)
    d.subexpression_reverse_values = zeros(num_subexpressions)
    for k in d.subexpression_order
        # Only load expressions which actually are used
        d.subexpression_forward_values[k] = NaN
        subex = _SubexpressionStorage(
            d.data.expressions[k],
            d.subexpression_linearity,
            moi_index_to_consecutive_index,
        )
        d.subexpressions[k] = subex
        d.subexpression_linearity[k] = subex.linearity
        if d.want_hess
            empty!(coloring_storage)
            _compute_gradient_sparsity!(coloring_storage, subex.nodes)
            # union with all dependent expressions
            for idx in _list_subexpressions(subex.nodes)
                union!(coloring_storage, subexpression_variables[idx])
            end
            subexpression_variables[k] = collect(coloring_storage)
            empty!(coloring_storage)
            linearity = _classify_linearity(
                subex.nodes,
                subex.adj,
                d.subexpression_linearity,
            )
            edgelist = _compute_hessian_sparsity(
                subex.nodes,
                subex.adj,
                linearity,
                coloring_storage,
                subexpression_edgelist,
                subexpression_variables,
            )
            subexpression_edgelist[k] = edgelist
        end
    end
    max_chunk = 1
    if d.data.objective !== nothing
        d.objective = _FunctionStorage(
            main_expressions[1],
            d.data.objective.values,
            N,
            coloring_storage,
            d.want_hess,
            d.subexpressions,
            individual_order[1],
            d.subexpression_linearity,
            subexpression_edgelist,
            subexpression_variables,
            moi_index_to_consecutive_index,
        )
        max_expr_length = max(max_expr_length, length(d.objective.nodes))
        max_chunk = max(max_chunk, size(d.objective.seed_matrix, 2))
    end
    for (k, (_, constraint)) in enumerate(d.data.constraints)
        idx = d.data.objective !== nothing ? k + 1 : k
        push!(
            d.constraints,
            _FunctionStorage(
                main_expressions[idx],
                constraint.expression.values,
                N,
                coloring_storage,
                d.want_hess,
                d.subexpressions,
                individual_order[idx],
                d.subexpression_linearity,
                subexpression_edgelist,
                subexpression_variables,
                moi_index_to_consecutive_index,
            ),
        )
        max_expr_length = max(max_expr_length, length(d.constraints[end].nodes))
        max_chunk = max(max_chunk, size(d.constraints[end].seed_matrix, 2))
    end
    # 10 is hardcoded upper bound to avoid excess memory allocation
    max_chunk = min(max_chunk, 10)
    if d.want_hess || want_hess_storage
        d.input_ϵ = zeros(max_chunk * N)
        d.output_ϵ = zeros(max_chunk * N)
        #
        len = max_chunk * max_expr_length
        d.forward_storage_ϵ = zeros(len)
        d.partials_storage_ϵ = zeros(len)
        d.reverse_storage_ϵ = zeros(len)
        #
        len = max_chunk * length(d.subexpressions)
        d.subexpression_forward_values_ϵ = zeros(len)
        d.subexpression_reverse_values_ϵ = zeros(len)
        #
        for k in d.subexpression_order
            len = max_chunk * length(d.subexpressions[k].nodes)
            d.subexpressions[k].forward_storage_ϵ = zeros(len)
            d.subexpressions[k].partials_storage_ϵ = zeros(len)
            d.subexpressions[k].reverse_storage_ϵ = zeros(len)
        end
        d.max_chunk = max_chunk
        if d.want_hess
            _compute_hessian_lagrangian_structure(d)
        end
    end
    return
end

function MOI.eval_objective(d::NLPEvaluator, x)
    if d.objective === nothing
        error("No nonlinear objective.")
    end
    _reverse_mode(d, x)
    return d.objective.forward_storage[1]
end

function MOI.eval_objective_gradient(d::NLPEvaluator, g, x)
    if d.objective === nothing
        error("No nonlinear objective.")
    end
    _reverse_mode(d, x)
    fill!(g, 0.0)
    _extract_reverse_pass(g, d, d.objective)
    return
end

function MOI.eval_constraint(d::NLPEvaluator, g, x)
    _reverse_mode(d, x)
    for i in 1:length(d.constraints)
        g[i] = d.constraints[i].forward_storage[1]
    end
    return
end

function MOI.jacobian_structure(d::NLPEvaluator)
    J = Tuple{Int64,Int64}[]
    for (row, constraint) in enumerate(d.constraints)
        for col in constraint.grad_sparsity
            push!(J, (row, col))
        end
    end
    return J
end

function MOI.eval_constraint_jacobian(d::NLPEvaluator, J, x)
    _reverse_mode(d, x)
    fill!(J, 0.0)
    offset = 0
    for ex in d.constraints
        for i in ex.grad_sparsity
            d.jac_storage[i] = 0.0
        end
        _extract_reverse_pass(d.jac_storage, d, ex)
        for (k, idx) in enumerate(ex.grad_sparsity)
            J[offset+k] = d.jac_storage[idx]
        end
        offset += length(ex.grad_sparsity)
    end
    return
end

function MOI.eval_constraint_jacobian_product(d::NLPEvaluator, y, x, w)
    fill!(y, 0.0)
    J_struct = MOI.jacobian_structure(d)
    nnz = length(J_struct)
    J = zeros(nnz)
    MOI.eval_constraint_jacobian(d, J, x)
    for (k, (i, j)) in enumerate(J_struct)
        y[i] += J[k] * w[j]
    end
    return
end

function MOI.eval_constraint_jacobian_transpose_product(
    d::NLPEvaluator,
    y::AbstractVector{Float64},
    x::AbstractVector{Float64},
    w::AbstractVector{Float64},
)
    fill!(y, 0.0)
    J_struct = MOI.jacobian_structure(d)
    nnz = length(J_struct)
    J = zeros(nnz)
    MOI.eval_constraint_jacobian(d, J, x)
    for (k, (i, j)) in enumerate(J_struct)
        y[j] += J[k] * w[i]
    end
    return y
end

function MOI.hessian_lagrangian_structure(d::NLPEvaluator)
    @assert d.want_hess
    return d.hessian_sparsity
end

function _compute_hessian_lagrangian_structure(d::NLPEvaluator)
    d.hessian_sparsity = Tuple{Int64,Int64}[]
    if d.objective !== nothing
        append!(d.hessian_sparsity, zip(d.objective.hess_I, d.objective.hess_J))
    end
    for c in d.constraints
        append!(d.hessian_sparsity, zip(c.hess_I, c.hess_J))
    end
    return
end

function MOI.eval_hessian_lagrangian(d::NLPEvaluator, H, x, σ, μ)
    @assert d.want_hess
    _reverse_mode(d, x)
    fill!(d.input_ϵ, 0.0)
    offset = 0
    if d.objective !== nothing
        offset += _eval_hessian(d, d.objective, H, σ, offset)::Int
    end
    for (i, ex) in enumerate(d.constraints)
        offset += _eval_hessian(d, ex, H, μ[i], offset)::Int
    end
    return
end

function MOI.eval_hessian_lagrangian_product(d::NLPEvaluator, h, x, v, σ, μ)
    _reverse_mode(d, x)
    fill!(h, 0.0)
    T = ForwardDiff.Partials{1,Float64}
    input_ϵ = reinterpret(T, d.input_ϵ)
    output_ϵ = reinterpret(T, d.output_ϵ)
    for i in 1:length(x)
        input_ϵ[i] = ForwardDiff.Partials((v[i],))
    end
    # forward evaluate all subexpressions once
    subexpr_forward_values_ϵ = reinterpret(T, d.subexpression_forward_values_ϵ)
    for i in d.subexpression_order
        subexpr = d.subexpressions[i]
        subexpr_forward_values_ϵ[i] = _forward_eval_ϵ(
            d,
            subexpr,
            reinterpret(T, subexpr.forward_storage_ϵ),
            reinterpret(T, subexpr.partials_storage_ϵ),
            input_ϵ,
            subexpr_forward_values_ϵ,
            d.data.operators,
        )
    end
    # we only need to do one reverse pass through the subexpressions as well
    subexpr_reverse_values_ϵ = reinterpret(T, d.subexpression_reverse_values_ϵ)
    fill!(subexpr_reverse_values_ϵ, zero(T))
    fill!(d.subexpression_reverse_values, 0.0)
    fill!(d.reverse_storage_ϵ, 0.0)
    fill!(output_ϵ, zero(T))
    if d.objective !== nothing
        _forward_eval_ϵ(
            d,
            d.objective,
            reinterpret(T, d.forward_storage_ϵ),
            reinterpret(T, d.partials_storage_ϵ),
            input_ϵ,
            subexpr_forward_values_ϵ,
            d.data.operators,
        )
        _reverse_eval_ϵ(
            output_ϵ,
            d.objective,
            reinterpret(T, d.reverse_storage_ϵ),
            reinterpret(T, d.partials_storage_ϵ),
            d.subexpression_reverse_values,
            subexpr_reverse_values_ϵ,
            σ,
            zero(T),
        )
    end
    for (i, con) in enumerate(d.constraints)
        _forward_eval_ϵ(
            d,
            con,
            reinterpret(T, d.forward_storage_ϵ),
            reinterpret(T, d.partials_storage_ϵ),
            input_ϵ,
            subexpr_forward_values_ϵ,
            d.data.operators,
        )
        _reverse_eval_ϵ(
            output_ϵ,
            con,
            reinterpret(T, d.reverse_storage_ϵ),
            reinterpret(T, d.partials_storage_ϵ),
            d.subexpression_reverse_values,
            subexpr_reverse_values_ϵ,
            μ[i],
            zero(T),
        )
    end
    for i in length(d.subexpression_order):-1:1
        j = d.subexpression_order[i]
        subexpr = d.subexpressions[j]
        _reverse_eval_ϵ(
            output_ϵ,
            subexpr,
            reinterpret(T, subexpr.reverse_storage_ϵ),
            reinterpret(T, subexpr.partials_storage_ϵ),
            d.subexpression_reverse_values,
            subexpr_reverse_values_ϵ,
            d.subexpression_reverse_values[j],
            subexpr_reverse_values_ϵ[j],
        )
    end
    for i in 1:length(x)
        h[i] += output_ϵ[i].values[1]
    end
    return
end


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff,
    BenchmarkTools, ProfileView

## Observed Data #############################################################

data_small = Feather.read("test/comparisons/data_unique_small.feather")
data_big = Feather.read("test/comparisons/data_unique_big.feather")
data_huge = Feather.read("test/comparisons/data_unique_huge.feather")

data_def_small = Feather.read("test/comparisons/data_def_unique_small.feather")
data_def_big = Feather.read("test/comparisons/data_def_unique_big.feather")
data_def_huge = Feather.read("test/comparisons/data_def_unique_huge.feather")

pars_small = Feather.read("test/comparisons/pars_unique_small.feather")
pars_big = Feather.read("test/comparisons/pars_unique_big.feather")
pars_huge = Feather.read("test/comparisons/pars_unique_huge.feather")

## small ############################################################

semobserved_small =
    SemObsCommon(data = Matrix(data_small); meanstructure = true)

diff_fin_small = SemFiniteDiff(BFGS(), Optim.Options())

@ModelingToolkit.variables x[1:7], m[1:2], load_t[1:4]

S =[x[1]  0     0     0     0     0
    0     x[2]  0     0     0     0
    0     0     x[3]  0     0     0
    0     0     0     x[4]  0     0
    0     0     0     0     x[5]  x[7]
    0     0     0     0     x[7]  x[6]]

F =[1.0 0 0 0 0 0
    0 1 0 0 0 0
    0 0 1 0 0 0
    0 0 0 1 0 0]

A =[0  0  0  0  1.0  load_t[1]
    0  0  0  0  1  load_t[2]
    0  0  0  0  1  load_t[3]
    0  0  0  0  1  load_t[4]
    0  0  0  0  0  0
    0  0  0  0  0  0]

M = [0.0, 0, 0, 0, m[1:2]...]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

data_def_small = Matrix(data_def_small)

start_val = [1.0, 1, 1, 1, 1, 1, 1, 1.0, 1.0]

imply_small = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def_small
)

loss_small = Loss([SemDefinition(semobserved_small, imply_small, 0.0, 0.0)])

model_fin_small = Sem(semobserved_small, imply_small, loss_small, diff_fin_small)

solution_small = sem_fit(model_fin_small)

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution_small.minimizer .- pars_small.Estimate[par_order]
        ) .< 0.001*abs.(pars_small.Estimate[par_order]))

@benchmark sem_fit(model_fin_small)

## big ############################################################

semobserved_big =
    SemObsCommon(data = Matrix(data_big); meanstructure = true)

diff_fin_big = SemFiniteDiff(LBFGS(), Optim.Options(x_tol = 1e-9))

@ModelingToolkit.variables x[1:18], m[1:2], load_t[1:15]

F = zeros(15, 17)

diag(F)
F[diagind(F)] .= 1.0

M = zeros(Expression, 17, 1)
M .= [zeros(15)..., m[1:2]...]
M = sparse(M)

#S
Ind = [1:17..., 16, 17]; J = [1:17..., 17, 16]; V = [x[1:17]..., x[18], x[18]]
S = sparse(Ind, J, V)

#F
F = sparse(F)

#A
Ind = [1:15..., 1:15...]; J = [fill(16, 15)..., fill(17, 15)...]; V = [ones(15)..., load_t...]
A = sparse(Ind, J, V, 17, 17)

data_def_big = Matrix(data_def_big)

start_val = ones(20) #[ones(15)..., 0.05, 0.05, 0.0, 1.0, 1.0]

imply_big = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:18]..., m[1:2]...],
    load_t,
    start_val,
    data_def_big
)

loss_big = Loss([SemDefinition(semobserved_big, imply_big, 0.0, 0.0)])

model_fin_big = Sem(semobserved_big, imply_big, loss_big, diff_fin_big)

solution_big = sem_fit(model_fin_big)

par_order = [collect(1:16); 18; 17; 19; 20]

all(
    abs.(solution_big.minimizer .- pars_big.Estimate[par_order]
        ) .< 0.001*abs.(pars_big.Estimate[par_order]))

@benchmark sem_fit(model_fin_big)

## huge ############################################################

semobserved_huge =
    SemObsCommon(data = Matrix(data_huge); meanstructure = true)

diff_fin_huge = SemFiniteDiff(BFGS(), Optim.Options())

@ModelingToolkit.variables x[1:33], m[1:2], load_t[1:30]

F = zeros(30, 32)

diag(F)
F[diagind(F)] .= 1.0

M = [zeros(30)..., m[1:2]...]

#S
Ind = [1:32..., 31, 32]; J = [1:32..., 32, 31];
V = [x[1:32]..., x[33], x[33]];

S = sparse(Ind, J, V)

#F
F = sparse(F)

#A
Ind = [1:30..., 1:30...]; J = [fill(31, 30)..., fill(32, 30)...];
V = [ones(30)..., load_t...];
A = sparse(Ind, J, V, 32, 32)

data_def_huge = Matrix(data_def_huge)

start_val = ones(35)#[ones(30)..., 0.05, 0.05, 0.0, 1.0, 1.0]

imply_huge = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:33]..., m[1:2]...],
    load_t,
    start_val,
    data_def_huge
)

loss_huge = Loss([SemDefinition(semobserved_huge, imply_huge, 0.0, 0.0)])

model_fin_huge = Sem(semobserved_huge, imply_huge, loss_huge, diff_fin_huge)

solution_huge = sem_fit(model_fin_huge)

par_order = [collect(1:31); 33; 32; 34; 35]

all(
    abs.(solution_huge.minimizer .- pars_huge.Estimate[par_order]
        ) .< 0.001*abs.(pars_huge.Estimate[par_order]))

@benchmark sem_fit(model_fin_huge)

## benchmark #######################################################
start_val = [ones(15)..., 0.05, 0.05, 0.0, 1.0, 1.0]

start_val = [ones(30)..., 0.05, 0.05, 0.0, 1.0, 1.0]

@code_warntype model_fin_big.imply(start_val)

@code_warntype model_fin_big(start_val)

@code_warntype model_fin_big.loss.functions[1](start_val, model_fin_big)

@benchmark model_fin_big.imply($start_val)

@benchmark model_fin_big(start_val)
ProfileView.@profview model_fin_big(start_val)

@benchmark model_fin_big.loss.functions[1](start_val, model_fin_big)

@benchmark model_fin_huge.imply($start_val)

@benchmark model_fin_huge($start_val)

@benchmark model_fin_big.loss.functions[1]($start_val, $model_fin_big)
#1.5ms

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end
# run once to trigger compilation (ignore this one)

@btime profile_test(10)
@btime to_prof(100, start_val, model_fin_big)

function to_prof(n, par, model)
    for i = 1:n
        model.loss.functions[1](par, model)
    end
end

ProfileView.@profview to_prof(1000, start_val, model_fin_big)

using LinearAlgebra, BenchmarkTools
BLAS.vendor()

using CUDA, MKL

a = rand(30,30); a = a'*a

BLAS.set_num_threads(8)

set_max_threads(n) = ccall((:mkl_set_num_threads, libmkl_rt), Cvoid, (Ptr{Int32},), Ref(Int32(n)));
set_max_threads(3)
MKL.typemin()

@benchmark inv(a)

a = CUDA.rand(Float64, 30, 30)
a32 = CUDA.rand(30, 30)

a32 isa StridedCuArray

inv(a32)

b = rand(30,30)
b32 = rand(Float32, 30 , 30)

@benchmark inv(b)

@benchmark inv(b32)

a = rand(5)

### bounds

a = rand(30,30); a = a*a'
matvec = [a for i = 1:500]
safe = copy(matvec)
pre = similar(matvec)

copyto!(safe, matvec)

cholvec = cholesky.(matvec)

function myf(pre, matvec)
    @inbounds @fastmath for i = 1:size(pre, 1)
        pre[i] = matvec[i]^2
    end
end

function myf2(pre, matvec)
    @inbounds @fastmath for i = 1:size(pre, 1)
        pre[i] = matvec[i]^2
    end
end

@benchmark myf(pre, cholvec)

@benchmark myf2(pre, matvec)

# 0.2720108 secs vs 100ms
# 36.52637 secs vs 30s (21s)
# 156.6648 secs vs 157s

ext: jl


text:
using PyPlot;

function plotPos(tSim,rSol, titleStr::String="Pos")
    figure(); clf();
    pygui(true);
    subplot(3,1,1);
    PyPlot.plot(tSim,rSol[:,1]); xlabel("t"); ylabel("x");
    subplot(3,1,2);
    PyPlot.plot(tSim,rSol[:,2]); xlabel("t"); ylabel("y")
    subplot(3,1,3);
    PyPlot.plot(tSim,rSol[:,3]); xlabel("t"); ylabel("z")
    PyPlot.suptitle(titleStr)
    tight_layout();
end

function plotVel(tSim,vSol, titleStr::String="Vel")
    figure(); clf();
    pygui(true);
    subplot(3,1,1);
    PyPlot.plot(tSim,vSol[:,1]); xlabel("t"); ylabel("u");
    subplot(3,1,2);
    PyPlot.plot(tSim,vSol[:,2]); xlabel("t"); ylabel("v")
    subplot(3,1,3);
    PyPlot.plot(tSim,vSol[:,3]); xlabel("t"); ylabel("w")
    PyPlot.suptitle(titleStr)
    tight_layout();

end

function plotQuat(tSim,βSol, titleStr::String="Quaternion")
    figure(); clf();
    pygui(true);
    subplot(4,1,1);
    PyPlot.plot(tSim,βSol[:,1]); xlabel("t"); ylabel("β0");
    subplot(4,1,2);
    PyPlot.plot(tSim,βSol[:,2]); xlabel("t"); ylabel("β1")
    subplot(4,1,3);
    PyPlot.plot(tSim,βSol[:,3]); xlabel("t"); ylabel("β2")
    subplot(4,1,4);
    PyPlot.plot(tSim,βSol[:,4]); xlabel("t"); ylabel("β3")
    PyPlot.suptitle(titleStr)
    tight_layout();

end

function plotQuatDot(tSim,βdotSol, titleStr::String="Quat Dot")
    figure(); clf();
    pygui(true);
    subplot(4,1,1);
    PyPlot.plot(tSim,βdotSol[:,1]); xlabel("t"); ylabel("̇β̇0");
    subplot(4,1,2);
    PyPlot.plot(tSim,βdotSol[:,2]); xlabel("t"); ylabel("β̇1")
    subplot(4,1,3);
    PyPlot.plot(tSim,βdotSol[:,3]); xlabel("t"); ylabel("β̇2")
    subplot(4,1,4);
    PyPlot.plot(tSim,βdotSol[:,4]); xlabel("t"); ylabel("β̇3")
    PyPlot.suptitle(titleStr)
    tight_layout();

end

function plotAngVel(tSim,ωSol, titleStr::String="Ang Vel")
    figure(); clf;
    pygui(true);
    subplot(3,1,1);
    PyPlot.plot(tSim,ωSol[:,1]); xlabel("t"); ylabel("ω1");
    subplot(3,1,2);
    PyPlot.plot(tSim,ωSol[:,2]); xlabel("t"); ylabel("ω2")
    subplot(3,1,3);
    PyPlot.plot(tSim,ωSol[:,3]); xlabel("t"); ylabel("ω3")
    PyPlot.suptitle(titleStr)
    tight_layout();

end

function plotErrNorm(tSim,βSol)
    errNormβ = zeros(length(tSim),1)
    for i=1:length(tSim)
        errNormβ[i] = norm(βSol[i,:],2)^2 - 1.0
    end
    errNorm = broadcast(abs,errNormβ)
    figure(); clf;
    pygui(true)
    plot(tSim,errNorm); title("Drift in ||q||")
    tight_layout();

end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
#make sure the clamped spline is working
x = 2*π*sort(rand(100)) .- 3
y = sin.(x)
ϕ = CubicSplineInterpolator(x, y, cos(x[1]), cos(x[end]))
xx = LinRange(minimum(x), maximum(x), 1000)
err = abs.(ϕ.(xx) - sin.(xx))
@test maximum(err) < 0.001

#test cubic hermite
f(x) = 3x^3 - 2x^2 + 5x - 1
df(x) = 9x^2 - 4x + 5
xa = -2
xb = 3.1
xx = LinRange(xa, xb, 100)
yy = cubichermite.(xx, xa, xb, f(xa), f(xb), df(xa), df(xb))
@test maximum(@. abs(yy - f.(xx)) ) < 1e-6

#test the type-flexible Bichebyshev functor
P = BichebyshevInterpolator((x,y)->sin(x)*cos(x), -π, π, 14, -π, π, 16)
for i in 1:10
    a, b = 2π*rand(2) .- π
    @test P(a,b) ≈ P(BigFloat(a),b)
end

#test cheby derivative
p = ChebyshevInterpolator(x->sin(x^2), -2π, 2π, 128)
d = chebyderiv(p)
for i = 1:10
    a = 4π*rand() - 2π
    @test d(a) - 2*a*cos(a^2) < 1e-6
end


text:
# ---
# title: 986. Interval List Intersections
# id: problem986
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers
# link: <https://leetcode.com/problems/interval-list-intersections/description/>
# hidden: true
# ---
#
# Given two lists of **closed** intervals, each list of intervals is pairwise
# disjoint and in sorted order.
#
# Return the intersection of these two interval lists.
#
# _(Formally, a closed interval`[a, b]` (with `a <= b`) denotes the set of real
# numbers `x` with `a <= x <= b`.  The intersection of two closed intervals is a
# set of real numbers that is either empty, or can be represented as a closed
# interval.  For example, the intersection of [1, 3] and [2, 4] is [2, 3].)_
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/30/interval1.png)**
#
#
#
#     Input: A = [[0,2],[5,10],[13,23],[24,25]], B = [[1,5],[8,12],[15,24],[25,26]]
#     Output: [[1,2],[5,5],[8,10],[15,23],[24,24],[25,25]]
#
#
#
#
# **Note:**
#
#   1. `0 <= A.length < 1000`
#   2. `0 <= B.length < 1000`
#   3. `0 <= A[i].start, A[i].end, B[i].start, B[i].end < 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module RocketPairwiseOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: pairwise()" begin

    println("Testing: operator pairwise()")

    run_proxyshowcheck("Pairwise", pairwise())

    run_testset([
        (
            source      = from(1:5) |> pairwise(),
            values      = @ts([ (1, 2), (2, 3), (3, 4), (4, 5), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = of(1) |> pairwise(),
            values      = @ts(c),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = completed() |> pairwise(),
            values      = @ts(c)
        ),
        (
            source      = faulted(Int, "e") |> pairwise(),
            values      = @ts(e("e"))
        ),
        (
            source      = never() |> pairwise(),
            values      = @ts()
        )
    ])

end

end


text:
struct KPDSupply{T<:Integer}
    item::String
    weight::T
    value::T
    quant::T
end

KPDSupply{T<:Integer}(itm::AbstractString, w::T, v::T, q::T=one(T)) = KPDSupply(itm, w, v, q)
Base.show(io::IO, kdps::KPDSupply) = print(io, kdps.quant, " ", kdps.item, " ($(kdps.weight) kg, $(kdps.value) €)")

using MathProgBase, Cbc
function solve(gear::Vector{<:KPDSupply}, capacity::Integer)
    w = getfield.(gear, :weight)
    v = getfield.(gear, :value)
    sol = mixintprog(-v, w', '<', capacity, :Bin, 0, 1, CbcSolver())
    gear[sol.sol .≈ 1]
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using RecipesBase

@recipe function f(img::AstroImage, header_number::Int)
    seriestype   := :heatmap
    aspect_ratio := :equal
    # Right now we only support single frame images,
    # gray scale is a good choice.
    color        := :grays
    img.data[header_number]
end

@recipe function f(img::AstroImage)
    seriestype   := :heatmap
    aspect_ratio := :equal
    # Right now we only support single frame images,
    # gray scale is a good choice.
    color        := :grays
    img.data[1]
end

@recipe function f(img::AstroImage, wcs::WCSTransform)
    seriestype   := :heatmap
    aspect_ratio := :equal
    color        := :grays
    xformatter   := x -> pix2world_xformatter(x, wcs)
    yformatter   := y -> pix2world_yformatter(y, wcs)
    xlabel       := labler_x(wcs)
    ylabel       := labler_y(wcs)
    img.data
end

function pix2world_xformatter(x, wcs::WCSTransform)
    res = round(pix_to_world(wcs, [float(x), float(x)])[1], digits = 2)
    if wcs.cunit[1] == "deg"       # TODO: add symbols for more units
        return string(res)*"°"
    else
        return res[1]
    end
end

function pix2world_yformatter(x, wcs::WCSTransform)
    res = round(pix_to_world(wcs, [float(x), float(x)])[2], digits = 2)
    if wcs.cunit[2] == "deg"       # TODO: add symbols for more units
        return string(res)*"°"
    else
        return res[1]
    end
end

function labler_x(wcs::WCSTransform)
    if length(wcs.ctype[1]) == 0
        return wcs.radesys
    elseif wcs.ctype[1][1:2] == "RA"
        return "Right Ascension (ICRS)"
    elseif wcs.ctype[1][1:4] == "GLON"
        return "Galactic Coordinate"
    elseif wcs.ctype[1][1:4] == "TLON"
        return "ITRS"
    else
        return wcs.ctype[1]
    end
end

function labler_y(wcs::WCSTransform)
    if length(wcs.ctype[2]) == 0
        return wcs.radesys
    elseif wcs.ctype[2][1:3] == "DEC"
        return "Declination (ICRS)"
    elseif wcs.ctype[2][1:4] == "GLAT"
        return "Galactic Coordinate"
    elseif wcs.ctype[2][1:4] == "TLAT"
        return "ITRS"
    else
        return wcs.ctype[2]
    end
end


text:
export printRow, drawLine!, setPoint!, setPixel!, addRow!

@deprecate printRow(io, canvas, row) printrow(io, canvas, row)

@deprecate drawLine!(args...; nargs...) lines!(args...; nargs...)

@deprecate setPoint!(args...; nargs...) points!(args...; nargs...)

@deprecate setPixel!(args...; nargs...) pixel!(args...; nargs...)

@deprecate addRow!(args...; nargs...) addrow!(args...; nargs...)


text:
using Documenter, DFWannier

makedocs(sitename="DFWannier")


text:
################################################################################
#
#  Commutator algebra
#
################################################################################

export is_GLZ_conjugate

# Given a square matrix A, determine the algebra C_A = {X | XA = AX }
# and its semisimple reduction
#
# This is mainly for dealing with the isomorphism C_A \cong \End_Q[x](Q^n)
# and the surjection to the semisimplification.
mutable struct CommutatorAlgebra
  A::fmpq_mat
  T::Generic.MatSpaceElem{fmpq_poly}
  Tinv::Generic.MatSpaceElem{fmpq_poly}
  el::Vector{fmpq_poly}
  invariant_factors::Vector{Vector{fmpq_poly}}
  invariant_factors_factored::Vector{Vector{fmpq_poly}}
  invariant_factors_grouped::Vector{Tuple{fmpq_poly,
                                          AnticNumberField,
                                          Vector{Tuple{Int, Int, Int}}}}
  invariant_factors_grouped_grouped::Vector{Vector{Tuple{Int, Vector{Tuple{Int, Int}}}}}
  irreducible_factors::Vector{Tuple{fmpq_poly,
                                    AnticNumberField,
                                    Vector{Tuple{Int, Int}}}}

  function CommutatorAlgebra(A)
    z = new()
    z.A = A
    return z
  end
end

matrix(C::CommutatorAlgebra) = C.A

dim(C::CommutatorAlgebra) = nrows(matrix(C))

function _compute_decomposition!(C::CommutatorAlgebra)
  A = matrix(C)
  Qx = Hecke.Globals.Qx
  x = gen(Qx)
  Ax = x - change_base_ring(Qx, A)
  S, _ , T = snf_with_transform(Ax)
  n = nrows(A)
  el = diagonal(S)
  Tinv = inv(T)
  C.el = el
  C.T = T
  C.Tinv = Tinv

  # Consistency check
  @hassert :Conjugacy 1 begin
    for i in 1:10
      _w = [rand(Qx, 1:5, 1:5) % C.el[i] for i in 1:n]
      v = Hecke._first_map_backward(_w, C)
      @assert Hecke._first_map_forward(v, C) == _w

      v = matrix(Hecke.Globals.QQ, 1, n, [rand(-10:10) for i in 1:n])
      w = Hecke._first_map_forward(v, C)
      @assert Hecke._first_map_backward(w, C) == v
    end
    true
  end

  invariant_factors = Vector{Vector{fmpq_poly}}()

  invariant_factors_factored = Vector{Vector{fmpq_poly}}()

  invariant_factors_grouped = Vector{Tuple{fmpq_poly,
                                           AnticNumberField,
                                           Vector{Tuple{Int, Int, Int}}}}()

  invariant_factors_grouped_grouped = Vector{Tuple{Int,
                                                   Vector{Tuple{Int, Int}}}}[]

  for i in 1:length(C.el)
    fac = factor(C.el[i])::Fac{fmpq_poly}
    inv_fac = Vector{fmpq_poly}()
    factored = Vector{fmpq_poly}()
    j = 1
    for (p, e) in fac
      push!(factored, inv(leading_coefficient(p)) * p^e)

      k = findfirst(x -> isequal(x[1], p), invariant_factors_grouped)
      if k isa Int
        push!(invariant_factors_grouped[k][3], (i, j, e))
      else
        K,  = NumberField(p, cached = false)
        push!(invariant_factors_grouped, (p, K, Tuple{Int, Int, Int}[(i, j, e)]))
      end

      push!(inv_fac, (inv(leading_coefficient(p)) * p)^e)
      j += 1
    end
    push!(invariant_factors_factored, factored)
  end

  for l in 1:length(invariant_factors_grouped)
    D = Tuple{Int, Vector{Tuple{Int, Int}}}[]
    M = invariant_factors_grouped[l][3]
    for (i, j, e) in M
      o = findfirst(x -> x[1] == e, D)
      if o == nothing
        push!(D, (e, Tuple{Int, Int}[(i, j)]))
      else
        push!(D[o][2], (i, j))
      end
    end
    push!(invariant_factors_grouped_grouped, D)
  end

  C.invariant_factors_factored = invariant_factors_factored

  C.invariant_factors_grouped = invariant_factors_grouped

  C.invariant_factors_grouped_grouped = invariant_factors_grouped_grouped

  @hassert :Conjugacy 1 begin
    for i in 1:10
      _w = fmpq_poly[rand(Qx, 1:5, 1:5) % C.el[i] for i in 1:n]
      v = Hecke._second_map_forward(_w, C)
      @assert Hecke._second_map_backward(v, C) == _w
    end

    for i in 1:10
      _w = [fmpq_poly[rand(Qx, 1:5, 1:5) % C.invariant_factors_factored[i][j]
               for j in 1:length(C.invariant_factors_factored[i])]
               for i in 1:dim(C)]
      _v = Hecke._third_map_forward(_w, C)
      @assert Hecke._third_map_backward(_v, C) == _w
    end
    true
  end

  return C
end

function _first_map_forward(w::fmpq_mat, C::CommutatorAlgebra)
  v = change_base_ring(Hecke.Globals.Qx, w) * C.T
  return elem_type(base_ring(v))[v[1, i] % C.el[i] for i in 1:dim(C)]
end

function _first_map_backward(v::Vector{fmpq_poly}, C::CommutatorAlgebra)
  A = matrix(C)
  n = dim(C)
  _w = matrix(Hecke.Globals.Qx, 1, n, v)
  w = _w * C.Tinv
  z = zero_matrix(Hecke.Globals.QQ, 1, n)
  for i in 1:n
    B = w[i](A)
    for j in 1:n
      z[1, j] += B[i, j]
    end
  end
  return z
end

function _second_map_forward(v::Vector{fmpq_poly}, C::CommutatorAlgebra)
  z = Vector{Vector{fmpq_poly}}()
  for i in 1:length(v)
    if length(C.invariant_factors_factored[i]) == 0
      push!(z, fmpq_poly[])
    else
      push!(z, fmpq_poly[v[i] % C.invariant_factors_factored[i][j]
                           for j in 1:length(C.invariant_factors_factored[i])])
    end
  end
  return z
end

function _second_map_backward(v::Vector{Vector{fmpq_poly}}, C::CommutatorAlgebra)
  n = dim(C)
  w = Vector{fmpq_poly}()
  for i in 1:n
    if length(C.invariant_factors_factored[i]) == 0
      push!(w, zero(Hecke.Globals.Qx))
    else
      push!(w, crt(v[i], C.invariant_factors_factored[i]))
    end
  end
  return w
end

function _third_map_forward(v::Vector{Vector{fmpq_poly}}, C::CommutatorAlgebra)
  z = Vector{Vector{fmpq_poly}}(undef, length(C.invariant_factors_grouped))
  for l in 1:length(C.invariant_factors_grouped)
    zz = fmpq_poly[]
    for (i, j) in C.invariant_factors_grouped[l][3]
      push!(zz, v[i][j])
    end
    z[l] = zz
  end
  return z
end

function _third_map_backward(v::Vector{Vector{fmpq_poly}}, C::CommutatorAlgebra)
  z = Vector{Vector{fmpq_poly}}(undef, dim(C))
  for i in 1:dim(C)
    z[i] = Vector{fmpq_poly}(undef, length(C.invariant_factors_factored[i]))
  end

  for l in 1:length(C.invariant_factors_grouped)
    k = 1
    for (i, j) in C.invariant_factors_grouped[l][3]
      z[i][j] = v[l][k]
      k += 1
    end
  end
  return z
end

function _std_basis_vector(C::CommutatorAlgebra, i::Int, j::Int)
  z = Vector{Vector{fmpq_poly}}()
  @assert 1 <= i <= length(C.invariant_factors_grouped)
  @assert 1 <= j <= length(C.invariant_factors_grouped[i][3])
  for k in 1:length(C.invariant_factors_grouped)
    push!(z, zeros(Hecke.Globals.Qx, length(C.invariant_factors_grouped[k][3])))
  end
  z[i][j] = one(Hecke.Globals.Qx)
  return z
end

function _from_std_basis(v::Vector{Vector{fmpq_poly}}, C::CommutatorAlgebra)
  return _first_map_backward(
            _second_map_backward(_third_map_backward(v, C), C), C)
end

function _to_std_basis(v::fmpq_mat, C::CommutatorAlgebra)
  return _third_map_forward(_second_map_forward(_first_map_forward(v, C), C), C)
end

function _decomposition_type(C::CommutatorAlgebra)
  l = length(C.invariant_factors_grouped)
  fields = AnticNumberField[]
  degs = Int[]
  for i in 1:l
    for (e, inds) in C.invariant_factors_grouped_grouped[i]
      push!(fields, C.invariant_factors_grouped[i][2])
      push!(degs, length(inds))
    end
  end

  return fields, degs
end

function _induce_action(C::CommutatorAlgebra, M)
  res = dense_matrix_type(fmpq_poly)[]
  @hassert :Conjugacy 1 M * matrix(C) == matrix(C) * M
  for i in 1:length(C.invariant_factors_grouped)
    n = length(C.invariant_factors_grouped[i][3])
    z = zero_matrix(Hecke.Globals.Qx, n, n)
    for j in 1:length(C.invariant_factors_grouped[i][3])
      s = _std_basis_vector(C, i, j)
      v = _from_std_basis(s, C)
      w = v * M
      ww = _to_std_basis(w, C)
      for k in 1:n
        z[j, k] = ww[i][k]
      end
    end
    push!(res, z)
  end
  return res
end

function _induce_action_mod(C::CommutatorAlgebra, N)
  res = dense_matrix_type(nf_elem)[]
  ac = _induce_action(C, N)
  for i in 1:length(C.invariant_factors_grouped)
    z = ac[i]
    k = 1
    for (e, inds) in C.invariant_factors_grouped_grouped[i]
      l = length(inds)
      zz = sub(z, (k):(k + l - 1), (k):(k + l - 1))
      k = k + l
      push!(res, map_entries(C.invariant_factors_grouped[i][2], zz))
    end
  end
  return res
end

################################################################################
#
#  General approach
#
################################################################################

@doc doc"""
    is_GLZ_conjugate(A::MatElem, B::MatElem) -> Bool, MatElem

Given two integral or rational matrices, determine whether there exists an
invertible integral matrix $T$ with $TA = BT$. If true, the second argument
is such a matrix $T$. Otherwise, the second argument is unspecified.

```jldoctest
julia> A = matrix(ZZ, 4, 4, [ 0, 1,  0, 0,
                             -4, 0,  0, 0,
                              0, 0,  0, 1,
                              0, 0, -4, 0]);

julia> B = matrix(ZZ, 4, 4,  [ 0, 1,  4,  0,
                              -4, 0,  0, -4,
                               0, 0,  0,  1,
                               0, 0, -4,  0]);

julia> fl, T = is_GLZ_conjugate(A, B);

julia> isone(abs(det(T))) && T * A == B * T
true
```
"""
is_GLZ_conjugate(A::Union{fmpz_mat, fmpq_mat}, B::Union{fmpz_mat, fmpq_mat})

function is_GLZ_conjugate(A::fmpz_mat, B::fmpz_mat)
  AQ = change_base_ring(FlintQQ, A)
  BQ = change_base_ring(FlintQQ, B)
  return _isGLZ_conjugate_integral(AQ, BQ)
end

function is_GLZ_conjugate(A::fmpq_mat, B::fmpq_mat)
  d = lcm(denominator(A), denominator(B))
  return _isGLZ_conjugate_integral(d*A, d*B)
end

function _isGLZ_conjugate_integral(A::fmpq_mat, B::fmpq_mat)
  if nrows(A) != nrows(B)
    return false, zero_matrix(FlintZZ, 0, 0)
  end

  if A == B
    return true, identity_matrix(FlintZZ, nrows(A))
  end

  CA, TA = rational_canonical_form(A)
  CB, TB = rational_canonical_form(B)

  if CA != CB
    return false, zero_matrix(FlintZZ, 0, 0)
  end

  Z = CommutatorAlgebra(A)
  _compute_decomposition!(Z)
  Ks, ns = _decomposition_type(Z)
  AA = _matrix_algebra(Ks, ns)
  O = _basis_of_integral_commutator_algebra(A, A)
  I = _basis_of_integral_commutator_algebra(A, B)
  @vprint :Conjugacy 1 "Algebra has dimension $(length(O))\n"
  @vprint :Conjugacy 1 "Semisimple quotient has dimension $(dim(AA))\n"
  ordergens = elem_type(AA)[]
  idealgens = elem_type(AA)[]
  dec = decompose(AA)

  _C = inv(TB) * TA

  @assert _C * A == B * _C

  invC = inv(_C)

  for bb in O
    b = _induce_action_mod(Z, bb)
    z = zero(AA)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      BB, mB = dec[i]::Tuple{AlgAss{fmpq},
                             AbsAlgAssMor{AlgAss{fmpq},AlgAss{fmpq},fmpq_mat}}
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = BB.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(BB))
    end
    push!(ordergens, z)
  end

  for bb in I
    b = _induce_action_mod(Z, invC * bb)
    z = zero(AA)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      BB, mB = dec[i]::Tuple{AlgAss{fmpq},
                             AbsAlgAssMor{AlgAss{fmpq},AlgAss{fmpq},fmpq_mat}}
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = BB.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(BB))
    end
    push!(idealgens, z)
  end

  OO = Order(AA, ordergens)
  OI = ideal_from_lattice_gens(AA, idealgens)
  @hassert :Conjugacy 1 OO == right_order(OI)
  @vprint :Conjugacy 1 "Testing if ideal is principal...\n"
  fl, y = _isprincipal(OI, OO, :right)::Tuple{Bool,
                                              AlgAssElem{fmpq,AlgAss{fmpq}}}

  if !fl
    return false, zero_matrix(FlintZZ, 0, 0)
  end

  @hassert :Conjugacy 1 y * OO == OI

  # I know invC * I maps surjectively onto OI
  # Let's lift the generator y

  d = denominator(OI, OO)

  Y = zero_matrix(FlintZZ, length(idealgens), dim(AA))
  for i in 1:length(idealgens)
    cc = coordinates(OO(d * idealgens[i]))
    @assert length(cc) == dim(AA)
    for j in 1:length(cc)
      Y[i, j] = cc[j]
    end
  end

  YY = matrix(FlintZZ, 1, dim(AA), coordinates(OO(d * y)))

  fl, vv = can_solve_with_solution(Y, YY, side = :left)
  @assert fl
  yy = zero_matrix(FlintQQ, nrows(A), nrows(A))
  for i in 1:length(vv)
    yy = yy + vv[1, i] * (invC * I[i])
  end

  T = _C * yy

  @assert abs(denominator(T)) == 1
  @assert T * A == B * T
  @assert abs(det(T)) == 1

  return fl, map_entries(FlintZZ, T)
end

################################################################################
#
#  Computation of commutator algebras
#
################################################################################

_basis_of_commutator_algebra(A) = _basis_of_commutator_algebra(A, A)

_basis_of_commutator_algebra(As::Vector) = _basis_of_commutator_algebra(A, A)

function _basis_of_integral_commutator_algebra(A::fmpq_mat, B::fmpq_mat)
  # Compute using saturation
  @assert isone(denominator(A))
  @assert isone(denominator(A))
  linind = transpose(LinearIndices(size(A)))
  n = nrows(A)
  z = zero_matrix(FlintQQ, n^2, n^2)
  for i in 1:n
    for j in 1:n
      for k in 1:n
        z[linind[i, j], linind[i, k]] = FlintZZ(z[linind[i, j], linind[i, k]] +
                                                                        A[k, j])
        z[linind[i, j], linind[k, j]] = FlintZZ(z[linind[i, j], linind[k, j]] -
                                                                        B[i, k])
      end
    end
  end
  r, K = right_kernel(z)
  KK = change_base_ring(FlintZZ, denominator(K) * K)
  KK = transpose(saturate(transpose(KK)))
  res = typeof(A)[]
  for k in 1:ncols(K)
    cartind = CartesianIndices(size(A))
    M = zero_matrix(base_ring(A), nrows(A), ncols(A))
    for l in 1:n^2
      i1, j1 = cartind[l].I
      M[j1, i1] = KK[l, k]
    end
    @assert M * A == B * M
    push!(res, M)
  end
  return res
end

function _basis_of_commutator_algebra(A::MatElem, B::MatElem)
  linind = transpose(LinearIndices(size(A)))
  n = nrows(A)
  z = zero_matrix(base_ring(A), n^2, n^2)
  for i in 1:n
    for j in 1:n
      for k in 1:n
        z[linind[i, j], linind[i, k]] = FlintZZ(z[linind[i, j], linind[i, k]] +
                                                                        A[k, j])
        z[linind[i, j], linind[k, j]] = FlintZZ(z[linind[i, j], linind[k, j]] -
                                                                        B[i, k])
      end
    end
  end
  r, K = right_kernel(z)
  res = typeof(A)[]
  for k in 1:ncols(K)
    cartind = cartesian_product_iterator([1:x for x in size(A)], inplace = true)
    M = zero_matrix(base_ring(A), nrows(A), ncols(A))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    @assert M * A == B * M
    push!(res, M)
  end
  return res
end

function _basis_of_commutator_algebra(As::Vector, Bs::Vector)
  linind = transpose(LinearIndices(size(As[1])))
  n = nrows(As[1])
  zz = zero_matrix(base_ring(As[1]), 0, n^2)
  for (A, B) in zip(As, Bs)
    z = zero_matrix(base_ring(A), n^2, n^2)
    for i in 1:n
      for j in 1:n
        for k in 1:n
          z[linind[i, j], linind[i, k]] = z[linind[i, j], linind[i, k]] +
                                                                        A[k, j]
          z[linind[i, j], linind[k, j]] = z[linind[i, j], linind[k, j]] -
                                                                        B[i, k]
        end
      end
    end
    zz = vcat(zz, z)
  end
  r, K = right_kernel(zz)
  res = eltype(As)[]
  for k in 1:ncols(K)
    cartind = cartesian_product_iterator([1:x for x in size(As[1])],
                                         inplace = true)
    M = zero_matrix(base_ring(As[1]), nrows(As[1]), ncols(As[1]))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    push!(res, M)
  end
  return res
end

function _basis_of_integral_commutator_algebra(As::Vector{fmpq_mat},
                                               Bs::Vector{fmpq_mat})
  # Compute using saturation
  @assert all(x -> isone(denominator(x)), As)
  @assert all(x -> isone(denominator(x)), Bs)
  linind = transpose(LinearIndices(size(As[1])))
  n = nrows(As[1])
  zz = zero_matrix(base_ring(As[1]), 0, n^2)
  for (A, B) in zip(As, Bs)
    z = zero_matrix(FlintQQ, n^2, n^2)
    for i in 1:n
      for j in 1:n
        for k in 1:n
          z[linind[i, j], linind[i, k]] = z[linind[i, j], linind[i, k]] + A[k, j]
          z[linind[i, j], linind[k, j]] = z[linind[i, j], linind[k, j]] - B[i, k]
        end
      end
    end
    zz = vcat(zz, z)
  end
  r, K = right_kernel(zz)
  KK = change_base_ring(FlintZZ, denominator(K) * K)
  KK = transpose(saturate(transpose(KK)))
  res = typeof(As[1])[]
  for k in 1:ncols(K)
    cartind = CartesianIndices(size(As[1]))
    M = zero_matrix(base_ring(As[1]), n, n)
    for l in 1:n^2
      i1, j1 = cartind[l].I
      M[j1, i1] = KK[l, k]
    end
    for i in 1:length(As)
      M * As[i] == Bs[i] * M
    end
    push!(res, M)
  end
  return res
end

################################################################################
#
#  Simultaneous conjugacy
#
################################################################################

#function _isconjugated_probabilistic(a::Vector{fmpz_mat}, b::Vector{fmpz_mat})
#  return isconjugated_probabilistic(map(x -> map(QQ, x), a),
#                                    map(x -> map(QQ, x), b))
#end
#
#function _isconjugated_probabilistic(a::Vector{fmpq_mat}, b::Vector{fmpq_mat})
#  B = _basis_of_commutator_algebra(a, b)
#  l = length(B)
#  for i in 1:50
#    c = sum(c * b for (c, b) in zip(rand(-1:1, l), B))
#    if !iszero(det(c))
#      return true, c
#    end
#  end
#  return false, zero_matrix(FlintQQ, 0, 0)
#end
#
#function _isGLZ_conjugate(A::Vector{fmpz_mat}, B::Vector{fmpz_mat})
#  return __isGLZ_conjugate(map(x -> change_base_ring(FlintQQ, x), A),
#                          map(x -> change_base_ring(FlintQQ, x), B))
#end
#
#function _isGLZ_conjugate(A::Vector{fmpq_mat}, B::Vector{fmpq_mat})
#  d1 = lcm(fmpz[denominator(x) for x in A])
#  d2 = lcm(fmpz[denominator(x) for x in B])
#  d = lcm(d1, d2)
#  return __isGLZ_conjugate(d .* A, d .* B)
#end
#
#function __isGLZ_conjugate(A::Vector{fmpq_mat}, B::Vector{fmpq_mat})
#
#  if A == B
#    return true, identity_matrix(FlintQQ, nrows(A[1]))
#  end
#  O = _basis_of_integral_commutator_algebra(A, A)
#  I = _basis_of_integral_commutator_algebra(A, B)
#  AA = matrix_algebra(FlintQQ, map(x -> map(FlintQQ, x), O))
#  ordergens = elem_type(AA)[]
#  idealgens = elem_type(AA)[]
#
#  fl, _C = _isconjugated_probabilistic(A, B)
#  if !fl
#    return false, zero_matrix(FlintQQ, 0, 0)
#  end
#
#  @assert all(_C * map(QQ, A[i]) == map(QQ, B[i]) * _C for i in 1:length(A))
#
#  invC = inv(_C)
#
#  @info "Dimension of the commutator algebra: $(dim(AA))"
#
#  OO = Order(AA, map(x -> AA(map(QQ, x)), O))
#  OI = ideal_from_lattice_gens(AA, map(x -> AA(invC * map(QQ, x)), I))
#  @assert OO == right_order(OI)
#
#  J = radical(AA)
#
#  # Now lets do this
#  if iszero(J)
#    # Semisimple case
#    if dim(AA) == nrows(A[1])^2
#      # this is a full matrix algebra
#      dec = decompose(AA)
#      B, mB = dec[1]
#      A.isomorphic_full_matrix_algebra = A, inv(mB)
#      fl, y = _isprincipal(OI, OO, :right)::Tuple{Bool, AlgAssElem{fmpq,typeof(AA)}}
#      yy = elem_in_algebra(y)
#    elseif is_commutative(AA)
#      @info "Algebra is commutative"
#      OI.order = OO
#      d = denominator(OI, OO)
#      # Fix this upstream!
#      fl, y = is_principal(d * OI)
#      if !fl
#        return false, zero_matrix(QQ, 0, 0)
#      end
#      yy = inv(QQ(d)) * elem_in_algebra(y)
#    else
#      error("Not implemented")
#    end
#  else
#    # Pass to the semisimple quotient
#    J = radical(AA)
#    S, AtoS = quo(AA, J)
#    @info "Semisimple quotient has dimension $(dim(S))"
#    !is_commutative(S) && error("Semisimple quotient must be commutative")
#    IS = ideal_from_lattice_gens(S, [AtoS(b) for b in basis(OI)])
#    OS = Order(S, [AtoS(elem_in_algebra(b)) for b in basis(OO)])
#    @info "Algebra is commutative"
#    IS.order = OS
#    d = denominator(IS, OS)
#    # Fix this upstream!
#    fl, yy = is_principal(d * IS)
#    if !fl
#      return false, zero_matrix(QQ, 0, 0)
#    end
#    yyy = inv(QQ(d)) * elem_in_algebra(yy)
#    # Now I have to lift this
#    # I cannot just lift, I need a preimage in OI
#    d = denominator(OI, OO)
#
#    Y = zero_matrix(FlintQQ, dim(AA), dim(S))
#    OIbasis = basis(OI)
#    for i in 1:dim(AA)
#      cc = coefficients(AtoS(OIbasis[i]))
#      for j in 1:length(cc)
#        Y[i, j] = cc[j]
#      end
#    end
#
#    YY = matrix(FlintQQ, 1, dim(S), coefficients(yyy))
#    # I look for a integral solution, but the matrices are rational ..
#    d = lcm(denominator(Y), denominator(YY))
#    fl, vv = can_solve_with_solution(map(ZZ, d*Y), map(ZZ, d*YY), side = :left)
#    @assert fl
#    yy = sum(vv[i] * OIbasis[i] for i in 1:dim(AA))
#  end
#
#  if !fl
#    return false, zero_matrix(FlintQQ, 0, 0)
#  end
#
#  @assert yy * OO == OI
#
#  D = _C * matrix(yy)
#  @assert all(D * map(QQ, A[i]) == map(QQ, B[i]) * D for i in 1:length(A))
#  @assert abs(det(D)) == 1
#
#  return fl, _C * matrix(yy)
#
#  ## I know invC * I maps surjectively onto OI
#  ## Let's let the generator y
#
#  #d = denominator(OI, OO)
#
#  #Y = zero_matrix(FlintZZ, length(idealgens), dim(AA))
#  #for i in 1:length(idealgens)
#  #  cc = coordinates(OO(d * idealgens[i]))
#  #  @assert length(cc) == dim(AA)
#  #  for j in 1:length(cc)
#  #    Y[i, j] = cc[j]
#  #  end
#  #end
#
#  #YY = matrix(FlintZZ, 1, dim(AA), coordinates(OO(d * y)))
#
#  #fl, vv = can_solve_with_solution(Y, YY, side = :left)
#  #@assert fl
#  #yy = zero_matrix(FlintQQ, nrows(A), nrows(A))
#  #for i in 1:length(vv)
#  #  yy = yy + vv[1, i] * (invC * I[i])
#  #end
#
#  #T = _C * yy
#
#  #@assert abs(denominator(T)) == 1
#  #@assert T * A == B * T
#  #@assert abs(det(T)) == 1
#
#  #return fl, T
#end


text:
# Copyright 2015, 2016 Martin Holters
# See accompanying license file.

using ACME
using Base.Test

tv, ti = ACME.topomat(sparse([1 -1 1; -1 1 -1]))
@test tv*ti'==spzeros(2,1)

# Pathological cases for topomat:
# two nodes, one loop branch (short-circuited) -> voltage==0, current arbitrary
@test ACME.topomat(spzeros(Int, 2, 1)) == (speye(1), spzeros(0, 1))
# two nodes, one branch between them -> voltage arbitrary, current==0
@test ACME.topomat(sparse([1,2], [1,1], [1,-1])) == (spzeros(0, 1), speye(1))

let circ = Circuit()
    model=DiscreteModel(circ, 1.)
    @test run!(model, zeros(0, 20)) == zeros(0, 20)
end

let circ = Circuit(), r = resistor(0)
    connect!(circ, r[1], r[2])
    model = DiscreteModel(circ, 1.)
    @test run!(model, zeros(0, 20)) == zeros(0, 20)
end

let circ = Circuit(), r = resistor(0), probe = currentprobe()
    connect!(circ, r[1], probe[:+])
    connect!(circ, r[2], probe[:-])
    orig_stderr = STDERR
    rd, wr = redirect_stderr()
    model = DiscreteModel(circ, 1.)
    # should warn because output is indeterminate
    @test !isempty(search(convert(ASCIIString, readavailable(rd)), "WARNING"))
    redirect_stderr(orig_stderr)
end

for num = 1:50
    let ps = rand(4, num)
        t = ACME.KDTree(ps)
        for i in 1:size(ps)[2]
            idx = ACME.indnearest(t, ps[:,i])[1]
            @test ps[:,i] == ps[:,idx]
        end
    end
end

let ps = rand(6, 10000)
    t = ACME.KDTree(ps)
    p = rand(6)
    best_p = ps[:,indmin(sumabs2(broadcast(-, ps, p),1))]
    idx = ACME.indnearest(t, p)[1]
    @test_approx_eq sumabs2(p - best_p) sumabs2(p - ps[:, idx])
end

let nleq = ACME.ParametricNonLinEq((res, J, Jp, p, z) ->
    begin
        res[1] = z[1]^2 - 1 + p[1]
        J[1,1] = 2*z[1]
        Jp[1,1] = 1
    end, 1, 1)
    solver = ACME.HomotopySolver{ACME.SimpleSolver}(nleq, [0.0], [1.0])
    ACME.solve(solver, [-0.5 + rand()])
    @test ACME.hasconverged(solver)
    ACME.solve(solver, [1.5 + rand()])
    @test !ACME.hasconverged(solver)
end

# simple circuit: resistor and diode in series, driven by constant voltage,
# chosen such that a prescribe current flows
let i = 1e-3, r=10e3, is=1e-12
    v_r = i*r
    v_d = 25e-3 * log(i/is+1)
    vsrc = voltagesource(v_r + v_d)
    r1 = resistor(r)
    d = diode(is=is)
    vprobe = voltageprobe()
    circ = Circuit()
    add!(circ, vsrc)
    add!(circ, r1, d)
    connect!(circ, vsrc[:+], :vcc)
    connect!(circ, vsrc[:-], :gnd)
    connect!(circ, r1[1], :vcc)
    connect!(circ, d[:-], vprobe[:-], :gnd)
    connect!(circ, r1[2], d[:+], vprobe[:+])
    model = DiscreteModel(circ, 1.)
    y = run!(model, zeros(0, 1))
    @test_approx_eq_eps y[1] v_d 1e-6
end

function checksteady!(model)
    x_steady = steadystate!(model)
    ACME.set_resabs2tol!(model.solver, 1e-25)
    run!(model, zeros(1, 1))
    @test_approx_eq model.x x_steady
end

include("../examples/sallenkey.jl")
let model=sallenkey()
    println("Running sallenkey")
    y = run!(model, sin(2π*1000/44100*(0:44099)'))
    @test size(y) == (1,44100)
    # TODO: further validate y

    # cannot check steady state: steadystate() does not work for matrix A having
    # eigenvalue 1
end

include("../examples/diodeclipper.jl")
let model=diodeclipper()
    println("Running diodeclipper")
    y = run!(model, sin(2π*1000/44100*(0:44099)'))
    @test size(y) == (1,44100)
    # TODO: further validate y
    checksteady!(model)
end

include("../examples/birdie.jl")
let model=birdie(vol=0.8)
    ACME.solve(model.solver, [0.003, -0.0002])
    @assert ACME.hasconverged(model.solver)
    println("Running birdie with fixed vol")
    y = run!(model, sin(2π*1000/44100*(0:44099)'))
    @test size(y) == (1,44100)
    # TODO: further validate y
    checksteady!(model)
end
let model=birdie()
    println("Running birdie with varying vol")
    y = run!(model, [sin(2π*1000/44100*(0:44099).'); linspace(1,0,44100).'])
    @test size(y) == (1,44100)
    # TODO: further validate y
end

include("../examples/superover.jl")
let model=superover(drive=1.0, tone=1.0, level=1.0)
    println("Running superover with fixed potentiometer values")
    y = run!(model, sin(2π*1000/44100*(0:44099)'))
    @test size(y) == (1,44100)
    # TODO: further validate y
    checksteady!(model)
end
let model=superover()
    println("Running superover with varying potentiometer values")
    y = run!(model, [sin(2π*1000/44100*(0:999)'); linspace(1,0,1000).'; linspace(0,1,1000).'; linspace(1,0,1000).'])
    @test size(y) == (1,1000)
    # TODO: further validate y
end


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

mutable struct LocationMetadata <: SwaggerModel
    regionType::Any # spec type: Union{ Nothing, String } # spec name: regionType
    regionCategory::Any # spec type: Union{ Nothing, String } # spec name: regionCategory
    geographyGroup::Any # spec type: Union{ Nothing, String } # spec name: geographyGroup
    longitude::Any # spec type: Union{ Nothing, String } # spec name: longitude
    latitude::Any # spec type: Union{ Nothing, String } # spec name: latitude
    physicalLocation::Any # spec type: Union{ Nothing, String } # spec name: physicalLocation
    pairedRegion::Any # spec type: Union{ Nothing, Vector{PairedRegion} } # spec name: pairedRegion

    function LocationMetadata(;regionType=nothing, regionCategory=nothing, geographyGroup=nothing, longitude=nothing, latitude=nothing, physicalLocation=nothing, pairedRegion=nothing)
        o = new()
        validate_property(LocationMetadata, Symbol("regionType"), regionType)
        setfield!(o, Symbol("regionType"), regionType)
        validate_property(LocationMetadata, Symbol("regionCategory"), regionCategory)
        setfield!(o, Symbol("regionCategory"), regionCategory)
        validate_property(LocationMetadata, Symbol("geographyGroup"), geographyGroup)
        setfield!(o, Symbol("geographyGroup"), geographyGroup)
        validate_property(LocationMetadata, Symbol("longitude"), longitude)
        setfield!(o, Symbol("longitude"), longitude)
        validate_property(LocationMetadata, Symbol("latitude"), latitude)
        setfield!(o, Symbol("latitude"), latitude)
        validate_property(LocationMetadata, Symbol("physicalLocation"), physicalLocation)
        setfield!(o, Symbol("physicalLocation"), physicalLocation)
        validate_property(LocationMetadata, Symbol("pairedRegion"), pairedRegion)
        setfield!(o, Symbol("pairedRegion"), pairedRegion)
        o
    end
end # type LocationMetadata

const _property_map_LocationMetadata = Dict{Symbol,Symbol}(Symbol("regionType")=>Symbol("regionType"), Symbol("regionCategory")=>Symbol("regionCategory"), Symbol("geographyGroup")=>Symbol("geographyGroup"), Symbol("longitude")=>Symbol("longitude"), Symbol("latitude")=>Symbol("latitude"), Symbol("physicalLocation")=>Symbol("physicalLocation"), Symbol("pairedRegion")=>Symbol("pairedRegion"))
const _property_types_LocationMetadata = Dict{Symbol,String}(Symbol("regionType")=>"String", Symbol("regionCategory")=>"String", Symbol("geographyGroup")=>"String", Symbol("longitude")=>"String", Symbol("latitude")=>"String", Symbol("physicalLocation")=>"String", Symbol("pairedRegion")=>"Vector{PairedRegion}")
Base.propertynames(::Type{ LocationMetadata }) = collect(keys(_property_map_LocationMetadata))
Swagger.property_type(::Type{ LocationMetadata }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LocationMetadata[name]))}
Swagger.field_name(::Type{ LocationMetadata }, property_name::Symbol) =  _property_map_LocationMetadata[property_name]

const _allowed_LocationMetadata_regionType = ["Physical", "Logical"]

const _allowed_LocationMetadata_regionCategory = ["Recommended", "Other"]

function check_required(o::LocationMetadata)
    true
end

function validate_property(::Type{ LocationMetadata }, name::Symbol, val)
    if name === Symbol("regionType")
        Swagger.validate_param(name, "LocationMetadata", :enum, val, _allowed_LocationMetadata_regionType)
    end
    if name === Symbol("regionCategory")
        Swagger.validate_param(name, "LocationMetadata", :enum, val, _allowed_LocationMetadata_regionCategory)
    end
end


text:
function run_opf(
      parent_folder,
      folder,
      objective,
      gen_areas,
      area_transfer,
      gen_zones,
      zone_transfer,
      results_dict,
      fix_f_override,
      direct_pq,
      master_subnet,
      suffix,
      print_results;
      override_param::Dict=Dict(),
      start_vals=Dict("sn"=>Dict()),
      no_converter_loss::Bool=false,
      regularize_f::Float64=0.0,
      ipopt_max_iter::Int64=10000,
      output_to_files=output_to_files
   )
   (result, res_summary, solution_pm, binding_cnstr) = multifrequency_opf(
      "$parent_folder/$(folder)", objective,
      gen_areas=gen_areas, area_interface=area_transfer,
      gen_zones=gen_zones, zone_interface=zone_transfer,
      print_results=print_results,
      fix_f_override=fix_f_override, direct_pq=direct_pq,
      master_subnet=1, suffix=suffix, override_param=override_param,
      start_vals=start_vals, no_converter_loss=no_converter_loss,
      output_to_files=output_to_files,
      regularize_f=regularize_f,
      ipopt_max_iter=ipopt_max_iter
      )
   results_dict_out = deepcopy(results_dict)
   for (label, value) in result
      # println("label: $label")
      if !(label in keys(results_dict_out))
         results_dict_out[label] = Dict{String, Any}()
      end
      # Comment '|| true' to prevent saving outputs when not solved to optimum
      if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED] || true
         # results_dict_out[label][rstrip(folder,'/')*suffix] = value
         results_dict_out[label][rstrip(folder,'/')] = value
      end
   end
   return (result, res_summary, solution_pm, results_dict_out)
end


text:
# Runs the model to SS, and issue the pertrubation -
function Simulation(time_start,time_stop,time_step_size,data_dictionary)

  # Run the model to steady-state -
  steady_state_array = EstimateSteadyState(0.01,data_dictionary);

  # Ok, reset the TRIGGER and CRAF -
  steady_state_array[53] = 1;
  steady_state_array[54] = 1;

  #@show steady_state_array[37:54]

  # Reset the IC's -
  data_dictionary["initial_condition_array"] = steady_state_array;

  # Solve the model equations -
  (T,X) = SolveBalances(time_start,time_stop,time_step_size,data_dictionary)

  # return my simulation time, and state -
  return (T,X);
end

function EstimateSteadyState(epsilon,data_dictionary)

  initial_condition_vector = data_dictionary["initial_condition_array"];
  ic_array = copy(data_dictionary["initial_condition_array"])
  number_of_states = length(ic_array)

  # Setup loop -
  EPSILON = epsilon;
  TSTART = 0.0;
  Ts = 1.0;
  TSTOP = 1000;
  did_reach_steady_state = false
  while (!did_reach_steady_state)

    # solve the balances -
    (TSIM,X1) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary)

    # Take a few additional steps -
    TNEXT_START = TSTOP+Ts;
    TNEXT_STOP = TNEXT_START+1.0;
    Ts = 0.1;

    # solve the balances again 0
    initial_condition_array = vec(X1[end,:])
    data_dictionary["initial_condition_array"] = initial_condition_array;
    (TSIM,X2) = SolveBalances(TNEXT_START,TNEXT_STOP,Ts,data_dictionary)

    # Find the difference -
    DIFF = norm((X2[end,:] - X1[end,:]));

    # Should we stop -or- go around again?
    if (DIFF<EPSILON)
      did_reach_steady_state = true;
      return (vec(X2[end,:]));
    else

      # No, we did *not* reach steady state ....
      TSTART = TSTOP+Ts
      TSTOP = 1.0 + TSTART;
      Ts = 0.1;

      initial_condition_array = vec(X2[end,:])
      data_dictionary["initial_condition_array"] = initial_condition_array;
    end
  end

  # return
  return XSS;
end


text:
##### SGS quadrature
using FastGaussQuadrature
abstract type SubdomainStatistics end

struct SubdomainMean{IT <: Int} <: SubdomainStatistics
    order::IT
end

struct GaussianQuadrature{IT <: Int} <: SubdomainStatistics
    order::IT
end

struct LogNormalQuadrature{IT <: Int} <: SubdomainStatistics
    order::IT
end

"""
    compute_subdomain_statistics!

at the moment only in the en subdomain
"""

function compute_subdomain_statistics! end

function compute_subdomain_statistics!(
    grid::Grid{FT},
    i,
    q,
    aux,
    params,
    model::SubdomainMean,
) where {FT, IT}
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params
    @inbounds for k in over_elems_real(grid)
        q_tot = q[:q_tot, k, en]
        ts = ActiveThermoState(param_set, q, aux, k, en)
        T = air_temperature(ts)
        q_liq = PhasePartition(ts).liq
        q_vap = q_tot - q_liq
        θ = dry_pottemp(ts)
        if q_liq > 0
            aux[:CF, k] = FT(1)
            aux[:θ_cloudy, k] = θ
            aux[:t_cloudy, k] = T
            aux[:q_tot_cloudy, k] = q_tot
            aux[:q_vap_cloudy, k] = q_vap
        else
            aux[:CF, k] = FT(0)
            aux[:θ_dry, k] = θ
            aux[:q_tot_dry, k] = q_tot
        end
    end
end

function compute_subdomain_statistics!(
    grid::Grid{FT},
    i,
    q,
    aux,
    params,
    model::GaussianQuadrature,
) where {FT, IT}
    gm, en, ud, sd, al = allcombinations(q)
    abscissas, weights = gausshermite(model.order)
    sqrt2 = sqrt(2.0)
    sqpi_inv = 1 / sqrt(pi)
    env_len = 10
    src_len = 6
    @unpack param_set = params

    inner_env = zeros(env_len)
    outer_env = zeros(env_len)
    inner_src = zeros(src_len)
    outer_src = zeros(src_len)
    i_Sqt = collect(1:src_len)
    i_SH = collect(1:src_len)
    i_Sqt_H = collect(1:src_len)
    i_Sqt_qt = collect(1:src_len)
    i_SH_H = collect(1:src_len)
    i_SH_qt = collect(1:src_len)
    i_ql = collect(1:env_len)
    i_T = collect(1:env_len)
    i_θl = collect(1:env_len)
    i_ρ = collect(1:env_len)
    i_cf = collect(1:env_len)
    i_qt_cld = collect(1:env_len)
    i_qt_dry = collect(1:env_len)
    i_T_cld = collect(1:env_len)
    i_T_dry = collect(1:env_len)
    i_rf = collect(1:env_len)

    @inbounds for k in over_elems_real(grid)
        if (
            q[:cv_q_tot, k, i] > eps(FT) &&
            q[:cv_θ_liq, k, i] > eps(FT) &&
            fabs(q[:cv_θ_liq_q_tot, k, i]) > eps(FT) &&
            q[:q_tot, k, i] > eps(FT) &&
            sqrt(q[:cv_q_tot, k, i]) < q[:q_tot, k, i]
        )
            σ_q = sqrt(q[:cv_q_tot, k, i])
            σ_h = sqrt(q[:cv_θ_liq, k, i])
            corr = max(
                min(q[:cv_θ_liq_q_tot, k, i] / max(σ_h * σ_q, 1e-13), 1),
                -1,
            )
            # limit σ_q to prevent negative qt_hat
            σ_qt_lim = (1e-10 - q[:q_tot, k, i]) / (sqrt2 * abscissas[0])
            # walking backwards to assure your q_t will not be smaller than 1e-10
            # TODO - check
            # TODO - change 1e-13 and 1e-10 to some epislon
            σ_q = min(σ_q, σ_qt_lim)
            qt_covar = σ_q * σ_q
            σ_θ_star = sqrt(max(1 - corr * corr, 0.0)) * σ_h

            # clean outer vectors
            outer_src = 0.0 * outer_src
            outer_env = 0.0 * outer_env

            for j_qt in 1:(model.order)
                qt_hat = q[:q_tot, k, i] + sqrt2 * σ_q * abscissas[j_qt]
                μ_θ_star =
                    q[:θ_lid, k, i] + sqrt2 * corr * σ_h * abscissas[j_qt]
                # clean innner vectors
                inner_src = 0.0 * inner_src
                inner_env = 0.0 * inner_env
                for j_θ in 1:(model.order)
                    θ_hat = sqrt2 * σ_θ_star * abscissas[j_θ] + μ_θ_star
                    ts = DummyThermoState(
                        param_set;
                        θ = θ_hat,
                        q_tot = qt_hat,
                        ρ_0 = aux[:ρ_0, k, i],
                        p_0 = aux[:p_0, k, i],
                    )
                    T = air_temperature(ts)
                    q_liq = PhasePartition(ts).liq
                    q_vap = q_tot - q_liq
                    θ = dry_pottemp(ts)

                    # microphysics rain should apply here to update qt, ql, thetal, T
                    # now collect the inner values of
                    inner_env[i_ql] += q_liq * weights[j_θ] * sqpi_inv
                    inner_env[i_T] += T * weights[j_θ] * sqpi_inv
                    inner_env[i_θl] += θ * weights[j_θ] * sqpi_inv
                    inner_env[i_ρ] += ρ * weights[j_θ] * sqpi_inv
                    # rain area fraction
                    if aux[:qr_src, k, i] > 0
                        inner_env[i_rf] += weights[j_θ] * sqpi_inv
                    end
                    # cloudy/dry categories for buoyancy in TKE
                    if q_liq > 0.0
                        inner_env[i_cf] += weights[j_θ] * sqpi_inv
                        inner_env[i_qt_cld] += qt_hat * weights[j_θ] * sqpi_inv
                        inner_env[i_T_cld] += T * weights[j_θ] * sqpi_inv
                    else
                        inner_env[i_qt_dry] += qt_hat * weights[j_θ] * sqpi_inv
                        inner_env[i_T_dry] += T * weights[j_θ] * sqpi_inv
                    end
                    # products for variance and covariance source terms
                    inner_src[i_Sqt] += -qr_src * weights[j_θ] * sqpi_inv
                    inner_src[i_SH] += θl_rain_src * weights[j_θ] * sqpi_inv
                    inner_src[i_Sqt_H] += -qr_src * θ * weights[j_θ] * sqpi_inv
                    inner_src[i_SH_H] +=
                        θl_rain_src * θ * weights[j_θ] * sqpi_inv
                    inner_src[i_Sqt_qt] +=
                        -qr_src * q_tot * weights[j_θ] * sqpi_inv
                    inner_src[i_SH_qt] +=
                        θl_rain_src * q_tot * weights[j_θ] * sqpi_inv
                end

                outer_env[i_ql] += inner_env[i_ql] * weights[j_qt] * sqpi_inv
                outer_env[i_T] += inner_env[i_T] * weights[j_qt] * sqpi_inv
                outer_env[i_θl] += inner_env[i_θl] * weights[j_qt] * sqpi_inv
                outer_env[i_ρ] += inner_env[i_ρ] * weights[j_qt] * sqpi_inv
                outer_env[i_cf] += inner_env[i_cf] * weights[j_qt] * sqpi_inv
                outer_env[i_qt_cld] +=
                    inner_env[i_qt_cld] * weights[j_qt] * sqpi_inv
                outer_env[i_qt_dry] +=
                    inner_env[i_qt_dry] * weights[j_qt] * sqpi_inv
                outer_env[i_T_cld] +=
                    inner_env[i_T_cld] * weights[j_qt] * sqpi_inv
                outer_env[i_T_dry] +=
                    inner_env[i_T_dry] * weights[j_qt] * sqpi_inv
                outer_env[i_rf] += inner_env[i_rf] * weights[j_qt] * sqpi_inv

                outer_src[i_Sqt] += inner_src[i_Sqt] * weights[j_qt] * sqpi_inv
                outer_src[i_SH] += inner_src[i_SH] * weights[j_qt] * sqpi_inv
                outer_src[i_Sqt_H] +=
                    inner_src[i_Sqt_H] * weights[j_qt] * sqpi_inv
                outer_src[i_SH_H] +=
                    inner_src[i_SH_H] * weights[j_qt] * sqpi_inv
                outer_src[i_Sqt_qt] +=
                    inner_src[i_Sqt_qt] * weights[j_qt] * sqpi_inv
                outer_src[i_SH_qt] +=
                    inner_src[i_SH_qt] * weights[j_qt] * sqpi_inv
            end

            # aux[:CF, k, i]           = outer_env[i_cf]
            # aux[:t_cloudy, k, i]     = outer_env[i_T]
            # aux[:q_tot_cloudy, k, i] = outer_env[i_qt_cld]
            # aux[:q_vap_cloudy, k, i] = outer_env[i_qt_cld] - outer_env[i_ql]
            # aux[:q_tot_dry, k, i]    = outer_env[i_qt_dry]
            # ts = TemperatureSHumEquil(param_set, aux[:t_cloudy, k], p_0=aux[:p_0, k, i], aux[:q_tot_cloudy, k, i])
            # aux[:θ_cloudy, k, i] = liquid_ice_pottemp(ts)
            # aux[:θ_dry, k, i] = dry_pottemp(ts)

            aux[:CF, k] = outer_env[i_cf]
            aux[:t_cloudy, k] = outer_env[i_T]
            aux[:q_tot_cloudy, k] = outer_env[i_qt_cld]
            aux[:q_vap_cloudy, k] = outer_env[i_qt_cld] - outer_env[i_ql]
            aux[:q_tot_dry, k] = outer_env[i_qt_dry]
            ts = TemperatureSHumEquil(
                param_set,
                aux[:t_cloudy, k],
                p_0 = aux[:p_0, k, i],
                aux[:q_tot_cloudy, k, i],
            )
            aux[:θ_cloudy, k] = liquid_ice_pottemp(ts)
            aux[:θ_dry, k] = dry_pottemp(ts)
        end
    end
end

function compute_subdomain_statistics!(
    grid::Grid{FT},
    i,
    q,
    aux,
    params,
    model::LogNormalQuadrature,
) where {FT, IT}
    gm, en, ud, sd, al = allcombinations(q)
    abscissas, weights = gausshermite(model.order)
    sqrt2 = sqrt(2.0)
    sqpi_inv = 1 / sqrt(pi)
    env_len = 10
    src_len = 6
    @unpack param_set = params

    inner_env = zeros(env_len)
    outer_env = zeros(env_len)
    inner_src = zeros(src_len)
    outer_src = zeros(src_len)
    i_Sqt = collect(1:src_len)
    i_SH = collect(1:src_len)
    i_Sqt_H = collect(1:src_len)
    i_Sqt_qt = collect(1:src_len)
    i_SH_H = collect(1:src_len)
    i_SH_qt = collect(1:src_len)
    i_ql = collect(1:env_len)
    i_T = collect(1:env_len)
    i_θl = collect(1:env_len)
    i_ρ = collect(1:env_len)
    i_cf = collect(1:env_len)
    i_qt_cld = collect(1:env_len)
    i_qt_dry = collect(1:env_len)
    i_T_cld = collect(1:env_len)
    i_T_dry = collect(1:env_len)
    i_rf = collect(1:env_len)

    @inbounds for k in over_elems_real(grid)
        if (
            q[:cv_q_tot, k, i] > eps(FT) &&
            q[:cv_θ_liq, k, i] > eps(FT) &&
            fabs(q[:cv_θ_liq_q_tot, k, i]) > eps(FT) &&
            q[:q_tot, k, i] > eps(FT) &&
            sqrt(q[:cv_q_tot, k, i]) < q[:q_tot, k, i]
        )
            σ_q = sqrt(
                log(q[:cv_q_tot, k, i] / q[:q_tot, k, i] / q[:q_tot, k, i]) + 1,
            )
            σ_h = sqrt(
                log(q[:cv_θ_liq, k, i] / q[:θ_liq, k, i] / q[:θ_liq, k, i]) + 1,
            )
            # Enforce Schwarz's inequality
            corr = max(
                min(q[:cv_θ_liq_q_tot, k, i] / max(σ_h * σ_q, 1e-13), 1),
                -1,
            )
            sd2_hq = log(
                corr * sqrt(q[:cv_q_tot, k, i] * q[:cv_θ_liq, k, i]) /
                (q[:θl, k, i] * q[:q_tot, k, i]) + 1.0,
            )
            σ_cond_θ_q = sqrt(max(σ_h * σ_h - sd2_hq * sd2_hq / σ_q / σ_q, 0.0))
            μ_q = log(
                q[:q_tot, k, i] * q[:q_tot, k, i] / sqrt(
                    q[:q_tot, k, i] * q[:q_tot, k, i] + q[:cv_q_tot, k, i],
                ),
            )
            μ_θ = log(
                q[:θ_liq, k, i] * q[:θ_liq, k, i] / sqrt(
                    q[:θ_liq, k, i] * q[:θ_liq, k, i] + q[:cv_θ_liq, k, i],
                ),
            )
            # clean outer vectors
            outer_src = 0.0 * outer_src
            outer_env = 0.0 * outer_env

            for j_qt in 1:(model.order)
                qt_hat = exp(μ_q + sqrt2 * σ_q * abscissas[j_q])
                μ_θ_star = μ_θ + sd2_hq / sd_q / sd_q * (log(qt_hat) - μ_q)
                # clean innner vectors
                inner_src = 0.0 * inner_src
                inner_env = 0.0 * inner_env
                for j_θ in 1:(model.order)
                    θ_hat = exp(μ_θ_star + sqrt2 * σ_cond_θ_q * abscissas[j_θ])
                    ts = DummyThermoState(
                        param_set;
                        θ = θ_hat,
                        q_tot = qt_hat,
                        ρ_0 = aux[:ρ_0, k, i],
                        p_0 = aux[:p_0, k, i],
                    )
                    T = air_temperature(ts)
                    q_liq = PhasePartition(ts).liq
                    q_vap = q_tot - q_liq
                    θ = dry_pottemp(ts)

                    # microphysics rain should apply here to update qt, ql, thetal, T
                    # now collect the inner values of
                    inner_env[i_ql] += q_liq * weights[j_θ] * sqpi_inv
                    inner_env[i_T] += T * weights[j_θ] * sqpi_inv
                    inner_env[i_θl] += θ * weights[j_θ] * sqpi_inv
                    inner_env[i_ρ] += ρ * weights[j_θ] * sqpi_inv
                    # rain area fraction
                    if aux[:qr_src, k, i] > 0
                        inner_env[i_rf] += weights[j_θ] * sqpi_inv
                    end
                    # cloudy/dry categories for buoyancy in TKE
                    if q_liq > 0.0
                        inner_env[i_cf] += weights[j_θ] * sqpi_inv
                        inner_env[i_qt_cld] += qt_hat * weights[j_θ] * sqpi_inv
                        inner_env[i_T_cld] += T * weights[j_θ] * sqpi_inv
                    else
                        inner_env[i_qt_dry] += qt_hat * weights[j_θ] * sqpi_inv
                        inner_env[i_T_dry] += T * weights[j_θ] * sqpi_inv
                    end
                    # products for variance and covariance source terms
                    inner_src[i_Sqt] += -qr_src * weights[j_θ] * sqpi_inv
                    inner_src[i_SH] += θl_rain_src * weights[j_θ] * sqpi_inv
                    inner_src[i_Sqt_H] += -qr_src * θ * weights[j_θ] * sqpi_inv
                    inner_src[i_SH_H] +=
                        θl_rain_src * θ * weights[j_θ] * sqpi_inv
                    inner_src[i_Sqt_qt] +=
                        -qr_src * q_tot * weights[j_θ] * sqpi_inv
                    inner_src[i_SH_qt] +=
                        θl_rain_src * q_tot * weights[j_θ] * sqpi_inv
                end

                outer_env[i_ql] += inner_env[i_ql] * weights[j_qt] * sqpi_inv
                outer_env[i_T] += inner_env[i_T] * weights[j_qt] * sqpi_inv
                outer_env[i_θl] += inner_env[i_θl] * weights[j_qt] * sqpi_inv
                outer_env[i_ρ] += inner_env[i_ρ] * weights[j_qt] * sqpi_inv
                outer_env[i_cf] += inner_env[i_cf] * weights[j_qt] * sqpi_inv
                outer_env[i_qt_cld] +=
                    inner_env[i_qt_cld] * weights[j_qt] * sqpi_inv
                outer_env[i_qt_dry] +=
                    inner_env[i_qt_dry] * weights[j_qt] * sqpi_inv
                outer_env[i_T_cld] +=
                    inner_env[i_T_cld] * weights[j_qt] * sqpi_inv
                outer_env[i_T_dry] +=
                    inner_env[i_T_dry] * weights[j_qt] * sqpi_inv
                outer_env[i_rf] += inner_env[i_rf] * weights[j_qt] * sqpi_inv

                outer_src[i_Sqt] += inner_src[i_Sqt] * weights[j_qt] * sqpi_inv
                outer_src[i_SH] += inner_src[i_SH] * weights[j_qt] * sqpi_inv
                outer_src[i_Sqt_H] +=
                    inner_src[i_Sqt_H] * weights[j_qt] * sqpi_inv
                outer_src[i_SH_H] +=
                    inner_src[i_SH_H] * weights[j_qt] * sqpi_inv
                outer_src[i_Sqt_qt] +=
                    inner_src[i_Sqt_qt] * weights[j_qt] * sqpi_inv
                outer_src[i_SH_qt] +=
                    inner_src[i_SH_qt] * weights[j_qt] * sqpi_inv
            end

            aux[:CF, k, i] = outer_env[i_cf]
            aux[:t_cloudy, k, i] = outer_env[i_T]
            aux[:q_tot_cloudy, k, i] = outer_env[i_qt_cld]
            aux[:q_vap_cloudy, k, i] = outer_env[i_qt_cld] - outer_env[i_ql]
            aux[:q_tot_dry, k, i] = outer_env[i_qt_dry]
            ts = TemperatureSHumEquil(
                param_set,
                aux[:t_cloudy, k],
                p_0 = aux[:p_0, k, i],
                aux[:q_tot_cloudy, k, i],
            )
            aux[:θ_cloudy, k, i] = liquid_ice_pottemp(ts)
            aux[:θ_dry, k, i] = dry_pottemp(ts)

        end
    end
end


text:
## Loading modules
using Distributed
using Plots

addprocs(6)

# Load modules in parallel
@everywhere begin
    include("RP_Cournot/modules/CourErr.jl")
    include("RP_Cournot/modules/MarkovChain.jl")
    include("RP_Cournot/modules/Testing.jl")
    include("modules/Parallel_mod.jl")
    using Random
    ## Setting seed
    Random.seed!(66636)
end

## Settin up directories

dir = @__DIR__ ;
datadir=dir*"/data/oildata.csv" ;
resultsdir=dir*"/results/"

## CourErr

# p, q = load_data(datadir)
pt, qt = MyModules.load_data_select(datadir)

## Markov chain

# Testing
p, q = pt[1:6], qt[1:6,:]
# eδ, eq, status = warm_start(p,q)
# w, δsim, qsim = jump(p,q,eδ,eq)

MyModules.jump_algo()
mchain, status = MyModules.gchain(p,q)
malg, salg =  MyModules.gchain_algo(p,q)

histogram(1:500000,malg[1,:])
histogram(1:500000,mchain[1,:])

mchain[:,1].==mchain[:,rand(1:500000)]
malg[:,1].==malg[:,rand(1:500000)]
rand(1:500000)
## Testing

objfun(collect(1:6), m=malg[:,1:1000])
W = objfun_bbo(collect(1:6), m=malg[:,1:1000])
@time f, d = my_dist(mchain[:,1:1000], collect(1:6))
@time my_dist(mchain[:,rand(1:500000,1000)], ones(6))
@time my_dist(malg[:,rand(1:500000,1000)], ones(6))
# f = ecdf(d)
# cdf(f,W)
f(W*Param.nfast)
1-f(W*Param.nfast)
# # TS_b=T*obj_val
quantile(d,[0.05, 0.5, 0.95])

quantile(d,0.05)
quantile(d,0.95)
pval = f(W)

# Testing the functions

# p, q = pt[1:12], qt[1:12,:]
# mchain, status = gchain(p,q)

@time TS, pval = RP_Cournot_test(mchain)

## Testing

# mchain_op = mchain[:,:,1:20000]
#
# objfun(collect(1:6))
#
# @time objfun_par(collect(1:6))
# @time objfun_s(collect(1:6))
#
# @time objfun_s_bbo(collect(1:6))
# @time objfun_bbo_par(collect(1:6))
# #
# @time TS, pval = RP_Cournot_test_sem(mchain)
# @time TS, pval = RP_Cournot_test_par(mchain)

## Rejection rate
acc_rate_table=SharedArray{Float64,2}(zeros(T*2))
@distributed for i=1:72
    my_dist(mchain[i,:,:], ones(6))
end

##

MyModules.acc_rate(mchain,ones(6),f=1)

##
MyModules.sim_err(malg,ones(6)*100)

@time RP_Cournot_test(malg)

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
using Documenter
using Flux
using NNlib # σ
using Zygote

include("contrib/html_writer.jl")

makedocs(
    sitename = "🦉",
    authors = "초보똥",
    clean = false,
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules=[Flux, NNlib, Zygote],
    pages = Any[
        "Home" => "index.md",
        "Flux ✅" => [
            "Flux 홈" => "Flux/index.md",
            "모델 만들기" =>
              ["기본적인 것" => "Flux/models/basics.md",
               "순환(Recurrence)" => "Flux/models/recurrence.md",
               "정규화(Regularisation)" => "Flux/models/regularisation.md",
               "모델 참조(Model Reference)" => "Flux/models/layers.md"],
            "모델 훈련시키기" =>
              ["최적화" => "Flux/training/optimisers.md",
               "훈련시키기" => "Flux/training/training.md"],
            "원-핫 인코딩" => "Flux/data/onehot.md",
            "GPU 지원" => "Flux/gpu.md",
            "저장 & 불러오기" => "Flux/saving.md",
            "커뮤니티" => "Flux/community.md"
        ],
        "DataFlow ✅" => [
            "DataFlow 버티스(vertices)" => "DataFlow/vertices.md",
        ],
        "Zygote ⏳" => [
            "Home" => "Zygote/index.md",
            "Custom Adjoints" => "Zygote/adjoints.md",
            "Utilities" => "Zygote/utils.md",
            "Complex Differentiation" => "Zygote/complex.md",
            "Flux" => "Zygote/flux.md",
            "Profiling" => "Zygote/profiling.md",
            "Internals" => "Zygote/internals.md",
            "Glossary" => "Zygote/glossary.md"
        ],
    ],
    format = Documenter.HTML(assets = ["assets/custom.css"], prettyurls = !("local" in ARGS))
)