text:
module QuadPencils

using LinearAlgebra, SparseArrays, SuiteSparse

export linearize, deflate, eigbasis

### PencilScaling ##########################################################################
# Compute the optimal scaling
############################################################################################
struct PencilScaling{T}
    γ::T
    δ::T
end

function pencilscaling(A0::M, A1::M, A2::M; scaling = nothing, kw...) where {T,M<:AbstractMatrix{T}}
    scaling === nothing && return PencilScaling(T(1), T(1))
    n0, n1, n2 = opnorm.(Matrix.((A0, A1, A2)))
    τ = n1 / sqrt(n0 * n2)
    if τ <= 1
        γ = sqrt(n0 / n2)
        δ = 2 / (n0 + γ * n1)
    elseif scaling === -
        # distinct tropical roots, choose γ+ to favor large eigenvalues
        γ = n1 / n2
        δ = 1 / max(n2 * γ^2, n1 * γ, n0)
    else
        # distinct tropical roots, choose γ- to favor small eigenvalues
        γ = n0 / n1
        δ = 1 / max(n2 * γ^2, n1 * γ, n0)
    end
    return PencilScaling(T(γ), T(δ))
end

function scale!(A0, A1, A2, s::PencilScaling)
    γ, δ = scalingfactors(s)
    if γ != 1 || δ != 1
        A0 .*= γ^2 * δ
        A1 .*= γ * δ
        A2 .*= δ
    end
    return A0, A1, A2
end

scalingfactors(s::PencilScaling) = s.γ, s.δ

### QuadPencil #############################################################################
# Encodes the quadratic pencil problem
############################################################################################
struct QuadPencil{T<:Complex,M<:AbstractMatrix{T}}
    A0::M
    A1::M
    A2::M
    scaling::PencilScaling{T}
end

function quadpencil(A0, A1, A2; kw...)
    size(A0) == size(A1) == size(A2) && size(A0, 1) == size(A0, 2) ||
        throw(DimensionMismatch("Expected square matrices of the same size"))
    ptype = complex(promote_type(eltype(A0), eltype(A1), eltype(A2)))
    A0´ = copy_to_eltype(A0, ptype)
    A1´ = copy_to_eltype(A1, ptype)
    A2´ = copy_to_eltype(A2, ptype)
    s = pencilscaling(A0´, A1´, A2´; kw...)
    scale!(A0´, A1´, A2´, s)
    return QuadPencil(A0´, A1´, A2´, s)
end

# This also ensures any Adjoint/Transpose wrapper is removed
copy_to_eltype(A, ::Type{T}) where {T} = copy!(similar(A, T, size(A)), A)

Base.size(q::QuadPencil, n...) = size(q.A0, n...)

### Linearized Pencils #####################################################################
# Build second companion linearization, or its Q*C2*V rotation
############################################################################################
struct Linearization{T,M<:AbstractMatrix{T},R}
    pencil::QuadPencil{T}
    A::M
    B::M
    P::M
    deflate_tol::R
end

linearize(As...; kw...) = linearize(quadpencil(As...; kw...); kw...)

function linearize(p::QuadPencil{T}; type = 1) where {T}
    A, B = AB(p, type)
    P = one(A)
    tol = default_tol(T)
    return Linearization(p, A, B, P, tol)
end

AB(p, type) = type == 1 ? C1(p) : type == 2 ? C2(p) :
    throw(ArgumentError("Unrecognized linearization `type`"))

C1(p) = [0I I; p.A0 p.A1], [I 0I; 0I -p.A2]
C2(p) = [p.A1 -I; p.A0 0I], [-p.A2 0I; 0I -I]

function Base.show(io::IO, l::Linearization{T,M}) where {T,M}
    print(io, summary(l), "\n",
"  Matrix size    : $(size(l.A, 1)) × $(size(l.A, 2))
  Matrix type    : $M
  Scalings γ, δ  : $(real.(scalingfactors(l)))
  Deflated       : $(deflationstring(l))")
end

Base.summary(l::Linearization{T}) where {T} =
    "Linearization{T}: second companion linearization of quadratic pencil"

deflationstring(l::Linearization) =
    isdeflated(l) ? "true ($(size(l.P, 1)) -> $(size(l.A, 1)))" : "false"

isdeflated(l) = size(l.P, 1) != size(l.A, 1)

Base.size(l::Linearization, n...) = size(l.A, n...)

scalingfactors(l::Linearization) = scalingfactors(l.pencil.scaling)

### deflate #############################################################################
# Deflates zero and infinite eigenvalues in a pencil
############################################################################################
deflate(A0, A1, A2; kw...) = deflate(linearize(A0, A1, A2; kw...); kw...)

function deflate(l::Linearization)
    Adef, Bdef, Pdef = deflate_zeros_infinites(l.A, l.B, l.P; atol = l.deflate_tol)
    return Linearization(l.pencil, Adef, Bdef, Pdef, l.deflate_tol)
end

function deflate_zeros_infinites(A, B, P = I; kw...)
    Adef0, Bdef0, Pdef0 = deflate_zeros(A, B; kw...)
    Adef, Bdef, Pdef∞ = deflate_infinites(Adef0, Bdef0; kw...)
    Pdef = P * Pdef0 * Pdef∞
    return Adef, Bdef, Pdef
end

deflate_zeros(A::AbstractMatrix, B; kw...) = deflate_zeros(pqr(A), B; kw...)

function deflate_zeros(qrA::Factorization{T}, B::AbstractMatrix; atol = default_tol(T)) where {T}
    Q´B = getQ´(qrA, B)
    Q´A = getRP´(qrA)
    n = size(qrA, 1)
    r = nonzero_rows(Q´A, atol)
    A1 = view(Q´A, 1:r, :)
    B1 = view(Q´B, 1:r, :)
    B2 = view(Q´B, r+1:n, :)
    qrB2´ = pqr(B2')
    P = getQ_rev(qrB2´)   # (Q' * B2')' = B2 * Q
    Pdef = view(P, :, 1:r)
    Adef = A1 * Pdef
    Bdef = B1 * Pdef
    # Qdef = Matrix(I, r, n) * getQ´(qrA)
    return Adef, Bdef, Pdef
end

function deflate_infinites(A, B; kw...)
    Bdef, Adef, Pdef = deflate_zeros(pqr(B), A; kw...)
    return Adef, Bdef, Pdef
end

default_tol(::Type{T}) where {T} = sqrt(eps(real(T)))

### Tools ##################################################################################

pqr(a::SparseMatrixCSC) = qr(a)
pqr(a::Adjoint) = qr(copy(a))
pqr(a) = qr(a, Val(true))

getQ_rev(qr::Factorization) = qr.Q * Idense(size(qr, 1), reverse(1:size(qr,1)))
getQ_rev(qr::SuiteSparse.SPQR.QRSparse) =
    Isparse(size(qr, 1), :, qr.prow) * sparse(qr.Q * Idense(size(qr, 1), reverse(1:size(qr,1))))

getQ´(qr::Factorization) = qr.Q' * Idense(size(qr, 1))
getQ´(qr::Factorization, B) = qr.Q' * B
getQ´(qr::SuiteSparse.SPQR.QRSparse) =
    sparse(qr.Q' * Idense(size(qr, 1))) * Isparse(size(qr,1), qr.prow, :)
getQ´(qr::SuiteSparse.SPQR.QRSparse, B) =
    qr.Q' * Matrix(B)[qr.prow, :]

getRP´(qr::Factorization) = qr.R * qr.P'
getRP´(qr::SuiteSparse.SPQR.QRSparse) = qr.R * Isparse(size(qr, 2), qr.pcol, :)

Idense(n) = Matrix(I, n, n)

function Idense(n, cols)
    m = zeros(Bool, n, length(cols))
    for (j, col) in enumerate(cols)
        m[col, j] = true
    end
    return m
end

# Equivalent to I(n)[rows, cols], but faster
Isparse(n, rows, cols) = Isparse(inds(rows, n), inds(cols, n))

function Isparse(rows, cols)
    rowval = Int[]
    nzval = Bool[]
    colptr = Vector{Int}(undef, length(cols) + 1)
    colptr[1] = 1
    for (j, col) in enumerate(cols)
        push_rows!(rowval, nzval, rows, col)
        colptr[j+1] = length(rowval) + 1
    end
    return SparseMatrixCSC(length(rows), length(cols), colptr, rowval, nzval)
end

inds(::Colon, n) = 1:n
inds(is, n) = is

function push_rows!(rowval, nzval, rows::AbstractUnitRange, col)
    if col in rows
        push!(rowval, 1 + rows[1 + col - first(rows)] - first(rows))
        push!(nzval, true)
    end
    return nothing
end

function push_rows!(rowval, nzval, rows, col)
    for (j, row) in enumerate(rows)
        if row == col
            push!(rowval, j)
            push!(nzval, true)
        end
    end
    return nothing
end

function nonzero_rows(m::AbstractMatrix{T}, atol = sqrt(eps(real(T)))) where {T}
    n = 0
    for row in eachrow(m)
        all(z -> abs(z) < atol, row) && break
        n += 1
    end
    return n
end

function chop!(A::AbstractArray{T}, atol = sqrt(eps(real(T)))) where {T}
    for (i, a) in enumerate(A)
        if abs(a) < atol
            A[i] = zero(T)
        elseif abs(a) > 1/atol || isnan(a)
            A[i] = T(Inf)
        end
    end
    return A
end

end


text:
using LiveServer

include("make.jl")

serve(dir="build")

ext: jl


text:
mutable struct Frame
    cl::ClosureObj
    ip::Int
    base_ptr::Int

    Frame(cl::ClosureObj, base_ptr::Int) = new(cl, 0, base_ptr)
end

instructions(f::Frame) = f.cl.fn.instructions
within_loop(f::Frame) = f.cl.fn.within_loop

mutable struct VM
    constants::Vector{Object}
    stack::Vector{Object}
    sp::Int64
    globals::Vector{Object}
    frames::Vector{Frame}
    input::IO
    output::IO

    VM(bc::ByteCode, globals::Vector{Object} = Object[]; input = stdin, output = stdout) =
        begin
            main_fn = CompiledFunctionObj(bc.instructions, 0, 0, false)
            main_closure = ClosureObj(main_fn, [])
            main_frame = Frame(main_closure, 0)
            frames = [main_frame]
            new(bc.constants, [], 1, globals, frames, input, output)
        end
end

Base.push!(vm::VM, obj::Object) = begin
    if vm.sp > length(vm.stack)
        append!(vm.stack, fill(_NULL, vm.sp - 1 - length(vm.stack)))
        push!(vm.stack, obj)
    else
        vm.stack[vm.sp] = obj
    end

    vm.sp += 1
end

push_frame!(vm::VM, frame::Frame) = push!(vm.frames, frame)

push_closure!(vm::VM, id::Int, free_count::Integer) = begin
    fn = vm.constants[id]
    free = vm.stack[vm.sp-free_count:vm.sp-1]
    vm.sp -= free_count
    closure = ClosureObj(fn, free)
    push!(vm, closure)
end

Base.pop!(vm::VM) = begin
    if vm.sp <= 1 || vm.sp > length(vm.stack) + 1
        error("stack underflow")
    end
    vm.sp -= 1
    return vm.stack[vm.sp]
end

pop_frame!(vm::VM) = pop!(vm.frames)

last_popped(vm::VM) = vm.sp > length(vm.stack) ? nothing : vm.stack[vm.sp]

current_frame(vm::VM) = vm.frames[end]

instructions(vm::VM) = instructions(current_frame(vm))

macro monkey_vm_str(code::String)
    quote
        run($(esc(Meta.parse("\"$(escape_string(code))\""))))
    end
end

run(code::String; input = stdin, output = stdout) = begin
    raw_program = parse(code; input, output)
    if !isnothing(raw_program)
        macro_env = Environment(; input, output)
        program = define_macros!(macro_env, raw_program)
        expanded = expand_macros(program, macro_env)

        syntax_check_result = analyze(expanded)
        if isa(syntax_check_result, ErrorObj)
            println(output, syntax_check_result)
            return syntax_check_result
        end

        c = Compiler()
        compile!(c, expanded)
        vm = VM(bytecode(c), Object[]; input, output)

        result = run!(vm)
        if isa(result, ErrorObj)
            println(output, result)
        end

        return result
    end
end

run!(vm::VM) = begin
    @debug "Executing instructions:\n$(string(instructions(vm)))"

    while current_frame(vm).ip < length(instructions(vm))
        current_frame(vm).ip += 1
        ip = current_frame(vm).ip
        ins = instructions(vm)
        op = OpCode(ins[ip])
        if op == OpConstant
            current_frame(vm).ip += 2
            const_id = read_uint16(ins, ip + 1) + 1
            if 1 <= const_id <= length(vm.constants)
                push!(vm, vm.constants[const_id])
            else
                return ErrorObj(
                    "bounds error: attempt to access $(length(vm.constants))-element vector at index [$const_id]",
                )
            end
        elseif OpAdd <= op <= OpGreaterThan
            right = pop!(vm)
            left = pop!(vm)
            err = execute_binary_operation!(vm, op, left, right)
            if isa(err, ErrorObj)
                return err
            end
        elseif OpMinus <= op <= OpBang
            right = pop!(vm)
            err = execute_unary_operation!(vm, op, right)
            if isa(err, ErrorObj)
                return err
            end
        elseif op == OpPop
            pop!(vm)
        elseif op == OpTrue
            push!(vm, _TRUE)
        elseif op == OpFalse
            push!(vm, _FALSE)
        elseif op == OpNull
            push!(vm, _NULL)
        elseif OpJump <= op <= OpJumpNotTruthy
            pos = read_uint16(ins, ip + 1)

            if op == OpJump
                current_frame(vm).ip = pos
            else
                current_frame(vm).ip += 2
                condition = pop!(vm)
                if !is_truthy(condition)
                    current_frame(vm).ip = pos
                end
            end
        elseif OpGetGlobal <= op <= OpSetGlobal
            current_frame(vm).ip += 2
            global_id = read_uint16(ins, ip + 1)

            if op == OpSetGlobal
                if global_id + 1 > length(vm.globals)
                    push!(vm.globals, pop!(vm))
                else
                    vm.globals[global_id+1] = pop!(vm)
                end
            else
                push!(vm, vm.globals[global_id+1])
            end
        elseif OpGetLocal <= op <= OpSetLocal
            current_frame(vm).ip += 1
            local_id = ins[ip+1]
            frame = current_frame(vm)

            if op == OpSetLocal
                vm.stack[frame.base_ptr+local_id] = pop!(vm)
            else
                push!(vm, vm.stack[frame.base_ptr+local_id])
            end
        elseif OpGetOuter <= op <= OpSetOuter
            current_frame(vm).ip += 3
            level = ins[ip+1]
            scope = Scope(ins[ip+2])
            id = ins[ip+3]

            frame = vm.frames[end-level]

            if op == OpGetOuter
                if scope == LocalScope
                    push!(vm, vm.stack[frame.base_ptr+id])
                elseif scope == FreeScope
                    push!(vm, frame.cl.free[id+1])
                elseif scope == FunctionScope
                    push!(vm, frame.cl)
                end
            else
                if scope == LocalScope
                    vm.stack[frame.base_ptr+id] = pop!(vm)
                elseif scope == FreeScope
                    frame.cl.free[id+1] = pop!(vm)
                end
            end
        elseif op == OpGetBuiltin
            current_frame(vm).ip += 1
            builtin_id = ins[ip+1] + 1
            builtin = BUILTINS[builtin_id].second
            push!(vm, builtin)
        elseif OpGetFree <= op <= OpSetFree
            current_frame(vm).ip += 1
            free_id = ins[ip+1] + 1

            if op == OpGetFree
                push!(vm, current_frame(vm).cl.free[free_id])
            else
                current_frame(vm).cl.free[free_id] = pop!(vm)
            end
        elseif op == OpArray
            current_frame(vm).ip += 2
            element_count = read_uint16(ins, ip + 1)
            arr = build_array!(vm, element_count)
            push!(vm, arr)
        elseif op == OpHash
            current_frame(vm).ip += 2
            element_count = read_uint16(ins, ip + 1)
            hash = build_hash!(vm, element_count)
            push!(vm, hash)
        elseif op == OpClosure
            current_frame(vm).ip += 3
            const_id = read_uint16(ins, ip + 1) + 1
            free_count = ins[ip+3]
            push_closure!(vm, const_id, free_count)
        elseif op == OpCurrentClosure
            push!(vm, current_frame(vm).cl)
        elseif op == OpIndex
            index = pop!(vm)
            left = pop!(vm)
            if isa(left, ArrayObj)
                if !isa(index, IntegerObj)
                    return ErrorObj("unsupported index type: $(type_of(index))")
                end

                i = index.value
                if 0 <= i < length(left.elements)
                    push!(vm, left.elements[i+1])
                else
                    push!(vm, _NULL)
                end
            elseif isa(left, HashObj)
                if index ∈ keys(left.pairs)
                    push!(vm, left.pairs[index])
                else
                    push!(vm, _NULL)
                end
            else
                return ErrorObj("index operator not supported: $(type_of(left))")
            end
        elseif op == OpCall
            current_frame(vm).ip += 1
            arg_count = ins[ip+1]
            callee = vm.stack[vm.sp-1-arg_count]
            err = call!(vm, callee, arg_count)
            if isa(err, ErrorObj)
                return err
            end
        elseif OpBreak <= op <= OpContinue
            current_frame(vm).ip += 1
            return_value = native_bool_to_boolean_obj(op == OpContinue)
            frame = pop_frame!(vm)
            vm.sp = frame.base_ptr - 1
            push!(vm, return_value)
        elseif OpReturnValue <= op <= OpReturn
            return_value = op == OpReturnValue ? pop!(vm) : _NULL
            frame = pop_frame!(vm)
            if op == OpReturnValue
                while within_loop(frame)
                    frame = pop_frame!(vm)
                end
            end
            if isempty(vm.frames)
                return return_value
            end
            vm.sp = frame.base_ptr - 1
            push!(vm, return_value)
        end
    end

    return last_popped(vm)
end

execute_unary_operation!(vm::VM, op::OpCode, right::Object) = begin
    if op == OpBang
        push!(vm, native_bool_to_boolean_obj(!is_truthy(right)))
    elseif op == OpMinus
        if isa(right, IntegerObj)
            push!(vm, IntegerObj(-right.value))
        else
            return ErrorObj("unknown operator: -$(type_of(right))")
        end
    end
end

execute_binary_operation!(vm::VM, op::OpCode, left::Object, right::Object) = begin
    if type_of(left) != type_of(right)
        return ErrorObj(
            "type mismatch: " *
            type_of(left) *
            " " *
            OPCODE_STRINGS[op] *
            " " *
            type_of(right),
        )
    end

    result = if op == OpEqual
        native_bool_to_boolean_obj(left == right)
    elseif op == OpNotEqual
        native_bool_to_boolean_obj(left != right)
    else
        return ErrorObj(
            "unknown operator: " *
            type_of(left) *
            " " *
            OPCODE_STRINGS[op] *
            " " *
            type_of(right),
        )
    end

    push!(vm, result)
end

execute_binary_operation!(vm::VM, op::OpCode, left::StringObj, right::StringObj) = begin
    lval = left.value
    rval = right.value

    result = if op == OpAdd
        StringObj(lval * rval)
    elseif op == OpEqual
        native_bool_to_boolean_obj(lval == rval)
    elseif op == OpNotEqual
        native_bool_to_boolean_obj(lval != rval)
    else
        return ErrorObj("unknown operator: STRING " * OPCODE_STRINGS[op] * " STRING")
    end

    push!(vm, result)
end

execute_binary_operation!(vm::VM, op::OpCode, left::IntegerObj, right::IntegerObj) = begin
    lval = left.value
    rval = right.value

    result = if op == OpAdd
        IntegerObj(lval + rval)
    elseif op == OpSub
        IntegerObj(lval - rval)
    elseif op == OpMul
        IntegerObj(lval * rval)
    elseif op == OpDiv
        if rval == 0
            return ErrorObj("divide error: division by zero")
        end
        IntegerObj(lval ÷ rval)
    elseif op == OpEqual
        native_bool_to_boolean_obj(lval == rval)
    elseif op == OpNotEqual
        native_bool_to_boolean_obj(lval != rval)
    elseif op == OpLessThan
        native_bool_to_boolean_obj(lval < rval)
    elseif op == OpGreaterThan
        native_bool_to_boolean_obj(lval > rval)
    else
        return ErrorObj("unknown operator: INTEGER " * string(op) * " INTEGER")
    end

    push!(vm, result)
end

build_array!(vm::VM, element_count::Integer)::ArrayObj = begin
    elements = vm.stack[vm.sp-element_count:vm.sp-1]

    vm.sp -= element_count

    return ArrayObj(elements)
end

build_hash!(vm::VM, element_count::Integer)::HashObj = begin
    elements = vm.stack[vm.sp-element_count:vm.sp-1]
    vm.sp -= element_count
    prs = Dict(elements[i] => elements[i+1] for i = 1:2:length(elements))
    return HashObj(prs)
end

call!(::VM, obj::Object, ::Integer) = return ErrorObj("not a function: $(type_of(obj))")

call!(vm::VM, cl::ClosureObj, arg_count::Integer) = begin
    if arg_count != cl.fn.param_count
        vm.sp += cl.fn.local_count - arg_count
        return ErrorObj("argument error: wrong number of arguments: got $(arg_count)")
    else
        frame = Frame(cl, vm.sp - arg_count)
        push_frame!(vm, frame)
        vm.sp += cl.fn.local_count - arg_count
    end
end

call!(vm::VM, builtin::Builtin, arg_count::Integer) = begin
    args = vm.stack[vm.sp-arg_count:vm.sp-1]
    result =
        builtin.fn(args...; env = Environment(; input = vm.input, output = vm.output))
    if isa(result, ErrorObj)
        return result
    end
    vm.sp -= arg_count + 1
    push!(vm, result)
end

native_bool_to_boolean_obj(b::Bool) = b ? _TRUE : _FALSE


text:
using Formatting
numbers = map(x -> parse(Int64,x),split(first(readlines(open("input"))), ","))
fuel_cost(pivot, numbers, distance) = sum(map(x -> distance(pivot, x), numbers))
minimum_fuel_cost(numbers, metric) = format(minimum(map(x -> fuel_cost(x, numbers, metric), 0:maximum(numbers))))
minimum_fuel_cost(numbers, (x, y) -> abs(x-y))
minimum_fuel_cost(numbers, (x, y) -> (abs(x-y)*(abs(x-y)+1)) / 2) # small gauss

max_stars_repo_name: Ximsa/aoc

max_issues_repo_name: Ximsa/aoc

max_forks_repo_name: Ximsa/aoc


text:
using DifferentialEquations
using SimpleDiffEq
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

function infection_rate(u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    β*c*I/N*S
end
function infection!(integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
infection_jump = ConstantRateJump(infection_rate,infection!);

function recovery_rate(u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    γ*I
end
function recovery!(integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
recovery_jump = ConstantRateJump(recovery_rate,recovery!);

tmax = 40.0
tspan = (0.0,tmax);

δt = 0.1
t = 0:δt:tmax;

u0 = [990,10,0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob = DiscreteProblem(u0,tspan,p)

prob_jump = JumpProblem(prob,Direct(),infection_jump,recovery_jump)

sol_jump = solve(prob_jump,SSAStepper());

out_jump = sol_jump(t);

df_jump = DataFrame(out_jump')
df_jump[!,:t] = out_jump.t;

@df df_jump plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_jump,FunctionMap())

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
"""
    SweepState

State of a DMRG sweep.
"""
mutable struct SweepState{L,T}
    "Wavefunction."
    psi::MPS{L,T}
    "Hamiltonian."
    H::MPO{L,T}

    "Sweep counter."
    sweep_num::Int
    "Site index."
    site::Int
    "Sweep direction."
    dir::Direction

    "Hamiltonian contractions."
    H_cntrctns::OffsetArray{Contraction{T,1},1}
    "Hamiltonian fluctuation contraction."
    H2_cntrctn::Contraction{T,2}

    # SweepDetails parameters.
    energy::Float64
    dH2::Float64
    middle_eigvals::Vector{Float64}
    mps_max_bond_dim::Int
    duration::Float64
end

"""
    SweepState(psi::MPS{L,T}, H::MPO{L,T})

Generate a blank `SweepState` for wavefunction `psi` and Hamiltonian `H`.
"""
function SweepState(psi::MPS{L,T}, H::MPO{L,T}) where {L,T}
    # Initialize the H contractions for the first sweep.
    H_cntrctns = OffsetArray{Contraction{T,1}}(undef, 0:(L+1))
    H_cntrctns[L+1] = cap_contraction(T, Left, L, 1)
    for i in L:-1:3
        H_cntrctns[i] = H_cntrctns[i+1] * (psi, H)
    end
    H_cntrctns[2] = dummy_contraction(T, 1)
    H_cntrctns[1] = dummy_contraction(T, 1)
    H_cntrctns[0] = cap_contraction(T, Right, L, 1)

    # Use dummy values for parameters that we don't know yet.
    SweepState(psi, H, -1, -1, Left, H_cntrctns, dummy_contraction(T, 2),
               NaN, NaN, Float64[], -1, NaN)
end

Base.broadcastable(x::SweepState) = Ref(x)


text:
include("JuliaOpenGL.jl")
main()

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


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

mutable struct IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus <: SwaggerModel
    currentCPUUtilizationPercentage::Any # spec type: Union{ Nothing, Int32 } # spec name: currentCPUUtilizationPercentage
    currentReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: currentReplicas
    desiredReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: desiredReplicas
    lastScaleTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastScaleTime
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration

    function IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus(;currentCPUUtilizationPercentage=nothing, currentReplicas=nothing, desiredReplicas=nothing, lastScaleTime=nothing, observedGeneration=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("currentCPUUtilizationPercentage"), currentCPUUtilizationPercentage)
        setfield!(o, Symbol("currentCPUUtilizationPercentage"), currentCPUUtilizationPercentage)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("currentReplicas"), currentReplicas)
        setfield!(o, Symbol("currentReplicas"), currentReplicas)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("desiredReplicas"), desiredReplicas)
        setfield!(o, Symbol("desiredReplicas"), desiredReplicas)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("lastScaleTime"), lastScaleTime)
        setfield!(o, Symbol("lastScaleTime"), lastScaleTime)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        o
    end
end # type IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus

const _property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus = Dict{Symbol,Symbol}(Symbol("currentCPUUtilizationPercentage")=>Symbol("currentCPUUtilizationPercentage"), Symbol("currentReplicas")=>Symbol("currentReplicas"), Symbol("desiredReplicas")=>Symbol("desiredReplicas"), Symbol("lastScaleTime")=>Symbol("lastScaleTime"), Symbol("observedGeneration")=>Symbol("observedGeneration"))
const _property_types_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus = Dict{Symbol,String}(Symbol("currentCPUUtilizationPercentage")=>"Int32", Symbol("currentReplicas")=>"Int32", Symbol("desiredReplicas")=>"Int32", Symbol("lastScaleTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("observedGeneration")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus)
    (getproperty(o, Symbol("currentReplicas")) === nothing) && (return false)
    (getproperty(o, Symbol("desiredReplicas")) === nothing) && (return false)
    true
end

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


text:
#src Introduction to DG methods
using Test: @test #src

# This tutorial is about how to set up a simple way to approximate the solution of a hyperbolic partial
# differential equation. First, we will implement a basic and naive algorithm. Then, we will use predefined
# features from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) to show how you can use Trixi.jl on your own.

# We will implement the scalar linear advection equation in 1D with the advection velocity $1$.
# ```math
# u_t + u_x = 0,\; \text{for} \;t\in \mathbb{R}^+, x\in\Omega=[-1,1]
# ```
# We define the domain $\Omega$ by setting the boundaries.
coordinates_min = -1.0 # minimum coordinate
coordinates_max = 1.0  # maximum coordinate

# We assume periodic boundaries and the following initial condition.
initial_condition_sine_wave(x) = 1.0 + 0.5 * sin(pi * x)

# ## The discontinuous Galerkin collocation spectral element method (DGSEM)
# ### i. Discretization of the physical domain
# To improve precision we want to approximate the solution on small parts of the physical domain.
# So, we split the domain $\Omega=[-1, 1]$ into elements $Q_l$ of length $dx$.

n_elements = 16 # number of elements

dx = (coordinates_max - coordinates_min) / n_elements # length of one element

# To make the calculation more efficient and storing less information, we transform each element
# $Q_l$ with center point $x_l$ to a reference element $E=[-1, 1]$
# ```math
# Q_l=\Big[x_l-\frac{dx}{2}, x_l+\frac{dx}{2}\Big] \underset{x(\xi)}{\overset{\xi(x)}{\rightleftarrows}} [-1, 1].
# ```
# So, for every element the transformation from the reference domain to the physical domain is defined by
# ```math
# x(\xi) = x_l + \frac{dx}{2} \xi,\; \xi\in[-1, 1]
# ```
# Therefore,
# ```math
# \begin{align*}
# u &= u(x(\xi), t) \\
# u_x &= u_\xi \frac{d\xi}{dx} \\[3pt]
# \frac{d\xi}{dx} &= (x_\xi)^{-1} = \frac{2}{dx} =: J^{-1}. \\
# \end{align*}
# ```
# Here, $J$ is the Jacobian determinant of the transformation.

# Using this transformation, we can transform our equation for all elements $Q_l$.
# ```math
# \frac{dx}{2} u_t^{Q_l} + u_\xi^{Q_l} = 0 \text{, for }t\in\mathbb{R}^+,\; \xi\in[-1, 1]
# ```

# ### ii. Polynomial approach
# Now, we want to approximate the solution in each element $Q_l$ by a polynomial of degree $N$. Since we transformed
# the equation, we can use the same polynomial approach for the reference coordinate $\xi\in[-1, 1]$ in every
# physical element $Q_l$. This saves a lot of ressources by reducing the amount of calculations needed
# and storing less information.

# For DGSEM we choose [Lagrange basis functions](https://en.wikipedia.org/wiki/Lagrange_polynomial)
# $\{l_j\}_{j=0}^N$ as our polynomial basis of degree $N$ in $[-1, 1]$.
# The solution in element $Q_l$ can be approximated by
# ```math
# u(x(\xi), t)\big|_{Q_l} \approx u^{Q_l}(\xi, t) = \sum_{j=0}^N u_j^{Q_l}(t) l_j(\xi)
# ```
# with $N+1$ coefficients $\{u_j^{Q_l}\}_{j=0}^N$.
# By contruction the Lagrange basis has some useful advantages. This basis is defined by $N+1$ nodes, which
# fulfill a Kronecker property at the exact same nodes. Let $\{\xi_i\}_{i=0}^N$ be these nodes.
# ```math
# l_j(\xi_i) = \delta_{i,j} =
# \begin{cases}
# 1, & \text{if } i=j \\
# 0, & \text{else.}
# \end{cases}
# ```
# Because of this property, the polynomial coefficients are exact the values of $u^{Q_l}$ at the nodes
# ```math
# u^{Q_l}(\xi_i, t) = \sum_{j=0}^N u_j^{Q_l}(t) \underbrace{l_j(\xi_i)}_{=\delta_{ij}} = u_i^{Q_l}(t).
# ```

# Next, we want to select the nodes $\{\xi_i\}_{i=0}^N$, which we use for the construction of the Lagrange
# polynomials. We choose the $N+1$ Gauss-Lobatto nodes, which are used for the
# [Gaussian-Lobatto quadrature](https://mathworld.wolfram.com/LobattoQuadrature.html).
# These always contain the boundary points at $-1$ and $+1$ and are well suited as interpolation nodes.
# The corresponding weights will be referred to as $\{w_j\}_{j=0}^N$.
# In Trixi.jl the basis with Lagrange polynomials on Gauss-Lobatto nodes is already defined.
using Trixi
polydeg = 3 #= polynomial degree = N =#
basis = LobattoLegendreBasis(polydeg)
# The Gauss-Lobatto nodes are
nodes = basis.nodes
# with the corresponding weights
weights = basis.weights

# To illustrate how you can integrate using numerical quadrature with this Legendre-Gauss-Lobatto nodes,
# we give an example for $f(x)=x^3$. Since $f$ is of degree $3$, a polynomial interpolation with $N=3$ is exact.
# Therefore, the integral on $[-1, 1]$ can be calculated by
# ```math
# \begin{align*}
# \int_{-1}^1 f(x) dx &= \int_{-1}^1 \Big( \sum_{j=0}^3 f(\xi_j)l_j(x) \Big) dx
# = \sum_{j=0}^3 f(\xi_j) \int_{-1}^1 l_j(x)dx \\
# &=: \sum_{j=0}^3 f(\xi_j) w_j
# = \sum_{j=0}^3 \xi_j^3 w_j
# \end{align*}
# ```
# Let's use our nodes and weights for $N=3$ and plug in
integral = sum(nodes.^3 .* weights)

# Using this polynomial approach leads to the equation
# ```math
# \frac{dx}{2} \dot{u}^{Q_l}(\xi, t) + u^{Q_l}(\xi, t)' = 0
# ```
# with $\dot{u}=\frac{\partial}{\partial t}u$ and $u'=\frac{\partial}{\partial x}u$.
# To approximate the solution, we need to get the polynomial coefficients $\{u_j^{Q_l}\}_{j=0}^N$
# for every element $Q_l$.

# After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0$
# for every node.
x = Matrix{Float64}(undef, length(nodes), n_elements)
for element in 1:n_elements
    x_l = coordinates_min + (element - 1) * dx + dx/2
    for i in 1:length(nodes)
        ξ = nodes[i] # nodes in [-1, 1]
        x[i, element] = x_l + dx/2 * ξ
    end
end

u0 = initial_condition_sine_wave.(x)

# To have a look at the initial sinus curve, we plot it.
using Plots
plot(vec(x), vec(u0), label="initial condition", legend=:topleft)

# ### iii. Variational formulation
# After defining the equation and initial condition, we want to implement an algorithm to
# approximate the solution.

# From now on, we only write $u$ instead of $u^{Q_l}$ for simplicity, but consider that all the following
# calculation only concern one element.
# Multiplying the new equation with the smooth Lagrange polynomials $\{l_i\}_{i=0}^N$ (test functions)
# and integrating over the reference element $E=[-1,1]$, we get the variational formulation of our
# transformed partial differential equation for $i=0,...,N$:
# ```math
# \begin{align*}
# \int_{-1}^1 \Big( \frac{dx}{2} \dot{u}(\xi, t) + u'(\xi, t) \Big) l_i(\xi)d\xi
#   &= \underbrace{\frac{dx}{2} \int_{-1}^1 \dot{u}(\xi, t) l_i(\xi)d\xi}_{\text{Term I}} + \underbrace{\int_{-1}^1 u'(\xi, t) l_i(\xi)d\xi}_{\text{Term II}} = 0
# \end{align*}
# ```

# We deal with the two terms separately. We write $\int_{-1, N}^1 \;\cdot\; d\xi$ for the approximation
# of the integral using numerical quadrature with $N+1$ basis points. We use the Gauss-Lobatto nodes
# again. The numerical scalar product $\langle\cdot, \cdot\rangle_N$ is defined by
# $\langle f, g\rangle_N := \int_{-1, N}^1 f(\xi) g(\xi) d\xi$.

# #### Term I:
# In the following calculation we approximate the integral numerically with quadrature on the Gauss-Lobatto
# nodes $\{\xi_i\}_{i=0}^N$ and then use the Kronecker property of the Lagrange polynomials. This approach
# of using the same nodes for the interpolation and quadrature is called collocation.
# ```math
# \begin{align*}
# \frac{dx}{2} \int_{-1}^1 \dot{u}(\xi, t) l_i(\xi)d\xi
# &\approx \frac{dx}{2} \int_{-1, N}^1 \dot{u}(\xi, t) l_i(\xi)d\xi \\
# &= \frac{dx}{2} \sum_{k=0}^N \underbrace{\dot{u}(\xi_k, t)}_{=\dot{u}_k(t)} \underbrace{l_i(\xi_k)}_{=\delta_{k,i}}w_k \\
# &= \frac{dx}{2} \dot{u}_i(t) w_i
# \end{align*}
# ```
# We define the Legendre-Gauss-Lobatto (LGL) mass matrix $M$ and by the Kronecker property follows:
# ```math
# M_{ij} = \langle l_j, l_i\rangle_N = \delta_{ij} w_j,\; i,j=0,...,N.
# ```
using LinearAlgebra
M = diagm(weights)
# Now, we can write the integral with this new matrix.
# ```math
# \frac{dx}{2} \int_{-1, N}^1 \dot{u}(\xi, t) \underline{l}(\xi)d\xi = \frac{dx}{2} M \underline{\dot{u}}(t),
# ```
# where $\underline{\dot{u}} = (\dot{u}_0, ..., \dot{u}_N)^T$ and $\underline{l}$ respectively.

# **Note:** Since the LGL quadrature with $N+1$ nodes is exact up to functions of degree $2N-1$ and
# $\dot{u}(\xi, t) l_i(\xi)$ is of degree $2N$, in general the following holds
# ```math
# \int_{-1}^1 \dot{u}(\xi, t) l_i(\xi) d\xi \neq \int_{-1, N}^1 \dot{u}(\xi, t) l_i(\xi) d\xi.
# ```
# With an exact integration the mass matrix would be dense. Choosing numerical integrating and quadrature
# with the exact same nodes (collocation) leads to the sparse and diagonal mass matrix $M$. This
# is called mass lumping and has the big advantage of an easy invertation of the matrix.

# #### Term II:
# We use spatial partial intergration for the second term:
# ```math
# \int_{-1}^1 u'(\xi, t) l_i(\xi) d\xi = [u l_i]_{-1}^1 - \int_{-1}^1 u l_i'd\xi
# ```
# The resulting integral can be solved exactly with LGL quadrature since the polynomial is now
# of degree $2N-1$.

# Again, we split the calculation in two steps.

# #### Surface term
# As mentioned before, we approximate the solution with a polynomial in every element. Therefore, in
# general the value of this approximation at the interfaces between two elements is not unique. To solve
# this problem we introduce the idea of the numerical flux $u^*$, which will give an exact value at
# the interfaces. One of many different approaches and definitions for the calculation of the
# numerical flux we will deal with in [4. Numerical flux](@ref numerical_flux).
# ```math
# [u l_i]_{-1}^1 = u^*\big|^1 l_i(+1) - u^*\big|_{-1} l_i(-1)
# ```
# Since the Gauss-Lobatto nodes contain the element boundaries $-1$ and $+1$, we can use the
# Kronecker property of $l_i$ for the calculation of $l_i(-1)$ and $l_i(+1)$.
# ```math
# [u \underline{l}]_{-1}^1 = u^*\big|^1 \left(\begin{array}{c} 0 \\ \vdots \\ 0 \\ 1 \end{array}\right)
# - u^*\big|_{-1} \left(\begin{array}{c} 1 \\ 0 \\ \vdots \\ 0\end{array}\right)
# = B \underline{u}^*(t)
# ```
# with the boundary matrix
# ```math
# B = \begin{pmatrix}
# -1 & 0 & \cdots & 0\\
# 0 & 0 & \cdots & 0\\
# \vdots & \vdots & 0 & 0\\
# 0 & \cdots & 0 & 1
# \end{pmatrix}
# \qquad\text{and}\qquad
# \underline{u}^*(t) = \left(\begin{array}{c} u^*\big|_{-1} \\ 0 \\ \vdots \\ 0 \\ u^*\big|^1\end{array}\right).
# ```
B = diagm([-1; zeros(polydeg - 1); 1])

# #### Volume term
# As mentioned before, the new integral can be solved exact since the function inside is of degree $2N-1$.
# ```math
# - \int_{-1}^1 u l_i'd\xi = - \int_{-1, N}^1 u l_i' d\xi
# = - \sum_{k=0}^N u(\xi_k, t) l_i'(\xi_k) w_k
# = - \sum_{k=0}^N u_k(t) D_{ki} w_k
# ```
# where $D$ is the derivative matrix defined by $D_{ki} = l_i'(\xi_k)$ for $i,k=0,...,N$.
D = basis.derivative_matrix

# To show why this matrix is called the derivative matrix, we go back to our example $f(x)=x^3$.
# We calculate the derivation of $f$ at the Gauss-Lobatto nodes $\{\xi_k\}_{k=0}^N$ with $N=8$.
# ```math
# f'|_{x=\xi_k} = \Big( \sum_{j=0}^8 f(\xi_j) l_j(x) \Big)'|_{x=\xi_k} = \sum_{j=0}^8 f(\xi_j) l_j'(\xi_k)
# = \sum_{j=0}^8 f(\xi_j) D_{kj}
# ```
# for $k=0,...,N$ and therefore, $\underline{f}' = D \underline{f}$.
basis_N8 = LobattoLegendreBasis(8)
plot(vec(x), x -> 3 * x^2, label="f'", lw=2)
scatter!(basis_N8.nodes, basis_N8.derivative_matrix * basis_N8.nodes.^3, label="Df", lw=3)

# Combining the volume term for every $i=0,...,N$ results in
# ```math
# \int_{-1}^1 u \underline{l'} d\xi = - D^T M \underline{u}(t)
# ```

# Putting all parts together we get the following equation for the element $Q_l$
# ```math
# \frac{dx}{2} M \underline{\dot{u}}(t) = - B \underline{u}^*(t) + D^T M \underline{u}(t)
# ```
# or equivalent
# ```math
# \underline{\dot{u}}^{Q_l}(t) = \frac{2}{dx} \Big[ - M^{-1} B \underline{u}^{{Q_l}^*}(t) + M^{-1} D^T M \underline{u}^{Q_l}(t)\Big].
# ```
# This is called the weak form of the DGSEM.

# **Note:** For every element $Q_l$ we get a system of $N+1$ ordinary differential equations to
# calculate $N+1$ coefficients. Since the numerical flux $u^*$ is depending on extern values at
# the interfaces, the equation systems of adjacent elements are weakly linked.

# ### [iv. Numerical flux](@id numerical_flux)
# As mentioned above, we still have to handle the problem of different values at the same point at
# the interfaces. This happens with the ideas of the numerical flux $f^*(u)=u^*$. The role of $f^*$
# might seem minor in this simple example, but is important for more complicated problems.
# There are two values at the same spatial coordinate. Let's say we are looking at the interface between
# the elements $Q_l$ and $Q_{l+1}$, while both elements got $N+1$ nodes as defined before. We call
# the first value of the right element $u_R=u_0^{Q_{l+1}}$ and the last one of the left element
# $u_L=u_N^{Q_l}$. So, for the value of the numerical flux on that interface the following holds
# ```math
# u^* = u^*(u_L, u_R).
# ```
# These values are interpreted as start values of a so-called Riemann problem. There are many
# different (approximate) Riemann solvers available and useful for different problems. We will
# use the local Lax-Friedrichs flux.
surface_flux = flux_lax_friedrichs

# The only missing ingredient is the flux calculation at the boundaries $-1$ and $+1$.
# ```math
# u^{{Q_{first}}^*}\big|_{-1} = u^{{Q_{first}}^*}\big|_{-1}(u^{bound}(-1), u_R)
# \quad\text{and}\quad
# u^{{Q_{last}}^*}\big|^1 = u^{{Q_{last}}^*}\big|^1(u_L, u^{bound}(1))
# ```
# The boundaries are periodic, which means that the last value of the last element $u^{Q_{last}}_N$
# is used as $u_L$ at the first interface and accordingly for the other boundary.

# Now, we implement a function, that calculates $\underline{\dot{u}}^{Q_l}$ for the given matrices,
# $\underline{u}$ and $\underline{u}^*$.
function rhs!(du, u, x, t)
    ## Reset du and flux matrix
    du .= zero(eltype(du))
    flux_numerical = copy(du)

    ## Calculate interface and boundary fluxes, $u^* = (u^*|_{-1}, 0, ..., 0, u^*|^1)^T$
    ## Since we use the flux Lax-Friedrichs from Trixi.jl, we have to pass some extra arguments.
    ## Trixi needs the equation we are dealing with and an additional `1`, that indicates the
    ## first coordinate direction.
    equations = LinearScalarAdvectionEquation1D(1.0)
    for element in 2:n_elements-1
        ## left interface
        flux_numerical[1, element] = surface_flux(u[end, element-1], u[1, element], 1, equations)
        flux_numerical[end, element-1] = flux_numerical[1, element]
        ## right interface
        flux_numerical[end, element] = surface_flux(u[end, element], u[1, element+1], 1, equations)
        flux_numerical[1, element+1] = flux_numerical[end, element]
    end
    ## boundary flux
    flux_numerical[1, 1] = surface_flux(u[end, end], u[1, 1], 1, equations)
    flux_numerical[end, end] = flux_numerical[1, 1]

    ## Calculate surface integrals, $- M^{-1} * B * u^*$
    for element in 1:n_elements
        du[:, element] -= (M \ B) * flux_numerical[:, element]
    end

    ## Calculate volume integral, $+ M^{-1} * D^T * M * u$
    for element in 1:n_elements
        flux = u[:, element]
        du[:, element] += (M \ transpose(D)) * M * flux
    end

    ## Apply Jacobian from mapping to reference element
    for element in 1:n_elements
        du[:, element] *= 2 / dx
    end

    return nothing
end

# Combining all definitions and the function that calculates the right-hand side, we define the ODE and
# solve it until `t=2` with OrdinaryDiffEq's `solve` function and the Runge-Kutta method `RDPK3SpFSAL49()`,
# which is optimized for discontinuous Galerkin methods and hyperbolic PDEs. We set some common
# error tolerances `abstol=1.0e-6, reltol=1.0e-6` and pass `save_everystep=false` to avoid saving intermediate
# solution vectors in memory.
using OrdinaryDiffEq
tspan = (0.0, 2.0)
ode = ODEProblem(rhs!, u0, tspan, x)

sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false)
@test maximum(abs.(u0 - sol.u[end])) < 5e-5 #src

plot(vec(x), vec(sol.u[end]), label="solution at t=$(tspan[2])", legend=:topleft, lw=3)

# ## Alternative Implementation based on Trixi.jl
# Now, we implement the same example. But this time, we directly use the functionality that Trixi.jl
# provides.

using Trixi, OrdinaryDiffEq, Plots

# First, define the equation with a advection_velocity of `1`.
advection_velocity = 1.0
equations = LinearScalarAdvectionEquation1D(advection_velocity)

# Then, create a DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux.
# The implementation of the basis and the numerical flux is now already done.
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# We will now create a mesh with 16 elements for the physical domain `[-1, 1]` with periodic boundaries.
# We use Trixi.jl's standard mesh [`TreeMesh`](@ref). Since it's limited to hypercube domains, we
# choose `2^4=16` elements. The mesh type supports AMR, that' why `n_cells_max` has to be set, even
# if we don't need AMR here.
coordinates_min = -1.0 # minimum coordinate
coordinates_max = 1.0  # maximum coordinate
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4, # number of elements = 2^4
                n_cells_max=30_000) # set maximum capacity of tree data structure (only needed for AMR)

# A semidiscretization collects data structures and functions for the spatial discretization.
# In Trixi, an initial condition has the following parameter structure and is of the type `SVector`.
initial_condition_sine_wave(x, t, equations) = SVector(1.0 + 0.5 * sin(pi * sum(x - equations.advection_velocity * t)))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sine_wave, solver)

# Again, combining all definitions and the function that calculates the right-hand side, we define the ODE and
# solve it until `t=2` with OrdinaryDiffEq's `solve` function and the Runge-Kutta method `RDPK3SpFSAL49()`.
tspan = (0.0, 2.0)
ode_trixi  = semidiscretize(semi, tspan)

sol_trixi = solve(ode_trixi, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false);

# We add a plot of the new approximated solution to the one calculated before.
plot!(sol_trixi, label="solution at t=$(tspan[2]) with Trixi.jl", legend=:topleft, linestyle=:dash, lw=2)
@test maximum(abs.(vec(u0) - sol_trixi.u[end])) ≈ maximum(abs.(u0 - sol.u[end])) #src

# ## Summary of the code
# To sum up, here is the complete code that we used.

# ### Raw implementation
## basis: Legendre-Gauss-Lobatto
using Trixi, LinearAlgebra, OrdinaryDiffEq, Plots
polydeg = 3 #= polynomial degree =#
basis = LobattoLegendreBasis(polydeg)
nodes = basis.nodes # Gauss-Lobatto nodes in [-1, 1]
D = basis.derivative_matrix
M = diagm(basis.weights) # mass matrix
B = diagm([-1; zeros(polydeg - 1); 1])

## mesh
coordinates_min = -1.0 # minimum coordinate
coordinates_max = 1.0  # maximum coordinate
n_elements      = 16   # number of elements

dx = (coordinates_max - coordinates_min) / n_elements # length of one element

x = Matrix{Float64}(undef, length(nodes), n_elements)
for element in 1:n_elements
    x_l = -1 + (element - 1) * dx + dx/2
    for i in 1:length(nodes) # basis points in [-1, 1]
        ξ = nodes[i]
        x[i, element] = x_l + dx/2 * ξ
    end
end

## initial condition
initial_condition_sine_wave(x) = 1.0 + 0.5 * sin(pi * x)
u0 = initial_condition_sine_wave.(x)

plot(vec(x), vec(u0), label="initial condition", legend=:topleft)

## flux Lax-Friedrichs
surface_flux = flux_lax_friedrichs

## rhs! method
function rhs!(du, u, x, t)
    ## reset du
    du .= zero(eltype(du))
    flux_numerical = copy(du)

    ## calculate interface and boundary fluxes
    equations = LinearScalarAdvectionEquation1D(1.0)
    for element in 2:n_elements-1
        ## left interface
        flux_numerical[1, element] = surface_flux(u[end, element-1], u[1, element], 1, equations)
        flux_numerical[end, element-1] = flux_numerical[1, element]
        ## right interface
        flux_numerical[end, element] = surface_flux(u[end, element], u[1, element+1], 1, equations)
        flux_numerical[1, element+1] = flux_numerical[end, element]
    end
    ## boundary flux
    flux_numerical[1, 1] = surface_flux(u[end, end], u[1, 1], 1, equations)
    flux_numerical[end, end] = flux_numerical[1, 1]

    ## calculate surface integrals
    for element in 1:n_elements
        du[:, element] -= (M \ B) * flux_numerical[:, element]
    end

    ## calculate volume integral
    for element in 1:n_elements
        flux = u[:, element]
        du[:, element] += (M \ transpose(D)) * M * flux
    end

    ## apply Jacobian from mapping to reference element
    for element in 1:n_elements
        du[:, element] *= 2 / dx
    end

    return nothing
end

## create ODE problem
tspan = (0.0, 2.0)
ode = ODEProblem(rhs!, u0, tspan, x)

## solve
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false)
@test maximum(abs.(vec(u0) - sol_trixi.u[end])) ≈ maximum(abs.(u0 - sol.u[end])) #src

plot(vec(x), vec(sol.u[end]), label="solution at t=$(tspan[2])", legend=:topleft, lw=3)

# ### Alternative Implementation based on Trixi.jl
using Trixi, OrdinaryDiffEq, Plots

## equation with a advection_velocity of `1`.
advection_velocity = 1.0
equations = LinearScalarAdvectionEquation1D(advection_velocity)

## create DG solver with flux lax friedrichs and LGL basis
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

## distretize domain with `TreeMesh`
coordinates_min = -1.0 # minimum coordinate
coordinates_max = 1.0 # maximum coordinate
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4, # number of elements = 2^4
                n_cells_max=30_000)

## create initial condition and semidiscretization
initial_condition_sine_wave(x, t, equations) = SVector(1.0 + 0.5 * sin(pi * sum(x - equations.advection_velocity * t)))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sine_wave, solver)

## solve
tspan = (0.0, 2.0)
ode_trixi  = semidiscretize(semi, tspan)
sol_trixi  = solve(ode_trixi, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false);

plot!(sol_trixi, label="solution at t=$(tspan[2]) with Trixi.jl", legend=:topleft, linestyle=:dash, lw=2)
@test maximum(abs.(vec(u0) - sol_trixi.u[end])) ≈ maximum(abs.(u0 - sol.u[end])) #src


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.06"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
n, N = 30, 10^6

# ╔═╡ d599e1d4-3a99-11eb-1ab2-637caef1df0d
begin
	dist1 = Uniform(1-sqrt(3),1+sqrt(3))
	dist2 = Exponential(1)
	dist3 = Normal(1,1)

	data1 = [mean(rand(dist1,n)) for _ in 1:N]
	data2 = [mean(rand(dist2,n)) for _ in 1:N]
	data3 = [mean(rand(dist3,n)) for _ in 1:N]
end

# ╔═╡ d59a153c-3a99-11eb-1cfe-f54834ffbc5c
begin
	stephist([data1 data2 data3], bins=100,
		c=[:blue :red :green], xlabel = "x", ylabel = "Density",
		label=["Average of Uniforms" "Average of Exponentials" "Average of Normals"],
		normed=true, xlims=(0,2), ylims=(0,2.5))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.06"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═d599e1d4-3a99-11eb-1ab2-637caef1df0d
# ╠═d59a153c-3a99-11eb-1cfe-f54834ffbc5c
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text: # Add queue definitions here


text:
@testset "Population dynamics" begin
    @test sum(P1) ≈ K
    @test all(vecstate(P1) .== vecstate(eqdist_ap))
    @test all(vecstate(Phalf1) .≥ vecstate(Phalf))
    @test sum(Phalf1) ≈ 51.25

    @test P1 isa PopState
    @test P1_pt isa PopState
    @test P1_uspt isa PopState
    @test P1_mspt isa PopState
end


text:
# Convenient include for the common functions
include("run_sim!.jl")
include("one_step!.jl")
include("approach_wall!.jl")
include("boundary_check.jl")
include("boundary_cross.jl")
include("evaluate_proposed.jl")

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
module DependencyTrees

export
    DependencyTree,
    deptree, deptoken,
    is_projective,
    Treebank,
    from_conllu,

    typed, untyped,

    labeled_accuracy, unlabeled_accuracy

import Base.==
import Base.getindex
import Base.reduce
import Base: iterate, IteratorSize, length

include("errors.jl")
include("tokens.jl")
include("trees.jl")
include("conllu.jl")

include("treebanks.jl")

include("transition_parsing/TransitionParsing.jl")
using .TransitionParsing

include("evaluation.jl")

end # module


text:
# EXPERIMENTAL, not to be used!
# To use, uncomment the code below and add NearestNeighbors.

#=

# ngram NN-model for faster approximate string matching
# So far it is not really workable, results are bad
# An example can be found in "./scripts/test_word_model.jl"
#TODO(Corneliu) Improve quality and performance of this model

# Build model
function build_nn_model(words; ngram_size::Int=2)
    #words = collect(keys(cptnet))
    # Filter
    words = words[isascii.(words)]
    # Build ngrams
    ngrams = get_ngrams(words, ngram_size)
    # Build a ngram dictionary
    m = length(ngrams)
    n = length(words)
    ngramdict = Dict(ngrams[i]=>i for i in 1:m)
    # One-hot encode ngrams
    encmatrix = zeros(Float32, m, n)
    @inbounds for i in 1:n  # words
        _wgrams = get_ngrams(words[i], ngram_size)
        for _wg in _wgrams
            encmatrix[ngramdict[_wg], i] += 1.0
        end
    end
    # The model is the whole tuple
    return (words, ngramdict, KDTree(encmatrix, leafsize=1000), ngram_size)
end

function get_similar_words(target, words, ngramdict, model, ngram_size)
    _keys = keys(ngramdict)
    _wgrams = [wg for wg in get_ngrams(target, ngram_size) if wg in _keys]
    wordvec = zeros(Float32, length(ngramdict))
    for _wg in _wgrams
        wordvec[ngramdict[_wg]] += 1.0
    end
    return words[knn(model, wordvec, 3, true)[1]]
end

function get_ngrams(word::S, n::Int=2) where S<:AbstractString
    l = length(word)
    if l<=n
        return [word]
    else
        sz = n*div(l,n)-(n-1)
        ngrams = Vector{S}(undef, sz)
        for i in 1:sz
            ngrams[i] = word[i:i+n-1]
        end
        #push!(ngrams, word[n*div(l,n)-n+2:end])
    end
    return ngrams
end

function get_ngrams(words::Vector{S}, n::Int=2) where S<:AbstractString
    ngrams = S[]
    for word in words
        l = length(word)
        if l<=n
            push!(ngrams, word)
        else
            for i in 1:n*div(l,n)-(n-1)
                push!(ngrams, word[i:i+n-1])
            end
            #push!(ngrams, word[n*div(l,n)-n+2:end])
        end
    end
    return unique(ngrams)
end

=#


text:
using BenchmarkTools, PEPS, ITensors, ITensorsGPU

suite = BenchmarkGroup()
Nx = 4
Ny = 4
χ  = 6
D  = 6
sites = siteinds("S=1/2",Nx*Ny)
A     = randomfPEPS(sites, Nx, Ny, mindim=D)
H     = PEPS.makeH_XXZ(Nx, Ny, 1.0)
is_cu = PEPS.is_gpu(A)
include("environments.jl")
include("ancillaries.jl")
include("gauge.jl")
include("hamiltonian.jl")

tune!(suite)
results = run(suite, verbose=true, seconds=20)
display(median(results))


text:
using NDTensors
using LinearAlgebra
using Test

const 𝟎 = NDTensors.EmptyNumber()

@testset "NDTensors.EmptyNumber" begin
  x = 2.3

  @test complex(𝟎) == 𝟎
  @test complex(NDTensors.EmptyNumber) == Complex{NDTensors.EmptyNumber}

  # Basic arithmetic
  @test 𝟎 + 𝟎 == 𝟎
  @test 𝟎 + x == x
  @test x + 𝟎 == x
  @test -𝟎 == 𝟎
  @test 𝟎 - 𝟎 == 𝟎
  @test x - 𝟎 == x
  @test 𝟎 * 𝟎 == 𝟎
  @test x * 𝟎 == 𝟎
  @test 𝟎 * x == 𝟎
  @test 𝟎 / x == 𝟎
  @test_throws DivideError() x / 𝟎 == 𝟎
  @test_throws DivideError() 𝟎 / 𝟎 == 𝟎

  @test float(𝟎) == 0.0
  @test float(𝟎) isa Float64
  @test norm(𝟎) == 0.0
  @test norm(𝟎) isa Float64
end


text:
sal = SAL()

push!(sal, (1,3,4,5))
push!(sal, (2,3,4,5))
push!(sal, (1,3,6))
println(values(sal.nodelist[1]))

println(collect(SimplexTrees.maximalsimplices(sal)))

@testset "SAL: DictSet" begin
    ds = SimplexTrees.DictSet{Int}()

    # Insertions.
    insert!(ds, 1, 1)
    insert!(ds, 1, 2)
    insert!(ds, 1, 3)
    @test ds[1] == IntSet([1, 2, 3])
    @test ds[2] == IntSet()
    @test ds[3] == IntSet()

    insert!(ds, 2, 2)
    @test ds[1] == IntSet([1, 2, 3])
    @test ds[2] == IntSet(2)
    @test ds[3] == IntSet()

    # Deletions, isempty
    delete!(ds, 1)
    @test ds[1] == IntSet()
    @test ds[2] == IntSet(2)

    insert!(ds, 1, 1)
    insert!(ds, 1, 2)
    insert!(ds, 1, 3)
    delete!(ds, 1, 2)
    @test ds[1] == IntSet([1, 3])
    @test ds[2] == IntSet(2)
    @test !isempty(ds)

    delete!(ds, 1, 1)
    delete!(ds, 1, 3)
    delete!(ds, 2)
    @test isempty(ds)

    # Deletions with sets
    insert!(ds, 1, 1)
    insert!(ds, 1, 2)
    insert!(ds, 1, 3)
    insert!(ds, 2, 3)
    insert!(ds, 3, 2)
    delete!(ds, IntSet([2, 3]))
    @test ds[1] == IntSet(1)
    @test ds[2] == IntSet()
    @test ds[3] == IntSet()

    delete!(ds, IntSet([1, 2, 3]))
    @test isempty(ds)
end

@testset "isface" begin
    @test isface((1, 2, 3), (1, 2, 3))
    @test isface([3, 6, 7], 1:10)
    @test isface([70, 80], 1:100)
    @test !isface([70, 80], [1])
    @test !isface([1, 6, 7], 1:6)
    @test !isface([1, 6, 7], [1, 2, 3, 4, 5, 7, 8, 9])
    @test !isface([1, 6, 7], [1, 7])
end

@testset "maximalsimplices" begin
    sal = SAL()
    push!(sal, (1, 2, 3, 4))
    push!(sal, (3, 4, 5))
    maxsxs = collect(values(SimplexTrees.maximalsimplices(sal)))
    @test [3, 4, 5] in maxsxs
    @test [1, 2, 3, 4] in maxsxs
    @test !([1, 2, 3] in maxsxs)
end

@testset "Labels" begin
    # Make sure lbl_vals ∩ lbl_holes == {}
end

@testset "SAL: in" begin
    sal = SAL()
    push!(sal, (1,3,4,5))
    push!(sal, (2,3,4,5))
    push!(sal, (1,3,6))

    @test (1, 3, 4, 5) in sal
    @test (2, 3, 4, 5) in sal
    @test (1, 3, 6) in sal

    @test (1, 3) in sal
    @test (1, 6) in sal
    @test (2, 3, 5) in sal

    @test !((2, 3, 6) in sal)
    @test !((1, 7) in sal)

    @test_throws ArgumentError (2, 1) in sal
    @test_throws ArgumentError (7, 2) in sal

    sal = SAL()
    push!(sal, (1, 2))
    push!(sal, (2, 3))
    push!(sal, (1, 3))
    @test !((1, 2, 3) in sal)
end

@testset "SAL: Insertion" begin
end

@testset "SAL: Removal" begin
    @testset "removemaximals" begin
        sal = SAL()
        push!(sal, (1, 2, 3))
        push!(sal, (2, 3, 4))
        push!(sal, (1, 2, 5))
        SimplexTrees.deletemaximals!(sal, IntSet([1, 2, 3]))

        @test isempty(sal)

        push!(sal, (1, 2, 3))
        push!(sal, (2, 3, 4))
        push!(sal, (1, 2, 5))

        SimplexTrees.deletemaximals!(sal, IntSet([1, 3]))
        @test length(SimplexTrees.maximalsimplices(sal)) == 1
    end
end


text:
# forked from FixedSizeStrings.jl

import Base: iterate, lastindex, getindex, sizeof, length, ncodeunits, codeunit, isvalid, read, write, setindex!, string, convert

struct Characters{N} <: AbstractString
    data::NTuple{N, UInt8}
    function Characters{N}(v::Vector{UInt8}) where N
        new(NTuple{N, UInt8}(v))
    end
    function Characters{N}(itr) where {N}
        isempty(itr) && return missing
        new(NTuple{N, UInt8}(rpad(itr, N)))
    end
end

function Characters{N}(v::Vector{UInt8}, v2) where N

    @simd for i in 1:min(N, length(v))
        @inbounds v2[i] = v[i]
    end
    @simd for i in length(v)+1:N
        @inbounds v2[i] = 0x20
    end

    Characters{N}(v2)
end
function Characters{N}(v::Vector{UInt8}, v2, st, en) where N
    o1 = min(N, en-st+1)
    copyto!(v2, 1, v, st, o1)
    @simd for i in o1+1:N
        @inbounds v2[i] = 0x20
    end

    Characters{N}(v2)
end

Characters(s::Characters) = s

function Characters(s::AbstractString)
    isempty(s) && return missing
    Characters{ncodeunits(s)}(collect(codeunits(s)))

end

macro c_str(str)
    Characters(str)
end

function Base.String(s::T) where T <: Characters
    len = ncodeunits(s)
    out = Base._string_n(len)
    ref = Ref{T}(s)
    GC.@preserve ref out begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_copyto!(pointer(out), ptr, len)
    end
    return out
end

function Base.print(io::IO, s::T) where T<:Characters
    print(io, String(s))
end
Base.string(s::Characters) = String(s)

function Base.:(==)(s1::Characters, s2::Characters)
    cmp(s1,s2) == 0
end

function Base.:(==)(s1::Characters, s2::AbstractString)
    return view(codeunits(s1), 1:length(s1)) == codeunits(s2)

end
Base.:(==)(s1::AbstractString, s2::Characters) = s2 == s1

Base.isequal(s1::Characters, s2::Characters) =  cmp(s1, s2) == 0#s1 == s2
function Base.isequal(s1::Characters, s2::AbstractString)
    return isequal(view(codeunits(s1), 1:length(s1)), codeunits(s2))

end
Base.isequal(s1::AbstractString, s2::Characters) = isequal(s2, s1)

function Base.isless(s1::Characters, s2::Characters)
    cmp(s1,s2)<0
end

function Base.isless(s1::Characters, s2::AbstractString)
    return isless(view(codeunits(s1), 1:length(s1)), codeunits(s2))
end
function Base.isless(s1::AbstractString, s2::Characters)
    return isless(codeunits(s1), view(codeunits(s2), 1:length(s2)))
end

function iterate(s::Characters{N}, i::Int = 1) where N
    i > length(s) && return nothing
    return (Char.(s.data[i]), i+1)
end

lastindex(s::Characters{N}) where {N} = length(s)

getindex(s::Characters, i::Int) = Char(s.data[i])

sizeof(s::Characters) = sizeof(s.data)
sizeof(::Type{Characters{N}}) where N = N

function length(s::Characters)
    s_end = length(s.data)
    @inbounds for i in length(s.data):-1:1
        s.data[i] == 0x20 ? s_end -= 1 : break
    end
    return s_end
end

ncodeunits(s::Characters) = length(s)

codeunit(::Type{Characters{N}}) where N = UInt8
codeunit(::Characters) = UInt8
codeunit(s::Characters, i::Integer) = s.data[i]

isvalid(s::Characters, i::Int) = checkbounds(Bool, s, i)

Characters(s::Symbol) = Characters(string(s))

Characters(::Missing) = missing
Characters{N}(::Missing) where N = missing

function read(io::IO, T::Type{Characters{N}}) where N
    return read!(io, Ref{T}())[]::T
end

function write(io::IO, s::Characters{N}) where N
    return write(io, Ref(s))
end

function Base.hash(s::Characters{N}, h::UInt) where N
    h += Base.memhash_seed
    ref = Ref(s.data)
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), ref, length(s), h % UInt32) + h
end

function Base.cmp(a::Characters, b::Characters)
    a === b && return 0
    a, b = Iterators.Stateful(a), Iterators.Stateful(b)
    for (c, d) in zip(a, b)
        c ≠ d && return ifelse(c < d, -1, 1)
    end
    isempty(a) && return ifelse(isempty(b), 0, -1)
    return 1
end


text:
#=
error_model.jl

    Error model specifications, with accompanying utility routines and update
    and initialization steps for the EM estimation routine.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/03/01
=#

"""
    AbstractErrorModel

Abstract type for error model specifications.
"""
abstract type AbstractErrorModel end

"""
    Independent <: AbstractErrorModel

Constructor for a independent error specification instance of the error model
type, with hyperparameters `σ` and `ω`.
"""
struct Independent{Tε, Tσ} <: AbstractErrorModel
    ε::Tε   # errors
    σ::Tσ   # variances
    ω::Tσ   # precisions
end
# Constructor
function Independent(ε::AbstractMatrix)
    # Get dims
    n= size(ε, 1)

    # covariance and precision vectors
    σ= similar(ε, n)
    ω= similar(σ)

    return Independent(ε, σ, ω)
end

"""
    Idiosyncratic <: AbstractErrorModel

Constructor for a idiosyncratic error specification instance of the error model
type, with hyperparameters `Σ` and `Ω`.
"""
struct Idiosyncratic{Tε, TΣ <: Symmetric} <: AbstractErrorModel
    ε::Tε   # errors
    Σ::TΣ   # covariance matrix
    Ω::TΣ   # precision matrix
end
# Constructor
function Idiosyncratic(ε::AbstractMatrix)
    # Get dims
    n= size(ε, 1)

    # covariance and precision matrices
    Σ= Symmetric(similar(ε, n, n))
    Ω= similar(Σ)

    return Idiosyncratic(ε, Σ, Ω)
end

"""
    SpatialErrorModel <: AbstractErrorModel

Constructor for a spatial error model specification instance of the error model
type, with hyperparameter `ρ`, idiosyncratic error specification `error`,
spatial weight matrix `W`, and group structure `groups`.
"""
struct SpatialErrorModel{Tε, TΣ, Tρ, Tf, TW, TG, Tg, Te} <: AbstractErrorModel
    ε::Tε       # spatial errors
    Σ::TΣ       # covariance matrix
    Ω::TΣ       # precision matrix
    ρ::Tρ       # spatial dependence
    ρ_max::Tf   # upper bound on absolute value of ρ
    W::TW       # spatial weight matrix
    G::TG       # spatial lag polynomial
    groups::Tg  # group structure
    error::Te   # idiosyncratic error specification
end
# Constructor
function SpatialErrorModel( ε::AbstractMatrix,
                            ρ::AbstractVector,
                            W::AbstractMatrix,
                            groups::AbstractVector,
                            error::AbstractErrorModel
                        )
    # dims
    n= sum(groups)

    # covariance and precision matrices
    Σ= Symmetric(similar(ε, n, n))
    Ω= similar(Σ)

    # spatial lag polynomial
    G= similar(ρ, n, n)

    # Upper bound on abs. val. of ρ to ensure G is row and column bounded
    ρ_max= max(inv(opnorm(W, 1)), inv(opnorm(W, Inf)))

    return SpatialErrorModel(ε, Σ, Ω, ρ, ρ_max, W, G, groups, error)
end

"""
    SpatialMovingAverageModel <: AbstractErrorModel

Constructor for a spatial moving average model specification instance of the
error model type, with hyperparameter `ρ`, idiosyncratic error specification
`error`, spatial weight matrix `W`, and group structure `groups`.
"""
struct SpatialMovingAverageModel{Tε, TΣ, Tρ, TW, TG, Tg, Te} <: AbstractErrorModel
    ε::Tε       # spatial errors
    Σ::TΣ       # covariance matrix
    Ω::TΣ       # precision matrix
    ρ::Tρ       # spatial dependence
    W::TW       # spatial weight matrix
    G::TG       # spatial MA polynomial
    groups::Tg  # group structure
    error::Te   # idiosyncratic error specification
end
# Constructor
function SpatialMovingAverageModel( ε::AbstractMatrix,
                                    ρ::AbstractVector,
                                    W::AbstractMatrix,
                                    groups::AbstractVector,
                                    error::AbstractErrorModel
                                )
    # dims
    n= sum(groups)

    # covariance and precision matrices
    Σ= Symmetric(similar(ε, n, n))
    Ω= similar(Σ)

    # spatial MA polynomial
    G= similar(ρ, n, n)

    return SpatialMovingAverageModel(ε, Σ, Ω, ρ, W, G, groups, error)
end

# Methods
nparams(model::Independent)= length(model.σ)
nparams(model::Idiosyncratic)= length(model.Σ) * (length(model.Σ) + 1) ÷ 2
nparams(model::SpatialErrorModel)= length(model.ρ) + nparams(model.error)
nparams(model::SpatialMovingAverageModel)= length(model.ρ) + nparams(model.error)

get_params!(ψ::AbstractVector, model::Independent)= ψ.= model.σ
function get_params!(ψ::AbstractVector, model::Idiosyncratic)
    k= 0
    @inbounds @fastmath for j in axes(model.Σ,2)
        for i in 1:j
            k+= 1
            ψ[idx+k]= model.Σ[i,j]
        end
    end

    return nothing
end
function get_params!(ψ::AbstractVector, model::SpatialErrorModel)
    get_params!(view(ψ,1:nparams(model.error)), model.error)
    ψ[nparams(model.error)+1:end].= model.ρ

    return nothing
end
function get_params!(ψ::AbstractVector, model::SpatialMovingAverageModel)
    get_params!(view(ψ,1:nparams(model.error)), model.error)
    ψ[nparams(model.error)+1:end].= model.ρ

    return nothing
end

store_params!(model::Independent, ψ::AbstractVector)= model.σ.= ψ
function store_params!(model::Idiosyncratic, ψ::AbstractVector)
    k= 0
    @inbounds @fastmath for j in axes(model.Σ,2)
        for i in 1:j
            k+= 1
            model.Σε[i,j]= ψ[idx+k]
        end
    end

    return nothing
end
function store_params!(model::SpatialErrorModel, ψ::AbstractVector)
    n= length(ψ)
    store_params!(model.error, view(ψ,1:nparams(model.error)))
    model.ρ.= view(ψ,nparams(model.error)+1:n)

    return nothing
end
function store_params!(model::SpatialMovingAverageModel, ψ::AbstractVector)
    n= length(ψ)
    store_params!(model.error, view(ψ,1:nparams(model.error)))
    model.ρ.= view(ψ,nparams(model.error)+1:n)

    return nothing
end

"""
    resid(model)

Retrieve residuals of error model

#### Arguments
  - `model::AbstractErrorModel` : error model
"""
resid(model::AbstractErrorModel)= model.ε

"""
    cov(model)

Retrieve covariance matrix of error model

#### Arguments
  - `model::AbstractErrorModel` : error model
"""
cov(model::AbstractErrorModel)= model.Σ
cov(model::Independent)= Diagonal(model.σ)

"""
    prec(model)

Retrieve precision matrix of error model

#### Arguments
  - `model::AbstractErrorModel` : error model
"""
prec(model::AbstractErrorModel)= model.Ω
prec(model::Independent)= Diagonal(model.ω)

"""
    spatial(model)

Retrieve spatial dependence parameter of error model

#### Arguments
  - `model::AbstractErrorModel` : error model
"""
spatial(model::AbstractErrorModel)= nothing
spatial(model::SpatialErrorModel)= model.ρ
spatial(model::SpatialMovingAverageModel)= model.ρ

"""
    init_ρ!(ρ, ε, W, groups, ρ_max)

Initialize spatial dependence parameter `ρ` using spatial lag regression of
residuals `ε`, storing the result in `ρ`.

#### Arguments
  - `ε::AbstractMatrix`     : residuals
  - `W::AbstractMatrix`     : spatial weight matrix
  - `groups::AbstractVector`: group structure
  - `ρ_max::Real`           : upper bound on abs. value

#### Returns
  - `ρ::AbstractVector` : spatial dependence
"""
function init_ρ!(   ρ::AbstractVector,
                    ε::AbstractMatrix,
                    W::AbstractMatrix,
                    groups::AbstractVector,
                    ρ_max::Real
                )
    # squared residuals
    ε_sq= ε * transpose(ε)

    idx= 1
    @inbounds @fastmath for i in 1:length(ρ)
        rng= idx:idx+groups[i]-1
        # numerator and denominator
        num= zero(eltype(ρ))
        denom= zero(eltype(ρ))
        for j in rng
            w= view(W,j,:)
            num+= dot(w, view(ε_sq,:,j))
            denom+= dot(w, ε_sq, w)
        end
        # ρᵢ
        ρ_tmp= num * inv(denom)
        ρ[i]= sign(ρ_tmp) * min(abs(ρ_tmp), .9 * ρ_max)
        idx+= groups[i]
    end

    return nothing
end

"""
    spatial_ar_polynomial!(G, ρ, W, groups)

Compute spatial autoregressive lag polynomial matrix `G` based on spatial
dependence `ρ` and spatial weights `W`, storing the result in `G`.

#### Arguments
  - `ρ::AbstractVector`     : spatial dependence
  - `W::AbstractMatrix`     : spatial weight matrix
  - `groups::AbstractVector`: group structure

#### Returns
  - `G::AbstractMatrix` : spatial AR lag polynomial
"""
function spatial_ar_polynomial!(G::AbstractMatrix,
                                ρ::AbstractVector,
                                W::AbstractMatrix,
                                groups::AbstractVector
                                )
    idx= 1
    @inbounds @fastmath for i in 1:length(ρ)
        rng= idx:idx+groups[i]-1
        G[rng,:].= -ρ[i] .* view(W,rng,:)
        idx+= groups[i]
    end
    @inbounds @fastmath for i in axes(G,1)
        G[i,i]+= one(eltype(G))
    end

    return nothing
end

"""
    spatial_ma_polynomial!(G, ρ, W, groups)

Compute spatial moving average lag polynomial matrix `G` based on spatial
dependence `ρ` and spatial weights `W`, storing the result in `G`.

#### Arguments
  - `ρ::AbstractVector`     : spatial dependence
  - `W::AbstractMatrix`     : spatial weight matrix
  - `groups::AbstractVector`: group structure

#### Returns
  - `G::AbstractMatrix` : spatial MA lag polynomial
"""
function spatial_ma_polynomial!(G::AbstractMatrix,
                                ρ::AbstractVector,
                                W::AbstractMatrix,
                                groups::AbstractVector
                                )
    idx= 1
    @inbounds @fastmath for i in 1:length(ρ)
        rng= idx:idx+groups[i]-1
        G[rng,:].= ρ[i] .* view(W,rng,:)
        idx+= groups[i]
    end
    @inbounds @fastmath for i in axes(G,1)
        G[i,i]+= one(eltype(G))
    end

    return nothing
end

"""
    init_error!(model)

Initialize the error model hyper parameters, storing the results in `model`.

#### Arguments
  - `model::AbstractErrorModel` : error model
"""
function init_error!(model::Independent, init::NamedTuple)
    # variances and precisions
    if haskey(init, :error)
        model.σ.= init.error.σ
        model.ω.= init.error.ω
    else
        T= size(model.ε,2)
        @inbounds @fastmath for i in axes(model.σ,1)
            ε_i= view(model.ε,i,:)
            model.σ[i]= inv(T) * sum(abs2, ε_i)
            model.ω[i]= inv(model.σ[i])
        end
    end

    return nothing
end
function init_error!(model::Idiosyncratic, init::NamedTuple)
    # variances and precisions
    if haskey(init, :error)
        model.Σ.= init.error.Σ
        model.Ω.= init.error.Ω
    else
        T= size(model.ε,2)
        mul!(model.Σ.data, model.ε, transpose(model.ε))
        lmul!(T, model.Σ.data)
        model.Ω.= model.Σ
        LinearAlgebra.inv!(cholesky!(model.Ω))
    end

    return nothing
end
function init_error!(model::SpatialErrorModel, init::NamedTuple)
    # spatial dependence
    haskey(init, :error) ? model.ρ.= init.error.ρ : init_ρ!(model.ρ, resid(model), model.W, model.groups, model.ρ_max)

    # spatial lag polynomial
    spatial_ar_polynomial!(model.G, model.ρ, model.W, model.groups)

    # idsiosyncratic errors
    mul!(resid(model.error), model.G, resid(model))
    init_error!(model.error, init)

    # variance and precision
    model.Ω.data.= transpose(model.G) * prec(model.error) * model.G
    model.Σ.= model.Ω
    LinearAlgebra.inv!(cholesky!(model.Σ))

    return nothing
end
function init_error!(model::SpatialMovingAverageModel, init::NamedTuple)
    # spatial dependence
    model.ρ.= haskey(init, :error) ? init.error.ρ : zero(eltype(model.ρ))

    # spatial MA polynomial
    spatial_ma_polynomial!(model.G, model.ρ, model.W, model.groups)

    # factorization
    fac= lu(model.G)

    # idsiosyncratic errors
    ldiv!(resid(model.error), fac, resid(model))
    init_error!(model.error, init)

    # variance and precision
    model.Σ.data.= model.G * cov(model.error) * transpose(model.G)
    model.Ω.= model.Σ
    LinearAlgebra.inv!(cholesky!(model.Ω))

    return nothing
end

# Estimation
# Log-likelihood
"""
    f_spatial(ρ, model, quad)

Compute objective function value `f` w.r.t. spatial dependence parameter ``ρ``
(average negative log-likelihood w.r.t. ``ρ``) for either a spatial error or
spatial moving average model.

#### Arguments
  - `ρ::AbstractVector`         : spatial dependence
  - `model::AbstractErrorModel` : error model
  - `quad::AbstractMatrix`      : quadratic smoother state space model component

#### Returns
  - `f::Real`   : objective function value
"""
function f_spatial( ρ::AbstractVector,
                    model::SpatialErrorModel,
                    quad::AbstractMatrix
                )
    # Get dims
    (n,T)= size(resid(model))

    # transform parameters back
    ρ.= logistic.(ρ; offset=model.ρ_max, scale=2 * model.ρ_max)
    # Spatial lag polynomial G
    spatial_ar_polynomial!(model.G, ρ, model.W, model.groups)
    # transform parameters
    ρ.= logit.(ρ; offset=model.ρ_max, scale=2 * model.ρ_max)

    # error component
    e= resid(model) * transpose(resid(model)) .+ quad

    # precision matrix
    tmp_nn= transpose(model.G) * prec(model.error) # buffer
    mul!(prec(model).data, tmp_nn, model.G)

    # objective function
    f= zero(eltype(ρ))
    @inbounds @fastmath for i in 1:n
        f+= dot(view(prec(model),:,i), view(e,:,i))
    end
    # logabsdet
    d,s= logabsdet(lu!(model.G))

    return f * inv(2 * T) - d - log(s)
end
function f_spatial( ρ::AbstractVector,
                    model::SpatialMovingAverageModel,
                    quad::AbstractMatrix
                )
    # Get dims
    (n,T)= size(resid(model))

    # transform parameters back
    ρ.= logistic.(ρ, offset=1., scale=2.)
    # Spatial lag polynomial G
    spatial_ma_polynomial!(model.G, ρ, model.W, model.groups)
    # transform parameters
    ρ.= logit.(ρ, offset=1., scale=2.)

    # error component
    e= resid(model) * transpose(resid(model)) .+ quad

    # variance and precision
    cov(model).data.= model.G * cov(model.error) * transpose(model.G)
    prec(model).= cov(model)
    LinearAlgebra.inv!(cholesky!(prec(model)))

    # objective function
    f= zero(eltype(ρ))
    @inbounds @fastmath for i in 1:n
        f+= dot(view(prec(model),:,i), view(e,:,i))
    end
    # logabsdet
    d,s= logabsdet(lu!(model.G))

    return f * inv(2 * T) + d + log(s)
end

# Gradient
"""
    ∇f_spatial!(∇f, ρ, model, quad)

Compute gradient `∇f` of objective function `f` w.r.t. spatial dependence ``ρ``
(gradient of average negative log-likelihood w.r.t. ``ρ``) for either a spatial
error or spatial moving average model.

#### Arguments
  - `ρ::AbstractVector`         : spatial dependence
  - `model::AbstractErrorModel` : error model
  - `quad::AbstractMatrix`      : quadratic smoother state space model component

#### Returns
  - `∇f::AbstractVector`: gradient
"""
function ∇f_spatial!(   ∇f::AbstractVector,
                        ρ::AbstractVector,
                        model::SpatialErrorModel,
                        quad::AbstractMatrix
                    )
    # Get dims
    T= size(resid(model),2)

    # transform parameters back
    ρ.= logistic.(ρ; offset=model.ρ_max, scale=2 * model.ρ_max)
    # Spatial lag polynomial G
    spatial_ar_polynomial!(model.G, ρ, model.W, model.groups)
    # transform parameters
    ρ.= logit.(ρ; offset=model.ρ_max, scale=2 * model.ρ_max)

    # error component
    tmp_nn= resid(model) * transpose(resid(model)) .+ quad
    e= tmp_nn * transpose(model.W)

    # variance component
    mul!(tmp_nn, transpose(model.G), prec(model.error))

    # log-determinant component
    ∇logdet= model.W / lu!(model.G)

    # gradient
    idx= 1
    @inbounds @fastmath for i in 1:length(ρ)
        rng= idx:idx+model.groups[i]-1
        ∇f_i= zero(eltype(∇f))
        for j in rng
            ∇f_i+= ∇logdet[j,j] - inv(T) * dot(view(e,:,j), view(tmp_nn,:,j))
        end
        jacob= 2 * model.ρ_max * logistic(ρ[i]) * (one(ρ[i]) - logistic(ρ[i]))
        ∇f[i]= ∇f_i * jacob
        idx+= model.groups[i]
    end

    return nothing
end
function ∇f_spatial!(   ∇f::AbstractVector,
                        ρ::AbstractVector,
                        model::SpatialMovingAverageModel,
                        quad::AbstractMatrix
                    )
    # Get dims
    T= size(resid(model),2)

    # transform parameters back
    ρ.= logistic.(ρ, offset=1., scale=2.)
    # Spatial lag polynomial G
    spatial_ma_polynomial!(model.G, ρ, model.W, model.groups)
    # transform parameters
    ρ.= logit.(ρ, offset=1., scale=2.)

    # variance and precision
    tmp_nn= cov(model.error) * transpose(model.G)
    mul!(cov(model).data, model.G, tmp_nn)
    prec(model).= cov(model)
    LinearAlgebra.inv!(cholesky!(prec(model)))

    # log-determinant component
    ∇logdet= model.W / lu!(model.G)

    # error component
    tmp_nn.= quad
    mul!(tmp_nn, resid(model), transpose(resid(model)), 1., 1.)
    e= tmp_nn * transpose(∇logdet)

    # gradient
    idx= 1
    @inbounds @fastmath for i in 1:length(ρ)
        rng= idx:idx+model.groups[i]-1
        ∇f_i= zero(eltype(∇f))
        for j in rng
            ∇f_i+= ∇logdet[j,j] - inv(T) * dot(view(e,:,j), view(prec(model),:,j))
        end
        jacob= 2 * logistic(ρ[i]) * (one(ρ[i]) - logistic(ρ[i]))
        ∇f[i]= ∇f_i * jacob
        idx+= model.groups[i]
    end

    return nothing
end

# EM
"""
    update_error!(model, quad, pen)

Update error model hyper parameters, storing the results in `model`.

#### Arguments
  - `model::AbstractErrorModel` : error model
  - `quad::AbstractMatrix`      : quadratic smoother state space model component
  - `pen::Penalization`         : penalization variables
"""
function update_error!(model::Independent, quad::AbstractMatrix, pen::Penalization)
    T= size(resid(model),2)
    @inbounds @fastmath for i in axes(model.σ,1)
        ε_i= view(resid(model),i,:)
        # variance
        model.σ[i]= inv(T) * ( sum(abs2, ε_i) + quad[i,i] )
        # precision
        model.ω[i]= inv(model.σ[i])
    end

    return nothing
end
function update_error!(model::Idiosyncratic, quad::AbstractMatrix, pen::Penalization)
    # variance
    Σ= cov(model).data
    mul!(Σ, resid(model), transpose(resid(model)))
    @inbounds @fastmath for j in axes(Σ,2)
        for i in 1:j
            Σ[i,j]+= quad[i,j]
        end
    end
    T= size(resid(model),2)
    lmul!(inv(T), Σ)

    # precision
    prec(model).= Σ
    LinearAlgebra.inv!(cholesky!(prec(model)))

    return nothing
end
function update_error!(model::SpatialErrorModel, quad::AbstractMatrix, pen::Penalization)
    # Idiosyncratic components
    mul!(resid(model.error), model.G, resid(model))
    quad_spat= model.G * quad * transpose(model.G)
    update_error!(model.error, quad_spat, pen)

    # Spatial dependence
    # Closure of functions
    f(x::AbstractVector)= f_spatial(x, model, quad)
    ∇f!(∇f::AbstractVector, x::AbstractVector)=  ∇f_spatial!(∇f, x, model, quad)

    # Initial value for proximal operator
    x0= logit.(model.ρ; offset=model.ρ_max, scale=2 * model.ρ_max)

    # Proximal operators
    prox_g!(x::AbstractVector, λ::Real)= prox!(x, λ, pen)
    prox_f!(x::AbstractVector, λ::Real)= smooth!(x, λ, f, ∇f!, x0)

    # Transform parameters
    model.ρ.= logit.(model.ρ; offset=model.ρ_max, scale=2 * model.ρ_max)

    # Penalized estimation via admm
    model.ρ.= admm!(model.ρ, prox_f!, prox_g!)

    # Transform parameters back
    model.ρ.= logistic.(model.ρ; offset=model.ρ_max, scale=2 * model.ρ_max)

    # Spatial lag polynomial G
    spatial_ar_polynomial!(model.G, model.ρ, model.W, model.groups)

    # variance and precision
    prec(model).data.= transpose(model.G) * prec(model.error) * model.G
    cov(model).= prec(model)
    LinearAlgebra.inv!(cholesky!(cov(model)))

    return nothing
end
function update_error!(model::SpatialMovingAverageModel, quad::AbstractMatrix, pen::Penalization)
    # inverse
    Ginv= inv(model.G)

    # idiosyncratic components
    mul!(resid(model.error), Ginv, resid(model))
    quad_spat= Ginv * quad * transpose(Ginv)
    update_error!(model.error, quad_spat, pen)

    # Spatial dependence
    # Closure of functions
    f(x::AbstractVector)= f_spatial(x, model, quad)
    ∇f!(∇f::AbstractVector, x::AbstractVector)=  ∇f_spatial!(∇f, x, model, quad)

    # Initial value for proximal operator
    x0= logit.(model.ρ, offset=1., scale=2.)

    # Proximal operators
    prox_g!(x::AbstractVector, λ::Real)= prox!(x, λ, pen)
    prox_f!(x::AbstractVector, λ::Real)= smooth!(x, λ, f, ∇f!, x0)

    # Transform parameters
    model.ρ.= logit.(model.ρ, offset=1., scale=2.)

    # Penalized estimation via admm
    model.ρ.= admm!(model.ρ, prox_f!, prox_g!)

    # Transform parameters back
    model.ρ.= logistic.(model.ρ, offset=1., scale=2.)

    # Spatial lag polynomial G
    spatial_ma_polynomial!(model.G, model.ρ, model.W, model.groups)

    # variance and precision
    cov(model).data.= model.G * cov(model.error) * transpose(model.G)
    prec(model).= cov(model)
    LinearAlgebra.inv!(cholesky!(prec(model)))

    return nothing
end


text:
d = readdlm("input.txt", ',')
steps = Dict("n" => [0, -1], "ne" => [1, -1], "se" => [1, 0],
             "s" => [0, 1], "sw" => [-1, 1], "nw" => [-1, 0])
loc = [0, 0]
dist(x) = Int((abs(x[1]) + abs(x[2]) + abs(sum(x))) / 2)
maxDist = 0
for x in d
    loc += steps[x]
    maxDist = max(maxDist, dist(loc))
end
println(dist(loc))
println(maxDist)


text:
using PyPlot
using ComplexEdu
using PlotEdu
if !isdefined(:colorCurrent1)
  using PlotColorTab20
end

close("all")
figure(figsize=(3.3,2.8))
rc("text",usetex=true);
rc("text.latex",preamble="\\usepackage{fourier}\\usepackage[scaled=0.95]{inconsolata}");
rc("font", family="serif")

rel=0.2
t1=collect(-1:0.01:-rel);
t2=collect(-rel:0.01:rel);
t3=collect(rel:0.01:2);

y1=fill(-1.0,size(t1));
y3=fill(+1.0,size(t3));
y2=t2*y3[1]/t3[1];

plot(t1,y1,color="blue",linestyle=lineStyle1,linewidth=lineWidth1)
plot(t2,y2,color="blue",linestyle=lineStyle1,linewidth=lineWidth1)
plot(t3,y3,color="blue",linestyle=lineStyle1,linewidth=lineWidth1)
xlim(-1,1)
ylim(-1.1,1.1)
yticks([-1,0,1])
xticks([-rel,0,rel],["",L"$\texttt{-v\_linear~~0~~v\_linear}$",""])
grid(true)
subplots_adjust(left=0.2, right=0.85, top=0.95, bottom=0.20)
xlabel(L"$\rightarrow~\texttt{v}$")
ylabel(L"$\rightarrow~\frac{\texttt{fRollingFriction}}{\texttt{f\_frictionSliding}}$")

save3fig("v_linear",dpi=96)

ylabel(L"$\rightarrow~\frac{\texttt{f\_constant}}{\texttt{f\_frictionSliding}}$")
save3fig("v_linearBase",dpi=96)

max_forks_repo_name: christiankral/ElectroMechanicalDrives


text:
function packmultihypo(fnc::CommonConvWrapper{T}) where {T<:AbstractFactor}
  @warn "packmultihypo is deprecated in favor of Vector only operations"
  fnc.hypotheses !== nothing ? string(fnc.hypotheses) : ""
end
function parsemultihypostr(str::AS) where {AS <: AbstractString}
  @warn "parsemultihypostr is deprecated in favor of Vector only operations"
  mhcat=nothing
  if length(str) > 0
    mhcat = convert(SamplableBelief, str)
  end
  return mhcat
end

## packing converters-----------------------------------------------------------
# heavy use of multiple dispatch for converting between packed and original data types during DB usage

function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: FactorOperationalMemory}
  return PackedFunctionNodeData(d.eliminated, d.potentialused, d.edgeIDs,
                                convert(P, _getCCW(d).usrfnc!),
                                d.multihypo, _getCCW(d).certainhypo, d.nullhypo,
                                d.solveInProgress, d.inflation)  # extract two values from ccw for storage -- ccw thrown away
end

## unpack converters------------------------------------------------------------

# FIXME see #1424
function convert(
            ::Type{GenericFunctionNodeData{CommonConvWrapper{F}}},
            packed::GenericFunctionNodeData{P} ) where {F <: AbstractFactor, P <: PackedInferenceType}
  #
  # TODO store threadmodel=MutliThreaded,SingleThreaded in persistence layer
  usrfnc = convert(F, packed.fnc)
  mhcat, nh = parseusermultihypo(packed.multihypo, packed.nullhypo)

  # TODO -- improve _prepCCW for hypotheses and certainhypo field recovery when deserializing
  # reconstitute from stored data
  # FIXME, add threadmodel=threadmodel
  # FIXME https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/590#issuecomment-776838053
  # FIXME dont know what manifolds to use in ccw
  ccw = _prepCCW(DFG.DFGVariable[], usrfnc, multihypo=mhcat, nullhypo=nh, inflation=packed.inflation)
  ccw.certainhypo = packed.certainhypo

  # CommonConvWrapper{typeof(usrfnc)}
  ret = FunctionNodeData{typeof(ccw)}(packed.eliminated, packed.potentialused, packed.edgeIDs, ccw,
                                                            packed.multihypo, packed.certainhypo, packed.nullhypo,
                                                            packed.solveInProgress, packed.inflation )
  #
  return ret
end

## Variables

"""
  $(SIGNATURES)
After deserializing a factor using decodePackedType, use this to
completely rebuild the factor's CCW and user data.
Dev Notes:
- TODO: We should only really do this in-memory if we can by without it (review this).
"""
function rebuildFactorMetadata!(dfg::AbstractDFG{SolverParams},
                                factor::DFGFactor,
                                neighbors = map(vId->getVariable(dfg, vId), getNeighbors(dfg, factor)) )
  #
  # Set up the neighbor data

  # Rebuilding the CCW
  fsd = getSolverData(factor)
  fnd_new = getDefaultFactorData( dfg,
                                  neighbors,
                                  getFactorType(factor),
                                  multihypo=fsd.multihypo,
                                  nullhypo=fsd.nullhypo,
                                  # special inflation override
                                  inflation=fsd.inflation,
                                  eliminated=fsd.eliminated,
                                  potentialused=fsd.potentialused,
                                  edgeIDs=fsd.edgeIDs,
                                  solveInProgress=fsd.solveInProgress)
  #

  factor_ = if typeof(fnd_new) != typeof(getSolverData(factor))
    # must change the type of factor solver data FND{CCW{...}}
    # create a new factor
    factor__ = DFGFactor(getLabel(factor),
                        getTimestamp(factor),
                        factor.nstime,
                        getTags(factor),
                        fnd_new,
                        getSolvable(factor),
                        Tuple(getVariableOrder(factor)))
    #

    # replace old factor in dfg with a new one
    deleteFactor!(dfg, factor, suppressGetFactor=true)
    addFactor!(dfg, factor__)

    factor__
  else
    setSolverData!(factor, fnd_new)
    # We're not updating here because we don't want
    # to solve cloud in loop, we want to make sure this flow works:
    # Pull big cloud graph into local -> solve local -> push back into cloud.
    # updateFactor!(dfg, factor)
    factor
  end

  #... Copying neighbor data into the factor?
  # JT TODO it looks like this is already updated in getDefaultFactorData -> _prepCCW
  # factormetadata.variableuserdata is deprecated, remove when removing deprecation
  # for i in 1:Threads.nthreads()
  #   ccw_new.fnc.cpt[i].factormetadata.variableuserdata = deepcopy(neighborUserData)
  # end

  return factor_
end

# import IncrementalInference: decodefg, loadjld

veeCategorical(val::Categorical) = val.p
veeCategorical(val::Union{Nothing, Vector{Float64}}) = val

function convert( ::Type{Tuple{ManifoldKernelDensity,Float64}},
                  p::TreeBelief )
  #
  (convert(ManifoldKernelDensity, p), p.infoPerCoord)
end

#


text:
function neural_ode(model,x,tspan,
                    args...;kwargs...)
  p = destructure(model)
  dudt_(u::TrackedArray,p,t) = restructure(model,p)(u)
  dudt_(u::AbstractArray,p,t) = Tracker.data(restructure(model,p)(u))
  prob = ODEProblem(dudt_,x,tspan,p)
  return diffeq_adjoint(p,prob,args...;u0=x,kwargs...)
end

function neural_ode_rd(model,x,tspan,
                       args...;kwargs...)
  dudt_(u,p,t) = model(u)
  _x = x isa TrackedArray ? x : param(x)
  prob = ODEProblem(dudt_,_x,tspan)
  # TODO could probably use vcat rather than collect here
  solve(prob, args...; kwargs...) |> Tracker.collect
end

function neural_dmsde(model,x,mp,tspan,
                      args...;kwargs...)
  dudt_(u,p,t) = model(u)
  g(u,p,t) = mp.*u
  prob = SDEProblem(dudt_,g,param(x),tspan,nothing)
  # TODO could probably use vcat rather than collect here
  solve(prob, args...; kwargs...) |> Tracker.collect
end


text:
# Autogenerated wrapper script for StarPU_jll for i686-w64-mingw32-libgfortran5
export libstarpu

using Hwloc_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("StarPU")
JLLWrappers.@declare_library_product(libstarpu, "libstarpu-1.3-5.dll")
function __init__()
    JLLWrappers.@generate_init_header(Hwloc_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libstarpu,
        "bin\\libstarpu-1.3-5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Bogosort

function isordered{T}(a::AbstractArray{T,1})
    if length(a) < 2
        return true
    end
    for i in 2:length(a)
        if a[i] < a[i-1]
            return false
        end
    end
    return true
end

function bogosort!{T}(a::AbstractArray{T,1})
    while !isordered(a)
        shuffle!(a)
    end
    return a
end

a = [rand(-10:10) for i in 1:10]
println("Before bogosort:")
println(a)
bogosort!(a)
println("\nAfter bogosort:")
println(a)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
Pkg.add("Gadfly")
using Gadfly
include("Jlo.jl")

x = linspace(-1,1,100)
g = Jlo.gaussian(x,0,1)

X,Y = Jlo.meshgrid(x,x)
G = Jlo.gaussian(X,0,0.1).*Jlo.gaussian(Y,0,0.1)

E0 = G
L_m = 10e-6
lambda_m = 1e-6
z_m = 30e-3
Eprop = Jlo.propTF(E0,L_m,lambda_m,z_m)
spy(abs(Eprop))

max_stars_repo_name: awojdyla/jlo

max_issues_repo_name: awojdyla/jlo

max_forks_repo_name: awojdyla/jlo


text:
"""
    IC methods
    zeros(e,d)
    rand(e,d)
    e => equation type
    d => domain type
"""
Base.zeros(eqs::AbstractEquations,d) = fill!(similar(eqs,d),0)
Random.rand(eqs::Union{NL,GQL},d::Domain{T},aη::T=1e-6) where T =  aη*exp.(im*rand!(Uniform(0,2π),similar(eqs,d)))
Random.rand(eqs::Union{CE2,GCE2},d::Domain{T},aη::T=1e-6) where T = convert(eqs,rand(NL(),d,aη),d)

function fullrank(eqs,d::Domain{T},aη::T=1e-6) where T
    @info "Setting a full rank (CE2) or random (NL) initial condition..."
    u = zeros(eqs,d)
    for m=1:d.nx-1
        for n=-d.ny+1:d.ny-1
            u.x[2][n+d.ny,n+d.ny,m] = aη
        end
    end
    u
end


text:
# tests for staggered grid

function test_staggered()

  fname = "input_vals_staggered.jl"
  mesh, sbp, eqn, opts = run_solver(fname)

  # test volume integrals when both grids are same
  opts["use_staggered_grid"] = false
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)

  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.calcVolumeIntegralsStaggered(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  # test face integrals
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.evalFaceIntegrals(mesh, sbp, eqn, opts)

  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.calcFaceIntegralsStaggered_nopre(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts, eqn.flux_func)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  # use different operators on the grids
  delete!(opts, "use_staggered_grid")
  opts["operator_type2"] = "SBPDiagonalE"
  opts["order2"] = 3
  opts["IC_name"] = "ICp1"
  opts["order"] = 2

  fname = make_input(opts, "input_vals_staggered2")

  mesh, sbp, eqn, opts = run_solver(fname)

  # check if the D operators are the same
  println("checking D")
  for d=1:2
    sbp2 = mesh.sbp2
    Hinv = inv(diagm(sbp2.w))
    D_tilde = mesh.I_F2S*Hinv*sbp2.Q[:, :, d]*mesh.I_S2F

    Hinv = inv(diagm(sbp.w))
    D_exact = Hinv*sbp.Q[:, :, d]
  end

  println("checking different operators")
  # because the operators are exact in this case, the residuals should be the
  # same
  opts["use_staggered_grid"] = false
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.evalVolumeIntegrals(mesh, sbp, eqn, opts)

  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.calcVolumeIntegralsStaggered(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  # test face integrals
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.evalFaceIntegrals(mesh, sbp, eqn, opts)

  res_orig = copy(eqn.res)
  fill!(eqn.res, 0.0)
  AdvectionEquationMod.calcFaceIntegralsStaggered_nopre(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts, eqn.flux_func)

  @test isapprox( vecnorm(eqn.res - res_orig), 0.0) atol=1e-13

  return nothing
end

add_func1!(AdvectionTests, test_staggered, [TAG_SHORTTEST, TAG_STAGGERED])


text:
# This file was generated, do not modify it. # hide
classif = machine(KNNClassifier(K=3), Xs, y)
fit!(classif, rows=train)
ŷ = predict_mode(classif, rows=test)

accuracy(ŷ, y[test]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials