text:
module LibSymspg
using Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibSymspg not installed properly, run Pkg.build(\"LibSymspg\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

export find_primitive, refine_cell,
        niggli_reduce!, delaunay_reduce!,
        ir_reciprocal_mesh,
        get_spacegroup, get_symmetry

include("version.jl")
include("symmetry-api.jl")
include("spacegroup-api.jl")
include("cell-reduce-api.jl")
include("latt-reduce-api.jl")
include("ir-mesh-api.jl")
include("utils.jl")

"""
return type: String, Int64 -> (symbol, spg_number)
"""
function get_spacegroup(lattice::Array{Float64, 2},
                         positions::Array{Float64, 2},
                         types::Array{Int64, 1},
                         num_atom::Int64,
                         symprec::Float64=1e-5)
    #
    positions = Array(transpose(positions))
    db = spg_get_dataset(lattice, positions, types, num_atom, symprec)
    return char2Str(db.international_symbol), Base.convert(Int64, db.spacegroup_number)
end

"""
get_symmetry
"""
function get_symmetry(lattice::Array{Float64, 2},
                      positions::Array{Float64, 2},
                      types::Array{Int64, 1},
                      num_atom::Int64,
                      symprec::Float64=1e-5)
    #
    positions = Array(transpose(positions))
    max_size = 48*num_atom
    rots = Array{Cint, 3}(undef, 3, 3, max_size)
    trans = Array{Float64, 2}(undef, 3, max_size)

    nop = spg_get_symmetry!(rots, trans, max_size, lattice, positions, types, num_atom, symprec)
    @assert nop ≠ 0

    op_rots = Vector{Array{Int, 2}}(undef, nop)
    op_trans = Vector{Vector{Float64}}(undef, nop)
    for i in 1:nop
        op_rots[i] = rots[:, :, i]
        op_trans[i] = trans[:, i]
    end

    op_rots, op_trans
end

function get_symmetry(lattice::Array{Float64, 2},
                      positions::Array{Float64, 2},
                      types::Array{Int64, 1},
                      symprec::Float64=1e-5)
    #
    @assert size(positions)[1] == size(types)[1]
    num_atom = size(types)[1]

    return get_symmetry(lattice, positions, types, num_atom, symprec)
end
"""
"""
function find_primitive(lattice::Array{Float64, 2},
               positions::Array{Float64, 2},
               types::Array{Int64, 1},
               symprec::Float64)

    positions = Array(transpose(positions))
    res_lattice, res_positions, res_types, res_N = spg_find_primitive(lattice, positions, types, symprec)

    res_lattice, Array(transpose(res_positions)), res_types, res_N
end

function find_primitive(lattice::Array{Float64, 2},
                positions::Array{Float64, 2},
                types::Array{Int64, 1},
                symprec::Float64,
                angle_tolerance::Float64)

    positions = Array(transpose(positions))
    res_lattice, res_positions, res_types, res_N = spgat_find_primitive(lattice, positions, types, symprec, angle_tolerance)

    res_lattice, Array(transpose(res_positions)), res_types, res_N
end

"""
"""
function refine_cell(lattice::Array{Float64, 2},
            positions::Array{Float64, 2},
            types::Array{Int64, 1},
            symprec::Float64)

    positions = Array(transpose(positions))
    res_lattice, res_positions, res_types, res_N = spg_refine_cell(lattice, positions, types, symprec)

    res_lattice, Array(transpose(res_positions)), res_types, res_N
end

function refine_cell(lattice::Array{Float64, 2},
            positions::Array{Float64, 2},
            types::Array{Int64, 1},
            symprec::Float64,
            angle_tolerance::Float64)

    positions = Array(transpose(positions))
    res_lattice, res_positions, res_types, res_N = spgat_refine_cell(lattice, positions, types, symprec, angle_tolerance)

    res_lattice, Array(transpose(res_positions)), res_types, res_N
end

"""
"""
niggli_reduce! = spg_niggli_reduce!

"""
"""
delaunay_reduce! = spg_delaunay_reduce!

"""
"""
function ir_reciprocal_mesh(mesh::Array{Int64, 1},
                            is_shift::Array{Bool, 1},
                            is_time_reversal::Bool,
                            lattice::Array{Float64, 2},
                            positions::Array{Float64, 2},
                            types::Array{Int64, 1},
                            num_atom::Int64,
                            symprec::Float64=1e-5)
    #
    positions = Array(transpose(positions))
    is_shift = [flag ? 1 : 0 for flag in is_shift]
    is_time_reversal = is_time_reversal ? 1 : 0

    return spg_get_ir_reciprocal_mesh(mesh, is_shift, is_time_reversal,
                                    lattice, positions, types, num_atom, symprec)
end

end #module LibSymspg

ext: jl


text:
function triangles(Ni, Nj, Nh)
    L = LinearIndices((1:Ni, 1:Nj))
    I = vec([
        (t == 1 ? L[i, j] : L[i + 1, j]) + Ni * Nj * (h - 1) for t in 1:2,
        i in 1:(Ni - 1), j in 1:(Nj - 1), h in 1:Nh
    ])
    J = vec([
        (t == 1 ? L[i + 1, j] : L[i + 1, j + 1]) + Ni * Nj * (h - 1) for
        t in 1:2, i in 1:(Ni - 1), j in 1:(Nj - 1), h in 1:Nh
    ])
    K = vec([
        (t == 1 ? L[i, j + 1] : L[i, j + 1]) + Ni * Nj * (h - 1) for
        t in 1:2, i in 1:(Ni - 1), j in 1:(Nj - 1), h in 1:Nh
    ])
    return (I, J, K)
end

function triangulate(space::SpectralElementSpace2D)
    Ni, Nj, _, _, Nh = size(space.local_geometry)
    return triangles(Ni, Nj, Nh)
end

function triangulate(space::ExtrudedFiniteDifferenceSpace)
    Ni, Nj, _, Nv, Nh = size(space.local_geometry)
    @assert Nj == 1 "triangulation only defined for 1D extruded fields"
    return triangles(Ni, Nv, Nh)
end


text:
# using SymPy

# x, y = @vars x y
# u = x*(1-x)*y*(1-y)
# v = x*(1-x)*y*(1-y)^2

# p = x*(1-x)*y*(1-y)

# ux = diff(u,x)
# uy = diff(u,y)
# vx = diff(v,x)
# vy = diff(v,y)
# px = diff(p,x)
# py = diff(p,y)
# f = u*ux + v*uy - (diff(ux,x)+diff(uy,y)) + px
# g = u*vx + v*vy - (diff(vx,x)+diff(vy,y)) + py
# h = diff(u, x) + diff(v, y)
# println(replace(replace(sympy.julia_code(simplify(f)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(g)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(h)), ".^"=>"^"), ".*"=>"*"))

function u_exact(x,y)
    if y == 0.0
        return 1.0
    else
        return 0.0
    end
end

function v_exact(x,y)
    0.0
end

function p_exact(x,y)
    0.0
end

function ffunc(x, y)
    #-x^2*y*(x - 1)^2*(y - 1)^2*(2*y - 1) + x*y^2*(x - 1)*(2*x - 1)*(y - 1)^2 + x*y*(y - 1) - 2*x*(x - 1) + y*(x - 1)*(y - 1) - 2*y*(y - 1)
    0.0
end
function gfunc(x, y)
    #x^2*y*(x - 1)^2*(y - 1)^3*(3*y - 1) - x*y^2*(x - 1)*(2*x - 1)*(y - 1)^3 + 3*x*y*(x - 1) + 5*x*(x - 1)*(y - 1) + 2*y*(y - 1)^2
    0.0
end

function hfunc(x,y)
    # (y - 1)*(-2*x*y*(x - 1) + x*y - x*(x - 1)*(y - 1) + y*(x - 1))
    0.0
end

# nn to estimate space varying viscosity
function nu_exact(x, y)
    1 + 6 * x^2 + x / (1 + 2 * y^2)
    # (1 + x^2 + x / (1 + 2 * y^2)) * 0.01
    # 1 + x^2
end

using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

m = 20
n = 20
h = 1/n
# nu = 0.01

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc, m, n, h), m, n, h)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc, m, n, h), m, n, h)
H = h^2*eval_f_on_fvm_pts(hfunc, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)
nu_fvm = constant(@. nu_exact(x, y))
nu_gauss = reshape(repeat(nu_fvm, 1, 4), (-1,))

Laplace = compute_fem_laplace_matrix1(nu_gauss, m, n, h)

# Laplace = constant(nu*compute_fem_laplace_matrix1(m, n, h))
function compute_residual(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    g5 = -F2
    G = g1 + g2 + g3 + g4 + g5

    H0 = -B * [u;v] + H

    R = [F;G;H0]
    return R
end

function compute_jacobian(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    J0 = [Fu Fv
          Gu Gv]
    J = [J0 -B'
        -B spdiag(zeros(size(B,1)))]
end

NT = 5
# S = zeros(m*n+2(m+1)*(n+1), NT+1)
# S[1:m+1, 1] = ones(m+1,1)

bd = bcnode("all", m, n, h)
# bd = [bd; bd .+ (m+1)*(n+1); ((1:m) .+ 2(m+1)*(n+1))]
bd = [bd; bd .+ (m+1)*(n+1); 2*(m+1)*(n+1)+(n-1)*m+1:2*(m+1)*(n+1)+(n-1)*m+2] # only apply Dirichlet to velocity

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    # residual[bd] .= 0.0
    residual = scatter_update(residual, bd, zeros(length(bd)))

    d = J\residual
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

# for i = 1:NT
#     residual = compute_residual(S[:,i])
#     J = compute_jacobian(S[:,i])

#     J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
#     residual[bd] .= 0.0

#     d = J\residual
#     S[:,i+1] = S[:,i] - d
#     @info i, norm(residual)
# end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, [u0; v0; p0])

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n))

sess = Session(); init(sess)
output = run(sess, S)

S = output
# out_v = output[:, 1:2*(m+1)*(n+1)]
# out_p = output[:, 2*(m+1)*(n+1)+1:end]

matwrite("steady_cavity_data.mat",
    Dict(
        # "V"=>output[end, 1:2*(m+1)*(n+1)]
        # "V"=>output[:, 1:2*(m+1)*(n+1)]
        "V"=>output
    ))

# figure(figsize=(20,10))
# subplot(321)
# visualize_scalar_on_fem_points(output[NT+1, 1:(m+1)*(n+1)], m, n, h)
# subplot(322)
# visualize_scalar_on_fem_points(u0, m, n, h)

# subplot(323)
# visualize_scalar_on_fem_points(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)], m, n, h)
# subplot(324)
# visualize_scalar_on_fem_points(v0, m, n, h)

# subplot(325)
# visualize_scalar_on_fvm_points(output[NT+1, 2*(m+1)*(n+1)+1:end], m, n, h)
# subplot(326)
# visualize_scalar_on_fvm_points(p0, m, n, h)

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

# prev_data = matread("steady_cavity_data.mat")["V"]

# final_u2=prev_data[NT+1, 1:(1+m)*(1+n)]
# final_v2=prev_data[NT+1, (1+m)*(1+n)+1:2*(m+1)*(n+1)]
# final_p2=prev_data[NT+1, 2*(m+1)*(n+1)+1:end]

# u12 = final_u2[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# u22 = final_u2[Int(n/2)+1:m+1:end]

# v12 = final_v2[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# v22 = final_v2[Int(n/2)+1:m+1:end]
####################################################################################

# final_u=output[NT+1, 1:(1+m)*(1+n)]
# final_v=output[NT+1, (1+m)*(1+n)+1:2*(m+1)*(n+1)]
# final_p=output[NT+1, 2*(m+1)*(n+1)+1:end]

# u1 = final_u[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# u2 = final_u[Int(n/2)+1:m+1:end]

# v1 = final_v[Int(n/2)*(m+1)+1: Int(n/2)*(m+1)+m+1]
# v2 = final_v[Int(n/2)+1:m+1:end]
# xx = 0:h:1

# figure();plot(xx, u1);#plot(xx, u12);
# savefig("u_horizontal.png")

# figure();plot(xx, u2);#plot(xx, u22);
# savefig("u_vertical.png")

# figure();plot(xx, v1);#plot(xx, v12);
# savefig("v_horizontal.png")

# figure();plot(xx, v2);#plot(xx, v22);
# savefig("v_vertical.png")

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

# p12 = final_p2[(Int(n/2)-1)*m+1: (Int(n/2)-1)*m+m]
# p22 = final_p2[Int(n/2)*m+1: Int(n/2)*m+m]
# p32 = 0.5 * (p12 .+ p22)

# p42 = final_p2[Int(n/2):m:end]
# p52 = final_p2[Int(n/2)+1:m:end]
# p62 = 0.5 * (p42 .+ p52)

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

# p1 = final_p[(Int(n/2)-1)*m+1: (Int(n/2)-1)*m+m]
# p2 = final_p[Int(n/2)*m+1: Int(n/2)*m+m]
# p3 = 0.5 * (p1 .+ p2)

# p4 = final_p[Int(n/2):m:end]
# p5 = final_p[Int(n/2)+1:m:end]
# p6 = 0.5 * (p4 .+ p5)

# xx = 0.005:h:1

# figure();plot(xx, p3);#plot(xx, p32);
# savefig("p_horizontal.png")

# figure();plot(xx, p6);#plot(xx, p62);
# savefig("p_vertical.png")


text:
module BDDSolver

    using Laplacians
    using StaticArrays
    using Arpack
    using Statistics

    using SparseArrays
    using LinearAlgebra

    using Random

    include("BlockCSC.jl")
    include("BDDTypes.jl")
    include("acBDD.jl")
    include("approxPQ.jl")
    include("generators.jl")
    include("solvers.jl")

    export approxCholBDD
    export BlockCSC, BlockVec, FlatVec
    export randOrthWt
end # module


text:
"""
# module Settings

# Examples

```jldoctest
julia>
```
"""
module Settings

__precompile__()

using CSV
using EquationsOfState
using EquationsOfState.Collections
using EquationsOfState.NonlinearFitting
import JSON
using MLStyle
import YAML

export load_settings, parse_settings, run_settings

function load_settings(file::AbstractString)
    open(file, "r") do io
        ext = lowercase(splitext(file)[2])
        dict = if ext == ".json"
            JSON.parse(io)
        elseif ext in (".yaml", ".yml")
            YAML.load(io)
        else
            error("Invalid extension!")
        end
        return dict
    end
end  # function load_settings

parse_settings(file::AbstractString) = parse_settings(load_settings(file))
function parse_settings(d::AbstractDict)
    dict = deepcopy(d)
    eos = dict["eos"]
    eos isa AbstractDict || error("The parameter 'eos' must be a YAML dict!")
    for (key, value) in eos
        dict["eos"][key] = eval(Symbol(key))(value...)
    end
    return dict
end  # function parse_settings

run_settings(file::AbstractString) = (run_settings ∘ parse_settings ∘ load_settings)(file)
function run_settings(dict::AbstractDict)
    eos = dict["eos"]
    eos isa AbstractDict || error("The parameter 'eos' must be a YAML dict!")
    any(!isa(e, EquationOfState) for e in values(eos)) && error("The equation of state is not an `EquationsOfState`!")
    dataframe = CSV.read(dict["input"])
    xdata, ydata = dataframe[1], dataframe[2]
    relation = @match dict["relation"] begin
        "E" => EnergyRelation
        "P" => PressureRelation
        "B" => BulkModulusRelation
    end
    result = Dict()
    for (key, value) in eos
        fitted = lsqfit(relation, value, collect(xdata), collect(ydata))
        println("The fitted equation of state is: ", fitted)
        result[key] = fitted
    end
    output = dict["output"]
    isfile(output) || touch(output)
    open(output, "r+") do io
        JSON.print(io, result)
    end
end  # function run_settings

end


text:
function evaluate{T1, T2}(polynomial::MultiPoly.MPoly{T1}, x::AbstractVector{T2})
    @assert length(x) == length(polynomial.vars)
    result = zero(promote_type(T1, T2))
    for (powers, coeff) in polynomial.terms
        term = one(T2)
        for i = 1:length(powers)
            if powers[i] != 0
                term *= x[i] ^ powers[i]
            end
        end
        result += term * coeff
    end
    result
end

type PolynomialScalarField{N, T} <: AbstractScalarField{N}
	polynomial::MultiPoly.MPoly{T}
end
PolynomialScalarField{T}(poly::MultiPoly.MPoly{T}) = PolynomialScalarField{length(poly.vars), T}(poly)
(field::PolynomialScalarField)(x::AbstractVector) = evaluate(field.polynomial, x)

type PolynomialVectorField{N, T} <: AbstractVectorField{N}
	partials::SVector{N, MultiPoly.MPoly{T}}
end
(field::PolynomialVectorField)(x::AbstractVector) = map(p -> evaluate(p, x), field.partials)

gradient{N, T}(field::PolynomialScalarField{N, T}) = PolynomialVectorField{N, T}(map(v -> MultiPoly.diff(field.polynomial, v), field.polynomial.vars))

function linear_fit{T}(coordinates::Array{T, 2}, data::Vector{T})
    dimension = size(coordinates, 1)
    num_points = size(coordinates, 2)
    A = hcat(coordinates', ones(num_points))
    v = A \ data
    coeffs = OrderedDict(zeros(dimension) => v[end])
    for j = 1:dimension
        powers = zeros(Int, dimension)
        powers[j] = 1
        coeffs[powers] = v[j]
    end
    return MultiPoly.MPoly{T}(coeffs, [:x, :y, :z][1:dimension])
end

destructure{T,N}(A::Array{T,N}) = reinterpret(eltype(T), A, (size(T)..., size(A)...))

function linear_fit{N, T}(coordinates::Vector{SVector{N, T}}, data::Vector{T})
    linear_fit(destructure(coordinates), data)
end

function linear_fit{N, T}(coordinates::Vector{SVector{N, T}}, data::Vector{SVector{N, T}})
    A = destructure(coordinates)
    MultiPoly.MPoly{T}[linear_fit(A, [d[i] for d in data]) for i in 1:N]
end


text:
module TestVarious

using Copulas
using Base.Test

## undefined testpcc
##------------------

@test_throws Exception Copulas.testpcc(187)

end


text:
using Test
using PETSc, MPI, LinearAlgebra, SparseArrays

m,n = 20,20
x = randn(n)
V = PETSc.VecSeq(x)

@test norm(x) ≈ norm(V) rtol=10eps()

S = sprand(m,n,0.1) + I
M = PETSc.MatSeqAIJ(S)

@test norm(S) ≈ norm(M) rtol=10eps()

w = M*x
@test w ≈ S*x

ksp = PETSc.KSP(M; ksp_rtol=1e-8, pc_type="jacobi", ksp_monitor=true)
#PETSc.settolerances!(ksp; rtol=1e-8)

@test PETSc.gettype(ksp) == "gmres" # default

pc = PETSc.PC(ksp)
@test PETSc.nrefs(pc) == 2
@test PETSc.gettype(pc) == "jacobi"

# create an extra handle, check ref count is incremented
pc_extra = PETSc.PC(ksp)
@test PETSc.nrefs(pc) == 3
# destroy extra handle, check ptr is set to null, ref count is decremented
PETSc.destroy(pc_extra)
@test pc_extra.ptr == C_NULL
@test PETSc.nrefs(pc) == 2

# safe to call destroy on null pointer
PETSc.destroy(pc_extra)

# set new pc, check ref counts are modified by ksp
pc_new = PETSc.PC{Float64}(MPI.COMM_SELF)
@test PETSc.nrefs(pc_new) == 1
PETSc.setpc!(ksp, pc_new)
@test PETSc.nrefs(pc_new) == 2
@test PETSc.nrefs(pc) == 1

PETSc.setpc!(ksp, pc)
@test PETSc.nrefs(pc_new) == 1
@test PETSc.nrefs(pc) == 2

y = ksp \ w
@test S*y ≈ w  rtol=1e-8

w = M'*x
@test w ≈ S'*x

y = ksp' \ w
@test S'*y ≈ w rtol=1e-8

f!(y,x) = y .= 2 .*x

M = PETSc.MatShell{Float64}(f!,10,10)

x = rand(10)

@test M*x ≈ 2x
@test PETSc.KSP(M) \ x ≈ x/2

function F!(fx, x)
  fx[1] = x[1]^2 + x[1]*x[2] - 3
  fx[2] = x[1]*x[2] + x[2]^2 - 6
end

J = zeros(2,2)
PJ = PETSc.MatSeqDense(J)
function updateJ!(x, args...)
    J[1,1] = 2x[1] + x[2]
    J[1,2] = x[1]
    J[2,1] = x[2]
    J[2,2] = x[1] + 2x[2]
end

S = PETSc.SNES{Float64}(MPI.COMM_SELF; ksp_rtol=1e-4, pc_type="none")
PETSc.setfunction!(S, F!, PETSc.VecSeq(zeros(2)))
PETSc.setjacobian!(S, updateJ!, PJ, PJ)
@test PETSc.solve!([2.0,3.0], S) ≈ [1.0,2.0] rtol=1e-4


text:
#--- show julia builds code

"""that can give you insight of types (of your code)

1.text Files read as a string
2.-julia lets Files gets parsed into Julia `Expr`
`Expr`: a structured representation of the code (you wrote)

3.it organizes structures (for you)- you can manipulate,
so that you don't do that (as a form of Raw Strings)

e.g.
Matching the paranthesises

- Expressions get lowered (@code_lowered) to this lowered representation llvm
e.g. you're v. familiar in seeing this message
the `Info: pre-compiling...` message, when you update a package (or install a new one)
Tip: this (llvm) is mosly what gets saved to a file when that happens

juila is parsing text files into Expressions
2. then lowering them
3. then julia saves the lowered code to a file
(on your harddrive)

it can save that work, next time you run julia

the package is asked for additional information (to be stored )

- There is more information into those pre-compile files - only most of what's happening (code also turned to be native )
#--------

when you call `f(args...)` either

-the native code gets looked up in the in-memory torage & then runs (typo @ 22:42)
OR:

1. f gets type-inferred for the argument types
2. th type-inferred code gets optimized by Julia (this gets stored in memory)
3. the result gets compiled by LLVM (this also gets stored in memory)
4. Julia runs the compiled code

when you use a package, it loads that past file - that has lowered code stored in it
whenever you make a function call, either
    if already compiled, julia has it in memory
        it will look it up , & run it
        or
    if never called it before,  then
    Julia will Perform type inference for specific types of arguments that
    you're passing in ,
    that type-inferred code gets optimized (optimize=true stage)
    result gets compiled by an LLVM
    finally, it will run the compile code
#------Dispatch runtime vs compile-time

Hoping you'd appreciate the importance of types to Julia's performance

Remember the addTwo
this leaves a question open

what is this + here?

+ is a function call , generic , 200 methods (more if we loaded additional packages)
(as packages can extend functions)

How does julia know which (of the methods)is the right one?

It chooses the `Most specific method` for the the types of an argument

types of the argument(here ): x[1] , x[2]
could've been written as:
+(x[1] , x[2]) #then, it'd be more appearant this is a function call , by plus
which is the applicable one here-
we can see the decision process using @code_typed (it hates the UndefVars )
give all arguments a value for  addTwo( x[1] , x[2])
    @code_typed: answers which of methods gets called
"""
addTwo(x) = x[1] + x[2]
"""adding 2 numbers"""
#we'll review this code
#@code_typed optimize = false addTwo(1, 2.0) #Any[] #no problem

#again rewriting in another way:
addTwo(x) = +(x[1], x[2])
"""
CodeInfo(
1 ─ %1 = Base.getindex(x, 1)::Int64
│   %2 = Base.getindex(x, 2)::Float64
│   %3 = (%1 + %2)::Float64
└──      return %3
) => Float64
"""
@code_typed optimize = false addTwo((1, 2.0)) #Any[]
#=
since it's a typed code, it can figure out which method is the one
(applicable for those 2 types):
while compiling this code,
    figure out which method is the one applicable
for those 2 types
-the point: it didn't have to run the code (to figure it out )

the Magic: type-infernce : when i run this code, what types will %1 & %2 be ?

If know with certainty, then I can figure out which method will
Be applicable to them  (at the time I'm building the code ) - not at the time I'm running

here's julia main Performance bit of  magic:
takeaway:
instead of running it with optimize=false  , can see having figured out which method is gonna be called, in optimize=true
=#
@code_typed optimize = true addTwo((1, 2.0)) #optimize=true performing inlining: inline_expansion

"""
 this ability to 'pear-in the future'
 & look-up the functions ahead of time
 is again a key component

 this is called
 "Compile-time Dispatch"

 meaning choosing hich method
 will get executed happens when the calling function is being compiled

alternative is called runtime Dispatch
we're gonna defeat the ability of julia's

takeaway:type-inference engine (will) work
by creating a Vector that contains "No Information about the types (of objects) stored in it  "

others like: vector int , Float64, nut this line will gonna create a vector of Any
"""

@code_typed optimize = true addTwo(Any[1, 2.0]) ## call irt on a Vector(Any)

#=
Literally any kind of Object, String , a Dictionary,an image, anything you want
Might be stored in the slots of this Vector
takeaway: Julia Cannot Predict the types of the Object that will be gen_call_with_extracted_types_and_kwargs  under any circumstances,

therefore eventhough we turned on optimize = true , you'd hope julia will do all types of inlining
-heere thre's nothing to inline - it has no clue what method
+ will be applicable in this circumstance

the first object of type Any - julia has no Idea what type it is
2nd tye: inferred to be a type Any , therefore + is also any
why: if it doesn't know the inputs gonna be , it can't figure out what the
output is gonna be

julia knows very very little here
 & it's not made any sort of progress on optimizing your code
 for  this particular caes

this ould have then to be executed by a runtime dispach
=#

"""
1.It knows this arrayref is gonna be called (that's progress of perform)
2.It Doesn't know what the result will be
that reason it knows

vector of any isn't the type
vector of any isn't the type (that julia understands)
that part is a concrete type can figureout
what's gonna be called here

but it doesn't know what the result will be

takeaway: the key decision has to be made, once called the + operator

(it depends on the arguments)

Only piece of information, julia has to go on  - It doesn't get very far

means: this method will be looked up & called when the function is running
-this is called run-time dispach:
schemeatuc of Runtime vs compile-time Dispatch

1. preparing the Arguments
2. Deciding which specifc  compuled blob to use.like looking up someone's
phonenumber in the phone book

julia literally scans through the method tables - living in a particular memory location
machine:
1.calling funcition invloves:
1.Arguments: getting them ready to be passed in
& then deciding which specific compile blob to call

like clicking-up somebody's phonenumber (in a phonebook) [me:Database Analogy]
-which page, column
-it may take time flipping pages** then look-up the phonenumber

* that is already been debated: for a phonebook, it's apriori ordered by the first letter
this includes a huge proning of the Search Space

#--- Schematic of a compiletime call in psuedo-julia

push!(execution_stack, args)
@goto compiled_blob_52383

the blob will retrieve the argument values by popping the execution execution_stack

-get the arguments ready to go
-transfer execution (to the blob of compile code) [super-fast]

schematic of a runtime in psuedo-Julia:
1.first: figure out which blob you're gonna get -for the types of arguments
(this happens during the running of a function)
2.look-up process: paging through the method-tables
trying to figure-out which is gonna be applicable to this
particular set of arguments is the slow process

# scan the method tables and their lists of compiked blobs for a match
# if the right blob hasn't been compiled yet, compile it now

"""
blob = get_blob_forargtypes(f, typesof(args)) #Error types of not defined

# the rest looks the same as a compiletime call:
push!(execution_stack, args)
goto(blob)

#---comparing performance of runtime vs compile-time dispatch
#=built-in @time macro
    handy for something, but it's a very fast function=#

using BenchmarkTools

@btime addTwo(z) setup = (z = rand(2))
""

"""
 0.001 ns (0 allocations: 0 bytes)
3.0
it's solved with a nano-second with the cpu codeblock (me:minimum)
"""

@btime addTwo((1, 2.0)) setup = (x = rand(1:5); y = rand())

"
we wanna generate  a random integer between 1 & 5 for x

 because these are not predictable*

 2.300 ns (0 allocations: 0 bytes) & a random floating number for y #minimum

2.0372456559896017
----------*footnote:
predictability depends on: 1. rngType (subsequently seed choice #TODO:check persistence of a seed) 2.rngSampler: the 'intended' expected shape

"""
@btime addTwo((x, y)) setup = (x = rand(1:5); y = rand())

"""
  2.300 ns (0 allocations: 0 bytes)
2.002507990680424
pretty much sure that is has been solved
"""
@btime addTwo((x, y)) setup = (x = rand(1:5); y = rand(); z = Any[x, y])

"""
this comes from fact julia has to do method look-up  , at run-time -rather than compile time

takeaway: make these benchmarks a handy-dandy tool for
Ballpark costs of runtime dispach (me: gives a rough sketch of performance)

1. single argument: 15-35ns [hardware specific ]
2. two arguments: ~100ns

cost of a single run-time dispatch

whether it's trivial or catastrophic, it totally depends on what you're Doing

1.Trivial:
If this is 1 operation , it will happen once
Over the course of your computation
completely utterly irrelevant

2. complicated
if it's in an inner-loop (of some computation)
that you're running a billion times
[billion x 1 op 15-35ns]

Many operations start out taking a large fraction a minute
After a few optimizations, it's a tiny fraction of a second

Eliminating run-time dispatch (for newbs) is a single most common thing they could do

Runtime dispatch isslow

#--- Union-splitting (new)

an intermediate case
-julia Can determine that there are only a few possible arguemnt types

instead of branching your code (checking each type & element it is )

here we know with certainty that
they are (a choice) [from a UnionSplitting] either:
1.Float32
2.Float64

what we
  4.600 ns (0 allocations: 0 bytes)
0.7754799894275148

Julia checks what type of arguemnt is
Then check 1 type of blob

If it's the other type then it will use the other compile blob

Warning: there are no blobs looked up at run-time

You just don't know which ones are applicable

It can check that at runtime
But especially if these types are  concrete types
There is only one such type  that fits that description

"""

@btime addTwo(z) setup = (x = rand(Float32); y = rand(); z = Union{Float32,Float64}[x, y])

argtypes = typeof(args)
push!(execution_stack, args)
if argtypes === Tuple{Float32}
    @goto compiled_blob_52383

else #he only other option is Tuple{Float64}
    @goto compiled_blob_52951 #demo #donotRun

#=
Note: the absence of the need to call `get_blob_got_argtypes`
Uion-splitting generalizes compiletime dispatch

now debug code make it working, let's start profiling

=#

max_stars_repo_path: Tutorial/Tim Holy/buildcode.jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_path: Tutorial/Tim Holy/buildcode.jl

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_path: Tutorial/Tim Holy/buildcode.jl

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
################################################################################
#
#   DEFINITION OF HAMILTONIAN
#   - label type L (of bond labels)
#   - unitcell type U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
#   - bond hamiltonian type HB<:AbstractBondHamiltonian{L,NS} where {NS}
#
#
#    NOTE: Hamiltonian defined here is in the level of Majoranas.
#
#
#    NOTE ABOUT GAUGE: Refer to sign_structure.jl
################################################################################

struct Hamiltonian{
        L,
        U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
        HB<:AbstractBondHamiltonian{L,NS} where {NS}
    }

    # the unitcell
    unitcell :: U

    # the bond Hamiltonian
    h_bond   :: HB

    # ...

end

# function to create such an object
function newHamiltonian(
            unitcell :: U,
            bond_hamiltonian :: HB
        ) :: Hamiltonian{L,U,HB} where {
             L,U<:AbstractUnitcell{S,B} where {
             S, B<:AbstractBond{L,N} where {N}},HB<:AbstractBondHamiltonian{L,NS} where {NS}}

    # return such an object
    return Hamiltonian{L,U,HB}(unitcell, bond_hamiltonian)
end

# delta function for unitcell and bond to compute the displacement vector of the bond
function delta(
        uc :: U,
        b  :: B
    ) :: Vector{Float64} where {L,N,S,B<:AbstractBond{L,N},U<:AbstractUnitcell{S,B}}

    # build up the difference vector within the unitcell
    d = point(site(uc,to(b))) .- point(site(uc,from(b)))
    # add all the bravais lattice vectors
    for i in 1:N
        d .+= wrap(b)[i] .* latticeVectors(uc)[i]
    end
    # return the vector
    return d
end

# function to obtain the Hamiltonian matrix at some k (nearest+next-nearest neighbor contribution)
function getMatrixAtK(
        hamiltonian :: Hamiltonian{L,U,HB},
        k :: Vector{<:Real},
        kappa :: Float64
    ) :: Matrix{Complex} where {L,NS,U,HB<:AbstractBondHamiltonian{L,NS}}

    # new matrix

    h = zeros(Complex, NS*numSites(hamiltonian.unitcell),NS*numSites(hamiltonian.unitcell))
    adj_delta_r = 0+0im
    all_bonds = bonds(hamiltonian.unitcell)
    # add all bonds to the matrix
    for b in bonds(hamiltonian.unitcell)
        # get the bond indices
        i_from  = from(b)
        i_to    = to(b)
        delta_r = delta(hamiltonian.unitcell, b)

        # get the interaction matrix and add it to the general matrix

    h[(i_from-1)*(NS)+1:(i_from)*(NS), (i_to-1)*(NS)+1:(i_to)*(NS)] .+= 0.5 .*im .* plaquette(hamiltonian.unitcell, b) .* bondtermhopping(hamiltonian.h_bond, b).* exp(im*dot(k,delta_r))

  end
    h_NN = getMatrixAtKNN(hamiltonian,k,kappa) #hamiltonian for next-nearest neighbor contribution

    # return the matrix
    return (h .+ h_NN)
end

# function determine the energy of a k_vector for a given hamiltonian
function energy(
            hamiltonian :: Hamiltonian{L,U,HB},
            k           :: Vector{Float64},
           kappa :: Float64
        ) :: Float64 where {L,NS,U,HB<:AbstractBondHamiltonian{L,NS}}

    # first compute the matrix at the given k
    M_k = getMatrixAtK(hamiltonian, k,kappa)

    # return the minimal eigenvalue
    return eigmin(M_k)
end

# function to obtain the Hamiltonian matrix at some k (only next-nearest neighbor contribution)
function getMatrixAtKNN(
        hamiltonian :: Hamiltonian{L,U,HB},
        k :: Vector{<:Real},
        kappa :: Float64
    ) :: Matrix{Complex} where {L,NS,U,HB<:AbstractBondHamiltonian{L,NS}}

    # new matrix

    h = zeros(Complex, NS*numSites(hamiltonian.unitcell),NS*numSites(hamiltonian.unitcell))

    all_bonds = getBondsNNNByBondDistance(hamiltonian.unitcell)
    # add all bonds to the matrix
    for b in all_bonds
        # get the bond indices
        i_from  = from(b)
        i_to    = to(b)
        delta_r = delta(hamiltonian.unitcell, b)

        # get the inteaction matrix and add it to the general matrix
          h[(i_from-1)*(NS)+1:(i_from)*(NS), (i_to-1)*(NS)+1:(i_to)*(NS)] .+= 0.5.* im.* plaquetteNNN(hamiltonian.unitcell, b) .* kappa .* exp(im*dot(k,delta_r))

           end
    # return the matrix
    return h
end


text:
#=
  本文件配置个性化参数，具体参见 initAlarmSegparas! 中的注释
=#

function initAlarmSegparas1!(SegAlarmParasDictObj::Dict,SrcDictObj::Dict )   #####################
    SegParaNameTuple1 =  ( segID = 0 ,
          microtime=Seg_MicroTimeParaStruct() ,
          microSpace=Seg_MicroSpaceParaStruct(),
          macroSpace= Seg_MacroSpaceParaStruct(),
          macroTime= Seg_MacroTimeParaStruct(),
          )

    SegParaNameTuple1.microtime.StaticThreshold = Int(2200)
    SegParaNameTuple1.microtime.DynamicThreshold_CalTimePeriod = UInt( 60 )
    SegParaNameTuple1.microtime.DynamicThresHold_calCFactor=Int( 130 )
    SegParaNameTuple1.microtime.minTh=2200
    SegParaNameTuple1.microtime.maxTh=2200
    #微观过滤和合并 时间
    SegParaNameTuple1.microtime.microMergeTimePeriod  = Int(2)
    SegParaNameTuple1.microtime.microFilterTimePeriod = Int(2)

    #微观过滤和合并 空间
    SegParaNameTuple1.microSpace.microMergeSpaceLength = UInt(1)
    SegParaNameTuple1.microSpace.microFilterSpaceLength = UInt(3)

    #宏观 空间 限制范围
    SegParaNameTuple1.macroSpace.macroMaxSpaceLength = UInt(100)
    SegParaNameTuple1.macroSpace.macroMinSpaceLength = UInt(4)
    #宏观 时间  限制范围
    SegParaNameTuple1.macroTime.macroMaxTimeInterval = UInt(30)
    SegParaNameTuple1.macroTime.macroMinTimeInterval = UInt(15)

    SegAlarmParasDictObj[SegParaNameTuple1.segID ] = SegParaNameTuple1
end

function initAlarmSegparas!(SegAlarmParasDictObj::Array{NamedTuple,1},SrcDictObj::Dict )   #####################
  a= keys(SrcDictObj)
  b= a
  # v=[]
  # for a1 in  b
  #   append!(v,a1)
  # end

  # for v1 in sort(v)
  #   println( v1)
  # end
  # exit()
  for segid in b
      if segid > 8192
          println( segid ,"---SEGID OUTof Range--$segid----")
          exit()
      end
      if segid < 0
        println( segid ,"-----SEGID OUTof Range--$segid---")
        exit()
      end
      # println( segid ,"---------")
      aobj = ( segID = segid , microtime=Seg_MicroTimeParaStruct() ,
                            microSpace=Seg_MicroSpaceParaStruct(),
                            macroSpace= Seg_MacroSpaceParaStruct(),
                            macroTime= Seg_MacroTimeParaStruct(),
                          )

      valVector = SrcDictObj[ segid ]
      #="通道,标牌坐标1,标牌坐标2,
      4 ~7 :动态阈值时长,动态阈值系数,动态阈值下限,动态阈值上限,
      8 ~11 : 事件合并时长,事件滤除时间,事件合并距离,事件滤除距离,
      12~15: 告警空间下限,告警空间上限,告警续时下限,告警续时上限"
      =#
     # println("exit out 1 !!!!"); exit();
     l1= length( valVector)
    #  for i in 1:l1
    #     println( valVector[i] )
    #  end
      aobj.microtime.StaticThreshold = typemax(Int)
      aobj.microtime.DynamicThreshold_CalTimePeriod = UInt( valVector[4-3] )

      aobj.microtime.DynamicThresHold_calCFactor = Int(  valVector[5-3]   )
      # println(aobj.microtime.DynamicThresHold_calCFactor );
      # println("---------------")
      aobj.microtime.minTh = Int(  valVector[6-3]   )
      aobj.microtime.maxTh = Int(  valVector[7-3]   )
      #微观过滤和合并 时间
      aobj.microtime.microMergeTimePeriod  = Int(  valVector[8-3]   )
      aobj.microtime.microFilterTimePeriod = Int(  valVector[9-3]   )

      #微观过滤和合并 空间
      aobj.microSpace.microMergeSpaceLength =  UInt( valVector[10-3]   )
      aobj.microSpace.microFilterSpaceLength = UInt( valVector[11-3]  )

      #宏观 空间 限制范围
      aobj.macroSpace.macroMaxSpaceLength = UInt(  valVector[13-3]   )
      aobj.macroSpace.macroMinSpaceLength = UInt(  valVector[12-3]   )
      #宏观 时间  限制范围
      #println("exit out 3 !!!!"); exit();
      aobj.macroTime.macroMaxTimeInterval = UInt(  valVector[15-3]   )
      aobj.macroTime.macroMinTimeInterval = UInt(  valVector[14-3]   )
    #  println("exit out 6 !!!!",typeof( segid)  , typeof(SegAlarmParasDictObj)  ); exit();
    # println( "-typeof(segid)--",segid, "=$typeof(aobj)== ",aobj)
      SegAlarmParasDictObj[ segid ] = aobj
    #  SegAlarmParasDictObj[aobj.segID ] = aobj
    # println("exit out 7 !!!!",typeof( segid)); # exit();
      # break
  end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using Documenter, OhMyREPL

makedocs(
    format = :html,
    sitename = "OhMyREPL",
    pages = Any[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Features" => Any[
            "features/syntax_highlighting.md",
            "features/bracket_highlighting.md",
            "features/bracket_complete.md",
            "features/prompt.md",
            "features/rainbow_brackets.md",
            "features/markdown_highlight.md",
            ],
        "Internals" => Any[
            "internals/passes.md"
            ]
    ]
)

deploydocs(
    repo = "github.com/KristofferC/OhMyREPL.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
"""
alpha,beta,vk,vkm1 = LanczosStep!(A::Function,vk,vkm1,beta;sigma=0.0,tol=1e-10)

Performs a Lanczos step

Implementation is based on Table 2.1 in

Choi, S.-C. T. (2006).
Iterative Methods for Singular Linear Equations and Least-squares Problems.
Phd thesis, Stanford University.
"""
function LanczosStep!(A::Function,vk,vkm1,beta;sigma=0.0,tol=1e-10)
    n     = length(vk)

    pk    = A(vk)
    if sigma!=0.0
        BLAS.axpy!(n,-sigma,vk,1,pk,1) #pk   -= sigma*vk
    end

    alpha = BLAS.dot(n,vk,1,pk,1)  # alpha = dot(vk,pk)
    BLAS.axpy!(n,-alpha,vk,1,pk,1) # pk   -= alpha*vk
    BLAS.axpy!(n,-beta,vkm1,1,pk,1)# pk   -= beta*vkm1
    vkm1  = BLAS.blascopy!(n,vk,1,vkm1,1) # vkm1 = copy(vk)
    vk    = pk
    beta  = BLAS.nrm2(n,vk,1)
    if beta>tol
         vk /= beta
    end
    return alpha, beta,vk,vkm1
end

"""
c,s,r = SymOrtho(a,b)

Computes a Givens rotation

Implementation is based on Table 2.9 in

Choi, S.-C. T. (2006).
Iterative Methods for Singular Linear Equations and Least-squares Problems.
Phd thesis, Stanford University.
"""
function symOrtho(a,b)
    c = 0.0; s = 0.0; r = 0.0
	if b==0
		s = 0.0
		r = abs(a)
		c = (a==0) ? c=1.0 : c = sign(a)
	elseif a == 0
		c = 0.0
		s = sign(b)
		r = abs(b)
	elseif abs(b) > abs(a)
		tau = a/b
		s   = sign(b)/sqrt(1+tau^2)
		c   = s*tau
		r   = b/s
	elseif abs(a) > abs(b)
		tau = b/a
		c   = sign(a)/sqrt(1+tau^2)
		s   = c*tau
		r   = a/c
	end
	return c,s,r
end


text:
function save_new_ply(filename::String, PC::PointCloud; comment="julia"::String, ascii=true::Bool)
   # per ora funziona solo in formato ascii BUG nella libreria
   ply = Ply()
   push!(ply, PlyComment(comment))

   points = PC.coordinates
   # position and colors
   if !isempty(PC.rgbs)
      rgbs = PC.rgbs
      vertex = PlyElement("vertex",
                          ArrayProperty("x", points[1,:]),
                          ArrayProperty("y", points[2,:]),
                          ArrayProperty("z", points[3,:]),
                          ArrayProperty("r", rgbs[1,:]),
                          ArrayProperty("g", rgbs[2,:]),
                          ArrayProperty("b", rgbs[3,:]))
   else
      vertex = PlyElement("vertex",
                          ArrayProperty("x", points[1,:]),
                          ArrayProperty("y", points[2,:]),
                          ArrayProperty("z", points[3,:]))
   end

   push!(ply, vertex)

   # For the sake of the example, ascii format is used, the default binary mode is faster.
   save_ply(ply, filename, ascii=true)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Epigraph of real or complex matrix spectral norm (i.e. maximum singular value)
for a matrix (stacked column-wise) of `nrows` rows and `ncols` columns with
`nrows ≤ ncols`.

    $(FUNCTIONNAME){T, R}(nrows::Int, ncols::Int, use_dual::Bool = false)
"""
mutable struct EpiNormSpectral{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d1::Int
    d2::Int
    is_complex::Bool

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_aux_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    W::Matrix{R}
    Z::Matrix{R}
    fact_Z
    Zi::Matrix{R}
    tau::Matrix{R}
    HuW::Matrix{R}
    Huu::T
    trZi2::T
    WtauI::Matrix{R}
    Zitau::Matrix{R}
    tempd1d2::Matrix{R}
    tempd1d2b::Matrix{R}
    tempd1d2c::Matrix{R}
    tempd1d2d::Matrix{R}
    tempd1d1::Matrix{R}
    tempd2d2::Matrix{R}
    tempd2d2b::Matrix{R}

    function EpiNormSpectral{T, R}(
        d1::Int,
        d2::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert 1 <= d1 <= d2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.is_complex = (R <: Complex)
        cone.dim = 1 + vec_length(R, d1 * d2)
        cone.d1 = d1
        cone.d2 = d2
        return cone
    end
end

reset_data(cone::EpiNormSpectral) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated =
    cone.hess_fact_updated = false)

# TODO only allocate the fields we use
function setup_extra_data!(
    cone::EpiNormSpectral{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    dim = cone.dim
    (d1, d2) = (cone.d1, cone.d2)
    cone.W = zeros(R, d1, d2)
    cone.Z = zeros(R, d1, d1)
    cone.Zi = zeros(R, d1, d1)
    cone.tau = zeros(R, d1, d2)
    cone.HuW = zeros(R, d1, d2)
    cone.WtauI = zeros(R, d2, d2)
    cone.Zitau = zeros(R, d1, d2)
    cone.tempd1d2 = zeros(R, d1, d2)
    cone.tempd1d2b = zeros(R, d1, d2)
    cone.tempd1d2c = zeros(R, d1, d2)
    cone.tempd1d2d = zeros(R, d1, d2)
    cone.tempd1d1 = zeros(R, d1, d1)
    cone.tempd2d2 = zeros(R, d2, d2)
    cone.tempd2d2b = zeros(R, d2, d2)
    return cone
end

get_nu(cone::EpiNormSpectral) = cone.d1 + 1

function set_initial_point!(
    arr::AbstractVector,
    cone::EpiNormSpectral{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

function update_feas(cone::EpiNormSpectral{T}) where T
    @assert !cone.feas_updated
    u = cone.point[1]

    if u > eps(T)
        @views vec_copyto!(cone.W, cone.point[2:end])
        copyto!(cone.Z, abs2(u) * I)
        mul!(cone.Z, cone.W, cone.W', -1, true)
        cone.fact_Z = cholesky!(Hermitian(cone.Z, :U), check = false)
        cone.is_feas = isposdef(cone.fact_Z)
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiNormSpectral{T}) where {T <: Real}
    u = cone.dual_point[1]
    if u > eps(T)
        W = @views vec_copyto!(cone.tempd1d2, cone.dual_point[2:end])
        return (u - sum(svdvals!(W)) > eps(T))
    end
    return false
end

function update_grad(cone::EpiNormSpectral)
    @assert cone.is_feas
    u = cone.point[1]
    Zi = cone.Zi

    ldiv!(cone.tau, cone.fact_Z, cone.W)
    inv_fact!(Zi, cone.fact_Z)
    copytri!(Zi, 'U', true)

    cone.grad[1] = -u * tr(Hermitian(Zi, :U))
    @views vec_copyto!(cone.grad[2:end], cone.tau)
    cone.grad .*= 2
    cone.grad[1] += (cone.d1 - 1) / u

    cone.grad_updated = true
    return cone.grad
end

function update_hess_aux(cone::EpiNormSpectral)
    @assert cone.grad_updated
    u = cone.point[1]
    tau = cone.tau
    Zitau = cone.Zitau
    WtauI = cone.WtauI

    copyto!(Zitau, tau)
    ldiv!(cone.fact_Z, Zitau)
    @. cone.HuW = -4 * u * Zitau
    cone.trZi2 = sum(abs2, cone.Zi)
    cone.Huu = 4 * abs2(u) * cone.trZi2 + (cone.grad[1] - 2 *
        (cone.d1 - 1) / u) / u
    copyto!(WtauI, I)
    mul!(WtauI, cone.W', tau, true, true)

    cone.hess_aux_updated = true
    return
end

function update_hess(cone::EpiNormSpectral)
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    d1 = cone.d1
    d2 = cone.d2
    Zi = cone.Zi
    tau = cone.tau
    WtauI = cone.WtauI

    # H_W_W part
    # TODO parallelize loops
    idx_incr = (cone.is_complex ? 2 : 1)
    r_idx = 2
    for i in 1:d2, j in 1:d1
        c_idx = r_idx
        @inbounds for k in i:d2
            taujk = tau[j, k]
            WtauIik = WtauI[i, k]
            lstart = (i == k ? j : 1)
            @inbounds for l in lstart:d1
                term1 = Zi[l, j] * WtauIik
                term2 = tau[l, i] * taujk
                spectral_kron_element!(H, r_idx, c_idx, term1, term2)
                c_idx += idx_incr
            end
        end
        r_idx += idx_incr
    end
    H .*= 2

    # H_u_W and H_u_u parts
    @views vec_copyto!(H[1, 2:end], cone.HuW)
    H[1, 1] = cone.Huu

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormSpectral,
    )
    cone.hess_aux_updated || update_hess_aux(cone)
    u = cone.point[1]
    W = cone.W
    tempd1d2 = cone.tempd1d2
    tempd1d1 = cone.tempd1d1

    @inbounds for j in 1:size(prod, 2)
        arr_1j = arr[1, j]
        @views vec_copyto!(tempd1d2, arr[2:end, j])
        prod[1, j] = cone.Huu * arr_1j + real(dot(cone.HuW, tempd1d2))
        mul!(tempd1d1, tempd1d2, W')
        @inbounds for k in 1:cone.d1
            @inbounds for i in 1:k
                tempd1d1[i, k] += tempd1d1[k, i]'
            end
            tempd1d1[k, k] -= 2 * u * arr_1j
        end
        mul!(tempd1d2, Hermitian(tempd1d1, :U), cone.tau, 2, 2)
        ldiv!(cone.fact_Z, tempd1d2)
        @views vec_copyto!(prod[2:end, j], tempd1d2)
    end

    return prod
end

function dder3(cone::EpiNormSpectral, dir::AbstractVector)
    @assert cone.hess_aux_updated

    u = cone.point[1]
    W = cone.W
    u_dir = dir[1]
    @views W_dir = vec_copyto!(cone.tempd1d2, dir[2:end])
    dder3 = cone.dder3

    Zi = Hermitian(cone.Zi, :U)
    tau = cone.tau
    Zitau = cone.Zitau
    WtauI = cone.WtauI
    tempd1d2b = cone.tempd1d2b
    tempd1d2c = cone.tempd1d2c
    tempd1d2d = cone.tempd1d2d
    tempd1d1 = cone.tempd1d1
    tempd2d2 = cone.tempd2d2
    tempd2d2b = cone.tempd2d2b

    mul!(tempd2d2b, W_dir', tau)
    ldiv!(tempd1d2d, cone.fact_Z, W_dir)
    mul!(tempd1d2b, tempd1d2d, WtauI)
    mul!(tempd1d2c, tempd1d2d, tempd2d2b')
    mul!(tempd1d1, tempd1d2d, W')
    mul!(tempd2d2, tempd2d2b, tempd2d2b)

    mul!(tempd2d2, W_dir', tempd1d2b, true, true)
    mul!(tempd1d2d, tau, tempd2d2)
    mul!(tempd1d2d, tempd1d2c, WtauI, true, true)
    mul!(tempd1d2d, tempd1d2b, tempd2d2b, true, true)

    ldiv!(cone.fact_Z, tempd1d2b)
    mul!(tempd1d2b, Zitau, tempd2d2b, true, true)

    mul!(tempd1d1, tau, W_dir', true, true)
    mul!(tempd1d2b, tempd1d1, Zitau, true, true)
    tempd1d2b .*= -2 * u

    const1 = 4 * u * u_dir * u
    @. tempd1d2c = const1 * Zitau - u_dir * tau
    ldiv!(cone.fact_Z, tempd1d2c)
    tempd1d2b .+= tempd1d2c

    axpby!(-2 * u_dir, tempd1d2b, -2, tempd1d2d)
    @views vec_copyto!(dder3[2:end], tempd1d2d)

    trZi3 = sum(abs2, ldiv!(tempd1d1, cone.fact_Z.L, Zi))
    @. tempd1d2b += 3 * tempd1d2c
    dder3[1] = -real(dot(W_dir, tempd1d2b)) - u * u_dir * (6 * cone.trZi2 -
        8 * u * trZi3 * u) * u_dir - (cone.d1 - 1) * abs2(u_dir / u) / u

    return dder3
end


text:
# this is the interface
macro sort(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Expr, Nothing})
    error("""\n
        The syntax is:
            sort <varlist>
    """)
end

# this is `sort <varlist>`
macro sort(t::Symbol, by::Nothing, sort::Nothing, arguments::Vector{Symbol}, filter::Nothing, use::Nothing, options::Nothing)
    return esc(
        quote
            Douglass.@sort!($t, $arguments)
        end
    )
end

macro sort!(t::Symbol, varlist::Vector{Symbol})
    esc(
        quote
            # check that all variables are present
            Douglass.@assert_vars_present($t, $varlist)
            sort!($t, $varlist)
        end
    )
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule XyceWrapper_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("XyceWrapper")
JLLWrappers.@generate_main_file("XyceWrapper", UUID("11346573-4739-5234-a9de-0f0819ae11a7"))
end  # module XyceWrapper_jll


text:
# Since Kronecker does not implement `TensorCore.:⊗` but instead exports its own function
# `Kronecker.:⊗`, only the module is imported and Kronecker.:⊗ and Kronecker.kronecker are
# called explicitly.
using .Kronecker: Kronecker

export kernelkronmat
export kronecker_kernelmatrix

@doc raw"""
    kernelkronmat(κ::Kernel, X::AbstractVector{<:Real}, dims::Int) -> KroneckerPower

Return a `KroneckerPower` matrix on the `D`-dimensional input grid constructed by ``\otimes_{i=1}^D X``,
where `D` is given by `dims`.

!!! warning

    Require `Kronecker.jl` and for `iskroncompatible(κ)` to return `true`.
"""
function kernelkronmat(κ::Kernel, X::AbstractVector{<:Real}, dims::Int)
    checkkroncompatible(κ)
    K = kernelmatrix(κ, X)
    return Kronecker.kronecker(K, dims)
end

@doc raw"""

    kernelkronmat(κ::Kernel, X::AbstractVector{<:AbstractVector}) -> KroneckerProduct

Returns a `KroneckerProduct` matrix on the grid built with the collection of vectors ``\{X_i\}_{i=1}^D``: ``\otimes_{i=1}^D X_i``.

!!! warning

    Requires `Kronecker.jl` and for `iskroncompatible(κ)` to return `true`.
"""
function kernelkronmat(κ::Kernel, X::AbstractVector{<:AbstractVector})
    checkkroncompatible(κ)
    Ks = kernelmatrix.(κ, X)
    return reduce(Kronecker.:⊗, Ks)
end

@doc raw"""
    iskroncompatible(k::Kernel)

Determine whether kernel `k` is compatible with Kronecker constructions such as [`kernelkronmat`](@ref)

The function returns `false` by default. If `k` is compatible it must satisfy for all ``x, x' \in \mathbb{R}^D`:
```math
k(x, x') = \prod_{i=1}^D k(x_i, x'_i).
```
"""
@inline iskroncompatible(κ::Kernel) = false # Default return for kernels

function checkkroncompatible(κ::Kernel)
    return iskroncompatible(κ) || throw(
        ArgumentError(
            "The chosen kernel is not compatible for Kronecker matrices (see [`iskroncompatible`](@ref))",
        ),
    )
end

function _kernelmatrix_kroneckerjl_helper(
    ::Type{<:MOInputIsotopicByFeatures}, Kfeatures, Koutputs
)
    return Kronecker.kronecker(Kfeatures, Koutputs)
end

function _kernelmatrix_kroneckerjl_helper(
    ::Type{<:MOInputIsotopicByOutputs}, Kfeatures, Koutputs
)
    return Kronecker.kronecker(Koutputs, Kfeatures)
end

"""
    kronecker_kernelmatrix(
        k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI, y::MOI
    ) where {MOI<:IsotopicMOInputsUnion}

Requires Kronecker.jl: Computes the `kernelmatrix` for the `IndependentMOKernel` and the
`IntrinsicCoregionMOKernel`, but returns a lazy kronecker product. This object can be very
efficiently inverted or decomposed. See also [`kernelmatrix`](@ref).
"""
function kronecker_kernelmatrix(
    k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI, y::MOI
) where {MOI<:IsotopicMOInputsUnion}
    x.out_dim == y.out_dim ||
        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
    Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
    Koutputs = _mo_output_covariance(k, x.out_dim)
    return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
end

function kronecker_kernelmatrix(
    k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI
) where {MOI<:IsotopicMOInputsUnion}
    Kfeatures = kernelmatrix(k.kernel, x.x)
    Koutputs = _mo_output_covariance(k, x.out_dim)
    return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
end

function kronecker_kernelmatrix(
    k::MOKernel, x::IsotopicMOInputsUnion, y::IsotopicMOInputsUnion
)
    return throw(
        ArgumentError("This kernel does not support a lazy kronecker kernelmatrix.")
    )
end

function kronecker_kernelmatrix(k::MOKernel, x::IsotopicMOInputsUnion)
    return kronecker_kernelmatrix(k, x, x)
end


text:
# ---
# title: 718. Maximum Length of Repeated Subarray
# id: problem718
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: Array, Hash Table, Binary Search, Dynamic Programming
# link: <https://leetcode.com/problems/maximum-length-of-repeated-subarray/description/>
# hidden: true
# ---
#
# Given two integer arrays `A` and `B`, return the maximum length of an subarray
# that appears in both arrays.
#
# **Example 1:**
#
#
#
#     Input:
#     A: [1,2,3,2,1]
#     B: [3,2,1,4,7]
#     Output: 3
#     Explanation:
#     The repeated subarray with maximum length is [3, 2, 1].
#
#
#
#
# **Note:**
#
#   1. 1 <= len(A), len(B) <= 1000
#   2. 0 <= A[i], B[i] < 100
#
#
#
#
## @lc code=start
using LeetCode

function longest_common_substring(itr1, itr2)
    m, n = length(itr1) + 1, length(itr2) + 1
    dp = fill(0, m, n)
    for i in 2: m, j in 2: n
        (itr1[i - 1] == itr2[j - 1]) && (dp[i, j] = dp[i - 1, j - 1] + 1)
    end
    return maximum(dp)
end

find_length718(nums1::Vector{Int}, nums2::Vector{Int}) = longest_common_substring(nums1, nums2)
## @lc code=end


text:
"""Functions for pretty-printing various objects and results"""

"Removes trailing zeros after a decimal place, and returns empty for 1.0"
function _strip_decimal_string(str)
    decimal_idx = findfirst('.', str)
    chop_idx = length(str)
    if !isnothing(decimal_idx)
        for i in length(str):-1:decimal_idx-1
            chop_idx = i
            cur_char = str[chop_idx]

            if (cur_char != '0' && cur_char != '.')
                break
            end
        end
    end
    # If we're left with just a 1 or -1, remove the 1
    if (chop_idx == decimal_idx-1 && str[chop_idx] == '1')
        chop_idx -= 1
    end
    return str[1:chop_idx]
end

"Converts a list of basis elements for a J matrix into a nice string summary"
function _coupling_basis_strings(coup_basis; digits, tol) :: Matrix{String}
    J = fill("", 3, 3)
    for (letter, basis_mat) in coup_basis
        for idx in eachindex(basis_mat)
            coeff = basis_mat[idx]
            if abs(coeff) > tol
                coeff = round(coeff; digits)
                if J[idx] == ""
                    float_str = @sprintf "%.4f" coeff
                else
                    float_str = @sprintf "%+.4f" coeff
                end
                float_str = _strip_decimal_string(float_str)
                J[idx] *= float_str * letter
            end
        end
    end
    for idx in eachindex(J)
        if J[idx] == ""
            J[idx] = "0"
        end
    end
    return J
end

function _print_allowed_coupling(basis_strs; prefix)
    max_len = maximum(length, basis_strs)
    minus_prefix = any(s -> startswith(s, '-'), basis_strs)
    if !minus_prefix
        max_len += 1
    end

    for i in 1:3
        print(i == 1 ? prefix : repeat(' ', length(prefix)))
        print('|')
        for j in 1:3
            elem = basis_strs[i, j]
            elem = repeat(' ', max_len-length(elem)) * elem
            print(elem * " ")
        end
        println('|')
    end
end

function print_allowed_coupling(cryst::Crystal, b::Bond; prefix="", digits=2, tol=1e-4)
    basis = basis_for_symmetry_allowed_couplings(cryst, b)
    basis_strs = _coupling_basis_strings(zip('A':'Z', basis); digits, tol)
    _print_allowed_coupling(basis_strs; prefix)
end

"""
    print_bond(cryst::Crystal, bond::Bond)
    print_bond(cryst::Crystal, i::Int)

Pretty-prints symmetry information for bond `bond` or atom index `i`.
"""
function print_bond(cryst::Crystal, b::Bond; digits=2, tol=1e-4)
    ri = cryst.positions[b.i]
    rj = cryst.positions[b.j] + b.n

    basis = basis_for_symmetry_allowed_couplings(cryst, b)
    basis_strs = _coupling_basis_strings(zip('A':'Z', basis); digits, tol)

    if b.i == b.j && iszero(b.n)
        # On site interaction
        class_i = cryst.classes[b.i]
        m_i = count(==(class_i), cryst.classes)
        if isempty(cryst.types[b.i])
            @printf "Atom %d, position [%.4g, %.4g, %.4g], multiplicity %d\n" b.i ri[1] ri[2] ri[3] m_i
        else
            @printf "Atom %d, type '%s', position [%.4g, %.4g, %.4g], multiplicity %d\n" b.i cryst.types[b.i] ri[1] ri[2] ri[3] m_i
        end

        _print_allowed_coupling(basis_strs; prefix="Allowed single-ion anisotropy or g-tensor: ")
        if any(J -> !(J ≈ J'), basis)
            println("""Note: The antisymmetric part is irrelevant to the single-ion anisotropy,
                       but could contribute meaningfully to the g-tensor.""")
        end
    else
        @printf "Bond(%d, %d, [%d, %d, %d])\n" b.i b.j b.n[1] b.n[2] b.n[3]
        (m_i, m_j) = (coordination_number(cryst, b.i, b), coordination_number(cryst, b.j, b))
        if m_i == m_j
            @printf "Distance %.4g, coordination %d\n" distance(cryst, b) m_i
        else
            @printf "Distance %.4g, coordination %d (from atom %d) and %d (from atom %d)\n" distance(cryst, b) m_i b.i m_j b.j
        end
        if isempty(cryst.types[b.i]) && isempty(cryst.types[b.j])
            @printf "Connects [%.4g, %.4g, %.4g] to [%.4g, %.4g, %.4g]\n" ri[1] ri[2] ri[3] rj[1] rj[2] rj[3]
        else
            @printf "Connects '%s' at [%.4g, %.4g, %.4g] to '%s' at [%.4g, %.4g, %.4g]\n" cryst.types[b.i] ri[1] ri[2] ri[3] cryst.types[b.j] rj[1] rj[2] rj[3]
        end
        _print_allowed_coupling(basis_strs; prefix="Allowed exchange matrix: ")
        antisym_basis_idxs = findall(J -> J ≈ -J', basis)
        if !isempty(antisym_basis_idxs)
            antisym_basis_strs = _coupling_basis_strings(collect(zip('A':'Z', basis))[antisym_basis_idxs]; digits, tol)
            println("Allowed DM vector: [$(antisym_basis_strs[2,3]) $(antisym_basis_strs[3,1]) $(antisym_basis_strs[1,2])]")
        end
    end
    println()
end

function print_bond(cryst::Crystal, i::Int; digits=2, tol=1e-4)
    print_bond(cryst, Bond(i, i, [0, 0, 0]); digits, tol)
end

"""
    print_bond_table(cryst::Crystal, max_dist)

Pretty-prints a table of bonds, one for each symmetry equivalence class, up to a
maximum bond length of `max_dist`. Equivalent to calling `print_bond(cryst, b)`
for every bond `b` in `reference_bonds(cryst, max_dist)`.
"""
function print_bond_table(cryst::Crystal, max_dist; digits=2, tol=1e-4)
    for b in reference_bonds(cryst, max_dist)
        print_bond(cryst, b; digits, tol)
    end
end


text:
@reexport module Plotters

using Plots
using GR
using Interact
using DataFrames

export fit!,transform!
export Plotter

import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload

using TSML.TSMLTypes
using TSML.ValDateFilters
using TSML.Utils

# setup plotting for publication
function setupplot(pdfoutput::Bool)
  Plots.gr();
  fntsm = Plots.font("sans-serif", 8);
  fntlg = Plots.font("sans-serif", 8);
  Plots.default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm);
  if pdfoutput == true
    Plots.default(size=(390,200)); #Plot canvas size
  else
    Plots.default(size=(500,300)); #Plot canvas size
  end
  return nothing
end

"""
    Plotter(
       Dict(
           :interactive => false,
           :pdfoutput => false
       )
    )

Plots a TS by default but performs interactive plotting if specified during instance creation.
- `:interactive` => boolean to indicate whether to use interactive plotting with `false` as default
- `:pdfoutput` => boolean to indicate whether ouput will be saved as pdf with `false` as default

Example:

    csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
    pltr = Plotter(Dict(:interactive => false))

    mpipeline = Pipeline(Dict(
         :transformers => [csvfilter,pltr]
       )
    )
    fit!(mpipeline)
    myplot = transform!(mpipeline)

Implements: `fit!`, `transform!`
"""
mutable struct Plotter <: Transformer
  model
  args
  function Plotter(args=Dict())
    default_args = Dict(
        :interactive => false,
        :pdfoutput => false
    )
    margs=mergedict(default_args, args)
    setupplot(margs[:pdfoutput])
    new(nothing,margs)
  end
end

"""
    fit!(pltr::Plotter, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}

Check validity of `features`: 2-column Date,Val data
"""
function fit!(pltr::Plotter, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  typeof(features) <: DataFrame || error("Outliernicer.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  pltr.model = pltr.args
end

"""
    transform!(pltr::Plotter, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Convert `missing` into `NaN` to allow plotting of discontinuities.
"""
function transform!(pltr::Plotter, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  features != [] || return DataFrame()
  typeof(features) <: DataFrame || error("Outliernicer.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  sum(names(features) .== (:Date,:Value))  == 2 || error("wrong column names")
  # covert missing to NaN
  df = deepcopy(features)
  df.Value = Array{Union{Missing, Float64,eltype(features.Value)},1}(missing,nrow(df))
  df.Value .= features.Value
  ndxmissing = findall(x->ismissing(x),df.Value)
  df.Value[ndxmissing] .= NaN
  setupplot(pltr.args[:pdfoutput])
  if pltr.args[:interactive] == true && pltr.args[:pdfoutput] == false
    interactiveplot(df)
  else
    pl=Plots.plot(df.Date,df.Value,xlabel="Date",ylabel="Value",legend=false,show=false);
    return pl
  end
end

function interactiveplot(df::Union{Vector,Matrix,DataFrame})
  mlength = length(df.Value)
  @manipulate for min in slider(1:mlength,label="min",value=1),max in slider(1:mlength,label="max",value=mlength)
     Plots.plot(df.Date[min:max],df.Value[min:max],xlabel="Date",ylabel="Value",legend=false,show=false);
  end
end

end


text:
abstract type AbstractConstraint end

mutable struct RessourceConstraint <: AbstractConstraint
    capacity::Float64
    variablecosts::IdDict{Any,Float64}
    utilization_log::Vector{Float64}
    setup_log::Vector{Int}
    name::String
end
RessourceConstraint(capacity, pairs; name = "") = RessourceConstraint(capacity, IdDict{Any,Float64}(pairs), zeros(0), zeros(0), name)
RessourceConstraint(capacity, pairs...; name = "") = RessourceConstraint(capacity, pairs, name = name)

function (rc::RessourceConstraint)()
    consumption = 0.0
    setups = 0
    for (element, cost) in rc.variablecosts
        q = sum(values(element.pull_orders))
        consumption += q*cost
        setups += q > 0
    end
    push!(rc.utilization_log, consumption/rc.capacity)
    push!(rc.setup_log, setups)
    if consumption > rc.capacity
        ratio = rc.capacity/consumption
        for (element, cost) in rc.variablecosts
            for destination in keys(element.pull_orders)
                element.pull_orders[destination] *= ratio
            end
        end
    end
end

function RLBase.reset!(rc::RessourceConstraint)
    empty!(rc.setup_log)
    empty!(rc.utilization_log)
end

function Base.show(io::IO, rc::RessourceConstraint)
    print(io, "Ressource $(rc.name): ")
    for (el, cost) in rc.variablecosts
        print(io, cost,"*",el.name,"+")
    end
    println(io, "\b≤",rc.capacity)
end
#=
struct RessourceConstraintSetup <: AbstractConstraint
    capacity::Float64
    variablecosts::Dict{Any,Float64}
    setupcosts::Dict{Any,Float64}
end

capacity constraint (inventories)
=#


text:
struct Pointer
    sym::String
    offset::UInt32
    pos::Char
    source::String
    target::String
end

Pointer(sym, offset, pos, source) = Pointer(
    sym, offset, pos, source[1:2], source[3:4]
)

is_semantic(ptr::Pointer) = ptr.source == "00" && ptr.target == "00"


text:
using Unitful

struct Atmosphere_test
    z::Vector{<:Unitful.Length}
    x::Vector{<:Unitful.Length}
    y::Vector{<:Unitful.Length}
    # temperature::Array{<:Unitful.Temperature, 3}
    # electron_density::Array{<:NumberDensity, 3}
    # hydrogen_populations::Array{<:NumberDensity, 3}
    nz::Integer
    nx::Integer
    ny::Integer
    Δx::Unitful.Length
    Δy::Unitful.Length
    function Atmosphere_test(z::Vector{<:Unitful.Length{T}},
                             x::Vector{<:Unitful.Length{T}},
                             y::Vector{<:Unitful.Length{T}}) where T <: AbstractFloat
                        # temperature::Array{<:Unitful.Temperature{T}, 3},
                        # electron_density::Array{<:NumberDensity{T}, 3},
                        # hydrogen_populations::Array{<:NumberDensity{T}, 3}) where T <: AbstractFloat
        #Funcy func
        nz = length(z)
        nx = length(x)
        ny = length(y)
        Δx = x[2] - x[1]
        Δy = y[2] - y[1]
        new{T}(z, x, y, nz, nx, ny, Δx, Δy)
               # temperature, electron_density, hydrogen_populations,
               # nz, nx, ny, Δx, Δy)
    end
end

z = collect(LinRange(0,1,11))u"m"
x = copy(z)
y = copy(z)

T = 100*rand(11, 3)u"K"

atmos = Atmosphere(z, x, y)

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
include("experiment_helpers.jl")

# Create Roadway
roadway_opts = MergingRoadwayOptions(merging_point=1.4, lane_width=0.28)
roadway = build_roadway(roadway_opts)

# Create players
T = Float64
p = 3
model = BicycleGame(p=p)
n = model.n
m = model.m

# Define the horizon of the problem
N = 30 # N time steps
dt = 0.10 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(3*SVector{model.ni[i],T}([1., 0.8, 1., 1.])) for i=1:p] # Quadratic state cost
R = [Diagonal(0.3*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2.0, 0.08, 0.0, 0]),
      SVector{model.ni[2],T}([2.0, 0.08, 0.0, 0]),
      SVector{model.ni[2],T}([2.0, 0.08, 0.0, 0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
radius = 0.20*ones(p)
μ = 10.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)
radius = 0.08
add_collision_avoidance!(game_con, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.2,     0.2,    0.2,
   -0.42,   -0.10,   0.16,
    0.45,    0.60,   0.60,
    0.15,    0.0,    0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-5
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

players = Vector{Player{T}}(undef, p)
players[1] = Player(model, roadway.lane[3])
players[2] = Player(model, roadway.lane[3])
players[3] = Player(model, roadway.lane[3])

# Create Scenario
sce = Scenario(model, roadway, players)

# Initialize visualizers
vis = Visualizer()
open(vis)

# Visualize trajectories
set_scenario!(vis, sce)
set_env!(vis, VehicleState(1.3, 0.0, 0.0, 0.0))
set_camera_birdseye!(vis, height=5.0)

build_waypoint!(vis, sce.player, N, key=0)
set_waypoint_traj!(vis, model, sce, prob.pdtraj.pr, key=0)
set_traj!(vis, model, sce, prob.pdtraj.pr)

display_ghost(vis, prob, sce)

# Verify the invertibility of the Jacobian of the KKT conditions.
residual_jacobian!(prob, prob.pdtraj)
null = nullspace(Matrix(prob.core.jac))
dim = size(null)[2]


text:
const dimensions = 3
const screen_z = 250

"""
Takes a point in 3D space and calculates the position of its "camera" image.
Returns nothing if the point is outside the camera's field of view.

#
- `point::Tuple3{Float32}`: Point in 3D space
- `screen_size::Tuple2{UInt16}`: Width & height of the snapshot screen
"""
function screen_image(point, screen_size)
    w, h = screen_size[1] ÷ 2, screen_size[2] ÷ 2
    # If the point is behind the screen, there can be no image
    if point[3] < screen_z
        return nothing
    end
    # Scale the vector such that it's z-coordinate matches the screen
    scaled = point ./ point[3] .* screen_z
    # Ensure the scaled vector is positioned within the borders of the screen
    if -w <= scaled[1] < w && -h <= scaled[2] < h
        return floor(Int32, scaled[1]), floor(Int32, scaled[2])
    end
    return nothing
end

"""
Determines if a point on a 2-sphere is visible to the "camera".

#
- `point::Tuple3{Float32}`: Point on a 2-sphere in 3D space
- `x::Integer`: Horizontal position/index of the image pixel
- `center::Tuple3{Float32}`: Center position of the sphere in 3D space
- `radius::Float32`: Radius of the sphere
- `screen_size::Tuple2{UInt16}`: Width & height of the snapshot screen
"""
function is_visible(point, center, radius, screen_size)
    # Important: We assume the point is on the sphere!
    if isnothing(screen_image(point, screen_size)) return false end
    # We interpret the point as a vector that can be arbitrarily scaled
    # to make a line through the origin an that point.
    # That way it serves as the parameter to the function
    # from R to R^3 in the form f(t) := t * point.
    # To get the quadratic equation of the intersection point of this line
    # with the sphere in the form 0 = a*t^2 + b*t + c, we calculate:
    a = sum(point.^2)
    b = -2 * sum(point .* center)
    c = sum(center.^2) - radius^2
    # The discriminant:
    d = b^2 - 4(a*c)
    # If `d` is less than 0, there is no point of intersection, i.e. the point
    # is not on the sphere at all.
    if d < 0 return false end
    # If `d` is 0, the line is tangent to the sphere,
    # i.e. there is exactly one point of intersection with the sphere,
    # which can only be our `point`.
    if d == 0 return true end
    # There are two solutions and thus two point of intersection,
    # one of them must be trivially 1, i.e. the `point` we are investigating.
    # The solutions to the quadratic equation:
    s1 = (-b + sqrt(d)) / 2a
    s2 = (-b - sqrt(d)) / 2a
    # Since we cannot simply compare with 1 due to float inaccuracy
    # we check which one of the two scalars is closer to 1:
    if abs(1 - s1) > abs(1 - s2)
        other_point = point .* s1
    else
        other_point = point .* s2
    end
    # We know the z-coordinate of our `point` is greater than or equal to the
    # z-coordinate of the screen, or the function would have returned early on.
    # If the z-coordinate of the `other_point` is in between the two,
    # it means the `other_point` is on the sphere's face closer to the screen,
    # so our `point` of interest is *not* visible.
    # Our `point` is only visible either if it is "in front" of the other,
    # or the other is behind the focus point (i.e. camera inside the sphere)
    return other_point[3] > point[3] || other_point[3] < 0
end

"""
Projects a pixel's position from its original 2D image plane
onto a 2-sphere in 3D space and returns the new sampled positions.

#
- `x::Integer`: Horizontal position/index of the image pixel
- `y::Integer`: Vertical position/index of the image pixel
- `w::Integer`: Width of the original image in pixels
- `h::Integer`: Height of the original image in pixels
- `center::Tuple3{Float32}`: Center position of the sphere in 3D space
- `radius::Float32`: Radius of the sphere
- `density::Integer`: Sampling density factor
"""
function samples(x, y, w, h, center, radius, density)
    if density == 1
        theta_phi = [(
            x * π/h,   # index 1 is theta
            y * 2π/w   # index 2 is phi
        )]
    elseif density > 1
        # Generate a `density x 2`-matrix of random floats in the range [0, 1[:
        deltas = rand(Float32, density, 2)
        # Construct an array of theta-phi-tuples with the length of `density`:
        theta_phi = [ (
            (x + deltas[i, 1]) * π/h,  # index 1 is theta with random delta from [0, 1[
            (y + deltas[i, 2]) * 2π/w  # index 2 is phi with random delta from [0, 1[
        ) for i=1:density ]
    else
        error("Invalid sampling density $density. Must be a positive integer.")
    end
    # Return the generated point positions on the sphere as an array of 3-tuples:
    return [ (
        center[1] + radius * sin(theta_phi[i][1]) * cos(theta_phi[i][2]),
        center[2] + radius * sin(theta_phi[i][1]) * sin(theta_phi[i][2]),
        center[3] + radius * cos(theta_phi[i][1])
    ) for i=1:density ]
end

"""
Takes RGBA encoded image data, projects it onto a sphere and takes a snapshot.

# Arguments
- `w::Integer`: Width of the original image
- `h::Integer`: Height of the original image
- `data::Array{Array4{UInt8}}`: Image data in a 2D array of integers 0-255
- `center::Tuple3{Float32}`: Center position of the sphere in 3D space
- `radius::Float32`: Radius of the sphere
- `density::Integer`: Sampling density factor
- `snapshot_size::Tuple2{UInt16}`: Width & height of the snapshot screen
"""
function snapshot_sphere(w, h, data, center, radius, density, snapshot_size)
    # Initialize 2D array with a length of height * witdth of the snapshot in
    # the first dimension and a length of 4 (for RGBA) in the second dimension.
    snap_w, snap_h = snapshot_size
    screen = zeros(UInt8, (snap_w * snap_h, 4))
    # Go through each row of the 2D data array, i.e. each pixel.
    for (i, pixel) in enumerate(eachrow(data))
        # Get the x-y-position of the pixel in the original image.
        x = 1 + (i - 1) % w
        y = 1 + (i - 1) ÷ w
        # Sample the the position onto the sphere.
        points = samples(x, y, w, h, center, radius, density)
        for point in points
            # If a point on the sphere is visible, get its screen image,
            # and save the corresponding pixel value in the screen array
            # at the appropriate index.
            if is_visible(point, center, radius, snapshot_size)
                img_x, img_y = screen_image(point, snapshot_size)
                idx = img_x + snap_w ÷ 2 + (img_y + snap_h ÷ 2) * snap_w
                screen[idx, :] = pixel
            end
        end
    end
    return screen
end

"""
Parse and convert command line arguments, read image data,
perform sanity checks, and start the calculations.
Image data must be supplied via standard input.
"""
function main()
    if length(ARGS) != 5
        error(
            "5 arguments required: " *
            "image_size, center, radius, density, snapshot_size"
        )
    end
    # Parse and convert string arguments to usable types:
    w, h = parse.(Int32, split(ARGS[1], ","))
    center = parse.(Float32, split(ARGS[2], ","))
    radius = parse(Float32, ARGS[3])
    density = parse(Int32, ARGS[4])
    snap_size = parse.(Int32, split(ARGS[5], ","))
    # Perform sanity checks, starting with the dimensions of the center point:
    if length(center) != dimensions
        error("Center point must be a $dimensions-tuple of numbers")
    end
    # Read the image data as a bytes sequence from standard input;
    # this assumes every pixel to be represented by exactly 4 bytes of data,
    # which correspond to its RGBA value:
    data = read(stdin, typemax(Int32))  # cap at an equivalent of 2 GB
    if sizeof(data) % 4 != 0
        error("Number of provided bytes must be a multiple of 4!")
    end
    # Make a vector of 4-byte-chunks out of the provided data:
    pixels = transpose(reshape(data, 4, :))
    # The number of 4-byte-chunks must correspond to the specified dimensions:
    num_pixels = size(pixels)[1]
    if num_pixels != w * h
        error(
            "Number of 4-byte strides $num_pixels " *
            "does not match provided width and height of $w x $h"
        )
    end
    # Start the calculations:
    data = snapshot_sphere(w, h, pixels, center, radius, density, snap_size)
    # Return to stdout:
    write(stdout, transpose(data))
end

main()

max_stars_repo_name: daniil-berg/spherify

max_issues_repo_name: daniil-berg/spherify

max_forks_repo_name: daniil-berg/spherify


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using BasicPOMCP
using POMCPOW
using ParticleFilters
using CPUTime
using DataFrames
using ARDESPOT
using Distributions
using DiscreteValueIteration

@everywhere begin
    using POMDPs
    using ParticleFilters
    using POMDPToolbox
    using ContinuousPOMDPTreeSearchExperiments
end

pomdp = SimpleLightDark()

function gen_sims(x::Vector{Float64}, n, k, seed)
    c = max(0.0, x[1])
    k_obs = x[2]
    @assert k_obs >= 1.0
    inv_alpha_obs = x[3]
    @assert inv_alpha_obs >= 0.1

    sims = []

    for i in 1:n
        rng = MersenneTwister(0)
        ro = ValueIterationSolver()
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(c),
                               final_criterion=MaxTries(),
                               max_depth=90,
                               max_time=1.0,
                               k_observation=k_obs,
                               alpha_observation=1/inv_alpha_obs,
                               estimate_value=FOValue(ro),
                               enable_action_pw=false,
                               check_repeat_obs=true,
                               # default_action=ReportWhenUsed(TagAction(false, 0.0)),
                               # default_action=TagAction(false, 0.0),
                               rng=rng
                              )

        planner = solve(solver, pomdp)
        filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(10_000),
                                           0.05,
                                           MersenneTwister(i+10000*k))

        srand(planner, i+40000*k)
        sim = Sim(deepcopy(pomdp),
                  deepcopy(planner),
                  filter,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=100,
                  metadata=Dict(:i=>i, :k=>k)
                 )
        push!(sims, sim)
    end

    return sims
end

start_mean = [115.0, 4.0, 14.0]
start_cov = diagm([25.0^2, 3.0^2, 5.0^2])
d = MvNormal(start_mean, start_cov)
K = 60  # 60 # number of parameter samples
n = 100 # 100 # number of evaluation simulations
m = 15  # 15 # number of elite samples
max_iters = 100

for i in 1:max_iters
    sims = []
    params = Vector{Vector{Float64}}(K)
    print("creating $K simulation sets")
    for k in 1:K
        p = rand(d)
        p[2] = max(1.0, p[2])
        p[3] = max(0.1, p[3])
        params[k] = p
        k_sims = gen_sims(p, n, k, i)
        print(".")
        append!(sims, k_sims)
    end
    println()
    results = run_parallel(sims)
    # results = run(sims)
    combined = by(results, :k) do df
        DataFrame(mean_reward=mean(df[:reward]))
    end
    @show mean(combined[:mean_reward])
    order = sortperm(combined[:mean_reward])
    elite = params[combined[:k][order[K-m:end]]]
    elite_matrix = Matrix{Float64}(length(start_mean), m)
    for k in 1:m
        elite_matrix[:,k] = elite[k]
    end
    try
        d = fit(typeof(d), elite_matrix)
    catch ex
        if ex isa Base.LinAlg.PosDefException
            println("pos def exception")
            d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
        else
            rethrow(ex)
        end
    end
    println("iteration $i")
    @show mean(d)
    @show det(cov(d))
    @show ev = eigvals(cov(d))
    for j in 1:length(ev)
        @show eigvecs(cov(d))[:,j]
    end
end


text:
using StatsFuns
using DistributionsAD
using Bijectors
using Bijectors: TransformedDistribution
using Random: AbstractRNG, GLOBAL_RNG
import Bijectors: bijector

update(d::TuringDiagMvNormal, μ, σ) = TuringDiagMvNormal(μ, σ)
update(td::TransformedDistribution, θ...) = transformed(update(td.dist, θ...), td.transform)
function update(td::TransformedDistribution{<:TuringDiagMvNormal}, θ::AbstractArray)
    μ, ω = θ[1:length(td)], θ[length(td) + 1:end]
    return update(td, μ, softplus.(ω))
end

"""
    bijector(model::Model; sym_to_ranges = Val(false))

Returns a `Stacked <: Bijector` which maps from the support of the posterior to ℝᵈ with `d`
denoting the dimensionality of the latent variables.
"""
function bijector(model::Model; sym_to_ranges::Val{sym2ranges} = Val(false)) where {sym2ranges}
    varinfo = Turing.VarInfo(model)
    num_params = sum([size(varinfo.metadata[sym].vals, 1)
                      for sym ∈ keys(varinfo.metadata)])

    dists = vcat([varinfo.metadata[sym].dists for sym ∈ keys(varinfo.metadata)]...)

    num_ranges = sum([length(varinfo.metadata[sym].ranges)
                      for sym ∈ keys(varinfo.metadata)])
    ranges = Vector{UnitRange{Int}}(undef, num_ranges)
    idx = 0
    range_idx = 1

    # ranges might be discontinuous => values are vectors of ranges rather than just ranges
    sym_lookup = Dict{Symbol, Vector{UnitRange{Int}}}()
    for sym ∈ keys(varinfo.metadata)
        sym_lookup[sym] = Vector{UnitRange{Int}}()
        for r ∈ varinfo.metadata[sym].ranges
            ranges[range_idx] = idx .+ r
            push!(sym_lookup[sym], ranges[range_idx])
            range_idx += 1
        end

        idx += varinfo.metadata[sym].ranges[end][end]
    end

    bs = inv.(bijector.(tuple(dists...)))

    if sym2ranges
        return Stacked(bs, ranges), (; collect(zip(keys(sym_lookup), values(sym_lookup)))...)
    else
        return Stacked(bs, ranges)
    end
end

"""
    meanfield(model::Model)
    meanfield(rng::AbstractRNG, model::Model)

Creates a mean-field approximation with multivariate normal as underlying distribution.
"""
meanfield(model::Model) = meanfield(GLOBAL_RNG, model)
function meanfield(rng::AbstractRNG, model::Model)
    # setup
    varinfo = Turing.VarInfo(model)
    num_params = sum([size(varinfo.metadata[sym].vals, 1)
                      for sym ∈ keys(varinfo.metadata)])

    dists = vcat([varinfo.metadata[sym].dists for sym ∈ keys(varinfo.metadata)]...)

    num_ranges = sum([length(varinfo.metadata[sym].ranges)
                      for sym ∈ keys(varinfo.metadata)])
    ranges = Vector{UnitRange{Int}}(undef, num_ranges)
    idx = 0
    range_idx = 1
    for sym ∈ keys(varinfo.metadata)
        for r ∈ varinfo.metadata[sym].ranges
            ranges[range_idx] = idx .+ r
            range_idx += 1
        end

        # append!(ranges, [idx .+ r for r ∈ varinfo.metadata[sym].ranges])
        idx += varinfo.metadata[sym].ranges[end][end]
    end

    # initial params
    μ = randn(rng, num_params)
    σ = softplus.(randn(rng, num_params))

    # construct variational posterior
    d = TuringDiagMvNormal(μ, σ)
    bs = inv.(bijector.(tuple(dists...)))
    b = Stacked(bs, ranges)

    return transformed(d, b)
end

"""
$(TYPEDEF)

Automatic Differentiation Variational Inference (ADVI) with automatic differentiation
backend `AD`.

# Fields

$(TYPEDFIELDS)
"""
struct ADVI{AD} <: VariationalInference{AD}
    "Number of samples used to estimate the ELBO in each optimization step."
    samples_per_step::Int
    "Maximum number of gradient steps."
    max_iters::Int
end

"""
    ADVI([samples_per_step=1, max_iters=1000])

Create an [`ADVI`](@ref) with the currently enabled automatic differentiation backend
`ADBackend()`.
"""
function ADVI(samples_per_step::Int=1, max_iters::Int=1000)
    return ADVI{ADBackend()}(samples_per_step, max_iters)
end

alg_str(::ADVI) = "ADVI"

function vi(model::Model, alg::ADVI; optimizer = TruncatedADAGrad())
    q = meanfield(model)
    return vi(model, alg, q; optimizer = optimizer)
end

function vi(model, alg::ADVI, q::TransformedDistribution{<:TuringDiagMvNormal}; optimizer = TruncatedADAGrad())
    Turing.DEBUG && @debug "Optimizing ADVI..."
    # Initial parameters for mean-field approx
    μ, σs = params(q)
    θ = vcat(μ, invsoftplus.(σs))

    # Optimize
    optimize!(elbo, alg, q, model, θ; optimizer = optimizer)

    # Return updated `Distribution`
    return update(q, θ)
end

function vi(model, alg::ADVI, q, θ_init; optimizer = TruncatedADAGrad())
    Turing.DEBUG && @debug "Optimizing ADVI..."
    θ = copy(θ_init)
    optimize!(elbo, alg, q, model, θ; optimizer = optimizer)

    # If `q` is a mean-field approx we use the specialized `update` function
    if q isa Distribution
        return update(q, θ)
    else
        # Otherwise we assume it's a mapping θ → q
        return q(θ)
    end
end

function optimize(elbo::ELBO, alg::ADVI, q, model, θ_init; optimizer = TruncatedADAGrad())
    θ = copy(θ_init)

    if model isa Model
        optimize!(elbo, alg, q, make_logjoint(model), θ; optimizer = optimizer)
    else
        # `model` assumed to be callable z ↦ p(x, z)
        optimize!(elbo, alg, q, model, θ; optimizer = optimizer)
    end

    return θ
end

# WITHOUT updating parameters inside ELBO
function (elbo::ELBO)(
    rng::AbstractRNG,
    alg::ADVI,
    q::VariationalPosterior,
    logπ::Function,
    num_samples
)
    #   𝔼_q(z)[log p(xᵢ, z)]
    # = ∫ log p(xᵢ, z) q(z) dz
    # = ∫ log p(xᵢ, f(ϕ)) q(f(ϕ)) |det J_f(ϕ)| dϕ   (since change of variables)
    # = ∫ log p(xᵢ, f(ϕ)) q̃(ϕ) dϕ                   (since q(f(ϕ)) |det J_f(ϕ)| = q̃(ϕ))
    # = 𝔼_q̃(ϕ)[log p(xᵢ, z)]

    #   𝔼_q(z)[log q(z)]
    # = ∫ q(f(ϕ)) log (q(f(ϕ))) |det J_f(ϕ)| dϕ     (since q(f(ϕ)) |det J_f(ϕ)| = q̃(ϕ))
    # = 𝔼_q̃(ϕ) [log q(f(ϕ))]
    # = 𝔼_q̃(ϕ) [log q̃(ϕ) - log |det J_f(ϕ)|]
    # = 𝔼_q̃(ϕ) [log q̃(ϕ)] - 𝔼_q̃(ϕ) [log |det J_f(ϕ)|]
    # = - ℍ(q̃(ϕ)) - 𝔼_q̃(ϕ) [log |det J_f(ϕ)|]

    # Finally, the ELBO is given by
    # ELBO = 𝔼_q(z)[log p(xᵢ, z)] - 𝔼_q(z)[log q(z)]
    #      = 𝔼_q̃(ϕ)[log p(xᵢ, z)] + 𝔼_q̃(ϕ) [log |det J_f(ϕ)|] + ℍ(q̃(ϕ))

    # If f: supp(p(z | x)) → ℝ then
    # ELBO = 𝔼[log p(x, z) - log q(z)]
    #      = 𝔼[log p(x, f⁻¹(z̃)) + logabsdet(J(f⁻¹(z̃)))] + ℍ(q̃(z̃))
    #      = 𝔼[log p(x, z) - logabsdetjac(J(f(z)))] + ℍ(q̃(z̃))

    # But our `forward(q)` is using f⁻¹: ℝ → supp(p(z | x)) going forward → `+ logjac`
    _, z, logjac, _ = forward(rng, q)
    res = (logπ(z) + logjac) / num_samples

    if q isa TransformedDistribution
        res += entropy(q.dist)
    else
        res += entropy(q)
    end

    for i = 2:num_samples
        _, z, logjac, _ = forward(rng, q)
        res += (logπ(z) + logjac) / num_samples
    end

    return res
end


text:
mutable struct EnergiesT
    Total::Float64
    Kinetic::Float64
    Ionic::Float64
    Hartree::Float64
    XC::Float64
    NN::Float64
end

function print_Energies( Energies::EnergiesT )
    @printf("Kinetic energy: %18.10f\n", Energies.Kinetic )
    @printf("Ionic   energy: %18.10f\n", Energies.Ionic )
    @printf("Hartree energy: %18.10f\n", Energies.Hartree )
    @printf("XC      energy: %18.10f\n", Energies.XC )
    @printf("NN      energy: %18.10f\n", Energies.NN )
    @printf("----------------------------------\n")
    @printf("Total   energy: %18.10f\n", Energies.Total )
end


text:
# ensure stdlib is in load path
function ensure_stdlib()
    if "JULIA_LOAD_PATH" in keys(ENV)
        _add_path = false
        comps = split(ENV["JULIA_LOAD_PATH"], ":")
        for needed in split("@:@v#.#:@stdlib", ":")
            if !(needed in comps)
                push!(comps, needed)
                _add_path = true
            end
        end
        if _add_path
            @show ENV["JULIA_LOAD_PATH"] = join(comps, ":")
        end
    end
    nothing
end

if !Sys.iswindows()
    test_script = joinpath(@__DIR__, "testprotoc.sh")
    protogen_path = joinpath(@__DIR__, "..", "plugin")
    path_env = "$(protogen_path):$(ENV["PATH"])"
    modified_env = Dict{String, String}()
    is_ci = "CI" in keys(ENV)
    if is_ci
        modified_env["COVERAGE"] = "--code-coverage=user --inline=no"
    end
    protoc_compilers = []
    for protoc_env in ("PROTOC2", "PROTOC3")
        if protoc_env in keys(ENV)
            push!(protoc_compilers, ENV[protoc_env])
        end
    end
    if isempty(protoc_compilers)
        @info("no protobuf compilers setup for testing compiler plugin, skipping compiler tests.")
    else
        @info("testing protobuf compiler plugin...")
        is_ci && @info("detected CI environment, will enable code coverage...")
        for protoc_compiler in protoc_compilers
            modified_env["PROTOC"] = protoc_compiler
            modified_env["PATH"] = path_env
            modified_env["JULIA"] = joinpath(Sys.BINDIR, Base.julia_exename())
            @info("testing protoc compiler plugin with ", protoc_compiler)
            run(setenv(`$test_script`, modified_env))
        end
    end
end

# test ProtoBuf.protoc
test_protos = joinpath(@__DIR__, "proto")
test_proto = joinpath(test_protos, "plugin.proto")
ensure_stdlib()

@info("testing ProtoBuf.protoc")
ProtoBuf.protoc(`-I=$test_protos --julia_out=/tmp $test_proto`)


text:
@registerN Charlier AbstractCDOP1 μ
export Charlier

"""
    Charlier{μ}

References: [Koekoek and Swarttouw §1.12](https://arxiv.org/pdf/math/9602214.pdf)
"""
Charlier

abcde(::Type{<:Charlier{μ}}) where {μ} = NamedTuple{(:a, :b, :c, :d, :e)}((0, 1, 0, -1, μ))

basis_symbol(::Type{<:Charlier{μ}}) where {μ} = "cᵘ"
Polynomials.domain(::Type{<:Charlier{μ}}) where {μ} = Polynomials.Interval(0, Inf)
weight_function(::Type{<:Charlier{μ}}) where {μ} = x -> μ^x / gamma(1 + x)

function k1k0(P::Type{<:Charlier{μ}}, n::Int) where {μ}
    -one(eltype(P)) / μ
end

function classical_hypergeometric(P::Type{<:Charlier{μ}}, n::Int, x) where {μ}
    pFq((-n, -x), (), -1 / μ)
end


text:
function iterate!(job::BasicMCJob, ::Type{SliceSampler}, ::Type{Univariate})
  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.loguprime = log(rand())+job.pstate.logtarget
  job.sstate.lstate.value = job.pstate.value
  job.sstate.rstate.value = job.pstate.value
  job.sstate.primestate.value = job.pstate.value

  job.sstate.runiform = rand()
  job.sstate.lstate.value = job.pstate.value-job.sstate.runiform*job.sampler.widths[1]
  job.sstate.rstate.value = job.pstate.value+(1-job.sstate.runiform)*job.sampler.widths[1]

  if job.sampler.stepout
    job.parameter.logtarget!(job.sstate.lstate)

    while job.sstate.lstate.logtarget > job.sstate.loguprime
      job.sstate.lstate.value -= job.sampler.widths[1]
      job.parameter.logtarget!(job.sstate.lstate)
    end

    job.parameter.logtarget!(job.sstate.rstate)

    while job.sstate.rstate.logtarget > job.sstate.loguprime
      job.sstate.rstate.value += job.sampler.widths[1]
      job.parameter.logtarget!(job.sstate.rstate)
    end
  end

  while true
    job.sstate.primestate.value =
      rand()*(job.sstate.rstate.value-job.sstate.lstate.value)+job.sstate.lstate.value
    job.pstate.logtarget = job.parameter.logtarget!(job.sstate.primestate)
    if job.pstate.logtarget > job.sstate.loguprime
      break
    else
      if job.sstate.primestate.value > job.pstate.value
        job.sstate.rstate.value = job.sstate.primestate.value
      elseif job.sstate.primestate.value < job.pstate.value
        job.sstate.lstate.value = job.sstate.primestate.value
      else
        @assert false "Shrunk to current position and still not acceptable"
      end
    end
  end

  job.pstate.value = job.sstate.primestate.value

  if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
    if job.tuner.verbose
      println("Burnin iteration ", job.fmt_iter(job.sstate.tune.totproposed), " of ", job.range.burnin)
    end

    job.sstate.tune.totproposed += job.sstate.tune.proposed
    job.sstate.tune.proposed = 0
  end
end

function iterate!(job::BasicMCJob, ::Type{SliceSampler}, ::Type{Multivariate})
  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  for i = 1:job.pstate.size
    job.sstate.loguprime = log(rand())+job.pstate.logtarget
    job.sstate.lstate.value = copy(job.pstate.value)
    job.sstate.rstate.value = copy(job.pstate.value)
    job.sstate.primestate.value = copy(job.pstate.value)

    job.sstate.runiform = rand()
    job.sstate.lstate.value[i] = job.pstate.value[i]-job.sstate.runiform*job.sampler.widths[i]
    job.sstate.rstate.value[i] = job.pstate.value[i]+(1-job.sstate.runiform)*job.sampler.widths[i]

    if job.sampler.stepout
      job.parameter.logtarget!(job.sstate.lstate)

      while job.sstate.lstate.logtarget > job.sstate.loguprime
        job.sstate.lstate.value[i] -= job.sampler.widths[i]
        job.parameter.logtarget!(job.sstate.lstate)
      end

      job.parameter.logtarget!(job.sstate.rstate)

      while job.sstate.rstate.logtarget > job.sstate.loguprime
        job.sstate.rstate.value[i] += job.sampler.widths[i]
        job.parameter.logtarget!(job.sstate.rstate)
      end
    end

    while true
      job.sstate.primestate.value[i] =
        rand()*(job.sstate.rstate.value[i]-job.sstate.lstate.value[i])+job.sstate.lstate.value[i]
      job.pstate.logtarget = job.parameter.logtarget!(job.sstate.primestate)
      if job.pstate.logtarget > job.sstate.loguprime
        break
      else
        if job.sstate.primestate.value[i] > job.pstate.value[i]
          job.sstate.rstate.value[i] = job.sstate.primestate.value[i]
        elseif job.sstate.primestate.value[i] < job.pstate.value[i]
          job.sstate.lstate.value[i] = job.sstate.primestate.value[i]
        else
          @assert false "Shrunk to current position and still not acceptable"
        end
      end
    end

    job.pstate.value[i] = job.sstate.primestate.value[i]
  end

  if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
    if job.tuner.verbose
      println("Burnin iteration ", job.fmt_iter(job.sstate.tune.totproposed), " of ", job.range.burnin)
    end

    job.sstate.tune.totproposed += job.sstate.tune.proposed
    job.sstate.tune.proposed = 0
  end
end


text:
using Bridge, StaticArrays, Bridge.Models
const R = ℝ
Random.seed!(2)
Random.seed!(3)
πH = 2000. # prior

t = 1.0
T = 5.00
n = 50001 # total imputed length
m = 100 # number of segments
M = div(n-1,m)
skippoints = 2
dt = (T-t)/(n-1)
tt = t:dt:T
si = 3.
# 10, 20, 8/3 srand(2)
P = Bridge.Models.Lorenz(ℝ{3}(10, 28, 8/3), ℝ{3}(si,si,si))
P2 = Psmooth = Bridge.Models.Lorenz(ℝ{3}(10, 28, 8/3), ℝ{3}(0,0,0))

x0 = Models.x0(P)
#x0 =  ℝ{3}(6, 0, 2)

W = sample(tt, Wiener{ℝ{3}}())
X = SamplePath(tt, zeros(ℝ{3}, length(tt)))
Bridge.solve!(Euler(), X, x0, W, P)
W = sample(tt, Wiener{ℝ{3}}())
X2 = SamplePath(tt, zeros(ℝ{3}, length(tt)))
Bridge.solve!(Euler(), X2, x0, W, P2)

Xtrue = copy(X)

# Observation scheme and subsample
_pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
SV = ℝ{3}
SM = typeof(one(Bridge.outer(zero(SV))))
if isdefined(:partial)
    if !partial
        L = I
        Σ = SDiagonal(1., 1., 1.)
        lΣ = chol(Σ)'
        RV = SV
        RM = SM

        V = SamplePath(collect(_pairs(Xtrue))[1:M:end])
        map!(y -> L*y + lΣ*randn(RV), V.yy, V.yy)
    else
        L = @SMatrix [0.0 1.0 0.0; 0.0 0.0 1.0]
        Σ = SDiagonal(1., 1.)
        lΣ = chol(Σ)'
        RV = ℝ{2}
        RM = typeof(one(Bridge.outer(zero(RV))))

        V_ = SamplePath(collect(_pairs(Xtrue))[1:M:end])
        V = SamplePath(V_.tt, map(y -> L*y + lΣ*randn(RV), V_.yy))
    end
else
    println("Info:  No observation scheme defined. Set `partial in [true, false]`")
end