text:
"""
    build_opf_powermodels((sp::JuMP.Model, data::Dict, params::Dict)

Network grid model builder.

Required parameters are:
-   sp is a JuMP model.
-   data is a dict with information of the subproblem.
-   param is a dict containing solution parameters.
"""
function build_opf_powermodels(sp::JuMP.Model, data::Dict, params::Dict)
    return PowerModels.build_model(data["powersystem"], params["model_constructor_grid"],
            params["post_method"], jump_model=sp, setting = params["setting"])
end

"""
    build_graph(data::Dict, params::Dict)

Graph description.
"""
function build_graph(params::Dict)
    # graph definition
    graph = SDDP.LinearGraph(0)
    SDDP.add_node(graph, 1)
    SDDP.add_edge(graph, 0 => 1, 1.0)
    for t in 2:params["stages"]
        SDDP.add_node(graph, t)
        SDDP.add_edge(graph, t-1 => t, params["discount_factor"])
    end
    return graph
end

ext: jl


text:
# # Global Climatology
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/global_ocean_circulation.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/global_ocean_circulation.ipynb)
#
# Advect particles with climatological monthly mean flow at selected depth level
# (e.g. `k=10` for 95 m) from a global ocean state estimate ([ECCO v4 r2](https://eccov4.readthedocs.io/en/latest/) ; see also <https://ecco-group.org>)
# which is here repeated for `ny` years. For additional documentation e.g. see :
# [1](https://JuliaClimate.github.io/MeshArrays.jl/dev/),
# [2](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),
# [3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),
# [4](https://en.wikipedia.org/wiki/Displacement_(vector))
#
# [![simulated particle movie (5m)](https://user-images.githubusercontent.com/20276764/84766999-b801ad80-af9f-11ea-922a-610ad8a257dc.png)](https://youtu.be/W5DNqJG9jt0)

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ## 1. Get Software & Iput Files
#
# - packages + helper functions
# - grid and velocity files

using IndividualDisplacements, OceanStateEstimation, DataFrames, Statistics, CSV

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/helper_functions.jl"))

OceanStateEstimation.get_ecco_velocity_if_needed();

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ## 2. Set Up Parameters & Inputs
#
# - select vertical level (k=1 by default; k=0 for 3D) & duration in years (ny=2 by default)
# - read grid variables
# - return FlowFields (𝑃) and ancillary variables etc (𝐷)
# - read & normalize velocities (𝐷.🔄)

𝑃,𝐷=global_ocean_circulation(k=1,ny=2);

𝐷.🔄(𝑃,𝐷,0.0)

fieldnames(typeof(𝑃))

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3. Main Computation Loop
#
# ### 3.1 Initialize Individuals & Solution
#
# - initial particle positions randomly over Global Ocean

np=100

#xy = init_global_randn(np,𝐷)
#df=DataFrame(x=xy[1,:],y=xy[2,:],f=xy[3,:])

p=dirname(pathof(IndividualDisplacements))
fil=joinpath(p,"../examples/worldwide/global_ocean_circulation.csv")
df=DataFrame(CSV.File(fil))

𝐼=Individuals(𝑃,df.x[1:np],df.y[1:np],df.f[1:np])
fieldnames(typeof(𝐼))

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# - initial integration from time 0 to 0.5 month

📌ini=deepcopy(𝐼.📌)
𝑇=(0.0,𝐼.𝑃.𝑇[2])
∫!(𝐼,𝑇)

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### 3.2 Iteration function example
#
# In addition, `step!` is defined to provide additional flexibility around `∫!` :
#
# - `𝐷.🔄(𝐼.𝑃,t_ϵ)` resets the velocity input streams to bracket t_ϵ=𝐼.𝑃.𝑇[2]+eps(𝐼.𝑃.𝑇[2])
# - `reset_📌!(𝐼)` randomly selects a fraction (`𝐷.frac`) of the particles and resets their positions before each integration period. This tends to maintain homogeneous coverage of the Global Ocean by particles.
# - `∫!(𝐼)` then solves for the individual trajectories over one month, with updated velocity fields (𝐼.𝑃.u0 etc), and adds diagnostics to the DataFrame used to record variables along the trajectory (𝐼.🔴).

function step!(𝐼::Individuals)
    t_ϵ=𝐼.𝑃.𝑇[2]+eps(𝐼.𝑃.𝑇[2])
    𝐷.🔄(𝐼.𝑃,𝐷,t_ϵ)
    reset_📌!(𝐼,𝐷.frac,📌ini)
    ∫!(𝐼)
end

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3.3 Iterate For `ny*12` Months
#

[step!(𝐼) for y=1:1, m=1:1]

add_lonlat!(𝐼.🔴,𝐷.XC,𝐷.YC);

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3.4 Compute summary statistics
#
# See [DataFrames.jl](https://juliadata.github.io/DataFrames.jl/latest/) documentation for detail and additinal functionalities.

gdf = groupby(𝐼.🔴, :ID)
sgdf= combine(gdf,nrow,:lat => mean)
sgdf[rand(1:size(sgdf,1),4),:]

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 4. Plot trajectories / individual positions
#
# ```
# using Plots
# p=plot(;xlims=(-180,180),ylims=(-90,90),legend=:none)
# p!(x,y)=scatter!(p,x,y,markersize=1.1,markerstrokewidth=0)
# [p!(gdf[i].lon,gdf[i].lat) for i in rand(collect(1:length(gdf)),10)]
# display(p)
# ```

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# Or select a background map (e.g. `lon`, `lat`, and `DL=log10(bottom depth)`)
# and a recipe to superimpose initial and final locations. Try:
#
# ```
# include(joinpath(dirname(pathof(IndividualDisplacements)),"../examples/recipes_plots.jl"))
# map(𝐼,OceanDepthLog(𝐷.Γ))
# ```


text: function print() end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
import Random: AbstractRNG, rand, Sampler

export GeneralTermSampler,
    large_terms,
    BoundedTermSampler

function S∞(z)
    sq = z^6 + 2z^5 - 5z^4 + 4z^3 - z^2 - 2z + 1
    num = z^3 - z^2 - z + 1 - √sq
    den = 2z^2 * (1 - z)
    num / den
end

p₁(x) = x^2 / (1 - x) / S∞(x)
p₂(x) = p₁(x) + x^2

const ρ = 0.509308127024237357194177485
const ρ² = ρ * ρ
const p₁ρ = (1 - ρ²) / 2
const p₂ρ = p₁ρ + ρ²

struct GeneralTermSampler <: Sampler{Term}
    x::Float64
    p₁::Float64
    p₂::Float64
end

GeneralTermSampler(x) = GeneralTermSampler(x, p₁(x), p₂(x))

const large_terms = GeneralTermSampler(ρ, p₁ρ, p₂ρ)

function rand_index(rng::AbstractRNG, s::GeneralTermSampler)
    result = 1

    while rand(rng) < s.x
        result += 1
    end

    return result
end

function rand(rng::AbstractRNG, s::GeneralTermSampler)
    p = rand(rng)

    if p < s.p₁
        return Var(rand_index(rng, s))
    elseif p < s.p₂
        return Abs(rand(rng, s))
    else
        App(rand(rng, s), rand(rng, s))
    end
end

struct BoundedTermSampler <: Sampler{Term}
    lower::Int
    upper::Int
end

function rand_ceiled(rng::AbstractRNG, s::BoundedTermSampler)::Union{Term, Nothing}
    p = rand(rng)

    if p < p₁ρ
        i = rand_index(rng, large_terms)
        return i < s.upper ? Var(i) : nothing
    elseif p < p₂ρ
        mbt = rand_ceiled(rng, s)
        if mbt !== nothing && length(mbt) + 2 ≤ s.upper
            return Abs(mbt)
        else
            return nothing
        end
    else
        mbt1, mbt2 = rand_ceiled(rng, s), rand_ceiled(rng, s)
        if mbt1 !== nothing && mbt2 !== nothing && length(mbt1) + length(mbt2) + 2 ≤ s.upper
            return App(mbt1, mbt2)
        else
            return nothing
        end
    end
end

function rand(rng::AbstractRNG, s::BoundedTermSampler)
    candidate = rand_ceiled(rng, s)

    while candidate == nothing || length(candidate) < s.lower
        candidate = rand_ceiled(rng, s)
    end

    return candidate
end


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

using MXNet

#--------------------------------------------------------------------------------
# Helper functions to construct larger networks

# basic Conv + BN + ReLU factory
function conv_factory(data, num_filter, kernel; stride=(1,1), pad=(0,0), act_type=:relu)
  conv = mx.Convolution(data, num_filter=num_filter, kernel=kernel, stride=stride, pad=pad)
  bn   = mx.BatchNorm(conv)
  act  = mx.Activation(bn, act_type=act_type)
  return act
end

# simple downsampling factory
function downsample_factory(data, ch_3x3)
  # conv 3x3
  conv = conv_factory(data, ch_3x3, (3,3), stride=(2,2), pad=(1,1))
  # pool
  pool = mx.Pooling(data, kernel=(3,3), stride=(2,2), pool_type=:max)
  # concat
  concat = mx.Concat(conv, pool)
  return concat
end

# a simple module
function simple_factory(data, ch_1x1, ch_3x3)
  # 1x1
  conv1x1 = conv_factory(data, ch_1x1, (1,1); pad=(0,0))
  # 3x3
  conv3x3 = conv_factory(data, ch_3x3, (3,3); pad=(1,1))
  # concat
  concat = mx.Concat(conv1x1, conv3x3)
  return concat
end

#--------------------------------------------------------------------------------
# Actual architecture
data    = mx.Variable(:data)
conv1   = conv_factory(data, 96, (3,3); pad=(1,1), act_type=:relu)
in3a    = simple_factory(conv1, 32, 32)
in3b    = simple_factory(in3a, 32, 48)
in3c    = downsample_factory(in3b, 80)
in4a    = simple_factory(in3c, 112, 48)
in4b    = simple_factory(in4a, 96, 64)
in4c    = simple_factory(in4b, 80, 80)
in4d    = simple_factory(in4b, 48, 96)
in4e    = downsample_factory(in4d, 96)
in5a    = simple_factory(in4e, 176, 160)
in5b    = simple_factory(in5a, 176, 160)
pool    = mx.Pooling(in5b, pool_type=:avg, kernel=(7,7), name=:global_pool)
flatten = mx.Flatten(pool, name=:flatten1)
fc      = mx.FullyConnected(flatten, num_hidden=10, name=:fc1)
softmax = mx.SoftmaxOutput(fc, name=:loss)

#--------------------------------------------------------------------------------
# Prepare data
filenames = mx.get_cifar10()
batch_size = 128
num_epoch  = 10
num_gpus   = 8

train_provider = mx.ImageRecordProvider(label_name=:loss_label,
        path_imgrec=filenames[:train], mean_img=filenames[:mean],
        rand_crop=true, rand_mirror=true, data_shape=(28,28,3),
        batch_size=batch_size, preprocess_threads=1)
test_provider = mx.ImageRecordProvider(label_name=:loss_label,
        path_imgrec=filenames[:test], mean_img=filenames[:mean],
        rand_crop=false, rand_mirror=false, data_shape=(28,28,3),
        batch_size=batch_size, preprocess_threads=1)

#--------------------------------------------------------------------------------
# Training model
gpus = [mx.Context(mx.GPU, i) for i = 0:num_gpus-1]
model = mx.FeedForward(softmax, context=gpus)

# optimizer
optimizer = mx.SGD(η=0.05, μ=0.9, λ=0.0001)

# fit parameters
mx.fit(model, optimizer, train_provider, n_epoch=num_epoch, eval_data=test_provider,
       initializer=mx.UniformInitializer(0.07), callbacks=[mx.speedometer()])

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
module INE

using HTTP, JSONStat
using Dates

export INEdataset

function INEdataset(id::String; from::Date=Date(2010, 12))
    resp = HTTP.get("https://servicios.ine.es/wstempus/jsstat/ES/DATASET/" *
                    id *"?date=" * Dates.format(from, "Ymmdd") *":").body

    return JSONStat.read(resp)
end

INEdataset(id::Int; from::Date=Date(2010, 12)) = INEdataset(string(id); from)

function INEmetadata() end

end # module

max_stars_repo_name: lsaenzt/INE

max_issues_repo_name: lsaenzt/INE

max_forks_repo_name: lsaenzt/INE


text:
# ---
# title: 678. Valid Parenthesis String
# id: problem678
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/valid-parenthesis-string/description/>
# hidden: true
# ---
#
# Given a string containing only three types of characters: '(', ')' and '*',
# write a function to check whether this string is valid. We define the validity
# of a string by these rules:
#
#   1. Any left parenthesis `'('` must have a corresponding right parenthesis `')'`.
#   2. Any right parenthesis `')'` must have a corresponding left parenthesis `'('`.
#   3. Left parenthesis `'('` must go before the corresponding right parenthesis `')'`.
#   4. `'*'` could be treated as a single right parenthesis `')'` or a single left parenthesis `'('` or an empty string.
#   5. An empty string is also valid.
#
# **Example 1:**
#
#
#
#     Input: "()"
#     Output: True
#
#
# **Example 2:**
#
#
#
#     Input: "(*)"
#     Output: True
#
#
# **Example 3:**
#
#
#
#     Input: "(*))"
#     Output: True
#
#
# **Note:**
#
#   1. The string size will be in the range [1, 100].
#
#
## @lc code=start
using LeetCode

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


text:
@testset "UTM functions" begin
    # Test zone prediction
    @test Geodesy.utm_zone(1.0, 1.0) == (31, true)
    @test Geodesy.utm_zone(-1.0, -1.0) == (30, false)
    @test Geodesy.utm_zone(89.0, -1.0) == (0, true)
    @test Geodesy.utm_zone(-89.0, 1.0) == (0, false)

    # Svalbard
    @test Geodesy.utm_zone(75.0, 8.0) == (31, true)
    @test Geodesy.utm_zone(83.9, 10.0) == (33, true)
    @test Geodesy.utm_zone(72.1, 20.0) == (33, true)
    @test Geodesy.utm_zone(75.0, 22.0) == (35, true)
    @test Geodesy.utm_zone(75.0, 32.0) == (35, true)
    @test Geodesy.utm_zone(75.0, 34.0) == (37, true)

    # Noway
    @test Geodesy.utm_zone(60.0, 2.9) == (31, true)
    @test Geodesy.utm_zone(60.0, 3.1) == (32, true)

    # From types
    @test Geodesy.utm_zone(LatLon(1.0, 1.0)) == (31, true)
    @test Geodesy.utm_zone(LLA(1.0, 1.0)) == (31, true)
    @test Geodesy.utm_zone(ECEF(LLA(1.0, 1.0), wgs84), wgs84) == (31, true)

    # central meridians
    @test Geodesy.utm_meridian(31) == 3
    @test_throws Exception Geodesy.utm_meridian(0)
    @test_throws Exception Geodesy.utm_meridian(61)

end


text:
# error to envoke warning message
Dict(1)


text:
module TestBasicOperations

using AbstractWrappedDataFrames
using DataFrames
using Test

struct WDF <: AbstractWrappedDataFrame
    df::DataFrame
end

struct NDF <: AbstractWrappedDataFrame
    df::DataFrame
    name::String
end

# This is required to support various DataFrames join operations
NDF(df::DataFrame) = NDF(df, "No name")

function test_operations(label::AbstractString, make_wdf::Function, args...)

    @testset "$label" begin

        df1 = make_wdf(DataFrame(x = [1,2], y = [4,5]), args...)
        df2 = make_wdf(DataFrame(x = [1,4], z = [7,8]), args...)

        @testset "Metadata" begin
            @test sort(names(df1)) == ["x", "y"]
            @test sort(string.(propertynames(df1))) == ["x", "y"]
            @test nrow(df1) == 2
            @test ncol(df1) == 2
        end

        @testset "Filter" begin
            @test first(df1, 1).x == [1]
            @test last(df1, 1).x == [2]
            @test filter(:x => ==(1), df1) |> nrow == 1
        end

        @testset "Add columns" begin
            sf3 = copy(df1)
            @test_nowarn sf3.z1 = 10
            @test_nowarn sf3.z2 = [10,11]
            @test_nowarn sf3.z3 = sf3.z2 .* 2
        end

        @testset "Summarize" begin
            @test select(df1, :x) |> ncol == 1
            @test transform(df1, :x => (x -> x .+ 1) => :q) |> ncol == 3
            @test combine(df1, :x => sum) |> nrow == 1
            @test groupby(df1, :x) |> keys |> length == 2
        end

        @testset "Joins" begin
            @test innerjoin(df1, df2, on = :x) isa DataFrame
            @test leftjoin(df1, df2, on = :x) isa DataFrame
            @test rightjoin(df1, df2, on = :x) isa DataFrame
            @test outerjoin(df1, df2, on = :x) isa DataFrame
            @test semijoin(df1, df2, on = :x) isa DataFrame
            @test antijoin(df1, df2, on = :x) isa DataFrame
        end

        @testset "Return types" begin
            @test select(df1, :x) isa DataFrame
            @test transform(df1, :x => (x -> x .+ 1) => :q) isa DataFrame
            @test combine(df1, :x => sum) isa DataFrame
            @test groupby(df1, :x) isa GroupedDataFrame

            sf3 = copy(df1)
            @test select!(sf3, :x) isa DataFrame
            @test transform!(sf3, :x => identity) isa DataFrame
        end
    end
end

function test()
    @testset "Basic Operations" begin
        test_operations("Simple", df -> WDF(df))
        test_operations("Complex", df -> NDF(df, "MyName"))
    end
end

end # module

using .TestBasicOperations
TestBasicOperations.test()


text:
export NeuralNetworkApproximator, ActorCritic

using Flux

"""
    NeuralNetworkApproximator(;kwargs)

Use a DNN model for value estimation.

# Keyword arguments

- `model`, a Flux based DNN model.
- `optimizer=nothing`
"""
Base.@kwdef struct NeuralNetworkApproximator{M,O} <: AbstractApproximator
    model::M
    optimizer::O = nothing
end

(app::NeuralNetworkApproximator)(x) = app.model(x)

# !!! watch https://github.com/FluxML/Functors.jl/blob/master/src/functor.jl#L2
Flux.functor(x::NeuralNetworkApproximator) =
    (model = x.model,), y -> NeuralNetworkApproximator(y.model, x.optimizer)

device(app::NeuralNetworkApproximator) = device(app.model)

RLBase.update!(app::NeuralNetworkApproximator, gs) =
    Flux.Optimise.update!(app.optimizer, params(app), gs)

Base.copyto!(dest::NeuralNetworkApproximator, src::NeuralNetworkApproximator) =
    Flux.loadparams!(dest.model, params(src))

Flux.testmode!(app::NeuralNetworkApproximator, mode = true) = testmode!(app.model, mode)

#####
# ActorCritic
#####

"""
    ActorCritic(;actor, critic, optimizer=ADAM())

The `actor` part must return logits (*Do not use softmax in the last layer!*), and the `critic` part must return a state value.
"""
Base.@kwdef struct ActorCritic{A,C,O} <: AbstractApproximator
    actor::A
    critic::C
    optimizer::O = ADAM()
end

Flux.functor(x::ActorCritic) =
    (actor = x.actor, critic = x.critic), y -> ActorCritic(y.actor, y.critic, x.optimizer)

RLBase.update!(app::ActorCritic, gs) = Flux.Optimise.update!(app.optimizer, params(app), gs)

function Flux.testmode!(app::ActorCritic, mode = true)
    testmode!(app.actor, mode)
    testmode!(app.critic, mode)
end


text:
using GLMakie
using LaTeXStrings
using JLD2
using Random
using Statistics
using Printf

"""
function mega_viz!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange;
    add_labels = false, title_string = "", colormap = :balance,
    heuristic = 1, random_seed = 12345)

# Description
Combines a contour plot with a heatmap and adds numbers to the contours.

"""
function contour_heatmap!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange;
    add_labels = false, colormap = :balance,
    heuristic = 1, random_seed = 12345)

    cplot = contour!(ax, ϕ, p_coord, slice_zonal, levels = contour_levels, color = :black, interpolate = true, show_axis = false)
    hm = heatmap!(ax, ϕ, p_coord, slice_zonal, colorrange = colorrange, colormap = colormap, interpolate = true)

    ax.limits = (extrema(ϕ)..., extrema(p_coord)...)

    ax.xlabel = "Latitude [ᵒ]"
    ax.ylabel = "Stretched Height"
    ax.xlabelsize = 25
    ax.ylabelsize = 25
    ax.xticks = ([-80, -60, -30, 0, 30, 60, 80], ["80S", "60S", "30S", "0", "30N", "60N", "80N"])
    pressure_levels = [1000, 850, 700, 550, 400, 250, 100, 10]
    ax.yticks = (pressure_levels .* 1e2, string.(pressure_levels))
    ax.yreversed = true

    # hack
    Random.seed!(random_seed)
    if add_labels
        list_o_stuff = []
        labeled_contours = contour_levels[1:1:end]
        for level in labeled_contours
            local fig_t, ax_t, cp_t = contour(ϕ, p_coord, slice_zonal, levels = [level], linewidth = 0)
            local segments = cp_t.plots[1][1][]
            local index_vals = []
            local beginnings = []
            for (i, p) in enumerate(segments)
                # the segments are separated by NaN, which signals that a new contour starts
                if isnan(p)
                    push!(beginnings, segments[i-1])
                    push!(index_vals, i)
                end
            end
            push!(list_o_stuff, (; segments, beginnings, index_vals))
        end

        for contour_index = 1:length(labeled_contours)

            local contour_val = labeled_contours[contour_index]
            local segments = list_o_stuff[contour_index].segments

            local indices = [0, list_o_stuff[contour_index].index_vals[1:end]...]
            for i = 1:length(indices)-1
                # heuristics for choosing where on line
                local index1 = rand(indices[i]+1:indices[i+1]-1) # choose random point on segment
                local index2 = round(Int, 0.5 * indices[i] + 0.5 * indices[i+1]) # choose point in middle
                β = (rand() - 0.5) * 0.9 + 0.5 # α ∈ [0,1]
                # α = (contour_index-1) / (length(labeled_contours)-1)
                α = contour_index % 2 == 0 ? 0.15 : 0.85
                α = rand([α, β])
                local index3 = round(Int, α * (indices[i] + 1) + (1 - α) * (indices[i+1] - 1)) # choose point in middle
                if heuristic == 3
                    local index = index3 # rand([index1, index2]) # choose between random point and point in middle
                elseif heuristic == 1
                    local index = index1
                elseif heuristic == 2
                    local index = index2
                end
                # end of heuristics
                local location = Point3(segments[index]..., 2.0f0)
                local sc = scatter!(ax, location, markersize = 20, align = (:center, :center), color = (:white, 0.1), strokecolor = :white)
                local anno = text!(ax, [("$contour_val", location)], align = (:center, :center), textsize = 20, color = :black)

                delete!(ax, sc)
                delete!(ax, cplot)
                delete!(ax, anno)

                push!(ax.scene, anno)
                push!(ax.scene, sc)
                push!(ax.scene, cplot)
            end
        end
    end # end of adding labels
end

function grab_state(s_string, jl_file)
    m_string = "moment_" * "$(length(s_string))"
    m_names = jl_file[m_string*"_names"]
    m_ind = findall(x -> x == s_string, m_names)[1]
    state = jl_file[m_string][1:end-1, :, :, m_ind]
    return mean(state, dims = 1)[1, :, :] ./ jl_file["times"]
end

function grab_grid(jl_file)
    lat_grid = jl_file["grid"]["latitude"]
    lon_grid = jl_file["grid"]["longitude"]
    rad_grid = jl_file["grid"]["radius"]

    # 6 is pressure
    eq_ind = argmin(abs.(lat_grid)) # pressure at equator
    p2 = mean(jl_file["moment_1"][1:end-1, eq_ind, :, 6], dims = 1)[1, :] ./ jl_file["times"]

    λ = lon_grid
    ϕ = lat_grid
    r = rad_grid

    p0 = 1e5 # pressure bottom
    pH = 800 # pressure top
    H = rad_grid[end] - rad_grid[1] # height of domain
    sH = -H / log(pH / p0) # scale height
    # roughly p = p0 * exp(-z / sH), so log(p / p0) * sH = z
    p = p0 * exp.(-(rad_grid .- rad_grid[1]) ./ sH)
    return λ, ϕ, r, p2
end

function eddy_variance(m2_string, jl_file)
    @assert length(m2_string) == 2
    string_1 = m2_string[1]
    string_2 = m2_string[end]

    state_1 = grab_state(string(string_1), jl_file)
    state_2 = grab_state(string(string_2), jl_file)
    state_12 = grab_state(m2_string, jl_file)
    state = state_12 .- state_1 .* state_2
    name = string_1 * "'" * string_2 * "'"
    return state, name
end

function plot_helper(s_string, slice_zonal)
    zonal_quantiles = extrema(slice_zonal[:])
    colorrange = zonal_quantiles
    # contour_levels = collect(0:20:120)
    if s_string == "u"
        contour_levels = collect(-8:4:28) # u-velocity
        latex_string = L"$\langle u$ $\rangle$"
    elseif s_string == "T"
        contour_levels = collect(190:10:300) # T
        latex_string = L"$\langle T$ $\rangle$"
    elseif s_string == "T'T'"
        contour_levels = collect(0:8:40) # TT
        latex_string = L"\langle T' T' \rangle"
    elseif s_string == "u'u'"
        # contour_levels = collect(0:40:160) # uu
        contour_levels = [0, 40, 80, 160, 240, 320]
        latex_string = L"\langle u' u' \rangle"
    elseif s_string == "v'T'"
        contour_levels = collect(-21:3:21) # vT
        latex_string = L"\langle v' T' \rangle"
    elseif s_string == "u'v'"
        contour_levels = collect(-30:10:30) # uv
        latex_string = L"\langle u' v' \rangle"
    else
        contour_levels = range(colorrange[1], colorrange[2], length = 10) # default
        latex_string = " "
    end
    return colorrange, contour_levels, latex_string
end

function latex_format(matrix)
    m, n = size(matrix)
    for i = 1:m
        for j = 1:n
            # @sprintf print(error_matrix[i, j], " ")
            tmp = matrix[i, j] * 100
            if j < n
                @printf " %2.0f\\%% &" tmp
            else
                @printf " %2.0f\\%% " tmp
            end
        end
        println(" ")
    end
end

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
"""
Reference Finite Element without change-of-basis, all the DOFs belong to the
polytope itself
"""
struct DiscRefFE{D,T} <: RefFE{D,T}
  polytope::Polytope{D}
  shfbasis::Basis{D,T}
  nfacedofs::Vector{Vector{Int}}
  # this type is unstable
end

"""
Constructor that given a vector of orders, generates high-order Lagrangian
RefFEs on n-cubes or n-tets.
"""
function DiscRefFE(::Type{T}, p::Polytope{D}, orders::Vector{Int}) where {D,T}
  @assert length(orders) == D
  @assert D > 0
  # @santiagobadia : This part is missing for anisotropic orders !
  basis = Gridap.RefFEs._monomial_basis(p,T,orders[1])
  # aux = zeros(Float64,numlocaldofs(dofsb),numlocaldofs(dofsb))
  # @assert numlocaldofs(dofsb) == length(basis)
  return DiscRefFE{D,T}(p, basis, nfacedofs)
end

"""
Version of the constructor for a scalar order
"""
function DiscRefFE(::Type{T}, p::Polytope{D}, order::Int) where {D,T}
  _order = order*ones(Int,D)
  return DiscRefFE(T,p,_order)
end


text:
using Gurobi,GLPK
using JuMP
using PMDPs
using POMDPSimulators
using POMDPPolicies
using POMDPs
using Suppressor

"""
Get maximal action value that is below the budget b.
"""
function optimal_price(b::Float64, actions::AbstractArray{<:Number})::PMDPs.Action
    best_action = actions[1]
    for a in actions
        if a > b
            return best_action
        else
            best_action=a
        end
    end
    return best_action
end

"""
Calculate hindsight optimal allocation of requests from history recorder simulation history.
Optimization goal can be either "revenue" or "utilization".
Return the optimal allocation and the value of the objective.
Usage:

mdp_mc = PMDPs.create_PMDP(PMDPg)
planner = PMDPs.get_MCTS_planner(mdp_mc)

rng = MersenneTwister(123)
hr = HistoryRecorder(max_steps=100, capture_exception=false, rng=rng)
h = simulate(hr, mdp_mc, planner)

MILP_hindsight_pricing(mdp_mc, h)
"""
function MILP_hindsight_pricing(mdp::PMDPs.PMDP, h::AbstractSimHistory;
            objective::Symbol=:revenue, verbose::Bool=false,
            kwargs...)

    kwargs = Dict(kwargs)

    # extract request trace from history
    trace = collect(eachstep(h, "s, info"))
    requests = [rec for rec in trace if rec.s.iₚ != PMDPs.empty_product_id(mdp)]

    if length(requests)==0
        result = (r = 0., u = 0., alloc = [], action_seq=[], requests=[])
    else
        # get data from trace
        request_resources = [[PMDPs.product(mdp, rec.s)...] for rec in requests]
        request_budgets = [rec.info.b for rec in requests]

        # prepare data
        """ E: Matrix that has requests as rows and columns as capacity edges
        R: List of maximum prices that users will accept
        e.g. matrix for 4 request in problem with 3 edges and capacity 3:
        E = [[1,1,0],
        [1,0,0],
        [1,1,1],
        [0,0,1]]
        R = [20, 10, 30, 10]
        capacity_constraints = [3,3,3]
        """
        E = request_resources
        R = [optimal_price(b, POMDPs.actions(mdp)) for b in request_budgets]

        capacity_constraints = [PMDPs.pp(mdp).c₀...]

        request_ind = 1:length(request_resources)
        capacity_ind = 1:length(capacity_constraints)

        local model
        @suppress_out begin
            if haskey(kwargs, :gurobi) && kwargs[:gurobi]
                if haskey(kwargs, :env)
                    model = Model(() -> Gurobi.Optimizer(kwargs[:env]))
                else
                    model = Model(Gurobi.Optimizer)
                end
                if ~verbose  set_optimizer_attribute(model, "OutputFlag", 0) end
            else
                model = Model(GLPK.Optimizer)
            end
        end
        # set_optimizer_attribute(model, "Presolve", 0)

        # Variables
        @variable(model, x[request_ind], Bin)

        # Constraints
        @constraint(model, [j in capacity_ind],
                sum(x[i]*E[i][j] for i in request_ind) <= capacity_constraints[j])

        # Objectives
        if objective == :revenue
            @objective(model, Max,
                sum(x[i]*R[i] for i in request_ind))
        elseif  objective == :utilization
            durations = [sum(rec) for rec in E]
            @objective(model, Max,
                sum(x[i]*durations[i] for i in request_ind))
        else
            error("unknown MILP optimization goal $optimization_goal")
        end

        # Optimize and capture std_out
        output = @capture_out optimize!(model)

        # Report results
        status = termination_status(model)
        obj_val = objective_value(model)
        optimal_alloc = JuMP.value.(x)

        # Calculate utilization
        utilization = sum(sum([x*p for (x,p) in zip(optimal_alloc.data, request_resources)]))

        if verbose
            print(output)
            println("Allocation: ", optimal_alloc.data)
        end

        action_seq = [(a!=0. ? a : PMDPs.REJECT_ACTION) for a in optimal_alloc.data.*R]

        result = (objective_val = obj_val,
                  u = utilization,
                  alloc = optimal_alloc.data,
                  action_seq = action_seq,
                  requests = requests)
    end
    return result
end

"""
Returns hindsight based policy for given history.

This is the preffered interface to getting the hindsight actions.
"""
function get_MILP_hindsight_policy(mdp::PMDPs.PMDP, h::AbstractSimHistory; kwargs...)
    (objective_val, u, alloc, action_seq, requests) = MILP_hindsight_pricing(mdp, h; objective=PMDPs.objective(mdp), kwargs...)

    timesteps = collect([req[:s].t for req in requests])
    pd = Dict(zip(timesteps, action_seq))

    hp = FunctionPolicy(s -> get(pd, s.t, PMDPs.REJECT_ACTION))
end


text:
#Should work with
#Julia 1.4
#Atom v0.12.10
#CUDAapi v4.0.0
#CUDAdrv v6.2.2
#CUDAnative v3.0.2
#CuArrays v2.0.1 #master (https://github.com/JuliaGPU/CuArrays.jl)
#DataFrames v0.20.2
#Flux v0.10.3 #master (https://github.com/FluxML/Flux.jl)
#ImageMagick v1.1.3
#Images v0.22.0
#Juno v0.8.1
#MLDatasets v0.4.1
#Metalhead v0.5.0
#NNlib v0.6.6
#RDatasets v0.6.1
#StatsBase v0.32.2
#Zygote v0.4.12
#
#
#
#Still has issues with speed and memory consumption.
#ConvTranspose slows the neural network down significantly.
#A more optimal way to feed images to the neural network probably can be found.
#resizing images in preprocessing should be faster but might lead to greater
#memory consumption. Resizing in preprocessing could be done with ConvTranspose,
#imresize, or by some other tool.

ENV["JULIA_CUDA_VERBOSE"] = true
ENV["JULIA_CUDA_MEMORY_POOL"] = "split"
ENV["JULIA_CUDA_MEMORY_LIMIT"] = 8000_000_000

using Random
using Statistics
using CuArrays
using Zygote
using Flux, Flux.Optimise
using Metalhead, Images
using Metalhead: trainimgs
using Images.ImageCore
using Flux: onehotbatch, onecold, logitcrossentropy, Momentum, @epochs
using Base.Iterators: partition
using Dates

CuArrays.allowscalar(false)

batch_size = 1

getarray(X) = Float32.(permutedims(channelview(X), (2, 3, 1)))

function make_minibatch(imgs,labels,batch_size)
  data_set = [(cat(imgs[i]..., dims = 4), labels[:,i])
          for i in partition(1:length(imgs), batch_size)]
  return data_set
end

X = trainimgs(CIFAR10)

train_idxs = 1:49000

train_imgs = [getarray(X[i].img) for i in train_idxs]
train_labels = Float32.(onehotbatch([X[i].ground_truth.class for i in train_idxs],1:10))
train_set = make_minibatch(train_imgs,train_labels,batch_size)
train_data = train_set |>
  x -> map(y->gpu.(y),x)

verify_idxs = 49001:50000
verify_imgs = cat([getarray(X[i].img) for i in verify_idxs]..., dims = 4)
verify_labels = Float32.(onehotbatch([X[i].ground_truth.class for i in verify_idxs],1:10))
verify_set = [(verify_imgs,verify_labels)]
verify_data = verify_set |>
  x -> map(y->gpu.(y),x)

tX = valimgs(CIFAR10)
test_idxs = 1:10000
test_imgs = [getarray(tX[i].img) for i in test_idxs]
test_labels = Float32.(onehotbatch([tX[i].ground_truth.class for i in test_idxs], 1:10))
test_set = make_minibatch(test_imgs,test_labels,batch_size)
test_data = test_set |>
  x -> map(y->gpu.(y),x)

identity_layer(n) = Chain(Conv((3,3), n=>n, pad = (1,1), stride = (1,1)),
                                  BatchNorm(n,relu),
                                  Conv((3,3), n=>n, pad = (1,1), stride = (1,1)),
                                  BatchNorm(n,relu))

convolution_layer(n) = Chain(Conv((3,3), n=> 2*n, pad = (1,1), stride = (2,2)),
                             BatchNorm(2*n,relu),
                             Conv((3,3), 2*n=>2*n, pad = (1,1), stride = (1,1)),
                             BatchNorm(2*n,relu))

simple_convolution(n) = Chain(Conv((1,1), n=>n, pad = (1,1), stride = (2,2)),
                              BatchNorm(n,relu))

m_filter(n) = Chain(
  Conv((3,3), n=>2*n, pad = (1,1), stride = (2,2)),
  BatchNorm(2*n,relu),
) |> gpu

struct Combinator
    conv::Chain
end |> gpu
Combinator(n) = Combinator(m_filter(n))# |> gpu

function (op::Combinator)(x, y)
  z = op.conv(y)
  return x + z
end

n = 7

m = Chain(
  ConvTranspose((n, n), 3 => 3, stride = n),
  Conv((7,7), 3=>64, pad = (3,3), stride = (2,2)),
  BatchNorm(64,relu),
  MaxPool((3,3), pad = (1,1), stride = (2,2)),
  SkipConnection(identity_layer(64), (variable_1, variable_2) -> variable_1 + variable_2),
  SkipConnection(identity_layer(64), (variable_1, variable_2) -> variable_1 + variable_2),
  SkipConnection(convolution_layer(64), Combinator(64)),
  SkipConnection(identity_layer(128), (variable_1, variable_2) -> variable_1 + variable_2),
  SkipConnection(convolution_layer(128), Combinator(128)),
  SkipConnection(identity_layer(256), (variable_1, variable_2) -> variable_1 + variable_2),
  SkipConnection(convolution_layer(256), Combinator(256)),
  SkipConnection(identity_layer(512), (variable_1, variable_2) -> variable_1 + variable_2),
  MeanPool((7,7)),
  x -> reshape(x, :, size(x,4)),
  Dense(512*1, 10),
  softmax,
) |> gpu

function accuracy(data_set)
  batch_size = size(data_set[1][1])[end]
  l = length(data_set)*batch_size
  s = 0f0
  for (x,y) in data_set
    s += sum((onecold(m(x|>gpu) |> cpu) .== onecold(y|>cpu)))
  end
  return s/l
end

loss(x, y) = sum(logitcrossentropy(m(x), y))
opt = Momentum(0.01)

number_of_epochs = 1

@epochs number_of_epochs train!(loss, params(m), train_data, opt, cb = Flux.throttle(() -> println("training... $(Dates.Time(Dates.now()))") , 10))

verify_acc = accuracy(verify_data)
@info "Verify accuracy : $(verify_acc)"

test_acc = accuracy(test_data)
@info "Test accuracy : $(test_acc)"

max_stars_repo_name: OTapio/model-zoo

max_issues_repo_name: OTapio/model-zoo

max_forks_repo_name: OTapio/model-zoo


text:
using SeisIO, SeisNoise
export partition

"""
    partition(A::Array{Float64,2}, offset::Int64, size::Int64)

Partition an array into two sub-arrays, centered on the midpoint.

# Arguments
- `A::Array{Float32,2},`    : Input array to be partitioned
- `offset::Int64`    : Number of samples away from the midpoint of the array to start the window.
- `size::Int64`    : Length of the window in number of samples.

# Output
- `A`::Array{Float32,2}    : Partitioned array.
"""
function partition(A::Array{Float32,2}, offset::Int64, size::Int64)
    # midpoint of the array
    zerolag = convert(Int64, (length(A[:,1])+1)/2)
    # define windows
    pos_win = zerolag+offset:zerolag+offset+size
    neg_win = zerolag-offset-size:zerolag-offset
    # stack subarrays
    A = hcat(A[neg_win, :], A[pos_win, :])

    return A
end

"""
    slide(C::CorrData, cc_len::Int64, cc_step::Int64)

Generate equal length sliding windows into a CorrData object for both positive and negative codas.

# Arguments
- `C::CorrData,`    : Input array to be partitioned
- `cc_len::Int64`    : Number of samples away from the midpoint of the array to start the window
- `cc_step::Int64`    : Length of the window in number of samples

# Output
- `corr_neg::Array{Float32,2}`    : Partitioned array of negative-lag data
- `corr_pos::Array{Float32,2}`    : Partitioned array of positive-lag data
"""
function slide_c3(C::CorrData, cc_len::Int64, cc_step::Int64)
    window_samples = Int(cc_len * C.fs)
    t=[1 0; length(C.corr[:, 1]) 0]
    su,eu = SeisIO.t_win(t, C.fs) * 1e-6
    starts = Array(range(su,stop=eu,step=cc_step))
    ends = starts .+ cc_len .- 1. / C.fs
    ind = findlast(x -> x <= eu,ends)
    starts = starts[1:ind]
    ends = ends[1:ind]

    # fill arrays with overlapping windows
    corr_neg = Array{Float64,2}(undef, window_samples,length(starts))
    corr_pos = Array{Float64,2}(undef, window_samples,length(starts))
    s = convert.(Int,round.((hcat(starts,ends) .- su) .* C.fs .+ 1.))
    for ii = 1:length(starts)
        corr_neg[:,ii] = C.corr[s[ii,1]:s[ii,2], 1]
        corr_pos[:,ii] = C.corr[s[ii,1]:s[ii,2], 2]
    end

    return corr_neg, corr_pos
end

"""
    window(A::Array{R,1}, win_len::Int64, win_step::Int64) where R<:Real

Generate equal length sliding windows into an array.

# Arguments
- `A::Array{R,1} where R<:Real,`    : Input array to be windowed
- `win_len::Int64`    : Number of samples in each window
- `win_step::Int64`    : Step size for successive windows

# Output
- `Awin::Array{R,2} where R<:Real`    : Windowed array
"""
function window(A::Array{R,1}, win_len::Int64, win_step::Int64) where R<:Real
    # left-most indices of each window
    minind = 1:win_len:length(A) - win_len

    # number of windows
    N = length(minind)

    # initialize array of size (win_len, nWindows)
    Awin = zeros(win_len, N)
    # fill matrix
    for i=1:N
        Awin[:,i] = A[minind[i]:minind[i]+win_len-1]
    end

    return Awin
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug, @sync, @context!, initialize_context

using GPUArrays

using LinearAlgebra

using BFloat16s: BFloat16

using CEnum: @cenum

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")
include("libcublas_deprecated.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

function math_mode!(handle, mode)
    flags = 0

    # https://github.com/facebookresearch/faiss/issues/1385
    if version() > v"11"
        flags = CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION
    end

    flags |= if mode == CUDA.PEDANTIC_MATH
        # prevent use of tensor cores
        if version() < v"11"
            CUBLAS_DEFAULT_MATH
        else
            CUBLAS_PEDANTIC_MATH
        end
    elseif mode == CUDA.DEFAULT_MATH
        # use tensor cores, but don't reduce precision
        if version() < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_DEFAULT_MATH
        end
    elseif mode == CUDA.FAST_MATH
        # we'll additionally select a compute-mode with reduced precision whenever possible
        if version() < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_TF32_TENSOR_OP_MATH
        end
    end

    cublasSetMathMode(handle, flags)

    return
end

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,cublasHandle_t}()
const idle_xt_handles = HandleCache{Any,cublasXtHandle_t}()

function handle()
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{handle::cublasHandle_t, stream::CuStream, math_mode::CUDA.MathMode}
    states = get!(task_local_storage(), :CUBLAS) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_handle = pop!(idle_handles, cuda.context) do
            cublasCreate()
        end

        finalizer(current_task()) do task
            push!(idle_handles, cuda.context, new_handle) do
                @context! skip_destroyed=true cuda.context cublasDestroy_v2(new_handle)
            end
        end

        cublasSetStream_v2(new_handle, cuda.stream)

        math_mode!(new_handle, cuda.math_mode)

        (; handle=new_handle, cuda.stream, cuda.math_mode)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    # update stream
    @noinline function update_stream(cuda, state)
        cublasSetStream_v2(state.handle, cuda.stream)
        (; state.handle, stream=cuda.stream, state.math_mode)
    end
    if state.stream != cuda.stream
        states[cuda.context] = state = update_stream(cuda, state)
    end

    # update math mode
    @noinline function update_math_mode(cuda, state)
        math_mode!(state.handle, cuda.math_mode)
        (; state.handle, state.stream, math_mode=cuda.math_mode)
    end
    if state.math_mode != cuda.math_mode
        states[cuda.context] = state = update_math_mode(cuda, state)
    end

    return state.handle
end

function xt_handle()
    cuda = CUDA.active_state()

    # every task maintains library state per set of devices
    LibraryState = @NamedTuple{handle::cublasXtHandle_t}
    states = get!(task_local_storage(), :CUBLASxt) do
        Dict{UInt,LibraryState}()
    end::Dict{UInt,LibraryState}

    # derive a key from the selected devices
    key = zero(UInt)
    for dev in devices()
        # we hash the device context to support device resets
        key = hash(context(dev), key)
    end

    # get library state
    @noinline function new_state(cuda)
        new_handle = pop!(idle_xt_handles, cuda.context) do
            cublasXtCreate()
        end

        finalizer(current_task()) do task
            push!(idle_xt_handles, cuda.context, new_handle) do
                # TODO: which context do we need to destroy this on?
                cublasXtDestroy(handle)
            end
        end

        devs = convert.(Cint, devices())
        cublasXtDeviceSelect(new_handle, length(devs), devs)

        (; handle=new_handle)
    end
    state = get!(states, key) do
        new_state(cuda)
    end

    return state.handle
end

function devices!(devs::Vector{CuDevice})
    task_local_storage(:CUBLASxt_devices, sort(devs; by=deviceid))
    return
end

devices() = get!(task_local_storage(), :CUBLASxt_devices) do
    # by default, select all devices
    sort(collect(CUDA.devices()); by=deviceid)
end::Vector{CuDevice}

ndevices() = length(devices())

## logging

const MAX_LOG_BUFLEN = UInt(1024*1024)
const log_buffer = Vector{UInt8}(undef, MAX_LOG_BUFLEN)
const log_cursor = Threads.Atomic{UInt}(0)
const log_cond = Ref{Base.AsyncCondition}()    # root

function log_message(ptr)
    # NOTE: this function may be called from unmanaged threads (by cublasXt),
    #       so we can't even allocate, let alone perform I/O.
    len = @ccall strlen(ptr::Cstring)::Csize_t
    old_cursor = log_cursor[]
    new_cursor = old_cursor + len+1
    if new_cursor >= MAX_LOG_BUFLEN
        # overrun
        return
    end

    @ccall memmove((pointer(log_buffer)+old_cursor)::Ptr{Nothing},
                   pointer(ptr)::Ptr{Nothing}, (len+1)::Csize_t)::Nothing
    log_cursor[] = new_cursor   # the consumer handles CAS'ing this value
    @ccall uv_async_send(log_cond[]::Ptr{Nothing})::Cint

    return
end

function _log_message(blob)
    # the message format isn't documented, but it looks like a message starts with a capital
    # and the severity (e.g. `I!`), and subsequent lines start with a lowercase mark (`!i`)
    #
    # lines are separated by a \0 if they came in separately, but there may also be multiple
    # actual lines separated by \n in each message.
    for message in split(blob, r"[\0\n]+(?=[A-Z]!)")
        code = message[1]
        lines = split(message[3:end], r"[\0\n]+[a-z]!")
        submessage = join(lines, '\n')
        if code == 'I'
            @debug submessage
        elseif code == 'W'
            @warn submessage
        elseif code == 'E'
            @error submessage
        elseif code == 'F'
            error(submessage)
        else
            @info "Unknown log message, please file an issue.\n$message"
        end
    end
    return
end

function __runtime_init__()
    # register a log callback
    log_cond[] = Base.AsyncCondition() do async_cond
        blob = ""
        while true
            message_length = log_cursor[]
            blob = unsafe_string(pointer(log_buffer), message_length)
            if Threads.atomic_cas!(log_cursor, message_length, UInt(0)) == message_length
                break
            end
        end
        _log_message(blob)
        return
    end
    if (isdebug(:init, CUBLAS) || Base.JLOptions().debug_level >= 2) &&
       !Sys.iswindows() # NVIDIA bug #3321130
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end


text:
module TestDiagrammaticPrograms
using Test

using Catlab.CategoricalAlgebra, Catlab.Programs.DiagrammaticPrograms
using Catlab.Programs.DiagrammaticPrograms: NamedGraph, MaybeNamedGraph
using Catlab.Graphs.BasicGraphs: TheoryGraph

# Graphs
########

para_parsed = @graph begin
  s
  t
  s → t
  s → t
end
para = @acset MaybeNamedGraph{Symbol} begin
  V = 2
  E = 2
  src = [1,1]
  tgt = [2,2]
  vname = [:s, :t]
  ename = [nothing, nothing]
end
@test para_parsed == para

para_parsed = @graph NamedGraph{Symbol} begin
  x, y
  (f, g): x → y
end
para = @acset NamedGraph{Symbol} begin
  V = 2
  E = 2
  src = [1,1]
  tgt = [2,2]
  vname = [:x, :y]
  ename = [:f, :g]
end
@test para_parsed == para

tri_parsed = @graph NamedGraph{Symbol} begin
  v0, v1, v2
  fst: v0 → v1
  snd: v1 → v2
  comp: v0 → v2
end
tri = @acset NamedGraph{Symbol} begin
  V = 3
  E = 3
  src = [1,2,1]
  tgt = [2,3,3]
  vname = [:v0, :v1, :v2]
  ename = [:fst, :snd, :comp]
end
@test tri_parsed == tri

# Categories
############

Δ¹_parsed = @category begin
  V, E
  (δ₀, δ₁): V → E
  σ₀: E → V

  σ₀ ∘ δ₀ == id(V)
  σ₀ ∘ δ₁ == id(V)
end
Δ¹_graph = @acset NamedGraph{Symbol} begin
  V = 2
  E = 3
  src = [1,1,2]
  tgt = [2,2,1]
  vname = [:V, :E]
  ename = [:δ₀, :δ₁, :σ₀]
end
Δ¹ = FinCat(Δ¹_graph, [ [1,3] => empty(Path, Δ¹_graph, 1),
                        [2,3] => empty(Path, Δ¹_graph, 1) ])
@test Δ¹_parsed == Δ¹

# Diagrams
##########

C = FinCat(TheoryGraph)
F_parsed = @diagram C begin
  v::V
  (e1, e2)::E
  t::tgt : e1 → v
  s::src : e2 → v
end
J = FinCat(@acset NamedGraph{Symbol} begin
  V = 3
  E = 2
  src = [2,3]
  tgt = [1,1]
  vname = [:v, :e1, :e2]
  ename = [:t, :s]
end)
@test dom(F_parsed) == J
F = FinFunctor([:V,:E,:E], [:tgt, :src], J, C)
@test F_parsed == F

F_parsed = @diagram C begin
  v => V
  (e1, e2) => E
  t: e1 → v => tgt
  s: e2 → v => src
end
@test F_parsed == F

end


text:
"""
    CallbackFunction()

Set a generic CPLEX callback function.

Note: before accessing `MOI.CallbackVariablePrimal`, you must call either
`callbackgetcandidatepoint(model::Optimizer, cb_data, cb_where)` or
`callbackgetrelaxationpoint(model::Optimizer, cb_data, cb_where)`.
"""
struct CallbackFunction <: MOI.AbstractOptimizerAttribute end

function MOI.set(model::Optimizer, ::CallbackFunction, f::Function)
    if MOI.get(model, MOI.NumberOfThreads()) != 1
        @warn(
            "When using callbacks, make sure to set `NumberOfThreads` to `1` " *
            "using `MOI.set(model, MOI.NumberOfThreads(), 1)`. Bad things can" *
            " happen if you don't! Only ignore this message if you know what " *
            "you are doing."
        )
    end
    model.has_generic_callback = true
    context_mask =
        CPX_CALLBACKCONTEXT_THREAD_UP |
        CPX_CALLBACKCONTEXT_THREAD_DOWN |
        CPX_CALLBACKCONTEXT_LOCAL_PROGRESS |
        CPX_CALLBACKCONTEXT_GLOBAL_PROGRESS |
        CPX_CALLBACKCONTEXT_CANDIDATE |
        CPX_CALLBACKCONTEXT_RELAXATION
    CPLEX.cbsetfunc(model.inner, context_mask, (cb_data, cb_where) -> begin
        model.callback_state = CB_GENERIC
        f(cb_data, cb_where)
        model.callback_state = CB_NONE
    end)
    return
end

"""
    callbackgetcandidatepoint(model::Optimizer, cb_data, cb_where)

Load the solution at a CPX_CALLBACKCONTEXT_CANDIDATE node so that it can be
accessed using `MOI.CallbackVariablePrimal`.
"""
function callbackgetcandidatepoint(model::Optimizer, cb_data, cb_where)
    N = length(model.variable_info)
    resize!(model.callback_variable_primal, N)
    CPLEX.cbgetcandidatepoint(
        cb_data, model.callback_variable_primal, Cint(0), Cint(N - 1), Ref{Float64}()
    )
    return
end

"""
    callbackgetrelaxationpoint(model::Optimizer, cb_data, cb_where)

Load the solution at a CB_MIPNODE node so that it can be accessed using
`MOI.CallbackVariablePrimal`.
"""
function callbackgetrelaxationpoint(model::Optimizer, cb_data, cb_where)
    N = length(model.variable_info)
    resize!(model.callback_variable_primal, N)
    CPLEX.cbgetrelaxationpoint(
        cb_data, model.callback_variable_primal, Cint(0), Cint(N - 1), Ref{Float64}()
    )
    return
end

# ==============================================================================
#    MOI callbacks
# ==============================================================================

function default_moi_callback(model::Optimizer)
    if MOI.get(model, MOI.NumberOfThreads()) != 1
        # The current callback system isn't thread-safe. As a work-around, set
        # the number of threads to 1, regardless of what the user intended.
        # We only do this for the MOI callbacks, since we assume the user knows
        # what they are doing if they use a solver-dependent callback.
        MOI.set(model, MOI.NumberOfThreads(), 1)
    end
    return (cb_data, cb_where) -> begin
    if cb_where == CPX_CALLBACKCONTEXT_CANDIDATE
        if cbcandidateispoint(cb_data) == 0
            return  # No candidate point available
        end
        callbackgetcandidatepoint(model, cb_data, cb_where)
        if model.lazy_callback !== nothing
            model.callback_state = CB_LAZY
            model.lazy_callback(cb_data)
        end
    elseif cb_where == CPX_CALLBACKCONTEXT_RELAXATION
        callbackgetrelaxationpoint(model, cb_data, cb_where)
            if model.lazy_callback !== nothing
                model.callback_state = CB_LAZY
                model.lazy_callback(cb_data)
            end
            if model.user_cut_callback !== nothing
                model.callback_state = CB_USER_CUT
                model.user_cut_callback(cb_data)
            end
            if model.heuristic_callback !== nothing
                model.callback_state = CB_HEURISTIC
                model.heuristic_callback(cb_data)
            end
        end
        model.callback_state = CB_NONE
    end
end

function MOI.get(
    model::Optimizer,
    ::MOI.CallbackVariablePrimal{CallbackContext},
    x::MOI.VariableIndex
)
    return model.callback_variable_primal[_info(model, x).column]
end

# ==============================================================================
#    MOI.LazyConstraint
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.LazyConstraintCallback, cb::Function)
    model.lazy_callback = cb
    return
end

function MOI.submit(
    model::Optimizer,
    cb::MOI.LazyConstraint{CallbackContext},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64}}
)
    if model.callback_state == CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    elseif model.callback_state == CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64, typeof(f), typeof(s)}(f.constant))
    end
    indices, coefficients = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    CPLEX.cbrejectcandidate(
        cb.callback_data,
        Cint(1),
        Cint(length(coefficients)),
        Float64[rhs],
        Cchar[sense],
        Cint[0],
        indices .- Cint(1),
        coefficients
    )
    return
end

# ==============================================================================
#    MOI.UserCutCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.UserCutCallback, cb::Function)
    model.user_cut_callback = cb
    return
end

function MOI.submit(
    model::Optimizer,
    cb::MOI.UserCut{CallbackContext},
    f::MOI.ScalarAffineFunction{Float64},
    s::Union{MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64}}
)
    if model.callback_state == CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == CB_HEURISTIC
        throw(MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), cb))
    elseif !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64, typeof(f), typeof(s)}(f.constant))
    end
    indices, coefficients = _indices_and_coefficients(model, f)
    sense, rhs = _sense_and_rhs(s)
    CPLEX.cbaddusercuts(
        cb.callback_data,
        Cint(1),
        Cint(length(coefficients)),
        [rhs],
        [sense],
        Cint[0],
        indices .- Cint(1),
        coefficients,
        [CPX_USECUT_FILTER],
        Cint[0]
    )
    return
end

# ==============================================================================
#    MOI.HeuristicCallback
# ==============================================================================

function MOI.set(model::Optimizer, ::MOI.HeuristicCallback, cb::Function)
    model.heuristic_callback = cb
    return
end

function MOI.submit(
    model::Optimizer,
    cb::MOI.HeuristicSolution{CallbackContext},
    variables::Vector{MOI.VariableIndex},
    values::MOI.Vector{Float64}
)
    if model.callback_state == CB_LAZY
        throw(MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), cb))
    elseif model.callback_state == CB_USER_CUT
        throw(MOI.InvalidCallbackUsage(MOI.UserCutCallback(), cb))
    end
    ind = Cint[_info(model, var).column - 1 for var in variables]
    CPLEX.cbpostheursoln(
        cb.callback_data, Cint(1), ind, values, NaN, CPXCALLBACKSOLUTION_SOLVE
    )
    return MOI.HEURISTIC_SOLUTION_UNKNOWN
end


text:
struct AvailabilityException <: Exception
    id::Union{Int, String}
    message::String
end

Base.showerror(io::IO, e::AvailabilityException) =
    print(io, "id=\"$(e.id)\", message=\"$(e.message)\" happened")

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
using FakePackage, ReTest
using ReTest: process_args, and
import ReTest: Test
include("../../setup.jl")

retest(FakePackage, recursive=false)
@test FakePackage.RUN == [1, 2]

FakePackage.runtests(r"begin-end")
@test FakePackage.RUN == [1, 2, 1]
FakePackage.runtests(r" for")
@test FakePackage.RUN == [1, 2, 1, 2]

retest(dry=true) # only to update InlineTest.TESTED_MODULES
@test FakePackage.Sub1 in ReTest.TESTED_MODULES
@test FakePackage.Sub2.SubSub in ReTest.TESTED_MODULES

retest(FakePackage, recursive=true)
@test FakePackage.RUN == [1, 2, 1, 2, 1, 2, 3, 4]
retest(FakePackage) # equivalent to recursive=true
@test FakePackage.RUN == [1, 2, 1, 2, 1, 2, 3, 4, 1, 2, 3, 4]

Test.@testset "test macros" begin
    empty!(FakePackage.FakeMacro.RUN)
    retest(FakePackage, "test macro", verbose=9)
    @test FakePackage.FakeMacro.RUN == [4, 2, 3]
end

# in the following, in FakePackage (not FakePackage.FakeMacro), the line `x = 3` is
# shown because statically, we don't know that 4!=2, which is the condition for not
# running the same macro on 3
check(FakePackage, "test macro", verbose=9, id=false, dry=true, [], output="""
FakePackage
  test macro
    inner
      macro begin x = 4
      macro for 1
        plus
      macro begin x = 3
      macro for 1
        plus

FakePackage.FakeMacro
  test macro
    macro begin x = 2
    macro for 1
      plus
    macro begin x = 3
    macro for 1
      plus
""")

module Included ############################################
using ReTest
RUN = []

@testset "include parent" begin
    include("included.jl")
    include(joinpath(@__DIR__, "included.jl"))
    if VERSION >= v"1.5"
        include(identity, "included.jl")
        include(identity, joinpath(@__DIR__, "included.jl"))
    end
end
end # Included ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

retest(Included)
@test Included.RUN == (VERSION >= v"1.5" ? [0, 0, 0, 0] : [0, 0])

ReTest.Test.@testset "retest: load" begin
    @test process_args(()).modules == [
        FakePackage => and(),
        FakePackage.FakeMacro => and(),
        FakePackage.Sub1 => and(),
        FakePackage.Sub2.SubSub => and(),
        Main.Included => and(),
    ]
    check([], dry=true, verbose=0, output = """
FakePackage
FakePackage.FakeMacro
FakePackage.Sub1
FakePackage.Sub2.SubSub
Main.Included
""")

    check([], dry=true, verbose=0, load=true, output = """
FakePackage
FakePackage.FakeMacro
FakePackage.Sub1
FakePackage.Sub2.SubSub
Main.Included
Main.FakePackageTests
Main.FakePackageTests.Sub
Main.FakePackageTests.Sub.SubSub
""")

    check(FakePackage, [], dry=true, verbose=0, load=true, recursive=false, output = """
FakePackage
Main.FakePackageTests
""")

    @test process_args((FakePackage, 2)).modules ==
        [
            FakePackage => and(2),
            FakePackage.FakeMacro => and(2),
            FakePackage.Sub1 => and(2),
            FakePackage.Sub2.SubSub => and(2),
        ]
    @test process_args((FakePackage, 2), load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(2),
        ]
    # cumulative patterns
    @test process_args((FakePackage, 2, FakePackageTests => 3),
                       load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3, 2),
        ]
    # specified pattern for test module
    @test process_args((FakePackage => 2, FakePackageTests => 3),
                       load=true, recursive=false).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3),
        ]
    @test process_args((FakePackage => 2, FakePackageTests => 3),
                       load=true, recursive=true).modules ==
        [
            FakePackage => and(2),
            FakePackageTests => and(3),
            FakePackage.FakeMacro => and(2),
            FakePackage.Sub1 => and(2),
            FakePackage.Sub2.SubSub => and(2),
            FakePackageTests.Sub => and(3),
            FakePackageTests.Sub.SubSub => and(3),
        ]
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
#
# Run:
#   julia hello.jl
#
using UnicodePlots
using Distributions
f() = rand(Normal())

println("hello, world")
@show f()
println(lineplot(1:10, (1:10).^2))

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


text:
#Define useful functions to handle parallelization

#Return the number of effective workers
#If more than 1 worker is present Julia will use nprocs()-1 workers
function NumberWorkers()
    if nprocs() >= 3
        return nprocs() - 1
    elseif nprocs() <3
        return 1
    else error("Number of workers must be a positive integer.")
    end
end

#Divide the total number of sample in nprocs() equal parts
#If more than 1 worker is present Julia will use nprocs()-1 workers
function partition_sample(nsample::Int64)
    if nprocs() == 1
        return nsample
    else
        return Int(round(nsample/(nprocs()-1)))
    end
end

#Divide number of entris of S-matrix into nprocs()-1 parts
#and sort them into a given array strip
#dim=number of variational parameters
function GenerateStrips(strip::Array{Array{Int64}}, dim::Int64)

    if nprocs() == 3
        dd = Int(round(dim/2))
        push!(strip,collect(1:dd))
        push!(strip,collect(dd+1:dim))
        return
    end

    dd = Int(round(dim/(2*(nprocs()-1))))

    for i=1:2
      push!(strip,collect(1+(i-1)*dd:i*dd))
    end

    for i=3:nprocs()-2
      push!(strip,collect(1+2*dd+(i-3)*2*dd:(i-2)*2*dd + 2*dd))
    end

    push!(strip,collect(1+strip[end][end]:dim))

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
#####################
# helper arithmetic to fuse complex high-arity operations
#####################

function __init__()
    set_zero_subnormals(true) # this is supposed to speed up floating point arithmetic on certain architectures
end

@inline safe_div(x::T, y::T) where T = iszero(y) ? 0 : x/y # seems to be much slower than replacing NaN after the fact

@inline replace_nan(x::T) where T = isnan(x) ? zero(T) : x

# need to define this because Julia says true + true = 2 instead of true
Base.@pure accumulator_op(::Type{<:Number}) = :(@fastmath +)
Base.@pure accumulator_op(::Type{Bool}) = :(|)
Base.@pure product_op(::Type{<:Number}) = :(@fastmath *)
Base.@pure product_op(::Type{Bool}) = :(&)

@inline assign(acc::AbstractArray, pr::AbstractArray) =
    acc .= pr # seems to be just as fast as copyto!

@inline @generated accumulate_val(acc::AbstractArray, x::AbstractArray) =
    :(@fastmath  acc .= $(accumulator_op(eltype(acc))).(acc, x))

# generate product expressions

const max_unroll_products = 15

"generate an expression for the loop-unrolled product of xs..."
function expand_product(num_factors, ::Type{E}, xs::Union{Symbol,Expr})  where E
    if num_factors<=0
        :()
    elseif num_factors==1
        :($xs[1])
    else
        :($(product_op(E)).($(expand_product(num_factors-1,E,xs)),$xs[$num_factors])) #strangely adding @inbounds here slows things down a lot
    end
end

"generate an expression for the loop-unrolled product of x1 and xs..."
function expand_product(num_products, ::Type{E}, x1::Union{Symbol,Expr}, xs::Union{Symbol,Expr}) where E
    if num_products == 0
        x1
    else
        @assert num_products <= max_unroll_products "Unrolling loop too many times ($num_products), the compiler will choke... implement a loopy version instead."
        #strangely adding @inbounds here slows things down a lot
        :($(product_op(E)).($(expand_product(num_products-1,E,x1,xs)),xs[$num_products]))
    end
end

# number array
const NArr = AbstractArray{<:Number}

# computer a product efficiently with almost no allocations

@inline @generated function prod_fast_unroll(x1::NArr, xs::NArr...)
    :(@fastmath $(expand_product(length(xs),eltype(x1),:x1,:xs)))
end

@inline function prod_fast(xs::AbstractArray{<:NArr})
    if length(xs) > max_unroll_products
        tmp = prod_fast(xs[max_unroll_products:end])
        prod_fast_unroll(tmp, xs[1:max_unroll_products-1]...)
    else
        prod_fast_unroll(xs...)
    end
end

@inline function prod_fast(x1::NArr, xs::AbstractArray{<:NArr})
    if length(xs) > max_unroll_products-1
        tmp = prod_fast(xs[max_unroll_products-1:end])
        prod_fast_unroll(x1, tmp, xs[1:max_unroll_products-2]...)
    else
        prod_fast_unroll(x1, xs...)
    end
end

# assign a product to an accumulator argument

@inline @generated function assign_prod_unroll(acc::NArr, x1::NArr, xs::NArr...)
    :(acc .= $(expand_product(length(xs),eltype(acc),:x1,:xs)))
end

@inline function assign_prod(acc::NArr, xs::AbstractVector{<:NArr})
    if length(xs) > max_unroll_products
        tmp = prod_fast(xs[max_unroll_products:end])
        assign_prod_unroll(acc, tmp, xs[1:max_unroll_products-1]...)
    else
        assign_prod_unroll(acc, xs...)
    end
end

@inline function assign_prod(acc::NArr, x1::NArr, xs::AbstractVector{<:NArr})
    if length(xs) > max_unroll_products-1
        tmp = prod_fast(xs[max_unroll_products-1:end])
        assign_prod_unroll(acc, tmp, x1, xs[1:max_unroll_products-2]...)
    else
        assign_prod_unroll(acc, x1, xs...)
    end
end

# accumulate a product into an accumulator argument

# Many implementations of `accumulate_prod` are very slow, and somehow don't properly fuse broadcasts':
# - acc .|= (&).(xs...) is terrible!
# - acc .|= unrolled_reduce((x,y) -> x .& y, x1, xs) is also bad, it still does memory allocations!
# Hence, meta-programming to the rescue.
@inline @generated function accumulate_prod_unroll(acc::NArr, x1::NArr, xs::NArr...)
    :(@fastmath acc .= $(accumulator_op(eltype(acc))).(acc, $(expand_product(length(xs),eltype(x1),:x1,:xs))))
end

# specialize for when there is a constant x0 involved, so that x1, xs... can still be multiplied as Bools. Make sure multiplication priorities are set correctly!
@inline @generated function accumulate_prod_unroll(acc::NArr, x0::Real, x1::NArr, xs::NArr...)
    :(@fastmath acc .= $(accumulator_op(eltype(acc))).(acc, x0 .* ($(expand_product(length(xs),eltype(x1),:x1,:xs)))))
end

@inline function accumulate_prod(acc::NArr, xs::AbstractVector{<:NArr})
    if length(xs) > max_unroll_products
        tmp = prod_fast(xs[max_unroll_products:end])
        accumulate_prod_unroll(acc, tmp, xs[1:max_unroll_products-1]...)
    else
        accumulate_prod_unroll(acc, xs...)
    end
end

@inline function accumulate_prod(acc::NArr, x1::Real, xs::NArr...)
    if length(xs) > max_unroll_products-1
        tmp = prod_fast(collect(xs)[max_unroll_products-1:end])
        accumulate_prod_unroll(acc, x1, tmp, collect(xs)[1:max_unroll_products-2]...)
    else
        accumulate_prod_unroll(acc, x1, xs...)
    end
end

# assign a product into an accumulator argument and divide by a normalizer `z`

@inline @generated function assign_prod_normalized_unroll(acc::NArr, z::NArr,
    x1::NArr, xs::NArr...)
    x1ze = (eltype(z) <: Bool) ? :x1 : :(replace_nan.(x1 ./ z)) # Bool z is always equal to 0 or 1, and hence does nothing useful for flows (assuming z > x1)
    # use fastmath for all but the division where we require proper NaN handling
    :(acc .= $(expand_product(length(xs),eltype(acc),x1ze,:xs)))
end

@inline function assign_prod_normalized(acc::NArr, z::NArr,
x1::NArr, xs::AbstractArray{<:NArr})
    if length(xs) > max_unroll_products
        tmp = prod_fast(xs[max_unroll_products:end])
        assign_prod_normalized_unroll(acc, z , x1, tmp, xs[1:max_unroll_products-1]...)
    else
        assign_prod_normalized_unroll(acc, z , x1, xs...)
    end
end

# accumulate a product into an accumulator argument and divide by a normalizer `z`

@inline @generated function accumulate_prod_normalized_unroll(acc::NArr, z::NArr,
                                                        x1::NArr, xs::NArr...)
    x1ze = (eltype(z) <: Bool) ? :x1 : :(replace_nan.(x1 ./ z)) # Bool z is always equal to 0 or 1, and hence does nothing useful for flows (assuming z > x1)
    # use fastmath for all but the division where we require proper NaN handling
    :(acc .= $(accumulator_op(eltype(acc))).(acc, $(expand_product(length(xs),eltype(acc),x1ze,:xs))))
end

@inline function accumulate_prod_normalized(acc::NArr, z::NArr,
                                        x1::NArr, xs::AbstractArray{<:NArr})
    if length(xs) > max_unroll_products
        tmp = prod_fast(xs[max_unroll_products:end])
        accumulate_prod_normalized_unroll(acc, z , x1, tmp, xs[1:max_unroll_products-1]...)
    else
        accumulate_prod_normalized_unroll(acc, z , x1, xs...)
    end
end

# Specialized version of sum_of_products for each 2-layer subcircuit
# Adds a lot of compile time but no noticable speedup over the 1-layer operators above
# test: spe3(Float16,Tuple{NTuple{3,Float16},NTuple{2,Float16},NTuple{2,Float16}})
@generated function sum_of_products(y, ps::Tuple)
    sum = expand_product(length(ps.parameters[1].parameters), Float32, :(ps[1]))
    for i in 2:length(ps.parameters)
        prod = expand_product(length(ps.parameters[i].parameters), Float32, :(ps[$i]))
        sum = :($(accumulator_op(eltype(y))).($sum, $prod))
    end
    # println(:(y .= $sum))
    :(y .= $sum)
end

#TODO: make this method rebust to very large xs sets as above
@inline @generated function count_conjunction(x1::BitVector, xs::BitVector...)::UInt32
    :(count($(expand_product(length(xs),eltype(x1),:x1,:xs))))
end

#TODO: make this method rebust to very large xs sets as above
# @inline @generated function sum_weighted_product(weights::NArr, x1::NArr, xs::NArr...)
#     :(sum($(expand_product(length(xs),eltype(x1),:x1,:xs)) .* weights))
# end
@inline @generated function sum_weighted_product(weights::NArr, x1::NArr, xs::NArr...)
    :(sum(weights[$(expand_product(length(xs),eltype(x1),:x1,:xs))]))
end


text:
# ## Basic commands
#
# This is a very brief and rough primer if you're new to Julia and wondering how to do simpe things that are relevant for data analysis.
#
# Defining a vector

x = [1, 3, 2, 5]
@show x
@show length(x)

# Operations between vectors

y = [4, 5, 6, 1]
z = x .+ y # elementwise operation

# Defining a matrix

X = [1  2; 3 4]

# You can also do that from a vector

X = reshape([1, 2, 3, 4], 2, 2)

# But you have to be careful that it fills the matrix by column; so if you want to get the same result as before let's just permute the dimensions

X = permutedims(reshape([1, 2, 3, 4], 2, 2))

# Function calls can be split with the `|>` operator so that the above can also be written

X = reshape([1,2,3,4], 2, 2) |> permutedims

# You don't have to do that of course but we will sometimes use it in these tutorials.
#
# There's a wealth of functions available for simple math operations

x = 4
@show x^2
@show sqrt(x)

# Element wise operations on a collection can be done with the dot syntax:

sqrt.([4, 9, 16])

# The packages `Statistics` and `StatsBase` offer a number of useful function for stats:

using Statistics, StatsBase

# **Note**: if you don't have `StatsBase`, you can add it using `using Pkg; Pkg.add("StatsBase")`.
# Right, let's now compute some simple statistics:

x = randn(1_000) # 500 points iid from a N(0, 1)
μ = mean(x)
σ = std(x)
@show (μ, σ)

# Indexing data starts at 1, use `:` to indicate the full range

X = [1 2; 3 4; 5 6]
@show X[1, 2]
@show X[:, 1]
@show X[1, :]
@show X[[1, 2], [1, 2]]

# `size` gives dimensions (nrows, ncolumns)

size(X)

# ## Loading data
#
# There are many ways to load data in Julia, one convenient one is via the CSV package.

using CSV

# Many datasets are available via the `RDatasets` package

using RDatasets

# And finally the `DataFrames` package allows to manipulate data easily

using DataFrames

# Let's load some data from RDatasets (the full list of datasets is available [here](http://vincentarelbundock.github.io/Rdatasets/datasets.html))

auto = dataset("ISLR", "Auto")
first(auto, 3)

# The `describe` function allows to get an idea for the data:

describe(auto, :mean, :median, :std)

# To retrieve column names, you can use `names`:

names(auto)

# Accesssing columns can be done in different ways:

mpg = auto.MPG
mpg = auto[:, 1]
mpg = auto[:, :MPG]
mpg |> mean

# To get dimensions you can use `size` and `nrow` and `ncol`

@show size(auto)
@show nrow(auto)
@show ncol(auto)

# For more detailed tutorials on basic data wrangling in Julia, consider
#
# * the [learn x in y](https://learnxinyminutes.com/docs/julia/) julia tutorial
# * the [`DataFrames.jl` docs](http://juliadata.github.io/DataFrames.jl/latest/)
# * the [`StatsBases.jl` docs](https://juliastats.org/StatsBase.jl/latest/)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using PoreMatMod
import Aqua

# ambiguity testing finds many "problems" outside the scope of this package
ambiguities=false

# to skip when checking for stale dependencies and missing compat entries
#   Aqua is added in a separate CI job, so (ironically) does not work w/ itself
stale_deps = (ignore=[:Aqua],)

Aqua.test_all(PoreMatMod;
    ambiguities=ambiguities,
    stale_deps=stale_deps
)


text:
#transform lhs ~ x + rhs to lhs + x ~ rhs
function _shift(@nospecialize(formula::FormulaTerm))
    lhs = formula.lhs
    rhs = formula.rhs
    if !(lhs isa Tuple)
        lhs = (lhs,)
    end
    if !(rhs isa Tuple)
        rhs = (rhs,)
    end
    i = findfirst(x -> x isa Term, rhs)
    FormulaTerm(tuple(lhs..., rhs[i]), Tuple(term for term in rhs if term != rhs[i]))
end

@recipe function f(::Type{Val{:linearfit}}, x, y, z)
    seriestype := :scatter
    @series begin
        x := x
        y := y
        seriestype := :scatter
        ()
    end
    X = hcat(ones(length(x)), x)
    yhat = X * (X'X \ X'y)
    @series begin
        seriestype := :path
        label := ""
        primary := false
        x := x
        y := yhat
        ()
    end
    primary := false
    ()
end


text:
using Polyhedra
using Documenter, Literate

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

const EXAMPLES = [
    "Convex hull and intersection",
    "Extended Formulation",
    "Minimal Robust Positively Invariant Set",
    "Convex hull of a set of points",
]

for example in EXAMPLES
    example_filepath = joinpath(EXAMPLES_DIR, example * ".jl")
    Literate.markdown(example_filepath, OUTPUT_DIR)
    Literate.notebook(example_filepath, OUTPUT_DIR)
end

# See https://juliadocs.github.io/Documenter.jl/v0.25/man/doctests/#Setup-Code
DocMeta.setdocmeta!(Polyhedra, :DocTestSetup, :(using Polyhedra); recursive=true)

makedocs(
    sitename = "Polyhedra",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/jump-dev/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Index" => "index.md",
        "Installation" => "installation.md",
        "Representation" => "representation.md",
        "Polyhedron" => "polyhedron.md",
        "Plot" => "plot.md",
        "Containment/Redundancy" => "redundancy.md",
        "Projection/Elimination" => "projection.md",
        "Optimization" => "optimization.md",
        "Utilities" => "utilities.md",
        "Examples" => map(EXAMPLES) do example
            example => "generated/$example.md"
        end
    ],
    # The following ensures that we only include the docstrings from
    # this module for functions define in Base that we overwrite.
    modules = [Polyhedra]
)

deploydocs(
    repo   = "github.com/JuliaPolyhedra/Polyhedra.jl.git",
)


text:
# ---
# title: 1451. Rearrange Words in a Sentence
# id: problem1451
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Sort
# link: <https://leetcode.com/problems/rearrange-words-in-a-sentence/description/>
# hidden: true
# ---
#
# Given a sentence `text` (A  _sentence_  is a string of space-separated words)
# in the following format:
#
#   * First letter is in upper case.
#   * Each word in `text` are separated by a single space.
#
# Your task is to rearrange the words in text such that all words are rearranged
# in an increasing order of their lengths. If two words have the same length,
# arrange them in their original order.
#
# Return the new text following the format shown above.
#
#
#
# **Example 1:**
#
#
#
#     Input: text = "Leetcode is cool"
#     Output: "Is cool leetcode"
#     Explanation: There are 3 words, "Leetcode" of length 8, "is" of length 2 and "cool" of length 4.
#     Output is ordered by length and the new first word starts with capital letter.
#
#
# **Example 2:**
#
#
#
#     Input: text = "Keep calm and code on"
#     Output: "On and keep calm code"
#     Explanation: Output is ordered as follows:
#     "On" 2 letters.
#     "and" 3 letters.
#     "keep" 4 letters in case of tie order by position in original text.
#     "calm" 4 letters.
#     "code" 4 letters.
#
#
# **Example 3:**
#
#
#
#     Input: text = "To be or not to be"
#     Output: "To be or to be not"
#
#
#
#
# **Constraints:**
#
#   * `text` begins with a capital letter and then contains lowercase letters and single space between words.
#   * `1 <= text.length <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
using Base.Test
using PortfolioModels

my_tests = ["preprocess.jl",
            "solve.jl",
            "build.jl"]

print_with_color(:cyan, "Running tests: ")
println("")

for my_test in my_tests
    print_with_color(:magenta, "**   ")
    print_with_color(:blue, "$my_test")
    println("")
    include(my_test)
end


text:
function do_max_pol_b0(method, model_key)

    # ep params
    me_params = lglob(iJR, :maxent, :params)
    @extract me_params: alpha epsconv maxiter
    @extract me_params: damp maxvar minvar

    # handle cache
    is_cached(;method) && return

    # setup
    model = iJR.load_model(model_key)

    lock(WLOCK) do
        @info("Doing... ",
            method, threadid()
        ); println()
    end

    # maxent
    epout = ChEP.maxent_ep(model;
        alpha, damp, maxiter,
        epsconv, maxvar, minvar,
    )

    # storing
    lock(WLOCK) do

        # Storing
        dat = Dict()
        dat[:exp_beta] = 0.0
        dat[:epouts] = Dict(0.0 => epout)
        dat[:model] = model |> ChU.compressed_model

        # caching
        datfile = dat_file(;method)
        sdat(iJR, dat, datfile)

        @info("Finished ", method); println()
    end
end


text:
# Test code for JULIA power files

using PowerModels
using Ipopt
using JuMP
using JSON
# using .PP2PM

function run_powermodels(json_path)
    # pm = PP2PM.load_pm_from_json(json_path)
    ipopt_solver = JuMP.with_optimizer(Ipopt.Optimizer, print_level=5)
    # data = PowerModels.parse_file("../../casefiles/case_SyntheticUSA.m")
    # data = PowerModels.parse_file("../../casefiles/case_ACTIVSg2000.m")
    data = PowerModels.parse_file("../../casefiles/case9target.m")
    net = JSON.parse(JSON.json(data))
    result = PowerModels.run_ac_opf(net, ipopt_solver)

    print(json(data["bus"], 4))
    #print(keys(data))
    #result = run_ac_opf(data, with_optimizer(Ipopt.Optimizer))
    #result = PowerModels.run_ac_opf(pm, ipopt_solver,
    #                                setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

run_powermodels(false)


text:
function MOI.delete(model::Optimizer, vi::MOI.VariableIndex)
    check_inbounds(model, vi)
    if vi in model.delvars
        throw(MOI.InvalidIndex(vi))
    end
    rhp_delete_var(model.ctx, vi.value-1)
    push!(model.delvars, vi)

end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{<: Union{SF, SF}, <: LS})
    rhp_delete_eqn(model.ctx, ci.value-1)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{VAF, <: VLS})
    if haskey(model.vaf_mapping, ci)
        for eidx in model.vaf_mapping[ci]
            rhp_delete_eqn(model.ctx, eidx)
        end
    else
        error("Unknown contraint $ci")
    end
    return
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}})
    ctx_setvarlb(model.ctx, ci.value-1, Inf64)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}})
    ctx_setvarub(model.ctx, ci.value-1, -Inf64)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.SingleVariable, Union{MOI.EqualTo{Float64}, MOI.Interval{Float64}}})
    ctx_setvarub(model.ctx, ci.value-1, Inf64)
    ctx_setvarlb(model.ctx, ci.value-1, -Inf64)
end


text:
module PlotsTests

# # don't let pyplot use a gui... it'll crash
# # note: Agg will set gui -> :none in PyPlot
# ENV["MPLBACKEND"] = "Agg"
# try
#     @eval import PyPlot
# catch err
#     warn("Couldn't import PyPlot: $err")
# end

# using Plots
# using FactCheck

# # note: wrap first include in a try block because of the ImageMagick init_deps bug
# try
#     include("imgcomp.jl")
# end
include("imgcomp.jl")

# don't actually show the plots
srand(1234)
default(show=false)

# note: we wrap in a try block so that the tests only run if we have the backend installed
# try
  # Pkg.installed("Gadfly")
  # gadfly()
  # backend()

img_eps = 5e-2

facts("Gadfly") do
    @fact gadfly() --> Plots.GadflyPackage()
    @fact backend() --> Plots.GadflyPackage()

    @fact typeof(plot(1:10)) --> Plots.Plot{Plots.GadflyPackage}

    # plot(x::AVec, y::AVec; kw...)              # one line (will assert length(x) == length(y))
    @fact plot(Int[1,2,3], rand(3)) --> not(nothing)
    # @fact_throws plot(1:5, 1:4)

    # plot(x::AVec, y::AMat; kw...)              # multiple lines (one per column of x), all sharing x (will assert length(x) == size(y,1))
    @fact plot(sort(rand(10)), rand(Int, 10, 3)) --> not(nothing)
    # @fact_throws(plot!(rand(10), rand(9,2)))

    # plot(x::AMat, y::AMat; kw...)              # multiple lines (one per column of x/y... will assert size(x) == size(y))
    @fact plot!(rand(10,3), rand(10,3)) --> not(nothing)

    image_comparison_facts(:gadfly, skip=[4,6,19,23,24], eps=img_eps)
end

facts("PyPlot") do
    @fact pyplot() --> Plots.PyPlotPackage()
    @fact backend() --> Plots.PyPlotPackage()
    image_comparison_facts(:pyplot, skip=[10,19,21,23], eps=img_eps)
end

# catch err
#     warn("Skipped Gadfly due to: ", string(err))
# end

# # note: we wrap in a try block so that the tests only run if we have the backend installed
# try
#   Pkg.installed("Qwt")
#   qwt()
#   backend()
#   facts("Qwt") do
#     @fact backend(:qwt) --> Plots.QwtPackage()
#     @fact backend() --> Plots.QwtPackage()
#     @fact typeof(plot(1:10)) --> Plots.Plot{Plots.QwtPackage}

#     # plot(y::AVec; kw...)                       # one line... x = 1:length(y)
#     @fact plot(1:10) --> not(nothing)
#     @fact length(current().o.lines) --> 1

#     # plot(x::AVec, f::Function; kw...)          # one line, y = f(x)
#     @fact plot(1:10, sin) --> not(nothing)
#     @fact current().o.lines[1].y --> sin(collect(1:10))

#     # plot(x::AMat, f::Function; kw...)          # multiple lines, yᵢⱼ = f(xᵢⱼ)
#     @fact plot(rand(10,2), sin) --> not(nothing)
#     @fact length(current().o.lines) --> 2

#     # plot(y::AMat; kw...)                       # multiple lines (one per column of x), all sharing x = 1:size(y,1)
#     @fact plot!(rand(10,2)) --> not(nothing)
#     @fact length(current().o.lines) --> 4

#     # plot(x::AVec, fs::AVec{Function}; kw...)   # multiple lines, yᵢⱼ = fⱼ(xᵢ)
#     @fact plot(1:10, Function[sin,cos]) --> not(nothing)
#     @fact current().o.lines[1].y --> sin(collect(1:10))
#     @fact current().o.lines[2].y --> cos(collect(1:10))

#     # plot(y::AVec{AVec}; kw...)                 # multiple lines, each with x = 1:length(y[i])
#     @fact plot([11:20 ; rand(10)]) --> not(nothing)
#     @fact current().o.lines[1].x[4] --> 4
#     @fact current().o.lines[1].y[4] --> 14
#   end
# catch err
#     warn("Skipped Qwt due to: ", string(err))
# end

# try
    # Pkg.installed("PyPlot")
    # pyplot()
    # backend()
    # facts("PyPlot") do
    #     @fact backend(:pyplot) --> Plots.PyPlotPackage()
    #     @fact backend() --> Plots.PyPlotPackage()
    #     @fact typeof(plot(1:10)) --> Plots.Plot{Plots.PyPlotPackage}

    #     # image_comparison_facts(:pyplot, skip=[19])
    # end
# catch err
#     warn("Skipped PyPlot due to: ", string(err))
# end

# try
#     Pkg.installed("UnicodePlots")
#     unicodeplots()
#     backend()
#     facts("UnicodePlots") do
#         @fact backend(:unicodeplots) --> Plots.UnicodePlotsPackage()
#         @fact backend() --> Plots.UnicodePlotsPackage()
#         @fact typeof(plot(1:10)) --> Plots.Plot{Plots.UnicodePlotsPackage}
#     end
# catch err
#     warn("Skipped UnicodePlots due to: ", string(err))
# end

# try
#     Pkg.installed("Winston")
#     winston()
#     backend()
#     facts("Winston") do
#         @fact backend(:winston) --> Plots.WinstonPackage()
#         @fact backend() --> Plots.WinstonPackage()
#         @fact typeof(plot(1:10)) --> Plots.Plot{Plots.WinstonPackage}
#     end
# catch err
#     warn("Skipped Winston due to: ", string(err))
# end

FactCheck.exitstatus()
end # module


text:
module SpringCollab2020RainbowSpinnerColorAreaController

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/RainbowSpinnerColorAreaController" RainbowSpinnerColorAreaController(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    colors::String="89E5AE,88E0E0,87A9DD,9887DB,D088E2", gradientSize::Number=280.0)

const placements = Ahorn.PlacementDict(
    "Rainbow Spinner Colour Area Controller (Spring Collab 2020)" => Ahorn.EntityPlacement(
        RainbowSpinnerColorAreaController,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::RainbowSpinnerColorAreaController) = 8, 8
Ahorn.resizable(entity::RainbowSpinnerColorAreaController) = true, true

Ahorn.selection(entity::RainbowSpinnerColorAreaController) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RainbowSpinnerColorAreaController, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.4, 0.4, 1.0, 0.4), (0.4, 0.4, 1.0, 1.0))
end

end


text:
@testset "112.path-sum.jl" begin
    @test has_path_sum(
        TreeNode{Int}([5, 4, 8, 11, nothing, 13, 4, 7, 2, nothing, nothing, nothing, 1]), 22
    ) == true
    @test has_path_sum(TreeNode{Int}([1, 2]), 1) == false
    @test has_path_sum(nothing, 0) == false
end


text:
"""
    abstract type MetabolicModel end

A helper supertype that wraps everything usable as a linear-like model for
COBREXA functions.

If you want your model type to work with COBREXA, add the `MetabolicModel` as
its supertype, and implement the accessor functions. Accessors
[`reactions`](@ref), [`metabolites`](@ref), [`stoichiometry`](@ref),
[`bounds`](@ref) and [`objective`](@ref) must be implemented; others are not
mandatory and default to safe "empty" values.
"""
abstract type MetabolicModel end

const SparseMat = SparseMatrixCSC{Float64,Int}
const SparseVec = SparseVector{Float64,Int}
const MatType = AbstractMatrix{Float64}
const VecType = AbstractVector{Float64}
const StringVecType = AbstractVector{String}

"""
    GeneAssociation = Vector{Vector{String}}

An association to genes, represented as a logical formula in a positive
disjunctive normal form (DNF). (The 2nd-level vectors of strings are connected
by "and" to form conjunctions, and the 1st-level vectors of these conjunctions
are connected by "or" to form the DNF.)
"""
const GeneAssociation = Vector{Vector{String}}

"""
    MetaboliteFormula = Dict{String,Int}

Dictionary of atoms and their abundances in a molecule.
"""
const MetaboliteFormula = Dict{String,Int}

"""
    Annotations = Dict{String,Vector{String}}

Dictionary used to store (possible multiple) standardized annotations of
something, such as a [`Metabolite`](@ref) and a [`Reaction`](@ref).

# Example
```
Annotations("PubChem" => ["CID12345", "CID54321"])
```
"""
const Annotations = Dict{String,Vector{String}}

"""
    Notes = Dict{String,Vector{String}}

Free-form notes about something (e.g. a [`Gene`](@ref)), categorized by
"topic".
"""
const Notes = Dict{String,Vector{String}}

_missing_impl_error = (m, a) -> throw(MethodError(m, a))

"""
    reactions(a::MetabolicModel)::Vector{String}

Return a vector of reaction identifiers in a model.
"""
function reactions(a::MetabolicModel)::Vector{String}
    _missing_impl_error(reactions, (a,))
end

"""
    metabolites(a::MetabolicModel)::Vector{String}

Return a vector of metabolite identifiers in a model.
"""
function metabolites(a::MetabolicModel)::Vector{String}
    _missing_impl_error(metabolites, (a,))
end

"""
    n_reactions(a::MetabolicModel)::Int

Get the number of reactions in a model.
"""
function n_reactions(a::MetabolicModel)::Int
    length(reactions(a))
end

"""
    n_metabolites(a::MetabolicModel)::Int

Get the number of metabolites in a model.
"""
function n_metabolites(a::MetabolicModel)::Int
    length(metabolites(a))
end

"""
    stoichiometry(a::MetabolicModel)::SparseMat

Get the sparse stoichiometry matrix of a model.
"""
function stoichiometry(a::MetabolicModel)::SparseMat
    _missing_impl_error(stoichiometry, (a,))
end

"""
    bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}

Get the lower and upper flux bounds of a model.
"""
function bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
    _missing_impl_error(bounds, (a,))
end

"""
    balance(a::MetabolicModel)::SparseVec

Get the sparse balance vector of a model (ie. the `b` from `S x = b`).
"""
function balance(a::MetabolicModel)::SparseVec
    return spzeros(n_metabolites(a))
end

"""
    objective(a::MetabolicModel)::SparseVec

Get the objective vector of a model.
"""
function objective(a::MetabolicModel)::SparseVec
    _missing_impl_error(objective, (a,))
end

"""
    coupling(a::MetabolicModel)::SparseMat

Get a matrix of coupling constraint definitions of a model. By default, there
is no coupling in the models.
"""
function coupling(a::MetabolicModel)::SparseMat
    return spzeros(0, n_reactions(a))
end

"""
    n_coupling_constraints(a::MetabolicModel)::Int

Get the number of coupling constraints in a model.
"""
function n_coupling_constraints(a::MetabolicModel)::Int
    size(coupling(a), 1)
end

"""
    coupling_bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}

Get the lower and upper bounds for each coupling bound in a model, as specified
by `coupling`. By default, the model does not have any coupling bounds.
"""
function coupling_bounds(a::MetabolicModel)::Tuple{SparseVec,SparseVec}
    return (spzeros(0), spzeros(0))
end

"""
    genes(a::MetabolicModel)::Vector{String}

Return identifiers of all genes contained in the model. By default, there are
no genes.

In SBML, these are usually called "gene products" but we write `genes` for
simplicity.
"""
function genes(a::MetabolicModel)::Vector{String}
    return []
end

"""
    n_genes(a::MetabolicModel)::Int

Return the number of genes in the model (as returned by [`genes`](@ref)). If
you just need the number of the genes, this may be much more efficient than
calling [`genes`](@ref) and measuring the array.
"""
function n_genes(a::MetabolicModel)::Int
    return length(genes(a))
end

"""
    reaction_gene_association(a::MetabolicModel, gene_id::String)::Maybe{GeneAssociation}

Returns the sets of genes that need to be present so that the reaction can work
(technically, a DNF on gene availability, with positive atoms only).

For simplicity, `nothing` may be returned, meaning that the reaction always
takes place. (in DNF, that would be equivalent to returning `[[]]`.)
"""
function reaction_gene_association(
    a::MetabolicModel,
    reaction_id::String,
)::Maybe{GeneAssociation}
    return nothing
end

"""
    metabolite_formula(
        a::MetabolicModel,
        metabolite_id::String,
    )::Maybe{MetaboliteFormula}

Return the formula of metabolite `metabolite_id` in `model`.
Return `nothing` in case the formula is not known or irrelevant.
"""
function metabolite_formula(
    model::MetabolicModel,
    metabolite_id::String,
)::Maybe{MetaboliteFormula}
    return nothing
end

"""
metabolite_charge(model::MetabolicModel, metabolite_id::String)::Maybe{Int}

Return the charge associated with metabolite `metabolite_id` in `model`.
Returns `nothing` if charge not present.
"""
function metabolite_charge(model::MetabolicModel, metabolite_id::String)::Maybe{Int}
    return nothing
end

"""
    reaction_annotations(a::MetabolicModel, reaction_id::String)::Annotations

Return standardized names that may help identifying the reaction. The
dictionary assigns vectors of possible identifiers to identifier system names,
e.g. `"Reactome" => ["reactomeID123"]`.
"""
function reaction_annotations(a::MetabolicModel, reaction_id::String)::Annotations
    return Dict()
end

"""
    metabolite_annotations(a::MetabolicModel, metabolite_id::String)::Annotations

Return standardized names that may help to reliably identify the metabolite. The
dictionary assigns vectors of possible identifiers to identifier system names,
e.g. `"ChEMBL" => ["123"]` or `"PubChem" => ["CID123", "CID654645645"]`.
"""
function metabolite_annotations(a::MetabolicModel, metabolite_id::String)::Annotations
    return Dict()
end

"""
    gene_annotations(a::MetabolicModel, gene_id::String)::Annotations

Return standardized names that identify the corresponding gene or product. The
dictionary assigns vectors of possible identifiers to identifier system names,
e.g. `"PDB" => ["PROT01"]`.
"""
function gene_annotations(a::MetabolicModel, gene_id::String)::Annotations
    return Dict()
end

"""
    reaction_notes(model::MetabolicModel, reaction_id::String)::Notes

Return the notes associated with reaction `reaction_id` in `model`.
"""
function reaction_notes(model::MetabolicModel, reaction_id::String)::Notes
    return Dict()
end

"""
    metabolite_notes(model::MetabolicModel, metabolite_id::String)::Notes

Return the notes associated with metabolite `reaction_id` in `model`.
"""
function metabolite_notes(model::MetabolicModel, metabolite_id::String)::Notes
    return Dict()
end

"""
    gene_notes(model::MetabolicModel, gene_id::String)::Notes

Return the notes associated with the gene `gene_id` in `model`.
"""
function gene_notes(model::MetabolicModel, gene_id::String)::Notes
    return Dict()
end

"""
    metabolite_compartment(model::MetabolicModel, metabolite_id::String)::Maybe{String}

Return the compartment of metabolite `metabolite_id` in `model` if it is assigned. If not,
return `nothing`.
"""
function metabolite_compartment(model::MetabolicModel, metabolite_id::String)::Maybe{String}
    return nothing
end

"""
    reaction_subsystem(model::MetabolicModel, reaction_id::String)::Maybe{String}

Return the subsystem of reaction `reaction_id` in `model` if it is assigned. If not,
return `nothing`.
"""
function reaction_subsystem(model::MetabolicModel, reaction_id::String)::Maybe{String}
    return nothing
end

"""
    reaction_stoichiometry(model::MetaboliteModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::MetabolicModel, rxn_id::String)::Dict{String,Float64}
    mets = metabolites(m)
    Dict(
        mets[k] => v for (k, v) in
        zip(findnz(stoichiometry(m)[:, first(indexin([rxn_id], reactions(m)))])...)
    )
end

"""
    precache!(a::MetabolicModel)::Nothing

Do whatever is feasible to get the model into a state that can be read from
as-quickly-as-possible. This may include e.g. generating helper index
structures and loading delayed parts of the model from disk. The model should
be modified "transparently" in-place. Analysis functions call this right before
applying modifications or converting the model to the optimization model using
[`make_optimization_model`](@ref); usually on the same machine where the
optimizers (and, generally, the core analysis algorithms) will run. The calls
are done in a good hope that the performance will be improved.

By default, it should be safe to do nothing.
"""
function precache!(a::MetabolicModel)::Nothing
    nothing
end