text:
###########################################################
#
#   Generic matrix generation functions
#
###########################################################

function matrix_from_adjpairs!(a::AbstractMatrix, g::AbstractGraph, gen; returnpermutation::Bool=false)
    @graph_requires g vertex_list vertex_map
    p = nothing
    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                a[ui, vi] = get(gen, g, u, v)
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)
                a[ui, vi] = val
                a[vi, ui] = val
            end
        end

    elseif implements_adjacency_list(g)
        p = zeros(Int,size(a,1))
        tempd = Dict{Int,Int}()
        verts = vertices(g)
        for i in 1:length(verts)
          p[i] = vertex_index(verts[i], g)
          tempd[p[i]] = i
        end
        i = 0
        for u in verts
            i+=1
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)
                a[i, tempd[vi]] = val # p[i] -> ui, tempd[vi]->i
            end
        end
    else
        error("g does not implement required interface.")
    end

    return returnpermutation ? (a, p) : a
end

matrix_from_adjpairs(g::AbstractGraph, gen; returnpermutation::Bool=false) =
    (n = num_vertices(g); matrix_from_adjpairs!(zeros(eltype(gen), n, n), g, gen, returnpermutation=returnpermutation))

function sparse_matrix_from_adjpairs(g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    n = num_vertices(g)
    m = num_edges(g)
    ne = is_directed(g) ? m : 2m
    I = Array{Int}(undef, ne)
    J = Array{Int}(undef, ne)
    vals = Array{eltype(gen)}(undef, ne)
    idx = 0

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = val
            end
        end

    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        end
    end

    @assert(idx == ne)
    sparse(I, J, vals, n, n)
end

function matrix_from_edges!(a::AbstractMatrix, g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
                a[vi, ui] = val
            end
        end

    elseif implements_incidence_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
            end
        end
    else
        error("g does not implement required interface.")
    end

    return a
end

matrix_from_edges(g::AbstractGraph, gen) =
    (n = num_vertices(g); matrix_from_edges!(zeros(eltype(gen), n, n), g, gen))

function sparse_matrix_from_edges(g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    n = num_vertices(g)
    m = num_edges(g)
    ne = is_directed(g) ? m : 2m
    I = Array{Int}(undef, ne)
    J = Array{Int}(undef, ne)
    vals = Array{eltype(gen)}(undef, ne)
    idx = 0

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = val
            end
        end

    elseif implements_incidence_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        end
    end

    @assert(idx == ne)
    sparse(I, J, vals, n, n)
end

###########################################################
#
#   Specific matrices from graphs
#
###########################################################

### adjacency matrix

mutable struct _GenUnit{T} end

Base.get(::_GenUnit{T}, g::AbstractGraph{V}, u::V, v::V) where {T,V} = one(T)
Base.eltype(::_GenUnit{T}) where {T} = T

adjacency_matrix(g::AbstractGraph, ::Type{T}; returnpermutation::Bool=false) where {T<:Number} = matrix_from_adjpairs(g, _GenUnit{T}(), returnpermutation=returnpermutation)
adjacency_matrix(g::AbstractGraph; returnpermutation::Bool=false) = adjacency_matrix(g, Bool, returnpermutation=returnpermutation)

adjacency_matrix_sparse(g::AbstractGraph, ::Type{T}) where {T<:Number} = sparse_matrix_from_adjpairs(g, _GenUnit{T}())
adjacency_matrix_sparse(g::AbstractGraph) = adjacency_matrix_sparse(g, Bool)

### weight matrix

mutable struct _GenEdgeWeight{Weights}
    weights::Weights
end

_GenEdgeWeight(weights::AbstractVector) = _GenEdgeWeight{typeof(weights)}(weights)

Base.get(gen::_GenEdgeWeight, g::AbstractGraph{V,E}, e::E) where {V,E} = gen.weights[edge_index(e, g)]
Base.eltype(gen::_GenEdgeWeight) = eltype(gen.weights)

weight_matrix(g::AbstractGraph, weights::AbstractVector) = matrix_from_edges(g, _GenEdgeWeight(weights))
weight_matrix_sparse(g::AbstractGraph, weights::AbstractVector) = sparse_matrix_from_edges(g, _GenEdgeWeight(weights))

### distance matrix

function init_distancemat(n::Int, dinf::T) where {T<:Number}
    a = fill(dinf, n, n)
    @inbounds for i = 1:n
        a[i,i] = zero(T)
    end
    return a
end

distance_matrix(g::AbstractGraph, dists::AbstractVector, dinf) =
    (n = num_vertices(g); matrix_from_edges!(init_distancemat(n, dinf), g, _GenEdgeWeight(dists)))

distance_matrix(g::AbstractGraph, dists::AbstractVector{T}) where {T<:Real} =
    distance_matrix(g, dists, typemax(T))

###########################################################
#
#   Laplacian matrix
#
###########################################################

laplacian_matrix(g::AbstractGraph) = laplacian_matrix(g, Float64)

function laplacian_matrix(g::AbstractGraph, ::Type{T}) where {T<:Number}
    @graph_requires g vertex_list vertex_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                a[ui, vi] -= 1
                a[vi, ui] -= 1
                a[ui, ui] += 1
                a[vi, vi] += 1
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                if ui < vi
                    a[ui, vi] -= 1
                    a[vi, ui] -= 1
                    a[ui, ui] += 1
                    a[vi, vi] += 1
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end
    return a
end

function laplacian_matrix(g::AbstractGraph, eweights::AbstractVector{T}) where {T<:Number}
    @graph_requires g vertex_list vertex_map edge_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                wi = eweights[edge_index(e, g)]
                a[ui, vi] -= wi
                a[vi, ui] -= wi
                a[ui, ui] += wi
                a[vi, vi] += wi
            end
        end
    else
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                if ui < vi
                    wi = eweights[edge_index(e, g)]
                    a[ui, vi] -= wi
                    a[vi, ui] -= wi
                    a[ui, ui] += wi
                    a[vi, vi] += wi
                end
            end
        end
    end
    return a
end

laplacian_matrix_sparse(g::AbstractGraph) = laplacian_matrix_sparse(g, Float64)

function laplacian_matrix_sparse(g::AbstractGraph, ::Type{T}) where {T<:Number}
    @graph_requires g vertex_list vertex_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)

    nnz = num_edges(g) * 2 + n
    I = Array{Int}(undef, nnz)
    J = Array{Int}(undef, nnz)
    vals = Array{T}(undef, nnz)
    idx = 0

    degs = zeros(T, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = -1

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = -1

                degs[ui] += 1
                degs[vi] += 1
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                if ui < vi
                    idx += 1
                    I[idx] = ui
                    J[idx] = vi
                    vals[idx] = -1

                    idx += 1
                    I[idx] = vi
                    J[idx] = ui
                    vals[idx] = -1

                    degs[ui] += 1
                    degs[vi] += 1
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end

    for i = 1:n
        di = degs[i]
        if isnz(di)
            idx += 1
            I[idx] = i
            J[idx] = i
            vals[idx] = di
        end
    end

    if idx < nnz
        I = I[1:idx]
        J = J[1:idx]
        vals = vals[1:idx]
    end
    sparse(I, J, vals, n, n)
end

function laplacian_matrix_sparse(g::AbstractGraph, eweights::AbstractVector{T}) where {T<:Number}
    @graph_requires g vertex_list vertex_map edge_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)

    nnz = num_edges(g) * 2 + n
    I = Array{Int}(undef, nnz)
    J = Array{Int}(undef, nnz)
    vals = Array{T}(undef, nnz)
    idx = 0

    degs = zeros(T, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                wi = eweights[edge_index(e, g)]
                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = -wi

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = -wi

                degs[ui] += wi
                degs[vi] += wi
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                if ui < vi
                    wi = eweights[edge_index(e, g)]
                    idx += 1
                    I[idx] = ui
                    J[idx] = vi
                    vals[idx] = -wi

                    idx += 1
                    I[idx] = vi
                    J[idx] = ui
                    vals[idx] = -wi

                    degs[ui] += wi
                    degs[vi] += wi
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end

    for i = 1:n
        di = degs[i]
        if isnz(di)
            idx += 1
            I[idx] = i
            J[idx] = i
            vals[idx] = di
        end
    end

    if idx < nnz
        I = I[1:idx]
        J = J[1:idx]
        vals = vals[1:idx]
    end
    sparse(I, J, vals, n, n)
end

###########################################################
#
#   adjacency list from matrix
#
###########################################################

function sparse2adjacencylist(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    colptr = A.colptr
    rowval = A.rowval
    n = size(A, 1)
    adjlist = Array{Array{Ti,1}}(undef, n)
    s = 0
    for j in 1:n
        adjj = Ti[]
        sizehint!(adjj, colptr[j+1] - colptr[j] - 1)
        for k in colptr[j]:(colptr[j+1] - 1)
            rvk = A.rowval[k]
            if rvk != j push!(adjj, rvk) end
        end
        s += length(adjj)
        adjlist[j] = adjj
    end
    GenericAdjacencyList{Ti, UnitRange{Ti}, Vector{Vector{Ti}}}(!ishermitian(A),
                                                             one(Ti):convert(Ti,n),
                                                             s, adjlist)
end


text:
# # HMM with discrete observations

using Distributions
using HMMBase
using Plots

# https://en.wikipedia.org/wiki/Viterbi_algorithm#Example

a = [0.6, 0.4]
A = [0.7 0.3; 0.4 0.6]
B = [Categorical([0.5, 0.4, 0.1]), Categorical([0.1, 0.3, 0.6])]
hmm = HMM(a, A, B)
#-

z, y = rand(hmm, 250, seq = true)
z_viterbi = viterbi(hmm, y);
#-

plot(y, linetype=:steppre, label="Observations", size=(600,200))
#-

plot(z, linetype=:steppre, label="True hidden state", size=(600,200))
#-

plot(z_viterbi, linetype=:steppre, label="Viterbi decoded hidden state", size=(600,200))

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ d80342ae-522b-11eb-19c0-cf931fb4ae81
using Statistics

# ╔═╡ 9c2ef32a-522d-11eb-29e6-29da29e4dec7
using Plots

# ╔═╡ 88b47ee6-522d-11eb-1a98-ed4dd43dfd11
Base.@kwdef mutable struct SampleAvg
    t::Int = 0
    avg::Float64 = 0.0
end

# ╔═╡ 8c94a810-522d-11eb-2024-01ca19bc1984
function (s::SampleAvg)(x)
    s.t += 1
    s.avg += (x - s.avg) / s.t
    s.avg
end

# ╔═╡ 8fde026e-522d-11eb-2aa1-f3e28c238ea6
function run_once(b)
    rms = Float64[]
    distribution = randn(b)
    expectation = mean(distribution)
    sample_avg = SampleAvg()

    for i in 1:2*b
        avg = sample_avg(distribution[rand(1:b)])
        push!(rms, abs(avg - expectation))
    end
    rms
end

# ╔═╡ 969b48aa-522d-11eb-2039-61916ddf87d4
begin
	n_runs = 1000
	p = plot(legend=:topright)

	for b in [2, 10, 100, 1000]
		rms = mean(run_once(b) for _ in 1:n_runs)
		xs = (1:2*b) ./ b
		plot!(p, xs, rms, label="b=$b")
	end

	p
end

# ╔═╡ Cell order:
# ╠═d80342ae-522b-11eb-19c0-cf931fb4ae81
# ╠═88b47ee6-522d-11eb-1a98-ed4dd43dfd11
# ╠═8c94a810-522d-11eb-2024-01ca19bc1984
# ╠═8fde026e-522d-11eb-2aa1-f3e28c238ea6
# ╠═9c2ef32a-522d-11eb-29e6-29da29e4dec7
# ╠═969b48aa-522d-11eb-2039-61916ddf87d4


text:
@test begin
    l = Segment(Point(0, 0), Point(3, 0))
    l(0) == l.O && l(1) == l.B && l(0.5) == 0.5*(l.O + l.B)
end

max_stars_repo_name: MiroK/fractional-dim


text:
using Markets, Test, JSON3
# implement your tests here
@test 1 == 1
using HTTP

trackFile="sampletrack.json"
url="http://localhost:8000/track"

open(trackFile, "r") do io
    trackList=JSON3.read(io)
end

HTTP.request("POST", url, [("Content-Type", "application/json")], trackFile)


text:
export FixedPointHomotopy, FixedPointHomotopyCache, gamma

"""
    FixedPointHomotopy(F, x₀; gamma=exp(i * 2π*rand()))

Construct the homotopy ``H(x, t) = (1-t)F(x) + γt(x-x₀)``.
"""
struct FixedPointHomotopy{S<:AbstractSystem,T<:AbstractVector} <: AbstractHomotopy
    F::S
    x₀::T
    gamma::Complex{Float64}

end
function FixedPointHomotopy(
    F::AbstractSystem,
    x₀::AbstractVector;
    gamma = randomish_gamma(),
)
    FixedPointHomotopy(F, x₀, gamma)
end

"""
    FixedPointHomotopyCache

An simple cache for `StartTargetHomotopy`s consisting of the caches for the
start and target system as well as a `Vector` and a `Matrix`.
"""
struct FixedPointHomotopyCache{SC} <: AbstractHomotopyCache
    F::SC
end

function cache(H::FixedPointHomotopy, x, t)
    F_cache = cache(H.F, x)
    FixedPointHomotopyCache(F_cache)
end

Base.size(H::FixedPointHomotopy) = size(H.F)

"""
    gamma(H::FixedPointHomotopy)

Obtain the gamma used in the FixedPointHomotopy.
"""
gamma(H::FixedPointHomotopy) = H.gamma

"""
    γ(H)

Obtain the gamma used in the FixedPointHomotopy.
"""
γ(H::FixedPointHomotopy) = gamma(H)

function evaluate!(u, H::FixedPointHomotopy, x, t, c::FixedPointHomotopyCache)
    evaluate!(u, H.F, x, c.F)

    u .= (γ(H) * t) .* (x .- H.x₀) .+ (1 - t) .* u

    u
end
(H::FixedPointHomotopy)(x, t, c = cache(H, x, t)) = evaluate(H, x, t, c)

function dt!(u, H::FixedPointHomotopy, x, t, c::FixedPointHomotopyCache)
    evaluate!(u, H.F, x, c.F)
    u .= γ(H) .* (x .- H.x₀) .- u
    u
end

function jacobian!(U, H::FixedPointHomotopy, x, t, c::FixedPointHomotopyCache)
    jacobian!(U, H.F, x, c.F)

    U .= (1 - t) .* U
    γt = γ(H) * t
    for i = 1:length(x)
        U[i, i] += γt
    end
    U
end

function evaluate_and_jacobian!(
    u,
    U,
    H::FixedPointHomotopy,
    x,
    t,
    c::FixedPointHomotopyCache,
)
    evaluate_and_jacobian!(u, U, H.F, x, c.F)

    u .= (γ(H) * t) .* (x .- H.x₀) .+ (1 - t) .* u
    U .= (1 - t) .* U
    γt = γ(H) * t
    for i = 1:length(x)
        U[i, i] += γt
    end
    U

    nothing
end

function jacobian_and_dt!(U, u, H::FixedPointHomotopy, x, t, c::FixedPointHomotopyCache)
    evaluate_and_jacobian!(u, U, H.F, x, c.F)

    U .= (1 - t) .* U
    γt = γ(H) * t
    for i = 1:length(x)
        U[i, i] += γt
    end
    u .= γ(H) .* (x .- H.x₀) .- u

    nothing
end


text:
using FunctionalStateMachine
using Test

## User state functions

function bar!(usrdata)
  println("do bar!")
  return FunctionalStateMachine.exitStateMachine
end

function foo!(usrdata)
  println("do foo!")
  return bar!
end

@testset "Test generic state machine..." begin

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, verbose=true); end

@test statemachine.next == emptyState

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, verbose=false); end

@test statemachine.next == emptyState

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, breakafter=foo!); end

@test statemachine.next == bar!

statemachine = StateMachine{Nothing}(next=foo!)
while statemachine(nothing, iterlimit=1); end

@test statemachine.next == bar!

statemachine = StateMachine{Nothing}(next=bar!)
while statemachine(nothing, verbose=true); end

@test statemachine.next == emptyState

end

#


text:
#== import modules here ==#
using PyCall
psi4 = pyimport("psi4")

#== build molecule ==#
psi4.set_memory("500 MB")

input_file = "/gpfs/jlse-fs0/users/davpoolechem/projects/Julia/Julia-Sandbox/QC-Programs/Psi4/S22_3/benzene_ammonia_2.xyz"

io = open(input_file, "r")
input_file_string = read(io, String)
close(io)

mol = psi4.geometry(input_file_string)
mol.reset_point_group("c1")

options = Dict(
  "REFERENCE" => "RHF",
  "BASIS" => "6-311++G(2d,2p)",
  "SCF_TYPE" => "DIRECT",
  "GUESS" => "CORE",
  "DF_SCF_GUESS" => false,
  "PUREAM" => false,
  "PROPERTIES" => [],
  "PRINT" => 0,
  "E_CONVERGENCE" => 1.0,
  "D_CONVERGENCE" => 2e-5,
  "MAXITER" => 200
)

psi4.set_options(options)
psi4.set_num_threads(56)

psi4.energy("scf/6-311++G(2d,2p)")


text:
"""
    parts(d::Verb)

Return the parts of the helper (this is usually a tuple of useful things,
e.g. a kernel generated from a query argument).
"""
parts(d::Verb) = tuple([ getfield(d, i) for i in 1:nfields(d) ]...)

# parts(helper::SelectHelper) = helper.res_field, helper.f, helper.arg_fields
# parts(helper::FilterHelper) = helper.f, helper.arg_fields
# parts(helper::GroupbyHelper) = helper.is_predicate, helper.f, helper.arg_fields
# parts(helper::SummarizeHelper) = (
#     helper.res_field,
#     helper.f,
#     helper.g,
#     helper.arg_fields
# )


text:
abstract type PetscVecBase <: AbstractVector{PetscScalar} end
"""
    A PETSc Vec wrapper.
"""
mutable struct PetscVec <: PetscVecBase

    " The pointer to the PETSc Vec object "
    vec::Ref{Vec}

    " Whether or not a size has been set by calling setSize() "
    sized::Bool

    " Whether or not the vector has ever been assembled (Note: the vector might currently NOT be assembled) "
    assembled::Bool

    name::String

    function PetscVec()
#        println("Creating PetscVec")
        vec = Ref{Vec}()
        ccall((:VecCreate, library), PetscErrorCode, (comm_type, Ref{Vec}), MPI.COMM_WORLD, vec)
        ccall((:VecSetType, library), PetscErrorCode, (Vec, VecType), vec[], VECMPI)
        new_vec = new(vec, false, false, "")

        function finalize(avec)
#            ccall((:VecDestroy, library), PetscErrorCode, (Ref{Vec},), avec.vec)
#            println(avec.name)
        end

        finalizer(finalize, new_vec)

        return new_vec
    end

    function PetscVec(name::String)
#        println("Creating ", name)
        new_vec = PetscVec()
        new_vec.name = name
        return new_vec
    end
end

"""
    A PETSc Vec wrapper for Ghosted Vectors.
"""
struct GhostedPetscVec <: PetscVecBase

    " The pointer to the PETSc Vec object "
    vec::Ref{Vec}

    " Whether or not a size has been set. "
    sized::Bool

    " Whether or not the vector has ever been assembled (Note: the vector might currently NOT be assembled) "
    assembled::Bool

    " Maps global indices to local indices "
    global_to_local_map::Dict{PetscInt, PetscInt}

    " The first index on this processor (1-based)"
    first_local_index::PetscInt

    " The last index on this processor (1-based)"
    last_local_index::PetscInt

    " Whether or not the raw array is present "
    raw_array_present::Bool

    " The raw local array.  This shouldn't be accessed directly!  Use [] to access values. (1-based)"
    raw_array::Vector{PetscScalar}

    " The local form.  Only present if the raw_array is.  Only for internal use!"
    local_form::PetscVec

    function GhostedPetscVec(ghosts::Vector{T};
                                n_local::PetscInt=PETSC_DECIDE, n_global::PetscInt=PETSC_DETERMINE) where T
        vec = Ref{Vec}()
        ghost_dofs = (PetscInt)[dof-1 for dof in ghosts]
        ccall((:VecCreateGhost, library), PetscErrorCode, (comm_type, PetscInt, PetscInt, PetscInt, Ptr{PetscInt}, Ref{Vec}),
              MPI.COMM_WORLD, n_local, n_global, length(ghost_dofs), ghost_dofs, vec)

        first = Ref{PetscInt}()
        last = Ref{PetscInt}()

        ccall((:VecGetOwnershipRange, library), PetscErrorCode, (Vec, Ref{PetscInt}, Ref{PetscInt}), vec[], first, last)

        # +1 is for 1-based indexing
        new_vec = new(vec, true, false, Dict{PetscInt, PetscInt}(), first[]+1, last[], false)

        # Need to set up the global_to_local map
        for i in 1:length(ghosts)
            new_vec.global_to_local_map[ghosts[i]] = i + (last[]-first[])
        end

        # This idea comes from libMesh
        # mapping = Ref{ISLocalToGlobalMapping}()

        # ccall((:VecGetLocalToGlobalMapping, library), PetscErrorCode, (Vec, Ref{ISLocalToGlobalMapping}), vec[], mapping)

        # println("mapping: ", mapping)

        # indices_ptr = Ref{Ptr{PetscInt}}()

        # ccall((:ISLocalToGlobalMappingGetIndices, library), PetscErrorCode, (ISLocalToGlobalMapping, Ref{Ptr{PetscInt}}), mapping[], indices_ptr)

        # println("indices_ptr: ", indices_ptr)

        # mapping_size = Ref{PetscInt}()

        # ccall((:ISLocalToGlobalMappingGetSize, library), PetscErrorCode, (ISLocalToGlobalMapping, Ref{PetscInt}), mapping[], mapping_size)

        # println("mapping_size: ", mapping_size)

        # indices = unsafe_wrap(Array, indices_ptr[], mapping_size[], false)

        # if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        #     println("indices: ", indices)
        # end

        # for i in 1:length(indices)
        #     if indices[i]+1 < new_vec.first_local_index || new_vec.last_local_index < indices[i]+1
        #         new_vec.global_to_local_map[indices[i]+1] = i
        #     end
        # end

        # if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        #     println(new_vec.global_to_local_map)
        # end

        # ccall((:ISLocalToGlobalMappingRestoreIndices, library), PetscErrorCode, (ISLocalToGlobalMapping, Ref{Ptr{PetscInt}}), mapping[], indices_ptr)

        return new_vec
    end
end

"""
    Set up the size of a PetscVecBase

    Note: This MUST be called before setting/getting values!
"""
function setSize!(vec::PetscVecBase; n_local::PetscInt=PETSC_DECIDE, n_global::PetscInt=PETSC_DETERMINE)
    @assert !vec.sized

    # Must provide _some_ size!
    @assert n_local != PETSC_DECIDE || n_global != PETSC_DETERMINE

    ccall((:VecSetSizes, library), PetscErrorCode, (Vec, PetscInt, PetscInt), vec.vec[], n_local, n_global)

    vec.sized = true
end

"""
    Must be called after setting/adding values to construct the vector
"""
function assemble!(vec::PetscVecBase)
    @assert vec.sized

    ccall((:VecAssemblyBegin, library), PetscErrorCode, (Vec,), vec.vec[])
    ccall((:VecAssemblyEnd, library), PetscErrorCode, (Vec,), vec.vec[])

    vec.assembled = true
end

"""
    Must be called after setting/adding values to construct the vector
"""
function assemble!(vec::GhostedPetscVec)
    @assert vec.sized

    _restoreArray(vec)

    ccall((:VecAssemblyBegin, library), PetscErrorCode, (Vec,), vec.vec[])
    ccall((:VecAssemblyEnd, library), PetscErrorCode, (Vec,), vec.vec[])

    ccall((:VecGhostUpdateBegin, library), PetscErrorCode, (Vec, InsertMode, ScatterMode), vec.vec[], INSERT_VALUES, SCATTER_FORWARD)
    ccall((:VecGhostUpdateEnd, library), PetscErrorCode, (Vec, InsertMode, ScatterMode), vec.vec[], INSERT_VALUES, SCATTER_FORWARD)

    vec.assembled = true
end

"""
    Use PETSc viewer to print out the vector
"""
function viewVec(vec::PetscVecBase)
    @assert vec.sized
    @assert vec.assembled

    viewer = ccall((:PETSC_VIEWER_STDOUT_, library), PetscViewer, (comm_type,), MPI.COMM_WORLD)
    ccall((:VecView, library), PetscErrorCode, (Vec, PetscViewer), vec.vec[], viewer)
end

"""
    Does vec[i] += v
"""
function plusEquals!(vec::PetscVecBase, v::Array{Float64}, i)
    i_ind = (PetscInt)[i_val-1 for i_val in i]

    @assert length(v) == length(i_ind)

    ccall((:VecSetValues, library), PetscErrorCode, (Vec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), vec.vec[], length(i_ind), i_ind, v, ADD_VALUES)
end

"""
    If the incoming type is not a float, make it so

    Does vec[i] += v
"""
function plusEquals!(vec::PetscVecBase, v::Array{T}, i) where T
    plusEquals!(vec, (Float64)[(Float64)(val) for val in v], i)
end

"""
    vec += v
"""
function plusEquals!(vec::PetscVecBase, v::PetscVecBase)
    ccall((:VecAXPY, library), PetscErrorCode, (Vec, PetscScalar, Vec), vec.vec[], 1.0, v.vec[])
end

"""
    vec = 0
"""
function zero!(vec::PetscVecBase)
    ccall((:VecZeroEntries, library), PetscErrorCode, (Vec,), vec.vec[])
end

"""
    vec = a*vec
"""
function scale!(vec::PetscVecBase, a::Real)
    ccall((:VecScale, library), PetscErrorCode, (Vec, PetscScalar), vec.vec[], (PetscScalar)(a))
end

import Base.copy!

"""
    vec = a
"""
function copy!(vec::PetscVecBase, a::PetscVecBase)
    # NOTE!  PETSc's calling sequence is _backwards_ from Julia!
    # The destination is the _second_ argument for PETSc
    ccall((:VecCopy, library), PetscErrorCode, (Vec, Vec), a.vec[], vec.vec[])
end

"""
    vec = a when "vec" is a ghosted vec
"""
function copy!(vec::GhostedPetscVec, a::PetscVecBase)
    # Basic idea here: copy over the purely local data directly
    # then use assemble!() to update the ghosts

    # This happens in five steps
    # 1. Get the local form and array for the destination (vec)
    # 2. Get the array for the src (a)
    # 3. Copy local values
    # 4. Restore everything
    # 5. Call assemble!() to update ghosts

    # 1:
    local_and_ghosted_data_array = _getArray(vec)

    # 2:
    src_local_array = _getArray(a)

    @assert length(local_and_ghosted_data_array) >= length(src_local_array)

    # 3:
    for i in 1:length(src_local_array)
        local_and_ghosted_data_array[i] = src_local_array[i]
    end

    # 4:
    _restoreArray(a, src_local_array)
    _restoreArray(vec)

    # 5:
    assemble!(vec)
end

import Base.similar
"""
    Creates a PetscVec() with the same storage as the passed in vector
"""
function similar(vec::PetscVec)
    new_vec = PetscVec()

    ccall((:VecDuplicate, library), PetscErrorCode, (Vec, Ref{Vec}), vec.vec[], new_vec.vec)

    return new_vec
end

"""
    L2 Norm
"""
function norm(vec::PetscVecBase)
    norm_value = Ref{PetscScalar}()

    ccall((:VecNorm, library), PetscErrorCode, (Vec, NormType, Ref{PetscScalar}), vec.vec[], NORM_2, norm_value)

    return norm_value[]
end

"""
    Serializes a parallel PETSc vec down into a Julia Array

    Only returns the array on processor 0
"""
function serializeToZero(vec::PetscVecBase)
    serialized_vec = PetscVec("serialized_vec")
    scatter = Ref{VecScatter}()

    ccall((:VecScatterCreateToZero, library), PetscErrorCode, (Vec, Ref{VecScatter}, Ref{Vec}), vec.vec[], scatter, serialized_vec.vec)

    ccall((:VecScatterBegin, library), PetscErrorCode, (VecScatter, Vec, Vec, InsertMode, ScatterMode), scatter[], vec.vec[], serialized_vec.vec[], INSERT_VALUES, SCATTER_FORWARD)
    ccall((:VecScatterEnd, library), PetscErrorCode, (VecScatter, Vec, Vec, InsertMode, ScatterMode), scatter[], vec.vec[], serialized_vec.vec[], INSERT_VALUES, SCATTER_FORWARD)

    # Have to copy the array because _getArray is reference memory that PETSc will destroy in a moment
    serialized_array = copy(_getArray(serialized_vec))

    ccall((:VecScatterDestroy, library), PetscErrorCode, (Ref{VecScatter},), scatter)
    ccall((:VecDestroy, library), PetscErrorCode, (Ref{Vec},), serialized_vec.vec)

    return serialized_array
end

#### AbstractArray Interface Definitions ###

"""
    PETSc Vectors are inherently 1D
"""
function linearindexing(vec::PetscVecBase)
    return Base.LinearFast()
end

import Base.size

"""
    Returns the _global_ size of the vector
"""
function size(vec::PetscVecBase)
    N = Ref{PetscInt}(0)

    ccall((:VecGetSize, library), PetscErrorCode, (Vec, Ref{PetscInt}), vec.vec[], N)

    return (N[])
end

import Base.setindex!

"""
    Sets the value at i
"""
function setindex!(vec::PetscVecBase, v, i)
    # Copy out the values
    val = (Float64)[v_val for v_val in v]

    # Call the specialization below
    setindex!(vec, val, i)
end

"""
    Sets the value at i,j.

    Specialization for when v is already an array of Float64 (faster because we don't need to copy it)
"""
function setindex!(vec::PetscVecBase, v::Array{Float64}, i)
    i_ind = (PetscInt)[i_val-1 for i_val in i]

    @assert length(v) == length(i_ind)

    ccall((:VecSetValues, library), PetscErrorCode, (Vec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), vec.vec[], length(i_ind), i_ind, v, INSERT_VALUES)
end

"""
    Don't do this
"""
function setindex!(vec::PetscVecBase, v, i, j)
    error("Attempt to index PetscVecBasetor using multiple dimensions!")
end

import Base.getindex

"""
    Don't do this either
"""
function getindex(vec::PetscVecBase, i, j)
    error("Attempt to index PetscVecBasetor using multiple dimensions!")
end

"""
    Proper getter for entries from the vector for integer indices
"""
function getindex(vec::PetscVecBase, i::T) where T<:Integer
    # Don't forget about 1-based indexing...
    i_ind = (PetscInt)[i-1]

    get_vals = Vector{Float64}(undef, 1)

    ccall((:VecGetValues, library), PetscErrorCode, (Vec, PetscInt, Ptr{PetscInt}, Ref{PetscScalar}), vec.vec[], 1, i_ind, get_vals)

    return get_vals[1]
end

"""
    Proper getter for entries from the vector
"""
function getindex(vec::PetscVecBase, i)
    i_ind = (PetscInt)[i_val-1 for i_val in i]

    get_vals = Vector{Float64}(undef, length(i_ind))

    ccall((:VecGetValues, library), PetscErrorCode, (Vec, PetscInt, Ptr{PetscInt}, Ref{PetscScalar}), vec.vec[], length(i_ind), i_ind, get_vals)

    return get_vals
end

"""
    Helper function for ghosted vector indices
"""
function _getindices(vec::GhostedPetscVec, indices)
    raw_indices = Array{PetscInt}(length(indices))

    for i in 1:length(indices)
        index = indices[i]
        if vec.first_local_index <= index && index <= vec.last_local_index # Within the local portion of the vector
            raw_indices[i] = (index - vec.first_local_index) + 1
        else # Within the ghosted part
            raw_indices[i] = vec.global_to_local_map[index]
        end
    end

    return raw_indices
end

"""
    Proper getter for entries from the vector for integer indices for Ghosted Vectors
"""
function getindex(vec::GhostedPetscVec, i::T) where T<:Integer
    if !vec.raw_array_present
        _getArray(vec)
    end

    return vec.raw_array[_getindices(vec, (PetscInt)[i])[1]]
end

"""
    Proper getter for entries from the vector for Ghosted Vectors
"""
function getindex(vec::GhostedPetscVec, i)
    if !vec.raw_array_present
        _getArray(vec)
    end

    return vec.raw_array[_getindices(vec, i)]
end

########## Private stuff

"""
    PRIVATE: Used internally.  Don't use.
"""
function _getArray(vec::PetscVecBase)
    raw_data = Ref{Ptr{PetscScalar}}()
    ccall((:VecGetArray, library), PetscErrorCode, (Vec, Ref{Ptr{PetscScalar}}), vec.vec[], raw_data)

    local_size = Ref{PetscInt}()
    ccall((:VecGetLocalSize, library), PetscErrorCode, (Vec, Ref{PetscInt}), vec.vec[], local_size)

    return unsafe_wrap(Array, raw_data[], local_size[], false)
end

"""
    PRIVATE: Used internally.  Don't use.
"""
function _getArray(vec::GhostedPetscVec)
    if !vec.raw_array_present
        vec.local_form = PetscVec("local_form")
        ccall((:VecGhostGetLocalForm, library), PetscErrorCode, (Vec, Ref{Vec}), vec.vec[], vec.local_form.vec)

        vec.raw_array = _getArray(vec.local_form)

        vec.raw_array_present = true
    end

    return vec.raw_array
end

import Base.unsafe_convert

"""
    PRIVATE: Used internally.  Don't use.
"""
function _restoreArray(vec::PetscVecBase, raw_data::Array{PetscScalar})
    # This mess is required because PETSc is expecting a PetscScalar** and Julia won't automatically
    # convert the array to that
    pointer = Ref{Ptr{PetscScalar}}()
    pointer[] = unsafe_convert(Ptr{PetscScalar}, raw_data)
    ccall((:VecRestoreArray, library), PetscErrorCode, (Vec, Ref{Ptr{PetscScalar}}), vec.vec[], pointer)
end

"""
    PRIVATE: Used internally.  Don't use.
"""
function _restoreArray(vec::GhostedPetscVec)
    if vec.raw_array_present
        _restoreArray(vec.local_form, vec.raw_array)
        ccall((:VecGhostRestoreLocalForm, library), PetscErrorCode, (Vec, Ref{Vec}), vec.vec[], vec.local_form.vec)
        vec.raw_array_present = false
    end
end


text:
"""
    mpiexec(fn)

A wrapper function for the MPI launcher executable. Calls `fn(cmd)`, where `cmd` is a `Cmd` object of the MPI launcher.

# Environment Variables

The behaviour of `mpiexec` can be controlled by the following environment variables:

- `JULIA_MPIEXEC`: the name or path of the launcher executable (set at compile time).
- `JULIA_MPIEXEC_ARGS`: additional arguments that are passed to the launcher. These are
  space seperated, supporting the same quoting rules as Julia `Cmd` objects. These can be
  modified at run time.

# Usage

```jldoctest
julia> mpiexec(cmd -> run(`\$cmd -n 3 echo hello world`));
hello world
hello world
hello world
```
"""
function mpiexec(fn)
    _mpiexec(cmd -> fn(`$cmd $(Base.shell_split(get(ENV, "JULIA_MPIEXEC_ARGS", "")))`))
end

const REFCOUNT = Threads.Atomic{Int}(-1)

"""
    refcount_inc()

Increment the MPI reference counter. This should be called at initialization of any object
which calls an MPI routine in its finalizer. A matching [`refcount_dec`](@ref) should be
added to the finalizer.

For more details, see [Finalizers](@ref).
"""
function refcount_inc()
    Threads.atomic_add!(REFCOUNT, 1)
end

"""
    refcount_dec()

Decrement the MPI reference counter. This should be added after an MPI call in an object
finalizer, with a matching [`refcount_inc`](@ref) when the object is initialized.

For more details, see [Finalizers](@ref).
"""
function refcount_dec()
    # refcount zero, all objects finalized, now finalize MPI
    if Threads.atomic_sub!(REFCOUNT, 1) == 1
        if !Finalized()
            # MPI can now be finalized, but MPI_Finalize is a collective and can act
            # like a barrier (this may be implementation specific), if we are terminating
            # due to a Julia exception, we should calling MPI_Finalize. We thus peek at the
            # current exception, and only if that field is nothing do we terminate.
            if ccall(:jl_current_exception, Any, ()) === nothing
                _Finalize()
            end
        end
    end
end

# Administrative functions
"""
    Init()

Initialize MPI in the current process.

All MPI programs must contain exactly one call to `MPI.Init` or
[`MPI.Init_thread`](@ref). In particular, note that it is not valid to call `MPI.Init` or
`MPI.Init_thread` again after calling [`MPI.Finalize`](@ref).

The only MPI functions that may be called before `MPI.Init`/`MPI.Init_thread` are
[`MPI.Initialized`](@ref) and [`MPI.Finalized`](@ref).

# External links
$(_doc_external("MPI_Init"))
"""
function Init()
    REFCOUNT[] == -1 || error("MPI.REFCOUNT in incorrect state: MPI may only be initialized once per session.")
    @mpichk ccall((:MPI_Init, libmpi), Cint, (Ptr{Cint},Ptr{Cint}), C_NULL, C_NULL)
    REFCOUNT[] = 1
    atexit(refcount_dec)

    for f in mpi_init_hooks
        f()
    end
end

"""
    ThreadLevel

An Enum denoting the level of threading support in the current process:

 - `MPI.THREAD_SINGLE`: Only one thread will execute.

 - `MPI.THREAD_FUNNELED`: The process may be multi-threaded, but the application must
   ensure that only the main thread makes MPI calls. See [`Is_thread_main`](@ref).

 - `MPI.THREAD_SERIALIZED`: The process may be multi-threaded, and multiple threads may
   make MPI calls, but only one at a time (i.e. all MPI calls are serialized).

 - `MPI.THREAD_MULTIPLE`: Multiple threads may call MPI, with no restrictions.

# See also

- [`Init_thread`](@ref)
- [`Query_thread`](@ref)
"""
@enum ThreadLevel begin
    THREAD_SINGLE     = MPI_THREAD_SINGLE
    THREAD_FUNNELED   = MPI_THREAD_FUNNELED
    THREAD_SERIALIZED = MPI_THREAD_SERIALIZED
    THREAD_MULTIPLE   = MPI_THREAD_MULTIPLE
end

"""
    Init_thread(required::ThreadLevel)

Initialize MPI and the MPI thread environment in the current process. The argument
specifies the required level of threading support, see [`ThreadLevel`](@ref).

The function will return the provided `ThreadLevel`, and values may be compared via
inequalities, i.e.

```julia
provided = Init_thread(required)
@assert provided >= required
```

All MPI programs must contain exactly one call to [`MPI.Init`](@ref) or
`MPI.Init_thread`. In particular, note that it is not valid to call `MPI.Init` or
`MPI.Init_thread` again after calling [`MPI.Finalize`](@ref).

The only MPI functions that may be called before `MPI.Init`/`MPI.Init_thread` are
[`MPI.Initialized`](@ref) and [`MPI.Finalized`](@ref).

# External links
$(_doc_external("MPI_Init_thread"))
"""
function Init_thread(required::ThreadLevel)
    REFCOUNT[] == -1 || error("MPI.REFCOUNT in incorrect state: MPI may only be initialized once per session.")
    r_provided = Ref{ThreadLevel}()
    # int MPI_Init_thread(int *argc, char ***argv, int required, int *provided)
    @mpichk ccall((:MPI_Init_thread, libmpi), Cint,
                  (Ptr{Cint},Ptr{Cvoid}, ThreadLevel, Ref{ThreadLevel}),
                  C_NULL, C_NULL, required, r_provided)
    provided = r_provided[]
    if provided < required
        @warn "Thread level requested = $required, provided = $provided"
    end

    REFCOUNT[] = 1
    atexit(refcount_dec)

    for f in mpi_init_hooks
        f()
    end
    return provided
end

"""
    Query_thread()

Query the level of threading support in the current process.
Returns a [`ThreadLevel`](@ref) value denoting

# External links
$(_doc_external("MPI_Query_thread"))
"""
function Query_thread()
    r_provided = Ref{ThreadLevel}()

    # int MPI_Query_thread(int *provided)
    @mpichk ccall((:MPI_Query_thread, libmpi), Cint,
                  (Ref{ThreadLevel},), r_provided)
    return r_provided[]
end

"""
    Is_thread_main()

Queries whether the current thread is the main thread according to MPI. This can be called
by any thread, and is useful for the  `THREAD_FUNNELED` [`ThreadLevel`](@ref).

# External links
$(_doc_external("MPI_Is_thread_main"))
"""
function Is_thread_main()
    r_flag = Ref{Cint}()
    # int MPI_Is_thread_main(int *flag)
    @mpichk ccall((:MPI_Is_thread_main, libmpi), Cint,
                  (Ref{Cint},), r_flag)
    return r_flag[] != 0
end

"""
    Finalize()

Marks MPI state for cleanup. This should be called after [`Init`](@ref), at most once, and
no further MPI calls (other than [`Initialized`](@ref) or [`Finalized`](@ref)) should be
made after it is called.

Note that this does not correspond exactly to `MPI_FINALIZE` in the MPI specification. In
particular:

- It may not finalize MPI immediately. Julia will wait until all MPI-related objects are
  garbage collected before finalizing MPI. As a result, [`Finalized()`](@ref) may return
  `false` after `Finalize()` has been called. See [Finalizers](@ref) for more details.

- It is optional: [`Init`](@ref) will automatically insert a hook to finalize MPI when
  Julia exits.

# External links
$(_doc_external("MPI_Finalize"))
"""
function Finalize()
    # calling atexit here is a bit silly, but it's to avoid a case where MPI is finalized
    # one object early, e.g.
    #
    # event         | REFCOUNT
    # ---------------------
    # Init()        |   1  : MPI_INIT
    # new object    |   2  : MPI_X_CREATE
    # Finalize()    |   1
    # atexit        |
    #  refcount_inc |   2  : relies on LIFO ordering
    #  refcount_dec |   1  : MPI_FINALIZE would otherwise be called here
    # finalizers    |
    #  refcount_dec |   0  : MPI_X_FREE, MPI_FINALIZE
    atexit(refcount_inc)
    refcount_dec()
end

function _Finalize()
    @mpichk ccall((:MPI_Finalize, libmpi), Cint, ())
end

"""
    Abort(comm::Comm, errcode::Integer)

Make a “best attempt” to abort all tasks in the group of `comm`. This function does not
require that the invoking environment take any action with the error code. However, a Unix
or POSIX environment should handle this as a return errorcode from the main program.

# External links
$(_doc_external("MPI_Abort"))
"""
function Abort(comm::Comm, errcode::Integer)
    @mpichk ccall((:MPI_Abort, libmpi), Cint, (MPI_Comm, Cint), comm, errcode)
end

"""
    Initialized()

Returns `true` if [`MPI.Init`](@ref) has been called, `false` otherwise.

It is unaffected by [`MPI.Finalize`](@ref), and is one of the few functions that may be
called before [`MPI.Init`](@ref).

# External links
$(_doc_external("MPI_Intialized"))
"""
function Initialized()
    flag = Ref{Cint}()
    @mpichk ccall((:MPI_Initialized, libmpi), Cint, (Ptr{Cint},), flag)
    flag[] != 0
end

"""
    Finalized()

Returns `true` if [`MPI.Finalize`](@ref) has completed, `false` otherwise.

It is safe to call before [`MPI.Init`](@ref) and after [`MPI.Finalize`](@ref).

# External links
$(_doc_external("MPI_Finalized"))
"""
function Finalized()
    flag = Ref{Cint}()
    @mpichk ccall((:MPI_Finalized, libmpi), Cint, (Ptr{Cint},), flag)
    flag[] != 0
end

function Wtick()
    @mpicall ccall((:MPI_Wtick, libmpi), Cdouble, ())
end

function Wtime()
    @mpicall ccall((:MPI_Wtime, libmpi), Cdouble, ())
end

"""
    MPI.has_cuda()

Check if the MPI implementation is known to have CUDA support. Currently only Open MPI
provides a mechanism to check, so it will return `false` with other implementations
(unless overriden).

This can be overriden by setting the `JULIA_MPI_HAS_CUDA` environment variable to `true`
or `false`.
"""
function has_cuda()
    flag = get(ENV, "JULIA_MPI_HAS_CUDA", nothing)
    if flag === nothing
        # Only Open MPI provides a function to check CUDA support
        @static if MPI_LIBRARY == OpenMPI
            # int MPIX_Query_cuda_support(void)
            return 0 != ccall((:MPIX_Query_cuda_support, libmpi), Cint, ())
        elseif MPI_LIBRARY == IBMSpectrumMPI
            return true
        else
            return false
        end
    else
        return parse(Bool, flag)
    end
end


text:
using ..Statistics
using ..Bases

typealias GPCDistribution ContinuousUnivariateDistribution

"A gpc pair, capturing the pairing of distributions with systems of orthogonal polymials"
immutable GPCPair
  dist::GPCDistribution
  polysys::PolynomialSystem
end

"The standard mapping from distribution to orthogonal polynomials"
map_dist_to_polysys = Dict{GPCDistribution, PolynomialSystem}()
map_dist_to_polysys[Uniform(-1,1)] = LegendrePolynomials()
map_dist_to_polysys[Normal(0,1)] = HermitePolynomials()

"Constructor for GPCPairs "
GPCPair(dist) = GPCPair(dist, map_dist_to_polysys[dist])

"The abstract base for gpc germs (maybe we extend that later to infinite germs)"
abstract GPCGermBase

"A gpc germ of fixed length (the standard case)"
immutable GPCGerm <: GPCGermBase
  pairs::Array{GPCPair}
end
function GPCGerm(dist::GPCDistribution, n::Int)
  GPCGerm(fill(GPCPair(dist), n))
end

import Base.length
length(germ::GPCGerm) = length(germ.pairs)

distributions(germ::GPCGerm) = [pair.dist for pair in germ.pairs]

import Base.rand
rand(germ::GPCGerm) = Float64[rand(dist) for dist=distributions(germ)]

import Distributions.pdf
function pdf{T<:Real}(germ::GPCGerm, x::Vector{T})
  prod([pdf(dist, xi) for (dist,xi) in zip(distributions(germ), x)])
end

import Distributions.cdf
function cdf{T<:Real}(germ::GPCGerm, x::Vector{T})
  prod([cdf(dist, xi) for (dist,xi) in zip(distributions(germ), x)])
end

export GPCGerm,
       GPCPair


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models: anomaly_score, generalized_anomaly_score_gpu
using BSON
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		required = true
		arg_type = Int
		help = "seed"
	"dataset"
		required = true
		arg_type = String
		help = "dataset"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

modelname ="DeepSVDD"

function sample_params()
	argnames = (
		:hdim,
		:zdim,
		:nlayers,
        :activation,
        :objective,
        :nu,
        :lr_svdd,
        :lr_ae,
		:batch_size,
		:init_seed,
	)
	par_vec = (
		2 .^(4:9),
		2 .^(3:8),
		3:4,
        ["relu", "swish", "tanh"],
        ["soft-boundary", "one-class"],
        [0.01f0, 0.1f0, 0.5f0, 0.99f0], # paper 0.1
        10f0 .^ (-4:-3),
        10f0 .^ (-4:-3),
		2 .^ (5:7),
		1:Int(1e8),
	)
	return NamedTuple{argnames}(map(x->sample(x,1)[1], par_vec))
end

function fit(data, parameters)
	all_parameters = merge(
		parameters,
		(
			iters=5000,
			check_every=30,
			ae_check_every = 30,
			patience=10,
			ae_iters=5000
		)
	)

    svdd = GenerativeAD.Models.ae_constructor(;all_parameters...) |> gpu

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(svdd, data, all_parameters)
	catch e
		println("Error caught => $(e).")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

    training_info = (
		fit_t = fit_t,
		model = info[2]|>cpu,
		history_svdd = info[1][1], # losses through time
		history_ae = info[1][2],
		npars = info[3], # number of parameters
		iters = info[4] # optim iterations of model
		)

    return training_info, [(x -> GenerativeAD.Models.anomaly_score_gpu(svdd, x), parameters)]
    # if there is no gpu on pc anomaly_score will automaticly run on cpu
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
    parameters = sample_params()

    for seed in 1:max_seed
		savepath = datadir("experiments/tabular$cont_string/$(modelname)/$(dataset)/seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

        @info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
		@info "Train/validation/test splits: $(size(data[1][1], 2)) | $(size(data[2][1], 2)) | $(size(data[3][1], 2))"
		@info "Number of features: $(size(data[1][1])[1])"

		# update parameter
		parameters = merge(parameters, (idim=size(data[1][1],1), ))
		# here, check if a model with the same parameters was already tested
		if GenerativeAD.check_params(savepath, parameters)
			#(X_train,_), (X_val, y_val), (X_test, y_test) = data
			training_info, results = fit(data, parameters)
			# saving model separately
			if training_info.model !== nothing
				tagsave(joinpath(savepath, savename("model", parameters, "bson")), Dict("model"=>training_info.model), safe = true)
				training_info = merge(training_info, (model = nothing,))
			end
			save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))
			# now loop over all anomaly score funs
			for result in results
				GenerativeAD.experiment(result..., data, savepath; save_entries...)
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, sampling new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
module ProfileSVG

using FlameGraphs, Colors, UUIDs, Profile
using Base.StackTraces: StackFrame
const FlameGraph = FlameGraphs.Node{FlameGraphs.NodeData}

export @profview

"""
    @profview f(args...)

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.
"""
macro profview(ex)
    return quote
        Profile.clear()
        @profile $(esc(ex))
        view()
    end
end

function flamegraph_kwargs(kwargs)
    gopts = copy(default_config.graph_options)

    keywords = (:lidict, :C, :combine, :recur, :norepl, :pruned, :filter)
    for (k, v) in kwargs
        if k in keywords
            gopts[k] = v
        end
    end
    gopts
end

struct FGConfig
    g::Union{FlameGraph, Nothing}
    graph_options::Dict{Symbol, Any}
    fcolor
    bgcolor::Symbol
    fontcolor::Symbol
    frameopacity::Float64
    yflip::Bool
    maxdepth::Int
    maxframes::Int
    width::Float64
    height::Float64
    roundradius::Float64
    font::String
    fontsize::Float64
    notext::Bool
    timeunit::Symbol
    delay::Float64
end

function FGConfig(g::Union{FlameGraph, Nothing} = nothing,
                  graph_options = nothing,
                  fcolor               = default_config.fcolor;
                  bgcolor::Symbol      = default_config.bgcolor,
                  fontcolor::Symbol    = default_config.fontcolor,
                  frameopacity::Real   = default_config.frameopacity,
                  yflip::Bool          = default_config.yflip,
                  maxdepth::Int        = default_config.maxdepth,
                  maxframes::Int       = default_config.maxframes,
                  width::Real          = default_config.width,
                  height::Real         = default_config.height,
                  roundradius::Real    = default_config.roundradius,
                  font::AbstractString = default_config.font,
                  fontsize::Real       = default_config.fontsize,
                  notext::Bool         = default_config.notext,
                  timeunit::Symbol     = default_config.timeunit,
                  delay::Real          = default_config.delay,
                  kwargs...)

    gopts = graph_options === nothing ? flamegraph_kwargs(kwargs) : graph_options
    delay = g === nothing || delay > 0 ? delay : last(Profile.init())
    FGConfig(g, gopts, fcolor,
             bgcolor, fontcolor, frameopacity,
             yflip, maxdepth, maxframes, width, height, roundradius,
             font, fontsize, notext, timeunit, delay)
end

"""
    ProfileSVG.init()

Initialize the settings.
"""
function init()
    global default_config = FGConfig(nothing,
                                     Dict{Symbol, Any}(),
                                     FlameGraphs.default_colors,
                                     bgcolor=:fcolor,
                                     fontcolor=:fcolor,
                                     frameopacity=1,
                                     yflip=false,
                                     maxdepth=50,
                                     maxframes=2000,
                                     width=960,
                                     height=0,
                                     roundradius=2,
                                     font="inherit",
                                     fontsize=12,
                                     notext=false,
                                     timeunit=:none,
                                     delay=0.0)
    nothing
end

"""
    ProfileSVG.set_default([fcolor]; kwargs...)

Set defult configurations of profiling results. See [`ProfileSVG.view`](@ref)
for details of `kwargs`.

# Examples
```julia
ProfileSVG.set_default(width=600, fontsize=9)
```
"""
function set_default(fcolor=default_config.fcolor; kwargs...)
    global default_config = FGConfig(nothing, nothing, fcolor; kwargs...)
end

include("svgwriter.jl")

"""
    ProfileSVG.view([fcolor,] data=Profile.fetch(); kwargs...)
    ProfileSVG.view([fcolor,] g::FlameGraph; kwargs...)

View profiling results.

# keywords for SVG style
- `bgcolor` (default: `:fcolor`)
  - The background style. One of `:fcolor`/`:classic`/`:transparent`.
- `fontcolor` (default: `:fcolor`)
  - The font color style. One of `:fcolor`/`:classic`/`:currentcolor`/`:bw`.
- `frameopacity` (default: `1`)
  - The opacity of frames in [0, 1].
- `yflip` (default: `false`)
  - If `true`, the "icicle" graph will be rendered.
- `maxdepth` (default: `50`)
  - The maximum number of the rendered rows.
- `maxframes` (default: `2000`)
  - The maximum number of the rendered frames.
- `width` (default: `960`)
  - The width of output SVG image in pixels.
- `height` (default: `0`)
  - The height of output SVG image in pixels. If you set the height to `0`, it
    will be calculated automatically according to the graph.
- `roundradius` (default: `2`)
  - The rounding radius of the corners of each frame in pixels.
- `font` (default: `"inherit"`)
  - The font family names for texts. This setting is used as the CSS
    `font-family` property, i.e. you can use a comma-separated list.
- `fontsize` (default: `12`)
  - The font size of texts for function information, in pixels (not points).
- `notext` (default: `false`)
  - If `true`, the texts overlaid on the frames will be hidden by the
    interactive feature.
- `timeunit` (default: `:none`)
  - If `:s`, `:ms`, `:us`, or `:µs` is specified, the duration of the block will
    be displayed in the lower right corner in the specified unit.
- `delay` (default: `0.0`)
  - The delay between backtraces, in seconds. If a non-positive number is
    specified, the current setting in `Profile.init()` will be used.

# keywords for `flamegraph`
- `lidict`
- `C`
- `combine`
- `recur`
- `pruned`
- `filter`
See [FlameGraphs](https://timholy.github.io/FlameGraphs.jl/stable/reference/#FlameGraphs.flamegraph)
for details.
"""
function view(fcolor, data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    gopts = flamegraph_kwargs(kwargs)
    g = flamegraph(data; gopts...)
    FGConfig(g, gopts, fcolor; kwargs...)
end
function view(data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    gopts = flamegraph_kwargs(kwargs)
    g = flamegraph(data; gopts...)
    FGConfig(g, gopts; kwargs...)
end
function view(fcolor, g::FlameGraph; kwargs...)
    FGConfig(g, nothing, fcolor; kwargs...)
end
function view(g::FlameGraph; kwargs...)
    FGConfig(g; kwargs...)
end

"""
    ProfileSVG.save([fcolor,] io::IO, data=Profile.fetch(); kwargs...)
    ProfileSVG.save([fcolor,] io::IO, g::FlameGraph; kwargs...)
    ProfileSVG.save([fcolor,] filename, data=Profile.fetch(); kwargs...)
    ProfileSVG.save([fcolor,] filename, g::FlameGraph; kwargs...)

Save profile results as an SVG file. See [`ProfileSVG.view`](@ref) for details
of `kwargs`.
"""
function save(fcolor, io::IO, data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    gopts = flamegraph_kwargs(kwargs)
    g = flamegraph(data; gopts...)
    show(io, MIME"image/svg+xml"(), FGConfig(g, gopts, fcolor; kwargs...))
end
function save(fcolor, io::IO, g::FlameGraph; kwargs...)
    show(io, MIME"image/svg+xml"(), FGConfig(g, nothing, fcolor; kwargs...))
end
function save(io::IO, data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    gopts = flamegraph_kwargs(kwargs)
    g = flamegraph(data; gopts...)
    show(io, MIME"image/svg+xml"(), FGConfig(g, gopts; kwargs...))
end
function save(io::IO, g::FlameGraph; kwargs...)
    show(io, MIME"image/svg+xml"(), FGConfig(g; kwargs...))
end

function save(fcolor, filename::AbstractString, data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    open(filename, "w") do file
        save(fcolor, file, data; kwargs...)
    end
    return nothing
end
function save(fcolor, filename::AbstractString, g::FlameGraph; kwargs...)
    open(filename, "w") do file
        save(fcolor, file, g; kwargs...)
    end
    return nothing
end
function save(filename::AbstractString, data::Vector{<:Unsigned}=Profile.fetch(); kwargs...)
    open(filename, "w") do file
        save(file, data; kwargs...)
    end
    return nothing
end
function save(filename::AbstractString, g::FlameGraph; kwargs...)
    open(filename, "w") do file
        save(file, g; kwargs...)
    end
    return nothing
end

Base.showable(::MIME"image/svg+xml", fg::FGConfig) = fg.g !== nothing

# If html embedding does not work, redefine the following method to return `false`.
Base.showable(::MIME"text/html", fg::FGConfig) = fg.g !== nothing

function Base.show(io::IO, mime::Union{MIME"image/svg+xml", MIME"text/html"}, fg::FGConfig)
    if mime isa MIME"text/html"
        print(io, """
            <!DOCTYPE html>
            <html>
            <body>
            """)
    else
        write_svgdeclaration(io)
    end

    show_flamegraph_body(io, fg)

    if mime isa MIME"text/html"
        print(io, """
            </body>
            </html>
            """)
    end
end

function extract_frameinfo(sf::StackFrame)
    file = string(sf.file)
    m = match(r"[^\\/]+$", file)
    if m !== nothing
        dirinfo = SubString(file, firstindex(file), m.offset - 1)
        basename = m.match
    else
        dirinfo = ""
        basename = file
    end
    shortinfo = "$(sf.func) in $basename:$(sf.line)"
    return shortinfo, dirinfo
end

function show_flamegraph_body(io::IO, fg::FGConfig)
    ncols, nrows = length(fg.g.data.span), FlameGraphs.depth(fg.g)
    if nrows > fg.maxdepth
        @warn """The depth of this graph is $nrows, exceeding the `maxdepth` (=$(fg.maxdepth)).
                 The deeper frames will be truncated."""
        nrows = fg.maxdepth
    end
    width = fg.width
    leftmargin = rightmargin = round(Int, width * 0.01)
    topmargin = botmargin = round(Int, max(width * 0.04, fg.fontsize * 3))

    idealwidth = width - (leftmargin + rightmargin)
    xstep = Float64(rationalize(idealwidth / ncols, tol = 1 / ncols))
    ystep = round(Int, fg.fontsize * 1.25)

    height = fg.height > 0.0 ? fg.height : ystep * nrows + botmargin * 2.0

    function flamerects(io::IO, g::FlameGraph, j::Int, nextidx::Vector{Int})
        j > fg.maxdepth && return
        nextidx[end] > fg.maxframes && return
        nextidx[end] += 1

        ndata = g.data
        color = fg.fcolor(nextidx, j, ndata)::Color
        bw = fg.fontcolor === :bw
        x = (first(ndata.span)-1) * xstep + leftmargin
        if fg.yflip
            y = topmargin + (j - 1) * ystep
        else
            y = height - j * ystep - botmargin
        end
        w = length(ndata.span) * xstep
        r = fg.roundradius
        shortinfo, dirinfo = extract_frameinfo(ndata.sf)
        write_svgflamerect(io, x, y, w, ystep, r, shortinfo, dirinfo, color, bw)

        for c in g
            flamerects(io, c, j + 1, nextidx)
        end
    end

    fig_id = string("fig-", replace(string(uuid4()), "-" => ""))

    write_svgheader(io, fig_id, width, height,
                    bgcolor(fg), fontcolor(fg), fg.frameopacity,
                    fg.font, fg.fontsize, fg.notext, xstep, fg.timeunit, fg.delay)

    nextidx = fill(1, nrows + 1) # nextidx[end]: framecount
    flamerects(io, fg.g, 1, nextidx)

    if nextidx[end] > fg.maxframes
        @warn """The maximum number of frames (`maxframes`=$(fg.maxframes)) is reached.
                 Some frames were truncated."""
    end

    write_svgfooter(io, fig_id)
end

function bgcolor(fg)
    fg.bgcolor === :fcolor && return "#" * hex(fg.fcolor(:bg))
    fg.bgcolor === :transparent && return "transparent"
    fg.bgcolor === :classic && return ""
    return "white"
end

function fontcolor(fg)
    fg.fontcolor === :fcolor && return "#" * hex(fg.fcolor(:font))
    fg.fontcolor === :currentcolor && return "currentcolor"
    fg.fontcolor === :bw && return ""
    return "black"
end

__init__() = init()

end # module


text:
## Load
@adddevice(Load, OnePort, LinearDispatchForm, :loads, (:L,), joinpath(@__DIR__, "attrs", "loads.csv"))

cost(d::Load) = 0.
function addtoexpr!(injection::AxisArray, ::NodalActivePower,
                    m::EnergyModel, d::Load)
    B, T = AxisArrays.axes(injection)
    D = axis(m, d)

    buses = get(d, :bus, D)
    p_set = get(d, :p_set, D, T)

    for (i, idx) in enumerate(indexin(buses, B)), t in T
        !isnothing(idx) && add_to_expression!(injection[idx,t], -p_set[i,t])
    end
end
addto!(jm::JuMP.AbstractModel, m::EnergyModel, d::Load) = nothing

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
module VisdomLog
export Visdom, report
using PyCall
using PyPlot: Figure

struct Visdom
  vis::PyObject
  d::Dict{Symbol, Pair{PyObject,Any}}
end

"Start a new visdom environment"
function Visdom(env::String)
  visdom = pyimport(:visdom)
  vd = visdom[:Visdom](env=env)
  Visdom(vd, Dict{Symbol, Pair{PyObject,Any}}())
end

# ABSTRACT this pattern!

"Add a value to a log of results"
function report(l, k::Symbol, y::Float64)
  if !haskey(l.d, k)
    win = l.vis[:line](X=[1], Y=[y], opts=Dict(:title=>String(k)))
    l.d[k] = (win => 2);
  else
    (win, x) = l.d[k]
    l.vis[:line](win=win, X=[x], Y=[y], update="append", opts=Dict(:ytype=>"log"))
    l.d[k] = (win => x + 1);
  end
end

"Display a histogram of current results"
function report(l, k::Symbol, vals::Array{Float64})
  if !haskey(l.d, k)
    win = l.vis[:histogram](X=vals[:], opts=Dict(:title=>String(k)))
    l.d[k] = win=>nothing;
  else
    (win, _) = l.d[k]
    l.vis[:histogram](win=win, X=vals[:], opts=Dict(:ytype=>"log", :title=>String(k)));
  end
end

"Plot a comparison of current results"
function report(l, k::Symbol, (xs, ys, legend)::Tuple{Vector,Matrix{Float64},Vector})
  if !haskey(l.d, k)
    win = l.vis[:line](X=xs, Y=ys, opts=Dict(:legend=>legend, :title=>String(k)))
    l.d[k] = win=>nothing;
  else
    (win, _) = l.d[k]
    l.vis[:line](win=win, X=xs, Y=ys, update="replace",
      opts=Dict(:legend=>legend, :title=>String(k)));
  end
end

"Plot a comparison of current results"
function report(l, k::Symbol, (xs, ys)::Tuple{Vector,Vector{Float64}})
  if !haskey(l.d, k)
    win = l.vis[:line](X=xs, Y=ys, opts=Dict(:title=>String(k)))
    l.d[k] = win=>nothing;
  else
    (win, _) = l.d[k]
    l.vis[:line](win=win, X=xs, Y=ys, update="replace",
      opts=Dict(:title=>String(k)));
  end
end

"Plot a matplotlib figure"
function report(l, k::Symbol, fig::Figure)
  if !haskey(l.d, k)
    win = l.vis[:matplot](plot=fig, opts=Dict(:title=>String(k)))
    l.d[k] = win=>nothing;
  else
    (win, _) = l.d[k]
    l.vis[:matplot](win=win, plot=fig, opts=Dict(:title=>String(k)));
  end
end

end

max_stars_repo_name: bogiebro/VisdomLog

max_issues_repo_name: bogiebro/VisdomLog

max_forks_repo_name: bogiebro/VisdomLog


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 230fdd70-c5ef-11eb-01d9-371bb366387d
begin
	using CSV
	using DataFrames
	using Plots
	using GLM
	using PlutoUI
end

# ╔═╡ 2d3fad52-f1ac-4383-ad78-782c5265148d
data = CSV.read("D:\\Workspace\\Machine Learning\\World GDP Analysis and Prediction\\Data\\World_GDP.csv",DataFrame)

# ╔═╡ 1a4f31a3-5b58-4f4f-8513-a416fd795ad5
scatter(data.year,data.GDP)

# ╔═╡ 79e34de2-59b5-4bdb-b16e-93c64d2010a2
reg = fit(LinearModel,@formula(GDP~year),data)

# ╔═╡ 35d8c82f-d32c-4c2c-8934-280b39761320
linear_regression = predict(reg)

# ╔═╡ 3d4580cb-eb49-4979-9063-60a03c4072b2
begin
	scatter(data.year,data.GDP)
	plot!(data.year,linear_regression)
end

# ╔═╡ 7b0a2af1-148f-42f7-87f6-3a237c00f8cf
md"Enter an year to predict : $@bind y TextField()"

# ╔═╡ 342439a5-e524-40a3-85a7-43a78cf07596
pyear = parse(Int64,y)

# ╔═╡ 75394533-5d15-47c2-9a34-5d98ffd25900
pGDP = coef(reg)[2] * pyear + coef(reg)[1]

# ╔═╡ Cell order:
# ╠═230fdd70-c5ef-11eb-01d9-371bb366387d
# ╠═2d3fad52-f1ac-4383-ad78-782c5265148d
# ╠═1a4f31a3-5b58-4f4f-8513-a416fd795ad5
# ╠═79e34de2-59b5-4bdb-b16e-93c64d2010a2
# ╠═35d8c82f-d32c-4c2c-8934-280b39761320
# ╠═3d4580cb-eb49-4979-9063-60a03c4072b2
# ╠═7b0a2af1-148f-42f7-87f6-3a237c00f8cf
# ╠═342439a5-e524-40a3-85a7-43a78cf07596
# ╠═75394533-5d15-47c2-9a34-5d98ffd25900

max_stars_repo_path: World GDP Analysis and Prediction/main.jl

max_stars_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_issues_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_forks_repo_name: chidhvilash/GDP-PerCapita-Analysis


text:
module Estimator
include("SimplePF.jl")
import .SimplePF
using JuMP

export nl_pf, add_objective

function build_sets(sys)
    sub_der = [Estimator.build_DER_set_buses(sys); sys.substation.bus]
    Ω = sys.buses
    bΩ = filter(bus -> bus ∉ sub_der, Ω)
    L = collect(1:length(sys.m_load))
    K = collect(1:reduce(*, length(der.scenario) for der in sys.dgs))
    D = collect(1:length(sys.dgs))
    S = collect(1:length(sys.m_new_dg))
    return (Ω, bΩ, L, K, D, S)
end

function build_DER_set_buses(sys)
    return [dg.bus for dg in sys.dgs]
end

function build_DER_scenario(sys)
    qnt_scenarios = reverse([length(der.scenario) for der in sys.dgs])
    values = zeros(Int, length(sys.dgs))
    μᴰᴱᴿ(k::Int, d::Int) = begin
        k_n = k - 1 # To turn value 0 as a frist index
        for (indx, qnt_scenario) in enumerate(qnt_scenarios)
            values[indx] = (k_n % qnt_scenario) + 1 # To back to value 1 as a first index
            k_n = div(k_n, qnt_scenario)
        end
        indx = reverse(values)[d]
        return sys.dgs[d].scenario[indx]
    end
    return μᴰᴱᴿ
end

function add_variables(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)

    @variable(model, -sys.VH ≤ V[[:Re, :Im], Ω, L, K, S] ≤ sys.VH)
    for b = Ω, l = L, k = K, s = S
        set_start_value.(V[:Re, b, l, k, s], 1.0)
        set_start_value.(V[:Im, b, l, k, s], 0.0)
    end
    @variable(model, sys.dgs[d].P_limit[1] ≤ pᴰᴱᴿ[d=D, L, K, S] ≤ sys.dgs[d].P_limit[2], start = 0.0)
    @variable(model, sys.dgs[d].Q_limit[1] ≤ qᴰᴱᴿ[d=D, L, K, S] ≤ sys.dgs[d].Q_limit[2], start = 0.0)
    @variable(model, 0.0 ≤ pᴴᶜ, start = 0.0)

    return model
end

function add_voltage_constraints(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    V = model[:V]

    @expression(model, V²[b=Ω, l=L, k=K, s=S],
        V[:Re, b, l, k, s]^2 + V[:Im, b, l, k, s]^2
    )

    @constraint(model, voltage_constraint[b=Ω, l=L, k=K, s=S],
        sys.VL^2 <= V²[b, l, k, s] <= sys.VH^2
    )

    @NLexpression(model, V_module[b=Ω, l=L, k=K, s=S],
        sqrt(V²[b, l, k, s])
    )

    return model
end

function add_I_V_relationship(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    G = real(sys.Y)
    B = imag(sys.Y)
    V = model[:V]
    @expression(model, I[z=[:Re, :Im], i=Ω, l=L, k=K, s=S],
        if z == :Re
            sum(
                G[i, j] * V[:Re, j, l, k, s] - B[i, j] * V[:Im, j, l, k, s]
                for j in Ω
            )
        elseif z == :Im
            sum(
                B[i, j] * V[:Re, j, l, k, s] + G[i, j] * V[:Im, j, l, k, s]
                for j in Ω
            )
        end
    )

    return model
end

function add_S_VI_relationship(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    V = model[:V]
    I = model[:I]

    @expression(model, P[b=Ω, l=L, k=K, s=S],
        SimplePF.mc_re(V[:Re, b, l, k, s], V[:Im, b, l, k, s], I[:Re, b, l, k, s], -I[:Im, b, l, k, s])
    )
    @expression(model, Q[b=Ω, l=L, k=K, s=S],
        SimplePF.mc_im(V[:Re, b, l, k, s], V[:Im, b, l, k, s], I[:Re, b, l, k, s], -I[:Im, b, l, k, s])
    )

    return model
end

function add_active_losses(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    P = model[:P]
    @expression(model, Ploss[l=L, k=K, s=S],
        sum(P[b, l, k, s] for b = Ω)
    )
    return model
end

function add_reactive_losses(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    Q = model[:Q]
    @expression(model, Qloss[l=L, k=K, s=S],
        sum(Q[b, l, k, s] for b = Ω)
    )
    return model
end

function add_substation_constraint(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    V = model[:V]
    I = model[:I]
    P = model[:P]
    Q = model[:Q]
    sub = sys.substation
    for l = L, k = K, s = S
        fix(V[:Im, sub.bus, l, k, s], 0.0, force=true)
        set_start_value(V[:Re, sub.bus, l, k, s], sub.voltage)
        set_lower_bound(V[:Re, sub.bus, l, k, s], 0.9)
        set_upper_bound(V[:Re, sub.bus, l, k, s], 1.1)
    end

    @constraint(model, sub_plimit[l=L, k=K, s=S],
        0 <= P[sub.bus, l, k, s] <= sub.P_limit
    )
    @constraint(model, sub_qlimit[l=L, k=K, s=S],
        0 <= Q[sub.bus, l, k, s] <= sub.Q_limit
    )

    @constraint(model, sub_current[l=L, k=K, s=S],
        I[:Re, sub.bus, l, k, s]^2 + I[:Im, sub.bus, l, k, s]^2 >= 0.0
    )

    return model
end

function add_line_constraint(model, sys)
    G = real(sys.Y)
    B = imag(sys.Y)

    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    F = sys.amp
    Ωᴮ = keys(F)

    V = model[:V]

    vij(z, i, j, l, k, s) = V[z, i, l, k, s] - V[z, j, l, k, s]
    @expression(model, f[z=[:Re, :Im], (i, j)=Ωᴮ, l=L, k=K, s=S],
        if z == :Re
            SimplePF.mc_re(vij(:Re, i, j, l, k, s), vij(:Im, i, j, l, k, s), G[i, j], B[i, j])
        elseif z == :Im
            SimplePF.mc_im(vij(:Re, i, j, l, k, s), vij(:Im, i, j, l, k, s), G[i, j], B[i, j])
        end
    )

    @constraint(model, line_limit[(i, j)=Ωᴮ, l=L, k=K, s=S],
        f[:Re, (i, j), l, k, s]^2 + f[:Im, (i, j), l, k, s]^2 ≤ (F[(i, j)])^2
    )

    return model
end

function add_power_injection_definition(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    B = Estimator.build_DER_set_buses(sys)
    PL = sys.PL
    QL = sys.QL
    P = model[:P]
    Q = model[:Q]
    pᴴᶜ = model[:pᴴᶜ]
    μᴸ = sys.m_load
    μᴴᶜ = sys.m_new_dg
    @constraint(model, p[b=bΩ, l=L, k=K, s=S],
        P[b, l, k, s] == μᴴᶜ[s] * pᴴᶜ - μᴸ[l] * PL[b])
    @constraint(model, q[b=bΩ, l=L, k=K, s=S],
        Q[b, l, k, s] == -μᴸ[l] * QL[b])

    αᴰᴱᴿ(d) = sys.dgs[d].alpha
    Sᴰᴱᴿ(d) = sys.dgs[d].S_limit
    Pᴰᴱᴿ(d) = (1.0 - αᴰᴱᴿ(d)) * Sᴰᴱᴿ(d)
    pᴰᴱᴿ = model[:pᴰᴱᴿ]
    qᴰᴱᴿ = model[:qᴰᴱᴿ]
    μᴰᴱᴿ = build_DER_scenario(sys)
    @constraint(model, p_wder[d=D, l=L, k=K, s=S],
        P[B[d], l, k, s] == μᴰᴱᴿ(k, d) * Pᴰᴱᴿ(d) + pᴰᴱᴿ[d, l, k, s] - μᴸ[l] * PL[B[d]])
    @constraint(model, q_wder[d=D, l=L, k=K, s=S],
        Q[B[d], l, k, s] == qᴰᴱᴿ[d, l, k, s] - μᴸ[l] * QL[B[d]])

    return model
end

function add_ders_limits(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)

    αᴰᴱᴿ(d) = sys.dgs[d].alpha
    Sᴰᴱᴿ(d) = sys.dgs[d].S_limit
    Pᴰᴱᴿ(d) = (1.0 - αᴰᴱᴿ(d)) * Sᴰᴱᴿ(d)
    μᴰᴱᴿ = build_DER_scenario(sys)
    pᴰᴱᴿ = model[:pᴰᴱᴿ]
    qᴰᴱᴿ = model[:qᴰᴱᴿ]
    Eᴰᴱᴿ(d) = sys.dgs[d].energy
    βᴰᴱᴿ(d) = sys.dgs[d].beta
    Tᶜᵘʳ = sys.time_curr
    @constraint(model, disco_der_limit[d=D, l=L, k=K, s=S],
        pᴰᴱᴿ[d, l, k, s]^2 + qᴰᴱᴿ[d, l, k, s]^2 ≤ (αᴰᴱᴿ(d) * Sᴰᴱᴿ(d))^2
    )

    @constraint(model, der_limit[d=D, l=L, k=K, s=S],
        (μᴰᴱᴿ(k, d) * Pᴰᴱᴿ(d) + pᴰᴱᴿ[d, l, k, s])^2 + qᴰᴱᴿ[d, l, k, s]^2 ≤ Sᴰᴱᴿ(d)^2
    )

    @constraint(model, energy_limit[d=D, l=L, k=K, s=S],
        -βᴰᴱᴿ(d) * Eᴰᴱᴿ(d) ≤ Tᶜᵘʳ * pᴰᴱᴿ[d, l, k, s] ≤ βᴰᴱᴿ(d) * Eᴰᴱᴿ(d)
    )

    return model
end

function nl_pf(model, sys)

    model = add_variables(model, sys)
    model = add_voltage_constraints(model, sys)
    model = add_I_V_relationship(model, sys)
    model = add_S_VI_relationship(model, sys)
    model = add_substation_constraint(model, sys)
    model = add_line_constraint(model, sys)
    model = add_power_injection_definition(model, sys)
    model = add_ders_limits(model, sys)

    return model
end

function add_objective(model, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    Nᴮ = length(bΩ)
    pᴴᶜ = model[:pᴴᶜ]
    set_objective(model, MAX_SENSE, pᴴᶜ * Nᴮ)

    return model
end

end


text:
# multi_output.jl
module MultiOutput
export foo!

# inputs: a, b; outputs: c, d
function foo!(a, b, c, d)
    c .= a + b
    d .= a - b
end

end


text:
using Abaco
using Test

Abaco.loginit()

# My time zone: Saturday, October 30, 2021 3:43:35 PM GMT+02:00 DST
metric_ts = 1635601415
metric_sn = "Pelmo"

interval = 5
ages = 4

values = Dict(
    "x" => 2.0,
    "y" => 4.0,
    "z" => 8.0
    )

expected_triggers = 2
actual_triggers = 0

function onresult(ts, ne, name, value, inputs)
    global actual_triggers
    @debug "age [$ts]: scope: [$ne] $name = $value"
    @test ne == metric_sn
    @test value == 14.0
    actual_triggers += 1
end

abaco = abaco_init(onresult, interval=interval, ages=ages)

# add a formula that is not evaluated because variable t is not feeded
formula(abaco, "r = (x + y) / exp(t)")

# this gets evaluated (2 + 4 + 8 = 14)
formula(abaco, "w = x + y + z")

ingest(abaco, metric_ts, metric_sn, values)
ingest(abaco, metric_ts+interval*ages, metric_sn, values)

@test actual_triggers === expected_triggers


text:
T=3
tsname="joNFFT"
@testset "$tsname" begin
for t=1:T # start test loop
    m=4^t
    n=3^t
    nodes=sort(rand(joFloat,n))

    verbose && println("$tsname ($m[,$m]) - not centered")
    @testset "$m [x $m]" begin
        @test isadjoint(joNFFT(m,nodes))[1]
        @test islinear(joNFFT(m,nodes))[1]
    end

    verbose && println("$tsname ($m[,$m]) - centered")
    @testset "$m [x $m]" begin
        @test isadjoint(joNFFT(m,nodes,centered=true))[1]
        @test islinear(joNFFT(m,nodes,centered=true))[1]
    end

end # end test loop
end


text:
import BangBang
import OnlineStats
using Documenter
using Transducers
# Not importing `Literate` here so that this file can be included from
# test.

EXAMPLE_PAGES = [
    "Tutorial: Missing values" => "examples/tutorial_missings.md",
    "Empty result handling" => "examples/empty_result_handling.md",
    "Parallel word count" => "examples/words.md",
    "Writing transducers" => "examples/transducers.md",
    "Writing reducibles" => "examples/reducibles.md",
]

transducers_rm_examples() =
    rm(joinpath(@__DIR__, "src/examples"); force=true, recursive=true)

function transducers_literate(;
        inputbase = joinpath(@__DIR__, "..", "examples"),
        outputbase = joinpath(@__DIR__, "src", "examples"),
        examples = EXAMPLE_PAGES,
        kwargs...)
    transducers_rm_examples()
    for (_, outpath) in examples
        name, = splitext(basename(outpath))
        inputfile = joinpath(inputbase, "$name.jl")
        outputdir = outputbase
        Literate.markdown(
            inputfile, outputdir;
            documenter = true,
            kwargs...)
    end
end

function transducers_makedocs(;
        examples = EXAMPLE_PAGES,
        kwargs...)
    if isempty(examples)
        # Make some dummy examples
        dummypage = joinpath("examples", "dummy.md")
        path = joinpath(@__DIR__, "src", dummypage)
        mkpath(dirname(path))
        write(
            path,
            """
            # Dummy page
            ## Empty result handling
            """
        )
        examples = dummypage
    end
    makedocs(;
        modules = [Transducers],
        pages = [
            "Home" => "index.md",
            "Manual" => "manual.md",
            "Interface" => "interface.md",
            "Examples" => examples,
            hide("Internals" => "internals.md"),
            hide("doctests/show_xf.md"),
            hide("doctests/show_rf.md"),
        ],
        repo = "https://github.com/tkf/Transducers.jl/blob/{commit}{path}#L{line}",
        sitename = "Transducers.jl",
        authors = "Takafumi Arakaki",
        root = @__DIR__,
        strict = true,
        kwargs...)
end


text:
mmHgToPa = 133.32;
MPaToPa = 1e6;
cmTom = 1e-2;
cm3Tom3 = 1e-6;
minTos = 60;

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using HydrogenicAtoms
using HydrogenicAtoms.DiracAnalytic: E, α
using AtomicLevels: @ro_str
using Test

@testset "HydrogenicAtoms.jl" begin
    @test E(0, 1, -1) isa Real
    @test E(0, 1, -1) ≈ α^-2
    @test E(1/α, 1, -1) isa Real
    @test E(1/α, 1, -1) ≈ 0.0
    @test_throws DomainError E(138, 1, -1)

    @test E(Complex, 0, 1, -1) isa Complex
    @test E(Complex, 0, 1, -1) ≈ α^-2
    @test E(Complex, 1/α, 1, -1) isa Complex
    @test_broken E(Complex, 1/α, 1, -1) ≈ 0.0
    @test E(Complex, 138, 1, -1) isa Complex
    @test E(Complex, 138, 1, -1) ≈ 0.0 + 2231.3523399043934im

    @test E(0, ro"1s") isa Real
    @test E(0, ro"1s") ≈ α^-2
    @test E(Complex, 0, ro"1s") isa Complex
    @test E(Complex, 0, ro"1s") ≈ α^-2

    @test E(100, ro"2s") ≈ E(Complex, 100, 2, -1)
    @test E(100, ro"3p-") ≈ E(Complex, 100, 3, 1)
    @test E(100, ro"5p") ≈ E(Complex, 100, 5, -2)
end


text:
# ## Simple linear regression
#
# `MLJ` essentially serves as a unified path to many existing Julia packages each of which provides their own functionalities and models, with their own conventions.
#
# The simple linear regression demonstrates this.
# Several packages offer it (beyond just using the backslash operator): here we will use `MLJLinearModels` but we could also have used `GLM`, `ScikitLearn` etc.
#
# To load the model from a given package use `@load ModelName pkg=PackageName`

using MLJ
MLJ.color_off() # hide

@load LinearRegressor pkg=MLJLinearModels

# Note: in order to be able to load this, you **must** have the relevant package in your environment, if you don't, you can always add it (``using Pkg; Pkg.add("MLJLinearModels")``).
#
# Let's load the _boston_ data set

import RDatasets: dataset
import DataFrames: describe, select, Not, rename!
boston = dataset("MASS", "Boston")
first(boston, 3)

# Let's get a feel for the data

describe(boston, :mean, :std, :eltype)

# So there's no missing value and most variables are encoded as floating point numbers.
# In MLJ it's important to specify the interpretation of the features (should it be considered as a Continuous feature, as a Count, ...?), see also [this tutorial section](/getting-started/choosing-a-model/#data_and_its_interpretation) on scientific types.
#
# Here we will just interpret the integer features as continuous as we will just use a basic linear regression:

data = coerce(boston, autotype(boston, :discrete_to_continuous));

# Let's also extract the target variable (`MedV`):

y = data.MedV
X = select(data, Not(:MedV));

# Let's declare a simple multivariate linear regression model:

mdl = LinearRegressor()

# First let's do a very simple univariate regression, in order to fit it on the data, we need to wrap it in a _machine_ which, in MLJ, is the composition of a model and data to apply the model on:

X_uni = select(X, :LStat) # only a single feature
mach_uni = machine(mdl, X_uni, y)
fit!(mach_uni)

# You can then retrieve the  fitted parameters using `fitted_params`:

fp = fitted_params(mach_uni)
@show fp.coefs
@show fp.intercept

# You can also visualise this

using PyPlot
ioff() # hide

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
Xnew = (LStat = collect(range(extrema(X.LStat)..., length=100)),)
plot(Xnew.LStat, predict(mach_uni, Xnew))

savefig(joinpath(@OUTPUT, "ISL-lab-3-lm1.svg")) # hide

# \figalt{Univariate regression}{ISL-lab-3-lm1.svg}

# The  multivariate case is very similar

mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
coefs = fp.coefs
intercept = fp.intercept
for (name, val) in coefs
    println("$(rpad(name, 8)):  $(round(val, sigdigits=3))")
end
println("Intercept: $(round(intercept, sigdigits=3))")

# You can use the `machine` in order to _predict_ values as well and, for instance, compute the root mean squared error:

ŷ = predict(mach, X)
round(rms(ŷ, y), sigdigits=4)

# Let's see what the residuals look like

figure(figsize=(8,6))
res = ŷ .- y
stem(res)

savefig(joinpath(@OUTPUT, "ISL-lab-3-res.svg")) # hide

# \figalt{Plot of the residuals}{ISL-lab-3-res.svg}

# Maybe that a histogram is more appropriate here

figure(figsize=(8,6))
hist(res, density=true)

savefig(joinpath(@OUTPUT, "ISL-lab-3-res2.svg")) # hide

# \figalt{Histogram of the residuals}{ISL-lab-3-res2.svg}

# ## Interaction and transformation
#
# Let's say we want to also consider an interaction term of `lstat` and `age` taken together.
# To do this, just create a new dataframe with an additional column corresponding to the interaction term:

X2 = hcat(X, X.LStat .* X.Age);

# So here we have a DataFrame with one extra column corresponding to the elementwise products between `:LStat` and `Age`.
# DataFrame gives this a default name (`:x1`) which we can change:

rename!(X2, :x1 => :interaction);

# Ok cool, now let's try the linear regression again

mach = machine(mdl, X2, y)
fit!(mach)
ŷ = predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)

# We get slightly better results but nothing spectacular.
#
# Let's get back to the lab where they consider regressing the target variable on `lstat` and `lstat^2`; again, it's essentially a case of defining the right DataFrame:

X3 = hcat(X.LStat, X.LStat.^2)
mach = machine(mdl, X3, y)
fit!(mach)
ŷ = predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

# which again, we can visualise:

Xnew = (LStat = Xnew.LStat, LStat2 = Xnew.LStat.^2)

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
plot(Xnew.LStat, predict(mach, Xnew))

savefig(joinpath(@OUTPUT, "ISL-lab-3-lreg.svg")) # hide

# \figalt{Polynomial regression}{ISL-lab-3-lreg.svg}
PyPlot.close_figs() # hide


text:
include("objectives.jl")
include("constraints.jl")


text:
# wrappers of low-level functionality

function curandCreateGenerator(typ)
  handle_ref = Ref{curandGenerator_t}()
  res = @retry_reclaim err->isequal(err, CURAND_STATUS_ALLOCATION_FAILED) ||
                            isequal(err, CURAND_STATUS_INITIALIZATION_FAILED) begin
    unsafe_curandCreateGenerator(handle_ref, typ)
  end
  if res != CURAND_STATUS_SUCCESS
    throw_api_error(res)
  end
  handle_ref[]
end

function curandGetProperty(property::libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(curandGetProperty(CUDA.MAJOR_VERSION),
                          curandGetProperty(CUDA.MINOR_VERSION),
                          curandGetProperty(CUDA.PATCH_LEVEL))


text:
function initialize_system(	particle_type::String,
							R::Array{Float64, 2},
							Lx::Float64,
							Ly::Float64,
							Lz::Float64,
							phi_initial::Float64,
							position_constraint_axis::String,
							position_constraint_lower::Float64,
							position_constraint_upper::Float64,
							orientation_axis::Array{Float64, 1},
							orientation_constraint_axis::Array{Float64, 1},
							orientation_constraint_lower::Float64,
							orientation_constraint_upper::Float64)

	number_of_particles::Int64 = size(R, 1)

	# Determine initial system size.
	if phi_initial > 0.0 # If phi_initial specified, compute (Lx, Ly, Lz), assuming cubic domain.
		if particle_type == "sphere"
			Lx = ( sum(4.0 * pi / 3.0 * R[:, 1].^3) / phi_initial )^(1/3)
		elseif particle_type == "ellipse"
			warn("Initial volume fraction specified but not used for particle type 'ellipse'.")
		elseif particle_type == "ellipsoid"
			Lx = ( sum(4.0 * pi / 3.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / phi_initial )^(1/3)
		elseif particle_type == "cuboid"
			Lx = ( sum(8.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / phi_initial )^(1/3)
		end
		Ly = Lx
		Lz = Lx
	else # If (Lx, Ly, Lz) specified, compute phi_initial.
		if particle_type == "sphere"
			phi_initial = sum(4.0 * pi / 3.0 * R[:, 1].^3 ) / (Lx * Ly * Lz)
		elseif particle_type == "ellipse"
			warn("Not computing initial volume fraction for particle type 'ellipse'.")
		elseif particle_type == "ellipsoid"
			phi_initial = sum(4.0 * pi / 3.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / (Lx * Ly * Lz)
		elseif particle_type == "cuboid"
			phi_initial = sum(8.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / (Lx * Ly * Lz)
		end
	end

	# Generate random positions. If there are position constraints, satisfy these.
	X = Array{Float64}(number_of_particles)
	Y = Array{Float64}(number_of_particles)
	Z = Array{Float64}(number_of_particles)
	if position_constraint_axis == "x"
		X = position_constraint_lower * Lx + (position_constraint_upper - position_constraint_lower) * Lx * rand(number_of_particles)
	else
		X = Lx * rand(number_of_particles)
	end
	if position_constraint_axis == "y"
		Y = position_constraint_lower * Ly + (position_constraint_upper - position_constraint_lower) * Ly * rand(number_of_particles)
	else
		Y = Ly * rand(number_of_particles)
	end
	if position_constraint_axis == "z"
		Z = position_constraint_lower * Lz + (position_constraint_upper - position_constraint_lower) * Lz * rand(number_of_particles)
	else
		Z = Lz * rand(number_of_particles)
	end

	# Generate random orientations.
	Q0::Array{Float64, 1} = zeros(number_of_particles)
	Q1::Array{Float64, 1} = zeros(number_of_particles)
	Q2::Array{Float64, 1} = zeros(number_of_particles)
	Q3::Array{Float64, 1} = zeros(number_of_particles)
	q0::Float64 = 0.0
	q1::Float64 = 0.0
	q2::Float64 = 0.0
	q3::Float64 = 0.0
	if particle_type != "sphere"
		orientation_axis_rotated::Array{Float64, 1} = zeros(3)
		angle_to_orientation_constraint_axis::Float64 = 0.0
		is_orientation_ok::Bool = false

		for current_particle = 1:number_of_particles
			is_orientation_ok = false
			while !is_orientation_ok
				(q0, q1, q2, q3) = generate_random_unit_quaternion()
				(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(q0, q1, q2, q3)
				orientation_axis_rotated[1] = a11 * orientation_axis[1] + a12 * orientation_axis[2] + a13 * orientation_axis[3]
				orientation_axis_rotated[2] = a21 * orientation_axis[1] + a22 * orientation_axis[2] + a23 * orientation_axis[3]
				orientation_axis_rotated[3] = a31 * orientation_axis[1] + a32 * orientation_axis[2] + a33 * orientation_axis[3]
				angle_to_orientation_constraint_axis = acos(orientation_axis_rotated[1] * orientation_constraint_axis[1] + orientation_axis_rotated[2] * orientation_constraint_axis[2] + orientation_axis_rotated[3] * orientation_constraint_axis[3])
				if orientation_constraint_lower <= angle_to_orientation_constraint_axis <= orientation_constraint_upper
					Q0[current_particle] = q0
					Q1[current_particle] = q1
					Q2[current_particle] = q2
					Q3[current_particle] = q3
					is_orientation_ok = true
				end
			end
		end
	end

	# Characteristic/rotation matrix entries.
	A11::Array{Float64, 1} = zeros(number_of_particles)
	A12::Array{Float64, 1} = zeros(number_of_particles)
	A13::Array{Float64, 1} = zeros(number_of_particles)
	A21::Array{Float64, 1} = zeros(number_of_particles)
	A22::Array{Float64, 1} = zeros(number_of_particles)
	A23::Array{Float64, 1} = zeros(number_of_particles)
	A31::Array{Float64, 1} = zeros(number_of_particles)
	A32::Array{Float64, 1} = zeros(number_of_particles)
	A33::Array{Float64, 1} = zeros(number_of_particles)

	a11::Float64 = 0.0
	a12::Float64 = 0.0
	a13::Float64 = 0.0
	a21::Float64 = 0.0
	a22::Float64 = 0.0
	a23::Float64 = 0.0
	a31::Float64 = 0.0
	a32::Float64 = 0.0
	a33::Float64 = 0.0
	if particle_type == "ellipse"
		for current_particle = 1:number_of_particles
			(a11, a12, a13, a21, a22, a23, a31, a32, a33) = characteristic_matrix_ellipse(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle], R[current_particle, 1], R[current_particle, 2])
			A11[current_particle] = a11
			A12[current_particle] = a12
			A13[current_particle] = a13
			A21[current_particle] = a21
			A22[current_particle] = a22
			A23[current_particle] = a23
			A31[current_particle] = a31
			A32[current_particle] = a32
			A33[current_particle] = a33
		end
	elseif particle_type == "ellipsoid"
		for current_particle = 1:number_of_particles
			(a11, a12, a13, a21, a22, a23, a31, a32, a33) = characteristic_matrix_ellipsoid(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle], R[current_particle, 1], R[current_particle, 2], R[current_particle, 3])
			A11[current_particle] = a11
			A12[current_particle] = a12
			A13[current_particle] = a13
			A21[current_particle] = a21
			A22[current_particle] = a22
			A23[current_particle] = a23
			A31[current_particle] = a31
			A32[current_particle] = a32
			A33[current_particle] = a33
		end
	elseif particle_type == "cuboid"
		for current_particle = 1:number_of_particles
			(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle])
			A11[current_particle] = a11
			A12[current_particle] = a12
			A13[current_particle] = a13
			A21[current_particle] = a21
			A22[current_particle] = a22
			A23[current_particle] = a23
			A31[current_particle] = a31
			A32[current_particle] = a32
			A33[current_particle] = a33
		end
	end

	return (Lx,
			Ly,
			Lz,
			phi_initial,
			X,
			Y,
			Z,
			Q0,
			Q1,
			Q2,
			Q3,
			A11,
			A12,
			A13,
			A21,
			A22,
			A23,
			A31,
			A32,
			A33)
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Computations
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatFiles

dir = "X:/VIVES/1-Personal/Florian/git/OECD_ICIO/src/"
dir_raw = "C:/Users/u0148308/Dropbox/STORE/3. GVCs/tasks/task4_MRIO_to_WIOD/output/" # location of raw data

# other scripts
include(dir * "functions.jl") # Script with functions to import and transform raw data

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

N = 43 + 3 # 43 countries + 3 regions of BE
S = 55 # industries

# Notes:
#   - import data in long format and reshape => far quicker!

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

# Intermediate demand
df = DataFrame(StatFiles.load(dir_raw * "MRIO_WIOT_isic4_" * "IG" * ".dta"))
sort!(df)
df.ID = df.use_country .* "_" .* df.use_isic
df = unstack(df[:, Not([:use_country, :use_isic])], :ID, :value)
Z = Matrix(df[:, 3:end]) # NS×NS

# Notes
#   - missing values in Z => transform to zeros instead
any(ismissing.(Z)) # check for missing values
Z .= ifelse.(ismissing.(Z), 0.0, Z)
Z = Float64.(Z)

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

# Final demand
df = DataFrame(StatFiles.load(dir_raw * "MRIO_WIOT_isic4_" * "FD" * ".dta"))
sort!(df)
df.ID = df.use_country .* "_" .* df.use_isic
df = unstack(df[:, Not([:use_country, :use_isic])], :ID, :value)
Y = Matrix(df[:, 3:end])

any(ismissing.(Y)) # check for missing values

Y = [sum(Y[i, (j-1)*4+1:j*4]) for i in 1:N*S, j in 1:N] # NS×N, sum final demand components

# Notes:
#   - some entries of Y are negative due to inventories
#       + inventory adjustment as in Antras and Chor (2018)?
any(Y .< 0.0) # check for negative values

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

# Value added
df = DataFrame(StatFiles.load(dir_raw * "MRIO_WIOT_isic4_" * "VA" * ".dta"))
sort!(df)
df = unstack(df, :supply_isic, :value)
transform!(df, [:TAX,:VA] => ByRow((t, v) -> t + v) => :NET)
W = Float64.(df.VA) # NS×1

# Notes:
#   - check wether we should use VA or NET?

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

# Gross output
X_1 = W .+ [sum(Z[:, j]) for j in 1:N*S] # NS×1, gross output computed from intermediate demand imports and VA
X_2 = [sum(Z[i, :]) + sum(Y[i, :]) for i in 1:N*S] # NS×1, gross output computed from intermediate and final demand exports

# Notes:
#   - X_1 and X_2 are pretty different
X = copy(X_2)

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

# Gross domestic product
GDP = [sum(W[(i-1)*S+1:i*S]) for i in 1:N] # N×1, sum of value added per country

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

# Domar weights
DW = X ./ repeat(GDP, inner=S) # NS×1
DW_WORLD = X ./ sum(GDP) # NS×1

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

# Technical coefficients
A = Z ./ repeat(X', N*S) # NS×NS, input coefficients
any(isnan.(A))
A .= ifelse.(isnan.(A), 0.0, A) # in case X = 0 (otherwise we cannot compute Leontief inverse since NaN in A)

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

# Leontief inverse
L = inv(I - A) # NS×NS, Leontief inverse
any(isnan.(L))
any(isinf.(L))

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

# Total requirements matrix
L * Z # NS×NS

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

# Country-industry final demand multipliers
FD_multi = [sum(L[(i-1)*S+1:i*S, j]) for i in 1:N, j in 1:N*S] # N×NS
FD_multi = reshape(FD_multi', (S, N*N)) # S×NN

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


text:
using Test, JuliaExamplePkg

@test solve_problem() == [5;7;9]