text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Linear algebra module. Provides array arithmetic,
matrix factorizations and other linear algebra related
functionality.
"""
module LinearAlgebra

import Base: \, /, *, ^, +, -, ==
import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, asec, asech,
    asin, asinh, atan, atanh, axes, big, broadcast, ceil, conj, convert, copy, copyto!, cos,
    cosh, cot, coth, csc, csch, eltype, exp, findmax, findmin, fill!, floor, getindex, hcat,
    getproperty, imag, inv, isapprox, isone, IndexStyle, kron, length, log, map, ndims,
    oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round, sec, sech,
    setindex!, show, similar, sin, sincos, sinh, size, size_to_strides, sqrt, StridedReinterpretArray,
    StridedReshapedArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec
using Base: hvcat_fill, iszero, IndexLinear, promote_op, promote_typeof,
    @propagate_inbounds, @pure, reduce, typed_vcat, has_offset_axes
using Base.Broadcast: Broadcasted

export
# Modules
    LAPACK,
    BLAS,

# Types
    Adjoint,
    Transpose,
    SymTridiagonal,
    Tridiagonal,
    Bidiagonal,
    Factorization,
    BunchKaufman,
    Cholesky,
    CholeskyPivoted,
    Eigen,
    GeneralizedEigen,
    GeneralizedSVD,
    GeneralizedSchur,
    Hessenberg,
    LU,
    LDLt,
    QR,
    QRPivoted,
    LQ,
    Schur,
    SVD,
    Hermitian,
    Symmetric,
    LowerTriangular,
    UpperTriangular,
    UnitLowerTriangular,
    UnitUpperTriangular,
    Diagonal,
    UniformScaling,

# Functions
    axpy!,
    axpby!,
    bunchkaufman,
    bunchkaufman!,
    cholesky,
    cholesky!,
    cond,
    condskeel,
    copyto!,
    copy_transpose!,
    cross,
    adjoint,
    adjoint!,
    det,
    diag,
    diagind,
    diagm,
    dot,
    eigen,
    eigen!,
    eigmax,
    eigmin,
    eigvals,
    eigvals!,
    eigvecs,
    factorize,
    givens,
    hessenberg,
    hessenberg!,
    isdiag,
    ishermitian,
    isposdef,
    isposdef!,
    issuccess,
    issymmetric,
    istril,
    istriu,
    kron,
    ldiv!,
    ldlt!,
    ldlt,
    logabsdet,
    logdet,
    lowrankdowndate,
    lowrankdowndate!,
    lowrankupdate,
    lowrankupdate!,
    lu,
    lu!,
    lyap,
    mul!,
    lmul!,
    rmul!,
    norm,
    normalize,
    normalize!,
    nullspace,
    ordschur!,
    ordschur,
    pinv,
    qr,
    qr!,
    lq,
    lq!,
    opnorm,
    rank,
    rdiv!,
    schur,
    schur!,
    svd,
    svd!,
    svdvals!,
    svdvals,
    sylvester,
    tr,
    transpose,
    transpose!,
    transpose_type,
    tril,
    triu,
    tril!,
    triu!,

# Operators
    \,
    /,

# Constants
    I

const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const BlasReal = Union{Float64,Float32}
const BlasComplex = Union{ComplexF64,ComplexF32}

if USE_BLAS64
    const BlasInt = Int64
else
    const BlasInt = Int32
end

# Check that stride of matrix/vector is 1
# Writing like this to avoid splatting penalty when called with multiple arguments,
# see PR 16416
"""
    stride1(A) -> Int

Return the distance between successive array elements
in dimension 1 in units of element size.

# Examples
```jldoctest
julia> A = [1,2,3,4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> LinearAlgebra.stride1(A)
1

julia> B = view(A, 2:2:4)
2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64:
 2
 4

julia> LinearAlgebra.stride1(B)
2
```
"""
stride1(x) = stride(x,1)
stride1(x::Array) = 1
stride1(x::DenseArray) = stride(x, 1)::Int

@inline chkstride1(A...) = _chkstride1(true, A...)
@noinline _chkstride1(ok::Bool) = ok || error("matrix does not have contiguous columns")
@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride1(A) == 1), B...)

"""
    LinearAlgebra.checksquare(A)

Check that a matrix is square, then return its common dimension.
For multiple arguments, return a vector.

# Examples
```jldoctest
julia> A = fill(1, (4,4)); B = fill(1, (5,5));

julia> LinearAlgebra.checksquare(A, B)
2-element Array{Int64,1}:
 4
 5
```
"""
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

function checksquare(A...)
    sizes = Int[]
    for a in A
        size(a,1)==size(a,2) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(a))"))
        push!(sizes, size(a,1))
    end
    return sizes
end

function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

"""
    ldiv!(Y, A, B) -> Y

Compute `A \\ B` in-place and store the result in `Y`, returning the result.

The argument `A` should *not* be a matrix.  Rather, instead of matrices it should be a
factorization object (e.g. produced by [`factorize`](@ref) or [`cholesky`](@ref)).
The reason for this is that factorization itself is both expensive and typically allocates memory
(although it can also be done in-place via, e.g., [`lu!`](@ref)),
and performance-critical situations requiring `ldiv!` usually also require fine-grained
control over the factorization of `A`.

# Examples
```jldoctest
julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];

julia> X = [1; 2.5; 3];

julia> Y = zero(X);

julia> ldiv!(Y, qr(A), X);

julia> Y
3-element Array{Float64,1}:
  0.7128099173553719
 -0.051652892561983674
  0.10020661157024757

julia> A\\X
3-element Array{Float64,1}:
  0.7128099173553719
 -0.05165289256198333
  0.10020661157024785
```
"""
ldiv!(Y, A, B)

"""
    ldiv!(A, B)

Compute `A \\ B` in-place and overwriting `B` to store the result.

The argument `A` should *not* be a matrix.  Rather, instead of matrices it should be a
factorization object (e.g. produced by [`factorize`](@ref) or [`cholesky`](@ref)).
The reason for this is that factorization itself is both expensive and typically allocates memory
(although it can also be done in-place via, e.g., [`lu!`](@ref)),
and performance-critical situations requiring `ldiv!` usually also require fine-grained
control over the factorization of `A`.

# Examples
```jldoctest
julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];

julia> X = [1; 2.5; 3];

julia> Y = copy(X);

julia> ldiv!(qr(A), X);

julia> X
3-element Array{Float64,1}:
  0.7128099173553719
 -0.051652892561983674
  0.10020661157024757

julia> A\\Y
3-element Array{Float64,1}:
  0.7128099173553719
 -0.05165289256198333
  0.10020661157024785
```
"""
ldiv!(A, B)

"""
    rdiv!(A, B)

Compute `A / B` in-place and overwriting `A` to store the result.

The argument `B` should *not* be a matrix.  Rather, instead of matrices it should be a
factorization object (e.g. produced by [`factorize`](@ref) or [`cholesky`](@ref)).
The reason for this is that factorization itself is both expensive and typically allocates memory
(although it can also be done in-place via, e.g., [`lu!`](@ref)),
and performance-critical situations requiring `rdiv!` usually also require fine-grained
control over the factorization of `B`.
"""
rdiv!(A, B)

copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

include("adjtrans.jl")
include("transpose.jl")

include("exceptions.jl")
include("generic.jl")

include("blas.jl")
include("matmul.jl")
include("lapack.jl")

include("dense.jl")
include("tridiag.jl")
include("triangular.jl")

include("factorization.jl")
include("qr.jl")
include("hessenberg.jl")
include("lq.jl")
include("eigen.jl")
include("svd.jl")
include("symmetric.jl")
include("cholesky.jl")
include("lu.jl")
include("bunchkaufman.jl")
include("diagonal.jl")
include("bidiag.jl")
include("uniformscaling.jl")
include("givens.jl")
include("special.jl")
include("bitarray.jl")
include("ldlt.jl")
include("schur.jl")
include("structuredbroadcast.jl")
include("deprecated.jl")

const ⋅ = dot
const × = cross
export ⋅, ×

function versioninfo(io::IO=stdout)
    if Base.libblas_name == "libopenblas" || BLAS.vendor() == :openblas || BLAS.vendor() == :openblas64
        openblas_config = BLAS.openblas_get_config()
        println(io, "BLAS: libopenblas (", openblas_config, ")")
    else
        println(io, "BLAS: ",Base.libblas_name)
    end
    println(io, "LAPACK: ",Base.liblapack_name)
end

function __init__()
    try
        BLAS.check()
        if BLAS.vendor() == :mkl
            ccall((:MKL_Set_Interface_Layer, Base.libblas_name), Cvoid, (Cint,), USE_BLAS64 ? 1 : 0)
        end
        Threads.resize_nthreads!(Abuf)
        Threads.resize_nthreads!(Bbuf)
        Threads.resize_nthreads!(Cbuf)
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module LinearAlgebra")
    end
end

end # module LinearAlgebra

max_stars_repo_name: amilsted/julia

max_issues_repo_name: amilsted/julia

max_forks_repo_name: amilsted/julia


text:
url = "https://www-nds.iaea.org/phsp/software/iaea_phsp_Sept2013.zip"
dir = joinpath(@__DIR__, "src")
rm(dir, force=true, recursive=true)
mkdir(dir)
path = joinpath(dir, "blob")
download(url, path)
run(Cmd(`unzip $path`, dir=dir))
rm(path)
#run(Cmd(`make all --just-print`, dir=dir))

cmds = [
`g++ -O2 -fPIC -c -o iaea_header.o iaea_header.cpp`,
`g++ -O2 -fPIC -c -o iaea_phsp.o iaea_phsp.cpp`,
`g++ -O2 -fPIC -c -o iaea_record.o iaea_record.cpp`,
`g++ -O2 -fPIC -c -o utilities.o utilities.cpp`,
`g++ -O2 -fPIC -c -o iaea_event_generator.o iaea_event_generator.cpp`,
`g++ -O2 -fPIC -shared -o libiaea_phsp.so iaea_header.o iaea_phsp.o iaea_record.o utilities.o iaea_event_generator.o -ldl`,
]
for cmd in cmds
    run(Cmd(cmd, dir=dir))
end

ext: jl


text:
@assert typeof(3) == Int64
@assert typeof(0x22) == UInt8
@assert typeof(0xFA31) == UInt16
@assert typeof(false) == Bool
@assert typeof('\u263a') == Char
@assert typeof("Hello") == String
@assert typeof(BigInt(72)^3897) == BigInt
@assert typeof(:hello) == Symbol
@assert typeof(r"\d+(\.\d+)") == Regex
@assert typeof(Float64) == DataType
@assert typeof(DataType) == DataType

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = [-5 0.1875 0.0625]

dense_stencil = @def_stencil_expression(
            @sum(i,-2,2,
                @sum(j,-2,2,
                    @sum(k,-2,2, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))

st_inst = NewStencilInstance(dense_stencil, 2, coefs)

radius = 2
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

function cpu_dense_stencil(data, radius, coefs)
    dim = size(data)
    out = zeros(Float32, dim)
    for x = (radius+1):(dim[1]-radius)
        for y = (radius+1):(dim[2]-radius)
            for z = (radius+1):(dim[3]-radius)
                res = coefs[1]*data[x,y,z]
                for ix = -radius:radius
                    for iy = -radius:radius
                        for iz = -radius:radius
                            if ix == 0 &&  iy == 0 &&  iz == 0
                                continue
                            end
                            ci = max(abs.([ix, iy, iz])...)
                            ci = coefs[ci+1]
                            res += ci*data[x+ix,y+iy,z+iz]
                        end
                    end
                end
                out[x,y,z] = res
            end
        end
    end
    return out
end

cout = cpu_dense_stencil(data, radius, coefs)
for i = 2:t_steps
    global cout
    cout = cpu_dense_stencil(cout, radius, coefs)
end

@assert(sum(abs.(cout .- o))./(reduce(*, size(data))) < 0.001)
mean_err = sum(abs.(cout .- o))./(reduce(*, size(data)))

#
# using Plots
# shift = -13
# heatmap(cout[:,:,Int(round(dz/2))+shift])
# heatmap(o[:,:,Int(round(dz/2))+shift])
#
# shift = -14
# heatmap(o[:,:,Int(round(dz/2))+shift] .- cout[:,:,Int(round(dz/2))+shift])
#
# diff =  abs.(cout .- o)#./(reduce(*, size(data)))
# findmax(diff)
# sum(diff)./(reduce(*, size(data)))


text:
using SententialDecisionDiagrams
const SDD = SententialDecisionDiagrams

# set up vtree and manager
vtree = SDD.read_vtree("$(@__DIR__)/input/opt-swap.vtree")
manager = SDD.sdd_manager(vtree)

println("reading sdd from file ...")
α = SDD.read_sdd("$(@__DIR__)/input/opt-swap.sdd", manager)
println("sdd size = $(SDD.size(α))")

# ref, perform the minimization, and then de-ref
SDD.ref(α)
println("minimizing sdd size ... ")
SDD.minimize(manager)  # see also SDD.minimize(m,limited=true)
println("done!")
println("sdd size = $(SDD.size(α))")
SDD.deref(α)

# augment the SDD
println("augmenting sdd ...")
β = α & (SDD.literal(4,manager) | SDD.literal(5,manager))
println("sdd size = $(SDD.size(β))")

# ref, perform the minimization again on new SDD, and then de-ref
SDD.ref(β)
println("minimizing sdd ... ")
SDD.minimize(manager)
println("done!")
println("sdd size = $(SDD.size(β))")
SDD.deref(β)


text:
function generate_file_template(parent_folder, file_type, file_name="", file_params=Dict())

  file_suffix = "jl"
  if contains(file_type, ".")
    file_type, file_suffix = map(s -> String(s), split(file_type, "."; limit=2))
  end
  if file_name == "" ; file_name = file_type ; end

  template_path = "$(dirname(@__FILE__))/../../../templates/$parent_folder/$file_type.$file_suffix"

  template = readstring(template_path)

  template_dictionary = Dict()

  template_dictionary["app"] = replace(rsplit(pwd(), "/"; limit=2)[2], ".jl", "")
  template_dictionary["user"] = LibGit2.getconfig("github.user", "")
  template_dictionary["name"] = file_name
  template_dictionary["title"] = title(file_name)
  template_dictionary["class"] = get_class_name(file_name)
  template_dictionary["params"] = join(file_params, ", ")
  template_dictionary["fields"] = join(file_params, "\n  ")

  render(template, template_dictionary)

end


text:
using Gtk.ShortNames, Gtk.GConstants

using Printf
using GR

x = randn(1000000)
y = randn(1000000)

function plot(ctx, w, h)
    global sl

    ENV["GKS_WSTYPE"] = "142"
    ENV["GKSconid"] = @sprintf("%lu", UInt64(ctx.ptr))

    plt = gcf()
    plt[:size] = (w, h)
    nbins = Int64(Gtk.GAccessor.value(sl))

    hexbin(x, y, nbins=nbins)
end

function draw(widget)
    ctx = Gtk.getgc(widget)
    w = Gtk.width(widget)
    h = Gtk.height(widget)

    Gtk.rectangle(ctx, 0, 0, w, h)
    Gtk.set_source_rgb(ctx, 1, 1, 1)
    Gtk.fill(ctx)

    plot(ctx, w, h)
end

function resize_event(widget)
    ctx = Gtk.getgc(widget)
    h = Gtk.height(widget)
    w = Gtk.width(widget)

    Gtk.paint(ctx)
end

function motion_notify_event(widget::Gtk.GtkCanvas, event::Gtk.GdkEventMotion)
    Gtk.GAccessor.text(lb, @sprintf("(%g, %g)", event.x, event.y))
end

function value_changed(widget::Gtk.GtkScale)
    global canvas
    draw(canvas)
end

win = Window("Gtk") |> (bx = Box(:v))
Gtk.set_gtk_property!(win, :double_buffered, false)
lb = Label("(-, -)")
sl = Scale(false, 10, 100, 1)
Gtk.GAccessor.value(sl, 30)
canvas = Canvas(600, 450)
push!(bx, lb, sl, canvas)

signal_connect(motion_notify_event, canvas, "motion-notify-event")
signal_connect(value_changed, sl, "value_changed")

canvas.resize = resize_event
canvas.draw = draw

Gtk.showall(win)

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
end


text:
abstract type AbstractJammer{T} <: AbstractEmitter{T} end

struct CWJammer{
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    C <: Union{Unitful.Gain, AbstractAmplitude},
} <: AbstractJammer{T}
    id::Int
    doppler::typeof(1.0Hz)
    phase::Float64
    jnr::C
    exists::E
    doa::D
    signal::StructArray{Complex{T},1,NamedTuple{(:re, :im),Tuple{Array{T,1},Array{T,1}}},Int}
end

struct NoiseJammer{
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    C <: Union{Unitful.Gain, AbstractAmplitude},
} <: AbstractJammer{T}
    id::Int
    jnr::C
    exists::E
    doa::D
    signal::StructArray{Complex{T},1,NamedTuple{(:re, :im),Tuple{Array{T,1},Array{T,1}}},Int}
end

function CWJammer(
    ::Type{T},
    id::Integer,
    jnr::J;
    doppler = 0.0Hz,
    phase = 0.0,
    exists::E = true,
    doa::D = SVector(0, 0, 1)
) where {
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    J <: Union{Unitful.Gain{Unitful.Decibel, :?, <:Real}, AbstractAmplitude}
}
    phase = phase / 2π
    signal = StructArray{Complex{T}}(undef, 0)
    CWJammer{T, D, E, J}(
        id,
        float(doppler),
        float(phase),
        jnr,
        exists,
        doa,
        signal
    )
end

function CWJammer(
    id::Integer,
    jnr::Union{Unitful.Gain{Unitful.Decibel, :?, <:Real}, AbstractAmplitude};
    doppler = 0.0Hz,
    phase = 0.0,
    exists = true,
    doa = SVector(0, 0, 1)
)
    CWJammer(
        Float32,
        id,
        jnr,
        doppler = doppler,
        phase = phase,
        exists = exists,
        doa = doa
    )
end

function gen_signal!(
    jammer::CWJammer{T},
    sampling_frequency,
    intermediate_frequency,
    n0,
    num_samples::Integer,
    rng
) where T
    resize!(jammer.signal, num_samples)
    gen_carrier!(
        jammer.signal,
        intermediate_frequency + get_carrier_doppler(jammer),
        sampling_frequency,
        get_carrier_phase_2pi(jammer),
        T(get_amplitude(jammer.jnr, n0, sampling_frequency))
    )
end

function NoiseJammer(
    ::Type{T},
    id,
    jnr::J;
    exists::E = true,
    doa::D = SVector(0, 0, 1)
) where {
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    J <: Union{Unitful.Gain{Unitful.Decibel, :?, <:Real}, AbstractAmplitude}
}
    signal = StructArray{Complex{T}}(undef, 0)
    NoiseJammer{T, D, E, J}(
        id,
        jnr,
        exists,
        doa,
        signal
    )
end

function NoiseJammer(
    id,
    jnr;
    exists = true,
    doa = SVector(0, 0, 1)
)
    NoiseJammer(
        Float32,
        id,
        jnr,
        exists = exists,
        doa = doa
    )
end

function gen_signal!(
    jammer::NoiseJammer{T},
    sample_frequency,
    intermediate_frequency,
    n0,
    num_samples::Integer,
    rng
) where T
    resize!(jammer.signal, num_samples)
    signal = randn!(rng, jammer.signal)
    jammer.signal .*= T(get_amplitude(jammer.jnr, n0, sample_frequency))
    jammer.signal
end

function propagate(
    jammer::CWJammer,
    num_samples,
    intermediate_frequency,
    sample_frequency,
    n0,
    rng
)
    carrier_phase = (intermediate_frequency + get_carrier_doppler(jammer)) * num_samples /
        sample_frequency + get_carrier_phase_2pi(jammer)
    modded_carrier_phase = mod(carrier_phase + 0.5, 1) - 0.5
    Δt = num_samples / sample_frequency
    exists = propagate(jammer.exists, Δt, rng)
    doa = propagate(jammer.doa, Δt, rng)
    jnr = propagate(jammer.jnr, n0, sample_frequency, Δt, rng)
    CWJammer(
        jammer.id,
        get_carrier_doppler(jammer),
        modded_carrier_phase,
        jnr,
        exists,
        doa,
        jammer.signal
    )
end

function propagate(
    jammer::NoiseJammer,
    num_samples,
    intermediate_frequency,
    sample_frequency,
    n0,
    rng
)
    Δt = num_samples / sample_frequency
    exists = propagate(jammer.exists, Δt, rng)
    doa = propagate(jammer.doa, Δt, rng)
    jnr = propagate(jammer.jnr, n0, sample_frequency, Δt, rng)
    NoiseJammer(jammer.id, jnr, exists, doa, jammer.signal)
end

@inline get_carrier_phase(jammer::CWJammer) = 2π * jammer.phase
@inline get_carrier_phase_2pi(jammer::CWJammer) = jammer.phase
@inline get_carrier_doppler(jammer::CWJammer) = jammer.doppler
@inline get_id(jammer::AbstractJammer) = jammer.id
@inline get_jammer_to_noise_ratio(jammer::AbstractJammer) = get_jammer_to_noise_ratio(jammer.jnr)


text:
using SnoopCompile

file = joinpath(dirname(@__FILE__),"..","examples","streaming.jl")
@eval SnoopCompile.@snoop "/tmp/psychotask_compiles.csv" begin
  ARGS = ["test"]
  include($file)
end


text:
"""
`f` from page `b.jl`.

Links:

- [`ccall`](@ref)
- [`while`](@ref)
- [`@time`](@ref)
- [`T`](@ref)
- [`f`](@ref)
- [`f(::Any)`](@ref)
- [`f(::Any, ::Any)`](@ref)
- [`f(::Any, ::Any, ::Any)`](@ref)

"""
f(x, y) = x + y


text:
# ---
# title: 409. Longest Palindrome
# id: problem409
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/longest-palindrome/description/>
# hidden: true
# ---
#
# Given a string `s` which consists of lowercase or uppercase letters, return
# _the length of the **longest palindrome**_  that can be built with those
# letters.
#
# Letters are **case sensitive** , for example, `"Aa"` is not considered a
# palindrome here.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abccccdd"
#     Output: 7
#     Explanation:
#     One longest palindrome that can be built is "dccaccd", whose length is 7.
#
#
# **Example 2:**
#
#
#
#     Input: s = "a"
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: s = "bb"
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 2000`
#   * `s` consits of lower-case **and/or** upper-case English letters only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/env bash
#=
exec julia --color=yes --startup-file=no "${BASH_SOURCE[0]}" "$@"
=#

@show ARGS  # put any Julia code here
println("arg1:")
println(ARGS[1])

#=
$> sh /Users/merlinr/repo/julia-repos/julia-examples/julia-script-arg-handling.jl hello world
ARGS = ["hello", "world"]
=#

max_stars_repo_name: mahiki/julia-is-great

max_issues_repo_name: mahiki/julia-is-great

max_forks_repo_name: mahiki/julia-is-great


text:
function maxdiff(x::AbstractArray, y::AbstractArray)
    return mapreduce((a, b) -> abs(a - b), max, x, y)
end


text:
export Object3D, Scene, Visual, Solid
export Box, Beam, Cylinder, Sphere, Ellipsoid
export Cone, Capsule, GearWheel, Grid, SpringShape
export CoordinateSystem, FileMesh
export Font, TextShape, ModelicaShape
export VisualMaterial
export MassProperties, MassPropertiesFromShape
export MassPropertiesFromShapeAndMass
export UniformGravityField
export RefPath, ptpJointSpace, scheduleReferenceMotion
export calculateRobotMovement
export getRefPathPosition, getRefPathInitPosition, getVariables
export multibodyResiduals!, setModiaJointVariables!
export Fix
export Revolute, RevoluteWithFlange
export Prismatic, PrismaticWithFlange
export J123, J132, J123or132, singularRem, FreeMotion, change_rotSequenceInNextIteration!

export buildModia3D

import Modia3D
using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")
include("$(ModiaLang.path)/models/Translational.jl")

include("$(Modia3D.path)/src/ModiaInterface/buildModia3D.jl")

Object3D(        ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Object3D), _path = true, kwargs...)
Scene(           ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Scene)                 , kwargs...)
Visual(          ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Visual)                , kwargs...)
Solid(           ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Solid)                 , kwargs...)
Box(             ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Box)                   , kwargs...)
Beam(            ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Beam)                  , kwargs...)
Cylinder(        ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Cylinder)              , kwargs...)
Sphere(          ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Sphere)                , kwargs...)
Ellipsoid(       ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Ellipsoid)             , kwargs...)
Capsule(         ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Capsule)               , kwargs...)
Cone(            ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Cone)                  , kwargs...)
SpringShape(     ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Spring)                , kwargs...)
GearWheel(       ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.GearWheel)             , kwargs...)
Grid(            ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Grid)                  , kwargs...)
VisualMaterial(  ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.VisualMaterial)        , kwargs...)
MassProperties(  ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.MassProperties)        , kwargs...)
CoordinateSystem(; kwargs...) = Par(; _constructor = :(Modia.Modia3D.CoordinateSystem)      , kwargs...)
FileMesh(        ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.FileMesh)              , kwargs...)
Font(            ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Font)                  , kwargs...)
TextShape(       ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.TextShape)             , kwargs...)
ModelicaShape(   ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.ModelicaShape)         , kwargs...)
Fix(             ; kwargs...) = Par(; _constructor = :(Modia.Modia3D.Fix)                   , kwargs...)

MassPropertiesFromShape()  = Par(; _constructor = :(Modia.Modia3D.MassPropertiesFromShape))
MassPropertiesFromShapeAndMass(;mass) = Par(; _constructor = :(Modia.Modia3D.MassPropertiesFromShapeAndMass), mass = mass)
UniformGravityField(; kwargs...) = Par(; _constructor = :(Modia.Modia3D.UniformGravityField), kwargs...)

RefPath(; kwargs...) = Modia3D.ReferencePath(; kwargs...)

ptpJointSpace(; kwargs...) = Modia3D.ptpJointSpace(; kwargs...)
scheduleReferenceMotion(; kwargs...) = Modia3D.scheduleReferenceMotion(; kwargs...)

calculateRobotMovement(args...) = Modia3D.calculateRobotMovement(args...)
getRefPathPosition(args...) = Modia3D.getRefPathPosition(args...)
getRefPathInitPosition(args...) = Modia3D.getRefPathInitPosition(args...)

getVariables(args...) = (args...,)

multibodyResiduals!(args...)     = Modia3D.multibodyResiduals!(args...)
setModiaJointVariables!(args...) = Modia3D.setModiaJointVariables!(args...)

Revolute(; obj1, obj2, axis=3, phi=Var(init=0.0), w=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia.Modia3D.Revolute), _path = true, ndof = 1),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    axis = Par(value = axis),
    canCollide = Par(value = canCollide),
    phi  = phi,
    w    = w,
    equations = :[
        w   = der(phi)
        qdd = der(w)   # standardized name for the generalized joint accelerations
        variables = getVariables(phi, w, 0.0) # standardized name for the generalized joint position, velocity, force
        ]
)

RevoluteWithFlange(; obj1, obj2, axis=3, phi=Var(init=0.0), w=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia.Modia3D.Revolute), _path = true, ndof = 1),
    obj1   = Par(value = obj1),
    obj2   = Par(value = obj2),
    axis   = Par(value = axis),
    canCollide = Par(value = canCollide),
    flange = Flange,
    phi    = phi,
    w      = w,
    equations = :[
        phi = flange.phi
        w   = der(phi)
        qdd = der(w)   # standardized name for the generalized joint accelerations
        variables = getVariables(phi, w, flange.tau) # standardized name for the generalized joint position, velocity, force
        ]
)

Prismatic(; obj1, obj2, axis=1, s=Var(init=0.0), v=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia.Modia3D.Prismatic), _path = true, ndof = 1),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    axis = Par(value = axis),
    canCollide = Par(value = canCollide),
    s    = s,
    v    = v,
    equations = :[
        v   = der(s)
        qdd = der(v)   # standardized name for the generalized joint accelerations
        variables = getVariables(s, v, 0.0) # standardized name for the generalized joint position, velocity, force
        ]
)

PrismaticWithFlange(; obj1, obj2, axis=1, s=Var(init=0.0), v=Var(init=0.0), canCollide=true) = Model(; _constructor = Par(value = :(Modia.Modia3D.Prismatic), _path = true, ndof = 1),
    obj1   = Par(value = obj1),
    obj2   = Par(value = obj2),
    axis   = Par(value = axis),
    canCollide = Par(value = canCollide),
    flange = TranslationalFlange,
    s      = s,
    v      = v,
    equations = :[
        s   = flange.s
        v   = der(s)
        qdd = der(v)   # standardized name for the generalized joint accelerations
        variables = getVariables(s, v, flange.f) # standardized name for the generalized joint position, velocity, force
        ]
)

"""
    J = J123(rot123::AbstractVector)

Return joint rot. kinematics matrix `J` for Cardan angles `rot123` (rotation sequence x-y-z).
"""
function J123(rot123::AbstractVector)

    (sbe, cbe) = sincos(rot123[2])
    (sga, cga) = sincos(rot123[3])
    return [     cga     -sga  0.0 ;
             cbe*sga  cbe*cga  0.0 ;
            -sbe*cga  sbe*sga  cbe ] / cbe

end

"""
    J = J132(rot132::AbstractVector)

Return joint rot. kinematics matrix `J` for Cardan angles `rot132` (rotation sequence x-z-y).
"""
function J132(rot132::AbstractVector)

    (sga, cga) = sincos(rot132[2])
    (sbe, cbe) = sincos(rot132[3])
    return [ cbe      0.0  sbe     ;
            -sbe*cga  0.0  cbe*cga ;
             cbe*sga  cga  sbe*sga ] / cga

end

"""
    next_isrot123 = change_rotSequenceInNextIteration!(rot::AbstractVector, isrot123::Bool, instantiatedModel::SimulationModel, x, rot_name)

Change rotation sequence of `rot` from `x-axis, y-axis, z-axis` to `x-axis, z-axis, y-axis` or visa versa in the next event iteration:

- If `isrot123 = true`, return `next_isrot123 = false` and `x[..] = rot132fromR(Rfromrot123(rot))`

- If `isrot123 = false`, return `next_isrot123 = true` and `x[..] = rot123fromR(Rfromrot132(rot))`
"""
function change_rotSequenceInNextIteration!(rot::AbstractVector, isrot123::Bool, instantiatedModel::SimulationModel, x, rot_name)::Bool
    if isrot123
        #println("        switch $rot_name 123 -> 132")
        next_rot      = Modia3D.rot132fromR(Modia3D.Rfromrot123(rot))
        next_isrot123 = false
    else
        #println("        switch $rot_name 132 -> 123")
        next_rot      = Modia3D.rot123fromR(Modia3D.Rfromrot132(rot))
        next_isrot123 = true
    end

    # Change x-vector with the next_rot values
    eqInfo = instantiatedModel.equationInfo
    startIndex = eqInfo.x_info[ eqInfo.x_dict[rot_name] ].startIndex
    x[startIndex]   = next_rot[1]
    x[startIndex+1] = next_rot[2]
    x[startIndex+2] = next_rot[3]
    return next_isrot123
end

singularRem(ang) = abs(rem2pi(ang, RoundNearest)) - 1.5  # is negative/positive in valid/singular angle range
J123or132(rot, isrot123) = isrot123 ? J123(rot) : J132(rot)

FreeMotion(; obj1, obj2, r=Var(init=zeros(3)), rot=Var(init=zeros(3)), v=Var(init=zeros(3)), w=Var(init=zeros(3))) = Model(; _constructor = Par(value = :(Modia.Modia3D.FreeMotion), _path = true, ndof = 6),
    obj1 = Par(value = obj1),
    obj2 = Par(value = obj2),
    r    = r,
    rot  = rot,
    v    = v,
    w    = w,

    next_isrot123 = Var(start=true),
    _rotName = "???",  # is changed by buildModia3D to the full path name of "rot"

    equations = :[
        der(r) = v

        isrot123 = pre(next_isrot123)
        rot2_singularity = positive(singularRem(rot[2]))
        next_isrot123 = if rot2_singularity; change_rotSequenceInNextIteration!(rot, isrot123, instantiatedModel, _x, _rotName) else isrot123 end
        der(rot) = J123or132(rot,isrot123) * w

        der(v) = qdd[1:3]
        der(w) = qdd[4:6]
        variables = getVariables(r, rot, v, w, isrot123)
        ]
)


text:
module GeoCoderUs
  include("geocommon.jl")
  using HTTPClient.HTTPC
  base_url = "http://rpc.geocoder.us/service"

  type GeoCoderClient
    username::String
    password::String
  end

  function geocode(location::String)
    resp = HTTPC.get("$base_url/csv",
                    RequestOptions(;query_params=[("address",location)]))
    return resp.body
  end

  function simple_geocode(location::String)
    csv = geocode(location)

    parts = split(strip(csv),",",3)
    lat = float(parts[1])
    lng = float(parts[2])
    place = parts[3]
    return GeoCommon.SimpleLocation(place,lat,lng)
  end

end

max_stars_repo_name: astrieanna/geojl

max_issues_repo_name: astrieanna/geojl

max_forks_repo_name: astrieanna/geojl


text:
"""
    v_func()

Lorem ipsum dolor sit amet.
"""
function v_func()
  9
end

"""
    v_func(cur_obj::GType)

Lorem ipsum dolor sit amet.
"""
function v_func(cur_obj::GType)
  8
end

"""
    v_func(cur_obj::AType)

Lorem ipsum dolor sit amet.
"""
function v_func(cur_obj::AType)
  7
end


text:
module Visualize

using GeometryBasics
using FileIO
using Plots

export Mesh3D, create_ply_file, show_closed_cross_section, show_multi_branch_cross_section

function Mesh3D(CrossSectionNodes, z, u, v, ϕ)

    #define planar cross-section geometry
    x=view(CrossSectionNodes,:,1)
    y=view(CrossSectionNodes,:,2)

    # write deformed cross-section in 3D coordinates, section by section along the z axis
    for i in 1:length(z)
        CrossSectionShape=hcat(x.+u[i]+x.*cos.(ϕ[i]).-y.*sin.(ϕ[i]),y.+v[i]+x.*sin.(ϕ[i]).+y.*cos.(ϕ[i]),ones(length(x))*z[i])
        if i==1
            global coordinates=CrossSectionShape   #start with first cross-section
        else
            coordinates=vcat(coordinates, CrossSectionShape)  #add following cross-sections along z
        end
    end

    # apply triangular mesh to 3D coordinates
    NumCrossSectionNodes=length(x)
    NumZNodes=length(z)

    #uppoly are triangles pointing towards +z, downpoly are triangles pointing towards -z
    uppoly=[[i i+1 i+1+NumCrossSectionNodes].+NumCrossSectionNodes*(j-1) for i=1:NumCrossSectionNodes-1 for j=1:NumZNodes-1]
    downpoly=[[i i+NumCrossSectionNodes i+NumCrossSectionNodes+1].+NumCrossSectionNodes*(j-1) for i=1:NumCrossSectionNodes-1 for j=1:NumZNodes-1]

    connectivity=vcat(uppoly,downpoly)  #combine all the triangles
    connectivity=permutedims(reshape(hcat(connectivity...), (length(connectivity[1]), length(connectivity))))  #convert into format to be read by Makie

    return coordinates, connectivity   #these are ready for Makie scenes, e.g., poly()

end

"""
    create_ply_file(coordinates, connectivity, filename)

Accepts triangular surface mesh `coordinates` and `connectivity` along with a `filename' and saves a .ply file.

https://en.wikipedia.org/wiki/PLY_(file_format)

The coordinates and connectivity are converted to GeometryBasics.jl primitives before saving with FileIO.jl.

"""

function create_ply_file(coordinates, connectivity, filename)

    points = [Point{3, Float64}(coordinates[i,1], coordinates[i,2], coordinates[i,3]) for i = 1:size(coordinates)[1]]

    facets = [TriangleFace{Cint}(connectivity[i, 1:3]) for i in eachindex(connectivity[:,1])]

    mesh = GeometryBasics.Mesh(points, facets)

    # using FileIO, MeshIO
    FileIO.save(filename, mesh)

end

function show_closed_cross_section(xcoords, ycoords, markershape, xlims, ylims)

    num_nodes = length(xcoords)

    plt = []

    for i = 1:num_nodes

        if i == 1

            plt = plot([xcoords[i], xcoords[i+1]], [ycoords[i], ycoords[i+1]], size = (600, 600), legend = false, linecolor = :black, markercolor = :black, markershape = markershape, seriestype = :line, xlims = xlims, ylims = ylims)

        elseif i<num_nodes

            plot!(plt, [xcoords[i], xcoords[i+1]], [ycoords[i], ycoords[i+1]], size = (600, 600), legend = false, linecolor = :black, markercolor = :black, markershape = markershape, seriestype = :line, xlims = xlims, ylims = ylims)

        elseif i==num_nodes

            plot!(plt, [xcoords[i], xcoords[1]], [ycoords[i], ycoords[1]], size = (600, 600), legend = false, linecolor = :black, markercolor = :black, markershape = markershape, seriestype = :line, xlims = xlims, ylims = ylims)

        end

    end

    return plt

end

function show_multi_branch_cross_section(xcoords, ycoords, element_connectivity, markershape, xlims, ylims)

    num_elem = size(element_connectivity)[1]

    plt = []

    for i = 1:num_elem

        node_i = Int(element_connectivity[i, 1])
        node_j = Int(element_connectivity[i, 2])

        if i == 1

            plt = plot([xcoords[node_i], xcoords[node_j]], [ycoords[node_i], ycoords[node_j]], size = (600, 600), legend = false, linecolor = :black, markercolor = :black, markershape = markershape, seriestype = :line, xlims = xlims, ylims = ylims)

        else

            plot!(plt, [xcoords[node_i], xcoords[node_j]], [ycoords[node_i], ycoords[node_j]], size = (600, 600), legend = false, linecolor = :black, markercolor = :black, markershape = markershape, seriestype = :line, xlims = xlims, ylims = ylims)

        end

    end

    return plt

end

end #module


text:
# Based on
# https://gitlab.com/agentydragon/solid-flask/-/blob/master/solid_flask_main.py

# This implementation is not secure.
# Many required cryptographic checks were ignored.

using HTTP
using JSON
using Base64
using Random
using SHA
using Sockets
using PyCall
using Dates

const pyjwk = pyimport("jwcrypto.jwk")
const pyjwt = pyimport("jwcrypto.jwt")

const _ISSUER        = "https://inrupt.net/"
const _PORT          = 3333
const _REDIRECT_PATH = "/oauth/callback"
const _REDIRECT_URI  = "http://localhost:" * string(_PORT) * _REDIRECT_PATH
const _AUTHORITY     = _ISSUER * ".well-known/openid-configuration"

const random_device = RandomDevice()
const state_storage = Dict{String,Dict{String,String}}()
const access_token  = Ref("")
const key           = Ref("")

function register_client(provider_info)
    url     = provider_info["registration_endpoint"]
    headers = ["content-type" => "application/json"]
    data    = Dict("redirect_uris" => [_REDIRECT_URI])
    result  = JSON.parse(String(HTTP.post(url, headers, JSON.json(data)).body))
    @info "registration response" result
    return result["client_id"]
end

function make_random_string()
    randstring(random_device, 50)
end

function make_verifier_challenge()
    code_verifier  = make_random_string()
    code_challenge = base64encode(sha256(code_verifier))
    code_challenge = replace(replace(replace(code_challenge, "+" => "-"), "/" => "_"), "=" => "")
    code_verifier, code_challenge
end

const provider_info = JSON.parse(String(HTTP.request("GET", _AUTHORITY).body))
@info "provider info" provider_info

const client_id = register_client(provider_info)

function template(web_id, resource, resource_content)
    web_str = web_id != "" ?
        "You are logged in as $(web_id)." :
        "You are not logged in."
    res_str = resource != "" ?
        "
$(resource_content)
" :
        "Use the form below to read a resource."

    """

Login status

    $(web_str)

Resource content

    $(res_str)

    <form action=/ method=GET>
      <input
          value="$(resource)"
          placeholder='https://you.solidcommunity.net/private/...'
          name='resource'>
      <input type=submit value=Read>
    </form>
    """
end

function make_token_for(pykeypair, uri, method)
    header = Dict(
        "typ" => "dpop+jwt",
        "alg" => "ES256",
        "jwk" => pykeypair.export(private_key=false, as_dict=true)
    )
    claims = Dict(
        "jti" => make_random_string(),
        "htm" => method,
        "htu" => uri,
        "iat" => Int(floor(datetime2unix(now())))
    )
    jwt = pyjwt.JWT(header=header, claims=claims)
    jwt.make_signed_token(pykeypair)
    return jwt.serialize()
end

function index(req::HTTP.Request)
    queryparams = HTTP.queryparams(HTTP.URI(req.target))
    tested_url  = get(queryparams, "resource", "")

    if access_token[] != "" && key[] != ""
        println("loading access token and key from session")
        keypair = pyjwk.JWK.from_json(key[])
        headers = [
            "authorization" => "DPoP " * access_token[],
            "dpop"          => make_token_for(keypair, tested_url, "GET")
        ]
        decoded_access_token = pyjwt.JWT(jwt=access_token[])
        web_id = JSON.parse(decoded_access_token.token.objects["payload"])["sub"]
    else
        headers = Pair{String, String}[]
        web_id = ""
    end

    if tested_url != ""
        resp = HTTP.get(tested_url, headers; status_exception = false)
        if resp.status == 401
            println("Got 401 trying to access " * tested_url)
            code_verifier, code_challenge = make_verifier_challenge()
            state = make_random_string()
            @assert !haskey(state_storage,state)
            state_storage[state] = Dict("code_verifier" => code_verifier, "redirect_url" => req.target)

            url = provider_info["authorization_endpoint"] * "?" * HTTP.escapeuri([
                "code_challenge"        => code_challenge,
                "state"                 => state,
                "response_type"         => "code",
                "redirect_uri"          => _REDIRECT_URI,
                "code_challenge_method" => "S256",
                "client_id"             => client_id,
                "scope"                 => "openid offline_access"
            ])
            return HTTP.Response(302,["Location" => url])
        elseif resp.status != 200
            error("Unexpected status code: $(resp.status) $(resp.body)")
        end
        resource_content = String(resp.body)
    else
        resource_content = ""
    end
    headers = ["Content-Type" => "text/html"]
    return HTTP.Response(200, headers; body=template(web_id,tested_url,resource_content))
end

function oauth_callback(req::HTTP.Request)
    queryparams = HTTP.queryparams(HTTP.URI(req.target))
    auth_code   = get(queryparams, "code", "")
    state       = get(queryparams, "state", "")

    @assert haskey(state_storage, state) "state $(state) not in state_storage?"

    pykeypair     = pyjwk.JWK.generate(kty="EC", crv="P-256")
    key[]         = pykeypair.export()

    state_info    = state_storage[state]
    code_verifier = state_info["code_verifier"]
    redirect_url  = state_info["redirect_url"]

    url     = provider_info["token_endpoint"]
    headers = [
        "dpop" => make_token_for(pykeypair, url, "POST"),
        "content-type" => "application/x-www-form-urlencoded"
    ]
    data    = HTTP.escapeuri([
        "grant_type"    => "authorization_code",
        "client_id"     => client_id,
        "redirect_uri"  => _REDIRECT_URI,
        "code"          => auth_code,
        "code_verifier" => code_verifier
    ])
    result = JSON.parse(String(HTTP.post(url, headers, data; redirect=false).body))
    access_token[] = result["access_token"]

    pop!(state_storage, state)
    return HTTP.Response(302, ["Location" => redirect_url])
end

const ROUTER = HTTP.Router()
HTTP.@register(ROUTER, "GET", "/", index)
HTTP.@register(ROUTER, "GET", "/oauth/callback", oauth_callback)
@async HTTP.serve(ROUTER, ip"127.0.0.1", _PORT)


text:
# ------------------------------------------------------------------------------------------
# ## Juliabox is based on JuliaRun
#
# The distributed compute facilities in JuliaBox are provided by
# [JuliaRun](https://juliacomputing.com/products/juliarun.html), Julia Computing's solution
# for scaling Julia processes. There are available in JuliaBox via the `JuliaRunClient`
# package, which allows you to run interactive or batch distributed processes.
#
# For the free edition of JuliaBox, users are limited to 2 CPU cores. For running larger
# clusters, please contact `juliabox@juliacomputing.com`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Understanding the cluster
#
# `JuliaRunClient` provides a `Context` object which is the entry point into the cluster.
# All operations need a reference to this object
# ------------------------------------------------------------------------------------------

using JuliaRunClient
ctx = Context()

# ------------------------------------------------------------------------------------------
# Let's see if we can connect to the cluster
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Running batch job
#
# JuliaRun is well suited to run batch processes in a distributed environment. We have below
# an example of running the Monte Carlo Pi computation as batch job.
#
# First, we create an instance of a `JuliaParBatch` object, with name _mcpi_
# ------------------------------------------------------------------------------------------

job = JuliaParBatch("mcpi")

# ------------------------------------------------------------------------------------------
# We then _submit_ the job to the *JuliaRun* system, providing it with the scripts to run on
# the master and worker nodes, as well as the attached filesystem volume names. Update
# `num_workers`, `master_cpu`, `master_mem`,`worker_cpu` and `worker_memory`  as
# appropriate.
# ------------------------------------------------------------------------------------------

num_workers=2
master_cpu=1
master_mem="2Gi"
worker_cpu=1
worker_memory="2Gi"

@result submitJob(ctx, job; start_script="/juliabox/scripts/mcpi_pmap_master.jl",
        run_volume="juliabox", image="juliabox", nworkers=num_workers,
        cpu=master_cpu, memory=master_mem,
        shell="/juliabox/scripts/master.sh", worker_shell="/juliabox/scripts/worker.sh",
        worker_cpu=worker_cpu, worker_memory=worker_memory,
        worker_start_script="/juliabox/scripts/mcpi_pmap_worker.jl")

# Wait for a while
sleep(25)

# ------------------------------------------------------------------------------------------
# While the job is running, we can query it's status.
# ------------------------------------------------------------------------------------------

@result getJobStatus(ctx, job)

# ------------------------------------------------------------------------------------------
# The logs for the job is visible on the shared filesystem
# ------------------------------------------------------------------------------------------

# Wait for a while to see the output
sleep(20)

;tail /mnt/juliabox/logs/output1

# ------------------------------------------------------------------------------------------
# Once the job is completed, we can clean it up and free it's resources. Again, this is
# important to do for cost reasons.
# ------------------------------------------------------------------------------------------

@result deleteJob(ctx, job; force=true)

;rm -f /mnt/juliabox/logs/output1

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# diagonal operators: joDiag

export joDiag
"""
    julia> op = joDiag(v;[makecopy=...,][DDT=...,][RDT=...,][name=...])

Diagonal matrix with elements from given vector

# Signature

    function joDiag(v::LocalVector{EDT};makecopy::Bool=true,
        DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),name::String="joDiag") where {EDT}

# Arguments

- `v`: vector of diagonal elements
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    joDiag(rand(m))

examples with DDT/RDT

    joDiag(rand(m); DDT=Float32)
    joDiag(rand(m); DDT=Float32,RDT=Float64)

"""
function joDiag(v::LocalVector{EDT};makecopy::Bool=true,
    DDT::DataType=EDT,RDT::DataType=promote_type(EDT,DDT),name::String="joDiag") where {EDT}

    vc= makecopy ? Base.deepcopy(v) : v

    return joMatrix(Diagonal(vc);DDT=DDT,RDT=RDT,name=name)
end


text:
export SolitaryWaveWhitham

"""
    `SolitaryWaveWhitham(param; kwargs...)`

Computes the Whitham solitary wave with prescribed velocity.

# Argument
- `param :: NamedTuple`: parameters of the problem containing velocity `c` and dimensionless parameters `ϵ` and `μ`, and mesh size `L` and number of collocation points `N`;

## Keywords (optional)
- `guess :: Vector{Real}`: initial guess for the surface deformation (if not provided, the exact formula for KdV is used);
- `x₀ :: Real`: center of solitary wave (if guess is not provided);
- `iterative :: Bool`: inverts Jacobian through GMRES if `true`, LU decomposition if `false`;
- `verbose :: Bool`: prints numerical errors at each step if `true`;
- `max_iter :: Int`: maximum number of iterations of the Newton algorithm;
- `tol :: Real`: general tolerance (default is `1e-10`);
- `ktol :: Real`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `gtol :: Real`: relative tolerance of the GMRES algorithm;
- `dealias :: Int`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `q :: Real`: Newton algorithm modified with
`u_{n+1}=q*u_{n+1}+(1-q)*u_n`
(default is `1`);
- `α :: Real`: adds `α` times spectral projection onto the Kernel to the Jacobian;
- `KdV :: Bool`: if `true` computes the KdV (instead of Whitham) solitary wave.
# Return values
`u :: Vector{Float64}` the solution
"""
function SolitaryWaveWhitham(
                param :: NamedTuple;
                guess = zeros(0) :: Vector{Float64},
                x₀ = 0 :: Real,
                iterative = false :: Bool,
                verbose = true :: Bool,
                max_iter = 20 :: Int,
                tol = 1e-14 :: Real,
                gtol = 1e-14 :: Real,
                ktol = 0 :: Real,
                dealias = 0 :: Int,
                q=1 :: Real,
                α=0 :: Real,
                KdV = false :: Bool)

        mesh=Mesh(param)
        c = param.c
        ϵ = param.ϵ
        μ = param.μ

        if guess == zeros(0)
                @info "Using the exact formula for the KdV solitary wave as initial guess"
                guess = 2*(c-1)/ϵ*sech.(sqrt(3*2*(c-1)/μ)/2*mesh.x.-x₀).^2
        end

        k  = mesh.k
        Dx =  1im * k
        if KdV == true
                @info string("Solving the KdV solitary wave with velocity c=",param.c)
                F₁ = 1 .-μ/6*k.^2
        else
                @info string("Solving the Whitham solitary wave with velocity c=",param.c)
                F₁ = sqrt.(tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k)))
                F₁[1]=1
        end

        if dealias == 0
                Π = ones(size(k))
        else
                Π = abs.(k) .< maximum(k) * (1-dealias/(dealias+2))
        end
        krasny(k) = (abs.(k).> ktol ).*k
        krasny!(k) = k[abs.(k).< ktol ].=0

        function filter( v :: Vector{Float64} )
                real.(ifft(krasny(Π.*fft(v))))
        end
        function filter( v :: Vector{Complex{Float64}} )
                ifft(krasny(Π.*fft(v)))
        end

        function F( u :: Vector{Float64} )
            -c*u+real.(ifft(F₁.*fft(u)))+3*ϵ/4*u.^2
        end
        function Fabs( u :: Vector{Float64} )
            c*abs.(u)+abs.(ifft(F₁.*fft(u)))+3*ϵ/4*u.^2
        end

        if iterative == false
                k = mesh.k
                x = mesh.x
                x₀ = mesh.x[1]
                FFT = exp.(-1im*k*(x.-x₀)');
                IFFT = exp.(1im*k*(x.-x₀)')/length(x);
                M = real.(IFFT*(Diagonal( F₁)*FFT))
                function JacF( v , dxv  )
                    Symmetric(M+Diagonal(3*ϵ/2*v .-c)+ α*dxv*dxv')
                end
        else
                function JacFfast( v ,dxv )
                        dF(φ) = -c*φ+real.(ifft(F₁.*fft(φ)))+3*ϵ/2*v.*φ+ α*dot(dxv,φ)*dxv
                        return LinearMap(dF, length(v); issymmetric=true, ismutating=false)
                end
        end

        flag=0
        iter = 0
        u = filter(guess)
        du = similar(u)
        fu = similar(u)
        dxu = similar(u)
        for i in range(0, stop=max_iter)
                dxu .= real.(ifft(Dx.*fft(u)))
                dxu ./= norm(dxu,2)
                fu .= F(u)
                relerr = norm(fu,Inf)/norm(Fabs(u),Inf)
                abserr = norm(fu,Inf)
                if relerr < tol
    			@info string("Converged : relative error ",relerr," in ",i," steps\n")
    			break
    		elseif verbose == true
                        print(string("absolute error at step ",i,": ",abserr,"\n"))
                        print(string("relative error at step ",i,": ",relerr,"\n"))

    		end
                if i == max_iter
                        @warn string("The algorithm did not converge after ",i," steps: final relative error is ",relerr,"\n")
                        break
                end

                if iterative == false
                        du .= JacF(u,dxu) \ fu
                else
                        du .= gmres( JacFfast(u,dxu) , fu ; verbose = verbose, tol = gtol )
                end
    		u .-= q*filter(du)
        end
        return u
end


text:
using PMDPs

PMDPs.run_small_instances()


text:
# PLOT example
using PyPlot
using PyCall
@pyimport numpy as np

x = np.random[:randn](200)
xm = mean(x)

plot(x)
title("α = $xm")
xlabel("time")

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
"""
Shock and Detonation Toolbox
"Postshock" module

Calculates CJ detonation speed and post-shock states for frozen and equilibrium cases.

This module defines the following functions:

    LSQ_CJspeed
    hug_fr
    hug_eq
    FHFP
    CJ_calc
    CJspeed
    PostShock_fr
    PostShock_eq
    shock_calc
    shk_eq_calc

################################################################################
Theory, numerical methods and applications are described in the following report:

    SDToolbox Numerical Tools for Shock and Detonation Wave Modeling,
    Explosion Dynamics Laboratory, Contributors:
    S. Browne, J. Ziegler, N. Bitter, B. Schmidt, J. Lawson and J. E. Shepherd,
    GALCIT Technical Report FM2018.001 Revised January 2021.
    California Institute of Technology, Pasadena, CA USA

Please cite this report and the website if you use these routines.

Please refer to LICENCE.txt or the above report for copyright and disclaimers.

http://shepherd.caltech.edu/EDL/PublicResources/sdt/

"""
module Postshock

export CJspeed, PostShock_eq, PostShock_fr

using PyCall
const ct = PyNULL()

function __init__()
    copy!(ct, pyimport_conda("cantera","cantera","cantera"))
    global R̄ = ct.gas_constant
end

# constants from sdtoolbox/config
const ERRFT = 1e-4
const ERRFV = 1e-4
const volumeBoundRatio = 5

"""
    LSQ_CJspeed(x,y)

Determines least squares fit of parabola to input data

FUNCTION SYNTAX:
[a,b,c,R2,SSE,SST] = LSQ_CJspeed(x,y)

INPUT:
    x = independent data points
    y = dependent data points

OUTPUT:
    a,b,c = coefficients of quadratic function (ax^2 + bx + c = 0)
    R2 = R-squared value
    SSE = sum of squares due to error
    SST = total sum of squares

"""
function LSQ_CJspeed(x,y)

    # Calculate Sums
    X = 0.0; X2 = 0.0; X3 = 0.0; X4 = 0.0;
    Y = 0.0; Y1 = 0.0; Y2 = 0.0;
    a = 0.0; b = 0.0; c = 0.0; R2 = 0.0
    n = length(x)

    for k = 1:n
        X += x[k]
        X2 += x[k]^2
        X3 += x[k]^3
        X4 += x[k]^4
        Y += y[k]
        Y1 += y[k]*x[k]
        Y2 += y[k]*x[k]^2
    end
    m = Y/n

    den = X3*n - X2*X
    temp = den*(X*X2-X3*n) + X2*X2*(X*X-n*X2) - X4*n*(X*X-X2*n)
    temp2 = den*(Y*X2-Y2*n) + (Y1*n-Y*X)*(X4*n-X2*X2)

    b = temp2/temp
    a = 1.0/den*(n*Y1 - Y*X - b*(X2*n-X*X))
    c = 1/n*(Y - a*X2 - b*X)

    SSE = 0.0; SST = 0.0;

    for k = 1:n
        f = a*x[k]^2 + b*x[k] + c
        SSE += (y[k] - f)^2
        SST += (y[k] - m)^2
    end
    R2 = 1 - SSE/SST

    return a,b,c,R2,SSE,SST
end

"""
    hug_fr(x,vb,h1,P1,v1,gas)

Computes difference in enthalpy computed from assumed (T, V)
state and fixed composition (frozen) hugoniot evaluation.  Used
with root solver such as 'fsolve' to compute frozen hugoniot as a function
of volume.  Input gas object is modified to correspond to input (T, V).

FUNCTION SYNTAX:
    diff = hug_fr(x,vb,h1,P1,v1,gas)

USAGE:
    fval = fsolve(hug_fr,Ta,args=(vb,h1,P1,v1,gas))
        = frozen Hugoniot temperature (K) corresponding to vb

INPUT:
    Ta = initial guess for frozen Hugoniot temperature (K)
    vb = desired frozen Hugoniot specific volume (m^3/kg)
    h1 = enthalpy at state 1 (J/kg)
    P1 = pressure at state 1 (Pa)
    v1 = specific volume at state 1 (m^3/kg)
    gas = working gas object

OUTPUT:
    enthalpy difference

"""
function hug_fr(x,vb,h1,P1,v1,gas)

    gas.TD = x, 1.0/vb
    hb1 = gas.enthalpy_mass
    Pb = R̄ * x/(gas.mean_molecular_weight*vb)

    hb2 = h1 + 0.5*(Pb-P1)*(vb+v1)
    return hb2-hb1
end

"""
    hug_eq(x,vb,h1,P1,v1,gas)

Computes difference in enthalpy computed from assumed (T, V)
state and equilibrium composition state hugoniot evaluation.  Used
with root solver such as 'fsolve' to compute equilibrium hugoniot as a function
of volume.  Input gas object is modified to correspond to input (T, V) and
an equilibrium composition.

FUNCTION SYNTAX:
    diff = hug_eq(x,vb,h1,P1,v1,gas)

USAGE:
    fval = fsolve(hug_eq,Ta,args=(vb,h1,P1,v1,gas))
        = equilibrium Hugoniot temperature (K) corresponding to vb

INPUT:
    Ta = initial guess for equilibrium Hugoniot temperature (K)
    vb = desired equilibrium Hugoniot specific volume (m^3/kg)
    h1 = enthalpy at state 1 (J/kg)
    P1 = pressure at state 1 (Pa)
    v1 = specific volume at state 1 (m^3/kg)
    gas = working gas object

OUTPUT:
    enthalpy difference

"""
function hug_eq(x,vb,h1,P1,v1,gas)

    gas.TD = x, 1.0/vb
    gas.equilibrate("TV")
    hb1 = gas.enthalpy_mass
    Pb = R̄ * x / (gas.mean_molecular_weight*vb)
    hb2 = h1 + 0.5*(Pb-P1)*(vb+v1)
    return hb2-hb1
end

"""
    FHFP(w1,gas2,gas1)

Uses the momentum and energy conservation equations to calculate
error in pressure and enthalpy given shock speed, upstream (gas1)
and downstream states (gas2).  States are not modified by these routines.

FUNCTION SYNTAX:
    [FH,FP] = FHFP(w1,gas2,gas1)

INPUT:
    w1 = shock speed (m/s)
    gas2 = gas object at working/downstream state
    gas1 = gas object at initial/upstream state

OUTPUT:
    FH,FP = error in enthalpy and pressure

"""
function FHFP(w1,gas2,gas1)

    P1 = gas1.P
    H1 = gas1.enthalpy_mass
    r1 = gas1.density
    P2 = gas2.P
    H2 = gas2.enthalpy_mass
    r2 = gas2.density
    w1s = w1^2
    w2s = w1s*(r1/r2)^2
    FH = H2 + 0.5*w2s - (H1 + 0.5*w1s)
    FP = P2 + r2*w2s - (P1 + r1*w1s)
    return FH, FP
end

"""
    CJ_calc(gas, gas1, ERRFT, ERRFV, x)

Calculates the Chapman-Jouguet wave speed using Reynolds' iterative method.

FUNCTION SYNTAX:
    [gas,w1] = CJ_calc(gas,gas1,ERRFT,ERRFV,x)

INPUT:
    gas = working gas object
    gas1 = gas object at initial state
    ERRFT,ERRFV = error tolerances for iteration
    x = density ratio

OUTPUT:
    gas = gas object at equilibrium state
    w1 = initial velocity to yield prescribed density ratio

"""
function CJ_calc(gas, gas1, ERRFT, ERRFV, x)

    T = 2000; r1 = gas1.density
    V1 = 1/r1
    i = 0; DT = 1000; DW = 1000;
    #PRELIMINARY GUESS
    V = V1/x; r = 1/V; w1 = 2000
    gas.TD = T, r
    gas.equilibrate("TV")
    H, P = gas.HP
    # START LOOP
    while (abs(DT) > ERRFT*T) || (abs(DW) > ERRFV*w1)
        i += 1
        if i == 500
            println("CJ_calc did not converge")
            return
        end
        # CALCULATE FH & FP FOR GUESS 1
        FH,FP = FHFP(w1,gas,gas1)
        # TEMPERATURE PERTURBATION
        DT = T*0.02; Tper = T + DT
        Vper = V; Rper = 1/Vper
        Wper = w1
        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(Wper,gas,gas1)
        #ELEMENTS OF JACOBIAN
        DFHDT = (FHX-FH)/DT; DFPDT = (FPX-FP)/DT
        # VELOCITY PERTURBATION
        DW = 0.02*w1; Wper = w1 + DW
        Tper = T; Rper = 1/V
        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(Wper,gas,gas1)
        #ELEMENTS OF JACOBIAN
        DFHDW = (FHX-FH)/DW; DFPDW = (FPX-FP)/DW
        # INVERT MATRIX
        J = DFHDT*DFPDW - DFPDT*DFHDW
        b = DFPDW, -DFHDW, -DFPDT, DFHDT
        a = -FH, -FP
        DT = (b[1]*a[1]+b[2]*a[2])/J; DW = (b[3]*a[1]+b[4]*a[2])/J
        # CHECK & LIMIT CHANGE VALUES
        # VOLUME
        DTM = 0.2*T
        if abs(DT) > DTM
            DT = DTM*DT/abs(DT)
        end
        # MAKE THE CHANGES
        T += DT; w1 += DW
        gas.TD = T, r
        gas.equilibrate("TV")
    end
    return w1
end

""" defaults to CEA Algorithm """
function CJspeed(P1, T1, q, mech; method = "CEA", kwargs...)
    if method == "CEA"
        _CJspeed_CEA(P1, T1, q, mech)
    elseif method == "umin"
        _CJspeed_umin(P1, T1, q, mech; kwargs...)
    elseif method == "aeq"
        _CJspeed_aeq(P1, T1, q, mech)
    end
end

"""
    _CJspeed_umin(P1, T1, q, mech; fullOutput=false)

Calculates CJ detonation velocity for a given pressure, temperature, and
composition.

FUNCTION SYNTAX:
    If only CJ speed required:
    cj_speed = CJspeed(P1,T1,q,mech)
    If full output required:
    [cj_speed,R2,plot_data] = CJspeed(P1,T1,q,mech,fullOutput=true)

INPUT:
    P1 = initial pressure (Pa)
    T1 = initial temperature (K)
    q = reactant species mole fractions in one of Cantera's recognized formats
    mech = cti file containing mechanism data (e.g. 'gri30.cti')

OPTIONAL INPUT:
    fullOutput = set true for R-squared value and pre-formatted plot data
                (the latter for use with sdtoolbox.utilities.CJspeed_plot)

OUTPUT
    cj_speed = CJ detonation speed (m/s)
    R2 = R-squared value of LSQ curve fit (optional)
    plot_data = tuple (rr,w1,dnew,a,b,c)
                rr = density ratio
                w1 = speed
                dnew = minimum density
                a,b,c = quadratic fit coefficients

"""
function _CJspeed_umin(P1, T1, q, mech; fullOutput=false)

    #DECLARATIONS
    numsteps = 20; maxv = 2.0; minv = 1.5
    w1 = zeros(Float64, numsteps+1)
    rr = zeros(Float64,numsteps+1)
    gas1 = ct.Solution(mech)
    gas  = ct.Solution(mech)
    #INTIAL CONDITIONS
    gas.TPX  = T1, P1, q
    gas1.TPX = T1, P1, q
    #INITIALIZE ERROR VALUES & CHANGE VALUES
    ERRFT = 1e-4;  ERRFV = 1e-4

    T1 = gas1.T; P1 = gas1.P
    counter = 1; R2 = 0.0; cj_speed = 0.0
    a = 0.0; b = 0.0; c = 0.0; dnew = 0.0
    while (counter <= 4) || (R2 < 0.99999)
        step = (maxv-minv)/numsteps; i = 1; x = minv
        while x <= maxv
            gas.TPX = T1, P1, q
            w1[i] = CJ_calc(gas, gas1, ERRFT, ERRFV, x)
            rr[i] = gas.density/gas1.density
            i += 1; x += step
        end
        a,b,c,R2,SSE,SST = LSQ_CJspeed(rr,w1)
        dnew = -b/(2.0*a)
        minv = dnew * (1 - 0.001)
        maxv = dnew * (1 + 0.001)
        counter += 1
        cj_speed = a*dnew^2 + b*dnew + c
    end

    if fullOutput
        # Optional output data for plotting (with sdtoolbox.utilities.CJspeed_plot)
        plot_data = rr,w1,dnew,a,b,c
        return cj_speed,R2,plot_data
    else
        return cj_speed
    end
end

"""
    PostShock_fr(U1, P1, T1, q, mech)

Calculates frozen post-shock state for a specified shock velocity and pre-shock state.

FUNCTION SYNTAX:
    gas = PostShock_fr(U1,P1,T1,q,mech)

INPUT:
    U1 = shock speed (m/s)
    P1 = initial pressure (Pa)
    T1 = initial temperature (K)
    q = reactant species mole fractions in one of Cantera's recognized formats
    mech = cti file containing mechanism data (e.g. 'gri30.cti')

OUTPUT:
    gas = gas object at frozen post-shock state

"""
function PostShock_fr(U1, P1, T1, q, mech)

    gas1 = ct.Solution(mech)
    gas  = ct.Solution(mech)
    # INTIAL CONDITIONS
    gas.TPX = T1, P1, q
    gas1.TPX = T1, P1, q
    # CALCULATES POST-SHOCK STATE
    gas = shk_calc(U1, gas, gas1, ERRFT, ERRFV)
end

"""
    PostShock_eq(U1, P1, T1, q, mech)

Calculates equilibrium post-shock state for a specified shock velocity and pre-shock state.

FUNCTION SYNTAX:
    gas = PostShock_eq(U1,P1,T1,q,mech)

INPUT:
    U1 = shock speed (m/s)
    P1 = initial pressure (Pa)
    T1 = initial temperature (K)
    q = reactant species mole fractions in one of Cantera's recognized formats
    mech = cti file containing mechanism data (e.g. 'gri30.cti')

OUTPUT:
    gas = gas object at equilibrium post-shock state

"""
function PostShock_eq(U1, P1, T1, q, mech)

    gas1 = ct.Solution(mech);
    gas  = ct.Solution(mech);
    # INTIAL CONDITIONS
    # workaround to avoid unsized object error when only one species in a .cti file
    # (flagged to be fixed in future Cantera version)
    if length(q) > 1
        gas.TPX = T1, P1, q
        gas1.TPX= T1, P1, q
    else
        gas.TP = T1, P1
        gas1.TP= T1, P1
    end
    # CALCULATES POST-SHOCK STATE
    gas = shk_eq_calc(U1, gas, gas1, ERRFT, ERRFV)
end

"""
    shk_calc(U1, gas, gas1, ERRFT, ERRFV)

Calculates frozen post-shock state using Reynolds' iterative method.

FUNCTION SYNTAX:
    gas = shk_calc(U1,gas,gas1,ERRFT,ERRFV)

INPUT:
    U1 = shock speed (m/s)
    gas = working gas object
    gas1 = gas object at initial state
    ERRFT,ERRFV = error tolerances for iteration

OUTPUT:
    gas = gas object at frozen post-shock state

"""
function shk_calc(U1, gas, gas1, ERRFT, ERRFV)

    r1 = gas1.density; V1 = 1/r1
    P1 = gas1.P; T1 = gas1.T
    i = 0
    deltaT = 1000; deltaV = 1000
    # PRELIMINARY GUESS
    Vg = V1/volumeBoundRatio; rg = 1/Vg
    Pg = P1 + r1*(U1^2)*(1-Vg/V1); Tg = T1*Pg*Vg/(P1*V1)
    gas.TD = Tg, rg
    Hg, Pg = gas.HP
    # SAVE STATE
    V = Vg; r = rg; P = Pg
    T = Tg; H = Hg
    # START LOOP
    while (abs(deltaT) > ERRFT*T) || (abs(deltaV) > ERRFV*V)
        i += 1
        if i == 500
            println("shk_calc did not converge for U = $U1")
            return gas
        end
        # CALCULATE FH & FP FOR GUESS 1
        FH,FP = FHFP(U1,gas,gas1)

        # TEMPERATURE PERTURBATION
        DT = T*0.02; Tper = T + DT;
        Vper = V; Rper = 1/Vper;
        gas.TD = Tper, Rper
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(U1,gas,gas1)
        # ELEMENTS OF JACOBIAN
        DFHDT = (FHX-FH)/DT; DFPDT = (FPX-FP)/DT

        # VOLUME PERTURBATION
        DV = 0.02*V; Vper = V + DV
        Tper = T; Rper = 1/Vper
        gas.TD = Tper, Rper
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(U1,gas,gas1)
        # ELEMENTS OF JACOBIAN
        DFHDV = (FHX-FH)/DV; DFPDV = (FPX-FP)/DV

        # INVERT MATRIX
        J = DFHDT*DFPDV - DFPDT*DFHDV
        b = DFPDV, -DFHDV, -DFPDT, DFHDT
        a = -FH, -FP
        deltaT = (b[1]*a[1]+b[2]*a[2])/J; deltaV = (b[3]*a[1]+b[4]*a[2])/J

        # CHECK & LIMIT CHANGE VALUES
        # TEMPERATURE
        DTM = 0.2*T
        if abs(deltaT) > DTM
            deltaT = DTM*deltaT/abs(deltaT)
        end
        # VOLUME
        V2X = V + deltaV
        if V2X > V1
            DVM = 0.5*(V1 - V)
        else
            DVM = 0.2*V
        end
        if abs(deltaV) > DVM
            deltaV = DVM*deltaV/abs(deltaV)
        end
        # MAKE THE CHANGES
        T += deltaT; V += deltaV; r = 1/V
        gas.TD = T, r
    end
    return gas
end

"""
    shk_eq_calc(U1, gas, gas1, ERRFT, ERRFV)

Calculates equilibrium post-shock state using Reynolds' iterative method.

FUNCTION SYNTAX:
    gas = shk_calc(U1,gas,gas1,ERRFT,ERRFV)

INPUT:
    U1 = shock speed (m/s)
    gas = working gas object
    gas1 = gas object at initial state
    ERRFT,ERRFV = error tolerances for iteration

OUTPUT:
    gas = gas object at equilibrium post-shock state

"""
function shk_eq_calc(U1, gas, gas1, ERRFT, ERRFV)

    r1 = gas1.density; V1 = 1/r1
    P1 = gas1.P; T1 = gas1.T
    i = 0
    deltaT = 1000; deltaV = 1000
    # PRELIMINARY GUESS
    V = V1/volumeBoundRatio; r = 1/V
    P = P1 + r1*(U1^2)*(1-V/V1); T = T1*P*V/(P1*V1)
    gas.TD = T,r
    H,P = gas.HP
    # START LOOP
    while (abs(deltaT) > ERRFT*T) || (abs(deltaV) > ERRFV*V)
        i += 1
        if i == 500
            println("shk_calc did not converge for U = $U1")
            return gas
        end
        # CALCULATE FH & FP FOR GUESS 1
        FH,FP = FHFP(U1,gas,gas1)

        # TEMPERATURE PERTURBATION
        DT = T*0.02; Tper = T + DT
        Vper = V; Rper = 1/Vper
        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(U1,gas,gas1)
        # ELEMENTS OF JACOBIAN
        DFHDT = (FHX-FH)/DT; DFPDT = (FPX-FP)/DT

        # VOLUME PERTURBATION
        DV = 0.02*V; Vper = V + DV
        Tper = T; Rper = 1/Vper
        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        # CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX = FHFP(U1,gas,gas1)
        # ELEMENTS OF JACOBIAN
        DFHDV = (FHX-FH)/DV; DFPDV = (FPX-FP)/DV

        # INVERT MATRIX
        J = DFHDT*DFPDV - DFPDT*DFHDV
        b = DFPDV, -DFHDV, -DFPDT, DFHDT
        a = -FH, -FP
        deltaT = (b[1]*a[1]+b[2]*a[2])/J; deltaV = (b[3]*a[1]+b[4]*a[2])/J

        # CHECK & LIMIT CHANGE VALUES
        # TEMPERATURE
        DTM = 0.2*T
        if abs(deltaT) > DTM
            deltaT = DTM*deltaT/abs(deltaT)
        end
        # VOLUME
        V2X = V + deltaV
        if V2X > V1
            DVM = 0.5*(V1 - V)
        else
            DVM = 0.2*V
        end
        if abs(deltaV) > DVM
            deltaV = DVM*deltaV/abs(deltaV)
        end
        # MAKE THE CHANGES
        T += deltaT; V += deltaV; r = 1/V
        gas.TD = T,r
        gas.equilibrate("TV")
    end

    return gas
end

"""
    _CJspeed_aeq(P1, T1, q, mech)

CJspeed_aeq
Calculates CJ detonation velocity and CJ state based on equilibrium sound speed

FUNCTION
SYNTAX
cj_speed,gas = CJspeed(P1,T1,q,mech)

INPUT
P1 = initial pressure (Pa)
T1 = initial temperature (K)
q = string of reactant species mole fractions
mech = cti file containing mechanism data (i.e. 'gri30.cti')

OUTPUT
cj_speed = CJ detonation speed (m/s)

"""
function _CJspeed_aeq(P1, T1, q, mech)

    gas2 = ct.Solution(mech)
    gas1 = ct.Solution(mech)
    gas  = ct.Solution(mech)

    #INTIAL CONDITIONS
    gas.TPX = T1, P1, q
    gas1.TPX = T1, P1, q

    #INITIALIZE ERROR VALUES & CHANGE VALUES
    ERRFT = 1e-4;  ERRFV = 1e-4

    r1 = gas1.density; V1 = 1/r1
    P1 = gas1.P; T1 = gas1.T
    i = 0
    #PRELIMINARY GUESS
    Vg = V1/5; rg = 1/Vg

    gas.TD = T1,rg
    gas.equilibrate("UV")
    Tg = gas.T
    gas2.TDX = Tg, rg, gas.X

    #SAVE STATE
    V = Vg; r = rg
    T = Tg
    deltaT = 1000; deltaV = 1000; cj_speed = 0
    #START LOOP
    while (abs(deltaT) > ERRFT*T) || (abs(deltaV) > ERRFV*V)
        i += 1
        if i == 500
            println("CJspeed_sound did not converge")
            return
        end

        #CALCULATE FH & FP FOR GUESS 1
        FH,FP,cj_speed = FHFP_CJ2(gas,gas1,gas2)

        #TEMPERATURE PERTURBATION
        DT = T*0.01; Tper = T + DT
        Vper = V; Rper = 1/Vper

        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        gas2.TDX = Tper, Rper, gas.X

        #CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX,cj_speed = FHFP_CJ2(gas,gas1,gas2)
        #ELEMENTS OF JACOBIAN
        DFHDT = (FHX-FH)/DT; DFPDT = (FPX-FP)/DT

        #VOLUME PERTURBATION
        DV = 0.01*V; Vper = V + DV
        Tper = T; Rper = 1/Vper

        gas.TD = Tper, Rper
        gas.equilibrate("TV")
        gas2.TDX = Tper, Rper, gas.X

        #CALCULATE FHX & FPX FOR "IO" STATE
        FHX,FPX,cj_speed = FHFP_CJ2(gas,gas1,gas2)
        #ELEMENTS OF JACOBIAN
        DFHDV = (FHX-FH)/DV; DFPDV = (FPX-FP)/DV

        #INVERT MATRIX
        J = DFHDT*DFPDV - DFPDT*DFHDV
        b = DFPDV, -DFHDV, -DFPDT, DFHDT
        a = -FH, -FP
        deltaT = (b[1]*a[1]+b[2]*a[2])/J; deltaV = (b[3]*a[1]+b[4]*a[2])/J

        #CHECK & LIMIT CHANGE VALUES
        #TEMPERATURE
        DTM = 0.2*T
        if abs(deltaT) > DTM
            deltaT = DTM*deltaT/abs(deltaT)
        end
        #VOLUME
        V2X = V + deltaV
        if V2X > V1
            DVM = 0.5*(V1 - V)
        else
            DVM = 0.2*V
        end
        if abs(deltaV) > DVM
            deltaV = DVM*deltaV/abs(deltaV)
        end
        #MAKE THE CHANGES
        T += deltaT; V += deltaV; r = 1/V
        gas.TD = T, r
        gas.equilibrate("TV")
        gas2.TDX = T, r, gas.X
    end

    FH,FP,cj_speed = FHFP_CJ2(gas,gas1,gas2)

    return cj_speed
end

"""
    FHFP_CJ2(gas,gas1,gas2)

FHFP_CJ2
Uses the momentum and energy conservation equations and the equilibrium sound speed to calculate error in current pressure and enthalpy guesses.  In this case, state 2 is in equilibrium.

FUNCTION
SYNTAX
[FH,FP,cj_speed] = FHFP_CJ2(gas,gas1,gas2)

INPUT
gas = working gas object
gas1 = gas object at initial state
gas2 = dummy gas object (for calculating numerical derivatives)

OUTPUT
FH,FP = error in enthalpy and pressure
cj_speed = CJ detonation speed (m/s)

"""
function FHFP_CJ2(gas,gas1,gas2)

    P1 = gas1.P
    H1 = gas1.enthalpy_mass
    r1 = gas1.density
    P2 = gas.P
    H2 = gas.enthalpy_mass
    r2 = gas.density

    speeds = equilSoundSpeeds(gas2)
    w2s = speeds[1]^2
    w1s = w2s*(r2/r1)^2
    FH = H2 + 0.5*w2s - (H1 + 0.5*w1s)
    FP = P2 + r2*w2s - (P1 + r1*w1s)
    return FH, FP, sqrt(w1s)
end

"""
    equilSoundSpeeds(gas)

Calculates equilibrium and frozen sound speeds. For the equilibrium sound speed, the gas is equilibrated holding entropy and specific volume constant.

FUNCTION
SYNTAX
[aequil,afrozen] = equilSoundSpeeds(gas)

INPUT
gas = working gas object (modified inside function)

OUTPUT
aequil = equilibrium sound speed (m/s)
afrozen = frozen sound speed (m/s)

"""
function equilSoundSpeeds(gas)

    # set the gas to equilibrium at its current T and P
    gas.equilibrate("TP")

    # save properties
    s0 = gas.entropy_mass
    p0 = gas.P
    r0 = gas.density

    # perturb the density
    r1 = r0*1.0001

    # set the gas to a state with the same entropy and composition but
    # the perturbed density
    gas.SV = s0, 1.0/r1

    # save the pressure for this case for the frozen sound speed
    pfrozen = gas.P

    # now equilibrate the gas holding S and V constant
    gas.equilibrate("SV")

    p1 = gas.P

    # equilibrium sound speed
    aequil = sqrt((p1 - p0)/(r1 - r0))

    # frozen sound speed
    afrozen = sqrt((pfrozen - p0)/(r1 - r0));
    return aequil, afrozen
end

"""
    _CJspeed_CEA(P1,T1,X1,mech)

This function calculates the CJ-detonation velocity for given initial
conditions.

It is based on (i.e. mostly copied from) the CJ-detonation part of the
NASA CEA Fortran95 code, available at
http://www.grc.nasa.gov/WWW/CEAWeb/ceaHome.htm,
which itself is based on "Calculation of Detonation Properties and Effect
of Independent Parameters on Gaseous Detonations" by Frank J. Zeleznik
and Sanford Gordon, ARS Journal, April 1962.

INPUT:
P1 = pressure in (Pa)
T1 = temperature in (K)
X1 = mole fraction column vector
mech = file name of kinetic mechanism, e.g. 'gri30.xml'

OUTPUT:
cj_speed = CJ detonation velocity in (m/s)
"""
function _CJspeed_CEA(P1,T1,X1,mech)

    gas = ct.Solution(mech)
    gas.TPX = T1,P1,X1

    M1 = gas.mean_molecular_weight
    h1 = gas.enthalpy_mass

    P = 15 * P1

    # Acquire initial guess for T2 by finding flame temperature according to
    # chapter 8.3 of CEA documentation
    h_init = h1 + 0.75 * R̄ * T1/M1 * P/P1
    gas.HP = h_init, P
    gas.equilibrate("HP")
    T_init = gas.T
    T = T_init

    # Further improving initial guesses for T and P by recursion formulas
    for i = 1:3
        gas.TPX = T,P,X1
        gas.equilibrate("TP")
        M2 = gas.mean_molecular_weight
        α = T1/T * M2/M1
        γₛ, cₚ_eq, Dlvtp, Dlvpt = partDeriv(T,P,X1,gas)
        P = P1 * (1+γₛ) / (2*γₛ*α) * (1 + (1- (4*γₛ*α)/(1+γₛ)^2)^0.5)
        r = α * P/P1
        T = T1 * T_init/T1 - 0.75 * R̄/(M1*cₚ_eq) * 15 + (R̄*γₛ)/(2*M1*cₚ_eq) * (r^2-1)/r * P/P1
    end

    # Main iteration loop
    cj_speed = 0
    I = 0
    while true
        I += 1
        gas.TPX = T,P,X1
        gas.equilibrate("TP")
        M2 = gas.mean_molecular_weight
        γₛ, cₚ_eq, Dlvtp, Dlvpt = partDeriv(T,P,X1,gas)
        α = T1/T * M2/M1
        r = α * P/P1
        a₁₁ = P1/P + γₛ*r*Dlvpt
        a₁₂ = γₛ*r*Dlvtp
        a₂₁ = 0.5*γₛ*(r^2 - 1 - Dlvpt*(1+r^2)) + Dlvtp - 1
        a₂₂ = -0.5*γₛ*Dlvtp*(r^2+1) - M2*cₚ_eq/R̄
        b₁ = P1/P - 1 + γₛ*(r-1)
        b₂ = M2*(gas.enthalpy_mass-h1)/(R̄*T) - 0.5*γₛ*(r^2-1)
        d = a₁₁*a₂₂ - a₁₂*a₂₁
        x₁ = (a₂₂*b₁-a₁₂*b₂)/d
        x₂ = (a₁₁*b₂-a₂₁*b₁)/d

        δ = 1
        temp = maximum((abs(x₁), abs(x₂)))
        if temp > 0.4054652
            δ = 0.4054652/temp
        end
        P = P*exp(x₁*δ)
        T = T*exp(x₂*δ)
        soundspeed_equil = sqrt(R̄*γₛ*T/M2)
        cj_speed = r*soundspeed_equil

        # Convergence test
        if I > 8
            println("Could not converge in 8 iterations")
            break
        elseif I <= 8 && temp < 0.5e-4
            break
        end
    end
    return cj_speed
end

"""
    partDeriv(T,P,X,gas)

Calculates all needed variables related to partial derivatives, i.e γₛ (isentropic exponent),
cₚ_eq (equilibrium specific heat at constant pressure), Dlvtp (d ln(v) / d ln(T) at constant
pressure), Dlvpt (d ln(v) / d ln(P) at constant temperature). Derivatives are approximated
via central differences.
"""
function partDeriv(T,P,X,gas)

    v = 1/gas.density
    ΔT = T*1e-4
    ΔP = P*1e-4

    v_P = zeros(2)
    v_T = zeros(2)
    h = zeros(2)

    sign = -1, 1

    for i = 1:2

        gas.TPX = T, P + sign[i]*ΔP, X
        gas.equilibrate("TP")
        v_P[i] = 1/gas.density

        gas.TPX = T + sign[i]*ΔT, P, X
        gas.equilibrate("TP")
        v_T[i] = 1/gas.density
        h[i] = gas.enthalpy_mass
    end

    Dlvtp = T/v * (v_T[2] - v_T[1]) / (2*ΔT)
    Dlvpt = P/v * (v_P[2] - v_P[1]) / (2*ΔP)
    cₚ_eq = (h[2]-h[1])/(2*ΔT)

    gas.TPX = T,P,X
    gas.equilibrate("TP")

    γₛ = -1/(Dlvpt + R̄ * Dlvtp^2 / (cₚ_eq * gas.mean_molecular_weight))

    return γₛ, cₚ_eq, Dlvtp, Dlvpt
end

end #module


text:
using Test

include("../src/inmet.jl")

using Main.inmet:listarEstacoes, coletarDadosEstacao, coletarDadosEstacoesCSV
using DataFrames
import Dates

@testset verbose = true "Coleta de dados" begin
    @testset "Listagem de estações INMET" begin
        listaEstacoes = listarEstacoes(true)
        colunasListaEstacoes = names(listaEstacoes)
        @test typeof(listaEstacoes) == DataFrame
        @test nrow(listaEstacoes) > 0
        @test "CD_SITUACAO" in colunasListaEstacoes
        @test "CD_ESTACAO" in colunasListaEstacoes
    end

    @testset "Coleta de dados de uma estação" begin
        formato = Dates.DateFormat("y-m-d")
        dataInicial = Dates.Date("2021-01-01", formato)
        dataFinal = Dates.Date("2021-01-02", formato)
        serie = coletarDadosEstacao("A756", dataInicial, dataFinal)
        @test typeof(serie) == DataFrame
        @test nrow(serie) > 0
        @test "CHUVA" in names(serie)
    end

    @testset "Persistência de dados em CSV" begin
        formato = Dates.DateFormat("y-m-d")
        dataInicial = Dates.Date("2021-01-01", formato)
        dataFinal = Dates.Date("2021-01-02", formato)
        coletarDadosEstacoesCSV(dataInicial, dataFinal, "/tmp/inmet", ["A756", "A223"])
        @test Base.isfile("/tmp/inmet/A756.csv")
        @test Base.isfile("/tmp/inmet/A223.csv")
    end
end;

max_stars_repo_name: ericmiguel/inmet-julia

max_issues_repo_name: ericmiguel/inmet-julia

max_forks_repo_name: ericmiguel/inmet-julia


text:
using Distributions
using MLBase

"Generate a correlation matrix of a given size and density of network edges."
function randcor(K, netDensity)
    IC = diagm(abs(randn(K)))
    for i in 1:K, j in 1:i-1
        IC[i,j] = rand() < netDensity ? randn()-1 : 0.0
        IC[j,i] = IC[i,j]
    end
    mineval = minimum(eig(IC)[1])
    if mineval < 0.01
        IC -= eye(K)*mineval*1.1
    end
    C = inv(IC)

    Base.cov2cor!(C, sqrt(diag(C)))
end

function upper(X::AbstractMatrix)
    x = Float64[]
    for i in 1:size(X)[1], j in i+1:size(X)[2]
        push!(x, X[i,j])
    end
    x
end
function area_under_pr(truth::AbstractVector, predictor::AbstractVector; resolution=4000)
    rocData = MLBase.roc(round(Int64, truth), float(invperm(sortperm(predictor))), resolution)
    vals = collect(map(x->(recall(x), -precision(x)), rocData))
    sort!(vals)
    xvals = map(x->x[1], vals)
    yvals = map(x->-x[2], vals)
    # println(xvals)
    # println(yvals)
    area_under_curve([0.0; xvals], [yvals[1]; yvals]) # make sure we extend all the way to zero
end
function area_under_curve(x, y) # must be sorted by increasing x
    area = 0.0
    lastVal = NaN
    for i in 2:length(x)
        v = (y[i-1]+y[i])/2 * (x[i]-x[i-1])
        if !isnan(v)
            area += v
            lastVal = v
        elseif !isnan(lastVal)
            area += lastVal
        end
    end
    area
end
function area_under_pr(truth::AbstractMatrix, predictor::AbstractMatrix)
    area_under_pr(abs(upper(truth)) .> 0.01, abs(upper(predictor)))
end


text:
using Hexagons

immutable HexagonalBinGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    function HexagonalBinGeometry(
            default_statistic::Gadfly.StatisticElement=Gadfly.Stat.identity())
        new(default_statistic)
    end

    function HexagonalBinGeometry(; xbincount=50, ybincount=50)
        new(Gadfly.Stat.hexbin(xbincount=xbincount, ybincount=ybincount))
    end
end

const hexbin = HexagonalBinGeometry

function default_statistic(geom::HexagonalBinGeometry)
    geom.default_statistic
end

function element_aesthetics(geom::HexagonalBinGeometry)
    [:x, :y, :xsize, :ysize, :color]
end

function render(geom::HexagonalBinGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    default_aes.xsize = [1.0]
    default_aes.ysize = [1.0]
    aes = inherit(aes, default_aes)

    Gadfly.assert_aesthetics_defined("Geom.hexbin", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.hexbin", aes, :x, :y)

    n = length(aes.x)

    forms = Array(Compose.Form, 0)
    for (i, c, xs, ys) in zip(1:n, cycle(aes.color), cycle(aes.xsize), cycle(aes.ysize))
        if !isna(c)
            form = compose(polygon(hexpoints((aes.x[i], aes.y[i]), xs, ys)...),
                           linewidth(0.1mm), fill(c), stroke(c))
            push!(forms, form)
        end
    end

    compose(combine(forms...), svgclass("geometry"))
end


text:
"""
    RegularizedSVDPseudoInverse(A; ϵ=1e-8)

Structure able to solve linear system `A` using a regularized SVD.

The threshold `ϵ` controls the regularization.
Let ``σ_{\\max}`` be the larger singular value.
A singular value ``σ`` such that ``σ \\leq ϵ σ_{\\max}`` is approximated by
zero when solving the least-squares problem.
"""
struct RegularizedSVDPseudoInverse
    A
    ϵ
    F
    I
    Uϵ
    iΣϵ
    Vϵ
    function RegularizedSVDPseudoInverse(A; ϵ=1e-8)
        F = svd(A)
        I = argmin(F.S .> ϵ * F.S[1]) - 1
        if I == 0 I = length(F.S) end
        iΣϵ = spdiagm(0 => 1 ./ F.S[1:I])
        return new(A, ϵ, F, I, F.U[:,1:I], iΣϵ, F.V[:,1:I])
    end
end

"""
    condition_number(B::RegularizedSVDPseudoInverse)

Computes condition number of the matrix for which `B` was constructed.

The condition number is defined as the ratio of the largest over the smallest
singular values of the matrix ``A``.
```math
κ := \\frac{σ_{\\max}}{σ_{\\min}}.
```
"""
function condition_number(B::RegularizedSVDPseudoInverse)
    return B.F.S[1] / B.F.S[end]
end

"""
    solve_via_regularizedSVD(B::RegularizedSVDPseudoInverse, b)

Solve the linear system ``Ax=b`` using the regularized SVD of ``A``.
"""
function solve_via_regularizedSVD(B::RegularizedSVDPseudoInverse, b)
    return B.Vϵ * ( B.iΣϵ * ( adjoint(B.Uϵ) * b ))
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
# To Array
"""
    to_array(draws::Array{Dict{Symbol,T},1}; labels::Array{Symbol,1}  = collect(keys(draws[1]))) where T<:Real
    to_array(dd::DataFrame; labels::Array{Symbol,1} =  Symbol.(names(dd)))
Convert draws or a `DataFrame` to an array and a vector of column labels.
"""
function to_array(draws::Array{Dict{Symbol,T},1}; labels::Array{Symbol,1}  = collect(keys(draws[1]))) where T<:Real
    X = Array{T,2}(undef,length(draws), length(labels))
    for i in 1:length(draws)
        for j in 1:length(labels)
            X[i,j] = draws[i][labels[j]]
        end
    end
    return X, labels
end
function to_array(dd::DataFrame; labels::Array{Symbol,1} =  Symbol.(names(dd)), float_type = typeof(AbstractFloat(1.0)))
    X = Array{float_type,2}(undef,size(dd)[1], length(labels))
    for c in 1:length(labels)
        X[:,c] = dd[!,labels[c]]
    end
    return X, labels
end

# To Draws
"""
    to_draws(X::Array{T,2}; labels::Array{Symbol,1} = Symbol.("x", 1:size(X)[2]))
    to_draws(dd::DataFrame; labels::Array{Symbol,1} =  Symbol.(names(dd)))
Convert array or dataframe to a vector of `Dict`s containing draws.
"""
function to_draws(X::Array{T,2}; labels::Array{Symbol,1} = Symbol.("x", 1:size(X)[2])) where T<:Real
    draws = Array{Dict{Symbol,T},1}()
    for i in 1:size(X)[1]
        draw = Dict{Symbol,T}()
        for j in 1:length(labels)
            draw[labels[j]] = X[i,j]
        end
        draws = vcat(draws, draw)
    end
    return draws
end
function to_draws(dd::DataFrame; labels::Array{Symbol,1} = Symbol.(names(dd)), float_type = typeof(AbstractFloat(1.0)))
    draws = Array{Dict{Symbol,float_type},1}()
    for i in 1:size(dd)[1]
        draw = Dict{Symbol,float_type}()
        for j in labels
            draw[j] = dd[!,j][i]
        end
        draws = vcat(draws, draw)
    end
    return draws
end

# To dataframe
"""
    to_dataframe(X::Array{T,2}; labels::Array{Symbol,1} = Symbol.("x", 1:size(X)[2]))
    to_dataframe(draws::Array{Dict{Symbol,T},1}; labels::Array{Symbol,1}  = collect(keys(draws[1])))
Convert `Array`s or `Vector`s of `Dict`s to a `DataFrame`.
"""
function to_dataframe(X::Array{T,2}; labels::Array{Symbol,1} = Symbol.("x", 1:size(X)[2])) where T<:Real
    dd = DataFrame()
        for j in 1:length(labels)
            dd[!,labels[j]] = X[:,j]
        end
    return dd
end
function to_dataframe(draws::Array{Dict{Symbol,T},1}; labels::Array{Symbol,1}  = collect(keys(draws[1]))) where T<:Real
    dd = DataFrame()
    X = Array{T,2}(undef,length(draws), length(labels))
    for lab in labels
        q = Array{T,1}(undef, length(draws))
        for j in 1:length(draws)
            q[j] = draws[j][lab]
        end
        dd[!,lab] = q
    end
    return dd
end


text:
function isenabled(controller::AbstractControl)
    return controller.isenabled
end

function enable!(controller::AbstractControl)
    controller.isenabled = true
end

function disable!(controller::AbstractControl)
    controller.isenabled  = false
end


text:
using Pkg

Pkg.develop(PackageSpec(path=joinpath(dirname(@__FILE__), "..")))
# when first running instantiate
Pkg.instantiate()

using Documenter, ProxAL

makedocs(
    sitename = "ProxAL.jl",
    format = Documenter.HTML(
        prettyurls = Base.get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.KaTeX()
    ),
    strict = true,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Formulation" => "man/formulation.md",
            "Algorithm" => "man/algorithm.md",
            "Usage" => "man/usage.md",
        ],
        "Library" => [
            "Model parameters" => "lib/modelparams.md",
            "Algorithm parameters" => "lib/algparams.md",
            "Main functions" => "lib/algorithm.md",
            "NLP blocks and backends" => "lib/backends.md",
            "Optimal Power Flow" => "lib/opf.md",
            "MPI communication" => "lib/mpi.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/exanauts/ProxAL.jl.git",
    target = "build",
    devbranch = "develop",
    devurl = "dev",
    push_preview = true,
)

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
# Array info
for (m, f) in ((:Base, :size), (:Base, :axes), (:Base, :firstindex), (:Base, :lastindex))
    @eval begin
        @inline $m.$f(A::AbstractDimArray, dims::AllDims) = $m.$f(A, dimnum(A, dims))
    end
end

# Reducing methods

# With a function arg version
for (m, f) in ((:Base, :sum), (:Base, :prod), (:Base, :maximum), (:Base, :minimum),
                     (:Base, :extrema), (:Statistics, :mean))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; dims=:, kw...) = $_f(A::AbstractDimArray, dims; kw...)
        $m.$f(f, A::AbstractDimArray; dims=:, kw...) =
            $_f(f, A::AbstractDimArray, dims; kw...)
        # Local dispatch methods
        # - Return a reduced DimArray
        $_f(A::AbstractDimArray, dims; kw...) =
            rebuild(A, $m.$f(parent(A); dims=dimnum(A, _astuple(dims)), kw...), reducedims(A, dims))
        $_f(f, A::AbstractDimArray, dims; kw...) =
            rebuild(A, $m.$f(f, parent(A); dims=dimnum(A, _astuple(dims)), kw...), reducedims(A, dims))
        # - Return a scalar
        $_f(A::AbstractDimArray, dims::Colon; kw...) = $m.$f(parent(A); kw...)
        $_f(f, A::AbstractDimArray, dims::Colon; kw...) =
            $m.$f(f, parent(A); dims=dims, kw...)
    end
end
# With no function arg version
for (m, f) in ((:Statistics, :std), (:Statistics, :var))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; corrected::Bool=true, mean=nothing, dims=:) =
            $_f(A, corrected, mean, dims)
        # Local dispatch methods - Returns a reduced array
        $_f(A::AbstractDimArray, corrected, mean, dims) =
            rebuild(A, $m.$f(parent(A); corrected=corrected, mean=mean, dims=dimnum(A, _astuple(dims))), reducedims(A, dims))
        # - Returns a scalar
        $_f(A::AbstractDimArray, corrected, mean, dims::Colon) =
            $m.$f(parent(A); corrected=corrected, mean=mean, dims=:)
    end
end
for (m, f) in ((:Statistics, :median), (:Base, :any), (:Base, :all))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; dims=:) = $_f(A, dims)
        # Local dispatch methods - Returns a reduced array
        $_f(A::AbstractDimArray, dims) =
            rebuild(A, $m.$f(parent(A); dims=dimnum(A, _astuple(dims))), reducedims(A, dims))
        # - Returns a scalar
        $_f(A::AbstractDimArray, dims::Colon) = $m.$f(parent(A); dims=:)
    end
end

# These are not exported but it makes a lot of things easier using them
function Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractDimArray, dims)
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, _astuple(dims))), reducedims(A, dims))
end
function Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractDimArray, dims::Colon)
    Base._mapreduce_dim(f, op, nt, parent(A), dims)
end
function Base._mapreduce_dim(f, op, nt, A::AbstractDimArray, dims)
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
end
@static if VERSION >= v"1.6"
    function Base._mapreduce_dim(f, op, nt::Base._InitialValue, A::AbstractDimArray, dims)
        rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
    end
    function Base._mapreduce_dim(f, op, nt::Base._InitialValue, A::AbstractDimArray, dims::Colon)
        Base._mapreduce_dim(f, op, nt, parent(A), dims)
    end
end

# TODO: Unfortunately Base/accumulate.jl kw methods all force dims to be Integer.
# accumulate wont work unless that is relaxed, or we copy half of the file here.
# Base._accumulate!(op, B, A, dims::AllDims, init::Union{Nothing, Some}) =
    # Base._accumulate!(op, B, A, dimnum(A, dims), init)

# Dimension dropping

function Base.dropdims(A::AbstractDimArray; dims)
    dims = DD.dims(A, dims)
    data = Base.dropdims(parent(A); dims=dimnum(A, dims))
    rebuildsliced(A, data, _dropinds(A, dims))
end

@inline _dropinds(A, dims::DimTuple) = dims2indices(A, map(d -> basetypeof(d)(1), dims))
@inline _dropinds(A, dim::Dimension) = dims2indices(A, basetypeof(dim)(1))

# Function application

@inline Base.map(f, A::AbstractDimArray) = rebuild(A, map(f, parent(A)))

function Base.mapslices(f, A::AbstractDimArray; dims=1, kw...)
    dimnums = dimnum(A, _astuple(dims))
    _data = mapslices(f, parent(A); dims=dimnums, kw...)
    rebuild(A, _data, reducedims(A, DimensionalData.dims(A, dimnums)))
end

# This is copied from base as we can't efficiently wrap this function
# through the kw with a rebuild in the generator. Doing it this way
# also makes it faster to use a dim than an integer.
function Base.eachslice(A::AbstractDimArray; dims=1, kw...)
    if dims isa Tuple && length(dims) != 1
        throw(ArgumentError("only single dimensions are supported"))
    end
    dim = first(dimnum(A, dims))
    dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
    return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
end

# Duplicated dims

for fname in (:cor, :cov)
    @eval function Statistics.$fname(A::AbstractDimArray{<:Any,2}; dims=1, kw...)
        newdata = Statistics.$fname(parent(A); dims=dimnum(A, dims), kw...)
        removed_idx = dimnum(A, dims)
        newrefdims = $dims(A)[removed_idx]
        newdim = $dims(A)[3 - removed_idx]
        rebuild(A, newdata, (newdim, newdim), (newrefdims,))
    end
end

# Rotations

struct _Rot90 end
struct _Rot180 end
struct _Rot270 end
struct _Rot360 end

Base.rotl90(A::AbstractDimMatrix) = rebuild(A, rotl90(parent(A)), _rot(_Rot90(), dims(A)))
Base.rotl90(A::AbstractDimMatrix, k::Integer) =
    rebuild(A, rotl90(parent(A), k), _rot(_rottype(k), dims(A)))

Base.rotr90(A::AbstractDimMatrix) = rebuild(A, rotr90(parent(A)), _rot(_Rot270(), dims(A)))
Base.rotr90(A::AbstractDimMatrix, k::Integer) =
    rebuild(A, rotr90(parent(A), k), _rot(_rottype(-k), dims(A)))

Base.rot180(A::AbstractDimMatrix) = rebuild(A, rot180(parent(A)), _rot(_Rot180(), dims(A)))

# Not type stable - but we have to lose type stability somewhere when
# dims are being swapped, by an Int value, so it may as well be here
function _rottype(k)
    k = mod(k, 4)
    if k == 1
        _Rot90()
    elseif k == 2
        _Rot180()
    elseif k == 3
        _Rot270()
    else
        _Rot360()
    end
end

_rot(::_Rot90, (dima, dimb)) = (flip(Relation, dimb), dima)
_rot(::_Rot180, dims) = map(d -> flip(Relation, d), dims)
_rot(::_Rot270, (dima, dimb)) = (dimb, flip(Relation, dima))
_rot(::_Rot360, dims) = dims

# Dimension reordering

for (pkg, fname) in [(:Base, :permutedims), (:Base, :adjoint),
                     (:Base, :transpose), (:LinearAlgebra, :Transpose)]
    @eval begin
        @inline $pkg.$fname(A::AbstractDimArray{<:Any,2}) =
            rebuild(A, $pkg.$fname(parent(A)), reverse(dims(A)))
        @inline $pkg.$fname(A::AbstractDimArray{<:Any,1}) =
            rebuild(A, $pkg.$fname(parent(A)), (AnonDim(Base.OneTo(1)), dims(A)...))
    end
end
for fname in [:permutedims, :PermutedDimsArray]
    @eval begin
        @inline function Base.$fname(A::AbstractDimArray, perm)
            rebuild(A, $fname(parent(A), dimnum(A, Tuple(perm))), sortdims(dims(A), Tuple(perm)))
        end
    end
end

# Concatenation
function Base._cat(catdims::Union{Int,Base.Dims}, Xin::AbstractDimArray...)
    Base._cat(dims(first(Xin), catdims), Xin...)
end
function Base._cat(catdims::AllDims, Xin::AbstractDimArray...)
    A1 = first(Xin)
    comparedims(Xin...)
    if all(hasdim(A1, catdims))
        # Concatenate an existing dim
        dnum = dimnum(A1, catdims)
        # cat the catdim, ignore others
        newdims = Tuple(_catifcatdim(catdims, ds) for ds in zip(map(dims, Xin)...))
    else
        # Concatenate a new dim
        add_dims = if (catdims isa Tuple)
            Tuple(d for d in catdims if !hasdim(A1, d))
        else
            (catdims,)
        end
        dnum = ndims(A1) + length(add_dims)
        newdims = (dims(A1)..., add_dims...)
    end
    newA = Base._cat(dnum, map(data, Xin)...)
    rebuild(A1, newA, formatdims(newA, newdims))
end

function _catifcatdim(catdims::Tuple, dims)
    anydimsmatch = any(map(d -> basetypeof(d) <: basetypeof(dims[1]), catdims))
    anydimsmatch ? vcat(dims...) : dims[1]
end
function _catifcatdim(catdim, dims)
    basetypeof(catdim) <: basetypeof(dims[1]) ? vcat(dims...) : dims[1]
end

function Base.vcat(dims::Dimension...)
    newmode = _vcat_modes(mode(dims)...)
    rebuild(dims[1], _vcat_index(newmode, index(dims)...), newmode)
end

# IndexModes may need adjustment for `cat`
_vcat_modes(modes::IndexMode...) = first(modes)
function _vcat_modes(modes::AbstractSampled...)
    _vcat_modes(sampling(first(modes)), span(first(modes)), modes...)
end
function _vcat_modes(::Any, ::Regular, modes...)
    _step = step(first(modes))
    map(modes) do mode
        step(span(mode)) == _step || error("Step sizes $(step(span(mode))) and $_step do not match ")
    end
    first(modes)
end
function _vcat_modes(::Intervals, ::Irregular, modes...)
    allbounds = map(bounds ∘ span, modes)
    newbounds = minimum(map(first, allbounds)), maximum(map(last, allbounds))
    rebuild(modes[1]; span=Irregular(newbounds))
end
_vcat_modes(::Points, ::Irregular, modes...) =
    rebuild(first(modes); span=Irregular(nothing, nothing))

# Index vcat depends on mode: NoIndex is always Colon()
_vcat_index(mode::NoIndex, A...) = OneTo(sum(map(length, A)))
# TODO: handle vcat OffsetArrays?
# Otherwise just vcat. TODO: handle order breaking vcat?
_vcat_index(mode::IndexMode, A...) = vcat(A...)

Base.inv(A::AbstractDimArray{T,2}) where T =
    rebuild(A, inv(parent(A)), reverse(map(d -> flip(IndexOrder, d), dims(A))))

# Index breaking

# TODO: change the index and traits of the reduced dimension and return a DimArray.
Base.unique(A::AbstractDimArray; dims::Union{DimOrDimType,Colon}=:) = _unique(A, dims)
Base.unique(A::AbstractDimArray{<:Any,1}) = unique(parent(A))

_unique(A::AbstractDimArray, dims::DimOrDimType) =
    unique(parent(A); dims=dimnum(A, dims))
_unique(A::AbstractDimArray, dims::Colon) = unique(parent(A); dims=:)

Base.diff(A::AbstractDimVector; dims=1) = _diff(A, dimnum(A, dims))
Base.diff(A::AbstractDimArray; dims) = _diff(A, dimnum(A, dims))

@inline function _diff(A::AbstractDimArray{<:Any,N}, dims::Integer) where {N}
    r = axes(A)
    # Copied from Base.diff
    r0 = ntuple(i -> i == dims ? UnitRange(1, last(r[i]) - 1) : UnitRange(r[i]), N)
    rebuildsliced(A, diff(parent(A); dims=dimnum(A, dims)), r0)
end


text:
using FractionalDiffEq, Plots, SpecialFunctions

fx0(x) = 0
function fgz(q)
    x = q[1, 2];t=q[1, 1];a=q[1, 3]
    f = exp(x)*t^(6-a)/gamma(7-a)*720
    return f
end

function f0t(k)
    x=k[1, 1]
    a=k[1, 2]
    f=x^6
    return f
end

function flt(k)
    x = k[1,1];a=k[1,2]
    f = exp(1)*x^6
    return f
end

result = solve(3, 0.5, 1, 20, 1, 20, 1, 1, fx0, fgz, f0t, flt, ADV_DIF())
XX, YY = meshgrid(0.05^2/6 .*(0:21), 0:0.05:1)
plotlyjs()
plot(XX, YY, result, st=:surface)