text:
const DEFAULT_REF_TYPE = Int32

# A wrapper only used for resolving method ambiguity when constructing ScaledArray
mutable struct RefArray{R}
    a::R
end

"""
    ScaledArray{T,R,N,RA,P} <: AbstractArray{T,N}

Array type that stores data as indices of a range.

# Fields
- `refs::RA<:AbstractArray{R,N}`: an array of indices.
- `pool::P<:AbstractRange`: a range that covers all possible values stored by the array.
- `invpool::Dict{T,R}`: a map from array elements to indices of `pool`.
"""
mutable struct ScaledArray{T,R,N,RA,P} <: AbstractArray{T,N}
    refs::RA
    pool::P
    invpool::Dict{T,R}
    function ScaledArray{T,R,N,RA,P}(rs::RefArray{RA}, pool::P, invpool::Dict{T,R}) where
            {T, R, N, RA<:AbstractArray{R, N}, P<:AbstractRange}
        eltype(P) == nonmissingtype(T) || throw(ArgumentError(
            "expect element type of pool being $(nonmissingtype(T)); got $(eltype(P))"))
        return new{T,R,N,RA,P}(rs.a, pool, invpool)
    end
end

ScaledArray(rs::RefArray{RA}, pool::P, invpool::Dict{T,R}) where {T,R,RA<:AbstractArray{R},P} =
    ScaledArray{T,R,ndims(RA),RA,P}(rs, pool, invpool)

const ScaledVector{T,R} = ScaledArray{T,R,1}
const ScaledMatrix{T,R} = ScaledArray{T,R,2}

const ScaledArrOrSub{T,R,N,RA,P} = Union{ScaledArray{T,R,N,RA,P},
    SubArray{<:Any, <:Any, <:ScaledArray{T,R,N,RA,P}}}

"""
    scale(sa::ScaledArrOrSub)

Return the step size of the `pool` of `sa`.
"""
scale(sa::ScaledArrOrSub) = step(DataAPI.refpool(sa))

Base.size(sa::ScaledArray) = size(sa.refs)
Base.IndexStyle(::Type{<:ScaledArray{T,R,N,RA}}) where {T,R,N,RA} = IndexStyle(RA)

function _validmin(min, xmin, isstart::Bool)
    if min === nothing
        min = xmin
    elseif min > xmin
        str = isstart ? "start =" : "stop ="
        throw(ArgumentError("$str $min is greater than the allowed minimum element $xmin"))
    end
    return min
end

function _validmax(max, xmax, isstart::Bool)
    if max === nothing
        max = xmax
    elseif max < xmax
        str = isstart ? "start =" : "stop ="
        throw(ArgumentError("$str $max is smaller than the allowed maximum element $xmax"))
    end
    return max
end

function validpool(x::AbstractArray, T::Type, start, step, stop, usepool::Bool)
    step === nothing && throw(ArgumentError("step cannot be nothing"))
    pool = DataAPI.refpool(x)
    xs = skipmissing(usepool && pool !== nothing ? pool : x)
    xmin, xmax = extrema(xs)
    applicable(+, xmin, step) || throw(ArgumentError(
        "step of type $(typeof(step)) does not match array with element type $(eltype(x))"))
    if xmin + step > xmin
        start = _validmin(start, xmin, true)
        stop = _validmax(stop, xmax, false)
    elseif xmin + step < xmin
        start = _validmax(start, xmax, true)
        stop = _validmin(stop, xmin, false)
    else
        throw(ArgumentError("step cannot be zero"))
    end
    start = convert(T, start)
    stop = convert(T, stop)
    return start:step:stop
end

function _scaledlabel!(labels::AbstractArray, invpool::Dict, xs::AbstractArray, start, step)
    z = zero(valtype(invpool))
    @inbounds for (i, x) in zip(eachindex(labels), xs)
        lbl = get(invpool, x, z)
        if lbl !== z
            labels[i] = lbl
        elseif ismissing(x)
            labels[i] = z
            invpool[x] = z
        else
            r = start:step:x
            lbl = length(r)
            labels[i] = lbl
            invpool[x] = lbl
        end
    end
end

function scaledlabel(xs::AbstractArray, stepsize,
        R::Type=DEFAULT_REF_TYPE, T::Type=eltype(xs);
        start=nothing, stop=nothing, usepool::Bool=true)
    pool = validpool(xs, T, start, stepsize, stop, usepool)
    T = Missing <: T ? Union{eltype(pool), Missing} : eltype(pool)
    start = first(pool)
    stepsize = step(pool)
    if R <: Integer
        while typemax(R) < length(pool)
            R = widen(R)
        end
    end
    # Array types with customized indexing could cause problems
    labels = similar(Array{R}, axes(xs))
    invpool = Dict{T,R}()
    _scaledlabel!(labels, invpool, xs, start, stepsize)
    return labels, pool, invpool
end

function ScaledArray(x::AbstractArray, reftype::Type, xtype::Type, start, step, stop, usepool::Bool=true)
    refs, pool, invpool = scaledlabel(x, step, reftype, xtype; start=start, stop=stop, usepool=usepool)
    return ScaledArray(RefArray(refs), pool, invpool)
end

function ScaledArray(sa::ScaledArray, reftype::Type, xtype::Type, start, step, stop, usepool::Bool=true)
    if step !== nothing && step != scale(sa)
        refs, pool, invpool = scaledlabel(sa, step, reftype, xtype; start=start, stop=stop, usepool=usepool)
        return ScaledArray(RefArray(refs), pool, invpool)
    else
        step = scale(sa)
        pool = validpool(sa, xtype, start, step, stop, usepool)
        T = Missing <: xtype ? Union{eltype(pool), Missing} : eltype(pool)
        refs = similar(sa.refs, reftype)
        invpool = Dict{T, reftype}()
        start0 = first(sa.pool)
        start = first(pool)
        stop = last(pool)
        if start == start0
            copy!(refs, sa.refs)
            copy!(invpool, sa.invpool)
        elseif start < start0 && start < stop || start > start0 && start > stop
            offset = length(start:step:start0) - 1
            refs .= sa.refs .+ offset
            for (k, v) in sa.invpool
                invpool[k] = v + offset
            end
        else
            offset = length(start0:step:start) - 1
            refs .= sa.refs .- offset
            for (k, v) in sa.invpool
                invpool[k] = v - offset
            end
        end
        return ScaledArray(RefArray(refs), pool, invpool)
    end
end

"""
    ScaledArray(x::AbstractArray, start, step[, stop]; reftype=Int32, usepool=true)
    ScaledArray(x::AbstractArray, step; reftype=Int32, start, stop, usepool=true)

Construct a `ScaledArray` from `x` given the `step` size.
If `start` or `stop` is not specified, it will be chosen based on the extrema of `x`.

# Keywords
- `reftype::Type=Int32`: the element type of field `refs`.
- `usepool::Bool=true`: find extrema of `x` based on `DataAPI.refpool`.
"""
ScaledArray(x::AbstractArray, start, step, stop=nothing;
    reftype::Type=DEFAULT_REF_TYPE, xtype::Type=eltype(x), usepool::Bool=true) =
        ScaledArray(x, reftype, xtype, start, step, stop, usepool)

ScaledArray(sa::ScaledArray, start, step, stop=nothing;
    reftype::Type=eltype(refarray(sa)), xtype::Type=eltype(sa), usepool::Bool=true) =
        ScaledArray(sa, reftype, xtype, start, step, stop, usepool)

ScaledArray(x::AbstractArray, step; reftype::Type=DEFAULT_REF_TYPE,
    start=nothing, stop=nothing, xtype::Type=eltype(x), usepool::Bool=true) =
        ScaledArray(x, reftype, xtype, start, step, stop, usepool)

ScaledArray(sa::ScaledArray, step=nothing; reftype::Type=eltype(refarray(sa)),
    start=nothing, stop=nothing, xtype::Type=eltype(sa), usepool::Bool=true) =
        ScaledArray(sa, reftype, xtype, start, step, stop, usepool)

Base.similar(sa::ScaledArrOrSub{T,R}, dims::Dims=size(sa)) where {T,R} =
    ScaledArray(RefArray(ones(R, dims)), DataAPI.refpool(sa), Dict{T,R}())

Base.similar(sa::ScaledArrOrSub, dims::Int...) = similar(sa, dims)

"""
    align(xs::AbstractArray, sa::ScaledArrOrSub)

Convert `xs` into a [`ScaledArray`](@ref) with a `pool`
that has the same first element and step size as the `pool` from `sa`.
"""
function align(xs::AbstractArray, sa::ScaledArrOrSub)
    pool = DataAPI.refpool(sa)
    invpool = DataAPI.invrefpool(sa)
    step = scale(sa)
    xmin, xmax = extrema(xs)
    start = first(pool)
    stop = last(pool)
    start < stop && xmin < start && throw(ArgumentError(
        "the minimum of xs $xmin is smaller than the minimum of pool $start"))
    start > stop && xmax > start && throw(ArgumentError(
        "the maximum of xs $xmax is greater than the maximum of pool $start"))
    refs = similar(DataAPI.refarray(sa), size(xs))
    _scaledlabel!(refs, invpool, xs, start, step)
    return ScaledArray(RefArray(refs), pool, invpool)
end

DataAPI.refarray(sa::ScaledArray) = sa.refs
DataAPI.refvalue(sa::ScaledArray, n::Integer) = getindex(DataAPI.refpool(sa), n)
DataAPI.refpool(sa::ScaledArray) = sa.pool
DataAPI.invrefpool(sa::ScaledArray) = sa.invpool

DataAPI.refarray(ssa::SubArray{<:Any, <:Any, <:ScaledArray}) =
    view(parent(ssa).refs, ssa.indices...)
DataAPI.refvalue(ssa::SubArray{<:Any, <:Any, <:ScaledArray}, n::Integer) =
    DataAPI.refvalue(parent(ssa), n)
DataAPI.refpool(ssa::SubArray{<:Any, <:Any, <:ScaledArray}) =
    DataAPI.refpool(parent(ssa))
DataAPI.invrefpool(ssa::SubArray{<:Any, <:Any, <:ScaledArray}) =
    DataAPI.invrefpool(parent(ssa))

@inline function Base.getindex(sa::ScaledArray, i::Int)
    refs = DataAPI.refarray(sa)
    @boundscheck checkbounds(refs, i)
    @inbounds n = refs[i]
    iszero(n) && return missing
    pool = DataAPI.refpool(sa)
    @boundscheck checkbounds(pool, n)
    return @inbounds pool[n]
end

Base.@propagate_inbounds function Base.getindex(sa::ScaledArrOrSub, I::AbstractVector)
    newrefs = DataAPI.refarray(sa)[I]
    pool = DataAPI.refpool(sa)
    invpool = DataAPI.invrefpool(sa)
    return ScaledArray(RefArray(newrefs), pool, invpool)
end

Base.@propagate_inbounds function Base.setindex!(sa::ScaledArray, val, ind::Int)
    invpool = DataAPI.invrefpool(sa)
    n = get(invpool, val, nothing)
    if n === nothing
        pool = DataAPI.refpool(sa)
        r = first(pool):step(pool):val
        last(r) > last(pool) && (sa.pool = r)
        n = length(r)
        invpool[val] = n
    end
    refs = DataAPI.refarray(sa)
    refs[ind] = n
    return sa
end

Base.@propagate_inbounds function Base.setindex!(sa::ScaledArray, ::Missing, ind::Int)
    refs = DataAPI.refarray(sa)
    z = zero(eltype(refs))
    invpool = DataAPI.invrefpool(sa)
    invpool[missing] = z
    refs[ind] = z
    return sa
end

allowmissing(sa::ScaledArray{T,R}) where {T,R} =
    ScaledArray(RefArray(sa.refs), sa.pool, convert(Dict{Union{T,Missing},R}, sa.invpool))

function disallowmissing(sa::ScaledArray{T,R}) where {T,R}
    T1 = nonmissingtype(T)
    any(x->iszero(x), sa.refs) && throw(ArgumentError("cannot convert missing to $T1"))
    delete!(sa.invpool, missing)
    return ScaledArray(RefArray(sa.refs), sa.pool, convert(Dict{T1,R}, sa.invpool))
end

function Base.:(==)(x::ScaledArrOrSub, y::ScaledArrOrSub)
    size(x) == size(y) || return false
    first(DataAPI.refpool(x)) == first(DataAPI.refpool(y)) &&
        scale(x) == scale(y) && return DataAPI.refarray(x) == DataAPI.refarray(y)
    eq = true
    for (p, q) in zip(x, y)
        # missing could arise
        veq = p == q
        isequal(veq, false) && return false
        eq &= veq
    end
    return eq
end


text:
function format_inner(inner)
    first_index = findfirst(isequal(' '), inner) + 1
    last_index = findlast(isequal(' '), inner) - 1
    return inner[first_index:last_index]
end

function parse_rules(s)
    rules = Dict{String,Vector{String}}()
    for line in readlines(IOBuffer(s))
        outer, inners = split(line, " bags contain ")
        if inners == "no other bags."
            continue
        end
        for inner in split(inners, ", ")
            formatted_inner = format_inner(inner)
            if formatted_inner in keys(rules)
                push!(rules[formatted_inner], outer)
            else
                rules[formatted_inner] = [outer,]
            end
        end
    end
    return rules
end
function fill_can_be_included_in(rules, bag, acc)
    push!(acc, bag)
    if bag in keys(rules)
        for next in rules[bag]
            if !(next in acc)
                fill_can_be_included_in(rules, next, acc)
            end
        end
    end
end

function run(s)
    # Your code here
    rules = parse_rules(s)

    acc = Set()
    fill_can_be_included_in(rules, "shiny gold", acc)
    return length(acc) - 1
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()

ext: jl


text:
"""
    alloc_mem(memsize::UnitPrefixedBytes; devs=(CUDA.device(),), dtype=Float32)
Allocates memory on the devices whose IDs are provided via `devs`.
Returns a vector of memory handles (i.e. `CuArray`s).

**Examples:**
```julia
alloc_mem(MiB(1024)) # allocate on the currently active device
alloc_mem(B(40_000_000); devs=(0,1)) # allocate on GPU0 and GPU1
```
"""
function alloc_mem(memsize::UnitPrefixedBytes; devs=(CUDA.device(),), dtype=Float32)
    all(map(isequal, devs, unique(devs))) ||
        throw(ArgumentError("No duplicates in `devs` allowed."))
    N = Int(bytes(memsize) ÷ sizeof(dtype))

    mem_handles = Vector{CuArray{dtype,1,CUDA.Mem.DeviceBuffer}}(undef, length(devs))
    for (i, dev) in pairs(devs)
        device!(dev)
        mem_handles[i] = CUDA.rand(N)
    end
    return mem_handles
end

"Reclaim the unused memory of the currently active GPU (i.e. `device()`)."
function clear_gpu_memory(device::CuDevice=CUDA.device(); gc=true)
    device!(device) do
        gc && GC.gc()
        CUDA.reclaim()
    end
    return nothing
end

"Reclaim the unused memory of all available GPUs."
function clear_all_gpus_memory(devices=CUDA.devices(); gc=true)
    gc && GC.gc()
    for dev in devices
        clear_gpu_memory(dev; gc=false)
    end
    return nothing
end

"""
Check if CUDA/GPU is available and functional.
If not, print some (hopefully useful) debug information.
"""
function functional(verbose=true)
    if CUDA.functional()
        verbose && @info("CUDA/GPU available.")
        hascuda = true
    else
        verbose && @info("No CUDA/GPU found.")
        hascuda = false
        if verbose
            # debug information
            @show Libdl.find_library("libcuda")
            @show filter(contains("cuda"), lowercase.(Libdl.dllist()))
            try
                @info("CUDA.versioninfo():")
                CUDA.versioninfo()
                @info("Successful!")
            catch ex
                @warn("Unsuccessful!", ex)
                println()
            end
        end
    end
    return hascuda
end

"""
    toggle_tensorcoremath([enable::Bool; verbose=true])
Switches the `CUDA.math_mode` between `CUDA.FAST_MATH` (`enable=true`) and `CUDA.DEFAULT_MATH` (`enable=false`).
For matmuls of `CuArray{Float32}`s, this should have the effect of using/enabling and not using/disabling tensor cores.
Of course, this only works on supported devices and CUDA versions.

If no arguments are provided, this functions toggles between the two math modes.
"""
function toggle_tensorcoremath(
    enable::Bool=CUDA.math_mode() == CUDA.DEFAULT_MATH; verbose=true
)
    if enable
        CUDA.math_mode!(CUDA.FAST_MATH)
        verbose && @info(
            "Tensor cores should be ACTIVE for matmuls (only involving `Float32` elements).",
            CUDA.math_mode()
        )
    else
        CUDA.math_mode!(CUDA.DEFAULT_MATH)
        verbose && @info(
            "Tensor cores should be NOT ACTIVE for matmuls (only involving `Float32` elements).",
            CUDA.math_mode()
        )
    end
    return nothing
end

_device2string(dev::CuDevice) = "GPU $(gpuid(dev)): $(CUDA.name(dev))"

L2_cachesize() = cachesize()[2]

"""
Tries to get the temperatures of the available CPUs (sockets not cores) in degrees Celsius.

Based on `cat /sys/class/thermal/thermal_zone*/temp`.
"""
function get_cpusocket_temperatures()
    temp_strs = split(
        read(`cat $(glob("thermal_zone*", "/sys/class/thermal"))/temp`, String)
    )
    return parse.(Float64, temp_strs) ./ 1000
end

"""
Get information about all cpu cores. Returns a vector of vectors. The outer index corresponds
to cpu cores. The inner vector contains the following information (in that order):

`user nice system idle iowait irq softirq steal guest ?`

See [proc(5)](https://man7.org/linux/man-pages/man5/proc.5.html) for more information.
"""
function get_cpu_stats()
    st = read("/proc/stat", String)
    cpulines = filter(l -> startswith(l, "cpu") && !startswith(l, "cpu "), split(st, '\n'))
    stats = Vector{Vector{Int64}}(undef, length(cpulines))
    for (i, line) in enumerate(cpulines)
        stats[i] = @views parse.(Int64, split(line)[2:end])
    end
    # user nice system idle iowait irq softirq steal guest ?
    return stats
end

"""
    get_cpu_utilizations(cores=0:Sys.CPU_THREADS-1; Δt=0.01)
Get the utilization (in percent) of the given cpu `cores` over a certain time interval `Δt`.

Based on [this](https://www.idnt.net/en-US/kb/941772).
"""
function get_cpu_utilizations(cores=0:(Sys.CPU_THREADS - 1); Δt=0.01)
    stats1 = get_cpu_stats()
    sleep(Δt)
    stats2 = get_cpu_stats()
    coreidcs = cores .+ 1 # + 1 for coreids -> coreindex mapping
    idles1, totals1 = _stats2idles_and_totals(stats1[coreidcs])
    idles2, totals2 = _stats2idles_and_totals(stats2[coreidcs])
    Δtotals = totals2 .- totals1
    Δidles = idles2 .- idles1
    return 100.0 .* (Δtotals .- Δidles) ./ Δtotals
end

"""
    get_cpu_utilization(core=getcpuid(); Δt=0.01)
Get the utilization (in percent) of the given cpu `core` over a certain time interval `Δt`.
"""
function get_cpu_utilization(core=getcpuid(); kwargs...)
    return only(get_cpu_utilizations([core]; kwargs...))
end

function _stats2idles_and_totals(stats=get_cpu_stats())
    idles = zeros(length(stats))
    totals = zeros(length(stats))
    for i in eachindex(stats)
        vals = stats[i]
        idles[i] = vals[4] + vals[5] # idle + iowait
        nonidle = vals[1] + vals[2] + vals[3] + vals[6] + vals[7] + vals[8] # user nice system irq softirq steal
        totals[i] = idles[i] + nonidle
    end
    return idles, totals
end

"""
Checks whether the given `CuDevice` has Tensor Cores.
"""
function hastensorcores(dev::CuDevice=CUDA.device())
    cap = CUDA.capability(dev)
    if cap == v"8.0.0"
        # A100
        return true
    elseif cap == v"7.0.0"
        # V100
        return true
    else
        return false
    end
end

function logspace(start, stop, length)
    return exp2.(range(log2(start), log2(stop); length=length))
end

function pkgversion(pkg::Module)
    return Pkg.Types.read_project(joinpath(Base.pkgdir(pkg), "Project.toml")).version
end

"""
MultiLogger struct which combines normal and error output streams.
Useful if you want your normal and error output that is printed to the terminal to also be saved to different files.

**Arguments:**
* `normal_file_name`: Path to normal output file.
* `error_file_name`: Path to error output file.
"""
struct MultiLogger
    normal_io::IO
    error_io::IO
    normal_logger::ConsoleLogger
    error_logger::ConsoleLogger

    function MultiLogger(normal_file_name::String, error_file_name::String)
        normal_io = open(normal_file_name, "w+")
        normal_logger = ConsoleLogger(normal_io)

        error_io = open(error_file_name, "w+")
        error_logger = ConsoleLogger(error_io)

        new(normal_io, error_io, normal_logger, error_logger)
    end
end

"""
Logging function for MultiLogger. Use this in combination with MultiLogger if you want
your normal and error output that is printed to the terminal to also be saved to different files.

**Arguments:**
* `MultiLogger`: Instance of MultiLogger struct.
* `text`: Text to be printed to terminal and written to file.
* `is_error` (default: `false`): Flag which decides whether `text` should be written into normal or error stream.
"""
function multi_log(multilogger::MultiLogger, text::String, is_error::Bool=false)

      @info(text)
      println("")

      if !is_error
          with_logger(multilogger.normal_logger) do
            @info(text)
            println("")
          end
          flush(multilogger.normal_io)
      else
          with_logger(multilogger.error_logger) do
            @error(text)
            println("")
          end
          flush(multilogger.error_io)
      end
end


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

mutable struct IoK8sApiAppsV1beta2RollingUpdateDaemonSet <: SwaggerModel
    maxUnavailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxUnavailable

    function IoK8sApiAppsV1beta2RollingUpdateDaemonSet(;maxUnavailable=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2RollingUpdateDaemonSet, Symbol("maxUnavailable"), maxUnavailable)
        setfield!(o, Symbol("maxUnavailable"), maxUnavailable)
        o
    end
end # type IoK8sApiAppsV1beta2RollingUpdateDaemonSet

const _property_map_IoK8sApiAppsV1beta2RollingUpdateDaemonSet = Dict{Symbol,Symbol}(Symbol("maxUnavailable")=>Symbol("maxUnavailable"))
const _property_types_IoK8sApiAppsV1beta2RollingUpdateDaemonSet = Dict{Symbol,String}(Symbol("maxUnavailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2RollingUpdateDaemonSet }) = collect(keys(_property_map_IoK8sApiAppsV1beta2RollingUpdateDaemonSet))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2RollingUpdateDaemonSet }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2RollingUpdateDaemonSet[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2RollingUpdateDaemonSet }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2RollingUpdateDaemonSet[property_name]

function check_required(o::IoK8sApiAppsV1beta2RollingUpdateDaemonSet)
    true
end

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


text:
mutable struct Todo
    Prefix::String
    Suffix::String
    Filename::String
    Body::String
    Line::Int64
    ID::Int64

    Todo() = new("", "", "", "", -1, -1) # empty initializer
    Todo(pre, suf, fn, body, ln, id) = new(pre, suf, fn, body, ln, id)
    Todo(body) = new("", "", "", body, -1, -1) # initialize only body

end

function Base.:isempty(t::Todo)
    flag = false
    for i in 1:4
        if isempty(getfield(t, i)) && !flag
            flag = true
            break
        end
    end # for

    return flag
end


text:
let

    for (runner_cmd, src_name) in [
            ["julia --startup-file=no ", "test_script.jl"],
            ["bash -c", "test_script.sh"],
        ]

        src_file = joinpath(@__DIR__, src_name)
        outfile = joinpath(tempname(), "out.log")
        GW._rm(outfile)
        GW._mkdir(outfile)
        GW._spawn_bash("$(runner_cmd) $(src_file) 2>&1 | tee $(outfile)")

        cat = ""
        done_flag = false
        for _ in 1:15
            println("-"^60)
            println(basename(src_file))
            if isfile(outfile)
                cat = read(outfile, String)
                println(cat)
            end
            done_flag = contains(cat, "DONE")
            done_flag && break
            sleep(1.0)
        end
        @assert done_flag
        GW._rm(outfile)
    end
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using DistributedEnvironments
using Test

@testset "DistributedEnvironments.jl" begin
    # Write your tests here.
end


text:
# Gradient tests for adjoint extended modeling
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: April 2020
#

using JUDI, SegyIO, JOLI, Flux, JUDI4Flux
using Test, ImageFiltering, LinearAlgebra, Random

Random.seed!(11)

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v[:, Int(round(end/2)):end] .= 4f0
v0 = imfilter(v, Float32.(Kernel.gaussian(10)))

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
time = 1000f0   # receiver recording time [ms]
dt = 1f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# setup wavelet
f0 = 0.01f0     # MHz
wavelet = ricker_wavelet(time, dt, f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

# return linearized data as Julia array
opt = Options(return_array=true, dt_comp=dt, sum_padding=true)

# Linear operators
Pr = judiProjection(info, recGeometry)
A_inv = judiModeling(info, model; options=opt)
Pw = judiLRWF(info, wavelet)
F = Pr*A_inv*adjoint(Pw)

#####################################################################################

G = ExtendedQForward(F)

misfit(x, m, y) = Flux.mse(G(x, m), y)

function loss(x, m, y)
    f = misfit(x, m, y)
    gs = gradient(() -> misfit(x, m, y), params(x, m))
    return f, gs[x], gs[m]
end

# Nonlinear modeling
m = reshape(m, n[1], n[2], 1, 1)
x = randn(Float32, n[1], n[2], 1, 1)
x0 = randn(Float32, n[1], n[2], 1, 1)
dx = x - x0
y = randn(Float32, recGeometry.nt[1], nxrec, 1, 1)

# Gradient test for extended modeling: weights
f0, g = loss(x0, m, y)[1:2]
h = .1f0
maxiter = 6

err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test extended source weights\n")
for j=1:maxiter
    f = loss(x0 + h*dx, m, y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dx, g))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for extended modeling: model
x = randn(Float32, n[1], n[2], 1, 1)
m = reshape(m, n[1], n[2], 1, 1)
m0 = reshape(m0, n[1], n[2], 1, 1)
dm = m - m0

f0, gw, gm = loss(x, m0, y)
h = .1f0
maxiter = 4

err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test extended source model\n")
for j=1:maxiter
    f = loss(x, m0 + h*dm, y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dm, gm))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
using Documenter, BasisFunctionExpansions

makedocs(doctest = false) # Due to lots of plots, this will just have to be run on my local machine

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-cinder"),
    repo = "github.com/baggepinnen/BasisFunctionExpansions.jl.git",
    julia  = "0.6",
    osname = "linux"
)


text:
@generated function gen_harmonic_mean(x...)
    sum = :(0)
    for i = 1:length(x)
        sum = :(1/x[$i] + $sum)
    end
    mean = :(length(x) / $sum)
    return mean
end

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
"""
Convert AbstractRange to list of tuples (start, end, grid)
"""
function ranges2tuples(param_ranges)
    param_ranges_ = Array{Tuple{Number,Number,Integer},1}(undef,length(param_ranges))
    for (i, ran) in enumerate(param_ranges)
        par_rang = (ran[1], ran[end], length(ran))
        param_ranges_[i] = par_rang
    end
    return param_ranges_
end


text:
using Bounds
using Test


text:
function AcmulB!(y::Array{Float64,1}, A::SparseMatrixCSC{Float64,Int64}, x::Array{Float64,1}; nthreads::Int64=4)

    (m,n) = size(A)
    if length(x) != m || length(y) != n
       throw(DimensionMismatch("length(x) != m || length(y) != n"))
    end
    p = ccall((:ac_mul_b_rr_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}),
                                                    &nthreads,  &m,         &n,         A.nzval,      A.rowval,   A.colptr,   x,            y             )
    return nothing
end

function AcmulB!(y::Array{Complex128,1}, A::SparseMatrixCSC{Float64,Int64}, x::Array{Complex128,1}; nthreads::Int64=4)

    (m,n) = size(A)
    if length(x) != m || length(y) != n
       throw(DimensionMismatch("length(x) != m || length(y) != n"))
    end
    p = ccall((:ac_mul_b_rc_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                                                       &nthreads,  &m,         &n,         A.nzval,      A.rowval,   A.colptr,   x,               y                )
    return nothing
end

function AcmulB!(y::Array{Complex128,1}, A::SparseMatrixCSC{Complex128,Int64}, x::Array{Complex128,1}; nthreads::Int64=4)

    (m,n) = size(A)
    if length(x) != m || length(y) != n
       throw(DimensionMismatch("length(x) != m || length(y) != n"))
    end
    p = ccall((:ac_mul_b_cc_, spmatveclib), Int64, (Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Int64}, Ptr{Int64}, Ptr{Complex128}, Ptr{Complex128}),
                                                       &nthreads,  &m,         &n,         A.nzval,         A.rowval,   A.colptr,   x,               y                )
    return nothing
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
module OceanTurbulenceParameterEstimation

export OneDimensionalTimeSeries, InverseProblem, FreeParameters,
        IdentityNormalization, ZScore, forward_map, observation_map,
        eki, lognormal_with_mean_std, iterate!, EnsembleKalmanInversion, UnscentedKalmanInversion,
        UnscentedKalmanInversionPostprocess, ConstrainedNormal

include("Observations.jl")
include("TurbulenceClosureParameters.jl")
include("InverseProblems.jl")
include("EnsembleKalmanInversions.jl")

using .Observations: OneDimensionalTimeSeries, ZScore
using .TurbulenceClosureParameters: FreeParameters
using .InverseProblems: InverseProblem, forward_map, observation_map
using .EnsembleKalmanInversions: iterate!, EnsembleKalmanInversion, UnscentedKalmanInversion,
       UnscentedKalmanInversionPostprocess,  ConstrainedNormal, lognormal_with_mean_std

end # module


text:
using Perceptrons
using Test

include("utils_test.jl")
include("linear_perceptron_test.jl")
include("kernel_perceptron_test.jl")
include("voted_perceptron_test.jl")
include("averaged_perceptron_test.jl")


text:
function Base.show(io::IO, ::MIME"text/plain", result::ConcreteResult)
    print(io, "TryExperimental.ConcreteResult ")
    value = result.value
    if value isa Ok
        printstyled(io, "(Ok)"; color = :green, bold = true)
        print(io, ": ")
        show(io, MIME"text/plain"(), Try.unwrap(value))
    else
        printstyled(io, "(Err)"; color = :red, bold = true)
        print(io, ": ")
        ex = Try.unwrap_err(err)
        backtrace = err.backtrace
        if backtrace === nothing
            showerror(io, ex)
        else
            showerror(io, ex, simplify_backtrace(err.backtrace))
        end
    end
end

# TODO: simplify arguments when they are too long
function Base.showerror(io::IO, ex::NotImplementedError)
    print(io, "Not Implemented: ")
    show(io, ex.f)
    print(io, '(')
    let isfirst = true
        for a in ex.args
            if isfirst
                isfirst = false
            else
                print(io, ", ")
            end
            show(IOContext(io, :compact => true, :limit => true), a)
        end
    end
    let isfirst = true
        for (k, v) in pairs(ex.kwargs)
            if isfirst
                isfirst = false
                print(io, "; ")
            else
                print(io, ", ")
            end
            print(io, k, " = ")
            show(IOContext(io, :compact => true, :limit => true), v)
        end
    end
    print(io, ')')
end


text:
include("../../src/SPICE.jl")
using .SPICE

# ---- State-change matrix
const ν = [1 0; -1 1; 0 -1]

# ---- Initial state-space
x0 = [50,50]

# ---- Highest orders of reaction for tau-leaping
const hor = [2,2]

# ---- Bounds on the initial search space for parameters
# ---- Real parameters are [0.5,0.0025,0.3]
const bounds = [[1e-6,1e-6,1e-6],[10.0,10.0,10.0]]

# ---- Hazard functions
function F(p)
    H(p)
    p.a .*= p.θ
end

function H(p)
    p.a[1] = p.x[1]
    p.a[2] = p.x[1] * p.x[2]
    p.a[3] = p.x[2]
end

system = System(Model(x0,F,H, ν,bounds,hor=hor,obsname=[:X1, :X2]), "./data/5",routine=CEM(ssa=:Tau, nElite = 10, nRepeat=1, nSamples=1000, maxIter=250, mSamples=20000, shoot=false, splitting=false, sampling=:log, tauOpt=TauOpt(ϵ=0.1)))
estimate(system, 1, "lotka")

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


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

using Markdown
using InteractiveUtils

# ╔═╡ 3395c983-7d14-4163-bebb-096c851cbde8
begin
	using Revise
	using Pkg
	Pkg.activate("../")

	using BayesianModels
	using DataFrames
	using DelimitedFiles
	using Downloads
	using ForwardDiff
	using LinearAlgebra
	using Plots
	using Zygote

	include("./plotting.jl")
end

# ╔═╡ e5883656-7aa1-11ec-2648-4f9569fddcf8
md"""
# Normal distribution

Learn a multivariate normal distribution
"""

# ╔═╡ 9b39f491-8836-4368-877b-6cf1eea0ff19
md"""
## Data
"""

# ╔═╡ 60f71c97-f4c5-4b4d-9665-4838d676793b
datapath = Downloads.download("https://www.cs.cmu.edu/afs/cs/project/ai-repository/ai/areas/speech/database/pb/pb.tgz", "data.tar.gz")

# ╔═╡ 52b2c2a9-01c2-4d39-b9d6-4bcd7a11dec1
run(`tar xf $datapath`)

# ╔═╡ 04ed61ae-a1ed-4114-a160-100e86f0769a
Text(read("PetersonBarney/HEADER", String))

# ╔═╡ ff11bde6-a2c8-47f5-a3bf-a4b6f9758747
# We remove the "*" marker.
open("PetersonBarney/verified_pb.data.filt", "w") do f
	write(f, replace(read("PetersonBarney/verified_pb.data", String), '*' => ""))
end

# ╔═╡ 4b071c89-67b6-4964-a4b1-c784f6d99f9c
begin
	rawX = readdlm("PetersonBarney/verified_pb.data.filt", '\t')

	# Normalize each dimension to have 0 mean and 1 standard deviation.
	local X = rawX[:, 5:end]
	local μ̄ = sum(X, dims=1) / size(X,1)
	local σ̄ = sqrt.( sum((X .- μ̄).^2, dims=1) / size(X,1) )
	rawX[:, 5:end] = (X .- μ̄) ./ σ̄

	data = DataFrame(
		rawX,
		[:type, :speaker, :phone_id, :phone_ascii, :f0, :f1, :f2, :f3]
	)
end

# ╔═╡ e2aa337d-4cdc-43d2-a8f6-5d0a607a3de0
arpa2ipa = Dict(
	"IY" => "i",
	"IH" => "I",
	"EH" => "ε",
	"AE" => "æ",
	"AH" => "Λ",
	"AA" => "ɑ",
	"AO" => "ɔ",
	"UH" => "ʊ",
	"UW" => "u",
	"ER" => "ɝ"
)

# ╔═╡ 570dc7f7-10dc-46ed-9958-2acc0ed76ad1
phones = Set(data[:, :phone_ascii])

# ╔═╡ c74b04e1-b19d-4398-b376-dda5fa88ca86
begin
	local p = plot()
	for phone in phones
		X = subset(data,  :phone_ascii => p -> p .== phone)
		X = subset(X, :type => t -> t .== 1)
		X = X[:, [:f1, :f2]]
		scatter!(X[:, 1], X[:, 2], label=phone)
	end
	p
end

# ╔═╡ 15c82689-be5d-411b-9594-1ba031fb2a0d
Xs = [
	Array{Float64}(
		subset(
			data,
			:type => t -> t .== 1,
			:phone_ascii => p -> p .== phone
		)[:, [:f0, :f1, :f2, :f3]])'
	for phone in phones
]

# ╔═╡ 6c8c5ecc-9dc2-446a-8726-699694677c42
md"""
## Basic Model

We start by learning a basic model: a set of independent Normal distributions each with a Normal-Wishart prior.
"""

# ╔═╡ 2d01cdec-0da2-48f7-9ae9-7cf2298e9d11
# Dimension of the observed space.
D = size(Xs[1], 1)

# ╔═╡ 3814457c-fb1b-4f49-9717-c9d9fcdcc070
prior = NormalWishart(zeros(4), 1., Matrix{Float64}(I, 4, 4), 4.)

# ╔═╡ 5fa65b07-fe1c-4119-93fa-85738d74af6c
models = [NormalModel(prior) for phone in phones]

# ╔═╡ c43a677b-58dd-4266-a93a-c914b7cb3ebc
begin
	qθs = [prior for phone in phones]
	local N = sum(size.(Xs, 2))
	L = []

	for t in 1:2000
		μs = getstats.(models, qθs)
		(l, (∇μs,)) = withgradient(μs -> sum(elbo.(models, Xs, qθs, μs))/N, μs)
		qθs = newposterior.(models, qθs, ∇μs; lrate=0.1)

		push!(L, l)
	end
	L
end

# ╔═╡ e8148b69-8fba-4d15-8c41-8e445d2cb5a0
plot(L, label=false)

# ╔═╡ 67ee5a37-0baa-4986-8558-77b720321771
begin
	local p = plot(legend=false)
	for (phone, X, qθ) in zip(phones, Xs, qθs)
		μ = qθ.μ[2:3]
		Σ = inv(qθ.ν*qθ.W)[2:3,2:3]
		plotnormal!(μ, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
		scatter!(X[2,:], X[3,:], label=phone)
	end
	p
end

# ╔═╡ d12fcd36-51b4-463e-b01c-807239522d6b
md"""
## GSM

Now, we build the GSM. We start by defining the parameter transformation from the "real space" ($\boldsymbol{\psi}$) to the "standard space" ($\boldsymbol{\Lambda}$, $\boldsymbol{\mu}$). The transformation is defined by:

```math
\begin{align}
	\boldsymbol{\psi} &= \begin{bmatrix}
		\boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
		\boldsymbol{\psi}^{(\mathbf{L})}_1 \\
		\boldsymbol{\psi}^{(\mathbf{L})}_2
	\end{bmatrix} \\
	\boldsymbol{\Lambda\mu} &= \boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
	\text{diag}(\mathbf{L}) &= \exp \{ \boldsymbol{\psi}^{(\mathbf{L})}_1 \} \\
	\text{vech}^{\dagger}(\mathbf{L}) &= \boldsymbol{\psi}^{(\mathbf{L})}_2 \\
\end{align}
```
where $\text{vech}^{\dagger}$ is the half-vectorization of a matrix **not including the diagonal**,
"""

# ╔═╡ 44f0dfd2-cde0-4cc5-b8e7-03d36932e335
function f(ψ::AbstractVector)
	L = length(ψ)
	D = Int((-(3/2) + sqrt((9/4) + 2L)) )

	Λμ = ψ[1:D]
	A = diagm(exp.(1e-3 .+ ψ[D+1:2*D]))
	L = A + CompressedLowerTriangular(D, 1, ψ[2*D+1:end])
	L⁻¹ = inv(L)

	Λ = L * L'
	μᵀΛμ = ((L⁻¹' * L⁻¹) * Λμ)' * Λμ
	logdetΛ = 2*sum(ψ[D+1:2*D])

	Λμ, μᵀΛμ, Λ, logdetΛ
end

# ╔═╡ fb6c11f7-e68c-4edf-965e-d36dd19cad2f
function f⁻¹((Λμ, μᵀΛμ, Λ, logdetΛ))
	C = cholesky(Symmetric(Λ))
	vcat(Λμ, log.(diag(C.L)), vech(C.L, 1))
end

# ╔═╡ 93109c2d-a784-408a-9466-c1d43dc8cd48
function f⁻¹(η::AbstractVector)
	L = length(η)
	D = Int((-(3/2) + sqrt((9/4) + 2L)))

	Λμ = η[1:D]
	diagΛ = η[D+1:2*D]
	vechΛ = η[2*D+1:end]
	Λ = diagm(diagΛ) + CompressedSymmetric(D, 1, vechΛ)

	C = cholesky(Λ)
	vcat(Λμ, log.(diag(C.L)), vech(C.L, 1))
end

# ╔═╡ 3d4c8ff8-1621-4710-a49e-7ccda892cd7c
md"""
 $P$ is the dimension of model's parameter space
"""

# ╔═╡ 4a09cd83-a32c-486d-afa0-af460a20d784
P = 2D + (D^2 - D) ÷ 2

# ╔═╡ a0cc6a05-6d8a-4b3f-9886-aa6139382a0c
md"""
 $Q$ is dimension of the subspace
"""

# ╔═╡ c2f69d52-7533-4bea-9457-f8ba4f324ee8
Q = 2

# ╔═╡ 0ae825a4-f21e-4ed9-a08c-2d781d2c9576
pW = JointNormalFixedCov(zeros(Q+1, P))

# ╔═╡ 0fd9bb0e-a8b0-4b55-9af8-8b3356fc53a5
ph = Normal(Q)

# ╔═╡ 37067855-9706-404e-ac02-8e15dc9a84d5
gsm = GSM(f, pW, [NormalModel(ph) for p in collect(phones)[1:end]])

# ╔═╡ 4274cf0e-0cc1-41f4-a0bf-ad46785fea7d
begin
	# Initialization of the variational posterior
	local qθ = (
		JointNormalFixedCov(randn(Q+1, P)),
		[Normal(Q) for model in gsm.models]
	)

	# Total number of observed data points
	local N = sum(size.(Xs, 2))

	Lg = []
	for t in 1:100000
		μ = getstats(gsm, qθ)
		(l, (∇μ,)) = withgradient(μ -> elbo(gsm, Xs, qθ, μ)/N, μ)
		qθ = newposterior(gsm, qθ, ∇μ; lrate=1e-2)

		push!(Lg, l)
	end
	qW, qH = qθ
	Lg
end

# ╔═╡ 35984926-37a0-4de4-a0b7-4eb93ed6bade
begin
	local p = plot(legend=false, xrange=(-3, 3), yrange=(-3,3))
	local k = 1
	local d1, d2 = 2, 3
	for (phone, X, qh) in zip(phones, Xs, qH)
		h = sample(qh)
		W = qW.M
		ψ = W' * vcat(h, 1)
		Λμ, μᵀΛμ, Λ, logdetΛ = f(ψ)
		Σ = inv(Λ)
		μ = Σ*Λμ
		μ = μ[d1:d2]
		Σ = Σ[d1:d2,d1:d2]
		plotnormal!(μ, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
		scatter!(X[d1,:], X[d2,:])
	end
	p
end

# ╔═╡ cfc2db2f-ec6f-4e7c-bda5-aee0a94ef15e
qW

# ╔═╡ 926070be-3e40-456f-95b5-53399d7c2644
μ(qH[1])

# ╔═╡ 70dad4a3-e470-4c57-93c6-840eb29e5317
begin
	local p = plot(legend=false)
	local k = 1
	local d1, d2 = 1, 2
	for (phone, X, qh) in zip(phones, Xs, qH)
		x, diagxxᵀ, vechxxᵀ =  unpack(qh, μ(qh))
		m = x[d1:d2]
		xxᵀ = diagm(diagxxᵀ) + CompressedSymmetric(size(x, 1), 1, vechxxᵀ)
		Σ = (xxᵀ - m*m')[d1:d2,d1:d2]
		annotate!(m[1], m[2], phone)
		plotnormal!(m, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ ebba18b7-1a20-45aa-b843-2baa1b60f468
plot(Lg[10000:10:end], label=false)

# ╔═╡ c6d3f311-7593-40cd-b01c-7b7f6b067a72
function gsm_η(gsm::GSM, qW, h)
	Λμ, μᵀΛμ, Λ, logdetΛ = f(qW.M' * vcat(h, 1))
	vcat(Λμ, diag(Λ), vech(Λ, 1))
end

# ╔═╡ c0a18a14-475f-4d99-8594-f8054e404f14
gsm_η(gsm, qW, zeros(2))

# ╔═╡ 3cf213b2-b9e6-4765-a444-76142e700d93
function gsm_μ(gsm::GSM, η)
	L = length(η)
	D = Int((-(3/2) + sqrt((9/4) + 2L)))

	Λμ = η[1:D]
	diagΛ = η[D+1:2*D]
	vechΛ = η[2*D+1:end]
	Λ = diagm(diagΛ) + CompressedSymmetric(D, 1, vechΛ)

	Σ = inv(Symmetric(Λ))
	μ = Σ * Λμ

	vcat(μ, vec(Σ + μ * μ'))
end

# ╔═╡ a05fb9e3-5ce9-4ccc-9418-178d3d4aa18f
function gsm_ξ(gsm::GSM, qW, η)
	W = qW.M
	W⁺ = inv(W*W') * W
	h = W⁺ * f⁻¹(η)
	h[1:end-1]
end

# ╔═╡ 45d5823c-2044-4815-a6e4-78e02c568f20
function gsm_A(gsm::GSM, η)
	L = length(η)
	D = Int((-(3/2) + sqrt((9/4) + 2L)))

	Λμ = η[1:D]
	diagΛ = η[D+1:2*D]
	vechΛ = η[2*D+1:end]
	Λ = diagm(diagΛ) + CompressedSymmetric(D, 1, vechΛ)
	Σ = inv(Λ)

	-(1/2) * logdet(Λ) + (1/2) * Λμ' * (Σ * Λμ)
end

# ╔═╡ 534dc966-f719-4049-aabc-bed102de1901
gsm_A(gsm, gsm_η(gsm, qW, zeros(2)))

# ╔═╡ b7c4d472-a001-41d8-b2ae-97ee9b095c47
ForwardDiff.hessian(η -> gsm_A(gsm, η), gsm_η(gsm, qW, zeros(2)))

# ╔═╡ 6eaac0ae-2153-4dec-9abc-8114f95b8a67
ForwardDiff.jacobian(h -> gsm_ξ(gsm, qW, h), gsm_η(gsm, qW, zeros(2)))

# ╔═╡ 9952e00f-0b61-4c12-a8f6-4c0384fdd34b
gsm_μ(gsm, gsm_η(gsm, qW, zeros(2)))

# ╔═╡ cee9c7c9-c356-482e-b90e-4c9ec869eb0a
gsm_ξ(gsm, qW, gsm_η(gsm, qW, zeros(2)))

# ╔═╡ 8488be6d-f08e-4d26-a41f-fa3350627d12
function trajectory(gsm, qW, A, B)
	traj = [A]
	while norm(traj[end] - B) > 1e-1
		η = gsm_η(gsm, qW, traj[end])
		H = ForwardDiff.hessian(η -> gsm_A(gsm, η), η)
		#J = ForwardDiff.jacobian(h -> gsm_ξ(gsm, qW, h), η)
		J = ForwardDiff.jacobian(h -> gsm_η(gsm, qW, h), traj[end])
		(∇h,) = Zygote.gradient(h -> -(h - B)' * (h - B), traj[end])

		#∇h = ∇h ./ norm(∇h)
		G = J' * H * J
		v = inv(G) * ∇h
		v = v ./ norm(v)

		push!(traj, traj[end] + 1e-3 * v)
	end
	hcat(traj...)
end

# ╔═╡ 839249a8-dae0-4b29-8717-278b5f495176
begin
	p = plot(xlims=(-3, 3), ylims=(-3, 3), legend=false)
	s, e = -4, 4
	for i = s:0.5:e
		println("i = $(i)")
		t = trajectory(gsm, qW, Float64[s, i], Float64[e, i])
		plot!(t[1, :], t[2, :], linecolor=:black)
	end

	for i = s:0.5:e
		t = trajectory(gsm, qW, [i, s], [i, e])
		plot!(t[1, :], t[2, :], linecolor=:black)
	end

	p
end

# ╔═╡ 6b77d67e-1690-4f60-812e-9687030e3aa2
begin
	local k = 1
	local d1, d2 = 1, 2
	for (phone, X, qh) in zip(phones, Xs, qH)
		x, diagxxᵀ, vechxxᵀ =  unpack(qh, μ(qh))
		m = x[d1:d2]
		xxᵀ = diagm(diagxxᵀ) + CompressedSymmetric(size(x, 1), 1, vechxxᵀ)
		Σ = (xxᵀ - m*m')[d1:d2,d1:d2]

		annotate!(m[1], m[2], arpa2ipa[phone])
		plotnormal!(m, Σ, σ=2, colorfillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ Cell order:
# ╟─e5883656-7aa1-11ec-2648-4f9569fddcf8
# ╠═3395c983-7d14-4163-bebb-096c851cbde8
# ╠═9b39f491-8836-4368-877b-6cf1eea0ff19
# ╠═60f71c97-f4c5-4b4d-9665-4838d676793b
# ╠═52b2c2a9-01c2-4d39-b9d6-4bcd7a11dec1
# ╠═04ed61ae-a1ed-4114-a160-100e86f0769a
# ╠═ff11bde6-a2c8-47f5-a3bf-a4b6f9758747
# ╠═4b071c89-67b6-4964-a4b1-c784f6d99f9c
# ╠═e2aa337d-4cdc-43d2-a8f6-5d0a607a3de0
# ╠═570dc7f7-10dc-46ed-9958-2acc0ed76ad1
# ╠═c74b04e1-b19d-4398-b376-dda5fa88ca86
# ╠═15c82689-be5d-411b-9594-1ba031fb2a0d
# ╟─6c8c5ecc-9dc2-446a-8726-699694677c42
# ╠═2d01cdec-0da2-48f7-9ae9-7cf2298e9d11
# ╠═3814457c-fb1b-4f49-9717-c9d9fcdcc070
# ╠═5fa65b07-fe1c-4119-93fa-85738d74af6c
# ╠═c43a677b-58dd-4266-a93a-c914b7cb3ebc
# ╠═e8148b69-8fba-4d15-8c41-8e445d2cb5a0
# ╠═67ee5a37-0baa-4986-8558-77b720321771
# ╟─d12fcd36-51b4-463e-b01c-807239522d6b
# ╠═44f0dfd2-cde0-4cc5-b8e7-03d36932e335
# ╠═fb6c11f7-e68c-4edf-965e-d36dd19cad2f
# ╠═93109c2d-a784-408a-9466-c1d43dc8cd48
# ╟─3d4c8ff8-1621-4710-a49e-7ccda892cd7c
# ╠═4a09cd83-a32c-486d-afa0-af460a20d784
# ╟─a0cc6a05-6d8a-4b3f-9886-aa6139382a0c
# ╠═c2f69d52-7533-4bea-9457-f8ba4f324ee8
# ╠═0ae825a4-f21e-4ed9-a08c-2d781d2c9576
# ╠═0fd9bb0e-a8b0-4b55-9af8-8b3356fc53a5
# ╠═37067855-9706-404e-ac02-8e15dc9a84d5
# ╠═4274cf0e-0cc1-41f4-a0bf-ad46785fea7d
# ╠═35984926-37a0-4de4-a0b7-4eb93ed6bade
# ╠═cfc2db2f-ec6f-4e7c-bda5-aee0a94ef15e
# ╠═926070be-3e40-456f-95b5-53399d7c2644
# ╠═70dad4a3-e470-4c57-93c6-840eb29e5317
# ╠═ebba18b7-1a20-45aa-b843-2baa1b60f468
# ╠═c6d3f311-7593-40cd-b01c-7b7f6b067a72
# ╠═c0a18a14-475f-4d99-8594-f8054e404f14
# ╠═3cf213b2-b9e6-4765-a444-76142e700d93
# ╠═a05fb9e3-5ce9-4ccc-9418-178d3d4aa18f
# ╠═45d5823c-2044-4815-a6e4-78e02c568f20
# ╠═534dc966-f719-4049-aabc-bed102de1901
# ╠═b7c4d472-a001-41d8-b2ae-97ee9b095c47
# ╠═6eaac0ae-2153-4dec-9abc-8114f95b8a67
# ╠═9952e00f-0b61-4c12-a8f6-4c0384fdd34b
# ╠═cee9c7c9-c356-482e-b90e-4c9ec869eb0a
# ╠═8488be6d-f08e-4d26-a41f-fa3350627d12
# ╠═839249a8-dae0-4b29-8717-278b5f495176
# ╠═6b77d67e-1690-4f60-812e-9687030e3aa2

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
const tips_cache = Dict{Integer, Matrix{Float64}}()

@inline function tips(iso_id::Integer, temp::T)::Float64 where T <: AbstractFloat
    data::Matrix{Float64} = get!(tips_cache, iso_id) do
        path = joinpath(module_path, "meta", "q", @sprintf("q%d.jld2", iso_id))
        if isfile(path) === false
            d = [c_T_ref isotopologue(iso_id)[1, :q_t0]; (c_T_ref + 1.0) isotopologue(iso_id)[1, :q_t0]]
        else
            load(path, "data")
        end
    end
    ind_lo = max(1, min(size(data)[1] - 1, searchsortedfirst(@view(data[:,1]), temp) - 1))
    ind_hi = min(size(data)[1], ind_lo + 1)

    # return linear interpolation
    return (data[ind_hi,2] - data[ind_lo,2]) / (data[ind_hi,1] - data[ind_lo,1]) * (temp - data[ind_lo,1]) + data[ind_lo,2]
end


text:
function Folds.reduce(op::OP, itr, ex::Executor; init = FoldsInit) where {OP}
    result = Transducers.fold(op, itr, ex; init = init)
    result isa FoldsInitOf && return reduce_empty(op, eltype(itr))
    return result
end

Folds.reduce(op::OP, itr; init = FoldsInit, kwargs...) where {OP} =
    Folds.reduce(op, itr, parallel_executor(bottom_foldable(itr); kwargs...); init = init)

Folds.mapreduce(f::F, op::OP, itr; init = FoldsInit, kwargs...) where {F,OP} =
    Folds.mapreduce(
        f,
        op,
        itr,
        parallel_executor(bottom_foldable(itr); kwargs...);
        init = init,
    )

function Folds.mapreduce(
    f::F,
    op::OP,
    itr,
    ex::Executor;
    init = FoldsInit,
) where {F,OP}
    result = Transducers.fold(op, itr |> Map(f), ex; init = init)
    result isa FoldsInitOf && return mapreduce_empty(f, op, eltype(itr))
    return result
end

mapreduce_check_no_kwargs(::NamedTuple{(),Tuple{}}) = nothing
@noinline function mapreduce_check_no_kwargs(kwargs)
    error(
        "`mapreduce(f, op, itrs..., executor)` only accepts `init` as keyword argument" *
        "\ngot:" *
        string(kwargs),
    )
end

function Folds.mapreduce(f, op, itr, itrs...; init = FoldsInit, kwargs...)
    args0, ex0 = de_snoc(itr, itrs...)
    if ex0 isa Executor
        xs = zip(args0...)
        ex = ex0
        mapreduce_check_no_kwargs((; kwargs...))
    else
        xs = zip(itr, itrs...)
        ex = parallel_executor(bottom_foldable(xs); kwargs...)
    end
    result = Transducers.fold(op, xs |> MapSplat(f), ex; init = init)
    result isa FoldsInitOf && return mapreduce_empty(f, op, eltype(xs))
    return result
end

@def_monoid return_nothing
@inline return_nothing(_) = nothing
@inline return_nothing(_, _) = nothing

foreach_check_no_kwargs(::NamedTuple{(),Tuple{}}) = nothing
@noinline function foreach_check_no_kwargs(kwargs)
    error(
        "`foreach(f, itrs..., executor)` accepts no keyword arguments" *
        "\ngot:" *
        string(kwargs),
    )
end

foreach_check_referenceables(_, _) = nothing
@noinline function foreach_check_referenceables(::DistributedEx, args)
    error(
        "`foreach(f, itrs..., ::DistributedEx)` is used with a referenceable. " *
        "Currently, this is not supported as the mutation will not be reflected " *
        "to the original collection."
    )
end

function Folds.foreach(f, itr, itrs...; kwargs...)
    args0, ex0 = de_snoc(itr, itrs...)
    if ex0 isa Executor
        xs = zip(args0...)
        ex = ex0
        foreach_check_no_kwargs((; kwargs...))
        Baselet.Specialized.any(isreferenceable, args0) &&
            foreach_check_referenceables(ex, args0)
    else
        xs = zip(itr, itrs...)
        ex = parallel_executor(bottom_foldable(xs); kwargs...)
    end
    transduce(MapSplat(f), return_nothing, nothing, xs, ex)
    return
end

function Folds.map!(f, dest, array, arrays...; kwargs...)
    Folds.foreach(referenceable(dest), array, arrays...; kwargs...) do y, xs...
        Base.@_inline_meta
        y[] = f(xs...)
    end
    return dest
end

"""
    ReducerFunctionAndFoldable

Reducer function with one or two arguments plus executor.
This function provides the following signature:

    reducer([f,] itr; init, executor_options...)
    reducer([f,] itr, exec; init)
"""
ReducerFunctionAndFoldable

function monoid_for end
default_init_for(reducer) = FoldsInit
default_simd_for(reducer) = Val(true)

"""
    exec(reducer, executor, f, itr, init)
"""
function exec(reducer, ex, f, itr, init)
    return Folds.mapreduce(
        f,
        monoid_for(reducer),
        itr,
        maybe_set_simd(ex, default_simd_for(reducer));
        init = default_to(init, default_init_for(reducer)),
    )
end

(reducer::ReducerFunctionAndFoldable)(itr; kwargs...) = reducer(identity, itr; kwargs...)
(reducer::ReducerFunctionAndFoldable)(f::F, itr; init = unspecified, kwargs...) where {F} =
    exec(reducer, parallel_executor(itr; kwargs...), f, itr, init)

(reducer::ReducerFunctionAndFoldable)(
    f::F,
    itr,
    ex::Executor;
    init = unspecified,
) where {F} = exec(reducer, ex, f, itr, init)

(reducer::ReducerFunctionAndFoldable)(itr, ex::Executor; init = unspecified) =
    exec(reducer, ex, identity, itr, init)

"""
    typeof(Folds.X :: Type{<:ReducerFunctionAndFoldable})
    typeof(Folds.X)

Equivalent to typeof(Folds.X)
"""
reducer_typeof
struct ReducerTypeGetter end
const reducer_typeof = ReducerTypeGetter()
function Base.getproperty(::ReducerTypeGetter, name::Symbol)
    f = getproperty(Folds, name)
    return typeof(f)::ReducerFunctionAndFoldable
end

monoid_for(::typeof(Folds.sum)) = add_sum
monoid_for(::typeof(Folds.prod)) = mul_prod
monoid_for(::typeof(Folds.maximum)) = max
monoid_for(::typeof(Folds.minimum)) = min

asbool(f) = x -> f(x)::Bool

default_init_for(::typeof(Folds.count)) = 0
exec(::typeof(Folds.count), ex, f, itr, init) =
    Folds.sum(ConvertTo(Int) ∘ asbool(f), itr, ex; init = init)
# TODO: Use `ConverTo(typeof(init))` if `init` is a "concrete" initial value.

# TODO: `any` and `all` should be done with "unordered" version
exec(::typeof(Folds.any), ex, f, itr, init) = Transducers.fold(
    right,  # no need to use `|`
    itr |> Map(asbool(f)) |> ReduceIf(identity),
    maybe_set_simd(ex, Val(true));
    init = default_to(init, false),
)

exec(::typeof(Folds.all), ex, f, itr, init) = Transducers.fold(
    right,  # no need to use `&`
    itr |> Map(asbool(f)) |> ReduceIf(!),
    maybe_set_simd(ex, Val(true));
    init = default_to(init, true),
)

exec(::typeof(Folds.findfirst), ex, f, array::AbstractArray, init) = Transducers.fold(
    right,
    pairs(array) |> ReduceIf(((i, x),) -> f(x)) |> Map(first),
    maybe_set_simd(ex, Val(true));
    init = default_to(init, nothing),
)

exec(::typeof(Folds.findlast), ex, f, array::AbstractArray, init) = Transducers.fold(
    right,
    Iterators.reverse(pairs(array)) |> ReduceIf(((i, x),) -> f(x)) |> Map(first),
    maybe_set_simd(ex, Val(true));
    init = default_to(init, nothing),
)

function exec(::typeof(Folds.findall), ex, f, array::AbstractArray, init)
    idxs = Folds.collect(pairs(array) |> Filter(((i, x),) -> f(x)) |> Map(first), ex)
    isempty(idxs) && return keytype(array)[]
    return idxs
end

_minmax((min0, max0), (min1, max1)) = (min(min0, min1), max(max0, max1))

exec(::typeof(Folds.extrema), ex, f, itr, init) = Transducers.fold(
    asmonoid(_minmax),
    itr |> Map(x -> (y = f(x); (y, y))),
    maybe_set_simd(ex, Val(true));
    init = default_to(init, DefaultInit),
)

struct PushUnique{F} <: Function
    f::F
end
PushUnique(::Type{T}) where {T} = PushUnique{Type{T}}(T)

function (f!::PushUnique)((ys, seen), x)
    fx = f!.f(x)
    return fx in seen ? (ys, seen) : (push!!(ys, x), push!!(seen, fx))
end

# TODO: do this with public API of Transducers or make it public
function Transducers.combine(f!::PushUnique, (ys1, seen1), (ys2, seen2))
    seen3 = setdiff!!(seen2, seen1)
    isempty(seen3) && return (ys1, seen1)
    return (append!!(Filter(x -> f!.f(x) in seen3), ys1, ys2), union!!(seen1, seen2))
end
# * Add an option to avoid re-compute `f(x)` in combine?
# * Iterate over `seen3` if `length(seen3) << length(ys2)`?

# Manually create a singleton callable since closure captures the
# types as `DataType` and causes type-instability:
struct InitUnique{X,Y} <: Function end
@inline function (::InitUnique{X,Y})() where {X,Y}
    if isbitstype(Y) || Base.isbitsunion(Y)
        return (X[], Set{Y}())
    else
        return (X[], Empty(Set))
    end
end

function exec(::typeof(Folds.unique), ex, f::F, itr, init) where {F}
    if Base.IteratorEltype(itr) isa Base.HasEltype
        X = eltype(itr)
    else
        X = Union{}
    end
    # Using inference as an optimization. The result of this inference
    # does not affect the returned value:
    Y = infer_eltype(f(x) for x in itr)
    ys, = Transducers.fold(PushUnique(f), itr, ex; init = OnInit(InitUnique{X,Y}()))
    if Base.IteratorEltype(itr) isa Base.HasEltype
        return ys::Vector{X}
    else
        return ys::Vector
    end
end

Folds.issorted(
    itr;
    lt = isless,
    by = identity,
    rev::Union{Bool,Nothing} = nothing,
    order::Base.Ordering = Base.Forward,
    kwargs...,
) = Folds.issorted(itr, Base.ord(lt, by, rev, order); kwargs...)

Folds.issorted(
    itr,
    ex::Executor;
    lt = isless,
    by = identity,
    rev::Union{Bool,Nothing} = nothing,
    order::Base.Ordering = Base.Forward,
    kwargs...,
) = Folds.issorted(itr, Base.ord(lt, by, rev, order), ex)

Folds.issorted(itr, order::Base.Ordering; kwargs...) =
    Folds.issorted(itr, order, parallel_executor(bottom_foldable(itr); kwargs...))

function make_geq(order)
    @inline geq((x1, x2),) = !Base.lt(order, x2, x1)
    return geq
end

#=
Folds.issorted(itr, order::Base.Ordering, ex::Executor) =
    Folds.all(make_geq(order), itr |> Consecutive(Val(2), Val(1)), ex)
=#

# TODO: check if using `zip` for arrays is better
Folds.issorted(itr::AbstractArray, order::Base.Ordering, ex::Executor) = Folds.all(
    make_geq(order),
    zip(
        view(itr, firstindex(itr):lastindex(itr)-1),
        view(itr, firstindex(itr)+1:lastindex(itr)),
    ),
    ex,
)


text:
# top level module for Pumi.  There are individual modules for each Pumi
# library/namespace

module PumiInterface

using MPI

push!(LOAD_PATH, dirname(@__FILE__))

using PumiConfig
using apf
using gmi

export apf, gmi

end # end module


text:
"""
    FEMMDeforLinearNICEModule

Formulation for the small displacement, small strain deformation
model for Nodally-Integrated Continuum Elements (NICE).

The approximation is  originally from Dohrmann et al IJNME 47 (2000).
The formulation was subsequently developed in Krysl, P. and Zhu, B.
Locking-free continuum displacement finite elements with nodal
integration, International Journal for Numerical Methods in Engineering,
76,7,1020-1043,2008.
"""
module FEMMDeforLinearNICEModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, FESetH8, FESetT4, manifdim, nodesperelem, gradN!
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: AbstractFEMMDeforLinear
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D
import FinEtoolsDeforLinear.MatDeforLinearElasticModule: AbstractMatDeforLinearElastic, tangentmoduli!, update!, thermalstrain!
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.FENodeToFEMapModule: FENodeToFEMap
import FinEtoolsDeforLinear.DeforModelRedModule: nstressstrain, nthermstrain, Blmat!
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, loc!, jac!, locjac!, adjugate3!
import FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: stiffness, nzebcloadsstiffness, mass, thermalstrainloads, inspectintegpoints
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtoolsDeforLinear.MatDeforModule: rotstressvec!
import LinearAlgebra: mul!, Transpose, UpperTriangular, eigvals
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, qr, diag, dot, cond, I
import Statistics: mean

"""
    AbstractFEMMDeforLinearNICE <: AbstractFEMMDeforLinear

Abstract FEMM type for Nodally Integrated Continuum Elements (NICE).
"""
abstract type AbstractFEMMDeforLinearNICE <: AbstractFEMMDeforLinear end

mutable struct _NodalBasisFunctionGradients
    gradN::FFltMat
    patchconn::FIntVec
    Vpatch::FFlt
end

"""
    FEMMDeforLinearNICEH8{MR<:AbstractDeforModelRed, S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearNICE

FEMM type for Nodally Integrated Continuum Elements (NICE) based on the eight-node hexahedron.
"""
mutable struct FEMMDeforLinearNICEH8{MR<:AbstractDeforModelRed, S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearNICE
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabfact::FFlt
    nodalbasisfunctiongrad::Vector{_NodalBasisFunctionGradients}
end

function FEMMDeforLinearNICEH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabfact = 0.05
    return FEMMDeforLinearNICEH8(mr, integdomain, mcsys, material, stabfact, _NodalBasisFunctionGradients[])
end

function FEMMDeforLinearNICEH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabfact = 0.05
    return FEMMDeforLinearNICEH8(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabfact, _NodalBasisFunctionGradients[])
end

function FEMMDeforLinearNICEH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M, stabfact::FFlt) where {MR<:AbstractDeforModelRed,  S<:FESetH8, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    return FEMMDeforLinearNICEH8(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabfact, _NodalBasisFunctionGradients[])
end

"""
    FEMMDeforLinearNICET4{MR<:AbstractDeforModelRed, S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearNICE

FEMM type for Nodally Integrated Continuum Elements (NICE) based on the 4-node tetrahedron.
"""
mutable struct FEMMDeforLinearNICET4{MR<:AbstractDeforModelRed, S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic} <: AbstractFEMMDeforLinearNICE
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabfact::FFlt
    nodalbasisfunctiongrad::Vector{_NodalBasisFunctionGradients}
end

function FEMMDeforLinearNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabfact = 0.015
    return FEMMDeforLinearNICET4(mr, integdomain, mcsys, material, stabfact, _NodalBasisFunctionGradients[])
end

function FEMMDeforLinearNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed,  S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabfact = 0.015
    return FEMMDeforLinearNICET4(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabfact, _NodalBasisFunctionGradients[])
end

function FEMMDeforLinearNICET4(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M, stabfact::FFlt) where {MR<:AbstractDeforModelRed,  S<:FESetT4, F<:Function, M<:AbstractMatDeforLinearElastic}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    return FEMMDeforLinearNICET4(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabfact, _NodalBasisFunctionGradients[])
end

function _buffers1(self::AbstractFEMMDeforLinearNICE, geom::NodalField, npts::FInt)
    fes = self.integdomain.fes
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    # Prepare buffers
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    adjJ = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradN = fill(zero(FFlt), nne, mdim);
    xl = fill(zero(FFlt), nne, mdim);
    return loc, J, adjJ, csmatTJ, gradN, xl
end

function _buffers2(self::AbstractFEMMDeforLinearNICE, geom::NodalField, u::NodalField, npts::FInt)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    Jac = fill(zero(FFlt), npts);
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    return loc, J, csmatTJ, Jac, D
end

function patchconn(fes, gl, thisnn)
    # Generate patch connectivity for a given node (thisnn)
    # from the connectivities of the finite elements attached to it.
    return vcat(collect(setdiff(Set([i for j=1:length(gl) for i in fes.conn[gl[j]]]), thisnn)), [thisnn])
end

function computenodalbfungrads(self, geom)
    # # Compute the nodal basis function gradients.
    # # Return the cell array of structures with attributes
    # %		 bfun_gradients{nix}.Nspd= basis function gradient matrix
    # #        bfun_gradients{nix}.Vpatch= nodal patch volume
    # #        bfun_gradients{nix}.patchconn= nodal patch connectivity

    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, adjJ, csmatTJ, gradN, xl = _buffers1(self, geom, npts)

    # Get the inverse map from finite element nodes to geometric cells
    fen2fe = FENodeToFEMap(fes.conn, nnodes(geom));
    # Initialize the nodal gradients, nodal patch, and patch connectivity
    bfungrads = fill(_NodalBasisFunctionGradients(fill(0.0, 0, 0), fill(0, 0), 0.0), nnodes(geom));
    # Now loop over all finite element nodes in the map
    lnmap = fill(0, length(fen2fe.map)); # Local node map: buffer to speed up operations
    for nix = 1:length(fen2fe.map)
        gl = fen2fe.map[nix];
        thisnn = nix; # We are at this node
        if !isempty(gl) # This node has an element patch in this block
            # establish local numbering of all nodes of the patch @ node thisnn
            p = patchconn(fes, gl, thisnn);
            np = length(p);
            lnmap[p] .= 1:np;# now store the local numbers
            c = reshape(geom.values[thisnn, :], 1, ndofs(geom))
            updatecsmat!(self.mcsys, c, J, 0);
            gradNavg = fill(0.0, np, ndofs(geom));# preallocate strain-displacement matrix
            Vpatch = 0.0;
            for k = 1:length(gl)
                i = gl[k]
                kconn = collect(fes.conn[i]);
                pci = findfirst(cx -> cx == thisnn, kconn);# at which node in the element are we with this quadrature point?
                @assert 1 <= pci <= nodesperelem(fes)
                # centered coordinates of nodes in the material coordinate system
                for cn = 1:length(kconn)
                    xl[cn, :] = (reshape(geom.values[kconn[cn], :], 1, ndofs(geom)) - c) * self.mcsys.csmat
                end
                jac!(J, xl, gradNparams[pci])
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                Jac = Jacobianvolume(self.integdomain, J, c, fes.conn[i], Ns[pci]);
                Vpatch += Jac * w[pci];
                sgradN = gradNparams[pci] * adjugate3!(adjJ, J);
                gradNavg[lnmap[kconn],:] += (w[pci] .* sgradN);
            end
            @assert Vpatch != 0
            gradNavg ./= Vpatch;
            bfungrads[nix] = _NodalBasisFunctionGradients(gradNavg, p, Vpatch);
            lnmap[p] .= 0; # Restore the buffer to pristine condition
        end
    end
    self.nodalbasisfunctiongrad = bfungrads
    return self
end

"""
    associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:AbstractFEMMDeforLinearNICE}
    return computenodalbfungrads(self, geom)
end

function Phi3(dim,np,lx,c)
    lx[:,1]=lx[:,1].-c[1];
    lx[:,2]=lx[:,2].-c[2];
    lx[:,3]=lx[:,3].-c[3];
    Phi = fill(0.0, dim*np, 12)
    Phi[1:dim:end-1, 1] .= 1;
    Phi[1:dim:end-1, 2]=lx[:,1]';
    Phi[1:dim:end-1, 3]=lx[:,2]';
    Phi[1:dim:end-1, 4]=lx[:,3]';
    Phi[2:dim:end, 5] .= 1;
    Phi[2:dim:end, 6]=lx[:,1]';
    Phi[2:dim:end, 7]=lx[:,2]';
    Phi[2:dim:end, 8]=lx[:,3]';
    Phi[3:dim:end, 9] .= 1;
    Phi[3:dim:end, 10]=lx[:,1]';
    Phi[3:dim:end, 11]=lx[:,2]';
    Phi[3:dim:end, 12]=lx[:,3]';
    return Phi
end

"""
    stiffness(self::AbstractFEMMDeforLinearNICE, assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::AbstractFEMMDeforLinearNICE, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, csmatTJ, Jac, D = _buffers2(self, geom, u, npts)
    tangentmoduli!(self.material, D, 0.0, 0.0, loc, 0)
    Dmod = sort(eigvals(D))
    stabDmod = mean(Dmod[1:2], dims=1)
    elmatsizeguess = 4*nodesperelem(fes)*ndofs(u)
    startassembly!(assembler, elmatsizeguess, elmatsizeguess, nnodes(u), u.nfreedofs, u.nfreedofs);
    for nix = 1:length(self.nodalbasisfunctiongrad)
        gradN = self.nodalbasisfunctiongrad[nix].gradN
        patchconn = self.nodalbasisfunctiongrad[nix].patchconn
        Vpatch = self.nodalbasisfunctiongrad[nix].Vpatch
        c = reshape(geom.values[nix, :], 1, ndofs(geom))
        updatecsmat!(self.mcsys, c, J, 0);
        nd = length(patchconn) * ndofs(u)
        Bnodal = fill(0.0, size(D, 1), nd)
        Blmat!(self.mr, Bnodal, Ns[1], gradN, c, self.mcsys.csmat);
        elmat = fill(0.0, nd, nd) # Can we SPEED it UP?
        DB = fill(0.0, size(D, 1), nd)
        add_btdb_ut_only!(elmat, Bnodal, Vpatch, D, DB)
        complete_lt!(elmat)
        if (self.stabfact > 0)
            Phi = Phi3(ndofs(u), length(patchconn), geom.values[patchconn,:], c);
            A1 = Phi* ((Phi'*Phi)\Phi');
            elmat += (self.stabfact*stabDmod).*(I-A1);
        end
        dofnums = fill(0, nd)
        gatherdofnums!(u, dofnums, patchconn); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

"""
    stiffness(self::AbstractFEMMDeforLinearNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::AbstractFEMMDeforLinearNICE, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom, u);
end

"""
nzebcloadsstiffness(self::AbstractFEMMDeforLinear,  assembler::A,
  geom::NodalField{FFlt},
  u::NodalField{T}) where {A<:AbstractSysvecAssembler, T<:Number}

Compute load vector for nonzero EBC for fixed displacement.
"""
function nzebcloadsstiffness(self::AbstractFEMMDeforLinearNICE,  assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:AbstractSysvecAssembler, T<:Number}
    error("Not implemented yet")
end

function nzebcloadsstiffness(self::AbstractFEMMDeforLinearNICE, geom::NodalField{FFlt}, u::NodalField{T}) where {T<:Number}
    assembler = SysvecAssembler()
    return  nzebcloadsstiffness(self, assembler, geom, u);
end

end


text:
module LibCEED

using StaticArrays, UnsafeArrays, Requires
using Base: RefValue

# import low-level C interface
include("C.jl")
import .C

export @interior_qf,
    @witharray,
    @witharray_read,
    Abscissa,
    AbscissaAndWeights,
    Basis,
    BasisCollocated,
    COPY_VALUES,
    Ceed,
    CeedDim,
    CeedInt,
    CeedScalar,
    CeedVector,
    CeedVectorActive,
    CeedVectorNone,
    Context,
    CopyMode,
    EVAL_CURL,
    EVAL_DIV,
    EVAL_GRAD,
    EVAL_INTERP,
    EVAL_NONE,
    EVAL_WEIGHT,
    ElemRestriction,
    ElemRestrictionNone,
    EvalMode,
    GAUSS,
    GAUSS_LOBATTO,
    HEX,
    LINE,
    MEM_DEVICE,
    MEM_HOST,
    MemType,
    NORM_1,
    NORM_2,
    NORM_MAX,
    NOTRANSPOSE,
    NormType,
    OWN_POINTER,
    Operator,
    PRISM,
    PYRAMID,
    QFunction,
    QFunctionNone,
    QUAD,
    QuadMode,
    RequestImmediate,
    RequestOrdered,
    STRIDES_BACKEND,
    TET,
    TRANSPOSE,
    TRIANGLE,
    Topology,
    TransposeMode,
    USE_POINTER,
    UserQFunction,
    add_input!,
    add_output!,
    apply!,
    apply_add!,
    apply,
    assemble,
    assemble_add_diagonal!,
    assemble_diagonal!,
    axpy!,
    ceedversion,
    create_composite_operator,
    create_elem_restriction,
    create_elem_restriction_strided,
    create_evector,
    create_h1_basis,
    create_identity_qfunction,
    create_interior_qfunction,
    create_lvector,
    create_tensor_h1_basis,
    create_tensor_h1_lagrange_basis,
    create_vectors,
    det,
    extract_array,
    extract_context,
    gauss_quadrature,
    get_preferred_memtype,
    getcompstride,
    getnumelements,
    getelementsize,
    getlvectorsize,
    getmultiplicity!,
    getmultiplicity,
    getdimension,
    getgrad,
    getgrad1d,
    getinterp,
    getinterp1d,
    getnumcomponents,
    getnumnodes,
    getnumnodes1d,
    getnumqpts,
    getnumqpts1d,
    getqref,
    getqweights,
    getresource,
    gettopology,
    getvoigt!,
    getvoigt,
    iscuda,
    isdeterministic,
    isrelease,
    lobatto_quadrature,
    norm,
    pointwisemult!,
    reciprocal!,
    scale!,
    set_context!,
    set_cufunction!,
    set_data!,
    set_field!,
    setarray!,
    setvalue!,
    setvoigt!,
    setvoigt,
    syncarray!,
    takearray!,
    witharray,
    witharray_read

include("Globals.jl")
include("Ceed.jl")
include("CeedVector.jl")
include("Basis.jl")
include("ElemRestriction.jl")
include("Quadrature.jl")
include("Context.jl")
include("UserQFunction.jl")
include("QFunction.jl")
include("Request.jl")
include("Operator.jl")
include("Misc.jl")

cuda_is_loaded = false

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("Cuda.jl")
    set_globals()
end

end # module

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
module RobustnessEvaluation

using Flux
using Statistics
using NeuralVerification
using Plots
using Measures
using Random
using ProgressMeter
using NeuralVerification: compute_output
using NeuralVerification: init_vars, BoundedMixedIntegerLP, encode_network!, _ẑᵢ₊₁, TOL
using GLPK
using JuMP
using LazySets
using StatsBase
import JuMP.MOI.OPTIMAL, JuMP.MOI.INFEASIBLE, JuMP.MOI.INFEASIBLE_OR_UNBOUNDED

struct RobustnessSummary
    samples
    errors_nominal
    mae_nominal
    errors_adversarial_true
    mae_adversarial_true
    errors_adversarial_pred
    mae_adversarial_pred
    saliency_maps
    ig_maps
end

export nominal_errors
include("nominal.jl")

export adversarial_errors
include("adversarial.jl")

export saliency_map, integrated_gradients, get_mask
include("explanations.jl")

export plot_summary, plot_nominal_errors, plot_adversarial_errors, plot_adversarial_perturbations, plot_explanations
include("plotting.jl")

export robustness_summary

function robustness_summary(network, X, Y; N_samples = 100, coeffs = [-0.74, -0.44],
    δ = 0.02, solver = GLPK.Optimizer, use_gurobi = false, PGD = false,
    step_factor = 0.1, iterations = 600, baseline = zeros(128))
    """ Runs all functions to get robustness summary for a given network and dataset
        Args:
            - network: neural network to analyze (flux model)
            - X: training examples (each column is an example)
            - Y: labels (each column is a label)
            ---------------
            - N_samples: number of samples to calculate on (actual indices will be chosen randomly)
            - coeffs: coefficients to get output (linear combination of network outputs)
            - δ: radius for hyperrectangle to check error within
            - solver: solver for mixed integer programs (ignored if PGD true)
            - use_gurobi: whether we are using gurobi (still need to put in solver)
            - PGD: whether or not to use PGD (if not, solve exactly with MIP Verify)
            - step_factor: step factor for PGD (ignored if PGD false)
            - iterations: number of steps for PGD (ignored if PGD false)
            - baseline: IG baseline (should be same size as input)
    """
    N = size(Y, 2)
    nnv_network = NeuralVerification.network(network)

    # Get samples to work with
    samples = randperm(N)[1:N_samples]
    X_samples = X[:, samples]
    Y_samples = Y[:, samples]

    # Nominal
    println("Calculating nominal errors...")
    errors_nominal = nominal_errors(nnv_network, X_samples, Y_samples, coeffs = coeffs)
    mae_nominal = mean(errors_nominal)

    # Adversarial
    println("Calculating adversarial errors...")
    errors_adversarial_true, errors_adversarial_pred = adversarial_errors(nnv_network,
        X_samples, Y_samples, coeffs = coeffs, δ = δ, solver = solver,
        use_gurobi = use_gurobi, PGD = PGD, step_factor = step_factor,
        iterations = iterations)
    mae_adversarial_true = mean(errors_adversarial_true)
    mae_adversarial_pred = mean(errors_adversarial_pred)

    # Explanations
    println("Getting explanations...")
    saliency_maps = zeros(size(X_samples))
    ig_maps = zeros(size(X_samples))
    for i = 1:N_samples
        saliency_maps[:, i] = saliency_map(network, X[:, i], coeffs = coeffs)
        ig_maps[:, i] = integrated_gradients(network, X[:, i], baseline, coeffs = coeffs)
    end

    return RobustnessSummary(samples, errors_nominal, mae_nominal, errors_adversarial_true,
        mae_adversarial_true, errors_adversarial_pred, mae_adversarial_pred,
        saliency_maps, ig_maps)
end

end # module


text:
#module B
#export bfunc

#include("a.jl")
#using ..A

function bfunc(::Astruct)
end

#end


text:
# Also see related functions in users.jl and personalization.jl

#### GET ####

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-current-user-follows
"""
    follow_check(item_type, ids)

**Summary**: Check to see if the current user is following one or more artists or other Spotify users.

# Arguments
- `item_type::String` _Required_: The ID type, either `artist` or `user`.\n
- `ids::String` _Required_: A comma separated list of the artist or user Spotify IDs to check. Maximum 50.\n

# Example
```julia-repl
julia> Spotify.follow_check("artist", "7fxBPUc2bTUgl7GLuqjajk")[1]
[ Info: We try the request without checking if current grant includes scope user-follow-read.
1-element JSON3.Array{Bool, Base.CodeUnits{UInt8, String}, Vector{UInt64}}:
 1
```
"""
function follow_check(item_type::String, ids::String)

    return spotify_request("me/following/contains?type=$item_type&ids=$ids"; scope = "user-follow-read")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-if-user-follows-playlist
"""
    follow_check_playlist(playlist_id::String, user_id::String)

**Summary**: Check to see if one or more Spotify users are following a specified playlist_id.\n

# Arguments
- `playlist_id::String` _Required_: The Spotify ID of the playlist.\n
- `user_id::String` _Required_: A comma separated list of the user Spotify IDS to check. Maximum 5.\n

# Example
```julia-repl
julia> Spotify.follow_check_playlist("37i9dQZF1DZ06evO3Khq6I", user_id)[1]
1-element JSON3.Array{Bool, Base.CodeUnits{UInt8, String}, Vector{UInt64}}:
 0
```
"""
function follow_check_playlist(playlist_id, user_id)

    return spotify_request("playlists/$playlist_id/followers/contains?ids=$user_id"; scope = "playlist-read-private")

end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-followed
"""
    follow_artists(item_type::String="artist", limit::Int64=20)

**Summary**: Get the current user's followed artists.

# Arguments
- `item_type::String` _Required_: The ID type. Currently only `artist` is supported. Default `artist`.\n
- `limit::Int64` _Optional_: The maximum number of items to return. Default 20, Minimum 1, Maximum 50. \n

# Example
```julia-repl
julia> Spotify.follow_artists()[1]["artists"]
[ Info: We try the request without checking if current grant includes scope user-follow-modify.
JSON3.Object{Base.CodeUnits{UInt8, String}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}} with 6 entries:
  :items   => JSON3.Object[{…
  :next    => nothing
  :total   => 2
  :cursors => {…
  :limit   => 20
  :href    => "https://api.spotify.com/v1/me/following?type=artist&limit=20"
```
"""
function follow_artists(item_type::String="artist", limit::Int64=20)

    return spotify_request("me/following?type=$item_type&limit=$limit"; scope = "user-follow-modify")

end

#### PUT ####

## https://developer.spotify.com/documentation/web-api/reference/follow/follow-artists-users/
@doc """
# Follow Artists or Users
**Summary**: Add the current user as a follower of one or more artists or other Spotify users.\n

`type` _Required_: The ID type: either `artist` or `user`. \n
`ids` _Required_: A comma-separated list of the artists or users Spotify IDs. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/follow/follow-artists-users/)
""" ->
function follow_artists_users(type, ids)
    return spotify_request("me/following?type=$type&ids=$ids", method="PUT"; scope = "user-follow-modify")
end

## https://developer.spotify.com/documentation/web-api/reference/follow/follow-playlist/
@doc """
# Follow a Playlist
**Summary**: Add the currend user as a follower of a playlist. \n

`playlist_id` _Required_: The Spotify ID of the playlist. Any playlist can be followed regardless of it's private/public status, as long as the ID is known.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/follow/follow-playlist/)
""" ->
function follow_playlist(playlist_id)
    return spotify_request("playlists/$playlist_id/followers", method="PUT"; scope = "user-follow-modify")
end

#### DELETE ####

## https://developer.spotify.com/documentation/web-api/reference/follow/unfollow-artists-users/
@doc """
# Unfollow Artists or Users
**Summary**: Remove the current user as a follower of one or more artists or other Spotify users.\n

`type` _Required_: The ID type: either `artist` or `user`. \n
`ids` _Required_: A comma-separated list of the artists or users Spotify IDs. Maximum 50.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/follow/unfollow-artists-users/)
""" ->
function unfollow_artists_users(type, ids)
    return spotify_request("me/following?type=$type&ids=$ids", method="DELETE"; scope = "user-follow-modify")
end

## https://developer.spotify.com/documentation/web-api/reference/follow/unfollow-playlist/
@doc """
# Unfollow a Playlist
**Summary**: Remove the current user as a follower of a playlist.\n

`playlist_id` _Required_: The Spotify ID of the playlist. Any playlist can be followed regardless of it's private/public status, as long as the ID is known.\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/follow/unfollow-playlist/)
""" ->
function unfollow_playlist(playlist_id)
    return spotify_request("playlists/$playlist_id/followers", method="DELETE"; scope = "playlist-modify-private")
end


text:
## Class shape transformation method
#==========================================================================================#

# Basic shape function
function shape_function(x, basis_func, coeffs, coeff_LE = 0)
    n     = length(coeffs)
    terms = basis_func.(x, n - 1, 0:n-1)
    dot(coeffs, terms) + coeff_LE * (x^0.5) * (1 - x)^(n - 0.5)
end

# Computing coordinates
CST_coordinates(class_func, basis_func, x, alphas, dz, coeff_LE, args...) = class_func(x) * shape_function(x, basis_func, alphas, coeff_LE) + x * dz

## Bernstein basis
#==========================================================================================#

bernstein_class(x, N1, N2) = x^N1 * (1 - x)^N2
bernstein_basis(x, n, k)   = binomial(n, k) * bernstein_class(x, k, n - k)

"""
    kulfan_CST(alpha_u, alpha_l,
               (Δz_u, Δz_l) = (0., 0.),
               (LE_u, LE_l) = (0., 0.),
               n            = 40)

Define a cosine-spaced foil with ``2n`` points using the Class Shape Transformation method on a Bernstein polynomial basis for the upper and lower coordinates.

The foil is defined by arrays of coefficients ``(α_u,~ α_l)`` for the upper and lower surfaces (not necessarily of the same lengths), trailing-edge displacement values ``(Δz_u,~ Δz_l)``, and coefficients for leading edge modifications on the upper and lower surfaces at the nose.
"""
function kulfan_CST(alpha_u, alpha_l, (dz_u, dz_l) = (0., 0.), (LE_u, LE_l) = (0., 0.), n :: Integer = 40, N1 = 0.5, N2 = 1.)
    # Cosine spacing for airfoil of unit chord length
    xs = cosine_spacing(0.5, 1, n)

    # λ-function for Bernstein polynomials
    bernie(x, alphas, dz, LE) = CST_coordinates(y -> bernstein_class(y, N1, N2), bernstein_basis, x, alphas, dz, LE)

    # Upper and lower surface generation
    upper_surf = [ bernie(x, alpha_u, dz_u, LE_u) for x ∈ xs ]
    lower_surf = [ bernie(x, alpha_l, dz_l, LE_l) for x ∈ xs ]

    # Counter-clockwise ordering
    @views Foil([ xs[end:-1:2] upper_surf[end:-1:2] ;
                  xs           lower_surf           ], "Kulfan CST")
end

"""
    camber_CST(α_c, α_t,
               (Δz_u, Δz_l) :: NTuple{2, Real},
               coeff_LE = 0.,
               n :: Integer = 40)

Define a cosine-spaced foil with ``2n`` points using the Class Shape Transformation method on a Bernstein polynomial basis for the camber and thickness coordinates.

The foil is defined by arrays of coefficients ``(α_c,~ α_t)`` for the upper and lower surfaces, trailing-edge spacing values ``(Δz_u,~Δz_l)``, and a coefficient for the leading edge modifications at the nose.
"""
function camber_CST(α_cam, α_thicc, dz_thicc = 0., coeff_LE = 0, n :: Integer = 40, N1 = 0.5, N2 = 1.)
    # Cosine spacing for airfoil of unit chord length
    xs = cosine_spacing(0.5, 1, n)

    # λ-function for Bernstein polynomials
    bernie(x, αs, dz = 0.) = CST_coordinates(y -> bernstein_class(y, N1, N2), bernstein_basis, x, αs, dz, coeff_LE)

    # Upper and lower surface generation
    cam   = [ bernie(x, α_cam) for x ∈ xs ]
    thicc = [ bernie(x, α_thicc, dz_thicc) for x ∈ xs ]

    Foil(camber_thickness_to_coordinates(xs, cam, thicc), "Camber-Thickness CST")
end

"""
    coordinates_to_CST(coords, num_dvs)

Convert coordinates to a specified number of CST variables by performing a least-squares solution.
"""
function coordinates_to_CST(coords, num_dvs)
    xs       = @views coords[:,1]
    S_matrix = reduce(hcat, @. bernstein_class(xs, 0.5, 1.0) * bernstein_basis(xs, num_dvs - 1, i) for i in 0:num_dvs - 1)
    alphas   = @views S_matrix \ coords[:,2]
end

"""
    camber_thickness_to_CST(coords, num_dvs)

Convert camber-thickness coordinates to a specified number of CST variables by performing a least-squares solution.
"""
function camber_thickness_to_CST(coords, num_dvs)
    xs, camber, thickness = (columns ∘ coordinates_to_camber_thickness)(coords)

    alpha_cam   = coordinates_to_CST([ xs camber ], num_dvs)
    alpha_thick = coordinates_to_CST([ xs thickness ], num_dvs)

    alpha_cam, alpha_thick
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module ScalingModel
using DifferentialEquations

function micro_1(du, u, parms, time)
    # PARAMETER DEFS
    # β transmition rate
    # r net population growth rate
    # μ hosts' natural mortality rate
    # Κ population size
    # α disease induced mortality rate

    β, r, μ, K, α = parms
    S, I = u
    dS = r*(1-S/K)*S - β*S*I
    dI = β*S*I-(μ+α)*I
    du = [dS,dI]
end

# +
# PARAMETER DEFS
# w and m are used to define the other parameters allometrically

function main()
    w = 1;
    m = 10;
    β = 0.0247*m*w^0.44;
    r = 0.6*w^-0.27;
    μ = 0.4*w^-0.26;
    K = 16.2*w^-0.7;
    α = (m-1)*μ;
    # -

    parms = [β,r,μ,K,α];
    init = [K,1.];
    tspan = (0.0,10.0);

    sir_prob = ODEProblem(micro_1,init,tspan,parms)

    sir_sol = solve(sir_prob);
    return sir_sol
end

# using Plots

# plot(sir_sol,xlabel="Time",ylabel="Number")

# m = [5,10,20,40]
# ws = 10 .^collect(range(-3,length = 601,3))
# βs = zeros(601,4)
# for i = 1:4
#     βs[:,i] = 0.0247*m[i]*ws.^0.44
# end
# plot(ws,βs,xlabel="Weight",ylabel="\\beta_min", xscale=:log10,yscale=:log10, label=["m = 5" "m = 10" "m = 20" "m = 40"],lw=3)

end


text:
using GZip

import Base.convert, Base.show

export STRAND_FORWARD, STRAND_REVERSE, STRAND_UNDEFINED

###################
# Strand definition
###################

immutable type Strand end
const STRAND_FORWARD = Strand()
const STRAND_REVERSE = Strand()
const STRAND_UNDEFINED = Strand()

function Base.convert(::Type{Strand}, c :: Char)
    if c == '+'
        STRAND_FORWARD
    elseif c == '-'
        STRAND_REVERSE
    else
        STRAND_UNDEFINED
    end
end

function Base.convert(::Type{Char}, s :: Strand)
    if s == STRAND_FORWARD
        '+'
    elseif s == STRAND_REVERSE
        '-'
    else
        '.'
    end
end

###########################################
# Generic API for individual Region objects
###########################################

abstract Region

function contig(gr :: Region)
    gr.contig
end

function range(gr :: Region)
    gr.range
end

function range_start (gr :: Region)
    range(gr)[1]
end

function range_end (gr :: Region)
    range(gr)[end]
end

function name(gr :: Region)
    gr.name
end

function score(gr :: Region)
    gr.score
end

function strand(gr :: Region)
    gr.strand
end

function isless(r1 :: Region, r2 :: Region)
    contig(r1) < contig(r2) \
        | (range_start(r1) < range_start(r2))
        | (range_end(r1) < range_end(r2))
end

####################
# Generic Region I/O
####################

type RegionStream{T <: Region}
    handle :: IO
end

function write{T}(rio :: RegionStream{T},  r :: Region)
    write(rio.handle, convert(T, r))
end

function advance{T}(it :: RegionStream{T})
    try
        read(it.handle, T)
    catch EOFError
        nothing
    end
end

Base.start{T}(it :: RegionStream{T}) = advance(it)
Base.next{T}(it :: RegionStream{T}, nxt) = (nxt, advance (it))
Base.done (it :: RegionStream, nxt) = ( nxt == nothing ) || eof (it.handle)

#########################################
# Mapping file extensions to Region types
#########################################

EXTENSIONS = Dict{String, Type}()

function register_file_extension!{T <: Region}(::Type{T}, ext :: String)
    EXTENSIONS[ext] = T
end

function infer_region_type_from_path(path :: String)
    for (ext,t) in EXTENSIONS
        if endswith(lowercase(path), ext)
            return t
        end
    end
    EXTENSIONS[".bed"]
end

#######################################
# Generic region I/O ( High-level API )
#######################################

export gropen, grload, grsave

function gropen(path :: String, mode :: String = "r")
    handle = endswith(path, ".gz") ? gzopen(path, mode) : open(path, mode)
    t = infer_region_type_from_path(path)
    RegionStream{t}(handle)
end

function grload(path :: String)
    collect(gropen(path))
end

function grsave(path :: String, regions)
    h = gropen(path, "w")
    for r in regions
        write(h.handle, r)
    end
    close(h.handle)
end


text:
using MINLPTests, JuMP, Ipopt, Juniper, Test

const OPTIMIZER = MINLPTests.JuMP.optimizer_with_attributes(
    Juniper.Optimizer, "nlp_solver"=>JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "sb" =>"yes"), "log_level" => []
)

@testset "MINLPTests" begin
    MINLPTests.test_nlp_mi(OPTIMIZER, exclude = [
        "005_011",  # Uses the function `\`
        "003_013",  # Bug in Juniper - handling of expression graph?
        "006_010"   # Bug in Juniper - handling of user-defined functions.
    ], objective_tol = 1e-5, primal_tol = 1e-5, dual_tol = NaN)
end


text:
using DRYlib
using Base.Test

# write your own tests here
@test 6 * 7 == 42


text:
using MarketEquilibria
using DelimitedFiles
using Random
using Plots
using Plots.PlotMeasures

form = [:linear, :CES, :CES, :CES, :cobb_douglas,
        :CES, :CES, :CES, :CES, :CES, :leontief]
# Rho value ignored for anything other than :CES
rho = [000, .75, .5, .25, 000, -0.5,
       -1., -2, -5, -10, 000]
rho_label = ["Linear", "0.75", "0.5", "0.25",
             "Cobb–\nDouglas", "-0.5", "-1.0", "-2.0",
             "-5.0", "-10.0", "Leontief"]
samp = size(rho)[1]

m, n = 4, 6

# endowments =  1 .+ randexp(m)
# A = rand(m, n)
# A ./= sum(A, dims=2)
# supplies = 1 .+ rand(n)
#
# writedlm("examples/Fisher/in/endowments.dat", endowments)
# writedlm("examples/Fisher/in/A.dat", A)
# writedlm("examples/Fisher/in/supplies.dat", supplies)

endowments = reshape(readdlm("examples/Fisher/in/endowments.dat"), :)
A          = readdlm("examples/Fisher/in/A.dat")
supplies   = reshape(readdlm("examples/Fisher/in/supplies.dat"), :)

prices = zeros(n, samp)
X      = zeros(m, n, samp)

for i in 1:samp
    global X
    global prices
    X[:, :, i], prices[:, i] = homogeneousfisher(endowments, A, supplies, form[i], rho[i])
    prices[:, i] .*= -1.
    signal = isapprox(endowments, X[:, :, i] * prices[:, i], atol=1e-3)
    @assert signal "Something went wrong; rho = $(rho[i])"
end

writedlm("examples/Fisher/out/X.dat", X)
writedlm("examples/Fisher/out/prices.dat", prices)

#A_bar = sum(A, dims=1) / m
A_bar = round.(A' * endowments, digits=4)

p = plot(1:samp,
         prices',
         xticks=(1:samp, rho_label),
         xlabel="ρ",
         ylabel="price",
         linestyle=[:dash :dot :solid],
         title="Equilibrium price for homogeneous CES utility, $m players, $n resources",
         titlefontsize=12,
         legend=:topleft,
         label=reshape(["π$i, ā = $(A_bar[i])" for i in 1:n], 1, :))

savefig(p, "examples/Fisher/EqPlot.png")
savefig(p, "examples/Fisher/EqPlot.pdf")

X1 = X[1, :, :] ./ sum(X[1, :, :], dims = 1)
X1A = round.(A[1, :], digits=4)

q = plot(1:samp,
         cumsum(X1[end:-1:1, :], dims=1)'[:,end:-1:1],
         fill=0,
         xlim=(1, samp),
         right_margin=5mm,
         ylim=(0, 1),
         xticks=(1:samp, rho_label),
         xlabel="ρ",
         ylabel="Player 1’s optimal allocation",
         ls=:solid,
         lc=:black, lw=.5,
         title="Equilibrium allocation for homogeneous CES utility",
         titlefontsize=12,
         legend=:top,
         label=reshape(["x$i, a = $(X1A[i])" for i in 1:n], 1, :))

savefig(q, "examples/Fisher/XPlot.png")
savefig(q, "examples/Fisher/XPlot.pdf")

max_stars_repo_name: maxkapur/MarketEquilibria

max_issues_repo_name: maxkapur/MarketEquilibria

max_forks_repo_name: maxkapur/MarketEquilibria


text:
@testset "Utils" begin
    @test isnothing(posinf(Int))
    @test posinf(Float64) == Inf
    @test posinf(Rational{Int}) == Inf

    @test isnothing(neginf(Int))
    @test neginf(Float64) == -Inf
    @test neginf(Rational{Int}) == -Inf

    @test !hasposinf(Int)
    @test hasposinf(Float64)
    @test hasposinf(Rational{Int})

    @test !hasneginf(Int)
    @test hasneginf(Float64)
    @test hasneginf(Rational{Int})

    @test !hasinf(Int)
    @test hasinf(Float64)
    @test hasinf(Rational{Int})

    @test !isposinf(1)
    @test isposinf(Inf)
    @test !isposinf(-Inf)

    @test !isneginf(1)
    @test !isneginf(Inf)
    @test isneginf(-Inf)
end