text:
abstract type Analysis end

struct FlowAnalysis <: Analysis
    reach::Dict
    ancestors::Dict
    sites::Set{Variable}
    addrs::Set{Any}
    map::Dict
    inv::Dict
    loops
    ir::IR
end

@inline function get_successors(ra::FlowAnalysis, var)
    haskey(ra.reach, var) || return []
    ra.reach[var]
end

@inline function get_ancestors(ra::FlowAnalysis, var)
    haskey(ra.ancestors, var) || return []
    ra.ancestors[var]
end

@inline function get_variable_by_address(ra::FlowAnalysis, addr)
    haskey(ra.inv, addr) || return (nothing, false)
    return (getfield(ra, :inv)[addr], true)
end

function Base.display(ra::FlowAnalysis)
    println(" __________________________________\n")
    println("              IR Reference\n")
    display(ra.ir)
    println("\n ------------ Loops ------------\n")
    display(ra.loops)
    println("\n ------------ Addresses ------------\n")
    println(ra.addrs)
    println()
    println(" ------------ Variable to address ------------\n")
    display(ra.map)
    println()
    println(" ------------ Address to variable ------------\n")
    display(ra.inv)
    println("\n ------------ Reachability ------------\n")
    for x in ra.sites
        haskey(ra.reach, x) ? println(" $x => $(ra.reach[x])") : println(" $x")
    end
    println("\n ------------ Ancestors ------------\n")
    for x in ra.sites
        haskey(ra.ancestors, x) ? println(" $x => $(ra.ancestors[x])") : println(" $x")
    end
    println(" __________________________________\n")
end

mutable struct CallGraph <: Analysis
    addresses::Set{Address}
    dependencies::Dict{Address, Set{Address}}
    CallGraph(addrs, map) = new(addrs, map)
    CallGraph() = new(Dict{Address, Set{Address}}())
end

function Base.display(cg::CallGraph)
    println("  __________________________________\n")
    println("             Dependencies\n")
    for x in cg.addresses
        haskey(cg.dependencies, x) ? println(" $x => $(cg.dependencies[x])") : println(" $x")
    end
    println("  __________________________________\n")
end

# Reaching analysis.
function reaching!(reach::Vector{Variable}, p, var, ir)
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                reaching!(reach, p, v, ir)
            end
        end
    end
end

function reaching(var::Variable, ir)
    reach = Variable[]
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                reaching!(reach, var, v, ir)
            end
        end
    end
    return reach
end

function transitive_closure!(work, reach, s)
    for (k, v) in reach
        if s in v && k != s
            push!(work, k)
            transitive_closure!(work, reach, k)
        end
    end
end

function flow_analysis(ir)
    sites = Set(Variable[])
    addrs = Set(Any[])
    var_addr_map = Dict{Variable, Any}()
    reach = Dict{Variable, Any}()
    for (v, st) in ir
        MacroTools.postwalk(st) do e
            @capture(e, call_(sym_, args__))
            if call isa GlobalRef && call.name == :trace
                push!(sites, v)
                if sym isa QuoteNode
                    sym = sym.value
                end
                push!(addrs, sym)
                var_addr_map[v] = sym
                reach[v] = Set(reaching(v, ir))
            elseif call == trace
                push!(sites, v)
                push!(addrs, unwrap(sym))
                var_addr_map[v] = unwrap(sym)
                reach[v] = Set(reaching(v, ir))
            else
                push!(sites, v)
                reach[v] = Set(reaching(v, ir))
            end
            e
        end
    end
    ancestors = Dict()
    for s in sites
        work = Set(Variable[])
        for (k, v) in reach
            s in v && begin
                push!(work, k)
                transitive_closure!(work, reach, k)
            end
        end
        ancestors[s] = work
    end
    return FlowAnalysis(reach, ancestors, sites, addrs, var_addr_map, Dict( v => k for (k, v) in var_addr_map), detectloops(ir), ir)
end

function dependency(a::Analysis)
    map = Dict{Any, Set{Any}}()
    for (k, vars) in a.reach
        depends = Any[]
        for v in vars
            haskey(a.map, v) && push!(depends, a.map[v].value)
        end
        map[a.map[k].value] = Set(depends)
    end
    addrs = [i.value for i in a.addrs]
    return CallGraph(Set(addrs), map)
end

# ------------ Driver ------------ #

# Returns the dependency analysis in call graph (tree) form.
function construct_graph!(parent, addr, call, type)
    ir = lower_to_ir(call, type)
    if control_flow_check(ir)
        graph = CallGraph()
    else
        analysis = flow_analysis(ir)
        dependencies = dependency(analysis)
        graph = CallGraph(dependencies)
    end
    parent[addr] = graph
end

function construct_graph(ir::IRTools.IR)
    if !control_flow_check(ir)
        graph = CallGraph()
    else
        analysis = flow_analysis(ir)
        dependencies = dependency(analysis)
        graph = dependencies
    end
    return graph
end

# Toplevel analysis driver.
function construct_graph(call::Function, type...)
    ir = lower_to_ir(call, type...)
    construct_graph(ir)
end

function construct_graph(call::Function)
    ir = lower_to_ir(call)
    construct_graph(ir)
end


text:
# This file was generated, do not modify it. # hide
codegen(post; ℓ=ℓ)

ext: jl

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
export Collect
export histmap
export binning
using StatsBase
using ArgCheck

struct Collect
    count::Int64
end

function (c::Collect)(iter)
    ret = eltype(iter)[]
    for (i,xi) in enumerate(iter)
        push!(ret, xi)
        i < c.count || break
    end
    return ret
end

function map_many_functions(getters, particles) where {N}
    map(getters) do f
        map(f, particles)
    end
end

function histmap(getters_particles...;
    edges=nothing,
    closed=:right,
    weight_function=p->p.weight,
    kw...)
    getters = Base.front(getters_particles)
    particles = Base.last(getters_particles)
    data = map_many_functions(getters, particles)
    weight = Weights(map(weight_function, particles))
    # if use_particle_weights
    #     weight = Weights(map(p->p.weight, particles))
    # else
    #     weight = Weights(fill(Float32(1), length(particles)))
    # end
    if edges == nothing
        fit(Histogram, data, weight; closed=closed, kw...)
    else
        fit(Histogram, data, weight, edges; closed=closed, kw...)
    end
end

function Base.filter(f, iter::AbstractPhspIterator; maxlength=10^7)
    if maxlength == nothing
        maxlength = -1
    end
    ret = eltype(iter)[]
    count = 0
    for p in iter
        if count == maxlength
            @warn("maxlength=$maxlength reached. Use `maxlength`=nothing to prevent this.")
            break
        end
        if f(p)
            count += 1
            push!(ret, p)
        else

        end
    end
    ret
end

const Edges{N} = NTuple{N, AbstractVector}
struct Binning{C <: AbstractArray,E <: Edges}
    edges::E
    content::C
    # Example:
    #
    # want to bin particles by their x,y coordinates
    # which are between -10 and 10 cm (-15 to 15 for y)
    # In that case we would have say
    # edges = (-10:10, -15:15)
    # content::Matrix{Vector{Particle}}
    # of size (20,30)
    # each cell of content stores the particles that have appropriate
    # x,y coordinates
end

function Base.map(f,b::Binning)
    Binning(b.edges, map(f,b.content))
end

function StatsBase.Histogram(b::Binning)
    closed = :right
    isdensity = false
    Histogram(b.edges, b.content, closed, isdensity)
end

function find_bin_index(edges, key::Number)
    find_bin_index(edges, (key,))
end

function find_bin_index(edges, key)
    map(find_bin_index, edges, key)
end

function find_bin_index(bdries::AbstractVector, key::Number)
    index = searchsortedfirst(bdries, key) - 1
    clamp(index, 1, length(bdries) - 1)
end

function binning_edges_keys_items(edges::NTuple{N}, keys, items) where {N}
    @argcheck length(keys) == length(items)
    T = eltype(items)
    C = Vector{T}
    dims = map(edges) do xs
        length(xs) - 1
    end
    content = Array{Vector{T},N}(undef, dims)
    for i in eachindex(content)
        content[i] = C()
    end
    for (key, item) in zip(keys, items)
        index = find_bin_index(edges, key)
        push!(content[index...], item)
    end
    Binning(edges, content)
end

function get_keys_edges(f, items, nbins::Nothing, edges::Nothing)
    get_keys_edges(f, items, 100, edges)
end

function get_keys_edges(f, items, nbins, edges::Nothing)
    keys = map(f, items)
    kmin,kmax = extrema(keys)
    edges = range(kmin,stop=kmax,length=nbins)
    keys, edges
end

function get_keys_edges(f, items, nbins::Nothing, edges)
    keys = map(f, items)
    keys, edges
end

function get_keys_edges(f, items, nbins, edges)
    @argcheck nbins == length(edges)
    get_keys_edges(f, items, nothing, edges)
end

function binning(f,items; nbins=nothing, edges=nothing)
    keys, edges = get_keys_edges(f, items, nbins, edges)
    binning_edges_keys_items(edges, keys, items)
end


text:
module jFEMTools

using SparseArrays
using LinearAlgebra
import Base:@propagate_inbounds
import Tensors
import Tensors: Vec,norm,dot,gradient,det, ⊗
import Base.ht_keyindex2!
import FastGaussQuadrature
import VoronoiDelaunay
import PlanarConvexHulls
import StaticArrays
import WriteVTK

# Abstract types
abstract type AbstractGalerkinElement end
abstract type AbstractVirtualElement <: AbstractGalerkinElement end
abstract type FiniteElement{dim,shape,FuncOrder,GeomOrder} <: AbstractGalerkinElement end
abstract type AbstractPolytopalMesh{dim,T} end

function Base.show(io::IO, mesh::AbstractPolytopalMesh{dim}) where {dim}
  println(io, "$dim D PolytopalMesh")
  println(io, "type: ", typeof(mesh))
  println(io, "Number of cells: ", getncells(mesh))
  println(io, "Number of vertices: ", getnvertices(mesh))
end

# Abstract Discrete Function Spaces
abstract type AbstractDiscreteFunctionSpace{dim,T} end
abstract type AbstractFEMFunctionSpace{dim,T,FE} <: AbstractDiscreteFunctionSpace{dim,T} end
abstract type AbstractVEMFunctionSpace{dim,T,FE} <: AbstractDiscreteFunctionSpace{dim,T} end

# Mesh related functions
export  rectangle_mesh, RectangleCell, TriangleCell, HexagonCell, HexahedronCell, TetrahedronCell,
        getncells, getverticesidx, getverticescoords,
        getnedges, cell_volume, cell_centroid, cell_diameter,
        mapToGlobalIdx, getvertexset, getvertexcoords,
        getnvertices, get_vertices_matrix, get_cell_connectivity_list,
  PolytopalMesh, unitSquareMesh, getnfacets
export FaceIndex, EdgeIndex, FacetIndex
export PolytopalMesh2, rectangle_mesh2, unitSquareMesh2
export hyper_rectagle_mesh2
# Assembler
export  start_assemble, assemble!,
        end_assemble
# TrialFunctions
export TrialFunction, TestFunction

# Element
export PoissonVirtualElement,LocalPoissonVirtualElement, VEMFunctionSpace

#VEM utils
export DofHandler
export VEMOperators
export assemble_stiffnessMat, assemble_massMat, assemble_load
export Dirichlet, apply!

#FEM
export CellIterator
export ContinuousLagrange
export FEMFunctionSpace

#IO
export vtk_grid

include("FEM/shapes.jl")
include("Quads/quadrature.jl")
include("tools.jl")
include("Quads/StrangQuad.jl")
include("Quads/GrundmannMoellerQuad.jl")
include("Quads/GaussQuadTetra.jl")
include("Quads/GaussQuadRect.jl")
include("Interpolations.jl")
include("mesh_operations.jl")
include("mesh.jl")
include("mesh2.jl")
include("mesh_generic_funcs.jl")
include("mesh_generators.jl")
include("mesh2_generators.jl")
include("functions.jl")
include("assembler.jl")
#include("plot_recipes.jl")
include("PoissonVirtualElement.jl")
include("DofHandler.jl")
include("VEMFunctionSpace.jl")
include("VEMOperators.jl")
include("boundary.jl")
include("FEM/basis.jl")
include("FEM/FiniteElement.jl")
include("FEM/LagrangeFE.jl")
include("FEM/FunctionSpaces.jl")
include("FEM/iterator.jl")
include("IO/vtk.jl")

# SnoopCompile output
include("precompile.jl")
_precompile_()

end # module

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
@testset "constructor" begin

    P1 = Person(age=24, field="Mechanical Engineering", courses = ["Artificial Intelligence", "Robotics"], id = 1, comment = "He is a genius")

    P2 = Person(age=18, field="Computer Engineering", GPA=4, courses=["Julia"], id = 2)

    U = University(name="Julia University", people=[P1, P2])

    U.people[2].GPA=4.2 # mutability support after Doc creation

    @test stripall(@capture_out(pprint(P1))) == stripall("""
    <person id="1">
      <age>24</age>
      <study-field>Mechanical Engineering</study-field>
      <GPA>4.5</GPA>
      <taken-courses>Artificial Intelligence</taken-courses>
      <taken-courses>Robotics</taken-courses>
      He is a genius
    </person>
    """)

    @test stripall(@capture_out(pprint(P2))) == stripall("""
    <person id="2">
      <age>18</age>
      <study-field>Computer Engineering</study-field>
      <GPA>4.2</GPA>
      <taken-courses>Julia</taken-courses>
    </person>
    """)

    @test stripall(@capture_out(pprint(U.aml))) == stripall("""
    <?xml version="1.0" encoding="UTF-8"?>
    <university university-name="Julia University">
      <person id="1">
        <age>24</age>
        <study-field>Mechanical Engineering</study-field>
        <GPA>4.5</GPA>
        <taken-courses>Artificial Intelligence</taken-courses>
        <taken-courses>Robotics</taken-courses>
        He is a genius
      </person>
      <person id="2">
        <age>18</age>
        <study-field>Computer Engineering</study-field>
        <GPA>4.2</GPA>
        <taken-courses>Julia</taken-courses>
      </person>
    </university>
    """)
end


text:
# low-memory way of accessing data by segment
function get_views(S::GphysData, inds::Array{Int64,1})
  L = Array{Int64,1}(undef,0)
  X = Array{SubArray,1}(undef,0)
  for i in inds
    n_seg = size(S.t[i],1)-1
    for k = 1:n_seg
      si = S.t[i][k,1]
      ei = S.t[i][k+1,1] - (k == n_seg ? 0 : 1)
      push!(X, view(S.x[i], si:ei))
      lx = ei-si+1
      push!(L, lx)
    end
  end
  ii = sortperm(L, rev=true)
  X = X[ii]
  L = L[ii]
  return L,X
end

function get_views(C::GphysChannel)
  L = Array{Int64,1}(undef,0)
  X = Array{SubArray,1}(undef,0)
  n_seg = size(C.t,1)-1
  for k = 1:n_seg
    si = C.t[k,1]
    ei = C.t[k+1,1] - (k == n_seg ? 0 : 1)
    push!(X, view(C.x, si:ei))
    lx = ei-si+1
    push!(L, lx)
  end
  ii = sortperm(L, rev=true)
  X = X[ii]
  L = L[ii]
  return L,X
end


text:
"""
"""
forwarddiff

function forwarddiff(bc, n)
    opn = _forwarddiff.(bc, n)
    eye = Eye{Bool}.(n)
    _forwarddiff(opn, eye)
end

_forwarddiff(::Dirichlet, n::Int) =
    Bidiagonal(-ones(n), ones(n-1), :U)

_forwarddiff(::Periodic, n::Int) =
    spdiagm(1-n => ones(1),
            0 => -ones(n),
            1 => ones(n-1))

_forwarddiff(opn::NTuple{1}, eye::NTuple{1}) = opn

_forwarddiff(opn::NTuple{2,Any}, eye::NTuple{2,Any}) =
    ApplyArray(kron, eye[2], opn[1]),
    ApplyArray(kron, opn[2], eye[1])

_forwarddiff(opn::NTuple{3,Any}, eye::NTuple{3,Any}) =
    ApplyArray(kron, eye[3], eye[2], opn[1]),
    ApplyArray(kron, eye[3], opn[2], eye[1]),
    ApplyArray(kron, opn[3], eye[2], eye[1])

"""
"""
backwarddiff

function backwarddiff(bc, n)
    opn = _backwarddiff.(bc, n)
    eye = Eye{Bool}.(n)
    _backwarddiff(opn, eye)
end

_backwarddiff(::Dirichlet, n::Int) =
    Bidiagonal(ones(n), -ones(n-1), :L)

_backwarddiff(::Periodic, n::Int) =
    spdiagm(-1 => -ones(n-1),
            0 => ones(n),
            n - 1 => -ones(1))

_backwarddiff(opn::NTuple{1}, eye::NTuple{1}) = opn

_backwarddiff(opn::NTuple{2,Any}, eye::NTuple{2,Any}) =
    ApplyArray(kron, eye[2], opn[1]),
    ApplyArray(kron, opn[2], eye[1])

_backwarddiff(opn::NTuple{3,Any}, eye::NTuple{3,Any}) =
    ApplyArray(kron, eye[3], eye[2], opn[1]),
    ApplyArray(kron, eye[3], opn[2], eye[1]),
    ApplyArray(kron, opn[3], eye[2], eye[1])


text:
# This script provides an implementation of
# the 4th-order Hermite's method for solving
# a system of differential equations and
# compares it to the exact solutions and
# other methods while solving a 1-body problem

println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra
using Printf

# Define the value of ϰ
const ϰ = -1

# Integrate equations of motion using the
# 4th-order Runge-Kutta's method, return
# the values of position and velocity on
# the last step
function rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F<:AbstractFloat,I<:Unsigned}
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Compute the solution
    for _ = 1:n
        k₁ = v
        a₁ = ϰ * r / norm(r)^3
        k₂ = v + h * a₁ / 2
        a₂ = ϰ * (r + h * k₁ / 2) / norm(r + h * k₁ / 2)^3
        k₃ = v + h * a₂ / 2
        a₃ = ϰ * (r + h * k₂ / 2) / norm(r + h * k₂ / 2)^3
        k₄ = v + h * a₃
        a₄ = ϰ * (r + h * k₃) / norm(r + h * k₃)^3
        r += h / 6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)
        v += h / 6 * (a₁ + 2 * a₂ + 2 * a₃ + a₄)
    end
    return r, v
end

# Integrate equations of motion using the
# 4th-order Hermite's method, return
# the values of position and velocity on
# the last step
function hermite(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F<:AbstractFloat,I<:Unsigned}
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Define the acceleration and jerk functions
    acc(r) = ϰ * r / norm(r)^3
    jerk(r, v) = ϰ * (v / norm(r)^3 - 3 * (r ⋅ v) .* r / norm(r)^5)
    # Compute the solution
    for _ = 1:n
        # Save old values
        rₒ, vₒ = copy(r), copy(v)
        aₒ, jₒ = acc(r), jerk(r, v)
        # Predict position and velocity
        r += v * h + aₒ * (h^2 / 2) + jₒ * (h^3 / 6)
        v += aₒ * h + jₒ * (h^2 / 2)
        # Predict acceleration and jerk
        a, j = acc(r), jerk(r, v)
        # Correct velocity and position
        v = vₒ + (aₒ + a) * (h / 2) + (jₒ - j) * (h^2 / 12)
        r = rₒ + (vₒ + v) * (h / 2) + (aₒ - a) * (h^2 / 12)
    end
    return r, v
end

println(" "^4, "> Integrating the 1-body problem...")

# Define the initial values of the position and velocity
r₀ = [1.0, 0.0]
v₀ = [0.0, 0.5]

# Calculate integrals and some of the orbit parameters
E₀ = 1 / 2 * norm(v₀)^2 - ϰ^2 / norm(r₀)

# Print the initial values
println(
    '\n',
    " "^6, "r₀: ", r₀[1], " ", r₀[2], '\n',
    " "^6, "v₀: ", v₀[1], " ", v₀[2], '\n',
    " "^6, "E₀: ", E₀,
)

# Define the output directories
const data_dir = "$(@__DIR__)/../data/1-body"
const files = joinpath.(
    data_dir,
    ["rk4.dat", "hermite.dat"]
)

# Prepare a list of methods
const methods = (rk4, hermite)

# Integrate equations of motion using all available
# methods and print the values of position and velocity
# on the last step
function integrate(h::Vector{F}, n::Vector{I}) where {F<:AbstractFloat,I<:Unsigned}
    # Open the data files
    io = open.(files, "a")

    # Prepare results vectors
    r = [Vector{Vector{F}}() for _ = 1:length(io)]
    v = [Vector{Vector{F}}() for _ = 1:length(io)]
    E = [Vector{F}() for _ = 1:length(io)]
    ΔE = [Vector{F}() for _ = 1:length(io)]

    # For each pair of parameters
    for i in eachindex(h)
        # Integrate and get the results of the last steps
        for j in eachindex(io)
            rᵢ, vᵢ = methods[j](r₀, v₀, h[i], n[i])

            push!(r[j], rᵢ)
            push!(v[j], vᵢ)

            Eᵢ = 1 / 2 * norm(vᵢ)^2 - ϰ^2 / norm(rᵢ)
            ΔEᵢ = abs(E₀ - Eᵢ)

            push!(E[j], Eᵢ)
            push!(ΔE[j], ΔEᵢ)
        end
    end

    # Get the format string based on the index of the iteration
    function get_format(i)
        Printf.Format(
            "\$ 10^{-$(length("$(n[i])") - 1)} \$ " *
            "& \$ $(n[i] % 10 == 0 ? "10^$(length("$(n[i])") - 1)" : n[i]) \$ " *
            "& \$ %.14f \$ " *
            "& \$ %.14f \$ " *
            "\\\\\n"
        )
    end

    # Print the data in the files in the LaTeX format

    for j in eachindex(io)
        println(io[j], "# Positions")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), r[j][i][1], r[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Velocities")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), v[j][i][1], v[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Total energy")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), E[j][i], ΔE[j][i])
        end
        println(io[j])
        close(io[j])
    end
end

# Truncate the previous results
open.(files; truncate = true)

# Integrate n = 2^m iterations
integrate(
    [1e-2, 1e-3, 1e-4, 1e-5], #, 1e-6, 1e-7],
    UInt.([10^i for i = 2:5]),
)

# Integrate the full cycle
integrate(
    [1e-2, 1e-3, 1e-4, 1e-5], #, 1e-6, 1e-7],
    UInt.([271, 2714, 27141, 271408]), #, 2714081, 27140809 ]),
)

println()


text:
# Functions for plotting molecular dynamics simulation results

PLOTTING_TIME_UNIT = u"ps"
PLOTTING_TEMPERATURE_UNIT = u"K"
PLOTTING_ENERGY_UNIT = u"hartree"

# Broadcastable function to extract a property from the system
function extract_property(result::MolecularDynamicsResult, func::Function, unit::Unitful.Units)
    t -> uconvert(unit, func(result, t))
end

"""
    plot_temperature(result::MolecularDynamicsResult, stride::Integer; time_unit = u"ps", temperature_unit = u"K")::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Converts time and temperature quantities to the specified units.
"""
function plot_temperature(result::MolecularDynamicsResult, stride::Integer;
    time_unit=PLOTTING_TIME_UNIT,
    temperature_unit=PLOTTING_TEMPERATURE_UNIT
)
    N = get_num_bodies(result)
    p = plot(title="Temperature during Simulation [n = $(N)]", xlab="Time", ylab="Temperature")
    plot_temperature!(p, result, stride;
        first_plot=true, time_unit=time_unit, temperature_unit=temperature_unit)
end
"""
    plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer; first_plot::Bool = false; first_plot::Bool = false, time_unit = u"ps", temperature_unit = u"K")::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Converts time and temperature quantities to the specified units.
Adds the new line to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, adds a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer;
    first_plot::Bool=false,
    time_unit=PLOTTING_TIME_UNIT,
    temperature_unit=PLOTTING_TEMPERATURE_UNIT
)
    time_range = 1:stride:length(result)
    times = extract_property(result, get_time, time_unit).(time_range)
    first_plot || vline!(p, [times[1]]; label=false, color=:black, linestyle=:dot, lw=2)

    plot!(p, times, extract_property(result, temperature, temperature_unit).(time_range);
        label=first_plot && "Simulation Temperature", color=1)

    reference_temp = uconvert(temperature_unit, reference_temperature(result))
    ismissing(reference_temp) || plot!(p, times, fill(reference_temp, length(times));
        label=first_plot && "Reference Temperature", color=2, linestyle=:dash, lw=2)
    p
end

"""
    plot_energy(result::MolecularDynamicsResult, stride::Integer; time_unit = u"ps", energy_unit = u"hartree")::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Converts time and energy quantities to the specified units.
"""
function plot_energy(result::MolecularDynamicsResult, stride::Integer;
    time_unit=PLOTTING_TIME_UNIT,
    energy_unit=PLOTTING_ENERGY_UNIT
)
    N = get_num_bodies(result)
    p = plot(title="Energy during Simulation [n = $(N)]", xlab="Time", ylab="Energy", legend=:right)
    plot_energy!(p, result, stride;
        first_plot=true, time_unit=time_unit, energy_unit=energy_unit)
end
"""
    plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer; first_plot::Bool = false, time_unit = u"ps", energy_unit = u"hartree")::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Converts time and energy quantities to the specified units.
Add the new lines to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, add a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer;
    first_plot::Bool=false,
    time_unit=PLOTTING_TIME_UNIT,
    energy_unit=PLOTTING_ENERGY_UNIT
)
    time_range = 1:stride:length(result)
    times = extract_property(result, get_time, time_unit).(time_range)
    first_plot || vline!(p, [times[1]]; label=false, color=:black, linestyle=:dot, lw=2)

    plot!(p, times, extract_property(result, potential_energy, energy_unit).(time_range);
        label=first_plot && "Potential Energy", color=1)
    plot!(p, times, extract_property(result, kinetic_energy, energy_unit).(time_range);
        label=first_plot && "Kinetic Energy", color=2)
    plot!(p, times, extract_property(result, total_energy, energy_unit).(time_range);
        label=first_plot && "Total Energy", color=3)
end

"""
    plot_rdf(result::MolecularDynamicsResult, σ::Real, start::Integer = 1, stop::Integer = length(result))::Plot

Plot the radial distribution function of a `MolecularDynamicsResult` averaging over the timesteps in `start:stop`
Use `σ` (from Lennard Jones) as a normalization factor for the radius.
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Real, start::Integer=1, stop::Integer=length(result))
    plot_rdf(result, σ * LENGTH_UNIT, start, stop)
end
"""
    plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, start::Integer = 1, stop::Integer = length(result))::Plot

    Plot the radial distribution function of a `MolecularDynamicsResult` averaging over the timesteps in `start:stop`
        Use `σ` (from Lennard Jones) as a normalization factor for the radius.
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, start::Integer=1, stop::Integer=length(result))
    N = get_num_bodies(result)
    T = length(result)
    r, g = rdf(result, start, stop)
    @assert length(r) == length(g)

    plot(austrip.(r / σ), g;
        title="Radial Distribution Function [n = $(N)] [T = $(T)]",
        xlab="Distance r/σ",
        ylab="Radial Distribution g(r)",
        legend=false
    )
end


text:
@testset "Test configure_logging" begin
    logfile = "testlog.txt"
    msg = "test log message"
    orig = Logging.global_logger()
    logger = configure_logging(; filename = logfile)

    try
        with_logger(logger) do
            @info msg
        end

        close(logger)

        @test isfile(logfile)
        open(logfile) do io
            lines = readlines(io)
            @test length(lines) == 2  # two lines per message
            @test occursin(msg, lines[1])
        end
    finally
        rm(logfile)
        Logging.global_logger(orig)
    end
end


text:
"""
Provides 2nd-level parsing of the intermediate representation
of a  **single network** in a pnml model.
See [`PnmlModel`](@ref).

$(TYPEDEF)

# Extended

The type parameter of a nets should map directly and simply
to subtypes of [`PnmlType`](@ref).

Additional constrants can be imposed. We want to run under the motto:
"syntax is not semantics, quack".

Since a PNML.Document can contain multiple networks it is possible that a higher-level
will create multiple PetriNet instances, each a different subtype.

Multiple [`Page`](@ref) can (are permitted) be merged into one page
by [`flatten_pages!`](@ref) without losing any Petri Net semantics.
Initial concrete `PetriNet`s are constructed by flattening to a single `Page`.
"""
abstract type PetriNet{T<:PnmlType} end

nettype(petrinet::N) where {T <: PnmlType, N <: PetriNet{T}} = T

nettype(::PnmlNet{T}) where {T <: PnmlType} = T
nettype(::Page{T}) where {T <: PnmlType} = T
nettype(::Place{T}) where {T <: PnmlType} = T
nettype(::Transition{T}) where {T <: PnmlType} = T
nettype(::Arc{T}) where {T <: PnmlType} = T

nettype(::Type{PnmlNet{T}}) where {T <: PnmlType} = T
nettype(::Type{Page{T}}) where {T <: PnmlType} = T
nettype(::Type{Place{T}}) where {T <: PnmlType} = T
nettype(::Type{Transition{T}}) where {T <: PnmlType} = T
nettype(::Type{Arc{T}}) where {T <: PnmlType} = T

#------------------------------------------------------------------
# Methods that should be implemented by concrete subtypes.
#------------------------------------------------------------------

pid(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("must implement id accessor")
pages(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")

places(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")
transitions(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")
arcs(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")
refplaces(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")
reftransitions(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = error("not implemented")

places(net::PnmlNet, page_idx=1) = places(pages(net)[page_idx])
transitions(net::PnmlNet, page_idx=1) = transitions(pages(net)[page_idx])
arcs(net::PnmlNet, page_idx=1) = arcs(pages(net)[page_idx])
refplaces(net::PnmlNet, page_idx=1) = refplaces(pages(net)[page_idx])
reftransitions(net::PnmlNet, page_idx=1) = reftransitions(pages(net)[page_idx])

places(page::Page) = page.places
transitions(page::Page) = page.transitions
arcs(page::Page) = page.arcs
refplaces(page::Page) = page.refPlaces
reftransitions(page::Page) = page.refTransitions

#------------------------------------------------------------------
# ARC, INSCRIPTION
#------------------------------------------------------------------

function has_arc(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    any(x -> pid(x) === id, arcs(petrinet))
end
function arc(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    getfirst(x -> pid(x) === id, arcs(petrinet))
end
arc_ids(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = map(pid, arcs(petrinet))
function all_arcs(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    filter(a -> source(a)===id || target(a)===id, arcs(petrinet))
end
function src_arcs(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    filter(a -> source(a)===id, arcs(petrinet))
end
function tgt_arcs(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    filter(a -> target(a)===id, arcs(petrinet))
end
function inscription(petrinet::N, arc_id::Symbol) where {T <: PnmlType, N <: PetriNet{T}}
    inscription(arc(petrinet, arc_id))
end

has_arc(net::PnmlNet, id::Symbol, page_idx=1) = has_arc(pages(net)[page_idx], id)
arc(net::PnmlNet, id::Symbol, page_idx=1) = arc(pages(net)[page_idx], id)
arc_ids(net::PnmlNet, page_idx=1) = arc_ids(pages(net)[page_idx])
all_arcs(net::PnmlNet, id::Symbol, page_idx=1) = all_arcs(pages(net)[page_idx], id)
src_arcs(net::PnmlNet, id::Symbol, page_idx=1) = src_arcs(pages(net)[page_idx], id)
tgt_arcs(net::PnmlNet, id::Symbol, page_idx=1) = tgt_arcs(pages(net)[page_idx], id)
inscription(net::PnmlNet, arc_id::Symbol, page_idx=1) = inscription(pages(net)[page_idx], arc_id)

has_arc(page::Page, id::Symbol) = any(x -> pid(x) === id, arcs(page))
arc(page::Page, id::Symbol) = getfirst(x -> pid(x) === id, arcs(page))
arc_ids(page::Page) = map(pid, arcs(page))
all_arcs(page::Page, id::Symbol) = filter(a -> source(a)===id || target(a)===id, arcs(page))
src_arcs(page::Page, id::Symbol) = filter(a -> source(a)===id, arcs(page))
tgt_arcs(page::Page, id::Symbol) = filter(a -> target(a)===id, arcs(page))
inscription(page::Page, arc_id::Symbol) = inscription(arc(page, arc_id))

#------------------------------------------------------------------
# REFERENCES
#------------------------------------------------------------------

function has_refP(petrinet::N, ref_id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    any(x -> pid(x) === ref_id, refplaces(petrinet))
end
function has_refT(petrinet::N, ref_id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    any(x -> pid(x) === ref_id, reftransitions(petrinet))
end

refplace_ids(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} =
    map(pid, refplaces(petrinet))
reftransition_ids(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} =
    map(pid, reftransitions(petrinet))

function refplace(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    refplace(petrinet.net, id)
end
function reftransition(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    reftransition(petrinet.net, id)
end

has_refP(net::PnmlNet, ref_id::Symbol, page_idx=1) = has_refP(pages(net)[page_idx], ref_id)
has_refT(net::PnmlNet, ref_id::Symbol, page_idx=1) = has_refP(pages(net)[page_idx], ref_id)
refplace_ids(net::PnmlNet, page_idx=1) = refplace_ids(net.pages[page_idx])
reftransition_ids(net::PnmlNet, page_idx=1) = reftransition_ids(net.pages[page_idx])
refplace(net::PnmlNet, id::Symbol, page_idx=1) = refplace(net.pages[page_idx], id)
reftransition(net::PnmlNet, id::Symbol, page_idx=1) = reftransition(net.pages[page_idx], id)

has_refP(page::Page, ref_id::Symbol) = any(x -> pid(x) === ref_id, refplaces(page))
has_refT(page::Page, ref_id::Symbol) = any(x -> pid(x) === ref_id, reftransitions(page))
refplace_ids(page::Page) = map(pid, page.refPlaces)
reftransition_ids(page::Page) = map(pid, page.refTransitions)
refplace(page::Page, id::Symbol) = getfirst(x -> pid(x) === id, refplaces(page))
reftransition(page::Page, id::Symbol) = getfirst(x -> pid(x) === id, reftransitions(page))

#------------------------------------------------------------------
# PLACES, MARKING,
#------------------------------------------------------------------

function has_place(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    any(x -> pid(x) === id, places(petrinet))
end
function place(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    getfirst(x -> pid(x) === id, places(petrinet))
end
function marking(petrinet::N, placeid::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    marking(place(petrinet, placeid)) #TODO specialization?
end
place_ids(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} = map(pid, places(petrinet))

function initialMarking end

initialMarking(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} =
    initialMarking(petrinet, place_ids(petrinet))

function initialMarking(petrinet::N, placeid_vec::Vector{Symbol}) where {T<:PnmlType, N<:PetriNet{T}}
    LVector( (; [p=>marking(petrinet, p) for p in placeid_vec]...))
end

has_place(net::PnmlNet, id::Symbol, page=1) = has_place(net.pages[page], id)
place(net::PnmlNet, id::Symbol, page_idx=1) = place(net.pages[page_idx], id)
place_ids(net::PnmlNet, page_idx=1) = place_ids(net.pages[page_idx])

marking(net::PnmlNet, placeid::Symbol, page_idx=1) = marking(net.pages[page_idx], placeid)

initialMarking(net::PnmlNet, page_idx=1) = initialMarking(net.pages[page_idx])

function initialMarking(net::PnmlNet, placeid_vec::Vector{Symbol})
    LVector( (; [p=>marking(net, p) for p in placeid_vec]...))
end

has_place(page::Page, id::Symbol) = any(x -> pid(x) === id, places(page))
place(page::Page, id::Symbol) = getfirst(x -> pid(x) === id, places(page))
place_ids(page::Page) = map(pid, places(page))

marking(page::Page,  placeid::Symbol) = marking(place(page, placeid))

#TODO Use marking (initialiezd to initialMarking in constructor).
# Return all places' marking as LVector
initialMarking(page::Page) = initialMarking(page, place_ids(page))
function initialMarking(page::Page, placeid_vec::Vector{Symbol})
    LVector( (; [p=>marking(page, p) for p in placeid_vec]...))
end

#------------------------------------------------------------------
# TRANSITIONS, CONDITIONS
#------------------------------------------------------------------
"""
$(TYPEDSIGNATURES)

Return function to be used like: any(ispid(sym), iterater_with_pid)
"""
ispid(x) = Base.Fix2(===, x)

function has_transition(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    #any(x -> pid(x) === id, transitions(petrinet))
    any(Base.Fix2(===,id), transition_ids(petrinet))
end

function transition(petrinet::N, id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    getfirst(x -> pid(x) === id, transitions(petrinet))
end
transition_ids(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} =
    map(pid, transitions(petrinet))

has_transition(net::PnmlNet, id::Symbol, page_idx=1) = has_transition(net.pages[page_idx], id)
transition(net::PnmlNet, id::Symbol, page_idx=1) = transition(net.pages[page_idx], id)
transition_ids(net::PnmlNet, page_idx=1) = transition_ids(net.pages[page_idx])

has_transition(page::Page, id::Symbol) = any(x -> pid(x) === id, page.transitions)
transition(page::Page, id::Symbol) = getfirst(x -> pid(x) === id, transitions(page))
transition_ids(page::Page) = map(pid, page.transitions)

#----------------------------------------
conditions(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}} =
    conditions(petrinet, transition_ids(petrinet))

function conditions(petrinet::N, idvec::Vector{Symbol}) where {T<:PnmlType, N<:PetriNet{T}}
    LVector( (; [t=>condition(petrinet, t) for t in idvec]...))
end

conditions(net::PnmlNet, page_idx=1) = conditions(net.pages[page_idx])

conditions(page::Page) = conditions(page, transition_ids(page))
conditions(page::Page, idvec::Vector{Symbol}) =
    LVector( (; [t=>condition(page, t) for t in idvec]...))

function condition(petrinet::N, trans_id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
    condition(transition(petrinet, trans_id))
end

condition(net::PnmlNet, trans_id::Symbol, page_idx=1) = condition(net.pages[page_idx], trans_id)
condition(page::Page, trans_id::Symbol) = condition(transition(page, trans_id))

#-----------------------------------------------------------------
"""
$(TYPEDSIGNATURES)

Transition function of a Petri Net.
Each transition has an input vector and an output vector.
Each labelled vector is indexed by the place on the other end of the arc.
Values are inscriptions.
"""
function transition_function end

function transition_function(petrinet::N) where {T<:PnmlType, N<:PetriNet{T}}
    transition_function(petrinet, transition_ids(petrinet))
end

function transition_function(petrinet::N, idvec::Vector{Symbol}) where {T<:PnmlType, N<:PetriNet{T}}
    LVector( (; [t=>in_out(petrinet,t) for t in idvec]...))
end

transition_function(net::PnmlNet, page_idx=1) = transition_function(net.pages[page_idx])
function transition_function(page::Page, idvec::Vector{Symbol})
    LVector( (; [t=>in_out(page, t) for t in idvec]...))
end

"""
$(TYPEDSIGNATURES)

Return tuple of input, output labelled vectors with key of place ids and
value of arc inscription's value for use as a transition function.
#TODO When do these get called "pre" and "post"?
"""
function in_out end

function in_out(petrinet::N, transition_id::Symbol) where {T<:PnmlType, N<:PetriNet{T}}
        (ins(petrinet, transition_id), outs(petrinet, transition_id))
end

in_out(net::PnmlNet, transition_id::Symbol, page_idx=1) =
    in_out(net.pages[page_idx], transition_id)

function in_out(page::Page, transition_id::Symbol)
    (ins(page, transition_id), outs(page, transition_id))
end

"""
$(TYPEDSIGNATURES)
Return arcs of `p` that have `transition_id` as the target.
"""
function ins(p, transition_id::Symbol)
    LVector( (; [source(a)=>inscription(a) for a in tgt_arcs(p, transition_id)]...))
end

"""
$(TYPEDSIGNATURES)
Return arcs of `p` that have `transition_id` as the source.
"""
function outs(p, transition_id::Symbol)
    #isempty(src_arcs(p, transition_id)) && @warn "no src_arcs for $transition_id"
    LVector( (; [target(a)=>inscription(a) for a in src_arcs(p, transition_id)]...))
end

#------------------------------------------------------------------
#
#------------------------------------------------------------------


text:
using JuMP, Cbc, NPZ

function Solve_func(year, tol)

	println("Solving weights for $year ...\n\n")

	array = npzread(string(year, "_input.npz"))

	A1 = array["A1"]
	A2 = array["A2"]
	b = array["b"]

	model = Model(Cbc.Optimizer)
	set_optimizer_attribute(model, "logLevel", 1)
	N = size(A1)[2]

	@variable(model, r[1:N] >= 0)
	@variable(model, s[1:N] >= 0)

	@objective(model, Min, sum(r[i] + s[i] for i in 1:N))

	# bound on top by tolerance
	@constraint(model, [i in 1:N], r[i] + s[i] <= tol)

	# Ax = b
	@constraint(model, [i in 1:length(b)], sum(A1[i,j] * r[j] + A2[i,j] * s[j]
		                          for j in 1:N) == b[i])

	optimize!(model)
	termination_status(model)

	r_vec = value.(r)
	s_vec = value.(s)

	npzwrite(string(year, "_output.npz"), Dict("r" => r_vec, "s" => s_vec))

	println("\n")

end

year_list = [x for x in 2014:2030]
tol = 0.70

# Run solver function for all years and tolerances (in order)
for i in year_list
	Solve_func(i, tol)
end

max_stars_repo_name: chusloj/taxdata

max_issues_repo_name: chusloj/taxdata

max_forks_repo_name: chusloj/taxdata


text:
# Test complete Lanczos factorization
@testset "Complete Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs
            A = rand(T,(n,n))
            v = rand(T,(n,))
            A = (A+A')
            iter = LanczosIterator(wrapop(A), wrapvec(v), orth)
            verbosity = 1
            fact = @constinferred initialize(iter; verbosity = verbosity)
            while length(fact) < n
                @constinferred expand!(iter, fact; verbosity = verbosity)
                verbosity = 0
            end

            V = hcat(unwrapvec.(basis(fact))...)
            H = rayleighquotient(fact)
            @test normres(fact) < 10*n*eps(real(T))
            @test V'*V ≈ I
            @test A*V ≈ V*H

            @constinferred initialize!(iter, deepcopy(fact); verbosity = 1)
            states = collect(Iterators.take(iter, n)) # collect tests size and eltype?
            @test rayleighquotient(last(states)) ≈ H
        end
    end
end

# Test complete Arnoldi factorization
@testset "Complete Arnoldi factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n))
            v = rand(T,(n,))
            iter = ArnoldiIterator(wrapop(A), wrapvec(v), orth)
            verbosity = 1
            fact = @constinferred initialize(iter; verbosity = verbosity)
            while length(fact) < n
                @constinferred expand!(iter, fact; verbosity = verbosity)
                verbosity = 0
            end

            V = hcat(unwrapvec.(basis(fact))...)
            H = rayleighquotient(fact)
            factor = (orth == cgs || orth == mgs ? 250 : 10)
            @test normres(fact) < factor*n*eps(real(T))
            @test V'*V ≈ I
            @test A*V ≈ V*H

            @constinferred initialize!(iter, deepcopy(fact); verbosity = 1)
            states = collect(Iterators.take(iter, n)) # collect tests size and eltype?
            @test rayleighquotient(last(states)) ≈ H
        end
    end
end

# Test incomplete Lanczos factorization
@testset "Incomplete Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            A = (A+A')
            iter = @constinferred LanczosIterator(wrapop(A), wrapvec(v), orth)
            krylovdim = n
            fact = @constinferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @constinferred expand!(iter, fact)

                Ṽ, H, r̃, β, e = fact
                V = hcat(unwrapvec.(Ṽ)...)
                r = unwrapvec(r̃)
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ V*H + r*e'
            end

            fact = @constinferred shrink!(fact, div(n,2))
            V = hcat(unwrapvec.(@constinferred basis(fact))...)
            H = @constinferred rayleighquotient(fact)
            r = unwrapvec(@constinferred residual(fact))
            β = @constinferred normres(fact)
            e = @constinferred rayleighextension(fact)
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ V*H + r*e'
        end
    end
end

# Test incomplete Arnoldi factorization
@testset "Incomplete Arnoldi factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr)
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            iter = @constinferred ArnoldiIterator(wrapop(A), wrapvec(v), orth)
            krylovdim = 3*n
            fact = @constinferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @constinferred expand!(iter, fact)

                Ṽ, H, r̃, β, e = fact
                V = hcat(unwrapvec.(Ṽ)...)
                r = unwrapvec(r̃)
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ V*H + r*e'
            end

            fact = @constinferred shrink!(fact, div(n,2))
            V = hcat(unwrapvec.(@constinferred basis(fact))...)
            H = @constinferred rayleighquotient(fact)
            r = unwrapvec(@constinferred residual(fact))
            β = @constinferred normres(fact)
            e = @constinferred rayleighextension(fact)
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ V*H + r*e'
        end
    end
end


text:
using Breadth
using Test

@testset "Baseline" begin
    # check for no errors basically for all functions
    x = Variable("x", Statements)
    y = Variable(:y, Statements)
    @test x == Variable(:x, Statements)
    println(x)
    println(x ∨ y)
    println(x ∧ y)
    println(~x)
    println(domain(x))
    println(codomain(x))
    println(domain(x ∨ x))
    println(codomain(x ∨ x))
    println(domain(~x))
    println(codomain(~x))
end

@testset "Theorems" begin
    x = Variable("x", Statements)
    y = Variable("y", Statements)
    println(isapplicable((x ∧ y) ∨ x, x ∧ y))
    println(rewrites(x ∨ x))
end

@testset "Proving" begin
    x = Variable("x", Statements)
    y = Variable("y", Statements)
    A = Variable("A", Sets)
    B = Variable("B", Sets)
    println(result(x ∈ (A ∩ B), ((x ∈ A) ∧ (x ∈ B))))
    println(result((x ∨ True) ∧ False))
    println(result(~(x ∨ y), ~x ∧ ~y))  # De Morgan
    println(result(~x ∧ ~y, ~(x ∨ y)))

    println(result(x ∈ (A ∩ B), ((x ∈ A) ∧ (x ∈ B))))
    println(result(x ∈ (A ∩ B), x ∈ (A ∪ B)))
end

@testset "Old proving" begin
    using Breadth: x, A, B, C, result, result_wp, ⊂, TRUE, FALSE, ⟹, ∧, ∨, P, Q, rewrites, rewrites_wp
    println(result_wp(x ∈ (A ∩ B), x ∈ (A ∪ B), 3))
    """
    x ∈ A ∩ B because Given
    (x ∈ A) ∧ (x ∈ B) because Definition of Element Of
    (x ∈ A) ∨ (x ∈ B) because Or is the weaker And
    x ∈ A ∪ B because Definition of Element Of
    """
    # println(result_wp(((A ∩ B) ⊂ (A ∪ B)) ∧ (x ∈ (A ∩ B) ⟹ x ∈ (A ∪ B)), TRUE, 4))
    A ∧ (~A ∨ B)
    (A ∧ ~A) ∨ (A ∨ B)
    result_wp(P ∧ (P ⟹ Q), Q, 5)
    result_wp(P ∧ (P ⟹ Q), P ∧ (~P ∨ Q), 5)
    result_wp(P ∧ (~P ∨ Q), FALSE ∨ (P ∧ Q), 5)
    result_wp(FALSE ∨ (P ∧ Q), P ∧ Q, 5)
    result_wp((P ∧ Q) ⟹ Q, TRUE, 5)
    result_wp((P ∧ Q) ⟹ Q, (~P ∨ ~Q) ∨ Q, 5)
    result_wp((~P ∨ ~Q) ∨ Q, ~P ∨ (Q ∨ ~Q),2 )
    result_wp(~P ∨ (Q ∨ ~Q), TRUE,2 )

    # println(result_wp((A ∩ B) ⊂ (A ∪ B), TRUE, 4))  # need to add the above proof
end


text:
function calc_grad( LF, Potentials, Focc, v::Array{Float64,2} )
  #
  Npoints = size(v)[1]
  Ncol    = size(v)[2]
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
  grad = zeros( Float64,Npoints,Ncol )
  #
  for ic = 1:Ncol
    Hv = op_H( LF, Potentials, v[:,ic] )
    grad[:,ic] = Hv # copy?
    for icc = 1:Ncol
      grad[:,ic] = grad[:,ic] - dot( v[:,icc], Hv ) * v[:,icc] * ΔV
    end
    grad[:,ic] = Focc[ic] * grad[:,ic]
  end
  return grad
end

function calc_grad( LF, ∇2::SparseMatrixCSC{Float64,Int64},
                Potentials, Focc, v::Array{Float64,2} )
  #
  Npoints = size(v)[1]
  Ncol    = size(v)[2]
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
  grad = zeros(Float64,Npoints,Ncol)
  for ic = 1:Ncol
    Hv = op_H( LF, ∇2, Potentials, v[:,ic] )
    grad[:,ic] = Hv # copy?
    for icc = 1:Ncol
      grad[:,ic] = grad[:,ic] - dot( v[:,icc], Hv ) * v[:,icc] * ΔV
    end
    grad[:,ic] = Focc[ic] * grad[:,ic]
  end
  return grad
end


text:
"""
SharedCircularDeque{T}(n)

Create a double-ended queue of maximum capacity `n`, implemented as a circular buffer. The element type is `T`.
The data buffers are shared memory segments, accessible across processes on the same physical node.

Instances can be serialized and deserialized or created afresh using identical path and capacity.

Methods on the datastructure itself are not locked
    - to avoid deadlock due to internal recursive calls
    - to enable efficient batched calls
However they can be locked by the caller using the `withlock(q.lck) do ... end` syntax.
"""
struct SharedCircularDeque{T}
    path::String
    buffer::SysVSharedArray{T,1}
    state::SysVSharedArray{Int,1} # capacity, n, first, last
    lck::NamedSemaphore # lock operations between processes
end

const CAP = 1
const LEN = 2
const FST = 3
const LST = 4

function SharedCircularDeque{T}(path::String, capacity::Int; create::Bool=false) where {T}
    buffer = SysVSharedVector((path,1), capacity, T)
    state = SysVSharedVector((path,2), 4, Int)
    if create
        S = state.A
        S[CAP] = capacity # capacity
        S[LEN] = 0        # n
        S[FST] = 1        # first
        S[LST] = capacity # last
    end
    lck = NamedSemaphore("/$(hash(path))")
    SharedCircularDeque{T}(path, buffer, state, lck)
end

function Base.hash(D::SharedCircularDeque, h::UInt)
    h += 0xe4fbea67fe10ce78 % UInt
    h = hash(D.path, h)
end

Base.length(D::SharedCircularDeque) = D.state.A[LEN]
Base.endof(D::SharedCircularDeque) = D.state.A[LEN]
Base.eltype(::Type{SharedCircularDeque{T}}) where {T} = T
capacity(D::SharedCircularDeque) = D.state.A[CAP]

function Base.empty!(D::SharedCircularDeque)
    S = D.state.A
    S[LEN] = 0
    S[FST] = 1
    S[LST] = S[CAP]
    D
end

function Base.delete!(D::SharedCircularDeque)
    delete!(D.state)
    delete!(D.buffer)
    delete!(D.lck)
    nothing
end

function Base.close(D::SharedCircularDeque)
    close(D.state)
    close(D.buffer)
    close(D.lck)
    nothing
end

Base.isempty(D::SharedCircularDeque) = D.state.A[LEN] == 0

@inline function front(D::SharedCircularDeque)
    S = D.state.A
    @boundscheck S[LEN] > 0 || throw(BoundsError())
    D.buffer.A[S[FST]]
end

@inline function back(D::SharedCircularDeque)
    S = D.state.A
    @boundscheck S[LEN] > 0 || throw(BoundsError())
    D.buffer.A[S[LST]]
end

@inline function Base.push!(D::SharedCircularDeque, v)
    S = D.state.A
    @boundscheck S[LEN] < S[CAP] || throw(BoundsError()) # prevent overflow
    S[LEN] += 1
    tmp = S[LST]+1
    S[LST] = ifelse(tmp > S[CAP], 1, tmp)  # wraparound
    @inbounds D.buffer.A[S[LST]] = v
    D
end

@inline function Base.pop!(D::SharedCircularDeque)
    S = D.state.A
    v = back(D)
    S[LEN] -= 1
    tmp = S[LST] - 1
    S[LST] = ifelse(tmp < 1, S[CAP], tmp)
    v
end

@inline function Base.unshift!(D::SharedCircularDeque, v)
    S = D.state.A
    @boundscheck S[LEN] < S[CAP] || throw(BoundsError())
    S[LEN] += 1
    tmp = S[FST] - 1
    S[FST] = ifelse(tmp < 1, S[CAP], tmp)
    @inbounds D.buffer.A[S[FST]] = v
    D
end

@inline function Base.shift!(D::SharedCircularDeque)
    S = D.state.A
    v = front(D)
    S[LEN] -= 1
    tmp = S[FST] + 1
    S[FST] = ifelse(tmp > S[CAP], 1, tmp)
    v
end

@inline function Base.splice!(D::SharedCircularDeque, idx)
    S = D.state.A
    C = S[CAP]
    j = S[FST] + idx - 1
    if j > C
        j -= C
    end
    A = D.buffer.A
    @inbounds ret = A[j]

    L = S[LST]
    S[LST] = L - 1
    S[LEN] -= 1

    if L < j
        while j < C
            @inbounds A[j] = A[j+1]
            j += 1
        end
        @inbounds A[j] = A[1]
        j = 1
    end

    while j < L
        @inbounds A[j] = A[j+1]
        j += 1
    end

    return ret
end

# getindex sans bounds checking
@inline function _unsafe_getindex(D::SharedCircularDeque, i::Integer)
    S = D.state.A
    j = S[FST] + i - 1
    if j > S[CAP]
        j -= S[CAP]
    end
    @inbounds ret = D.buffer.A[j]
    return ret
end

@inline function Base.getindex(D::SharedCircularDeque, i::Integer)
    @boundscheck 1 <= i <= D.state.A[LEN] || throw(BoundsError())
    return _unsafe_getindex(D, i)
end

@inline getindex(D::SharedCircularDeque, R::UnitRange) = [D[idx] for idx in R]

# Iteration via getindex
@inline Base.start(d::SharedCircularDeque) = 1
@inline Base.next(d::SharedCircularDeque, i) = (_unsafe_getindex(d, i), i+1)
@inline Base.done(d::SharedCircularDeque, i) = i == d.state.A[LEN] + 1

function Base.show(io::IO, D::SharedCircularDeque{T}) where T
    print(io, "SharedCircularDeque{$T}([")
    for i = 1:length(D)
        print(io, D[i])
        i < length(D) && print(io, ',')
    end
    print(io, "])")
end


text:
using Plots

function1(point)=point/2
function2(point)=point/2+[1;0]
function3(point)=point/2+[1/2;sqrt(3)/2]

num=500000
p=100
finalPoints=[]
for number_of_points in 1:num
    point=[rand();rand()]
    for number_of_iteration in 1:p
        point=rand((function1(point), function2(point), function3(point)))
    end
    push!(finalPoints, point)
end
scatter(hcat(finalPoints...)[1,:], hcat(finalPoints...)[2,:],markersize=0.7, color=:purple, xlabel="x", ylabel="y")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter2\\report\\Fig\\SierpinskiRandomTriangle.png")