text:
using BinaryBuilder

include("../../../fancy_toys.jl")

name = "CUDA_full"
version = v"11.0.3"

sources_linux = [
    FileSource("https://developer.download.nvidia.com/compute/cuda/11.0.3/local_installers/cuda_11.0.3_450.51.06_linux.run",
               "b079c4e408adf88c3f1ffb8418a97dc4227c37935676b4bf4ca0beec6c328cc0", "installer.run")
]
sources_linux_ppc64le = [
    FileSource("https://developer.download.nvidia.com/compute/cuda/11.0.3/local_installers/cuda_11.0.3_450.51.06_linux_ppc64le.run",
               "4775b21df004b1433bafff9b48a324075c008509f4c0fe28cd060d042d2e0794", "installer.run")
]
sources_win10 = [
    FileSource("http://developer.download.nvidia.com/compute/cuda/11.0.3/local_installers/cuda_11.0.3_451.82_win10.exe",
               "a639e0c097717c5f544f7a81b7162b4cb49243a30fd892f3267bc5532c8e2584", "installer.exe")
]

script = raw"""
cd ${WORKSPACE}/srcdir

# use a temporary directory to avoid running out of tmpfs
temp=${WORKSPACE}/tmpdir
mkdir ${temp}

mkdir ${prefix}/cuda
if [[ ${target} == *-linux-gnu ]]; then
    # don't run the embedded installer (libc issue, doesn't work on ppc),
    # but instruct makeself to uncompress the embedded archive.
    sh installer.run --noexec --target "${temp}"
    cd ${temp}/builds
    find .

    cp cuda_documentation/EULA.txt ${prefix}/cuda

    for project in cuda_cudart cuda_cuobjdump cuda_cupti cuda_gdb cuda_memcheck \
                   cuda_nvcc cuda_nvdisasm cuda_nvml_dev cuda_nvprof cuda_nvprune \
                   cuda_nvrtc cuda_nvtx cuda_sanitizer_api \
                   libcublas libcufft libcurand libcusolver libcusparse \
                   libnpp libnvjpeg; do
        [[ -d ${project} ]] || { echo "${project} does not exist!"; exit 1; }
        cp -a ${project}/* ${prefix}/cuda
    done

    cp -a integration/Sanitizer/* ${prefix}/cuda/bin
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    apk add p7zip

    7z x installer.exe -o${temp}
    cd ${temp}
    find .

    mv cuda_documentation/Doc/EULA.txt ${prefix}/cuda

    for project in cuda_cudart cuda_cuobjdump cuda_cupti cuda_memcheck \
                   cuda_nvcc cuda_nvdisasm cuda_nvml_dev cuda_nvprof cuda_nvprune \
                   cuda_nvrtc cuda_nvtx cuda_sanitizer_api \
                   libcublas libcufft libcurand libcusolver libcusparse  \
                   libnpp libnvjpeg; do
        [[ -d ${project} ]] || { echo "${project} does not exist!"; exit 1; }
        cp -a ${project}/*/* ${prefix}/cuda
    done

    # NVIDIA Tools Extension Library
    7z x "nsight_nvtx/nsight_nvtx/NVIDIA NVTX Installer.x86_64".*.msi -o${temp}/nvtx_installer
    find nvtx_installer
    for file in nvtx_installer/*.*_*; do
        mv $file $(echo $file | sed 's/\.\(\w*\)_.*/.\1/')
    done
    mv nvtx_installer/*.dll ${prefix}/cuda/bin
    mv nvtx_installer/*64_*.lib ${prefix}/cuda/lib/x64
    mv nvtx_installer/*32_*.lib ${prefix}/cuda/lib/Win32
    mv nvtx_installer/*.h ${prefix}/cuda/include

    # fixup
    chmod +x ${prefix}/cuda/bin/*.exe

    # clean-up
    rm ${prefix}/cuda/*.nvi
fi

rm -rf ${temp}

cd ${prefix}/cuda
install_license EULA.txt
"""

products = Product[
    # this JLL isn't meant for use by Julia packages, but only as build dependency
]

dependencies = []

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux, script,
                   [Linux(:x86_64)], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("powerpc64le-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_ppc64le, script,
                   [Linux(:powerpc64le)], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("x86_64-w64-mingw32")
    build_tarballs(ARGS, name, version, sources_win10, script,
                   [Windows(:x86_64)], products, dependencies;
                   skip_audit=true)
end


text:
using LinearAlgebra

const companion = CommutativeRings.companion

@testset "vector spaces" begin

    Z = ZZ/7

    xnull = VectorSpace(Matrix{Z}(undef, 5, 0))
    xall = VectorSpace(diagm(ones(Z, 5)))
    @test intersect(xnull, xnull) == xnull
    @test intersect(xnull, xall) == xnull
    @test intersect(xall, xall) == xall
    @test sum(xnull, xnull) == xnull
    @test sum(xall, xnull) == xall
    @test sum(xall, xall) == xall

    va = VectorSpace(Z[1 0; 0 1; 0 0; 0 0; 0 0])
    vb = VectorSpace(Z[0 0; 0 0; 0 0; 0 1; 1 0])

    vc = sum(va, vb)
    @test Matrix(vc) == [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 1 0]

    vc = intersect(va, vb)
    @test Matrix(vc) == Matrix{Z}(undef,5,0)

    vc = complement(vb)
    @test Matrix(vc) == Z[0 0 1; 0 1 0; 1 0 0; 0 0 0; 0 0 0]

    wa = VectorSpace(Z[1, 2, 3, 4, 5])
    wb = VectorSpace(Z[1, 2, 4, 4, 0])
    wc = VectorSpace(Z[0, 0, 2, 1, 0])
    xa = sum(wa, wc)
    xb = sum(wb, wc)

    @test intersect(wa, xa) == wa
    @test intersect(xa, wc) == wc
    @test intersect(wc, xb) == wc
    @test intersect(xb, wb) == wb

    xint = intersect(xa, xb)
    @test xint == wc

    xsum = sum(xa, xb)
    @test issubset(wa, xsum)
    @test issubset(wb, xsum)
    @test issubset(wc, xsum)
    @test issubset(xa, xsum)
    @test issubset(xb, xsum)

    @test intersect(xnull, xsum) == xnull
    @test intersect(xall, xsum) == xsum
    @test sum(xnull, xsum) == xsum
    @test sum(xall, xsum) == xall
end

@testset "determinant" begin
    G = ZZ/89
    A = G.([1 2 3; 4 5 6; 7 8 9])
    B = G.([1 2 3; 4 5 6; 7 8 8])
    @test iszero(det(A))
    @test iszero(adjugate(A) * A)
    @test !iszero(det(B))
    @test adjugate(B) ./ det(B) == inv(B)

    p = characteristic_polynomial(B)
    @test iszero(p(B))
    @test characteristic_polynomial(companion(p)) == p
end

ext: jl


text:
using Krotov:
    SquareParametrization,
    TanhParametrization,
    TanhSqParametrization,
    LogisticParametrization,
    LogisticSqParametrization

using Plots
Plots.default(
    linewidth               = 3,
    size                    = (550, 300),
    legend                  = :top,
    foreground_color_legend = nothing,
    background_color_legend = RGBA(1, 1, 1, 0.8),
)

function plot_positive_parametrization_comparison()

    u_vals = collect(range(-3, 3, length=101))
    ϵ_vals = collect(range(0, 1, length=101))
    ϵ_max = 1.0

    pnl1 = plot(
        u_vals,
        abs.(u_vals);
        linestyle=:dash,
        color="black",
        label="",
        xlabel="u",
        ylabel="ϵ",
        legend=false
    )
    plot!(pnl1, u_vals, TanhSqParametrization(ϵ_max).epsilon_of_u.(u_vals), label="TanhSq")
    plot!(
        pnl1,
        u_vals,
        LogisticSqParametrization(ϵ_max).epsilon_of_u.(u_vals),
        label="LogisticSq(k=1)"
    )
    plot!(
        pnl1,
        u_vals,
        LogisticSqParametrization(ϵ_max, k=4.0).epsilon_of_u.(u_vals),
        label="LogisticSq(k=4)"
    )
    plot!(pnl1, u_vals, SquareParametrization().epsilon_of_u.(u_vals), label="Square")
    ylims!(pnl1, (0, 1.2))

    pnl2 = plot(
        ϵ_vals,
        ϵ_vals;
        linestyle=:dash,
        color="black",
        label="",
        xlabel="ϵ",
        ylabel="u"
    )
    plot!(pnl2, ϵ_vals, TanhSqParametrization(ϵ_max).u_of_epsilon.(ϵ_vals), label="TanhSq")
    plot!(
        pnl2,
        ϵ_vals,
        LogisticSqParametrization(ϵ_max).u_of_epsilon.(ϵ_vals),
        label="LogisticSq(k=1)"
    )
    plot!(
        pnl2,
        ϵ_vals,
        LogisticSqParametrization(ϵ_max, k=4.0).u_of_epsilon.(ϵ_vals),
        label="LogisticSq(k=4)"
    )
    plot!(pnl2, ϵ_vals, SquareParametrization().u_of_epsilon.(ϵ_vals), label="Square")
    ylims!(pnl2, (0, 3))

    pnl3 = plot(
        u_vals,
        sign.(u_vals);
        linestyle=:dash,
        color="black",
        label="",
        xlabel="u",
        ylabel="∂ϵ/∂u",
        legend=false
    )
    plot!(
        pnl3,
        u_vals,
        TanhSqParametrization(ϵ_max).de_du_derivative.(u_vals),
        label="TanhSq"
    )
    plot!(
        pnl3,
        u_vals,
        LogisticSqParametrization(ϵ_max).de_du_derivative.(u_vals),
        label="LogisticSq(k=1)"
    )
    plot!(
        pnl3,
        u_vals,
        LogisticSqParametrization(ϵ_max, k=4.0).de_du_derivative.(u_vals),
        label="LogisticSq(k=4)"
    )
    plot!(pnl3, u_vals, SquareParametrization().de_du_derivative.(u_vals), label="Square")
    ylims!(pnl3, (-2, 2))

    plot(
        pnl1,
        pnl2,
        pnl3,
        layout=(1, 3),
        size=(1000, 300),
        left_margin=20Plots.px,
        bottom_margin=20Plots.px
    )

end

if abspath(PROGRAM_FILE) == @__FILE__
    gui(plot_positive_parametrization_comparison())
    readline()
end


text:
reload("src/NNN.jl")

const batch_size  = 20
const input_size  = 2
const output_size = 1
const hidden_size = 40

function toy_data(n)
    x = zeros(Float32, input_size,  n*batch_size)
    x[1,:] = rand(n*batch_size)
    y = zeros(Float32, output_size, n*batch_size)

    for i in 0:batch_size:((n-1)*batch_size)
        a = floor(Int64, rand() * batch_size)+1
        b = floor(Int64, rand() * batch_size)+1
        total = 0.0f0
        for j in 1:batch_size
            pos = i + j
            if j == a || j == b
                x[2, pos] = 1.0f0
                total += x[1, pos]
            end
            y[1,pos] = total
        end
    end

    (x, y)
end

function rnn_train(xtrn, ytrn, xtst, ytst, iter)
    net = NNN.Network(batch_size)
    NNN.push!(net,
             NNN.SimpleRecurrentLayer(batch_size, input_size, hidden_size;
                                     activation=NNN.LeakyRelu()),
             NNN.SGD(;λ=0.0,
                    ε=0.05,
                    gradient_clip=10.0f0),
             NNN.SGD(;λ=0.0,
                    ε=0.05,
                    gradient_clip=10.0f0))
    NNN.push!(net,
             NNN.OutputLayer(batch_size, hidden_size, output_size;
                            activation=NNN.IdentityAct(),
                            loss=NNN.rnn_quadratic_cost),
             NNN.SGD(;λ=0.000001,
                    ε=0.05,
                    gradient_clip=10.0f0))

    @time NNN.train(xtrn, ytrn, net; iter=iter, batch_size=batch_size)
    errors = NNN.test(xtst,  ytst, net; batch_size=batch_size)

    (net, errors)
end

function profile_rnn(xtrn, ytrn, xtst, ytst, iter)
    Profile.clear()
    Profile.init()

    net, errors = @profile rnn_train(xtrn, ytrn, xtst, ytst, iter)

    open("/tmp/prof.dat", "w") do io
        Profile.print(io)
    end

    net
end

xtrn, ytrn = toy_data(20000);
xtst, ytst = toy_data(1000);

net = profile_rnn(xtrn, ytrn, xtst, ytst, 1); "OK"

# NNN.train(xtrn, ytrn, net; iter=10, batch_size=batch_size)


text:
abstract type MessagePassing <: GraphNet end

"""
    message(mp::MessagePassing, x_i, x_j, e_ij)

Message function for the message-passing scheme,
returning the message from node `j` to node `i` .
In the message-passing scheme. the incoming messages
from the neighborhood of `i` will later be aggregated
in order to [`update`](@ref) the features of node `i`.

By default, the function returns `x_j`.
Layers subtyping [`MessagePassing`](@ref) should
specialize this method with custom behavior.

# Arguments

- `mp`: message-passing layer.
- `x_i`: the features of node `i`.
- `x_j`: the features of the nighbor `j` of node `i`.
- `e_ij`: the features of edge (`i`, `j`).

See also [`update`](@ref).
"""
function message end

@inline message(mp::MessagePassing, x_i, x_j, e_ij) = x_j
@inline message(mp::MessagePassing, i::Integer, j::Integer, x_i, x_j, e_ij) = x_j

"""
    update(mp::MessagePassing, m, x)

Update function for the message-passing scheme,
returning a new set of node features `x′` based on old
features `x` and the incoming message from the neighborhood
aggregation `m`.

By default, the function returns `m`.
Layers subtyping [`MessagePassing`](@ref) should
specialize this method with custom behavior.

# Arguments

- `mp`: message-passing layer.
- `m`: the aggregated edge messages from the [`message`](@ref) function.
- `x`: the node features to be updated.

See also [`message`](@ref).
"""
function update end

@inline update(mp::MessagePassing, m, x) = m
@inline update(mp::MessagePassing, i::Integer, m, x) = m

@inline update_edge(mp::MessagePassing, e, vi, vj, u) = GeometricFlux.message(mp, vi, vj, e)
@inline update_vertex(mp::MessagePassing, ē, vi, u) = GeometricFlux.update(mp, ē, vi)


text:
@impl begin
    struct DistMakeFactors
        numpartitions::Int64 = 100
    end

    function make_factors(sf::Dist{T},
            range::Vector{T},
            parranges::NTuple{N,Vector},
            id,
            parids::Tuple)::Tuple{Vector{<:Scruff.Utils.Factor}, Vector{<:Scruff.Utils.Factor}} where {T,N}

        (lowers, uppers) = bounded_probs(sf, range, ())
        keys = (id,)
        dims = (length(range),)
        return ([Factor(dims, keys, lowers)], [Factor(dims, keys, uppers)])
    end
end

@impl begin
    struct DistSendLambda end

    function send_lambda(sf::Dist{T},
                       lambda::Score{<:T},
                       range::__OptVec{<:T},
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_idx::Integer)::Score where {N,T}

        SoftScore(Float64[], Float64[])
    end
end


text:
"""
logplots.jl contains code to help gain insight into belief propagation.
It creates a Logger called `BPLogger` that plots the results, for
a single variable, of the following methods in the three pass BP algorithm:

* compute_pi
* compute_lambda
* compute_bel
* incoming_pi
* outgoing_pis

To use `BPLogger`, first install `Plots` and use `pyplot()`.  Then,

```
julia> include("src/utils/logplots.jl")
histogram_plot_discrete (generic function with 1 method)

julia> logger = BPLogger(:x1) # :x1 is the name of the variable being plotted
BPLogger(:x1)

julia> with_logger(logger) do
         include("mybpcall.jl") # mybpcall.jl contains a call to threepassbp()
       end
```

"""

export BPLogger

import Logging: handle_message, shouldlog, min_enabled_level

using Plots
using Logging

struct BPLogger <: AbstractLogger
    varname::Vector{Symbol}
end
BPLogger(varname::Symbol...) = BPLogger(collect(varname))

function handle_message(
        logger::BPLogger,
        level,
        message,
        _module,
        group,
        id,
        filepath,
        line; kwargs...)

    if isempty(kwargs) || !in(:type, keys(kwargs))
        with_logger(global_logger()) do
            @logmsg level message kwargs...
        end
        return
    end

    d = Dict(kwargs...)
    type = get(d, :type, :cont)
    if (type != :cont && type != :discrete)
        write_missing("Value of :type must be either :cont or :discrete, not $(type)",
            level, _module, filepath, line; kwargs...)
        return
    end

    if (type == :cont && !in(:numBins, keys(kwargs)))
        write_missing("Continuous missing a required parameter [:numBins]",
            level, _module, filepath, line; kwargs...)
        return
    end

    if !issubset([:range,:prob,:varname,:fname,:name], keys(kwargs))
        write_missing("Missing a required parameter [:range,:prob,:varname,:fname,:name]",
            level, _module, filepath, line; kwargs...)
        return
    end

    if (d[:varname] in logger.varname)
        if (type == :cont)
            histogram_plot(
                d[:range],
                d[:prob],
                d[:numBins],
                d[:varname],
                d[:fname],
                d[:name])
        else
            histogram_plot_discrete(
                d[:range],
                d[:prob],
                d[:varname],
                d[:fname],
                d[:name])
        end
    end
end

function write_missing(msg, level, _module, filepath, line; kwargs...)
    buf = IOBuffer()
    iob = IOContext(buf, stderr)
    levelstr = level == Logging.Warn ? "Warning" : string(level)
    msglines = split(chomp(string(msg)::String), '\n')
    println(iob, "┌ ", levelstr, ": ", msglines[1])
    for i in 2:length(msglines)
        println(iob, "│ ", msglines[i])
    end
    for (key, val) in kwargs
        println(iob, "│   ", key, " = ", val)
    end
    println(iob, "└ @ ", _module, " ", filepath, ":", line)
    write(stderr, take!(buf))
    nothing
end

function shouldlog(logger::BPLogger, level, _module, group, id)
    true
    # group == :threepassbp
end

function min_enabled_level(logger::BPLogger)
    Logging.Debug
end

"""
    histogram_plot(range, prob, numBins::Int64, name::String)

A utility function that, if the Plots module is loaded, returns a histogram
bar graph; otherwise it returns a string with the values that would have
been a histogram

"""
function histogram_plot(range, prob, numBins::Int64, varname::Symbol, fname, name::String)
    @debug("histogram_plot", range=range,
                prob=prob,numBins=numBins,varname=varname,
                fname=fname,name=name)
    if isempty(prob)
        prob = zeros(length(range))
    end
    range_min = minimum(range)
    range_max = maximum(range)
    if(range_min==range_max)
        range_max = range_min+1
    end
    numBins = max(min(length(range), numBins),1)
    stepsize = (range_max - range_min)/ numBins
    samplespace = range_min : stepsize : range_max
    bins = zeros(numBins)
    for (i,lb) in enumerate(samplespace)
        curr_bin_min = samplespace[i]
        curr_bin_max = samplespace[i] + stepsize
        idx = findall(x -> curr_bin_min <= range[x] < curr_bin_max, 1:length(range))
        if(!isempty(idx))
            if(i==length(samplespace)) # last elements
                bins[i-1] += sum(prob[idx])
            else
                bins[i] = sum(prob[idx])
            end
        end
    end
    x_min = samplespace[1]
    x_max = samplespace[end] + stepsize

    @debug "bin_centers=$(collect(samplespace)[1:numBins] .+ (stepsize/2)), bins=$bins"
    gui(bar(samplespace[1:numBins].+(stepsize/2), bins; reuse=false,title="$(varname).$(fname)", label=name, xlims=(x_min-2, x_max+2), hover = samplespace[1:numBins].+(stepsize/2), legend=true))
end

function histogram_plot_discrete(range, prob:: Array{Float64,1}, varname, fname, name::String)
    gui(bar(string.(range),prob; reuse=false,title="$(varname).$(fname)",label=name, hover = string.(range), legend=true))
end


text:
# abstract type NLP end

# struct ControlProblem <: NLP

# end


text:
@generated function __new__(T, args...)
  quote
    Base.@_inline_meta
    $(Expr(:new, :T, [:(args[$i]) for i = 1:length(args)]...))
  end
end

@generated function __splatnew__(T, args)
  quote
    Base.@_inline_meta
    $(Expr(:splatnew, :T, :args))
  end
end

literal_getindex(x, ::Val{i}) where i = getindex(x, i)
literal_indexed_iterate(x, ::Val{i}) where i = Base.indexed_iterate(x, i)
literal_indexed_iterate(x, ::Val{i}, state) where i = Base.indexed_iterate(x, i, state)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module LazIO

using Libdl
using FileIO
using LasIO
using FixedPointNumbers
using LASzip_jll

# Module initialization function
function __init__()
    # temporary _ until LasIO defers this key
    add_format(format"LAZ_", (), ".laz", [:LazIO])
end

include("utils.jl")
include("laszip_h.jl")
include("laszip.jl")
include("fileio.jl")
include("convert.jl")
include("dataset.jl")
include("write.jl")
include("table.jl")

export return_number
export number_of_returns
export scan_direction
export edge_of_flight_line
export classification
export synthetic
export key_point
export withheld

end # module

max_stars_repo_name: evetion/LazIO

max_issues_repo_name: evetion/LazIO

max_forks_repo_name: evetion/LazIO


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/19/2014

module CommonInterfaces

export
    addObserver,
    clearObservers!,

    initialize,
    step

function addObserver()
end

function clearObservers!()
end

function initialize()
end

function update()
end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
const Symbolic = Union{Expr, Symbol}
const Value = Union{Symbolic, Number}

struct Rule
    pat
    rpat
    check
end


text:
using DifferentialEquations, Test

lags = [.2]
f = function (u,h,p,t)
  out = -h(p,t-.2) + u
end
h = (p,t) -> 0.0

prob = DDEProblem(f,1.0,h,(0.0,10.0),constant_lags = lags)

sol = solve(prob)

@test typeof(sol.alg) <: CompositeAlgorithm
@test typeof(sol.alg.algs[1]) <: Tsit5
@test typeof(sol.alg.algs[2]) <: Rosenbrock23

max_stars_repo_name: ChrisRackauckas/DiffEq


text:
using JuMP
using CPLEX
include("parser.jl")
include("relaxed.jl")
include("base.jl")

export DW

function DW(fileName :: String, X_i::Array{Float64,2},  X_fi::Array{Float64,3}, X_ikf::Array{Float64,4}, E::Array{Float64,5})

	#Data acquisition, model definition
	cout_ouverture, Fct_commod, func_cost, func_capacity, nb_nodes, nb_arcs, nb_commodities, latency, node_capacity, commodity, nb_func, exclusion = read_instance(fileName)
	score_init=sum(cout_ouverture*X_i[1,i]+sum(func_cost[f,i]*X_fi[1,f,i] for f in 1:nb_func) for i in 1:nb_nodes)
	max_func_cap=maximum(func_capacity)

	#state sub-problem definition
	m0 = Model(CPLEX.Optimizer)
	set_silent(m0)
	@variable(m0, x_i[1:nb_nodes] >= 0)
	@variable(m0, x_fi[1:nb_func,1:nb_nodes]>=0)
	@constraint(m0,[i in 1:nb_nodes], x_i[i] <= 1) # Relaxed binary constraint
	@constraint(m0,[i in 1:nb_nodes], sum(x_fi[f,i] for f in 1:nb_func) <= node_capacity[i], base_name = "cap" ) # (4) capacité de noeud

	#flow sub problem definition

	m_k=Dict{Int,JuMP.Model}()
	m_k[0]=m0

	for k in 1:nb_commodities

		#model definition
		mk = Model(CPLEX.Optimizer)
		set_silent(mk)
		@variable(mk, x_ikf[1:nb_nodes,1:nb_func]>=0, base_name="x_ikf_"*string(k))
		@variable(mk, e[1:nb_nodes,1:nb_nodes,1:nb_func+1] >= 0,base_name="e_"*string(k))

		for i in 1:nb_nodes
			@constraint(mk,[f in 1:nb_func], x_ikf[i,f] <= 1) # Relaxed binary constraint
			for j in 1:nb_nodes
				@constraint(mk,[f in 1:nb_func+1], e[i,j,f] <= 1) # Relaxed binary constraint
			end
		end

		#Constraints of exclusion and possible arc
		for i in 1:nb_nodes
			if exclusion[k,1]>0
				@constraint(mk,sum(x_ikf[i,w] for w in exclusion[k,:]) <= 1, base_name = "exclusion_"*string(k)*"_"*string(i)) #exclusion
			end
			for j in 1:nb_nodes
				if latency[i,j] == 0
					@constraint(mk,[f in 1:nb_func+1],e[i,j,f]==0, base_name = "no_arc_"*string(i)*"_"*string(j)) #absence d'arcs
				end
			end
		end

		#Contrainte de flots successifs
		size_fk=length( findall( y -> y > 0, Fct_commod[k,:]))
		tab_fk=sortperm( Fct_commod[k,:])[1:size_fk]
		@constraint(mk,sum(e[commodity[k,1],j,tab_fk[1]] - e[j,commodity[k,1],tab_fk[1]] for j in 1:nb_nodes)-x_ikf[commodity[k,1],tab_fk[1]]==1, base_name = "init_flot_"*string(k)) #initialisation du flot à la source
		@constraint(mk,[i in 1:nb_nodes ; i!=commodity[k,1]],sum(e[i,j,tab_fk[1]] - e[j,i,tab_fk[1]]  for j in 1:nb_nodes)+x_ikf[i,tab_fk[1]]==0, base_name = "1er_flot_"*string(k)) #flot pour traiter la première fonction de puis la source
		for t in 2:size_fk
			@constraint(mk,[i in 1:nb_nodes],sum(e[i,j,tab_fk[t]] - e[j,i,tab_fk[t]] for j in 1:nb_nodes)-x_ikf[i,tab_fk[t-1]]+x_ikf[i,tab_fk[t]]==0, base_name = "mid_flow_"*string(t)*"_"*string(k)) #flot d'un traitement au suivant
		end
		@constraint(mk,[i in 1:nb_nodes ; i!=commodity[k,2]],sum(e[i,j,nb_func+1] - e[j,i,nb_func+1] for j in 1:nb_nodes)-x_ikf[i,tab_fk[end]]==0, base_name = "last_flow_"*string(k)) #flot du dernier traitement vers le puit
		@constraint(mk,sum(e[j,commodity[k,2],end] for j in 1:nb_nodes)+x_ikf[commodity[k,2],tab_fk[end]]==1, base_name = "end_flow_"*string(k)) #fin de flot sur le puit

		for f in 1:nb_func
			@constraint(mk,sum(x_ikf[i,f] for i in 1:nb_nodes)<=1, base_name = "concatenation_"*string(k)*"_"*string(f)) #concatenation des problemes de flot
		end

		#Contraintes de latence
		@constraint(mk, sum( sum( sum(e[i,j,f]*latency[i,j] for f in 1:nb_func+1) for j in 1:nb_nodes) for i in 1:nb_nodes) <= commodity[k,4], base_name = "latence_"*string(k))

		m_k[k]=mk
	end

	#Modele maitre dual
	m = Model(CPLEX.Optimizer)
	set_silent(m)
	@variable(m,mu_2[1:nb_nodes,1:nb_commodities,1:nb_func]>=0)
	@variable(m,mu_3[1:nb_nodes,1:nb_func]>=0)
	@variable(m,eta)
	@constraint(m, [ i in 1:nb_nodes*nb_func], mu_3[i] <= score_init/max_func_cap)
	@constraint(m, [ i in 1:nb_nodes*nb_func*nb_commodities], mu_2[i] <= score_init/max_func_cap)
	@objective(m, Max, eta )

	#Column generation
	not_opti=true
	n=0
	while not_opti
		n+=1

		@constraint(m, sum(X_i[end,i]*(cout_ouverture-sum(sum(mu_2[i,k,f] for f in 1:nb_func) for k in 1:nb_commodities))+sum(X_fi[end,f,i]*(func_cost[f,i]-mu_3[i,f]*func_capacity[f]) for f in 1:nb_func)  for i in 1:nb_nodes) - sum( sum(sum( (-mu_3[i,f]*commodity[k,3] - mu_2[i,k,f])*X_ikf[end,i,k,f] for i in 1:nb_nodes) for f in 1:nb_func) for k in 1:nb_commodities) - eta >=0, base_name="con_"*string(n))
		optimize!(m)

		#Get dual values and solutions
		eta_sol=JuMP.value(eta)
		mu_2_sol=JuMP.value.(mu_2)
		mu_3_sol=JuMP.value.(mu_3)

		#setup for next column to add on last index
		min_cr=-eta_sol #Minimum ecart complementaire
		X_i=vcat(X_i,X_i[1:1,:])
		X_fi=vcat(X_fi,X_fi[1:1,:,:])
		X_ikf=vcat(X_ikf,X_ikf[1:1,:,:,:])
		E=vcat(E,E[1:1,:,:,:,:])

		#state sub-problem resolution
		@objective(m_k[0], Min, sum(x_i[i]*(cout_ouverture-sum(sum(mu_2_sol[i,k,f] for f in 1:nb_func) for k in 1:nb_commodities))+sum(x_fi[f,i]*(func_cost[f,i]-mu_3_sol[i,f]*func_capacity[f]) for f in 1:nb_func)  for i in 1:nb_nodes) )
		optimize!(m_k[0])
		min_cr=min_cr+objective_value(m_k[0])
		X_i[end,:]=JuMP.value.(x_i)
		X_fi[end,:,:]=JuMP.value.(x_fi)

		#flow sub-problem resolution
		for k in 1:nb_commodities
			x_ikf=object_dictionary(m_k[k])[:x_ikf]
			e=object_dictionary(m_k[k])[:e]
			@objective(m_k[k], Min, - sum(sum( (-mu_3_sol[i,f]*commodity[k,3] - mu_2_sol[i,k,f])*x_ikf[i,f] for i in 1:nb_nodes) for f in 1:nb_func) )
			optimize!(m_k[k])
			min_cr = min_cr + objective_value(m_k[k])
			X_ikf[end,:,k,:]=JuMP.value.(x_ikf)
			E[end,:,:,k,:]=JuMP.value.(e)
		end

		#Check optimality
		println(min_cr)
		if min_cr>=-1e-7
			println(solution_summary(m))
			not_opti=false
			lbd_sol=Vector{Float64}(zeros(n))
			for i in 1:n
				lbd_sol[i]=dual(constraint_by_name(m,"con_"*string(i)))
			end
			x_i=sum(lbd_sol[w].*X_i[w,:] for w in 1:n)
			x_fi=sum(lbd_sol[w].*X_fi[w,:,:] for w in 1:n)
			x_ikf=sum(lbd_sol[w].*X_ikf[w,:,:,:] for w in 1:n)
			e=sum(lbd_sol[w].*E[w,:,:,:,:] for w in 1:n)

			return x_i, x_fi, x_ikf, e
		end
	end

	return 0

end

ov, x_i, x_fi, x_ikf, e = Static("di-yuan/di-yuan_1/",true)
X_i=Array{Float64,2}(zeros(1,size(x_i,1)))
X_i[1,:]=x_i
X_fi=Array{Float64,3}(zeros(1,size(x_fi,1),size(x_fi,2)))
X_fi[1,:,:]=x_fi
X_ikf=Array{Float64,4}(zeros(1,size(x_ikf,1),size(x_ikf,2),size(x_ikf,3)))
X_ikf[1,:,:,:]=x_ikf
E=Array{Float64,5}(zeros(1,size(e,1),size(e,2),size(e,3),size(e,4)))
E[1,:,:,:,:]=e

x_i, x_fi, x_ikf, e=DW_mp("di-yuan/di-yuan_1/",X_i,X_fi,X_ikf,E)
println("fin")
#x_i, x_fi, x_ikf, e=Relaxed("di-yuan/di-yuan_1/")

max_stars_repo_name: VSpitzer/RORT

max_issues_repo_name: VSpitzer/RORT

max_forks_repo_name: VSpitzer/RORT


text:
using Zygote: @adjoint
### helper functions for working with array types
addparent(A::Type{<:AbstractArray}, P′::DataType = A) = A
addparent(A::Type{CuArray{T,N,Nothing}}, P′::DataType = A) where {T,N} = CuArray{T,N,P′}
addparent(A::Type{CuArray{T,N,P}}, P′::DataType) where {T,N,P} = CuArray{T,N,P′}

vectorof(::Type{<:Array{T}}) where T = Vector{T}
matrixof(::Type{<:Array{T}}) where T = Matrix{T}
tensor₃of(::Type{<:Array{T}}) where T = Array{T,3}

vectorof(::Array{T}) where T = Vector{T}
matrixof(::Array{T}) where T = Matrix{T}
tensor₃of(::Array{T}) where T = Array{T,3}

vectorof(::Type{CuArray{T,N,P}}) where {T,N,P} = CuVector{T,P}
matrixof(::Type{CuArray{T,N,P}}) where {T,N,P} = CuMatrix{T,P}
tensor₃of(::Type{CuArray{T,N,P}}) where {T,N,P} = CuArray{T,3,P}

vectorof(::CuArray{T,N,P}) where {T,N,P} = CuVector{T,P}
matrixof(::CuArray{T,N,P}) where {T,N,P} = CuMatrix{T,P}
tensor₃of(::CuArray{T,N,P}) where {T,N,P} = CuArray{T,3,P}

### Work around for github.com/FluxML/Flux.jl/issues/1084
function Flux.dropout(x, p)
   q = 1 - p
   y = rand!(similar(x))
   y .= Flux._dropout_kernel.(y, p, q)
   x .* y
end

@adjoint function Flux.dropout(x, p)
   q = 1 - p
   y = rand!(similar(x))
   y .= Flux._dropout_kernel.(y, p, q)
   return x .* y, Δ -> (Δ .* y, nothing)
end

### Custom adjoint definitions
@adjoint function Base.reduce(::typeof(hcat), As::AbstractVector{<:AbstractVecOrMat})
   cumsizes = cumsum(size.(As, 2))
   return reduce(hcat, As), Δ -> (nothing, map((sz, A) -> Zygote.pull_block_horz(sz, Δ, A), cumsizes, As))
end

@adjoint function Base.reduce(::typeof(vcat), As::AbstractVector{<:AbstractVecOrMat})
   cumsizes = cumsum(size.(As, 1))
   return reduce(vcat, As), Δ -> (nothing, map((sz, A) -> Zygote.pull_block_vert(sz, Δ, A), cumsizes, As))
end

veclength(::Tuple) = 1
veclength(x::AbstractVector{<:Tuple}) = length(x)
pull_block_vert(idx, ȳ, x::Tuple) = ȳ[idx]
pull_block_vert(lastidx, ȳ, x::AbstractVector{<:Tuple}) = ȳ[(lastidx-length(x)+1):lastidx]

@adjoint function Base.vcat(xs::Union{AbstractVector{T},T}...) where T <: Tuple
   lastidx = Ref(0)
   lastidxs = ntuple(length(xs)) do i
      @inbounds lastidx[] += veclength(xs[i])
   end
   return vcat(xs...), function (ȳ)
      (map((lastidx, x) -> pull_block_vert(lastidx, ȳ, x), lastidxs, xs)...,)
   end
end

pull_block_horz(lastidx, ȳ, x::AbstractArray{<:Number,3}) = ȳ[:, (lastidx-size(x,2)+1):lastidx, :]

@adjoint function Base.hcat(xs::T₃...) where T₃ <: AbstractArray{<:Number,3}
   lastidx = Ref(0)
   lastidxs = ntuple(length(xs)) do i
      @inbounds lastidx[] += size(xs[i], 2)
   end
   return hcat(xs...), function (ȳ)
      (map((lastidx, x) -> pull_block_horz(lastidx, ȳ, x), lastidxs, xs)...,)
   end
end

Zygote.refresh()


text:
export create_samples, dist_create_samples, save_results_csv, save_stats,
	   save_polytope, results_to_array
export sample_polytope_cprnd, sample_uniform, sample_uniform_w_pf

import Distributed
import Distributed: nprocs, addprocs, rmprocs
export Distributed, nprocs, addprocs, rmprocs

import JuMP: optimizer_with_attributes
export JuMP, optimizer_with_attributes

import Ipopt
export Ipopt

import PowerModels
const PM = PowerModels
export PM, PowerModels

import DelimitedFiles
export DelimitedFiles  # For testing io functions

# The following code exports symbols used throughout OPFLearn. If you don't want all of these
# symbols in your environment, then use `import OPFLearn` instead of `using OPFLearn`.
include_all = false
if include_all
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    println("$(sym)")
    @eval export $sym
end
end


text:
function LongitudeExtremumDetector(arg0::OneAxisEllipsoid)
    return LongitudeExtremumDetector((OneAxisEllipsoid,), arg0)
end

function LongitudeExtremumDetector(arg0::jdouble, arg1::jdouble, arg2::OneAxisEllipsoid)
    return LongitudeExtremumDetector((jdouble, jdouble, OneAxisEllipsoid), arg0, arg1, arg2)
end

function g(obj::LongitudeExtremumDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_body(obj::LongitudeExtremumDetector)
    return jcall(obj, "getBody", BodyShape, ())
end


text:
######### StanSample Bernoulli example  ###########

using StanSample

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(@__DIR__, "tmp")

sm = SampleModel("bernoulli", bernoulli_model;
  method = StanSample.Sample(
    save_warmup=false,                           # Default
    thin=1,
    adapt = StanSample.Adapt(delta = 0.85)),
  tmpdir = tmpdir
);

rc = stan_sample(sm; data=bernoulli_data);

if success(rc)
  (samples, cnames) = read_samples(sm, output_format=:array)
end


text:
#aproximates liquid volume at a known pressure and t,
#by using isothermal compressibility
function volume_compress(model,p,T,z=SA[1.0];v0=vcompress_v0(model,p,T,z),max_iters=100)
    logv0 = log(v0)
    function f_fixpoint(_v)
        _v = exp(_v)
        _p,dpdv = p∂p∂v(model,_v,T,z)
        β = -1/_v*dpdv^-1
        _Δ =  -(p-_p)*β
        sign_Δ = sign(_Δ)
        Δ = abs(_Δ)
        vv = _v*exp(sign_Δ*Δ^(1-_Δ))#^((1+Δ)^4)
        return log(vv)
    end
        res = Solvers.fixpoint(f_fixpoint,logv0,Solvers.SimpleFixPoint(),rtol = 1e-12,max_iters=max_iters)
        return exp(res)
end

function vcompress_v0(model,p,T,z=SA[1.0])
    lb_v   = lb_volume(model,z)
    v0 = 1.1*lb_v
    return v0
end

function vcompress_v0(model::SAFTVRMieModel,p,T,z=SA[1.0])
    lb_v   = lb_volume(model,z)
    v0 = 1.5*lb_v
    return v0
end

function volume_virial(model,p,T, z=SA[1.] )
    B = second_virial_coefficient(model,T,z)
    a = p/(R̄*T)
    b = -1
    c = -B
    Δ = b*b-4*a*c
    n = sum(z)
    if Δ <= 0
        #virial approximation could not be calculated
        #return value at spinodal
        return -2*B
    end
    #only the left root has physical meaning
    return (-b + sqrt(b*b-4*a*c))/(2*a)
end

#(z = pv/rt)
#(RT/p = v/z)

function volume(model::EoSModel,p,T,z=SA[1.0];phase=:unknown,threaded=true)

    fp(_v) = log(pressure(model,_v,T,z)/p)

#Threaded version
    phase = Symbol(phase)
    if phase != :unknown
        v0 = x0_volume(model,p,T,z,phase=phase)
        #return Solvers.ad_newton(fp,vg0)
        return volume_compress(model,p,T,z,v0=v0)
    end

    if threaded
        _vg = Threads.@spawn begin
            vg0 = x0_volume(model,$p,$T,$z,phase=:v)
            volume_compress(model,$p,$T,$z;v0=vg0)
            #Solvers.ad_newton(fp,vg0)
        end
        vl0 = x0_volume(model,p,T,z,phase=:l)
        _vl = Threads.@spawn volume_compress(model,$p,$T,$z;v0=$vl0)
        #fp(_v) = pressure(model,_v,T,z) - p
        vg = fetch(_vg)
        vl = fetch(_vl)
    else
        vg0 = x0_volume(model,p,T,z,phase=:v)
        vl0 = x0_volume(model,p,T,z,phase=:l)

        vg =  volume_compress(model,p,T,z,v0=vg0)
        #vg = Solvers.ad_newton(fp,vg0)
        vl =  volume_compress(model,p,T,z,v0=vl0)
    end

# Serial version
#=
    vg0 = volume_virial(model,p,T,z)
    vg =Solvers.ad_newton(fp,vg0,rtol=1e-08)

    vl =  volume_compress(model,p,T,z)
=#
    function gibbs(v)
        _df,f =  ∂f(model,v,T,z)
        dv,dt = _df
        if abs((p+dv)/p) > 0.03
            return Inf
        else
            return f  +p*v
        end
    end
    #this catches the supercritical phase as well

    if isnan(vl)
        return vg
    end
    if isnan(vg)
        return vl
    end
    if vl ≈ vg
        return vl
    end
        gg = gibbs(vg)
        gl = gibbs(vl)
        #@show vg,vl
        #@show gg,gl
        return ifelse(gg<gl,vg,vl)

end

export volume


text:
import LinearAlgebra

export ucase,
    ucrss,
    uddf,
    union,
    unitim,
    unload,
    unorm,
    unormg,
    utc2et

function _ucase(in)
    n = length(in) + 1
    out = Array{SpiceChar}(undef, n)
    ccall((:ucase_c, libcspice), Cvoid, (Cstring, SpiceInt, Ref{SpiceChar}),
          in, n, out)
    chararray_to_string(out)
end

"""
    ucase(in)

!!! warning "Deprecated"
    Use `uppercase(in)` instead.
"""
ucase

@deprecate ucase uppercase

function _ucrss(v1, v2)
    vout = Array{SpiceDouble}(undef, 3)
    ccall((:ucrss_c ,libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          v1, v2, vout)
    vout
end

"""
    ucrss(v1, v2)

!!! warning "Deprecated"
    Use `LinearAlgebra.normalize(LinearAlgebra.cross(v1, v2))` instead.
"""
ucrss

@deprecate ucrss(v1, v2) LinearAlgebra.normalize(LinearAlgebra.cross(v1, v2))

"""
    uddf(udfunc, x, dx)

Routine to calculate the first derivative of a caller-specified function using
a three-point estimation.

### Arguments ###

- `udfunc`: A callable that computes the scalar value of interest,
    e.g. `f(x::Float64) -> Float64`
- `x`: Independent variable of `udfunc`
- `dx`: Interval from `x` for derivative calculation

### Output ###

Returns the approximate derivative of `udfunc` at `x`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/uddf_c.html)
"""
function uddf(udfunc, x, dx)
    function _udfunc(et::SpiceDouble, value::Ptr{SpiceDouble})
        value_ = unsafe_wrap(Array{SpiceDouble}, value, 1)
        value_[1] = udfunc(et)
        nothing
    end
    func = @cfunction($_udfunc, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    deriv = Ref{SpiceDouble}()
    ccall((:uddf_c, libcspice), Cvoid,
          (Ref{Cvoid}, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          func, x, dx, deriv)
    deriv[]
end

"""
    union(a::T, b::T) where T <: SpiceCell

Compute the union of two sets of any data type to form a third set.

### Arguments ###

- `a`: First input set
- `b`: Second input set

### Output ###

Returns a cell containing the union of `a` and `b`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/union_c.html)
"""
function Base.union(a::SpiceCell{T}, b::SpiceCell{T}) where T
    size = max(a.cell.size, b.cell.size)
    if T <: SpiceChar
        length = max(a.cell.length, b.cell.length)
        out = SpiceCell{T}(size, length)
    else
        out = SpiceCell{T}(size)
    end
    a_data = a.data
    b_data = b.data
    out_data = out.data
    GC.@preserve a_data b_data out_data begin
        ccall((:union_c, libcspice), Cvoid, (Ref{Cell{T}}, Ref{Cell{T}}, Ref{Cell{T}}),
              a.cell, b.cell, out.cell)
    end
    handleerror()
    out
end

"""
    unitim(epoch, insys, outsys)

Transform time from one uniform scale to another.

### Arguments ###

- `epoch`: An epoch to be converted
- `insys`: The time scale associated with the input epoch
- `outsys`: The time scale associated with the function value

The uniform time scales are:

- `:TAI`
- `:TDT`
- `:TDB`
- `:ET`
- `:JED`
- `:JDTDB`
- `:JDTDT`

### Output ###

Returns the time in the system specified by `outsys` that is equivalent to the
`epoch` in the `insys` time scale.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/unitim_c.html)
"""
function unitim(epoch, insys, outsys)
    out = ccall((:unitim_c, libcspice), SpiceDouble,
                (SpiceDouble, Cstring, Cstring),
                epoch, string(insys), string(outsys))
    handleerror()
    out
end

"""
    unload(file)

Unload a SPICE kernel.

### Arguments ###

- `file`: The file name of a kernel to unload

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/unload_c.html)
"""
function unload(file)
    ccall((:unload_c, libcspice), Cvoid, (Cstring,), file)
    handleerror()
end

function _unorm(v1)
    vout = Array{SpiceDouble}(undef, 3)
    vmag = Ref{SpiceDouble}()
    ccall((:unorm_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          v1, vout, vmag)
    vout, vmag[]
end

"""
    unorm(v1)

!!! warning "Deprecated"
    Use `(LinearAlgebra.normalize(v1), LinearAlgebra.norm(v1))` instead.
"""
unorm

@deprecate unorm(v) (normalize(v), norm(v))

function _unormg(v1)
    ndim = length(v1)
    vout = Array{SpiceDouble}(undef, ndim)
    vmag = Ref{SpiceDouble}()
    ccall((:unormg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}),
          v1, ndim, vout, vmag)
    vout, vmag[]
end

"""
    unormg(v1)

!!! warning "Deprecated"
    Use `(LinearAlgebra.normalize(v1), LinearAlgebra.norm(v1))` instead.
"""
unormg

@deprecate unormg(v) (normalize(v), norm(v))

"""
    utc2et(utcstr)

Convert an input time from Calendar or Julian Date format, UTC, to ephemeris
seconds past J2000.

### Arguments ###

- `utcstr`: Input time string, UTC

### Output ###

Returns the equivalent of utcstr, expressed in ephemeris seconds past J2000.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/utc2et_c.html)
"""
function utc2et(utcstr)
    et = Ref{SpiceDouble}()
    ccall((:utc2et_c, libcspice), Cvoid, (Cstring, Ref{SpiceDouble}),
          string(utcstr), et)
    handleerror()
    et[]
end


text:
using Test
using MPIClusterManagers
using Distributed

import MPI

# Start workers via `mpiexec` that communicate among themselves via MPI;
# communicate with the workers via TCP
nprocs = parse(Int, ARGS[1])

mgr = MPIManager(np=nprocs)
addprocs(mgr)

refs = []
for w in workers()
    push!(refs, @spawnat w MPI.Comm_rank(MPI.COMM_WORLD))
end
ids = falses(nworkers())
for r in refs
    id = fetch(r)
    @test !ids[id+1]
    ids[id+1] = true
end
for id in ids
    @test id
end

s = @distributed (+) for i in 1:10
    i^2
end
@test s == 385

@mpi_do mgr n = 2
@mpi_do mgr n + 2


text:
# Do not modify this file
module Assignment02

using Base: Bool
export question1,
       question2,
       question3,
       question4,
       question5,
       generate_sequence,
       bookprice

include("assignment.jl")

end


text:
using ITensors
using ITensorImpurity
using UnicodePlots
using Random

Random.seed!(1234)

# Number of sites
# N = 16, 32, 64, 128
N = 8
Nᴸ = N ÷ 2
Nᴿ = N ÷ 2 - 1
@assert iseven(Nᴸ)
@assert isodd(Nᴿ)
@assert Nᴸ + Nᴿ + 1 == N

@show N
@show Nᴸ
@show Nᴿ

# Bath hopping
t = 1.0

# Time steps
# Typical time steps: δτ = 0.01 - 0.1
# Can use larger time steps, up to δτ ≈ 0.4,
# deeper in the Kondo regime (large U/t′²)
# Reference: https://arxiv.org/pdf/0807.0581.pdf
δτ = 0.05t
τf = 6.0t

# Hybridization paramater
U = t
t′ = 0.4
# Source: https://arxiv.org/pdf/0807.0581.pdf
# U / Γ = 3.125
Γ = 2 * t′ ^ 2 / t

# Bias potential
V = 0.005U

# Impurity potential
# Particle-hole symmetric point is: Vᵍ = -U / 2
Vᵍ = -U / 2

@show t, U, Γ, U / Γ, t′, V, Vᵍ

s = siteinds("Electron", N)
# Without bias potential
ℋ₀ = anderson(N; Nᴸ, t, V=0.0, t′, U, Vᵍ)
# With bias potential
ℋ = anderson(N; Nᴸ, t, V, t′, U, Vᵍ)

H₀ = MPO(ℋ₀, s)
H = MPO(ℋ, s)

ψ0 = randomMPS(s; linkdims=4)

sweeps = Sweeps(10)
setmaxdim!(sweeps, 50, 100)
setcutoff!(sweeps, 1e-8)
e, ψ0 = dmrg(H₀, ψ0, sweeps)

ntot = expect(ψ0, "ntot")
println("⟨ψ|ntot|ψ⟩ = ")
display(ntot)
@show N
@show sum(ntot)

expH = exp(Val{:trotter_order_2}(), -im * δτ, ℋ, s)
U = MPO(s, "Id")
U = apply(expH, U; cutoff=1e-6)

if N ≤ 4
  # Compare trotter decomposition with ED
  ∏U = contract_slice(linkinds, U)
  ∏H = contract_slice(linkinds, H)
  @show norm(∏U - exp(-im * δt * ∏H))
end

τi = 0.0t
τ⃗ = τi:δτ:τf
nlayers = length(τ⃗)
f(ψτ) = ITensorImpurity.J(ψτ, t′; Nᴸ) / V
ψτ, J⃗ = apply_layers([expH for _ in 1:(nlayers - 1)], ψ0; f, cutoff=1e-8)

p = lineplot(τ⃗, real(J⃗))
display(p)

nothing

max_stars_repo_name: mtfishman/ITensorImpurity

max_issues_repo_name: mtfishman/ITensorImpurity

max_forks_repo_name: mtfishman/ITensorImpurity


text:
using SafeTestsets

@safetestset "util" begin include("util.jl") end
@safetestset "tnmean" begin include("tnmean.jl") end
@safetestset "tnmom2" begin include("tnmom2.jl") end
@safetestset "tnmom1c" begin include("tnmom1c.jl") end
@safetestset "tnmom2c" begin include("tnmom2c.jl") end
@safetestset "tnvar" begin include("tnvar.jl") end
@safetestset "truncnorm" begin include("truncnorm.jl") end
@safetestset "rejection" begin include("rejection.jl") end


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

mutable struct GatewayRouteListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{GatewayRoute} } # spec name: value

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

const _property_map_GatewayRouteListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_GatewayRouteListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{GatewayRoute}")
Base.propertynames(::Type{ GatewayRouteListResult }) = collect(keys(_property_map_GatewayRouteListResult))
Swagger.property_type(::Type{ GatewayRouteListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_GatewayRouteListResult[name]))}
Swagger.field_name(::Type{ GatewayRouteListResult }, property_name::Symbol) =  _property_map_GatewayRouteListResult[property_name]

function check_required(o::GatewayRouteListResult)
    true
end

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


text: return "boom!"; // Expect: "Cannot return from top-level code."


text:
import Libdl

if VERSION < v"1.3" || haskey(ENV, "SCIPOPTDIR")
    # User-provided SCIP library
    const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
    if !isfile(depsjl_path)
        error("SCIP was not built properly, please run Pkg.build(\"SCIP\")")
    end
    include(depsjl_path)
else
    # Artifact from BinaryBuilder package
    using SCIP_jll: libscip
end

function __init__()
    major = SCIPmajorVersion()
    minor = SCIPminorVersion()
    patch = SCIPtechVersion()
    current = VersionNumber("$major.$minor.$patch")
    required = VersionNumber("6")
    upperbound = VersionNumber("8")
    if current < required || current >= upperbound
        error("SCIP is installed at version $current, " *
              "supported are $required up to (excluding) $upperbound.")
    end
end


text:
# test simple blob detector
img_gray = OpenCV.imread(joinpath(test_dir, "shared", "pic1.png"), OpenCV.IMREAD_GRAYSCALE)

detector = OpenCV.SimpleBlobDetector_create()

# Compare centers of keypoints and se how many of them match,
kps = OpenCV.detect(detector, img_gray)

kps_expect = [OpenCV.Point{Float32}(174.9114f0, 227.75146f0),OpenCV.Point{Float32}(106.925545f0, 179.5765f0)]
for kp in kps
    closest_match = 100000
    for kpe in kps_expect
        dx = kpe.x - kp.pt.x
        dy = kpe.y - kp.pt.y
        if sqrt(dx*dx+dy*dy) < closest_match
            closest_match = sqrt(dx*dx+dy*dy)
        end
    end

    @test closest_match < 10
end

println("feature2d test passed")


text:
module DAinfTests

using Test
using LightPropagation

struct DAsemiinf_test{T <: AbstractFloat} <: DiffusionParameters
    ρ::T                                 # distance away from isotropic point source (cm)
    μa::T                                # absorption coefficient (cm⁻¹)
    μsp::T                               # reduced scattering coefficient (cm⁻¹)
    n_med::T                             # medium's index of refraction
    n_ext::T                             # external medium's index of refraction
    ω::T                                 # modulation frequency (1/ns)
    z::T                                 # z depth of detector within medium (cm)
    t::Union{T, AbstractVector{T}}       # time vector (ns)

    # do not need to provide these arguments (calculated from previous inputs)
    D::T                                 # Diffusion coefficient
    ν::T                                 # Speed of light in medium (cm/ns)
    z0::T                                # isotroptic source depth
    zb::T                                # Extrapolation length
end

# this manually changes the definition of the isotropic source depth to bury deep within a semi-infinite medium
# this approximates an infinite medium when we are close to source to compare the infinite and semi-infinite solutions

SI_test = DAsemiinf_test(0.0, 0.1, 10.0, 1.0, 1.0, 0.0, 20.0, 0.1:0.1:2.5, 1.0 / (3.0 * 10.0), 29.9792458, 19.0, 2.0 / (3.0 * 10.0))
@test fluence_DA_semiinf_CW(SI_test) ≈ fluence_DA_inf_CW(1.0, 0.1, 10.0)

t = 0.1:0.1:2.5
@test fluence_DA_semiinf_TD(SI_test) ≈ fluence_DA_inf_TD(t, 1.0, 0.1, 10.0)

# Hard code DA inf solution to compare other geometries

# CW
@test fluence_DA_inf_CW(0.1, 0.1, 10.0) ≈ 20.076563644369305
@test fluence_DA_inf_CW(0.1, 0.2, 15.0) ≈ 26.52859839465854
@test fluence_DA_inf_CW(60.0, 0.2, 15.0) ≈ 4.007233568620553e-80

# TD
d = [5.87809115624927e-14, 0.0019712791801801076, 2.10256922566133e-6, 3.038559267148949e-9]
@test fluence_DA_inf_TD(0.01:1.0:4.0, 1.0, 0.2, 15.0) ≈ d

end # module

max_issues_repo_name: heltonmc/LightPropagation


text:
"""
This example learns to classify hand-written digits from the
[MNIST](http://yann.lecun.com/exdb/mnist) dataset.  There are 60000
training and 10000 test examples. Each input x consists of 784 pixels
representing a 28x28 image.  The pixel values are normalized to
[0,1]. Each output y is converted to a ten-dimensional one-hot vector
(a vector that has a single non-zero component) indicating the correct
class (0-9) for a given image.  10 is used to represent 0.
You can run the demo using `julia mnist.jl` on the command line or
`julia> include(KnetLayers.dir("examples/mnist.jl"))` at the Julia prompt.
"""

using Knet, KnetLayers
import Knet: Data
#Data
include(Knet.dir("data","mnist.jl"))
dtrn,dtst = mnistdata(xsize=(784,:)); # dtrn and dtst = [ (x1,y1), (x2,y2), ... ] where xi,yi are

#Model
(m::MLP)(x,y) = nll(m(x),y)
(m::MLP)(d::Data) = mean(m(x,y) for (x,y) in d)
HIDDEN_SIZES = [100,50]
model = MLP(784,HIDDEN_SIZES...,10)

#Train
EPOCH=10
progress!(sgd(model,repeat(dtrn,EPOCH)))

#Test
@show 100accuracy(model, dtst)


text:
using DiscreteMarkovChains
using Test
using SafeTestsets

@time begin
    @time @safetestset "Utils" begin include("Utils.jl") end
    @time @safetestset "DMC" begin include("MCCore.jl") end
    @time @safetestset "CMC" begin include("CMC.jl") end
    @time @safetestset "Advanced" begin include("Advanced.jl") end
end


text:
# Create a Range object

struct Range{T}
    start::T
    step::T
    stop::T
end

function Base.getindex(r::Range, i::Int)
    start = r.start
    target = start + i*r.step
    if target <= r.stop
        return target
    else
        throw(DomainError("exceeding range limit $(target) > $(r.stop)"))
    end
end

r = Range(0.0,1.0,10.0)
for i in 1:10
    ri = r[i]
    println("r[$(i)] is $(ri)")
end

# Extension of the Exercise
Base.length(r::Range) = Int((r.stop - r.start) / r.step)

function Base.iterate(r::Range, state = 1)
target = r.start + state * r.step
if target <= r.stop
    return target, state+1
else
    return nothing
end
end

for i in r
    println("Next in range is $(i)")
end

# Whackadoodle

mutable struct MyIterator
    n::Int
    prev::Float64
end

aardvark = MyIterator(10,0.0)

println("The aardvark iterator has $(aardvark.n) elements")

function Base.iterate(i::MyIterator, s = 1)
s > i.n && return nothing
if s==1
    i.prev = rand()
    return i.prev, 2
elseif isodd(s)
    i.prev = i.prev/rand()
    return i.prev, s+1
else
    i.prev = i.prev*rand()
    return i.prev,s+1
end
end

for i in aardvark
    println("Next in aardvark is $(i)")
end

w= MyIterator(500,0.0)
using Statistics: mean, std
std(w)

max_stars_repo_name: AdamVStephen/learning-julia

max_issues_repo_name: AdamVStephen/learning-julia

max_forks_repo_name: AdamVStephen/learning-julia


text:
"""
    heatmap(data)

Creates an `EChart` heatmap.

## Methods
```julia
heatmap(h::StatsBase.Histogram)
```

## Arguments
* `legend::Bool = false` : display legend?
* `show::Bool = true` : show slider denoting color range
* `calculable::Bool = true` : allow user to slide endpoints to change data displayed
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `area`, such as displaying a legend when two or more series present.

## Examples
```julia
using RDatasets
h2 = fit(Histogram, (df[:Price], df[:Carat]), closed = :left)
hs2 = heatmap(h2)
```
"""
function heatmap(h::StatsBase.Histogram;
                    legend::Bool = false,
                    show::Bool = true,
                    calculable::Bool = true,
                    kwargs...)

    #Get labels
    xedges = collect(h.edges[1])
    yedges = collect(h.edges[2])

    h.closed == :left ? left = "[" : left = "("
    left == "[" ? right = ")" : right = "]"

    xlabels = []
    for i in 1:length(xedges)-1
        push!(xlabels, string(left, xedges[i], " - ", xedges[i + 1], right))
    end

    ylabels = []
    for i in 1:length(yedges)-1
        push!(ylabels, string(left, yedges[i], " - ", yedges[i + 1], right))
    end

    #Format data, subtract 1 bc JavaScript!
    xlen, ylen = size(h.weights)
    data = []
    for i in 1:xlen
        for j in 1:ylen
            push!(data, [i - 1, j - 1, h.weights[i,j]])
        end
    end

    #Build heatmap
    ec = ECharts.newplot(kwargs, ec_charttype = "heatmap")
    ec.yAxis = [Axis(_type = "category", data = ylabels)]
    ec.xAxis = [Axis(_type = "category", data = xlabels)]
    ec.series = [HeatmapSeries(name = "heatmap", data = data)]

    minval, maxval = extrema(h.weights)

    ec.visualMap = ECharts.VisualMap(show = show,
                                     min = minval,
                                     max = maxval,
                                     calculable = calculable,
                                     orient = "vertical",
                                     right = "2%",
                                     top = "middle"
                                     )

    #Show legend if requested, though doesn't make sense in this context
    legend ? legend!(ec) : nothing

    return ec

end


text:
using Base.Test

function solve(d)
    n = length(d)
    i = 1
    configs = Dict()

    while true
        v, idx = findmax(d)
        j = idx
        d[j] = 0
        for _ in 1:v
            j = (j) % n + 1
            d[j] += 1
        end
        h = hash(d)
        if h in keys(configs)
            return i, i - configs[h]
        else
            configs[h] = i
        end
        i += 1
    end
end

@test solve([0, 2, 7, 0]) == (5, 4)

d = readdlm("input.txt", Int)
println(solve(d))


text:
#Adapted from https://github.com/github/developer.github.com/blob/master/Rakefile#L21-L48

#Run this file from its current directory to build the documentation and push the result into gh-pages

last_commit=readchomp(`git --no-pager log -1 --pretty=format:"%h:%s"`)

ENV["GIT_DIR"]=abspath(chomp(read(`git rev-parse --git-dir`, String)))

old_sha = chomp(read(`git rev-parse refs/remotes/origin/gh-pages`, String))

#run(`julia make.jl`)

cd("build") do

    gif="/tmp/dev.gh.i"
    ENV["GIT_INDEX_FILE"]=gif
    ENV["GIT_WORK_TREE"]=pwd()
    run(`git add -A`)
    tsha=chomp(read(`git write-tree`, String))
    mesg="Deploy docs for master@$last_commit"

    if length(old_sha) == 40
        csha = chomp(read(`git commit-tree $tsha -p $old_sha -m $(mesg)`, String))
    else
        csha = chomp(read(`git commit-tree $tsha -m $(mesg)`, String))
    end

     print("Created commit $csha")

     run(`git --no-pager show $csha --stat`)

     run(`git update-ref refs/heads/gh-pages $csha `)

     run(`git push origin gh-pages `)

end


text:
# Apply default faceting look to a grid of `AxisEntries`

function facet_wrap!(fig, aes::AbstractMatrix{AxisEntries}; facet)

    scale = get(aes[1].categoricalscales, :layout, nothing)
    isnothing(scale) && return

    # Link axes and hide decorations if appropriate
    attr = clean_facet_attributes(aes, facet)
    link_axes!(aes; attr.linkxaxes, attr.linkyaxes)
    hideinnerdecorations!(aes, attr.hidexdecorations, attr.hideydecorations)

    # delete empty axes
    deleteemptyaxes!(aes)

    # add facet labels
    panel_labels!(fig, aes, scale)

    # span axis labels if appropriate
    nonempty_aes = get_nonempty_aes(aes)
    is2d = isaxis2d(first(nonempty_aes))

    if is2d && consistent_ylabels(aes)
        span_ylabel!(fig, aes)
    end
    if is2d && consistent_xlabels(aes)
        span_xlabel!(fig, aes)
    end

    return
end

function facet_grid!(fig, aes::AbstractMatrix{AxisEntries}; facet)
    M, N = size(aes)
    row_scale, col_scale = map(sym -> get(aes[1].categoricalscales, sym, nothing), (:row, :col))
    all(isnothing, (row_scale, col_scale)) && return

    # Link axes and hide decorations if appropriate
    attr = clean_facet_attributes(aes, facet)
    link_axes!(aes; attr.linkxaxes, attr.linkyaxes)
    hideinnerdecorations!(aes, attr.hidexdecorations, attr.hideydecorations)

    # span axis labels if appropriate
    nonempty_aes = get_nonempty_aes(aes)
    is2d = isaxis2d(first(nonempty_aes))

    if !isnothing(row_scale) && consistent_ylabels(aes)
        is2d && span_ylabel!(fig, aes)
        row_labels!(fig, aes, row_scale)
    end
    if !isnothing(col_scale) && consistent_xlabels(aes)
        is2d && span_xlabel!(fig, aes)
        col_labels!(fig, aes, col_scale)
    end
    return
end

function clean_facet_attributes(aes, facet)
    linkxaxes = get(facet, :linkxaxes, automatic)
    linkyaxes = get(facet, :linkyaxes, automatic)
    hidexdecorations = get(facet, :hidexdecorations, automatic)
    hideydecorations = get(facet, :hideydecorations, automatic)

    if linkxaxes ∉ [:all, :colwise, :minimal, :none, true, false, automatic]
        @warn "Replaced invalid keyword linkxaxes = $linkxaxes by automatic"
        linkxaxes = automatic
    end

    if linkyaxes ∉ [:all, :rowwise, :minimal, :none, true, false, automatic]
        @warn "Replaced invalid keyword linkyaxes = $linkyaxes by automatic"
        linkyaxes = automatic
    end

    if hidexdecorations ∉ [true, false, automatic]
        @warn "Replaced invalid keyword hidexdecorations = $hidexdecorations by automatic"
        hidexdecorations = automatic
    end

    if hideydecorations ∉ [true, false, automatic]
        @warn "Replaced invalid keyword hideydecorations = $hideydecorations by automatic"
        hideydecorations = automatic
    end

    if linkxaxes === automatic
        if consistent_xlabels(aes)
            linkxaxes = :all
        elseif colwise_consistent_xlabels(aes)
            linkxaxes = :colwise
        else
            linkxaxes = :none
        end
    end

    if linkxaxes == :minimal
        if colwise_consistent_xlabels(aes)
            linkxaxes = :colwise
        else
            linkxaxes = :none
        end
    end

    if linkyaxes === automatic
        if consistent_ylabels(aes)
            linkyaxes = :all
        elseif rowwise_consistent_ylabels(aes)
            linkyaxes = :rowwise
        else
            linkyaxes = :none
        end
    end

    if linkyaxes == :minimal
        if rowwise_consistent_ylabels(aes)
            linkyaxes = :rowwise
        else
            linkyaxes = :none
        end
    end

    if linkxaxes == true
        linkxaxes = :all
    elseif linkxaxes == false
        linkxaxes = :none
    end

    if linkyaxes == true
        linkyaxes = :all
    elseif linkyaxes == false
        linkyaxes = :none
    end

    if hidexdecorations === automatic
        hidexdecorations = (linkxaxes != :none)
    end
    if hideydecorations === automatic
        hideydecorations = (linkyaxes != :none)
    end

    (; linkxaxes, linkyaxes, hidexdecorations, hideydecorations)
end

# link axes

function link_axes!(aes; linkxaxes, linkyaxes)
    if linkxaxes == :all
        linkxaxes!(aes...)
    elseif linkxaxes == :colwise
        link_cols!(aes)
    end

    if linkyaxes == :all
        linkyaxes!(aes...)
    elseif linkyaxes == :rowwise
        link_rows!(aes)
    end
end

function link_rows!(aes)
    M, N = size(aes)
    for i in 1:M
        linkyaxes!(aes[i,:]...)
    end
end

function link_cols!(aes)
    M, N = size(aes)
    for i in 1:N
        linkxaxes!(aes[:,i]...)
    end
end

# facet labels

function col_labels!(fig, aes, scale)
    zipped_scale = zip(plotvalues(scale), datavalues(scale))

    titlegap, attributes = facetlabelattributes(first_nonempty_axis(aes))

    facetlabelpadding = lift(titlegap) do gap
        return (0f0, 0f0, gap, 0f0)
    end

    for (index, label) in zipped_scale
        Label(fig[1, index, Top()], to_string(label);
              padding=facetlabelpadding, attributes...)
    end
end

function row_labels!(fig, aes, scale)
    _, N = size(aes)
    zipped_scale = zip(plotvalues(scale), datavalues(scale))

    titlegap, attributes = facetlabelattributes(first_nonempty_axis(aes))

    facetlabelpadding = lift(titlegap) do gap
        return (gap, 0f0, 0f0, 0f0)
    end

    for (index, label) in zipped_scale
        Label(fig[index, N, Right()], to_string(label);
              rotation=-π/2, padding=facetlabelpadding,
              attributes...)
    end
end

function panel_labels!(fig, aes, scale)
    zipped_scale = zip(plotvalues(scale), datavalues(scale))

    titlegap, attributes = facetlabelattributes(first_nonempty_axis(aes))

    facetlabelpadding = lift(titlegap) do gap
        return (0f0, 0f0, gap, 0f0)
    end

    for (index, label) in zipped_scale
        Label(fig[index..., Top()], to_string(label);
              padding=facetlabelpadding, attributes...)
    end
end

function facetlabelattributes(ax)
    titlegap = ax.titlegap

    attributes = (
        color=ax.titlecolor,
        font=ax.titlefont,
        textsize=ax.titlesize,
    )

    (; titlegap, attributes)
end

# consistent axis labels

function consistent_xlabels(aes)
    nonempty_aes = get_nonempty_aes(aes)
    ax = first(nonempty_aes).axis
    return all(ae -> ae.axis.xlabel[] == ax.xlabel[], nonempty_aes)
end

function consistent_ylabels(aes)
    nonempty_aes = get_nonempty_aes(aes)
    ax = first(nonempty_aes).axis
    return all(ae -> ae.axis.ylabel[] == ax.ylabel[], nonempty_aes)
end

function colwise_consistent_xlabels(aes)
    _, N = size(aes)
    all(consistent_xlabels(aes[:,i]) for i in 1:N)
end

function rowwise_consistent_ylabels(aes)
    M, _ = size(aes)
    all(consistent_ylabels(aes[i,:]) for i in 1:M)
end

# spanned axis labels

function span_xlabel!(fig, aes)
    M, N = size(aes)

    for ae in aes
        ae.axis.xlabelvisible[] = false
    end
    protrusion = lift(
        (xs...) -> maximum(x -> x.bottom, xs),
        (MakieLayout.protrusionsobservable(ae.axis) for ae in aes[M, :])...
    )

    ax = first_nonempty_axis(aes)

    xlabelpadding = lift(protrusion, ax.xlabelpadding) do val, p
        return (0f0, 0f0, 0f0, val + p)
    end
    xlabelcolor = ax.xlabelcolor
    xlabelfont = ax.xlabelfont
    xlabelsize = ax.xlabelsize
    xlabelattributes = (
        color=xlabelcolor,
        font=xlabelfont,
        textsize=xlabelsize,
    )
    Label(fig[M, :, Bottom()], ax.xlabel;
        padding=xlabelpadding, xlabelattributes...)
end

function span_ylabel!(fig, aes)
    M, N = size(aes)

    for ae in aes
        ae.axis.ylabelvisible[] = false
    end

    protrusion = lift(
        (xs...) -> maximum(x -> x.left, xs),
        (MakieLayout.protrusionsobservable(ae.axis) for ae in aes[:, 1])...
    )
    # TODO: here and below, set in such a way that one can change padding after the fact?
    ax = first_nonempty_axis(aes)

    ylabelpadding = lift(protrusion, ax.ylabelpadding) do val, p
        return (0f0, val + p, 0f0, 0f0)
    end

    ylabelcolor = ax.ylabelcolor
    ylabelfont = ax.ylabelfont
    ylabelsize = ax.ylabelsize
    ylabelattributes = (
        color=ylabelcolor,
        font=ylabelfont,
        textsize=ylabelsize,
    )
    Label(fig[:, 1, Left()], ax.ylabel;
        rotation=π/2, padding=ylabelpadding, ylabelattributes...)
end

empty_ae(ae) = isempty(ae.entries)
get_nonempty_aes(aes) = filter(!empty_ae, aes)
first_nonempty_axis(aes) = first(get_nonempty_aes(aes)).axis

function facet!(fig, aes::AbstractMatrix{AxisEntries}; facet)
    update(fig) do f
        facet_wrap!(f, aes; facet)
        facet_grid!(f, aes; facet)
    end
    return
end

function facet!(fg::FigureGrid; facet)
    facet!(fg.figure, fg.grid; facet)
    return fg
end

## Layout helpers

isaxis2d(::Axis) = true
isaxis2d(::Axis3) = false
isaxis2d(ax::AxisSpec) = ax.type <: Axis
isaxis2d(ae::AxisEntries) = isaxis2d(ae.axis)
isaxis2d(ae::AxisSpecEntries) = isaxis2d(ae.axis)

for sym in [:hidexdecorations!, :hideydecorations!, :hidedecorations!]
    @eval function $sym(ae::AxisEntries; kwargs...)
        axis = ae.axis
        isaxis2d(axis) && $sym(axis; kwargs...)
    end
end

for sym in [:linkxaxes!, :linkyaxes!, :linkaxes!]
    @eval function $sym(ae::AxisEntries, aes::AxisEntries...)
        axs = filter(isaxis2d, map(ae->ae.axis, (ae, aes...)))
        isempty(axs) || $sym(axs...)
    end
end

function hideinnerdecorations!(aes, hidexdecorations, hideydecorations)
    I, J = size(aes)

    if hideydecorations
        hideydecorations!.(aes[:, 2:end], grid = false)
    end

    if hidexdecorations
        for i in 1:I, j in 1:J
            # don't hide x decorations if axis below is empty
            below_empty = (i < I) && empty_ae(aes[i+1,j])

            if (i < I) && !below_empty
                hidexdecorations!(aes[i,j],
                    grid = false,
                    ticks = hidexdecorations,
                    ticklabels = hidexdecorations
                )
            end

            if (i < I) && below_empty
                # improve appearance with empty axes
                aes[i,j].axis.alignmode = Mixed(bottom = MakieLayout.GridLayoutBase.Protrusion(0))
            end
        end
    end
end

function deleteemptyaxes!(aes::Matrix{AxisEntries})
    for ae in aes
        if isempty(ae.entries)
            delete!(ae.axis)
        end
    end
end

function resizetocontent!(fig::Figure)
    figsize = size(fig.scene)
    sz = map((Col(), Row()), figsize) do dir, currentsize
        inferredsize = determinedirsize(fig.layout, dir)
        return ceil(Int, something(inferredsize, currentsize))
    end
    sz == figsize || resize!(fig.scene, sz)
    return fig
end

function resizetocontent!(fg::FigureGrid)
    resizetocontent!(fg.figure)
    return fg
end


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

# example 0
true_d = MixtureModel(Normal, [(5.0, 2), (10.0, 2), (20.0, 1), (-2, 1)], [1.0/3, 1.0/4, 1.0/4, 1.0/6])
# true_d = Gamma(2.0, 3.0)
ll(x) = logpdf(true_d, x)

init_approx_d =  Normal(0, 40)
x = linspace(-10,30,500)
max_iter = 50

figure()
subplot(3, 2, 1)
plot(x, exp(ll(x)), label="true")
plot(x, pdf(init_approx_d, x), label="init approx")
for (col, method) in enumerate(["inclusive", "exclusive"])
    approx_d, trace = BVB.infer(ll, init_approx_d, method=method, max_iter=max_iter,
                                eval_cross_entropy=true, eval_ELBO=true, show_plot=true)
    subplot(3, 2, 1)
    plot(x, pdf(approx_d, x), label=method); legend(loc=0)

    subplot(3, 2, 2)
    plot(trace["iter"], trace["cross_entropy"], ".-", label=method)
    xlabel("x"); ylabel("cross entropy"); legend(loc=0); grid("on")

    subplot(3, 2, 3)
    plot(trace["iter"], trace["ELBO"], ".-", label=method)
    xlabel("x"); ylabel("ELBO"); legend(loc=0); grid("on")

    subplot(3, 2, 4)
    vlines(collect(1:max_iter) + 0.5 * (col-1), 0., trace["α"], label=method, colors=(col==1 ? "blue" : "green"))
    xlabel("iter"); ylabel("alpha"); legend(loc=0); grid("on")

    subplot(3, 2, 5)
    plot(1:length(probs(approx_d)), sort(probs(approx_d), rev=true), ".-", label=method)
    xlabel("component"); ylabel("weight"); legend(loc=0); grid("on")

    subplot(3, 2, 6)
    plot(x, ll(x) - logpdf(approx_d,x), "-", label=method)
    axhline(y=0., ls="--", color="red")
    xlabel("x"); ylabel("log L - log p"); legend(loc=0); grid("on")
end

show()

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
#=
Implement a queue using two stacks. Recall that a queue is a FIFO (first-in, first-out) data structure with the following methods: enqueue, which inserts an element into the queue, and dequeue, which removes it.
=#

using Test
include("Solutions/problem53_queue_from_two_stack.jl")

queue= Queue([],[])
enqueue(queue,1)
enqueue(queue,2)
enqueue(queue,3)
enqueue(queue,4)
enqueue(queue,5)

@test dequeue(queue) == 1
@test dequeue(queue) == 2
enqueue(queue, 6)
@test dequeue(queue) == 3

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module TestUtils
    using Base.Test
    using DataFrames

    macro testid(s::String, b::Bool)
        s = esc(s)
        b = esc(b)
        quote
            @test DataFrames.isidentifier($s) == $b
            @test DataFrames.is_valid_identifier(symbol($s)) == $b
        end
    end

    # 1-character

    @testid "a" true
    @testid "1" false
    @testid "!" false
    @testid "_" true

    # 2-characters, initial alphabetical character

    @testid "aa" true
    @testid "a1" true
    @testid "a!" true
    @testid "a_" true

    # 2-characters, initial underscore

    @testid "_a" true
    @testid "_1" true
    @testid "_!" true
    @testid "__" true

    # 3-characters, initial alphabetical character

    @testid "aaa" true
    @testid "aa1" true
    @testid "aa!" true
    @testid "aa_" true

    @testid "a1a" true
    @testid "a11" true
    @testid "a1!" true
    @testid "a1_" true

    @testid "a!a" true
    @testid "a!1" true
    @testid "a!!" true
    @testid "a!_" true

    @testid "a_a" true
    @testid "a_1" true
    @testid "a_!" true
    @testid "a__" true

    # 3-characters, initial underscore

    @testid "_aa" true
    @testid "_a1" true
    @testid "_a!" true
    @testid "_a_" true

    @testid "_1a" true
    @testid "_11" true
    @testid "_1!" true
    @testid "_1_" true

    @testid "_!a" true
    @testid "_!1" true
    @testid "_!!" true
    @testid "_!_" true

    @testid "__a" true
    @testid "__1" true
    @testid "__!" true
    @testid "___" true

    # False

    @testid "1a" false
    @testid "!a" false
    @testid "a*" false
    @testid "begin" false
    @testid "end" false

    @test DataFrames.is_valid_identifier(symbol("a\u212b")) == false

    @test DataFrames.makeidentifier("%_B*\tC*") == "B_C"
    @test DataFrames.makeidentifier("2a") == "x2a"
    @test DataFrames.makeidentifier("!") == "x!"
    @test DataFrames.makeidentifier("\t_*") == "x"
    @test DataFrames.makeidentifier("begin") == "_begin"
    @test DataFrames.makeidentifier("end") == "_end"

    @test DataFrames.make_unique([:x, :x, :x_1, :x2]) == [:x, :x_2, :x_1, :x2]

    # Check that reserved words are up to date
    f = "$JULIA_HOME/../../src/julia-parser.scm"
    if isfile(f)
        r1 = r"define reserved-words '\(([^)]+)"
        r2 = r"define \(parse-block s\) \(parse-Nary s parse-eq '\([^(]+\(([^)]+)"
        body = readall(f)
        m1, m2 = match(r1, body), match(r2, body)
        if m1 == nothing || m2 == nothing
            error("Unable to extract keywords from 'julia-parser.scm'.")
        else
            rw = split(m1.captures[1]*" "*m2.captures[1], r"\W+")
            @test rw == DataFrames.RESERVED_WORDS
        end
    else
        warn("Unable to find validate reserved words against parser. ",
             "Expected if Julia was not built from source.")
    end
end


text:
mutable struct MongoCursor
    _wrap_::Ptr{Void}

    MongoCursor(_wrap_::Ptr{Void}) = begin
        cursor = new(_wrap_)
        finalizer(cursor, destroy)
        return cursor
    end
end
export MongoCursor

# Iterator

Base.start(cursor::MongoCursor) = nothing

Base.next(cursor::MongoCursor, state::Void) =
    (BSONObject(ccall(
        (:mongoc_cursor_current, libmongoc),
        Ptr{Void}, (Ptr{Void},),
        cursor._wrap_
        ), Union{}), state)

function Base.done(cursor::MongoCursor, state::Void)
    return !ccall(
        (:mongoc_cursor_next, libmongoc),
        Bool, (Ptr{Void}, Ptr{Ptr{Void}}),
        cursor._wrap_,
        Array{Ptr{Void}}(1)
        )
end

Base.iteratorsize(::Type{MongoCursor}) = Base.SizeUnknown()
Base.eltype(::Type{MongoCursor}) = BSONObject

destroy(collection::MongoCursor) =
    ccall(
        (:mongoc_cursor_destroy, libmongoc),
        Void, (Ptr{Void},),
        collection._wrap_
        )


text:
# todo
# paddedSequence(x...) = Sequence(trim.(x)...)
trimhv(x) = CombinedParsers.trim(x;whitespace=Atomic(Repeat(CharIn(horizontal_space_char,vertical_space_char))))

@generated function separatedSequence(x...; whitespace=Atomic(Repeat(CharIn(horizontal_space_char,vertical_space_char))))
    f = Expr(Symbol("->"), :v, Expr(:tuple, Any[Expr(:ref, :v, i) for i in 2*1:length(x)]))
    a = Iterators.flatten([ [:i, QuoteNode(whitespace)] for i in x ])
    quote
        (
        # map($f,
            Sequence(whitespace, $a...))
    end
end

function separatedTriple(f::Function, open,body,close; whitespace=Atomic(Repeat(CharIn(horizontal_space_char,vertical_space_char))))
    Sequence(whitespace, open, whitespace, body, whitespace, close, whitespace) do v
        f(v[4])
    end
end

# https://upload.wikimedia.org/wikipedia/commons/0/0c/Ebnf-syntax-diagram.png
begin
    letter                              = CharIn('A':'Z','a':'z')
    decimal_digit                       = CharIn('0':'9')
    concatenate_symbol                  = ','
    defining_symbol                     = '='
    definition_separator_symbol         = CharIn("|/!")
    start_comment_symbol                = "(*"
    start_group_symbol                  = "("
    end_comment_symbol                  = "*)"
    end_group_symbol                    = ")"
    end_option_symbol                   = Either(']', "/)")
    end_repeat_symbol                   = Either('}', ":)")

    except_symbol                       = '-'
    first_quote_symbol                  = '''
    repetition_symbol                   = '*'
    second_quote_symbol                 = '"'
    special_sequence_symbol             = '?'

    start_option_symbol                 = Either('[', "(/")
    start_repeat_symbol                 = Either('{', "(:")
    terminator_symbol                   = CharIn(";.")

    other_character                     = CharIn(" :+_%@&#\$<>\\^`~") # is it `?
    space_character                     = ' '
    horizontal_tabulation_character     = horizontal_space_char#  "\t" #^-?
    new_line                            = '\n' #?
    vertical_tabulation_character       = vertical_space_char# Never() #?
    form_feed                           = '\f' # Never() #?

    terminal_character                  = !Atomic(Either(
        letter,
        decimal_digit,
        concatenate_symbol,
        defining_symbol,
        definition_separator_symbol,
        end_comment_symbol,
        end_group_symbol,
        end_option_symbol,
        end_repeat_symbol,
        except_symbol,
        first_quote_symbol,
        repetition_symbol,
        second_quote_symbol,
        special_sequence_symbol,
        start_comment_symbol,
        start_group_symbol,
        start_option_symbol,
        start_repeat_symbol,
        terminator_symbol,
        other_character
    ))

    first_terminal_character            = Sequence(NegativeLookahead(first_quote_symbol),
                                        terminal_character)

    second_terminal_character           = Sequence(NegativeLookahead(second_quote_symbol),
                                         terminal_character)

    @with_names terminal_string         = Either(Sequence(2,first_quote_symbol,
                                                  !!Repeat1(first_terminal_character),
                                                  first_quote_symbol),
                                         Sequence(2,second_quote_symbol,
                                                  !!Repeat1(second_terminal_character),
                                                  second_quote_symbol))

    gap_separator                       = CharIn(space_character, horizontal_tabulation_character,
                           new_line, vertical_tabulation_character, form_feed)

    @with_names _integer                = CombinedParsers.Numeric(Int) #!Repeat1(decimal_digit)

    meta_identifier_character           = CharIn(letter, decimal_digit, "-_ ") # optimize!
    @with_names meta_identifier         = map(Symbol,!join(!Sequence(letter, Lazy(Repeat(meta_identifier_character))), gap_separator))
    meta_identifier("abc")

    special_sequence_character          = Sequence(NegativeLookahead(special_sequence_symbol),
                                          terminal_character)
    @with_names special_sequence        = Sequence(2,
                                                   special_sequence_symbol, trimhv(!!Lazy(Repeat(special_sequence_character))), special_sequence_symbol)
    special_sequence("? ASCII char 32 ?")

    commentless_symbol                  = Either(
        Sequence(
            NegativeLookahead(Either(
                letter,
                decimal_digit,
                first_quote_symbol,
                second_quote_symbol,
                start_comment_symbol,
                end_comment_symbol,
                special_sequence_symbol,
                other_character)),
            terminal_character
        ),
        meta_identifier,
        _integer,
        terminal_string,
        special_sequence
    )

    comment_symbol                      = Either{Any}(Any[
        other_character,
        commentless_symbol
    ])
    @with_names bracket_textual_comment = Sequence(start_comment_symbol, !Repeat(comment_symbol), end_comment_symbol)
    pushfirst!(comment_symbol, bracket_textual_comment)

    empty_sequence                      = parser(Always() => Always())

    struct Substitution<:CombinedParser{Nothing,Nothing}
        name::Symbol
    end
    syntactic_primary                   = Either{CombinedParser}(
        Any[map(n->Substitution(n), meta_identifier),
            map(parser,terminal_string),
            map(s->Never(),special_sequence),
            empty_sequence])

    syntactic_primary("\"1b\"")

    syntactic_factor                    = Sequence(
        Optional(
            Sequence(1, _integer, trimhv(repetition_symbol));
        default                         = 1),
        syntactic_primary) do v
            v[1]                        == 1 ? v[2] : Repeat(v[1], v[2])
        end

    syntactic_factor("\"1b\"")

    syntactic_exception                 = syntactic_factor
    syntactic_term                      = Sequence(1, syntactic_factor,
                              # todo: handle exceptions
                              Optional(Sequence(except_symbol, syntactic_exception)))

    @with_names single_definition       = map(p -> sSequence(p...)::CombinedParser,
                                        join(syntactic_term, trimhv(concatenate_symbol)))

    definitions_list                    = map(p -> sEither(p...)::CombinedParser,
                           join(single_definition, trimhv(definition_separator_symbol)))

    definitions_list("3 * \"aa\", \"B\"")

    @syntax for optional_sequence in syntactic_primary
        separatedTriple(start_option_symbol, definitions_list, end_option_symbol) do v
            Optional(v)
        end
    end
    @syntax for repeated_sequence in syntactic_primary
        separatedTriple(start_repeat_symbol, definitions_list, end_repeat_symbol) do v
            Repeat(v)
        end
    end
    @syntax for grouped_sequence in syntactic_primary
        separatedTriple(identity, start_group_symbol, definitions_list, end_group_symbol)
    end

    @syntax syntax_rule                         = Sequence(meta_identifier, trimhv(defining_symbol), definitions_list, trimhv(terminator_symbol)) do v
        with_name(v[1], v[3])
    end;
    match(syntax_rule,"twelve                 = \"1\", \"2\" ;")

    syntax                              = Repeat1(syntax_rule);
end;

using BenchmarkTools
@btime match(syntax,
    """
    digit excluding zero = "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9" ;
    digit                = "0" | digit excluding zero ;
    natural number = digit excluding zero, { digit } ;
    integer = "0" | [ "-" ], natural number ;
    """);

syntax(
    """
    twelve a                          = "1", "2" ;
    """)

syntax(
    """
    twelve                          = "1", "2" ;
    """, log=true)

syntax(
    """
    twelve                          = "1", "2" ;    two hundredone                 = "2", "0", "1" ;
    """, log=true)


text:
## Canonical Form of Normal distribution

immutable NormalCanon <: ContinuousUnivariateDistribution
    η::Float64       # σ^(-2) * μ
    λ::Float64    # σ^(-2)
    μ::Float64       # μ

    function NormalCanon(η::Real, λ::Real)
        @check_args(NormalCanon, λ > zero(λ))
    	new(η, λ, η / λ)
    end
    NormalCanon() = new(0.0, 1.0, 0.0)
end

@distr_support NormalCanon -Inf Inf

## conversion between Normal and NormalCanon

Base.convert(::Type{Normal}, d::NormalCanon) = Normal(d.μ, 1.0 / sqrt(d.λ))
Base.convert(::Type{NormalCanon}, d::Normal) = (λ = 1.0 / σ^2; NormalCanon(λ * d.μ, λ))
canonform(d::Normal) = convert(NormalCanon, d)

#### Parameters

params(d::NormalCanon) = (d.η, d.λ)

#### Statistics

mean(d::NormalCanon) = d.μ
median(d::NormalCanon) = mean(d)
mode(d::NormalCanon) = mean(d)

skewness(d::NormalCanon) = 0.0
kurtosis(d::NormalCanon) = 0.0

var(d::NormalCanon) = 1.0 / d.λ
std(d::NormalCanon) = sqrt(var(d))

entropy(d::NormalCanon) = 0.5 * (log2π + 1.0 - log(d.λ))

#### Evaluation

pdf(d::NormalCanon, x::Float64) = (sqrt(d.λ) / sqrt2π) * exp(-0.5 * d.λ * abs2(x - d.μ))
logpdf(d::NormalCanon, x::Float64) = 0.5 * (log(d.λ) - log2π - d.λ * abs2(x - d.μ))

zval(d::NormalCanon, x::Float64) = (x - d.μ) * sqrt(d.λ)
xval(d::NormalCanon, z::Float64) = d.μ + z / sqrt(d.λ)

cdf(d::NormalCanon, x::Float64) = normcdf(zval(d,x))
ccdf(d::NormalCanon, x::Float64) = normccdf(zval(d,x))
logcdf(d::NormalCanon, x::Float64) = normlogcdf(zval(d,x))
logccdf(d::NormalCanon, x::Float64) = normlogccdf(zval(d,x))

quantile(d::NormalCanon, p::Float64) = xval(d, norminvcdf(p))
cquantile(d::NormalCanon, p::Float64) = xval(d, norminvccdf(p))
invlogcdf(d::NormalCanon, lp::Float64) = xval(d, norminvlogcdf(lp))
invlogccdf(d::NormalCanon, lp::Float64) = xval(d, norminvlogccdf(lp))

#### Sampling

rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.λ)
rand!{T<:Real}(cf::NormalCanon, r::AbstractArray{T}) = rand!(convert(Normal, cf), r)


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

mutable struct ManagementPolicyRule <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    definition::Any # spec type: Union{ Nothing, ManagementPolicyDefinition } # spec name: definition

    function ManagementPolicyRule(;enabled=nothing, name=nothing, type=nothing, definition=nothing)
        o = new()
        validate_property(ManagementPolicyRule, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(ManagementPolicyRule, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ManagementPolicyRule, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(ManagementPolicyRule, Symbol("definition"), definition)
        setfield!(o, Symbol("definition"), definition)
        o
    end
end # type ManagementPolicyRule

const _property_map_ManagementPolicyRule = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("name")=>Symbol("name"), Symbol("type")=>Symbol("type"), Symbol("definition")=>Symbol("definition"))
const _property_types_ManagementPolicyRule = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("name")=>"String", Symbol("type")=>"String", Symbol("definition")=>"ManagementPolicyDefinition")
Base.propertynames(::Type{ ManagementPolicyRule }) = collect(keys(_property_map_ManagementPolicyRule))
Swagger.property_type(::Type{ ManagementPolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicyRule[name]))}
Swagger.field_name(::Type{ ManagementPolicyRule }, property_name::Symbol) =  _property_map_ManagementPolicyRule[property_name]

const _allowed_ManagementPolicyRule_type = ["Lifecycle"]

function check_required(o::ManagementPolicyRule)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("definition")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ManagementPolicyRule }, name::Symbol, val)
    if name === Symbol("type")
        Swagger.validate_param(name, "ManagementPolicyRule", :enum, val, _allowed_ManagementPolicyRule_type)
    end
end