text:
#!/usr/bin/env julia
# Script to check if registered packages have changed url's

using Requests, URIParser
fix = true && !("--dry-run" in ARGS)
tofix = []
asyncmap(readdir(Pkg.dir("METADATA"))) do pkg
    urlfile = Pkg.dir("METADATA", pkg, "url")
    if isfile(urlfile)
        url = readchomp(urlfile)
        req = get(replace(replace(url, "git://", "https://"),
            ".jl.git", ".jl"); allow_redirects = false)
        if !(endswith(url, "/$pkg.jl.git") #= || endswith(url, "/$pkg.jl") =# )
            warn("Unexpected end of url for $pkg: $url")
        end
        if statuscode(req) == 404
            info("$pkg Not Found")
            push!(tofix, (pkg, urlfile, url, req))
        elseif statuscode(req) == 301
            info("$pkg relocated from $url to $(req.headers["Location"])")
            push!(tofix, (pkg, urlfile, url, req))
        elseif statuscode(req) != 200
            warn("$pkg unexpected status: $req")
        end
    elseif isdir(Pkg.dir("METADATA", pkg)) && pkg != ".test" && pkg != ".git"
        println("No url file for $pkg")
    end
end
if fix
    info("Fixing redirects")
    for (pkg, urlfile, url, req) in tofix
        uri = URI(url)
        if uri.host != "github.com"
            warn("Skipping non-github package at $url")
            continue
        end
        status = statuscode(req)
        if endswith(uri.path, "/$pkg.jl.git")
            if status == 404
                newurl = "$(uri.scheme)://github.com/JuliaPackageMirrors/$pkg.jl.git"
            elseif status == 301
                reluri = URI(req.headers["Location"])
                if reluri.host != "github.com"
                    warn("Skipping non-github redirect at $reluri")
                    continue
                end
                if endswith(reluri.path, "/$pkg.jl")
                    newurl = "$(uri.scheme)://github.com$(reluri.path).git"
                else
                    warn("Not changing $url to $reluri because of unexpected end of new url")
                    continue
                end
            end
            newreq = get(replace(replace(newurl, "git://", "https://"),
                ".jl.git", ".jl"); allow_redirects = false)
            if statuscode(newreq) == 200
                info("Replacing $url with $newurl")
                open(urlfile, "w") do f
                    println(f, newurl)
                end
            else
                warn("Package from $url returned $status but not found at $newurl")
            end
        else
            warn("Not fixing $status from $url because of unexpected end of url")
        end
    end
end


text:
@testset "compute_int_I.jl" begin

end

ext: jl


text:
function moment0(spectrum::Vector{<:Real}, channel_width::Real)
    return sum(spectrum) * channel_width;
end

function moment0(cube::AbstractSpectralCube)
    m0 = mapslices(i->moment0(i, cube.z_resolution), cube.data; dims=[3])[:, :, 1]
    header = header2d_of(cube, cube.header["LINE"], "Moment0")
    return create_map(m0, header)
end


text:
# Functions to estimate the derivative(s) of the mean spectrum (probably could do better w/ spline or GP)
function calc_deriv_simple{T}(x::AbstractArray{T,1})
  @assert length(x)>=3
  dx = similar(x)
  dx[1] = x[2]-x[1]
  dx[end] = x[end]-x[end-1]
  for i in 2:(length(x)-1)
     dx[i] = (x[i+1]-x[i-1])/2
  end
  return dx
end

function calc_deriv2_simple{T}(x::AbstractArray{T,1})
  @assert length(x)>=3
  dx2 = similar(x)
  dx2[1] = x[3]-2*x[2]+x[1]
  dx2[end] = x[end]-2*x[end-1]+x[end-2]
  for i in 2:(length(x)-1)
     dx2[i] = (x[i+1]+x[i-1]-2*x[i])
  end
  return dx2
end

function calc_doppler_component_simple{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,1})
  @assert length(lambda) == length(flux)
  dlambdadpix = calc_deriv_simple(lambda);
  dfluxdpix = calc_deriv_simple(flux);
  doppler_basis = dfluxdpix.*(lambda./dlambdadpix)
end

function calc_doppler_component_simple{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,2})
  doppler_basis = calc_doppler_component_simple(lambda,vec(mean(flux,2)))
end

function calc_doppler_quadratic_term_simple{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,1})
  @assert length(lambda) == length(flux)
  dlambdadpix = calc_deriv_simple(lambda);
  d2fluxdpix2 = calc_deriv2_simple(flux);
  doppler_quad_term = 0.5*d2fluxdpix2.*(lambda.^2./dlambdadpix.^2)
end

function calc_doppler_quadratic_term_simple{T}(lambda::AbstractArray{T,1}, flux::AbstractArray{T,2})
  doppler_quad_term = calc_doppler_quadratic_term_simple(lambda,vec(mean(flux,2)))
end

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
const AbstractArray4{T} = AbstractArray{T,4}

δkron(T::Type{<:Number}, i::Integer, j::Integer) = i == j ? one(T) : zero(T)

struct Id2{dim,T<:Number} <: AbstractMatrix{T} end
@pure Base.size(::Id2{dim}) where {dim} = (dim, dim)
Base.getindex(::Id2{dim,T}, i::Integer, j::Integer) where {dim,T} = δkron(T, i, j)
function Base.replace_in_print_matrix(::Id2, i::Integer, j::Integer, s::AbstractString)
    i == j ? s : Base.replace_with_centered_mark(s)
end

isidentity(a::AbstractMatrix{T}) where {T} = a ≈ Id2{size(a, 1),T}()
isidentity(a::AbstractMatrix{Sym}) = a == Id2{size(a, 1),Sym}()

isdiagonal(a::AbstractMatrix{T}) where {T} = norm(a - Diagonal(a)) <= eps(T)
isdiagonal(a::AbstractMatrix{Sym}) = isdiag(a)

simplifyif(x) = x
simplifyif(x::Sym) = simplify(x)
simplifyif(m::Matrix{Sym}) = simplify.(m)
simplifyif(m::Symmetric{Sym}) = Symmetric(simplify.(m))

@inline LinearAlgebra.issymmetric(t::Tensor{2, 2, T}) where {T <: Real} = @inbounds t[1,2] ≈ t[2,1]

@inline function LinearAlgebra.issymmetric(t::Tensor{2, 3, T}) where {T <: Real}
    return @inbounds t[1,2] ≈ t[2,1] && t[1,3] ≈ t[3,1] && t[2,3] ≈ t[3,2]
end

function Tensors.isminorsymmetric(t::Tensor{4, dim, T}) where {dim, T <: Real}
    @inbounds for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
        if !(t[i,j,k,l] ≈ t[j,i,k,l]) || !(t[i,j,k,l] ≈ t[i,j,l,k])
            return false
        end
    end
    return true
end

function Tensors.ismajorsymmetric(t::FourthOrderTensor{dim, T}) where {dim, T <: Real}
    @inbounds for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
        if !(t[i,j,k,l] ≈ t[k,l,i,j])
            return false
        end
    end
    return true
end

@inline function Tensors.majortranspose(S::SymmetricTensor{4, dim}) where {dim}
    SymmetricTensor{4, dim}(@inline function(i, j, k, l) @inbounds S[k,l,i,j]; end)
end

function Tensors.otimes(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    ec1 = ntuple(i -> i, order1)
    ec2 = ntuple(i -> order1 + i, order2)
    ec3 = ntuple(i -> i, order1 + order2)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

function contract(t::AbstractArray{T,order}, i::Integer, j::Integer) where {T,order}
    m = min(i, j)
    M = max(i, j)
    ec1 = ntuple(k -> k == j ? i : k, order)
    ec2 = (Tuple(1:m-1)..., Tuple(m+1:M-1)..., Tuple(M+1:order)...)
    return einsum(EinCode((ec1,), ec2), (AbstractArray{T}(t),))
end

contract(t::AbstractArray{T,2}, ::Integer, ::Integer) where {T} = tr(t)

function Tensors.dcontract(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    newc = order1 + order2
    ec1 = (ntuple(i -> i, order1 - 2)..., newc, newc + 1)
    ec2 = (newc, newc + 1, ntuple(i -> order1 - 2 + i, order2 - 2)...)
    ec3 = ntuple(i -> i, order1 + order2 - 4)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

Tensors.dcontract(t1::AbstractArray{T1,2}, t2::AbstractArray{T2,2}) where {T1,T2} =
    dot(AbstractArray{T1}(t1), AbstractArray{T2}(t2))

function Tensors.dotdot(
    v1::AbstractArray{T1,order1},
    S::AbstractArray{TS,orderS},
    v2::AbstractArray{T2,order2},
) where {T1,TS,T2,order1,orderS,order2}
    newc = order1 + orderS
    ec1 = (ntuple(i -> i, order1 - 1)..., newc)
    ecS = (newc, ntuple(i -> order1 - 1 + i, orderS - 1)...)
    ec3 = ntuple(i -> i, order1 + orderS - 2)
    v1S = einsum(EinCode((ec1, ecS), ec3), (AbstractArray{T1}(v1), AbstractArray{TS}(S)))
    newc += order2
    ecv1S = (ntuple(i -> i, order1 + orderS - 3)..., newc)
    ec2 = (newc, ntuple(i -> order1 + orderS - 3 + i, order2 - 1)...)
    ec3 = ntuple(i -> i, newc - 4)
    return einsum(EinCode((ecv1S, ec2), ec3), (v1S, AbstractArray{T2}(v2)))
end

function qcontract(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    newc = order1 + order2
    ec1 = (ntuple(i -> i, order1 - 4)..., newc, newc + 1, newc + 2, newc + 3)
    ec2 = (newc, newc + 1, newc + 2, newc + 3, ntuple(i -> order1 - 4 + i, order2 - 4)...)
    ec3 = ntuple(i -> i, order1 + order2 - 8)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

qcontract(t1::AbstractArray{T1,4}, t2::AbstractArray{T2,4}) where {T1,T2} =
    dot(AbstractArray{T1}(t1), AbstractArray{T2}(t2))

function Tensors.otimesu(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    ec1 = (ntuple(i -> i, order1 - 1)..., order1 + 1)
    ec2 = (order1, ntuple(i -> order1 + 1 + i, order2 - 1)...)
    ec3 = ntuple(i -> i, order1 + order2)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

function Tensors.otimesl(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    ec1 = (ntuple(i -> i, order1 - 1)..., order1 + 2)
    ec2 = (order1, order1 + 1, ntuple(i -> order1 + 2 + i, order2 - 2)...)
    ec3 = ntuple(i -> i, order1 + order2)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

otimesul(t1::AbstractArray{T1}, t2::AbstractArray{T2}) where {T1,T2} =
    (otimesu(t1, t2) + otimesl(t1, t2)) / promote_type(T1, T2)(2)

function sotimes(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    ec1 = ntuple(i -> i, order1)
    ec2 = ntuple(i -> order1 + i, order2)
    ec3 = ntuple(i -> i, order1 + order2)
    t3 = einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
    ec1 = (ntuple(i -> i, order1 - 1)..., order1 + 1)
    ec2 = (order1, ntuple(i -> order1 + 1 + i, order2 - 1)...)
    ec3 = ntuple(i -> i, order1 + order2)
    t4 = einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
    return (t3 + t4) / promote_type(T1, T2)(2)
end

@inline function sotimes(S1::Vec{dim}, S2::Vec{dim}) where {dim}
    return SymmetricTensor{2,dim}(@inline function (i, j)
        @inbounds (S1[i] * S2[j] + S1[j] * S2[i]) / 2
    end)
end

@inline function sotimes(S1::SecondOrderTensor{dim}, S2::SecondOrderTensor{dim}) where {dim}
    TensorType = Tensors.getreturntype(
        otimes,
        Tensors.get_base(typeof(S1)),
        Tensors.get_base(typeof(S2)),
    )
    TensorType(@inline function (i, j, k, l)
        @inbounds (S1[i, j] * S2[k, l] + S1[i, k] * S2[j, l]) / 2
    end)
end

Tensors.otimes(α::Number, t::AbstractArray) = α * t
Tensors.otimes(t::AbstractArray, α::Number) = α * t

sotimes(α::Number, t::AbstractArray) = α * t
sotimes(t::AbstractArray, α::Number) = α * t

const ⊙ = qcontract
const ⊠ = otimesu
const ⊠ˢ = otimesul
const ⊗ˢ = sotimes
const sboxtimes = otimesul

export isidentity, contract, qcontract, otimesu, otimesul, sboxtimes, sotimes, ⊙, ⊠, ⊠ˢ, ⊗ˢ
export ⋅, ⊡, ⊗


text:
@inline squared(x::Real) = x*x
@inline squared(x::Vector) = x.*x
@inline maxzero(x::Real) = max(x,zero(x))


text:
"""
    read_surfaces(filename)

Read surface triangulation from file.
"""
function read_surfaces(filename)

    println("Loading surface triangulation from " * filename)

    # Get list of nodes from .node file
    points = open(filename * ".node", "r") do io
        readline(io)
        readline(io)
        npoint = parse(Int, split(readline(io))[1])
        readline(io)
        points = zeros(3, npoint)
        for ipoint = 1:npoint
            tmp = parse.(Float64, split(readline(io)))
            points[:, ipoint] = tmp[2:4]
        end
        points
    end

    # Read .poly file
    facets, facetmarkers, regions = open(filename * ".poly", "r") do io

        # Read list of holes (refer to separate file)
        readline(io)
        readline(io)
        readline(io)

        # Read list of facets
        readline(io)
        readline(io)
        nfacet = parse(Int, split(readline(io))[1])
        readline(io)
        facets = zeros(Int, 3, nfacet)
        facetmarkers = zeros(Int, nfacet)
        for ifacet = 1:nfacet
            tmp = parse.(Int, split(readline(io)))
            facetmarkers[ifacet] = tmp[3]
            tmp = parse.(Int, split(readline(io)))
            facets[:, ifacet] = tmp[2:4]
        end

        # Read list of holes (empty)
        readline(io)
        readline(io)

        # Read list of interior points with their corresponding compartment
        readline(io)
        nregion = parse(Int, split(readline(io))[1])
        regions = zeros(3, nregion)
        for iregion = 1:nregion
            tmp = parse.(Float64, split(readline(io)))
            regions[:, iregion] = tmp[2:4]
        end

        facets, facetmarkers, regions
    end

    (; points, facets, facetmarkers, regions)
end


text:
α = 0.025

@testset "multiprocessing tests" begin
  addprocs(2)
  eval(Expr(:toplevel, :(@everywhere using DatagenCopulaBased)))
  eval(Expr(:toplevel, :(@everywhere import DatagenCopulaBased: nestedfrankgen)))
  eval(Expr(:toplevel, :(@everywhere using Random)))
  eval(Expr(:toplevel, :(@everywhere Random.seed!(44))))
  @test nestedfrankgen(5., 3., [1, 1, 2]) == [9, 54, 63]
  x = nestedarchcopulagen(250000, [3, 2],  [8., 10.], 2., "frank", 2)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  @test corkendall(x)[1:4,1] ≈ [1.0, 0.60262, 0.60262, 0.2139] atol=1.0e-2
  @test corkendall(x)[3:5,4] ≈ [0.2139, 1.0, 0.6658] atol=1.0e-2
  @test corkendall(x)[6:7,6] ≈ [1.0, 0.2139] atol=1.0e-2
  @test tail(x[:,4], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,4], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "r", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
end


text:
@testset "FiniteDiff comparison" begin
    Random.seed!(100)
    x0 = srand(3)
    v0 = srand(3)
    q0 = rand(QuatRotation).q
    ϕ0 = srand(3)
    timestep0 = 0.01
    x1 = Dojo.next_position(x0, v0, timestep0)
    q1 = Dojo.next_orientation(q0, ϕ0, timestep0)

    # ∇v
    FD∂i∂v = FiniteDiff.finite_difference_jacobian(v0 -> Dojo.next_position(x0, v0, timestep0), v0)
    @test norm(FD∂i∂v - Dojo.linear_integrator_jacobian_velocity(timestep0), Inf) < 1.0e-6
    # ∇ϕ
    FD∂i∂ϕ = FiniteDiff.finite_difference_jacobian(ϕ0 -> Dojo.vector(Dojo.next_orientation(q0, ϕ0, timestep0)), ϕ0)
    @test norm(FD∂i∂ϕ - Dojo.rotational_integrator_jacobian_velocity(q0, ϕ0, timestep0), Inf) < 1.0e-6

    # ∇x
    FD∂i∂x = FiniteDiff.finite_difference_jacobian(x0 -> Dojo.next_position(x0, v0, timestep0), x0)
    @test norm(FD∂i∂x - Dojo.linear_integrator_jacobian_position(), Inf) < 1.0e-6
    # ∇q
    FD∂i∂q = FiniteDiff.finite_difference_jacobian(q0 ->
        Dojo.vector(Dojo.next_orientation(Quaternion(q0...), ϕ0, timestep0)),
        Dojo.vector(q0))
    @test norm(FD∂i∂q - Dojo.rotational_integrator_jacobian_orientation(q0, ϕ0, timestep0, attjac=false), Inf) < 1.0e-6

    # ∇vϕ
    vel0 = [v0; ϕ0]
    FDvel = FiniteDiff.finite_difference_jacobian(vel0 ->
        [Dojo.next_position(x0, SVector{3}(vel0[1:3]), timestep0);
        Dojo.vector(Dojo.next_orientation(q0, SVector{3}(vel0[4:6]), timestep0))], vel0)
    @test norm(FDvel - Dojo.integrator_jacobian_velocity(q0, ϕ0, timestep0), Inf) < 1.0e-6

    # ∇xq
    con0 = [x0; Dojo.vector(q0)]
    FDcon = FiniteDiff.finite_difference_jacobian(con0 ->
        [Dojo.next_position(SVector{3}(con0[1:3]), v0, timestep0);
        Dojo.vector(Dojo.next_orientation(Quaternion(con0[4:7]...), ϕ0, timestep0))], con0)
    @test norm(FDcon - Dojo.integrator_jacobian_configuration(q0, ϕ0, timestep0, attjac=false), Inf) < 1e-8
    @test norm(FDcon * cat(I(3), Dojo.LVᵀmat(q0), dims=(1,2)) - Dojo.integrator_jacobian_configuration(q0, ϕ0, timestep0, attjac=true), Inf) < 1e-8
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

Base.prod(x::NDArray; dims = :) = _prod(x, dims)
@_remap _prod(x::NDArray, ::Colon) prod(x)
@_remap _prod(x::NDArray, dims)    prod(x; axis = 0 .- dims, keepdims = true)

Base.maximum(x::NDArray; dims = :) = _nd_maximum(x, dims)
@_remap _nd_maximum(x::NDArray, ::Colon) max(x)
@_remap _nd_maximum(x::NDArray, dims)    max(x; axis = 0 .- dims, keepdims = true)

Base.minimum(x::NDArray; dims = :) = _nd_minimum(x, dims)
@_remap _nd_minimum(x::NDArray, ::Colon) min(x)
@_remap _nd_minimum(x::NDArray, dims)    min(x; axis = 0 .- dims, keepdims = true)

###############################################################################
# min/max
###############################################################################

import Base: min, max

broadcasted(::typeof(max), x::NDArray{T}, y::NDArray{T}) where {T} =
  _broadcast_maximum(x, y)

broadcasted(::typeof(min), x::NDArray{T}, y::NDArray{T}) where {T} =
  _broadcast_minimum(x, y)

###############################################################################
# argmin/argmax
###############################################################################

# TODO: support CartesianIndex ?
"""
    argmax(x::NDArray; dims) -> indices

Note that `NaN` is treated as greater than all other values in `argmax`.

## Examples

```julia-repl
julia> x = NDArray([0. 1 2; 3 4 5])
2×3 NDArray{Float64,2} @ CPU0:
 0.0  1.0  2.0
 3.0  4.0  5.0

julia> argmax(x, dims = 1)
1×3 NDArray{Float64,2} @ CPU0:
 2.0  2.0  2.0

julia> argmax(x, dims = 2)
2×1 NDArray{Float64,2} @ CPU0:
 3.0
 3.0
```

See also [`argmin`](@ref mx.argmin).
"""
Base.argmax(x::NDArray; dims = :) = _argmax(x, dims) .+ 1
@_remap _argmax(x::NDArray, ::Colon) argmax(x)
@_remap _argmax(x::NDArray, dims)    argmax(x; axis = 0 .- dims, keepdims = true)

"""
    argmin(x::NDArray; dims) -> indices

Note that `NaN` is treated as less than all other values in `argmin`.

## Examples

```julia-repl
julia> x = NDArray([0. 1 2; 3 4 5])
2×3 NDArray{Float64,2} @ CPU0:
 0.0  1.0  2.0
 3.0  4.0  5.0

julia> argmax(x, dims = 1)
1×3 NDArray{Float64,2} @ CPU0:
 2.0  2.0  2.0

julia> argmax(x, dims = 2)
2×1 NDArray{Float64,2} @ CPU0:
 3.0
 3.0
```

See also [`argmax`](@ref mx.argmax).
"""
Base.argmin(x::NDArray; dims = :) = _argmin(x, dims) .+ 1
@_remap _argmin(x::NDArray, ::Colon) argmin(x)
@_remap _argmin(x::NDArray, dims)    argmin(x; axis = 0 .- dims, keepdims = true)

################################################################################
# remapping to solving type unstablility
################################################################################

@_remap _broadcast_maximum(x::NDArray, y::NDArray)  broadcast_maximum(x, y)
@_remap _broadcast_maximum!(x::NDArray, y::NDArray) broadcast_maximum(x, y)

@_remap _broadcast_minimum(x::NDArray, y::NDArray)  broadcast_minimum(x, y)
@_remap _broadcast_minimum!(x::NDArray, y::NDArray) broadcast_minimum(x, y)


text:
using Plots, Plots.PlotMeasures
import Fontconfig, Cairo
# import Reel, Compose
using LightGraphs, GraphPlot
using Random

"""
    plotscenarios(scen::Array{Float64,2}; savepath::String ="",
        save::Bool = false, fileformat::String = "png", kwargs...)

Plots a set of scenarios.

Parameters:
-   scen        : Scenarios matrix (Stages x Scenarious).
-   quants      : Lower Quantiles to plot (their simmetric will be taken into account)
-   save        : Bool to indicate if figure is to be saved.
-   savepath    : Path save figure.
-   fileformat  : Figure file format.
-   kwargs      : Aditional keyword arguments for plot function.
"""
function plotscenarios(scen::Array{Float64,2}; quants=collect(0.05:0.1:0.25),
        savepath::String ="",
        save::Bool = false, fileformat::String = "png", kwargs...)

    med_scen = Statistics.median(scen;dims=2)
    nstag,nscen = size(scen)

    # plot if all scenarious are equal

    if prod(isapprox.(med_scen-quantile_scen(scen,0.0), 0.0, atol=1e-1))
        p1 = plot(round.(med_scen,digits =5),
            color = "Steel Blue",
            xticks = (collect(1:Int(floor(nstag/4)):nstag), [string(i) for  i in collect(1:Int(floor(nstag/4)):nstag)]),
            label = "Max and Min";
            kwargs...)

        if save
            savefig(p1, savepath*"$fileformat")
            return nothing
        end
        return p1
    end

    # plot
    p1 = plot(med_scen, ribbon=(med_scen-quantile_scen(scen,0.0),quantile_scen(scen,1.0)-med_scen),
                color = "gray",
                xticks = (collect(1:Int(floor(nstag/4)):nstag), [string(i) for  i in collect(1:Int(floor(nstag/4)):nstag)]),
                label = "Median";
                kwargs...)
    for q in quants
        plot!(p1, med_scen, ribbon=(med_scen-quantile_scen(scen,q),quantile_scen(scen,1-q)-med_scen),
                            color = "gray", label = "")
    end
    plot!(p1, maximum(scen; dims=2), label = "Max and Min", color = "Steel Blue")
    plot!(p1, minimum(scen; dims=2), label = "", color = "Steel Blue")
    plot!(p1, mean(scen; dims=2), label = "Mean", color = "Black")
    if save
        savefig(p1, savepath*"$fileformat")
        return nothing
    end

    return p1
end

"""
    HydroPowerModels.plotresults(results::Dict;nc::Int = 3)

Common Plots.

Parameters:
-   results        : Simulations output.
-   nc             : Number of figures per line.

"""
function plotresults(results::Dict;nc::Int = 3)
    plt_total = Array{Plots.Plot}(undef,10000)
    nplots = 0
    nsim = length(results[:simulations])
    nstages = length(results[:simulations][1])

    # Thermal Generation
    ngen = length(results[:data][1]["powersystem"]["gen"])
    idxhyd = idx_hydro(results[:data][1])
    idxgen = setdiff(collect(1:ngen),idxhyd)
    baseMVA =  [results[:simulations][i][j][:powersystem]["solution"]["baseMVA"] for i=1:nsim, j=1:nstages]'
    scen_gen = [[results[:simulations][i][j][:powersystem]["solution"]["gen"]["$gen"]["pg"] for i=1:nsim, j=1:nstages]'.*baseMVA for gen =1:ngen]

    plt =   [plotscenarios(scen_gen[gen], title  = "Thermal Generation $gen",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for gen in idxgen
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Thermal Reactive Generation
    if results[:params]["model_constructor_grid_forward"] == PowerModels.ACPPowerModel
        scen_qgen = [[results[:simulations][i][j][:powersystem]["solution"]["gen"]["$gen"]["qg"] for i=1:nsim, j=1:results[:params]["stages"]]'.*baseMVA for gen =1:ngen]
        plt =   [plotscenarios(scen_qgen[gen], title  = "Thermal Reactive Generation $gen",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for gen in idxgen
        ]
        plt_total[nplots+1:nplots+length(plt)] = plt
        nplots +=length(plt)
    end

    # circuit MVA
    baseMVA = results[:data][1]["powersystem"]["baseMVA"]

    # Branch flow

    nbrc = length(results[:data][1]["powersystem"]["branch"])
    idxbrc = collect(1:nbrc)
    scen_branch = [[results[:simulations][i][j][:powersystem]["solution"]["branch"]["$brc"]["pf"] for i=1:nsim, j=1:nstages]'.*baseMVA for brc =1:nbrc]

    plt =   [plotscenarios(scen_branch[brc], title  = "Branch Flow $brc",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for brc in idxbrc
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Branch Reactive flow

    if results[:params]["model_constructor_grid_forward"] == PowerModels.ACPPowerModel
        scen_branch_qf = [[results[:simulations][i][j][:powersystem]["solution"]["branch"]["$brc"]["qf"] for i=1:nsim, j=1:results[:params]["stages"]]'.*baseMVA for brc =1:nbrc]
        plt =   [plotscenarios(scen_branch_qf[brc], title  = "Branch Reactive Flow $brc",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for brc in idxbrc
        ]
        plt_total[nplots+1:nplots+length(plt)] = plt
        nplots +=length(plt)
    end

    # Voltage angle
    try
        nbus = length(results[:data][1]["powersystem"]["bus"])
        idxbus = collect(1:nbus)
        scen_va = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["va"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])

        plt =   [plotscenarios(scen_va[bus], title  = "Voltage angle $bus",
                    ylabel = "Radians",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for bus in idxbus
        ]
        plt_total[nplots+1:nplots+length(plt)] = plt
        nplots +=length(plt)
    catch
    end

    # Nodal price
    try
    nbus = length(results[:data][1]["powersystem"]["bus"])
    idxbus = collect(1:nbus)
    scen_pld = convert(Array{Array{Float64,2},1},[[-results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["lam_kcl_r"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])/baseMVA

    plt =   [plotscenarios(scen_pld[bus], title  = "Nodal price bus $bus",
                ylabel = "\$/MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for bus in idxbus
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)
    catch
    end
    # Deficit
    try
    nbus = length(results[:data][1]["powersystem"]["bus"])
    idxbus = collect(1:nbus)
    scen_def = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["deficit"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])

    plt =   [plotscenarios(scen_def[bus].*baseMVA, title  = "Deficit bus $bus",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
            for bus in idxbus
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)
    catch
    end

    # Hydro Generation
    nHyd = results[:data][1]["hydro"]["nHyd"]

    plt =   [   plotscenarios(scen_gen[gen], title  = "Hydro Generation $gen",
                    ylabel = "MW",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for gen in idxhyd
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Reservoir Outflow
    scen_turn = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:outflow][res] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])
    plt =   [   plotscenarios(scen_turn[res], title  = "Reservoir Outflow $res",
                    ylabel = "m³/s",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for res = 1:nHyd
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Hydro Spill
    scen_spill = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:spill][res] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    plt =   [   plotscenarios(scen_spill[res], title  = "Hydro Spill $res",
                    ylabel = "Hm³",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for res = 1:nHyd
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Reservoir Volume
    scen_voume = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:reservoir][res].out for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    plt =   [   plotscenarios(scen_voume[res], title  = "Volume Reservoir $res",
                    ylabel = "Hm³",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for res = 1:nHyd
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)

    # Inflows
    scen_inflows = convert(Array{Array{Float64,2},1},[[results[:data][1]["hydro"]["Hydrogenerators"][res]["inflow"][cidx(j,results[:data][1]["hydro"]["size_inflow"][1]),results[:simulations][i][j][:noise_term]] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    plt =   [   plotscenarios(scen_inflows[res], title  = "Inflows Reservoir $res",
                    ylabel = "m³/s",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
                for res = 1:nHyd
    ]
    plt_total[nplots+1:nplots+length(plt)] = plt
    nplots +=length(plt)
    if mod(nplots,nc) > 0 && floor(Int,nplots/nc) > 0
        l = @layout [ Plots.grid(floor(Int,nplots/nc),nc);  Plots.grid(1,mod(nplots,nc))]
        nlines = floor(Int,nplots/nc)+1
        l.heights = Plots.grid(2,1,heights=[floor(Int,nplots/nc)/nlines;1/nlines]).heights
    elseif floor(Int,nplots/nc) > 0
        l = @layout Plots.grid(floor(Int,nplots/nc),nc)
        nlines = floor(Int,nplots/nc)
        l.heights = Plots.grid(nlines,1,heights=[1/nlines for n = 1:nlines]).heights
    else
        l = @layout Plots.grid(1,mod(nplots,nc))
        nlines = 1
        l.heights = Plots.grid(1,1,heights=[1]).heights
    end

    return plot(plt_total[1:nplots]...,layout=l,size = (nc*400,400*ceil(Int,nplots/nc)),legend=false)
end

"""Common descriptive statistics"""
function descriptivestatistics_results(results::Dict;nitem::Int = 3,quants::Array{Float64}=[0.25;0.5;0.75])

    dcp_stats = Dict()

    nsim = length(results[:simulations])
    nstages = length(results[:simulations][1][1][:powersystem]["solution"])

    # Thermal Generation first nitem gen

    dcp_stats["pg"] = Dict()

    ngen = length(results[:data][1]["powersystem"]["gen"])
    idxhyd = idx_hydro(results[:data][1])
    idxgen = setdiff(collect(1:min(ngen,nitem)),idxhyd)
    baseMVA =  [results[:simulations][i][j][:powersystem]["solution"]["baseMVA"] for i=1:nsim, j=1:nstages]'
    scen_gen = [[results[:simulations][i][j][:powersystem]["solution"]["gen"]["$gen"]["pg"] for i=1:nsim, j=1:nstages]'.*baseMVA for gen =1:ngen]

    for i = 1:size(idxgen,1)
        gen = idxgen[i]
        dcp_stats["pg"]["$i"] = quantile_scen(scen_gen[gen], quants, output_dict=true)
    end

    # Branch flow first nitem brc

    dcp_stats["pf"] = Dict()

    nbrc = length(results[:data][1]["powersystem"]["branch"])
    idxbrc = collect(1:min(nbrc,nitem))
    scen_branch = [[results[:simulations][i][j][:powersystem]["solution"]["branch"]["$brc"]["pf"] for i=1:nsim, j=1:nstages]'.*baseMVA for brc =1:nbrc]

    for i = 1:size(idxbrc,1)
        brc = idxbrc[i]
        dcp_stats["pf"]["$i"] = quantile_scen(scen_branch[brc], quants, output_dict=true)
    end

    # Voltage angle first nitem bus

    if results[:params]["model_constructor_grid"] != PowerModels.AbstractPowerModel{PowerModels.SOCWRForm}

        dcp_stats["va"] = Dict()

        nbus = length(results[:data][1]["powersystem"]["bus"])
        idxbus = collect(1:min(nbus,nitem))
        scen_va = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["va"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])

        for i = 1:size(idxbus,1)
            bus = idxbus[i]
            dcp_stats["va"]["$i"] = quantile_scen(scen_va[bus], quants, output_dict=true)
        end

    end
    # Hydro Generation first nitem Hydro

    dcp_stats["hpg"] = Dict()

    for i = 1:size(idxhyd[1:min(results[:data][1]["hydro"]["nHyd"],nitem)],1)
        gen = idxhyd[1:min(results[:data][1]["hydro"]["nHyd"],nitem)][i]
        dcp_stats["hpg"]["$i"] = quantile_scen(scen_gen[gen], quants, output_dict=true)
    end

    # Reservoir Volume first nitem Hydro

    dcp_stats["volume"] = Dict()

    scen_voume = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:reservoir][res].out for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    for res = 1:min(results[:data][1]["hydro"]["nHyd"],nitem)
        dcp_stats["volume"]["$res"] = quantile_scen(scen_voume[res], quants, output_dict=true)
    end

    return dcp_stats
end

"""
    HydroPowerModels.plot_grid(data::Dict;path=nothing,size_fig = [15cm, 15cm],node_label=false,nodelabeldist=4.5)

Plot Grid installed Power.

Paremeters:
-   data         : HydroPowerModel single stage data.
-   path         : Path to save grid plot.
-   size_fig     : Size figure.
-   node_label   : Plot nodel label on grid.
-   nodelabeldist: Nodel label distance from node.
"""
function plot_grid(data::Dict;path=nothing,size_fig = [15cm, 15cm],node_label=false,nodelabeldist=4.5)

    gatherusefulinfo!(data)

    nbus = length(data["powersystem"]["bus"])

    nNodes = nbus

    g = Graph(nbus)

    for brc in values(data["powersystem"]["branch"])
        add_edge!(g, brc["f_bus"], brc["t_bus"])
    end

    # nodes gen
    thermal_nodes = fill(0.0,nbus)
    hydro_nodes = fill(0.0,nbus)

    idxhyd = idx_hydro(data)

    for i in 1:length(data["powersystem"]["gen"])
        bus_i = data["powersystem"]["gen"]["$i"]["gen_bus"]
        if i in idxhyd
            hydro_nodes[bus_i] += data["powersystem"]["gen"]["$i"]["pmax"]*data["powersystem"]["gen"]["$i"]["mbase"]
        else
            thermal_nodes[bus_i] += data["powersystem"]["gen"]["$i"]["pmax"]*data["powersystem"]["gen"]["$i"]["mbase"]
        end
    end

    # nodes loads
    load_nodes = fill(0.0,nbus)

    for i in 1:length(data["powersystem"]["load"])
        bus_i = data["powersystem"]["load"]["$i"]["load_bus"]
        load_nodes[bus_i] += data["powersystem"]["load"]["$i"]["pd"]*data["powersystem"]["baseMVA"]
    end

    # number of nodes
    num_nodes = nbus+ sum(load_nodes .> 0)+sum(hydro_nodes .> 0)+sum(thermal_nodes .> 0)

    # node size
    nodesize = fill(0.0,num_nodes)

    # nodes membership (3: Hydro, 2: Thermal)
    membership = fill(1,num_nodes)

    # node label
    if node_label
        nodelabel = [1:nbus;fill("",sum(hydro_nodes .> 0));fill("",sum(thermal_nodes .> 0));fill("",sum(load_nodes .> 0))]
    else
        nodelabel = nothing
    end
    # create nodes
    for bus_i in 1:length(data["powersystem"]["bus"])
        if hydro_nodes[bus_i] > 0
            add_vertex!(g)
            nNodes +=1
            add_edge!(g, nNodes, bus_i)
            membership[nNodes] = 3
            nodesize[nNodes] = log(hydro_nodes[bus_i])
        end

        if thermal_nodes[bus_i] > 0
            add_vertex!(g)
            nNodes +=1
            add_edge!(g, nNodes, bus_i)
            membership[nNodes] = 2
            nodesize[nNodes] = log(thermal_nodes[bus_i])
        end

        if load_nodes[bus_i] > 0
            add_vertex!(g)
            nNodes +=1
            add_edge!(g, nNodes, bus_i)
            membership[nNodes] = 4
            nodesize[nNodes] = log(load_nodes[bus_i])
        end

    end

    for n in 1:num_nodes
        nodesize[n] = nodesize[n] == 0 ? unique(sort(nodesize))[2] : nodesize[n]
    end

    nodecolor = ["black", "red", "blue", "orange"]

    # membership color
    nodefillc = nodecolor[Int64.(membership)]

    if path != nothing
        draw(PDF(path, size_fig...), gplot(g, nodefillc=nodefillc, nodesize=nodesize, nodelabel=nodelabel,nodelabeldist=nodelabeldist))
    else
        gplot(g, nodefillc=nodefillc, nodesize=nodesize, nodelabel=nodelabel,nodelabeldist=nodelabeldist)
    end
end

""" Plot Hydro Grid installed volume"""
function plot_hydro_grid(data::Dict;path=nothing,size_fig = [12cm, 12cm],node_label=false,nodelabeldist=8.5)

    gatherusefulinfo!(data)

    nHyd = data["hydro"]["nHyd"]

    nNodes = nHyd

    g = DiGraph(nNodes)

    # nodes bus
    node2bus = []

    # hydro_size
    hydro_size = fill(0.0,nHyd)

    for i=1:nHyd
        hydro = data["hydro"]["Hydrogenerators"][i]
        hydro_size[i] = hydro["max_volume"]

        for hyd in hydro["downstream_turn"]
            j = findall(x->x["index"]==hyd,data["hydro"]["Hydrogenerators"])
            if !isempty(j)
                add_edge!(g, i, j[1])
            end
        end
        for hyd in hydro["downstream_spill"]
            j = findall(x->x["index"]==hyd,data["hydro"]["Hydrogenerators"])
            if !isempty(j)
                add_edge!(g, i, j[1])
            end
        end

        if hydro["index_grid"] != nothing
            bus_i = data["powersystem"]["gen"]["$(hydro["i_grid"])"]["gen_bus"]
            node = findall(x->x==bus_i,node2bus)
            if isempty(node)
                append!(node2bus,bus_i)
                add_vertex!(g)
                nNodes +=1
                node = [size(node2bus,1)]
            end
            add_edge!(g, i, Int64(nHyd+node[1]))
        end
    end

    # node size
    nodesize = fill(maximum(log.(hydro_size.*1000,1.002))*0.5,nNodes)

    nodesize[1:nHyd] = log.(hydro_size.*1000,1.002)
    nodecolor = ["black", "blue"]

    # membership color
    membership = fill(1,nNodes)
    membership[1:nHyd] .= 2
    nodefillc = nodecolor[membership]

    # node label
    if node_label
        nodelabel = [fill("",nHyd);node2bus]
    else
        nodelabel = nothing
    end

    if path != nothing
        draw(PDF(path, size_fig...), gplot(g, nodefillc=nodefillc, nodesize=nodesize, nodelabel=nodelabel,nodelabeldist=nodelabeldist, arrowlengthfrac=0.005))
    else
        gplot(g, nodefillc=nodefillc, nodesize=nodesize, nodelabel=nodelabel,nodelabeldist=nodelabeldist, arrowlengthfrac=0.005)
    end
end

# """ Plot Grid Stage Dispatched Power"""
# function plot_grid_dispatched_stage(results::Dict,t::Int;quant::Float64=0.5,size_fig = [15cm, 15cm],node_label=false,nodelabeldist=4.5)

#     data = results[:data][1]

#     nbus = length(data["powersystem"]["bus"])

#     nNodes = nbus

#     g = Graph(nbus)

#     for brc in values(data["powersystem"]["branch"])
#         add_edge!(g, brc["f_bus"], brc["t_bus"])
#     end

#     # nodes gen
#     thermal_nodes = fill(0.0,nbus)
#     hydro_nodes = fill(0.0,nbus)

#     idxhyd = idx_hydro(data)

#     for i in 1:length(data["powersystem"]["gen"])
#         bus_i = data["powersystem"]["gen"]["$i"]["gen_bus"]
#         if i in idxhyd
#             hydro_nodes[bus_i] += Statistics.quantile([results[:simulations][s][t][:powersystem]["solution"]["gen"]["$i"]["pg"] for s in 1:length(results[:simulations])], quant)*data["powersystem"]["gen"]["$i"]["mbase"]
#         else
#             thermal_nodes[bus_i] += Statistics.quantile([results[:simulations][s][t][:powersystem]["solution"]["gen"]["$i"]["pg"] for s in 1:length(results[:simulations])], quant)*data["powersystem"]["gen"]["$i"]["mbase"]
#         end
#     end

#     # nodes loads
#     load_nodes = fill(0.0,nbus)

#     for i in 1:length(data["powersystem"]["load"])
#         bus_i = data["powersystem"]["load"]["$i"]["load_bus"]
#         load_nodes[bus_i] += data["powersystem"]["load"]["$i"]["pd"]*data["powersystem"]["baseMVA"]
#     end

#     # number of nodes
#     num_nodes = nbus+ sum(load_nodes .> 0)+sum(hydro_nodes .> 0)+sum(thermal_nodes .> 0)

#     # node size
#     nodesize = fill(0.1,num_nodes)

#     # nodes membership (3: Hydro, 2: Thermal)
#     membership = fill(1,num_nodes)

#     # node label
#     if node_label
#         nodelabel = [1:nbus;fill("",sum(hydro_nodes .> 0));fill("",sum(thermal_nodes .> 0));fill("",sum(load_nodes .> 0))]
#     else
#         nodelabel = nothing
#     end
#     # create nodes
#     for bus_i in 1:length(data["powersystem"]["bus"])
#         if hydro_nodes[bus_i] > 0
#             add_vertex!(g)
#             nNodes +=1
#             add_edge!(g, nNodes, bus_i)
#             membership[nNodes] = 3
#             nodesize[nNodes] = log(hydro_nodes[bus_i])
#         end

#         if thermal_nodes[bus_i] > 0
#             add_vertex!(g)
#             nNodes +=1
#             add_edge!(g, nNodes, bus_i)
#             membership[nNodes] = 2
#             nodesize[nNodes] = log(thermal_nodes[bus_i])
#         end

#         if load_nodes[bus_i] > 0
#             add_vertex!(g)
#             nNodes +=1
#             add_edge!(g, nNodes, bus_i)
#             membership[nNodes] = 4
#             nodesize[nNodes] = log(load_nodes[bus_i])
#         end

#     end

#     for n in 1:num_nodes
#         nodesize[n] = nodesize[n] == 0 ? unique(sort(nodesize))[2] : nodesize[n]
#     end

#     nodecolor = ["black", "red", "blue", "orange"]

#     # membership color
#     nodefillc = nodecolor[Int64.(membership)]

#     gplot(g, nodefillc=nodefillc, nodesize=nodesize, nodelabel=nodelabel,nodelabeldist=nodelabeldist)
# end

# """ Plot Grid Dispatched Power"""
# function plot_grid_dispatched(results::Dict;seed=1111,quant::Float64=0.5,size_fig = [15cm, 15cm],node_label=false,nodelabeldist=4.5)
#     duration=length(results[:simulations][1])
#     roll(fps=1, duration=duration) do t, dt
#         Random.seed!(seed)
#         plot_grid_dispatched_stage(results,Int64(t+1);quant=quant,size_fig = size_fig,node_label=node_label,nodelabeldist=nodelabeldist)
#     end
# end

"""
    HydroPowerModels.plot_aggregated_results(results::Dict;nc::Int=3)

Plot Aggregated Results. Figures are of aggregated quantities, but the methods used to aggregate were chosen in order to help analysis. For example: The final nodal price is an average of nodal prices weighted by the contribution of local loads to the total demand; Reservoir volume was grouped weighted by the amount of energy that could be produced by the stored water (as was the inflow of water).

Paremeter:
-   results: Simulation results.

"""
function plot_aggregated_results(results::Dict;nc::Int=3)
    plt_total = Array{Plots.Plot}(undef,20)
    nplots = 0
    nsim = length(results[:simulations])
    nstages = length(results[:simulations][1])

    # Thermal Generation
    ngen = length(results[:data][1]["powersystem"]["gen"])
    idxhyd = idx_hydro(results[:data][1])
    idxgen = setdiff(collect(1:ngen),idxhyd)
    baseMVA =  [results[:simulations][i][j][:powersystem]["solution"]["baseMVA"] for i=1:nsim, j=1:nstages]'
    scen_gen_all = [[results[:simulations][i][j][:powersystem]["solution"]["gen"]["$gen"]["pg"] for i=1:nsim, j=1:nstages]'.*baseMVA for gen =1:ngen]

    scen_gen = deepcopy(scen_gen_all[idxgen[1]])
    scen_gen .=0
    for gen in idxgen
        scen_gen = scen_gen .+ scen_gen_all[gen]
    end
    plt = plotscenarios(scen_gen, title  = "Thermal Generation",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
    plt_total[nplots+1] = plt
    nplots += 1

    # circuit MVA
    baseMVA = results[:data][1]["powersystem"]["baseMVA"]

    # nodes loads
    nbus = length(results[:data][1]["powersystem"]["bus"])
    load_nodes = fill(0.0,nstages,nbus)
    for t in 1:nstages
        data = results[:data][min(t,length(results[:data]))]
        for i in 1:length(data["powersystem"]["load"])
            bus_i = data["powersystem"]["load"]["$i"]["load_bus"]
            load_nodes[t,bus_i] += data["powersystem"]["load"]["$i"]["pd"]*data["powersystem"]["baseMVA"]
        end
    end

    # Nodal price
    idxbus = collect(1:nbus)
    try
    scen_pld_all = convert(Array{Array{Float64,2},1},[[-results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["lam_kcl_r"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])/baseMVA

    scen_pld = deepcopy(scen_pld_all[idxbus[1]])
    scen_pld .=0
    for bus in idxbus
        scen_pld = scen_pld .+ scen_pld_all[bus].*hcat(fill(load_nodes[:,bus],nsim)...)
    end
    for t=1:nstages
        scen_pld[t,:] /= sum(load_nodes[t,bus] for bus in idxbus)
    end
    plt = plotscenarios(scen_pld, title  = "Load Weighted Average Nodal price ",
                ylabel = "\$/MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
    plt_total[nplots+1] = plt
    nplots += 1
    catch
    end
    # Deficit
    try
    scen_def_all = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:powersystem]["solution"]["bus"]["$bus"]["deficit"] for i=1:nsim, j=1:nstages]' for bus =1:nbus])
    scen_def = deepcopy(scen_def_all[idxbus[1]])
    scen_def .=0
    for bus in idxbus
        scen_def = scen_def .+ scen_def_all[bus]
    end
    plt = plotscenarios(scen_def.*baseMVA, title  = "Deficit",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
    plt_total[nplots+1] = plt
    nplots += 1
    catch
    end

    # Hydro Generation
    scen_gen = deepcopy(scen_gen_all[idxhyd[1]])
    scen_gen .=0
    for gen in idxhyd
        scen_gen = scen_gen .+ scen_gen_all[gen]
    end
    plt = plotscenarios(scen_gen, title  = "Hydro Generation",
                ylabel = "MW",
                xlabel = "Stages",
                bottom_margin = 10mm,
                right_margin = 10mm,
                left_margin = 10mm
                )
    plt_total[nplots+1] = plt
    nplots += 1

    # Reservoir Outflow
    nHyd = results[:data][1]["hydro"]["nHyd"]
    scen_turn_all = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:outflow][res] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    scen_turn = deepcopy(scen_turn_all[1])
    scen_turn .=0
    for res = 1:nHyd
        scen_turn = scen_turn .+ scen_turn_all[res]
    end

    plt = plotscenarios(scen_turn, title  = "Reservoir Outflow",
                    ylabel = "m³/s",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
    plt_total[nplots+1] = plt
    nplots += 1

    # Reservoir Spill
    scen_spill_all = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:spill][res] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])./(0.0036*results[:params]["stage_hours"])

    scen_spill = deepcopy(scen_spill_all[1])
    scen_spill .=0
    for res = 1:nHyd
        scen_spill = scen_spill .+ scen_spill_all[res]
    end

    plt = plotscenarios(scen_spill, title  = "Reservoir Spill",
                    ylabel = "m³/s",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
    plt_total[nplots+1] = plt
    nplots += 1

    # Reservoir Volume

    water_energy!(results[:data][1])

    scen_voume_all = convert(Array{Array{Float64,2},1},[[results[:simulations][i][j][:reservoirs][:reservoir][res].out for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    scen_voume = deepcopy(scen_voume_all[1])
    scen_voume .=0
    for res = 1:nHyd
        scen_voume = scen_voume .+ (scen_voume_all[res].*results[:data][1]["hydro"]["Hydrogenerators"][res]["water_energy"])./(0.0036*results[:params]["stage_hours"])
    end

    plt = plotscenarios(scen_voume, title  = "Volume Reservoir",
                    ylabel = "MW",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
    plt_total[nplots+1] = plt
    nplots += 1

    # Inflows

    scen_inflows_all = convert(Array{Array{Float64,2},1},[[results[:data][1]["hydro"]["Hydrogenerators"][res]["inflow"][cidx(j,results[:data][1]["hydro"]["size_inflow"][1]),results[:simulations][i][j][:noise_term]] for i=1:nsim, j=1:nstages]' for res = 1:results[:data][1]["hydro"]["nHyd"]])

    scen_inflows = deepcopy(scen_inflows_all[1])
    scen_inflows .=0
    for res = 1:nHyd
        scen_inflows = scen_inflows .+ (scen_inflows_all[res].*results[:data][1]["hydro"]["Hydrogenerators"][res]["water_energy"])
    end

    plt = plotscenarios(scen_inflows, title  = "Inflows",
                    ylabel = "MW",
                    xlabel = "Stages",
                    bottom_margin = 10mm,
                    right_margin = 10mm,
                    left_margin = 10mm
                    )
    plt_total[nplots+1] = plt
    nplots += 1

    if mod(nplots,nc) > 0 && floor(Int,nplots/nc) > 0
        l = @layout [ Plots.grid(floor(Int,nplots/nc),nc);  Plots.grid(1,mod(nplots,nc))]
        nlines = floor(Int,nplots/nc)+1
        l.heights = Plots.grid(2,1,heights=[floor(Int,nplots/nc)/nlines;1/nlines]).heights
    elseif floor(Int,nplots/nc) > 0
        l = @layout Plots.grid(floor(Int,nplots/nc),nc)
        nlines = floor(Int,nplots/nc)
        l.heights = Plots.grid(nlines,1,heights=[1/nlines for n = 1:nlines]).heights
    else
        l = @layout Plots.grid(1,mod(nplots,nc))
        nlines = 1
        l.heights = Plots.grid(1,1,heights=[1]).heights
    end

     return plot(plt_total[1:nplots]...,layout=l,size = (4*400, 500*ceil(Int,nplots/nc)),legend=false)
end

"""
    HydroPowerModels.plot_bound(m)

Plots the SDDP outer bound per iteration.
"""
function plot_bound(m)
    niter = length(m.policygraph.most_recent_training_results.log)
    val = round.([m.policygraph.most_recent_training_results.log[iter].bound for iter in 1:niter],digits =5)

    xticks = unique!([collect(1:Int(floor(niter/4)):niter);niter])
    plot(val,
        xticks = (xticks, [string(i) for  i in xticks]),
        label = "Bound",
        title = "Bound x iterations",
        bottom_margin = 10mm,
        right_margin = 10mm,
        left_margin = 10mm
        )
end


text:
"""
    mesh_connectivity_2D(cells, points)

Compute connectivity of 2D unstructured mesh
"""
function mesh_connectivity_2D(cells, points::AbstractMatrix{T}) where {T<:Real}
    cellid = extract_cell(cells)

    facePoints, faceCells = mesh_face_connectivity_2D(cellid)
    cellNeighbors = mesh_cell_neighbor_2D(cellid, facePoints, faceCells)
    cellFaces = mesh_cell_face(cellid, faceCells)
    cellType = mesh_cell_type(cellNeighbors)
    cellArea = mesh_cell_area_2D(points, cellid)
    cellCenter = mesh_cell_center(points, cellid)
    cellNormals = mesh_cell_normals_2D(points, cellid, cellCenter)
    faceCenter = mesh_face_center(points, facePoints)
    faceType = mesh_face_type(faceCells, cellType)
    faceArea = mesh_face_area_2D(points, facePoints)

    return cellid,
        cellType,
        cellNeighbors,
        cellFaces,
        cellCenter,
        cellArea,
        cellNormals,
        facePoints,
        faceCells,
        faceCenter,
        faceType,
        faceArea
end

# ------------------------------------------------------------
# Face-indexed connectivity
# ------------------------------------------------------------

"""
    mesh_face_connectivity_2D(cells::AbstractArray{<:Integer,2})

Compute affliated points and neighbor cells of faces
"""
function mesh_face_connectivity_2D(cells::T) where {T<:AbstractArray{<:Integer,2}}
    nNodesPerCell = size(cells, 2)
    nCells = size(cells, 1)
    nEdgesMax = nNodesPerCell * nCells

    tmpEdgeNodes = -ones(Int, nEdgesMax, 2)
    tmpEdgeCells = -ones(Int, nEdgesMax, 2)

    counter = 0
    for i = 1:nCells, k = 1:nNodesPerCell
        isNewEdge = true
        for j = 1:counter
            if tmpEdgeNodes[j, :] == [cells[i, k], cells[i, k%nNodesPerCell+1]] ||
               tmpEdgeNodes[j, :] == [cells[i, k%nNodesPerCell+1], cells[i, k]]
                isNewEdge = false
                tmpEdgeCells[j, 2] = i
            end
        end
        if isNewEdge
            counter += 1
            tmpEdgeNodes[counter, 1] = cells[i, k]
            tmpEdgeNodes[counter, 2] = cells[i, k%nNodesPerCell+1]
            tmpEdgeCells[counter, 1] = i
        end
    end

    nEdges = counter
    edgeNodes = tmpEdgeNodes[1:nEdges, :]
    edgeCells = tmpEdgeCells[1:nEdges, :]

    return edgeNodes, edgeCells
end

"""
    mesh_face_center(
        nodes::X,
        edgeNodes::Y,
    ) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Integer,2}}

Compute central points of cell faces
"""
function mesh_face_center(
    nodes::X,
    edgeNodes::Y,
) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Integer,2}}

    edgeCenter = zeros(size(edgeNodes, 1), size(nodes, 2))
    for i in axes(edgeCenter, 1)
        pids = edgeNodes[i, :]

        tmp = zeros(size(nodes, 2))
        for i in eachindex(pids)
            tmp .+= nodes[pids[i], :]
        end
        tmp ./= length(pids)

        @. edgeCenter[i, :] = tmp
    end

    return edgeCenter

end

"""
    function mesh_face_type(
        edgeCells::X,
        cellType::Y,
    ) where {X<:AbstractArray{<:Integer,2},Y<:AbstractArray{<:Integer,1}}

Compute type of faces

- 0: inner
- 1: boundary
"""
function mesh_face_type(
    edgeCells::X,
    cellType::Y,
) where {X<:AbstractArray{<:Integer,2},Y<:AbstractArray{<:Integer,1}}
    edgeType = zeros(eltype(edgeCells), size(edgeCells, 1))

    for i in axes(edgeCells, 1)
        i1 = edgeCells[i, 1]
        i2 = edgeCells[i, 2]
        if cellType[i1] == 0 && cellType[i2] == 0
            edgeType[i] = 0
        else
            edgeType[i] = 1
        end
    end

    return edgeType
end

"""
    mesh_face_area_2D(points, facePoints)

Compute area of faces
"""
function mesh_face_area_2D(points, facePoints)
    faceArea = zeros(size(facePoints, 1))

    for i in eachindex(faceArea)
        faceArea[i] = norm(points[facePoints[i, 1], 1:2] .- points[facePoints[i, 2], 1:2])
    end

    return faceArea
end

# ------------------------------------------------------------
# Cell-indexed connectivity
# ------------------------------------------------------------

"""
    mesh_cell_neighbor_2D(cells::AbstractMatrix{T}, edgeNodes::AbstractMatrix{T}, edgeCells::AbstractMatrix{T}) where {T<:Integer}

Compute neighbor cells of cells
"""
function mesh_cell_neighbor_2D(
    cells::AbstractMatrix{T},
    edgeNodes::AbstractMatrix{T},
    edgeCells::AbstractMatrix{T},
) where {T<:Integer}
    nCells = size(cells, 1)
    nNodesPerCell = size(cells, 2)
    nEdges = size(edgeNodes, 1)

    cellNeighbors = -ones(Int, nCells, nNodesPerCell)
    for i = 1:nCells, k = 1:nNodesPerCell, j = 1:nEdges
        if length(intersect(edgeNodes[j, :], [cells[i, k], cells[i, k%nNodesPerCell+1]])) == 2 # shared face
            idx = findall(x->x!=i, edgeCells[j, :]) |> first # idx takes value 1 or 2
            cellNeighbors[i, k] = edgeCells[j, idx]
        end
    end

    # cellNeighbors are already regularized, i.e.,
    # face 1 -> neighbor cell 1
    # face 2 -> neighbor cell 2
    # face 3 -> neighbor cell 3
    # so we can get rid of this part
    #=
    cellNeighbors_regular = -ones(Int, nCells, nNodesPerCell)
    for i = 1:nCells
        cids = cellNeighbors[i, :]

        fpids = [cells[i, j:j+1] for j = 1:nNodesPerCell-1]
        push!(fpids, [cells[i, nNodesPerCell], cells[i, 1]])

        for j = 1:nNodesPerCell
            fpid = fpids[j]

            for k in eachindex(cids)
                if cids[k] != -1 && length(intersect(fpid, cells[cids[k], :])) == 2
                    cellNeighbors_regular[i, j] = cids[k]
                end
            end
        end
    end
    =#

    return cellNeighbors
end

"""
    mesh_cell_type(cellNeighbors::T) where {T<:AbstractArray{<:Integer,2}}

Compute types of elements

- 0: inner
- 1: boundary
"""
function mesh_cell_type(cellNeighbors::T) where {T<:AbstractArray{<:Integer,2}}
    cellid = zeros(eltype(cellNeighbors), size(cellNeighbors, 1))
    for i in axes(cellNeighbors, 1)
        if -1 in cellNeighbors[i, :]
            cellid[i] = 1
        end
    end

    return cellid
end

"""
    mesh_area_2D(nodes::AbstractArray{<:AbstractFloat,2}, cells::AbstractArray{<:Int,2})

Compute areas of 2D elements
"""
function mesh_cell_area_2D(
    nodes::X,
    cells::Y,
) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Integer,2}}

    ΔS = zeros(size(cells, 1))

    if size(cells, 2) == 3 # triangular mesh
        for i in eachindex(ΔS)
            ΔS[i] = abs(
                (
                    nodes[cells[i, 1], 1] *
                    (nodes[cells[i, 2], 2] - nodes[cells[i, 3], 2]) +
                    nodes[cells[i, 2], 1] *
                    (nodes[cells[i, 3], 2] - nodes[cells[i, 1], 2]) +
                    nodes[cells[i, 3], 1] * (nodes[cells[i, 1], 2] - nodes[cells[i, 2], 2])
                ) / 2,
            )
        end
    elseif size(cells, 2) == 4 # quadrilateral mesh
        for i in eachindex(ΔS)
            d1 = [
                nodes[cells[i, 1], 1] - nodes[cells[i, 2], 1],
                nodes[cells[i, 1], 2] - nodes[cells[i, 2], 2],
            ]
            d2 = [
                nodes[cells[i, 2], 1] - nodes[cells[i, 3], 1],
                nodes[cells[i, 2], 2] - nodes[cells[i, 3], 2],
            ]
            d3 = [
                nodes[cells[i, 3], 1] - nodes[cells[i, 4], 1],
                nodes[cells[i, 3], 2] - nodes[cells[i, 4], 2],
            ]
            d4 = [
                nodes[cells[i, 4], 1] - nodes[cells[i, 1], 1],
                nodes[cells[i, 4], 2] - nodes[cells[i, 1], 2],
            ]

            a = sqrt(d1[1]^2 + d1[2]^2)
            b = sqrt(d2[1]^2 + d2[2]^2)
            c = sqrt(d3[1]^2 + d3[2]^2)
            d = sqrt(d4[1]^2 + d4[2]^2)
            T = 0.5 * (a + b + c + d)

            alpha = acos((d4[1] * d1[1] + d4[2] * d1[2]) / (a * d))
            beta = acos((d2[1] * d3[1] + d2[2] * d3[2]) / (b * c))

            ΔS[i] = sqrt(
                (T - a) * (T - b) * (T - c) * (T - d) -
                a * b * c * d * cos(0.5 * (alpha + beta)) * cos(0.5 * (alpha + beta)),
            )
        end
    end

    return ΔS

end

"""
    mesh_center_2D(nodes::AbstractArray{<:AbstractFloat,2}, cells::AbstractArray{<:Integer,2})

Compute central points of elements
"""
function mesh_cell_center(
    nodes::X,
    cells::Y,
) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Integer,2}}

    cellMidPoints = zeros(size(cells, 1), size(nodes, 2))
    for i in axes(cellMidPoints, 1) # nCells
        for j in axes(cells, 2) # nNodesPerCell
            cellMidPoints[i, :] .+= nodes[cells[i, j], :]
        end
    end
    cellMidPoints ./= size(cells, 2)

    return cellMidPoints

end

"""
    mesh_cell_face(
        cells::X,
        edgeCells::Y,
    ) where {X<:AbstractArray{<:Integer,2},Y<:AbstractArray{<:Integer,2}}

Compute surrounding faces of cell
"""
function mesh_cell_face(
    cells::X,
    edgeCells::Y,
) where {X<:AbstractArray{<:Integer,2},Y<:AbstractArray{<:Integer,2}}
    ncell = size(cells, 1)
    vv = [Int[] for i = 1:ncell]
    for i in axes(edgeCells, 1)
        if edgeCells[i, 1] != -1
            push!(vv[edgeCells[i, 1]], i)
        end
        if edgeCells[i, 2] != -1
            push!(vv[edgeCells[i, 2]], i)
        end
    end

    cellEdges = zero(cells)
    for i in axes(cellEdges, 1)
        cellEdges[i, :] .= vv[i]
    end

    return cellEdges
end

"""
    mesh_cell_normals_2D(points, cells, cellCenter)

Compute unit normal vectors of cells
"""
function mesh_cell_normals_2D(points, cells, cellCenter)
    ncell = size(cells, 1)
    np = size(cells, 2)
    cell_normal = zeros(ncell, np, 2)

    for i in 1:ncell
        pids = [cells[i, j:j+1] for j = 1:np-1]
        push!(pids, [cells[i, np], cells[i, 1]])

        for j = 1:np
            cell_normal[i, j, :] .= unit_normal(points[pids[j][1], :], points[pids[j][2], :])
        end

        p = [(points[pids[j][1], :] .+ points[pids[j][2], :]) ./ 2 for j = 1:np]

        for j = 1:np
            if dot(cell_normal[i, j, :], p[j][1:2] - cellCenter[i, 1:2]) < 0
                cell_normal[i, j, :] .= -cell_normal[i, j, :]
            end
        end
    end

    return cell_normal
end


text:
# Autogenerated wrapper script for tblis_jll for x86_64-apple-darwin
export tblis, tci

using Hwloc_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("tblis")
JLLWrappers.@declare_library_product(tblis, "@rpath/libtblis.0.dylib")
JLLWrappers.@declare_library_product(tci, "@rpath/libtci.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Hwloc_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        tblis,
        "lib/libtblis.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        tci,
        "lib/libtci.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Pkg.add("JSON") ... an external library http://docs.julialang.org/en/latest/packages/packagelist/
using JSON

sample = Dict()
sample["blue"] = [1, 2]
sample["ocean"] = "water"

@show sample jsonstring = json(sample)
@show jsonobj = JSON.parse(jsonstring)

@assert jsonstring == "{\"ocean\":\"water\",\"blue\":[1,2]}"
@assert jsonobj == Dict("ocean" => "water", "blue" => [1, 2])
@assert typeof(jsonobj) == Dict{String, Any}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CodeTracking, Markdown

"""
    does_thmhold(hold)

Check if `hold` is true.

# Examples

```jldoctest
julia> does_thmhold(true)
Theorem holds! 😀
```

```jldoctest
julia> does_thmhold(false)
Theorem does not holds! 😱
```
"""
function does_thmhold(hold)
if(hold == true)
    println("Theorem holds! 😀")
else
    println("Theorem does not holds! 😱")
end
end

"""
    @code_md(expr)

Return the source code for evaluating `expr` in markdown form.
"""
macro code_md(expr)
    quote
        Markdown.parse(
        """
        ```julia
        $(@code_string $expr)
        ```
        """)
    end
end


text:
# implementation of the self-consistent-field (SCF) method

```
Form the density matrix R from the occupied columns of U. C contains (as columns) the MO
coefficients and the R-matrix is formed from the first nocc columns.

C - The eigenvectors are stored as columns of length m
R - The charge and bond order matrix generated is stored in R by columns (m x m)
m - The dimension of the matrices (basis set size)
nocc - The number of occupied orbitals (half of the number of eletrons in the closed-shell case)

```
function scfr(C, R, m::Int64, nocc::Int64)  # C and R : Matrix
  for i in 1:m
    for j in 1:i
      sum = 0
      for k in 1:nocc
        kk = m*(k-1)
        ik = kk + i
        jk = kk + j
        sum = sum + C[ik]*C[jk]
      end
      ij = m*(j - 1) + i
      ji = m*(i - 1) + j
      R[ij] = sum;
      R[ji] = sum;
    end
  end
end

export scfr


text:
export supported_platforms

"""
    downloads_dir(postfix::String = "")

Builds a path relative to the `downloads_cache`.
"""
function downloads_dir(postfix::String = "")
    global downloads_cache
    return joinpath(downloads_cache, postfix)
end

"""
    rootfs_dir(postfix::String = "")

Builds a path relative to the `rootfs_cache`.
"""
function rootfs_dir(postfix::String = "")
    global rootfs_cache
    return joinpath(rootfs_cache, postfix)
end

"""
    shards_dir(postfix::String = "")

Builds a path relative to the `shards_cache`.
"""
function shards_dir(postfix::String = "")
    global shards_cache
    return joinpath(shards_cache, postfix)
end

shard_path_squashfs(shard_name) = downloads_dir("rootfs-$(shard_name).squashfs")
rootfs_path_squashfs() = shard_path_squashfs("base")

function ccache_dir()
    global ccache_override

    if ccache_override != ""
        return ccache_override
    else
        return shards_dir("ccache")
    end
end

function sandbox_path()
    global sandbox_override

    if sandbox_override != ""
        return sandbox_override
    else
        return rootfs_dir("sandbox")
    end
end

"""
    get_shard_url(target::String = "base"; squashfs::Bool = use_squashfs)

Returns the URL from which a rootfs image (tarball/squashfs) can be downloaded
"""
function get_shard_url(target::String = "base"; squashfs::Bool = use_squashfs)
    # These constants are what should be updated for a new rootfs build:
    rootfs_urlroot = "https://julialangmirror-s3.julialang.org/binarybuilder"
    rootfs_version = "2018-05-26"

    shard_name = "rootfs-$(target)"
    ext = squashfs ? "squashfs" : "tar.gz"

    return "$(rootfs_urlroot)-rootfs-$(target)-$(rootfs_version).$(ext)"
end

"""
    get_shard_url(target::String = "base"; squashfs::Bool = use_squashfs)

Returns the sha256 hash for a rootfs image (tarball/squashfs).
"""
function get_shard_hash(triplet::String = "base"; squashfs::Bool = use_squashfs)
    # These constants are what should be updated for a new rootfs build.
    # You can get these dictionaries spat out by running `make buildall && make push-rootfs`
    # in `julia-docker/crossbuild`
    squashfs_hashes = Dict(
        "aarch64-linux-gnu" => "05ee1834cc5d5d4f5d2864b9d47aea894ccfe984f151f6cef802c707aae14aff",
        "aarch64-linux-musl" => "631c6acd70a32b07d6c202369bea24ece17a1a09f85257f3a88a3dc6cd9ff281",
        "arm-linux-gnueabihf" => "e74ec696b2396267f4ca9b4ac04abe94d092f994435234ac4503d13a20aaaa04",
        "arm-linux-musleabihf" => "1c9d3c6e3f082459b70ea49370641ee8f00c5b058732d140bf10e445d65fbca5",
        "base" => "154a39efdefe67d680184869a15c9322594e93ce4c03621c81a35b17bd392e4b",
        "i686-linux-gnu" => "75cad3c1e5566edd6e162aadee5e54c45db49ae984abeed165a3c0fe36463eb6",
        "i686-linux-musl" => "7488e2c6e810b6c738d3af07a4c43e49d0018dd5ab613bc838de9692b8e63543",
        "i686-w64-mingw32" => "7eeeee6959f5d411941b17539b0bee07d2a1523930a1ec6e17dc45d08ae1a940",
        "powerpc64le-linux-gnu" => "36bb3449e5ec810905f08f68f9ea68a3f082e4dcd63505c469121108b8d2d47c",
        "x86_64-apple-darwin14" => "bc33d23a7c5974c6c2d310ef01eb9098cbf49618a07f6488bae89b8c80d1105b",
        "x86_64-linux-gnu" => "60d42526c76c19066d67b716173ac64fe29f93222e36c295c04e5ee5e9ab317a",
        "x86_64-linux-musl" => "725a3fbe875b068b3ee4a2536ed1c7bfda4a2bde1ff8aaa3422aa9bbd51e5556",
        "x86_64-unknown-freebsd11.1" => "538c077fab1977514aa49de54e4b6ddcde6ce036668d68cbc3cc14cc164c710d",
        "x86_64-w64-mingw32" => "c4c8015d8f8ed84ae5162902072093037ff9364d372ba03e33fcaf1938c9f6af",
    )
    tarball_hashes = Dict(
        "aarch64-linux-gnu" => "4b6a9d4fce1dfc1c4bfc6a439ea377225e493b29b4cbc4b5f312e9841611156f",
        "aarch64-linux-musl" => "3b0b832c7a498b6c778045038efedcc8f06778f0d653e093fc29d6fa4feb3f58",
        "arm-linux-gnueabihf" => "2b4f224c3c7f43b4ecdc63edd37ab73cd42411b9288bf3cc0d01599ea53541e1",
        "arm-linux-musleabihf" => "04b2a8a26721fc07045a67e4f867f715fe48a7d2e0e29636ece62e2e243fd3d9",
        "base" => "ecc87a85e79d3a64e853634678ce8e182a12367fdc062dbd046164f2ba0dbd3d",
        "i686-linux-gnu" => "848726ca7acc7dc23f5b43f9992d8235761b694260ae4336ff177b6952e3f759",
        "i686-linux-musl" => "b7a429b48cafd5fada36af52b24492f08ff67579211ee7c6854489ee214858df",
        "i686-w64-mingw32" => "7c64955cc3d1bc70bdeed05cdd0d4df140d2ec68ebb3e4c907db74f7370d1616",
        "powerpc64le-linux-gnu" => "34a1f375464d9a7269a8ff3eefb95e8f4c2711ab851dec38c959685adb7724b9",
        "x86_64-apple-darwin14" => "896a8b7ff4c423c189448041cbda6b7124076893bc66158eb78b79541f8647db",
        "x86_64-linux-gnu" => "757a695d8cbe2d62ba3b75f9a0ff568d172bf63db8ff6fde3ad757b86aee1f8e",
        "x86_64-linux-musl" => "f7b23fa8d799f553679d5bd5cb8149975e1a930740cadb2ceb721d360d490482",
        "x86_64-unknown-freebsd11.1" => "34ac437cea654fc92b64fea713673ba2dd47c721e0cd30df037d74a48fff32a1",
        "x86_64-w64-mingw32" => "35b5bbbab349e14590454dfec32d9b4bff1ea150d7cbb7d1d7ce83622876ee31",
    )
	if squashfs
        return squashfs_hashes[triplet]
    else
        return tarball_hashes[triplet]
    end
end

"""
    supported_platforms()

Return the list of supported platforms as an array of `Platform`s.  These are the platforms we
officially support building for, if you see a mapping in `get_shard_hash()` that isn't
represented here, it's probably because that platform is still considered "in beta".
"""
function supported_platforms()
    return [
        # glibc Linuces
        Linux(:i686),
        Linux(:x86_64),
        Linux(:aarch64),
        Linux(:armv7l),
        Linux(:powerpc64le),

        # musl Linuces
        Linux(:i686, :musl),
        Linux(:x86_64, :musl),
        Linux(:aarch64, :musl),
        Linux(:armv7l, :musl),

        # BSDs
        MacOS(:x86_64),
        FreeBSD(:x86_64),

        # Windows
        Windows(:i686),
        Windows(:x86_64),
    ]
end

"""
    getuid()

Wrapper around libc's `getuid()` function
"""
function getuid()
    return ccall(:getuid, Cint, ())
end

# Note: produce these values by #including squashfs_fs.h from linux in Cxx.jl
# and running the indicated command
const offsetof_id_table_start = 0x30    # offsetof(struct suqashfs_super_block, id_table_start)
const offsetof_no_ids = 0x1a            # offsetof(struct suqashfs_super_block, no_ids)

# From squashfs_fs.h
const SQUASHFS_COMPRESSED_BIT = UInt16(1) << 15
const SQUASHFS_MAGIC = 0x73717368

"""
    rewrite_squashfs_uids(path, new_uid)

In order for the sandbox to work well, we need to have the uids of the squashfs
images match the uid of the current unpriviledged user. Unfortunately there is
no mount-time option to do this for us. However, fortunately, squashfs is simple
enough that if the id table is uncompressed, we can just manually patch the uids
to be what we need. This functions performs this operation, by rewriting all
uids/gids to new_uid.
"""
function rewrite_squashfs_uids(path, new_uid)
    open(path, "r+") do file
       # Check magic
       if read(file, UInt32) != SQUASHFS_MAGIC
           error("`$path` is not a squashfs file")
       end
       # Check that the image contains only one id (which we will rewrite)
       seek(file, offsetof_no_ids)
       if read(file, UInt16) != 1
           error("`$path` uses more than one uid/gid")
       end
       # Find the index table
       seek(file, offsetof_id_table_start)
       offset = read(file, UInt64)
       seek(file, offset)
       # Find the correct metdata block
       index = read(file, UInt64)
       seek(file, index)
       # Read the metadata block
       size = read(file, UInt16)
       # Make sure it's uncompressed (yes, I know that flag is terribly
       # named - it indicates that the data is uncompressed)
       if ((size & SQUASHFS_COMPRESSED_BIT) == 0)
           error("Metadata block is compressed")
       end
       p = position(file)
       uid = read(file, UInt32)
       if uid != 0
           error("Expected all uids to be 0 (not $(uid)) inside the image")
       end
       seek(file, p)
       write(file, UInt32(new_uid))
    end
    return nothing
end

_sudo_cmd = nothing
function sudo_cmd()
    global _sudo_cmd

    # Use cached value if we've already run this
    if _sudo_cmd != nothing
        return _sudo_cmd
    end

    if getuid() == 0
        # If we're already root, don't use any kind of sudo program
        _sudo_cmd = ``
    elseif success(`sudo -V`)
        # If `sudo` is available, use that
        _sudo_cmd = `sudo`
    else
        # Fall back to `su` if all else fails
        _sudo_cmd = `su root -c`
    end
    return _sudo_cmd
end

"""
    update_rootfs(triplets::Vector{AbstractString};
                  automatic::Bool = automatic_apple, verbose::Bool = true,
                  squashfs::Bool = use_squashfs)

Updates the stored rootfs containing all cross-compilers and other compilation
machinery for the given triplets.  If `automatic` is set, when downloading Apple
SDKs, you will automatically accept the Apple license agreement and download the
macOS SDK for usage in targeting macOS.  See the help for `download_osx_sdk()`
for more details on this.
"""
function update_rootfs(triplets::Vector{S}; automatic::Bool = automatic_apple,
                       verbose::Bool = false, squashfs::Bool = use_squashfs,
                       mount::Bool=use_squashfs && Compat.Sys.islinux()) where {S <: AbstractString}
    # Check to make sure we have the latest version of both the base and the
    # given shard downloaded properly, and extracted if it's not a squashfs.
    for shard_name in ["base", triplets...]
        url = get_shard_url(shard_name; squashfs=squashfs)
        hash = get_shard_hash(shard_name; squashfs=squashfs)

        if shard_name == "base"
            dest_dir = rootfs_dir()
        else
            dest_dir = shards_dir(shard_name)
        end

        if squashfs
            squashfs_path = shard_path_squashfs(shard_name)
            file_existed = isfile(squashfs_path)

            # If squashfs, verify/redownload the squashfs image
            if !(download_verify(
                url,
                hash,
                squashfs_path;
                verbose = verbose,
                force = true
            ) && file_existed)
                if mount && Compat.Sys.islinux()
                    # Unmount the mountpoint. It may point to a previous version
                    # of the file. Also, we're about to mess with it
                    unmount_shard(dest_dir)

                    # Patch squashfs files for current uid
                    rewrite_squashfs_uids(squashfs_path, getuid())

                    # Touch SHA file to prevent re-verifcation from blowing the
                    # fs away
                    touch(string(squashfs_path,".sha256"))
                end
            end

            # Then mount it, if it hasn't already been mounted:
            if mount && Compat.Sys.islinux() && !success(`mountpoint $(dest_dir)`)
                mkpath(dest_dir)
                Compat.@info("Mounting $(squashfs_path) to $(dest_dir)")
                run(`$(sudo_cmd()) mount $(squashfs_path) $(dest_dir) -o ro,loop`)
            end
        else
            # If it has been mounted previously, unmount here
            unmount_shard(dest_dir)

             # If tarball, verify/redownload/reunpack the tarball
            download_verify_unpack(
                url,
                hash,
                dest_dir;
                tarball_path = downloads_dir("rootfs-$(shard_name).tar.gz"),
                verbose = verbose,
                force = true,
            )
        end
    end

    # If we're targeting the macOS SDK here, make sure it's ready to go.
    if any(triplets .== "x86_64-apple-darwin14")
        download_osx_sdk(;automatic=automatic, verbose=verbose)
    end
end

# Helper for when you're only asking for a single triplet
update_rootfs(triplet::AbstractString; kwargs...) = update_rootfs([triplet]; kwargs...)

# Helper to unmount any shards that may be mounted, so as not to exhaust the number of loopback devices
function unmount_shard(dest_dir::AbstractString; fail_on_error::Bool = false)
    # This function only matters on Linux
    @static if !Compat.Sys.islinux()
        return
    end

    if success(`mountpoint $(dest_dir)`)
        Compat.@info("Unmounting $(dest_dir)`")
        try
            run(pipeline(cmd, stdin=devnull, stdout=devnull, stderr=devnull))
        catch e
            # By default we don't error out if this unmounting fails
            if fail_on_error
                rethrow(e)
            end
        end
    end
    return
end

function unmount_all_shards(;fail_on_error::Bool = false)
    # This function only matters on Linux
    @static if !Compat.Sys.islinux()
        return
    end

    dest_dirs = shards_dir.(triplet.(supported_platforms()))
    push!(dest_dirs, rootfs_dir())

    for dest_dir in dest_dirs
        unmount_shard(dest_dir; fail_on_error = fail_on_error)
    end
    return
end

"""
    download_osx_sdk(;automatic::Bool = automatic_apple, verbose::Bool = false,
                      version::AbstractString = "10.10")

Apple restricts distribution and usage of the macOS SDK, a necessary component
to build software for macOS targets.  Please read the Apple and Xcode SDK
agreement for more information on the restrictions and legal terms you agree to
when using the SDK to build software for Apple operating systems:
https://images.apple.com/legal/sla/docs/xcode.pdf.

If `automatic` is set, this method will automatically agree to the Apple usage
terms and download the macOS SDK, enabling building for macOS.

To set this on an environment level, set the `BINARYBUILDER_AUTOMATIC_APPLE`
environment variable to `"true"`.
"""
function download_osx_sdk(;automatic::Bool = automatic_apple,
                           verbose::Bool = false,
                           version::AbstractString = "10.10")
    urlbase = "https://github.com/phracker/MacOSX-SDKs/releases/download/10.13"

    # Right now, we only support one version, but in the future, we may need
    # to support multiple macOS SDK versions
    sdk_metadata = Dict(
        "10.10" => (
            "$(urlbase)/MacOSX10.10.sdk.tar.xz",
            "4a08de46b8e96f6db7ad3202054e28d7b3d60a3d38cd56e61f08fb4863c488ce",
            shards_dir("MacOSX10.10.sdk"),
        ),
    )

    if !haskey(sdk_metadata, version)
        error("Unknown macOS version $(version); cannot download SDK!")
    end
    url, hash, dest = sdk_metadata[version]

    # If it already exists, return out
    if isdir(dest)
        if verbose
            Compat.@info("macOS SDK $(dest) already exists")
        end
        return
    end

    # Ask if we should download it, unless we're automated
    if !automatic
        if !isinteractive()
            msg = strip("""
            This is not an interactive Julia session, so we will not prompt you
            to download and install the macOS SDK, see the docstring for the
            `download_osx_sdk()` method for more details.
            """)
            Compat.@warn(msg)
            return
        end
        msg = strip("""
        macOS SDK not yet downloaded!

        Apple restricts distribution and usage of the macOS SDK, a necessary
        component to build software for macOS targets.  Please read the Apple
        and Xcode SDK agreement for more information on the restrictions and
        legal terms you agree to when using the SDK to build software for Apple
        operating systems: https://images.apple.com/legal/sla/docs/xcode.pdf.
        """)
        printstyled(msg, bold=true)
        println()
        while true
            print("Would you like to download and use the macOS SDK? [y/N]: ")
            answer = lowercase(strip(readline(stdin)))
            if answer == "y" || answer == "yes"
                break
            elseif answer == "n" || answer == "no"
                return
            else
                println("Unrecognized answer. Answer `y` or `n`.")
            end
        end
    else
        if verbose
            Compat.@info("Automatic macOS SDK install initiated")
        end
    end

    download_verify_unpack(
        url,
        hash,
        dest;
        tarball_path=downloads_dir(basename(url)),
        verbose=verbose,
        force=true
    )

    # These macOS tarballs have a nasty habit of nesting, let's fix that:
    dir_in_dir = joinpath(dest, basename(dest))
    if isdir(dir_in_dir)
        mv(dir_in_dir, "$(dirname(dir_in_dir))2")
        mv("$(dirname(dir_in_dir))2", dest; remove_destination=true)
    end
end


text:
function ContinuousScheduler(arg0::MeasurementBuilder, arg1::DatesSelector)
    return ContinuousScheduler((MeasurementBuilder, DatesSelector), arg0, arg1)
end

function generate(obj::ContinuousScheduler, arg0::List)
    return jcall(obj, "generate", SortedSet, (List,), arg0)
end


text:
using Plots;
using Printf;
gr();
anim = Animation();
f(x) = 1/2*cos(8pi*(x-1/2));

function forward(u_old::Array, i::Int64)
	return -(2u_old[ir[i]] + 3u_old[i] - 6u_old[il[i]] + u_old[il2[i]])C / 6 + (u_old[ir[i]] - 2u_old[i] + u_old[il[i]])C^2 / 2 - (u_old[ir[i]] - 3u_old[i] + 3u_old[il[i]] - u_old[il2[i]])C^3 / 6;
end

function backward(u_old::Array, i::Int64)
	return (2u_old[il[i]] + 3u_old[i] - 6u_old[ir[i]] + u_old[ir2[i]])*(-C) / 6 + (u_old[il[i]] - 2u_old[i] + u_old[ir[i]])*(-C)^2 / 2 + (u_old[il[i]] - 3u_old[i] + 3u_old[ir[i]] - u_old[ir2[i]])*(-C)^3 / 6;
end

function main()
#=********Numerical Set Up*********=#
	global c = 1.0;
	global nx = 101;
	global lx = 100;
	global x = range(0.0, stop=lx, length=nx);
	global dx = lx /(nx-1);
	global u1 = zeros(Float64, nx);
	global u2 = zeros(Float64, nx);
	@inbounds for i = 1:nx
		if 40.0 < x[i] < 60.0;
			u1[i] = f(x[i]);
		end
	end
	u2 = copy(u1);
	u = zeros(Float64, nx);
	global t = 0.0;
	global dt = 5e-1;
	global tlims = 300;
	global C = dt * c / dx;
	global ir = zeros(Int64, nx);
	global ir2 = zeros(Int64, nx);
	global il = zeros(Int64, nx);
	global il2 = zeros(Int64, nx);
	@inbounds for i = 1:nx
		ir[i] = i + 1;
		ir2[i] = i + 2;
		il[i] = i - 1;
		il2[i] = i - 2;
	end
	ir[end] = 1;
	ir2[end - 1] = 1;
	ir2[end] = 2;
	il[begin] = nx;
	il2[begin + 1] = nx;
	il2[begin] = nx - 1;
#=******Finish Numerical Set Up*************=#

	while t < tlims
		u_old1 = copy(u1);
		u_old2 = copy(u2);
		@inbounds for i = 1:nx
			u1[i] = u_old1[i] + forward(u_old1, i);
			u2[i] = u_old2[i] + backward(u_old2, i);
		end
		u1[begin] = u1[end];
		u2[begin] = u2[end];
		u = (u1 + u2) ./ 2;

		plt = plot(x, u, xlims=(0, lx), xticks=0:lx/5:lx, ylims=(-0.2, 1.2), yticks=-0.2:0.2:1.2,
		xlabel="x", ylabel="u", color="blue", title="QUICKEST with euler", label=@sprintf("t=%.1f", t));
		frame(anim, plt);
		t += dt;
		@printf("t = %.4f\n", t);
	end
	gif(anim, "quickest.gif", fps=5);
end

main();