text:
"""
    MLP(;M=4, L=4,)

Multi level Picard algorithm.

# Arguments
* `L`: number of Picard iterations (Level),
* `M`: number of Monte Carlo integrations (at each level `l`, `M^(L-l)`integrations),
"""
struct MLP{T, MCS} <: HighDimPDEAlgorithm where {T <: Int, MCS <: MCSampling}
    M::T # nb of MC integrations
    L::T # nb of levels
    K::T # nb MC integration non local term
    mc_sample!::MCS
end

MLP(; M=4, L=4, K=10, mc_sample = NoSampling()) = MLP(M ,L, K, mc_sample) #Note: mc_sample mutates its first argument but for the user interface we hide this technicality

# function DiffEqBase.__solve(
function solve(
        prob::PIDEProblem,
        alg::MLP;
        multithreading=true,
        verbose=false,
        )

    # unbin stuff
    x = prob.x
    neumann_bc = prob.neumann_bc
    K = alg.K
    M = alg.M
    L = alg.L
    mc_sample! = alg.mc_sample!
    g, f = prob.g, prob.f

    # errors
    !isnothing(prob.u_domain) ? error("`MLP` algorithm cannot be solved on a domain, i.e with argument `u_domain`.") : nothing
    isnothing(x) ? error("`MLP` algorithm needs a grid 'x'") : nothing
    # if !isnothing(neumann_bc)
    #     if !(typeof(mc_sample!) <: UniformSampling || typeof(mc_sample!) <: NoSampling)
    #         error("Only `UniformSampling` or  is covered with Neumann BC")
    #     end
    # end
    if multithreading
        usol = _ml_picard_mlt(M, L, K, x, prob.tspan[1], prob.tspan[2], mc_sample!, g, f, verbose, prob, neumann_bc)
    else
        usol = _ml_picard(M, L, K, x, prob.tspan[1], prob.tspan[2], mc_sample!, g, f, verbose, prob, neumann_bc)
    end
    return x, prob.tspan, [g(x),usol]
    # sol = DiffEqBase.build_solution(prob,alg,ts,usol)
    # save_everystep ? iters : u0(X0)[1]

end

function _ml_picard(M, # monte carlo integration
                    L, # level
                    K, # non local term monte carlo
                    x::xType, # initial point
                    s, # time
                    t::tType, # time
                    mc_sample!,
                    g,
                    f,
                    verbose,
                    prob,
                    neumann_bc
                    ) where {xType, tType}

    elxType = eltype(xType)
    if L == 0
        return zero(elxType)
    end

    x2 = similar(x)
    _integrate(mc_sample!) ? x3 = similar(x) : x3 = nothing
    p = prob.p

    a = zero(elxType)
    for l in 0:(min(L- 1, 1))
        verbose && println("loop l")
        b = zero(elxType)
        num = M^(L - l) # ? why 0.5 in sebastian code?
        for k in 1:num
            verbose && println("loop k")
            r = s + (t - s) * rand(tType)
            _mlt_sde_loop!(x2, x, s, r, prob, neumann_bc)
            b2 = _ml_picard(M, l, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b3 = zero(elxType)
            # non local integration
            for h in 1:K
                verbose && println("loop h")
                mc_sample!(x3, x2)
                b3 += f(x2, x3, b2, _ml_picard(M, l, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) #TODO:hardcode, not sure about t
            end
            b += b3 / K
        end
        a += (t - s) * b / num
    end

    for l in 2:(L-1)
        b = zero(elxType)
        num = M^(L - l)
        for k in 1:num
            r = s + (t - s) * rand(tType)
            _mlt_sde_loop!(x2, x, s, r, prob, neumann_bc)
            b2 = _ml_picard(M, l, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b4 = _ml_picard(M, l - 1, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b3 = zero(elxType)
                # non local integration
            for h in 1:K
                mc_sample!(x3, x2)
                b3 += f(x2, x3, b2, _ml_picard(M, l, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) - f(x2, x3, b4, _ml_picard(M, l - 1, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) #TODO:hardcode, not sure about t
            end
            b += b3 / K
        end
        a += (t - s) * b / num
    end

    num = M^(L)
    a2 = zero(elxType)
    for k in 1:num
        verbose && println("loop k3")
        _mlt_sde_loop!(x2, x, s, t, prob, neumann_bc)
        a2 += g(x2)
    end
    a2 /= num

    return a + a2
end

_ml_picard(M::Int, L::Int, K::Int, x::Nothing, s::Real, t::Real, mc_sample!, g, f, verbose::Bool, prob, neumann_bc) = nothing

function _ml_picard_mlt(M, # monte carlo integration
                        L, # level
                        K, # non local term monte carlo
                        x::xType, # initial point
                        s, # time
                        t, # time
                        mc_sample!,
                        g,
                        f,
                        verbose,
                        prob,
                        neumann_bc) where {xType}

    elxType = eltype(xType)

    # distributing tasks
    NUM_THREADS = Threads.nthreads()
    tasks = [Threads.@spawn(_ml_picard_call(M, L, K, x, s, t, mc_sample!, g, f, verbose, NUM_THREADS, thread_id, prob, neumann_bc)) for thread_id in 1:NUM_THREADS]

    # first level
    num = M^(L)
    x2 = similar(x)
    a2 = zero(elxType)
    for k in 1:num
        verbose && println("loop k3")
        _mlt_sde_loop!(x2, x, s, t, prob, neumann_bc)
        a2 += g(x2)
    end
    a2 /= num

    # fetching tasks
    a = sum([fetch(t) for t in tasks])

    return a + a2
end

function _ml_picard_call(M, # monte carlo integration
                        L, # level
                        K, # non local term monte carlo
                        x::xType, # initial point
                        s, # time
                        t::tType, # time
                        mc_sample!,
                        g,
                        f,
                        verbose,
                        NUM_THREADS,
                        thread_id,
                        prob,
                        neumann_bc
                        ) where {xType, tType}

    x2 = similar(x)
    _integrate(mc_sample!) ? x3 = similar(x) : x3 = nothing
    p = prob.p
    elxType = eltype(xType)

    a = zero(elxType)
    for l in 0:(min(L - 1, 1))
        b = zero(elxType)
        num = M^(L - l)
        loop_num = _get_loop_num(M, num, thread_id, NUM_THREADS)
        for k in 1:loop_num
            verbose && println("loop k")
            r = s + (t - s) * rand(tType)
            _mlt_sde_loop!(x2, x, s, r, prob, neumann_bc)
            b2 = _ml_picard(M, l, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b3 = zero(elxType)
            for h in 1:K # non local integration
                verbose && println("loop h")
                mc_sample!(x3, x2)
                b3 += f(x2, x3, b2, _ml_picard(M, l, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) #TODO:hardcode, not sure about t
            end
        b += b3 / K
        end
        a += (t - s) * b / num
    end

    for l in 2:(L-1)
        b = zero(elxType)
        num = M^(L - l)
        loop_num = _get_loop_num(M, num, thread_id, NUM_THREADS)
        for k in 1:loop_num
            r = s + (t - s) * rand(tType)
            _mlt_sde_loop!(x2, x, s, r, prob, neumann_bc)
            b2 = _ml_picard(M, l, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b4 = _ml_picard(M, l - 1, K, x2, r, t, mc_sample!, g, f, verbose, prob, neumann_bc)
            b3 = zero(elxType)
                # non local integration
            for h in 1:K
                mc_sample!(x3, x2)
                b3 += f(x2, x3, b2, _ml_picard(M, l, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) - f(x2, x3, b4, _ml_picard(M, l - 1, K, x3, r, t, mc_sample!, g, f, verbose, prob, neumann_bc), nothing, nothing, p, t) #TODO:hardcode, not sure about t
            end
        b += b3 / K
        end
        a += (t - s) * b / num
    end

    return a

end

function _get_loop_num(M, num, thread_id, NUM_THREADS) #decides how many iteration given thread id and num
    if num < NUM_THREADS
        # each thread only goes once through the loop
        loop_num = thread_id > num ? 0 : 1
    else
        remainder =  M % num
        if (remainder > 0) && (thread_id <= remainder) # remainder > 0 iff num == M or num == 1
            # each thread goes  num / NUM_THREADS + the remainder
            loop_num = num / NUM_THREADS + 1
        else
            loop_num = num / NUM_THREADS
        end
    end
end

function _mlt_sde_loop!(x2,
                        x,
                        s,
                        t,
                        prob,
                        neumann_bc)
    #randn! allows to save one allocation
    dt = t - s
    randn!(x2)
    x2 .= x + (prob.μ(x, prob.p, t) .* dt .+ prob.σ(x, prob.p, t) .* sqrt(dt) .* x2)
    if !isnothing(neumann_bc)
        x2 .= _reflect(x, x2, neumann_bc...)
    end
end


text:
using Avalon
import Avalon: batchiter
using MLDatasets
using BenchmarkTools

mutable struct Net
    conv1::Conv2d
    conv2::Conv2d
    conv3::Conv2d
    fc::Linear
end

Net() = Net(
    Conv2d(1, 16, 3; padding=(1,1)),
    Conv2d(16, 32, 3; padding=(1,1)),
    Conv2d(32, 32, 3; padding=(1,1)),
    Linear(288, 10)
)

function (m::Net)(x::AbstractArray)
    x = maxpool2d(relu.(m.conv1(x)), (2, 2))
    x = maxpool2d(relu.(m.conv2(x)), (2, 2))
    x = maxpool2d(relu.(m.conv3(x)), (2, 2))
    x = reshape(x, 288, :)
    x = relu.(m.fc(x))
    x = softmax(x)
    return x
end

function get_mnist_data(train::Bool; device=CPU())
    X, Y = train ? MNIST.traindata() : MNIST.testdata()
    X = convert(Array{Float64}, reshape(X, 28, 28, 1, :)) |> device
    # replace class label like "0" with its position like "1"
    Y = Y .+ 1
    # transform to one-hot encoding
    Y1h = zeros(10, length(Y))
    for (j, y) in enumerate(Y)
        Y1h[y, j] = 1.0
    end
    Y1h = Y1h |> device
    return X, Y1h
end

function main(device=best_available_device())
    # instantiate the model
    m = Net() |> device
    # load training data
    X_trn, Y_trn = get_mnist_data(true);
    X_trn = device(X_trn)
    # set loss function and optimizer, then fit the model
    loss_fn = CrossEntropyLoss()
    opt = Adam(lr=1e-3)
    @time fit!(m, X_trn, Y_trn, loss_fn; n_epochs=10, opt=opt, batch_size=100, device=device)

    x = X_trn[:,:,:,1:100] |> device
    @benchmark m(x)
end

# function main2()
#     x = X_trn[:,:,:,1:100] |> device
#     y = Y_trn[:, 1:100] |> device

#     _, g = grad((m, x, y) -> loss_fn(m(x), y), m, x, y)
# end

# using Yota
# import Avalon: batchiter

# function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat)
#   return -sum(y .* log.(ŷ)) * 1 // size(y, 2)
# end

function Avalon.fit!(m::Net, X::AbstractArray, Y::AbstractArray, loss_fn;
              n_epochs=10, batch_size=100, opt=SGD(1e-3), device=CPU(), report_every=1)
    f = (m, x, y) -> loss_fn(m(x), y)
    num_batches = size(X)[end] // batch_size
    for epoch in 1:n_epochs
        epoch_loss = 0
        @time for (i, (x, y)) in enumerate(batchiter((X, Y), sz=batch_size))
            x = to_device(device, copy(x))
            y = to_device(device, copy(y))
            loss, g = grad(f, m, x, y)
            update!(opt, m, g[2])
            epoch_loss += loss
        end
        if epoch % report_every == 0
            println("Epoch $epoch: avg_loss=$(epoch_loss / num_batches)")
        end
    end
    return m
end

# using Printf
# using Yota: Tape, exec!, play!

# function profile_tape!(tape::Tape, args...; use_compiled=true, debug=false)
#     for (i, val) in enumerate(args)
#         @assert(tape[i] isa Input, "More arguments than the original function had")
#         tape[i].val = val
#     end
#     timings = Vector{Float64}(undef, 0)
#     if use_compiled && tape.compiled != nothing
#         Base.invokelatest(tape.compiled)
#     else
#         for op in tape
#             if debug
#                 println(op)
#             end
#             t = @elapsed exec!(tape, op)
#             push!(timings, t)
#         end
#     end
#     for (t, op) in zip(timings, tape)
#         @printf "%.3f\t%s\n" (t * 10^6) op
#     end
#     return tape[tape.resultid].val
# end

ext: jl


text:
module Eive

include("cga.jl")
include("estimator.jl")

export CGA
export Estimator

import .Estimator: eive

export eive

end # module


text:
using go
using ArgParse

# Some flags are used for different purposes for h5 and loading sgf
function parseargs()
    s = ArgParseSettings()
    @add_arg_table s begin
        "model_name"
          help = "Model name (saved in models dir)"
          arg_type = AbstractString
          required = true
        "--h5path"
          help = "Path to h5 file to load training data from"
          arg_type = AbstractString
        "--kgs"
          help = "Load data from kgs sgf files"
          action = :store_true
        "--cgs"
          help = "Load data from cgs sgf files"
          action = :store_true
        "--rating"
          help = "Minimum elo rating (for cgs)"
          default = 2500
          arg_type = Int
        "--epochs"
          help = "Number training epochs"
          arg_type = Int
          default = 10
        "--filename_substring"
          help = "Filter down to sgf files which contain this substring"
          arg_type = AbstractString
          default = ".sgf"
        "--n_train"
          help = "Maximum number of sgf files to load OR how many moves to load from h5 file"
          arg_type = Int
          default = 10000000  # Basically no limit
        "--n_validation"
          help = "How many moves to use for validation in h5"
          arg_type = Int
          default = 10000
        "--load"
          help = "Load weights when trained model already exists"
          action = :store_true
    end
    parse_args(s)
end

pargs = parseargs()

# Set the feature extractors
feats = go.LIBERTIES

LOAD_H5 = pargs["h5path"] != nothing

if !LOAD_H5
    cgs = pargs["cgs"]
    kgs = pargs["kgs"]
    if (kgs + cgs) != 1
        println(STDERR, "Exactly one of h5path, cgs, or kgs must be set")
        exit()
    end
    DATASOURCE = cgs ? "cgs" : "kgs"
    println(STDERR, "Finding SGF files...")
    files = go.find_sgf("../data/$(DATASOURCE)/");
    println(STDERR, "Found $(length(files)) sgf files")

    # Filter to filenames containing regex (e.g. month/year)
    if "filename_substring" in keys(pargs)
        matcher = Regex(pargs["filename_substring"])
        filtered = Vector{AbstractString}()
        for filename in files
            if match(matcher, filename) != nothing
                push!(filtered, filename)
            end
        end
        files = filtered
        println(STDERR, "Contains filtered down to $(length(files))")
    end

    if cgs
        println(STDERR, "Finding SGF games with highly ranked players...")
        files =  filter(x -> go.players_over_rating(x, pargs["rating"]), files)
    end
    # if kgs...

    println("Found $(length(files)) matching sgf files")

    files = files[randperm(length(files))][1:min(end, pargs["n_train"])]

    println(STDERR, "Extracting features from $(length(files)) games...")
    examples = go.generate_training_data(files; progress_update=500, features=feats)
    X = go.batch_training_examples([x[1] for x in examples])
    Y = go.batch_training_examples([x[2][:] for x in examples])
end

println(STDERR, "Initializing model...")
clf = go.CLF_LINEAR(feats)

println(STDERR, "Training...")
tm = time()
MODEL_NAME = pargs["model_name"]
checkpoint_path = "../models/$(MODEL_NAME)_checkpoint.hf5"
go.save_model(clf, "../models/", MODEL_NAME, save_weights=false, save_yaml=true)

if pargs["load"]
    model.load_weights("../models/$(MODEL_NAME).hf5")
end

if LOAD_H5
    go.keras_train_h5(clf, pargs["h5path"],
                      epochs=pargs["epochs"],
                      n_train=pargs["n_train"],
                      n_validation=pargs["n_validation"],
                      checkpoint_path=checkpoint_path)
else
    go.train_model(clf, X, Y, epochs=pargs["epochs"], validation_split=0.05, checkpoint_path=checkpoint_path)
end
println("Took $(time() - tm) seconds")

println(STDERR, "Saving model...")
go.save_model(clf, "../models/", MODEL_NAME, save_weights=true, save_yaml=false, overwrite=pargs["load"])


text:
#Pkg.add("HyperbolicPlane")

using HyperbolicPlane, Plots

"""
`pencil(s::Real, nlines::Int, ncycles::Int)` draws a pencil of `nlines`
lines with common ideal point `exp(im*s)` and `ncycles` Horocycles that
are orthogonal to those lines.
"""
function pencil(s::Real=1, nlines::Int=5, ncycles::Int=5)
    Plots.plot()

    # draw a pencil of lines
    tlist = [ 2*pi*k/(nlines+1) for k=1:nlines ]
    for t in tlist
        L = HLine(s, s+t)
        set_color(L,:navyblue)
        draw(L)
    end

    # draw a pencil of horocycles
    xlist = [ 2*t/(ncycles+1) - 1 for t in 1:ncycles]

    plist = [ HPoint(x+0im) for x in xlist ]
    f = rotation(s)
    for P in plist
        H = Horocycle(f(P), s)
        set_color(H,:royalblue)
        draw(H)
    end
    draw(HPlane())
    finish()
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using ChemistryFeaturization
using Test

const testdir = dirname(@__FILE__)

tests = [
    "module_tests",
    "utils/ElementFeatureUtils_tests",
    "utils/GraphBuilding_tests",
    "atoms/AtomGraph_tests",
    "features/ElementFeature_tests",
    "features/OrbitalFeature_tests",
    "utils/OrbitalFeatureUtils_tests",
    "featurizations/GraphNodeFeaturization_tests",
    "featurizedatoms_tests",
    # TODO: add Weave stuff
    # TODO: add SpeciesFeature tests
]

@testset "ChemistryFeaturization" begin
    for t in tests
        tp = abspath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
# Examples from Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis

using IntervalArithmetic, IntervalRootFinding, StaticArrays

A = [4..5 -1..1 1.5..2.5; -0.5..0.5 -7.. -5 1..2; -1.5.. -0.5 -0.7.. -0.5 2..3]
sA = SMatrix{3}{3}(A)
mA = MMatrix{3}{3}(A)

b = [3..4, 0..2, 3..4]
sb = SVector{3}(b)
mb = MVector{3}(b)

p = fill(-1e16..1e16, 3)

rts = gauss_seidel_interval!(p, A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
rts = gauss_seidel_interval!(p, sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
rts = gauss_seidel_interval!(p, mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default

rts = gauss_seidel_interval(A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
rts = gauss_seidel_interval(sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
rts = gauss_seidel_interval(mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default

rts = gauss_seidel_contractor!(p, A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
rts = gauss_seidel_contractor!(p, sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
rts = gauss_seidel_contractor!(p, mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default

rts = gauss_seidel_contractor(A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
rts = gauss_seidel_contractor(sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
rts = gauss_seidel_contractor(mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default

rts = gauss_elimination_interval!(p, A, b, precondition=true) # Gaussian Elimination; precondition=true by default
rts = gauss_elimination_interval(A, b, precondition=true) # Gaussian Elimination; precondition=true by default


text:
@testset "type level" begin
    @case S{A, B}(a :: A, b :: B)
    s = S(1, "2")
    @test @match s begin
        ::S{Int, String} => true
    end

    @testset "type level for func pattern" begin
        f = (Int ⇒ Int)(x -> x)
        Feature.@activate TypeLevel
        function()
                 @match f begin
                ::(T ⇒ G)  => T == G
                _          => false
                end
        end |> it -> @test it()

    end

    @testset "type level for type parameter destruction" begin
        @case S{A, B}(a :: A, b :: B)
        s = S(1, "2")
        function ()
            @match s begin
                ::Kind{A, B} => Kind === S && A === Int && B === String
            end
        end |> it -> @test it()
    end

end


text:
"""
    `linearize(t::NewType) -> Vector{NewType}`

Return ordered list of all types, which type `t` can inherit from.
The fist element of the list is `Any`, the last element is `t`.
If `i` <= `j` then list[j] isa list[i]
"""

function linearize(x::NewType)
    _linearize(NewType[], x)
end
function _linearize(list::Vector{NewType}, x::NewType)
    if x ∉ list
        for y in supertypes(x)
            _linearize(list, y)
        end
        push_new!(list, x, is_pushed(typeof(x)))
    end
    list
end

export linearize

is_pushed(::Type) = false
is_pushed(::Type{NewDataType}) = true
is_pushed(::Type{NewUnionAll}) = true

function push_new!(x::Vector{NewType}, y::NewType, doit::Bool)
    if doit && y ∉ x && y != x
        push!(x, y)
    end
    x
end

push_new!(x::Vector{NewType}, y::NewUnionAll, doit::Bool) = push_new!(x, y.body, doit)

function linearize(x::Type)
    _linearize(Type[], x)
end
function _linearize(list::Vector{Type}, x::Type)
    if x ∉ list
        for y in supertypes(x)
            _linearize(list, y)
        end
        push_new!(list, x)
    end
    list
end

export linearize

function push_new!(x::Vector{Type}, y::DataType)
    if y ∉ x && y != x
        push!(x, y)
    end
    x
end

push_new!(x::Vector{Type}, y::UnionAll) = push_new!(x, y.body)
push_new!(x::Vector{Type}, y::Type) = x

supertypes(x::DataType) = x == Any ? Type[] : [supertype(x)]
supertypes(x::UnionAll) = supertypes(x.body)


text:
"""
**Rewiring graph**

(Based on Staniczenko et al., 2010)
This function calculate the rewiring graph R from the adjacency matrix A.
The rewiring graph indicates overlap species:
R[i,j] = 1 if i and j share at least one predator AND i has at least one predator
that does not prey on j (i is in that case an overlap species).
Used internally by rewire().
"""
function rewiring_graph(parameters)
  S = size(parameters[:A], 1)
  R = convert(Array{Int64,2}, zeros(S, S)) #pre-allocate
  Atmp = deepcopy(parameters[:A])
  links = findall(!iszero, Atmp)
  consumers = [links[x][1] for x = 1:length(links)]
  unique_consumers = unique(consumers)
  resources = [links[x][2] for x = 1:length(links)]
  unique_resources = unique(resources)
  for r1 in unique_resources #for each unique resource
    c1 = consumers[resources .== r1]
    if length(c1) > 1
      for r2 in unique_resources
        c2 = consumers[resources .== r2]
        shared_cons = [x ∈ c2 for x = c1]
        (sum(shared_cons) > 0) & (sum(.!shared_cons) > 0) ? R[r1,r2] = 1 : nothing
      end
    end
  end
  return R
end

"""
**Potential predators**

(Based on Staniczenko et al., 2010)
This function identify the potential predators from the rewiring graph R.
Used internally by rewire().
"""
function potential_newlinks(A, R, parameters)
  S = size(A, 1)
  R_slices = [R[:, i] for i in 1:S]
  #identify "overlap species" = species that share at least one
  #consumer but also have 1 unshared consumer (results in consumers competition)
  overlap_species = map(x -> findall(x .== 1), R_slices)
  sp = convert(BitArray{1}, [overlap_species[x] != Int64[] for x in 1:size(overlap_species,1)])
  #identify overlap species consumers
  potential_predators = map(x -> A[:,x], overlap_species)[sp]
  is_link = map(x -> unique(findall(!iszero, x)), potential_predators)
  id_predators = map(x -> (i->i[1]).(x), is_link)
  #identify species consumers
  predators = map(x -> unique(findall(!iszero, x)), [A[:, i] for i in 1:S])[sp]
  #new consumers are those who consume the overlap species but not the species under focus
  new_pred = map((x,y) -> filter(x -> x ∉ y, x), id_predators, predators)
  #pre-allocate empty array
  P = convert(Array{Int64,2}, zeros(S, S))
  #return the identified new potential consumers
  map((x,y) -> P[x,y] .= 1, new_pred, findall(sp))
  TR = parameters[:trophic_rank]
  idx = findall(!iszero, P)
  idx_c = (i->i[1]).(idx)
  idx_r = (i->i[2]).(idx)
  map((x,y) -> if TR[x] <= TR[y]
                 P[x,y]=0
               end,
       idx_c, idx_r
  )
  return P
end

"""
**Rewire**

(Based on Staniczenko et al., 2010)
This function identify the resources that are released following an extinction event
and sample a new predator for each of them (when possible) from the matrix of
potential predators (from `potential_predators()`).
"""
function Staniczenko_rewire(parameters)
  S = size(parameters[:A], 1)
  A = parameters[:A]
  #identify extinct species
  Ɛ = parameters[:extinctions]
  #identify potential new links
  R = rewiring_graph(parameters)
  P = potential_newlinks(A, R, parameters)
  #keep those that contain one of the released preys
  released_links = findall(!iszero, A[Ɛ,:])
  released_preys = unique((i->i[2]).(released_links))
  [filter!(x -> x .!= i, released_preys) for i in Ɛ]
  all_possible_predators = map(x -> findall(!iszero, x), [P[:,i] for i in released_preys])
  trm = findall(x -> x != Int64[], all_possible_predators)
  filter!(x -> x != Int64[], all_possible_predators)
  released_preys = released_preys[trm]
  #draw a new predator for each released preys
  new_predator = map(sample, all_possible_predators)
  A_rewired = deepcopy(A)
  A_rewired[Ɛ, :] .= 0
  A_rewired[:, Ɛ] .= 0
  map((x,y) -> A_rewired[x,y] = 1, new_predator, released_preys)
  return A_rewired
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
function getLinkLengthDict(adj_mtx)
    link_length_dict = Dict()
    for i=1:size(adj_mtx,1)
        for j=1:size(adj_mtx,2)
            if adj_mtx[i,j]==1
                link_length_dict[i,j] = 1
            end
        end
    end

    return link_length_dict
end

function getLinkLengthDict(start_node, end_node, link_length)
    link_length_dict = Dict()
    for k=1:length(start_node)
        link_length_dict[ start_node[k] , end_node[k] ] = link_length[k]
    end

    return link_length_dict
end

function getPathLength(path, link_length_dict)

    if length(link_length_dict) == 0
        return 0.0
    else
        path_length = 0
        for k=1:length(path)-1
            i = path[k]
            j = path[k+1]

            path_length += link_length_dict[i,j]
        end

        return path_length
    end

end

function getNoNode(start_node, end_node)
    no_node = max( maximum(start_node), maximum(end_node) )
    return no_node
end

function getAdjacency(start_node, end_node)

    no_node = max( maximum(start_node), maximum(end_node) )

    adj = zeros(Int, no_node, no_node)

    for k=1:length(start_node)
        i = start_node[k]
        j = end_node[k]
        adj[i,j] = 1

        # if option==:undirected
        #     adj[j,i] = 1
        # end
    end

    return adj
end


text:
using Pkg
Pkg.add(["JuMP", "AmplNLWriter", "Bonmin_jll"])
using JuMP, AmplNLWriter, Bonmin_jll

model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe))
@variable(model, 1 <= n <= 1000, Int)
@variable(model, 1 <= sx <= 76, Int)
@variable(model, 1 <= sy <= 162, Int)
@NLconstraint(model, c2, 56 <= (sx+1)*sx/2 <= 76)
@NLconstraint(model, c3, -162 <= (2*sy-n+1)*n/2 <= -134)
@NLconstraint(model, c4, n <= 1000)
@NLobjective(model, Max, sy*(sy + 1)/2)

optimize!(model)

print(solution_summary(model, verbose=true))


text:
using ArgParse
include("./test.jl")

"""
    parse_commandline()

Parse arguments from command line.

"""
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--graph-file"
            help = "graph file path"
            required = true
        "--pop-file"
            help = "population file path"
            required = true
        "-p"
            help = "plot results"
            action = :store_true
        "--enum-file"
            help = "all possible enumeration of original graph"
        "--benchmark"
            help = "type of benchmark to calculate. Default type is cut edges"
            default = "cut_edges"
        "--metric-kl"
            help = "Compute the Kullback-Leibler divergence between two distributions"
            default = 0
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    graph = load_graph_from_edge_list(parsed_args["graph-file"], parsed_args["pop-file"])
    test_enumeration = calculate_benchmark(graph, parsed_args["benchmark"])
    spanning_tree = spanning_tree_distribution(graph, parsed_args["enum-file"])
    if parsed_args["metric-kl"] != 0
        kl = calculate_Kullback_Leibler(spanning_tree, test_enumeration)
        if kl > parse(Int64, parsed_args["metric-kl"])
            println(kl)
            error("Maximum KL distance exceeded")
        end
    end
    if parsed_args["p"]
        fig = figure(figsize = (10, 5))
        plot_distribution(test_enumeration)
        plot_distribution(spanning_tree, "blue")
        show()
    end
end

main()

max_stars_repo_name: mggg/benchmarks

max_issues_repo_name: mggg/benchmarks

max_forks_repo_name: mggg/benchmarks


text:
# Deprecated functionality

macro profile(ex)
    Base.depwarn("`CUDAnative.@profile` is deprecated, use `CUDAdrv.@profile` instead", :profile)
    quote
        CUDAdrv.@profile begin
            $(esc(ex))
        end
    end
end

@deprecate nearest_warpsize(dev::CuDevice, threads::Integer) nextwarp(dev, threads)


text:
using AlgebraicMultigrid
import IterativeSolvers: cg
import Metis

#--------------------------------------------------------------------#

function SolvePoisson(mesh, K, f, solver)
    # Solve for hybrid variable
    if solver == 0
        @time Th   = K\f
    elseif solver == 1
        println("Direct solve:")
        PC  = 0.5.*(K.+K')
        # @time permetis, iperm = Metis.permutation(PC)
        # PCc = cholesky(PC, perm=convert(Vector{Int64},permetis))
        PCc = cholesky(PC)
        Th  = zeros(mesh.nf)
        dTh = zeros(mesh.nf,1)
        r   = zeros(mesh.nf,1)
        r  .= f - K*Th
        # @time Th   = K\f
        for rit=1:5
            r    .= f - K*Th
            println("It. ", rit, " - Norm of residual: ", norm(r)/length(r))
            if norm(r)/length(r) < 1e-10
                break
            end
            dTh  .= PCc\r
            Th  .+= dTh[:]
        end
    elseif solver == 2
        println("AMG preconditionned CG solver:")
        # ml = ruge_stuben(K)
        @time ml = smoothed_aggregation(K) #pas mal
        @time p  = aspreconditioner(ml)
        @time Th = cg(K, f, Pl = p, reltol=1e-6)
    end
    return Th
    end


text:
const HC = HomotopyContinuation

ImplicitFunction(f::HC.ModelKit.Expression) = ImplicitFunction(HC.System([f]))
ImplicitFunction(f::HC.ModelKit.System) = ImplicitFunction(HC.CompiledSystem(f))
function ImplicitFunction(F::HC.AbstractSystem)
    ImplicitFunction{HC.nvariables(F)}((x,y) -> first(F(SVector(x,y))))
end

max_issues_repo_name: saschatimme/ImplicitPlots

max_forks_repo_name: saschatimme/ImplicitPlots


text:
using Base.Test

module X1
    using Reexport
    @reexport module Y1
        const Z1 = 1
        export Z1
    end
end

@test union!(Set(), names(X1)) == union!(Set(), [:X1, :Y1, :Z1])
@test X1.Z1 == 1

module Y2
    const Z2 = 2
    export Z2
end
module X2
    using Reexport
    @reexport using Y2
end
@test union!(Set(), names(X2)) == union!(Set(), [:X2, :Y2, :Z2])
@test X2.Z2 == 2

module X3
    using Reexport
    module Y3
        const Z3 = 3
        export Z3
    end
    module Y4
        const Z4 = 4
        export Z4
    end
    @reexport using .Y3, .Y4
end
@test union!(Set(), names(X3)) == union!(Set(), [:X3, :Y3, :Y4, :Z3, :Z4])
@test X3.Z3 == 3
@test X3.Z4 == 4

module X4
    using Reexport
    @reexport importall Y2
end
@test union!(Set(), names(X4)) == union!(Set(), [:X4, :Y2, :Z2])
@test X4.Z2 == 2

module X5
    using Reexport
    module Y3
        const Z3 = 3
        export Z3
    end
    module Y4
        const Z4 = 4
        export Z4
    end
    @reexport importall .Y3, .Y4
end
@test union!(Set(), names(X3)) == union!(Set(), [:X3, :Y3, :Y4, :Z3, :Z4])
@test X3.Z3 == 3
@test X3.Z4 == 4


text:
# Copyright (c) 2019 Yury Dvorkin, Miles Lubin, and contributors                #src
#                                                                               #src
# This work is licensed under a Creative Commons Attribution-ShareAlike 4.0     #src
# International License. URL: http://creativecommons.org/licenses/by-sa/4.0.    #src

# # Power Systems

# **Originally Contributed by**: Yury Dvorkin and Miles Lubin

# This tutorial demonstrates how to formulate basic power systems engineering
# models in JuMP using a 3 bus example.

# We will consider basic "economic dispatch" and "unit commitment" models
# without taking into account transmission constraints.

# ## Illustrative example

# In the following notes for the sake of simplicity, we are going to use a three
# bus example mirroring the interface between Western and Eastern Texas. This
# example is taken from R. Baldick, "[Wind and Energy Markets: A Case Study of Texas](http://dx.doi.org/10.1109/JSYST.2011.2162798),"
# IEEE Systems Journal, vol. 6, pp. 27-34, 2012.

# ![Power Systems](../../assets/power_systems.png)

# For this example, we set the following characteristics of generators,
# transmission lines, wind farms and demands:

# | Quantity         | Generator 1 | Generator 2 |
# |:-----------------|:-----------:|------------:|
# | $g_{min}$, MW    | 0           | 300         |
# | $g_{max}$, MW    | 1000        | 1000        |
# | $c^g$, \$/MWh    | 50          | 100         |
# | $c^{g0}$, \$/MWh | 1000        | 0           |

# | Quantity      | Line 1 | Line 2 |
# |:--------------|:------:|-------:|
# | $f^{max}$, MW | 100    | 1000   |
# | $x$, p.u.     | 0.001  | 0.001  |

# | Quantity        | Wind farm 1 | Wind farm 2 |
# |:----------------|:-----------:|------------:|
# | $w^{f}$, MW     | 150         | 50          |
# | $c^{w}$, \$/MWh | 50          | 50          |

# | Quantity | Bus 1 | Bus 2 | Bus 3 |
# |:---------|:-----:|:-----:|------:|
# | $d$, MW  | 0     | 0     | 15000 |

# ## Economic dispatch

# Economic dispatch (ED) is an optimization problem that minimizes the cost of
# supplying energy demand subject to operational constraints on power system
# assets. In its simplest modification, ED is an LP problem solved for an
# aggregated load and wind forecast and for a single infinitesimal moment.

# Mathematically, the ED problem can be written as follows:

# ```math
# \min \sum_{i \in I} c^g_{i} \cdot g_{i} + c^w \cdot w,
# ```

# where $c_{i}$ and $g_{i}$ are the incremental cost (\$/MWh) and power output
# (MW) of the $i^{th}$ generator, respectively, and $c^w$ and $w$ are the
# incremental cost (\$/MWh) and wind power injection (MW), respectively.

# Subject to the constraints:

# * Minimum ($g^{\min}$) and maximum ($g^{\max}$) limits on power outputs of
#   generators:
#   $g^{\min}_{i} \leq g_{i} \leq g^{\max}_{i}.$
# * Constraint on the wind power injection:
#   $0 \leq w \leq w^f,$
#   where $w$ and $w^f$ are the wind power injection and wind power forecast,
#   respectively.
# * Power balance constraint:
#   $\sum_{i \in I} g_{i} + w = d^f,$
#   where $d^f$ is the demand forecast.

# Further reading on ED models can be found in A. J. Wood, B. F. Wollenberg, and
# G. B. Sheblé, "Power Generation, Operation and Control", Wiley, 2013.

# ## JuMP Implementation of Economic Dispatch

using JuMP
import DataFrames
import GLPK

# Define some input data about the test system.

## Maximum power output of generators
g_max = [1000, 1000]
## Minimum power output of generators
g_min = [0, 300]
## Incremental cost of generators
c_g = [50, 100]
## Fixed cost of generators
c_g0 = [1000, 0]
## Incremental cost of wind generators
c_w = 50
## Total demand
d = 1500
## Wind forecast
w_f = 200
nothing #hide

# Create a function solve_ed, which solves the economic dispatch problem for a
# given set of input parameters.
function solve_ed(g_max, g_min, c_g, c_w, d, w_f)
    ## Define the economic dispatch (ED) model
    ed = Model(GLPK.Optimizer)
    ## Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) ## power output of generators
    @variable(ed, 0 <= w <= w_f) ## wind power injection
    ## Define the objective function
    @objective(ed, Min, c_g' * g + c_w * w)
    ## Define the constraint on the maximum and minimum power output of each
    ## generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) ## maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) ## minimum
    ## Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)
    ## Define the power balance constraint
    @constraint(ed, sum(g) + w == d)
    ## Solve statement
    optimize!(ed)
    ## return the optimal value of the objective function and its minimizers
    return value.(g), value(w), w_f - value(w), objective_value(ed)
end

# Solve the economic dispatch problem

(g_opt, w_opt, ws_opt, obj) = solve_ed(g_max, g_min, c_g, c_w, d, w_f);

println("Dispatch of Generators: ", g_opt, " MW")
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

# ### Economic dispatch with adjustable incremental costs

# In the following exercise we adjust the incremental cost of generator G1 and
# observe its impact on the total cost.

c_g_scale_df = DataFrames.DataFrame(
    Symbol("Dispatch of Generator 1(MW)") => Float64[],
    Symbol("Dispatch of Generator 2(MW)") => Float64[],
    Symbol("Dispatch of Wind(MW)") => Float64[],
    Symbol("Spillage of Wind(MW)") => Float64[],
    Symbol("Total cost(\$)") => Float64[],
)
for c_g1_scale = 0.5:0.1:3.0
    ## Update the incremental cost of the first generator at every iteration.
    c_g_scale = [c_g[1] * c_g1_scale, c_g[2]]
    ## Solve the ed problem with the updated incremental cost
    g_opt, w_opt, ws_opt, obj = solve_ed(g_max, g_min, c_g_scale, c_w, d, w_f)
    push!(c_g_scale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

#-

ENV["COLUMNS"]=250 #hide

c_g_scale_df

# ## Modifying the JuMP model in place

# Note that in the previous exercise we entirely rebuilt the optimization model
# at every iteration of the internal loop, which incurs an additional
# computational burden. This burden can be alleviated if instead of re-building
# the entire model, we modify a specific constraint(s) or the objective
# function, as it shown in the example below.

# Compare the computing time in case of the above and below models.

function solve_ed_inplace(c_w_scale)
    start = time()
    obj_out = Float64[]
    w_out = Float64[]
    g1_out = Float64[]
    g2_out = Float64[]
    ed = Model(GLPK.Optimizer)
    ## Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) ## power output of generators
    @variable(ed, 0 <= w <= w_f ) ## wind power injection
    ## Define the objective function
    @objective(ed, Min, c_g' * g + c_w * w)
    ## Define the constraint on the maximum and minimum power output of each generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) ## maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) ## minimum
    ## Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)
    ## Define the power balance constraint
    @constraint(ed, sum(g) + w == d)
    optimize!(ed)
    for c_g1_scale = 0.5:0.01:3.0
        @objective(
            ed,
            Min,
            c_g1_scale * c_g[1] * g[1] + c_g[2] * g[2] + c_w_scale * c_w * w,
        )
        optimize!(ed)
        push!(obj_out, objective_value(ed))
        push!(w_out, value(w))
        push!(g1_out, value(g[1]))
        push!(g2_out, value(g[2]))
    end
    elapsed = time() - start
    print(string("elapsed time: ", elapsed, " seconds"))
    return obj_out, w_out, g1_out, g2_out
end

solve_ed_inplace(2.0)

# Adjusting specific constraints and/or the objective function is faster than
# re-building the entire model.

# ## A few practical limitations of the economic dispatch model

# ### Inefficient usage of wind generators

# The economic dispatch problem does not perform commitment decisions and, thus,
# assumes that all generators must be dispatched at least at their minimum power
# output limit. This approach is not cost efficient and may lead to absurd
# decisions. For example, if $ d = \sum_{i \in I} g^{\min}_{i}$, the wind power
# injection must be zero, i.e. all available wind generation is spilled, to meet
# the minimum power output constraints on generators.

# In the following example, we adjust the total demand and observed how it
# affects wind spillage.

demandscale_df = DataFrames.DataFrame(
    Symbol("Dispatch of Generators(MW)") => Float64[],
    Symbol("Dispatch of Generator 2(MW)") => Float64[],
    Symbol("Dispatch of Wind(MW)") => Float64[],
    Symbol("Spillage of Wind(MW)") => Float64[],
    Symbol("Total cost(\$)") => Float64[],
)

for demandscale = 0.2:0.1:1.5
    g_opt,w_opt,ws_opt,obj = solve_ed(g_max, g_min, c_g, c_w, demandscale*d, w_f)
    push!(demandscale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

demandscale_df

# This particular drawback can be overcome by introducing binary decisions on
# the "on/off" status of generators. This model is called unit commitment and
# considered later in these notes.

# For further reading on the interplay between wind generation and the minimum
# power output constraints of generators, we refer interested readers to R.
# Baldick, "Wind and Energy Markets: A Case Study of Texas," IEEE Systems
# Journal, vol. 6, pp. 27-34, 2012.

# ### Transmission-infeasible solution

# The ED solution is entirely market-based and disrespects limitations of the
# transmission network. Indeed, the flows in transmission lines would attain the
# following values:

# ```math
# f_{1-2} = 150 MW \leq f_{1-2}^{\max} = 100 MW
# ```

# ```math
# f_{2-3} = 1200 MW \leq f_{2-3}^{\max} = 1000 MW
# ```

# Thus, if this ED solution was enforced in practice, the power flow limits on
# both lines would be violated. Therefore, in the following section we consider
# the optimal power flow model, which amends the ED model with network
# constraints.

# The importance of the transmission-aware decisions is emphasized in E.
# Lannoye, D. Flynn, and M. O'Malley, "Transmission, Variable Generation, and
# Power System Flexibility," IEEE Transactions on Power Systems, vol. 30, pp.
# 57-66, 2015.

# ## Unit Commitment model

# The Unit Commitment (UC) model can be obtained from ED model by introducing
# binary variable associated with each generator. This binary variable can
# attain two values: if it is "1", the generator is synchronized and, thus, can
# be dispatched, otherwise, i.e. if the binary variable is "0", that generator
# is not synchronized and its power output is set to 0.

# To obtain the mathematical formulation of the UC model, we will modify the
# constraints of the ED model as follows:

# ```math
# g^{\min}_{i} \cdot u_{t,i} \leq g_{i} \leq g^{\max}_{i} \cdot u_{t,i},
# ```
# where $u_{i} \in \{0,1\}.$ In this constraint, if $u_{i} = 0$, then
# $g_{i}  = 0$. On the other hand, if $u_{i} = 1$, then
# $g^{max}_{i} \leq g_{i} \leq g^{min}_{i}$.

# For further reading on the UC problem we refer interested readers to G.
# Morales-Espana, J. M. Latorre, and A. Ramos, "Tight and Compact MILP
# Formulation for the Thermal Unit Commitment Problem," IEEE Transactions on
# Power Systems, vol. 28, pp. 4897-4908, 2013.

# In the following example we convert the ED model explained above to the UC
# model.

function solve_uc(g_max, g_min, c_g, c_w, d, w_f)
    ## Define the unit commitment (UC) model
    uc = Model(GLPK.Optimizer)
    ## Define decision variables
    @variable(uc, 0 <= g[i=1:2] <= g_max[i]) ## power output of generators
    @variable(uc, u[i = 1:2], Bin) ## Binary status of generators
    @variable(uc, 0 <= w <= w_f) ## wind power injection
    ## Define the objective function
    @objective(uc, Min, c_g' * g + c_w * w)
    ## Define the constraint on the maximum and minimum power output of each
    ## generator.
    @constraint(uc, [i = 1:2], g[i] <= g_max[i] * u[i]) ## maximum
    @constraint(uc, [i = 1:2], g[i] >= g_min[i] * u[i]) ## minimum
    ## Define the constraint on the wind power injection
    @constraint(uc, w <= w_f)
    ## Define the power balance constraint
    @constraint(uc, sum(g) + w == d)
    ## Solve statement
    optimize!(uc)
    status = termination_status(uc)
    if status != MOI.OPTIMAL
        return status, zeros(length(g)), 0.0, 0.0, zeros(length(u)), Inf
    end
    return status, value.(g), value(w), w_f - value(w), value.(u), objective_value(uc)
end

# Solve the economic dispatch problem
status, g_opt, w_opt, ws_opt, u_opt, obj =
    solve_uc(g_max, g_min, c_g, c_w, d, w_f)

println("Dispatch of Generators: ", g_opt[:], " MW")
println("Commitments of Generators: ", u_opt[:])
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

# ### Unit Commitment as a function of demand

# After implementing the UC model, we can now assess the interplay between the
# minimum power output constraints on generators and wind generation.

uc_df = DataFrames.DataFrame(
    Symbol("Commitment of Generator 1(MW)") => Float64[],
    Symbol("Commitment of Generator 2(MW)") => Float64[],
    Symbol("Dispatch of Generator 1(MW)") => Float64[],
    Symbol("Dispatch of Generator 2(MW)") => Float64[],
    Symbol("Dispatch of Wind(MW)") => Float64[],
    Symbol("Spillage of Wind(MW)") => Float64[],
    Symbol("Total cost(\$)") => Float64[],
)

for demandscale = 0.2:0.1:1.5
    status, g_opt, w_opt, ws_opt, u_opt, obj =
        solve_uc(g_max, g_min, c_g, c_w, demandscale*d, w_f)
    if status == MOI.OPTIMAL
        push!(uc_df, (u_opt[1], u_opt[2], g_opt[1], g_opt[2], w_opt, ws_opt, obj))
    else
        println("Status: $status for demandscale = $demandscale \n")
    end
end

#-

uc_df

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# A function that replaces a value at a specific place (say the third )
function Base.replace(x::String, sep::String, rep::String, pos::Integer)
  a = split(x, sep)
  if (length(a)<pos+1) return x end
  join([a[1:(pos-1)]..., join(a[pos:(pos+1)], rep), a[(pos+2):end]...], sep)
end

replace("a1 is a2 is a3", " ", ".", 2)
replace("a1 is a2 is a3", "-", ".", 2)
replace("a1 is a2 is a3", "2", ".", 2)
replace("a1 is a2 is a3", "2", ".", 1)

cd("C:/Users/francis.smart.ctr/GitDir/SoftDatesJL/SoftDatesJL")

using Pkg
Pkg.activate(".")
Pkg.add("Distributions")
using Pkg, Dates, DataFrames, Distributions, CSV

# used lipsum.com to grab 2000 year old document (https://www.lipsum.com/)
lipsum = read("sampledata/Lipsum.txt", String)

lip_par = [replace(i, r"\r|\n" => "") for i in split(lipsum, "\n\r") if i!=""]
lip_sen  = [replace(i, r"\r|\n|^ " => "") for i in split(lipsum, ".") if i!=""]
lip_wrd  = [replace(i, r"\r|\n|^ " => "") for i in split(lipsum, r"\n\r|\s+") if i!=""]

# Units from wikipedia (https://en.wikipedia.org/wiki/United_States_customary_units)
units = split(read("sampledata/Units.txt", String), "\r\n")

sampleset = [sample(lip_sen, 200, replace = false)..., units...]

# DataGenerators
# Type1
# Number of dates to generate
n = 250

dtrange = (Date(now()) - Year(5)):Day(1):Date(now())

startdt = firstdayofweek.(rand(dtrange, n)) + Day.(rand([-2,-1,0,0,0,0,1], n))
enddt  = startdt .+ Day.(rand([4,5,6,6,7,7,7,7,7,8], n))

# Number of days covered from the start to the end
span = Dates.value.(enddt-startdt)

# The reporting rate by date -> I set at .9
coverage = [rand(Binomial(span[i], .9),1)[1] for i in 1:n]

daydates = [sort(sample(startdt[i]:Day(1):enddt[i], coverage[i], replace = false)) for i in 1:n]

notes =  [join(sample(lip_sen, rand(Binomial(3, .25),1)[1]), " ") for i in 1:n]

txtfield = [[join(sample(sampleset, rand([1,2,2,3,4,5,6])),
             sample([". ", "; ", " ", " - ", ", ", ".\n\r", "\n\r", "\n\r\n\r"], 1)[1])
             for j in 1:length(daydates[i])]  for i in 1:n]

sep = sample(split(raw"./- ","") , n)

# txtblocks = [join((notes[i], , ), )]

# We have start dates and end dates and where things get tricky is that we
# have individual records that start at Start and end at End. They are listed
# at the beginning of lines.

# Some lines start with a note before the first date is entered.
# Not every date is covered

# Information to be spread accross the dates
V0 = collect(1:n)
V1 = [join(sample(lip_wrd, rand(Binomial(9, .25),1)[1]), " ") for i in 1:n]
V2 = [sample(lip_sen, 1)[1] for i in 1:n]
V3 = [sample(1:(n*5), 1)[1] for i in 1:n]

# Generate the text fields
[ begin
   x = [join([Dates.format.(daydates[i][j],"dd$(sep[i])mm$(sep[i])yyyy"), txtfield[i][j]]," ") for j in 1:length(daydates[i])]
   replace(join([notes[i], x...], "\n\r"), r"^\n\r" => "")
   end
   for i = 1 : n ]

# Generate date type 1
fmlist = [
  ["dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e mm$(sep[i])yyyy" for i in 1:n] ,
  ["e mm$(sep[i])yy"   for i in 1:n] ,
  ["yy$(sep[i])mm$(sep[i])dd"   for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yyyy" for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yy"   for i in 1:n] ,
  ["U dd, yyyy" for i in 1:n] ,
  ["U dd, yy"   for i in 1:n] ,
  ["u dd, yyyy" for i in 1:n] ,
  ["u dd, yy"   for i in 1:n] ]

txtblocks =
    [begin
     [ begin
     x = [join([Dates.format.(daydates[i][j], fmlist[ii][i]), txtfield[i][j]]," ")
       for j = 1:length(daydates[i])]
     replace(join([notes[i], x...], "\n\r"), r"^\n\r" => "")
     end
     for i = 1 : n ]
    end
    for ii = 1:length(fmlist)]

dt = DataFrame(
  type = 1,
  Start = startdt ,
  End = enddt,
  Dates = daydates,
  Text = txtblocks[1],
  V0=V0, V1=V1, V2=V2, V3=V3)

for i in 2:length(fmlist)
    append!(dt,
     DataFrame(type = i,
      Start = startdt ,
      End = enddt,
      Dates = daydates,
      Text = txtblocks[i],
      V0=V0, V1=V1, V2=V2, V3=V3))
    end

CSV.write("sampledata/SampleDates.csv", dt)

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

# Save starting text
dt[:Text0] = dt[:Text]

# Input Errors

    sloppy_rt = .10

    sloppy_entries = sample(1:n, Int(round(n*sloppy_rt)), replace = false) |> sort

    dt[:sloppy] = [dt[i, :V0] ∈ sloppy_entries for i in 1:size(dt)[1]]

    symbolset = split(raw"./- ","")

    txsloppy = dtrange[dt[:sloppy], :Text]

    # Replace some symbols in the sloppy dates with alternative symbols
    for i in 1:size(txsloppy)[1]
      position = sample(1:2,1)[1]
      for j in 1:length(symbolset)
        txsloppy[i] = replace(txsloppy[i],
          string(symbolset[j % length(symbolset) + 1]),
          string(symbolset[j]),
          position)
      end
    end

    dt[dt[:sloppy], :Text] = txsloppy

# There are also some entries which have the wrong date - 15-19 -> +- 1
    error_rt = .10

    error_entries = sample(1:n, Int(round(n*error_rt)), replace = false) |> sort

    dt[:error] = [dt[i, :V0] ∈ error_entries for i in 1:size(dt)[1]]

    txerror = dt[dt[:error], :Text]

    txerror = [replace(txerror[i], r"1[6-8]\b" => (x -> "$(parse(Int, x)-1)")) for i in 1:size(txerror)[1]]

    # 10% of years are off
    dt[dt[:error], :Text] = txerror

# Drop zeros for some portion of dates
    drop0_rt = .10

    drop0_entries = sample(1:n, Int(round(n*drop0_rt)), replace = false) |> sort

    dt[:drop0] = [dt[i, :V0] ∈ drop0_entries for i in 1:size(dt)[1]]

    txdrop0 = dt[dt[:drop0], :Text]

    txdrop0 = [replace(txdrop0[i], "0"=>"",count = sample([1,1,1,1,2,3],1)[1]) for i in 1:size(txdrop0)[1]]

    dt[dt[:drop0], :Text] = txdrop0

# Rand Δ of numbers are replaced with other random numbers 10 -> 30
    rand_rt = .10

    rand_entries = sample(1:n, Int(round(n*rand_rt)), replace = false) |> sort

    dt[:rand] = [dt[i, :V0] ∈ rand_entries for i in 1:size(dt)[1]]

    txrand = dt[dt[:rand], :Text]

    txrand = [replace(txrand[i], string(sample(0:9,1)[1]), string(sample(0:9,1)[1]), 1) for i in 1:size(txrand)[1]]

    dt[dt[:rand], :Text] = txrand

# Non-standard abbreviation of days of the week: Mon->M or Mo, Tuesday -> Tuesda
    abbrev_rt = .30

    abbrev_entries = sample(1:n, Int(round(n*abbrev_rt)), replace = false) |> sort

    dt[:abbrev] = [dt[i, :V0] ∈ abbrev_entries for i in 1:size(dt)[1]]

    txabbrev = dt[dt[:abbrev], :Text]

    daynames = join(dayname.(now() .+ Day.(0:6)), "|")
    dayabbrs = join(dayabbr.(now() .+ Day.(0:6)), "|")

    randcut(x) = x[1:sample(1:length(x)-1)[1]]

    txabbrev = [replace(txabbrev[i], Regex("\\b($daynames|$dayabbrs)\\b") => randcut ) for i in 1:size(txabbrev)[1]]

    dt[dt[:abbrev], :Text] = txabbrev

CSV.write("sampledata/SampleDatesErrors.csv", dt)

########################################################################
dayranges =
  [begin
    y = daydates[i]
    x = sort(sample(1:length(y), minimum([2,length(y)]), replace = false))
    (length(x) == 2) ? ([y[1:(x[1]-1)], [y[x[1]],y[x[2]]], y[(x[2]+1):end]]) : [y, y[1:0], y[1:0]]
   end
   for i in 1:n]

seprange = sample( ["-", "thru", "through", "till", "to"] , n)

# Generate date type 1
fm2list = [
  ["dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e mm$(sep[i])yyyy" for i in 1:n] ,
  ["e mm$(sep[i])yy"   for i in 1:n] ,
  ["yyyy$(sep[i])mm$(sep[i])dd" for i in 1:n] ,
  ["yy$(sep[i])mm$(sep[i])dd"   for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yyyy" for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yy"   for i in 1:n] ,
  ["U dd, yyyy" for i in 1:n] ,
  ["U dd, yy"   for i in 1:n] ]

fr1list = [
    ["dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
    ["dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
    ["E dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
    ["E dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
    ["e dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
    ["e dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
    ["e mm$(sep[i])yyyy" for i in 1:n] ,
    ["e mm$(sep[i])yy"   for i in 1:n] ,
    ["yyyy$(sep[i])mm$(sep[i])dd" for i in 1:n] ,
    ["yy$(sep[i])mm$(sep[i])dd"   for i in 1:n] ,
    ["mm$(sep[i])dd$(sep[i])yyyy" for i in 1:n] ,
    ["mm$(sep[i])dd$(sep[i])yy"   for i in 1:n] ,
    ["mm$(sep[i])dd" for i in 1:n] ,
    ["mm$(sep[i])dd" for i in 1:n]]

fr2list = [
  ["dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["E dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yyyy" for i in 1:n] ,
  ["e dd$(sep[i])mm$(sep[i])yy"   for i in 1:n] ,
  ["e mm$(sep[i])yyyy" for i in 1:n] ,
  ["e mm$(sep[i])yy"   for i in 1:n] ,
  ["yyyy$(sep[i])mm$(sep[i])dd" for i in 1:n] ,
  ["yy$(sep[i])mm$(sep[i])dd"   for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yyyy" for i in 1:n] ,
  ["mm$(sep[i])dd$(sep[i])yy"   for i in 1:n] ,
  ["dd$(sep[i])yyyy" for i in 1:n] ,
  ["dd$(sep[i])yyyy" for i in 1:n] ]

dtrange = dt[1:0,]

for k in 1:length(fm2list)

  textblock = [ begin
     dtr = dayranges[i]
      x1 =
        if length(dtr[1]) > 0
        join([join([Dates.format(dtr[1][j],fm2list[k][i]), txtfield[i][j]]," ")
              for j in 1:length(dtr[1])], "\n\r")
        else ""
        end

      x2 = if length(dtr[2]) == 2
        d1 = Dates.format(dtr[2][1], fr1list[k][i])
        sepr = seprange[i]
        d2 = Dates.format(dtr[2][2], fr2list[k][i])
        t1 = txtfield[i][length(dtr[1])+1]
        "$d1 $sepr $d2 $t1"
        else ""
        end

      x3 =
       if length(dtr[3]) > 0
         join([join([Dates.format(dtr[3][j],fm2list[k][i]), txtfield[i][j]]," ")
           for j in 1:length(dtr[3])], "\n\r")
       else ""
       end

      replace(join([x1, x2, x3], "\n\r"), r"^\n\r"=>"")
     end
     for i = 1 : n ]

    df_temp = DataFrame(
        type = k,
        Start = startdt ,
        End = enddt,
        Dates = daydates,
        Text = textblock,
        V0=V0, V1=V1, V2=V2, V3=V3)

    append!(dtrange, df_temp)
end

dtrange
dtrange[:Text]

CSV.write("sampledata/SampleDateRanges.csv", dtrange)

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

# Save starting text
dtrange[:Text0] = dtrange[:Text]

# Input Errors

    sloppy_rt = .10

    sloppy_entries = sample(1:n, Int(round(n*sloppy_rt)), replace = false) |> sort

    dtrange[:sloppy] = [dtrange[i, :V0] ∈ sloppy_entries for i in 1:size(dtrange)[1]]

    symbolset = split(raw"./- ","")

    txsloppy = dtrange[dtrange[:sloppy], :Text]

    # Replace some symbols in the sloppy dates with alternative symbols
    for i in 1:size(txsloppy)[1]
      position = sample(1:2,1)[1]
      for j in 1:length(symbolset)
        txsloppy[i] = replace(txsloppy[i],
          string(symbolset[j % length(symbolset) + 1]),
          string(symbolset[j]),
          position)
      end
    end

    dtrange[dtrange[:sloppy], :Text] = txsloppy

# There are also some entries which have the wrong date - 15-19 -> +- 1
    error_rt = .10

    error_entries = sample(1:n, Int(round(n*error_rt)), replace = false) |> sort

    dtrange[:error] = [dtrange[i, :V0] ∈ error_entries for i in 1:size(dtrange)[1]]

    txerror = dtrange[dtrange[:error], :Text]

    txerror = [replace(txerror[i], r"1[6-8]\b" => (x -> "$(parse(Int, x)-1)")) for i in 1:size(txerror)[1]]

    # 10% of years are off
    dtrange[dtrange[:error], :Text] = txerror

# Drop zeros for some portion of dates
    drop0_rt = .10

    drop0_entries = sample(1:n, Int(round(n*drop0_rt)), replace = false) |> sort

    dtrange[:drop0] = [dtrange[i, :V0] ∈ drop0_entries for i in 1:size(dtrange)[1]]

    txdrop0 = dtrange[dtrange[:drop0], :Text]

    txdrop0 = [replace(txdrop0[i], "0"=>"",count = sample([1,1,1,1,2,3],1)[1]) for i in 1:size(txdrop0)[1]]

    dtrange[dtrange[:drop0], :Text] = txdrop0

# Rand Δ of numbers are replaced with other random numbers 10 -> 30
    rand_rt = .10

    rand_entries = sample(1:n, Int(round(n*rand_rt)), replace = false) |> sort

    dtrange[:rand] = [dtrange[i, :V0] ∈ rand_entries for i in 1:size(dtrange)[1]]

    txrand = dtrange[dtrange[:rand], :Text]

    txrand = [replace(txrand[i], string(sample(0:9,1)[1]), string(sample(0:9,1)[1]), 1) for i in 1:size(txrand)[1]]

    dtrange[dtrange[:rand], :Text] = txrand

# Non-standard abbreviation of days of the week: Mon->M or Mo, Tuesday -> Tuesda
    abbrev_rt = .30

    abbrev_entries = sample(1:n, Int(round(n*abbrev_rt)), replace = false) |> sort

    dtrange[:abbrev] = [dtrange[i, :V0] ∈ abbrev_entries for i in 1:size(dtrange)[1]]

    txabbrev = dtrange[dtrange[:abbrev], :Text]

    daynames = join(dayname.(now() .+ Day.(0:6)), "|")
    dayabbrs = join(dayabbr.(now() .+ Day.(0:6)), "|")

    randcut(x) = x[1:sample(1:length(x)-1)[1]]

    txabbrev = [replace(txabbrev[i], Regex("\\b($daynames|$dayabbrs)\\b") => randcut ) for i in 1:size(txabbrev)[1]]

    dtrange[dtrange[:abbrev], :Text] = txabbrev

CSV.write("sampledata/SampleDatesRangeErrors.csv", dtrange)


text:
using MacroTools
using MacroTools: postwalk

@testset "rewrite_struct" begin
    def = Meta.parse("struct Sequence <: BNFNode
                  elements::Tuple{Vararg{<:BNFNode}}
                  function Sequence(elements...)
                      new(elements)
                  end
              end")
    got = postwalk(rmlines, AnotherParser.rewrite_struct(def))
    want = postwalk(rmlines,
                    Meta.parse("Base.@__doc__ struct Sequence <: BNFNode
                                    source::Union{Nothing, LineNumberNode}
                                    elements::Tuple{Vararg{<:BNFNode}}
                                    function Sequence(elements...)
                                        new(nothing, elements)
                                    end
                                    function Sequence(source::Union{Nothing, LineNumberNode},
                                                      elements...)
                                        new(source, elements)
                                    end
                                end"))
    # ??? Why does Expr == Expr fail?
    # @test got == want
    @test string(got) == string(want)
end


text:
# ---
# title: 1358. Number of Substrings Containing All Three Characters
# id: problem1358
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/number-of-substrings-containing-all-three-characters/description/>
# hidden: true
# ---
#
# Given a string `s` consisting only of characters _a_ , _b_ and _c_.
#
# Return the number of substrings containing **at least**  one occurrence of all
# these characters _a_ , _b_ and _c_.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abcabc"
#     Output: 10
#     Explanation: The substrings containing at least one occurrence of the characters  _a_ ,  _b_  and  _c are "_abc _" , "_abca _" , "_abcab _" , "_abcabc _" , "_bca _" , "_bcab _" , "_bcabc _" , "_cab _" , "_cabc _" _and _ "_abc _" _( **again** ) _._
#
#
# **Example 2:**
#
#
#
#     Input: s = "aaacb"
#     Output: 3
#     Explanation: The substrings containing at least one occurrence of the characters  _a_ ,  _b_  and  _c are "_aaacb _" , "_aacb _" _and _ "_acb _"._ __
#
#
# **Example 3:**
#
#
#
#     Input: s = "abc"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `3 <= s.length <= 5 x 10^4`
#   * `s` only consists of  _a_ , _b_ or _c  _characters.
#
#
## @lc code=start
using LeetCode

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


text:
# additional math functions & functors

sqr(x::Number) = x * x
rcp(x::Number) = one(x) / x

rsqrt(x::Number) = one(x) / sqrt(x)
rcbrt(x::Real) = one(x) / cbrt(x)

xlogx(x::AbstractFloat) = x > zero(x) ? x * log(x) : zero(x)
xlogx(x::Real) = xlogx(float(x))

xlogy{T<:AbstractFloat}(x::T, y::T) = x > zero(T) ? x * log(y) : zero(x)
xlogy{T<:Real}(x::T, y::T) = xlogy(float(x), float(y))
xlogy(x::Real, y::Real) = xlogy(promote(x, y)...)

logistic(x::AbstractFloat) = rcp(one(x) + exp(-x))
logistic(x::Real) = logistic(float(x))

logit(x::AbstractFloat) = log(x / (one(x) - x))
logit(x::Real) = logit(float(x))

softplus(x::AbstractFloat) = x <= 0 ? log1p(exp(x)) : x + log1p(exp(-x))
softplus(x::Real) = softplus(float(x))

invsoftplus(x::AbstractFloat) = log(exp(x) - one(x))
invsoftplus(x::Real) = invsoftplus(float(x))

logsumexp{T<:AbstractFloat}(x::T, y::T) = x > y ? x + log1p(exp(y - x)) : y + log1p(exp(x - y))
logsumexp{T<:Real}(x::T, y::T) = logsumexp(float(x), float(y))
logsumexp(x::Real, y::Real) = logsumexp(promote(x, y)...)

@vectorize_1arg Number sqr
@vectorize_1arg Number rcp
@vectorize_1arg Real rsqrt
@vectorize_1arg Real rcbrt

@vectorize_1arg Real xlogx
@vectorize_2arg Real xlogy
@vectorize_1arg Real logistic
@vectorize_1arg Real logit
@vectorize_1arg Real softplus
@vectorize_1arg Real invsoftplus

const sigmoid = logistic


text:
using MacroTools: @forward
export AdjointRegister

"""
    AdjointRegister{B, T, RT} <: AbstractRegister{B, T}

Lazy adjoint for a quantum register.
"""
struct AdjointRegister{B, RT <: AbstractRegister{B}} <: AbstractRegister{B}
    parent::RT
end

Base.parent(reg::AdjointRegister) = reg.parent

"""
    adjoint(register) -> register

Lazy adjoint for quantum registers.
"""
Base.adjoint(reg::AbstractRegister) = AdjointRegister(reg)
Base.adjoint(reg::AdjointRegister) = parent(reg)

viewbatch(reg::AdjointRegister, i::Int) = adjoint(viewbatch(parent(reg), i))

@forward AdjointRegister.parent nqubits, nremain, nactive

function Base.summary(io::IO, reg::AdjointRegister{B, RT}) where {B, RT}
    print(io, "adjoint(", summary(reg.parent), ")")
end


text:
using NLsolve
using Plots

# Parameters
const f = 2.0 # depensation in piscivores (B)
const m = 0.01 # yr⁻¹ natural piscivore mortality
const a = 0.37 # B recruitment assymptote (shape parameter)
const r = 8.7 # mg⋅m⁻²⋅d⁻¹ maximum recycling of P
const q = 8 # Or 2, steepness coefficient of R
const k = 90 # mg⋅m⁻² half saturation constant of P recycling #h in matlab
const s = 0.001 # d⁻¹ flushing rate # matlab q1
const b = 0.5 # adult density when recruitment = a/2 (B recruitment shape parameter)
const α = 0.8 # mg⋅m⁻²⋅d⁻¹ maximum growth rate of A per P input rate
const v = 0.07 # m⋅d⁻¹ sinking rate per algal cell
const g = 0.1 # d⁻¹ C grazing coefficient

# Drivers
# H yr⁻¹ harvest rate
#H = [0.1, 0.4]
# DOC mg/m2 dissolved organic carbon
DOC = [6.0, 6.0, 18.0, 18.0]
# L mg⋅m⁻²⋅d⁻¹ P loading
Ploading = 0.0:0.1:6.0

# State variables
# B g/m2 piscivore population biomass (pike)
## W g/m2 planktivore population biomass (bream)
#W = @. 3.2*exp(-3.6926*B)
## C mm zooplankton community length
#C = W < 2.58 ? 0.5 : 0.35 # mm
C = [0.5, 0.35, 0.5, 0.35] # For steady state, we assume a low or high harvest rate.
## Z m depth of thermocline
z = @. 13.1826/(1 + DOC)^0.61 # log10(z) = 1.12 - 0.61log10(DOC+1)
## P mg/m2 phosphorus
## R mg/m2 recycled phosphorus
## A mg/m2 chlorophyll

function beisnernl!(F, x, L, z, C)
    R = (r*(x[1]^q))/(k^q + x[1]^q)
    F[1] = α*x[2]*x[1] - ((v/z) + g*C + s)*x[1]
    F[2] = L + R - α*x[2]*x[1] - s*x[2]
end

handle = Vector{Plots.Plot}(undef,4)

for i in 1:4
    stable = Vector{Tuple{Float64,Float64}}(undef,0)
    unstable = Vector{Tuple{Float64,Float64}}(undef,0)
    for L in Ploading
        # Lower bound
        sol = nlsolve((F,x) -> beisnernl!(F,x, L, z[i], C[i]), [1.0;0.0])
        if converged(sol)
            push!(stable, (L, sol.zero[1]))
        end
        # Upper bound
        sol = nlsolve((F,x) -> beisnernl!(F,x, L, z[i], C[i]), [300.0;0.0])
        if converged(sol)
            push!(stable, (L, sol.zero[1]))
        end
        # Instability region
        sol = nlsolve((F,x) -> beisnernl!(F,x, L, z[i], C[i]), [70.0;0.0])
        if converged(sol)
            push!(unstable, (L, sol.zero[1]))
        end
    end
    handle[i] = scatter(stable; legend = false, markersize = 3)
    xlims!((0,6))
    ylims!((0, 300))
    scatter!(handle[i], unstable; legend = false, color = :red, markershape=:hline)
end

plot(handle..., layout = (2,2))


text:
# Samuel J Palmer
# Based on "A New Method of Interpolation and Smooth Curve Fitting Based on Local Parameters", Akima, 1970

module CubicSplines

function extrapolate(xs, ys)
    # Quadratic extrapolation of (x4,y4), (x5,y5)
    x1, x2, x3 = xs
    y1, y2, y3 = ys
    x4 = x3-x1+x2
    x5 = x3-x1+x3
    y4 = y3 + 2*(x4-x3)*(y3-y2)/(x3-x2) - (x4-x3)*(y2-y1)/(x2-x1)
    y5 = y4 + 2*(x5-x4)*(y4-y3)/(x4-x3) - (x5-x4)*(y3-y2)/(x3-x2)
    return x4, x5, y4, y5
end

struct CubicSpline
    xs
    ys
    ts
    logx
    function CubicSpline(xs, ys; logx=false)
        @assert length(xs) == length(ys)
        @assert issorted(xs)
        if !issorted(xs)
            idx = sortperm(xs)
            xs = xs[idx]
            ys = ys[idx]
        end
        if logx
            xs = log.(xs)
        end
        xb, xa, yb, ya = extrapolate(xs[3:-1:1], ys[3:-1:1])
        xy, xz, yy, yz = extrapolate(xs[end-2:end], ys[end-2:end])
        xs = [xa; xb; xs; xy; xz]
        ys = [ya; yb; ys; yy; yz]
        ms = zeros(eltype(ys), length(xs)-1)
        for i in 1:length(ms)-1
            ms[i] = (ys[i+1]-ys[i]) ./ (xs[i+1]-xs[i])
        end
        ts = zeros(eltype(ys), length(xs))
        for i in 3:length(xs)-2
            m1, m2, m3, m4 = ms[i-2:i+1]
            if m1 == m2 && m3 == m4
                ts[i] = (m2+m3)/2
            else
                numer = abs(m4-m3)*m2 + abs(m2-m1)*m3
                denom = abs(m4-m3)    + abs(m2-m1)
                ts[i] = numer ./ denom
            end
        end
        return new(xs[3:end-2], ys[3:end-2], ts[3:end-2], logx)
    end
end

function Base.getindex(spline::CubicSpline, x::Real)
    if spline.logx
        x = log(x)
    end
    if x < spline.xs[1]
        throw("x too small ($x < $(spline.xs[1]))")
    elseif x > spline.xs[end]
        throw("x too big ($x > $(spline.xs[end]))")
    end
    i = 1
    while !(spline.xs[i] <= x <= spline.xs[i+1])
        i += 1
    end
    x1, x2 = spline.xs[i:i+1]
    y1, y2 = spline.ys[i:i+1]
    t1, t2 = spline.ts[i:i+1]
    p0 = y1
    p1 = t1
    p2 = (3(y2-y1)/(x2-x1)-2t1-t2)/(x2-x1)
    p3 = (t1+t2-2(y2-y1)/(x2-x1))/(x2-x1)^2
    return p0 + p1*(x-x1) + p2*(x-x1)^2 + p3*(x-x1)^3
end

function Base.getindex(spline::CubicSpline, xs::AbstractArray{<:Real,1})
    return [spline[x] for x in xs]
end

function Base.firstindex(spline::CubicSpline)
    return spline.xs[1]
end

function Base.lastindex(spline::CubicSpline)
    return spline.xs[end]
end

export CubicSpline

end # module


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

using Test
using Future
using SparseArrays

@testset "Future.copy! for AbstractSet" begin
    for S = (Set, BitSet)
        s = S([1, 2])
        for a = ([1], UInt[1], [3, 4, 5], UInt[3, 4, 5])
            @test s === Future.copy!(s, Set(a)) == S(a)
            @test s === Future.copy!(s, BitSet(a)) == S(a)
        end
    end
end

@testset "Future.copy! for AbstractDict" begin
    s = Dict(1=>2, 2=>3)
    for a = ([3=>4], [0x3=>0x4], [3=>4, 5=>6, 7=>8], Pair{UInt,UInt}[3=>4, 5=>6, 7=>8])
        @test s === Future.copy!(s, Dict(a)) == Dict(a)
        if length(a) == 1 # current limitation of Base.ImmutableDict
            @test s === Future.copy!(s, Base.ImmutableDict(a[])) == Dict(a[])
        end
    end
end

@testset "Future.copy! for AbstractVector" begin
        s = Vector([1, 2])
        for a = ([1], UInt[1], [3, 4, 5], UInt[3, 4, 5])
            @test s === Future.copy!(s, Vector(a)) == Vector(a)
            @test s === Future.copy!(s, SparseVector(a)) == Vector(a)
        end
end

@testset "Future.copy! for AbstractArray" begin
    @test_throws ArgumentError Future.copy!(zeros(2, 3), zeros(3, 2))
    s = zeros(2, 2)
    @test s === Future.copy!(s, fill(1, 2, 2)) == fill(1, 2, 2)
    @test s === Future.copy!(s, fill(1.0, 2, 2)) == fill(1.0, 2, 2)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using MechanismGeometries
using RigidBodyDynamics
const rbd = RigidBodyDynamics
using GeometryBasics
using StaticArrays: SVector, SDiagonal
using CoordinateTransformations: AffineMap, IdentityTransformation, LinearMap, Translation, transform_deriv
using ColorTypes: RGBA, BGR
using ValkyrieRobot
using Test
using LinearAlgebra

function homog(::IdentityTransformation)
    Matrix(1.0I, 4, 4)
end

function homog(m::AffineMap)
    H = Matrix(1.0I, 4, 4)
    H[1:3, 1:3] .= transform_deriv(m, SVector(0., 0., 0.))
    H[1:3, 4] .= m(SVector(0., 0., 0.))
    H
end

homog(m::LinearMap) = homog(m ∘ Translation(0, 0, 0))

homog(t::Translation) = homog(AffineMap(Matrix(1.0I, 3, 3), t(SVector(0., 0., 0.))))

@testset "MechanismGeometries" begin
    @testset "skeletons" begin
        @testset "unit sphere" begin
            frame = CartesianFrame3D("com")
            inertia = SpatialInertia(frame, SDiagonal(1., 1, 1), SVector(0., 0, 0), 1.0)
            # I = diag([1, 1, 1]), M = 1
            # so the equivalent ellipsoid must have the same moment of inertia
            # Ixx = 1/5 * M * (b^2 + c^2)
            # so 1 = 1/5 * (b^2 + c^2) and b = c
            # so 1 = 2/5 b^2
            # and b^2 = 5/2 -> b = sqrt(5/2)
            element = MechanismGeometries.inertial_ellipsoid(inertia)
            @test element.frame === frame
            @test homog(element.transform) ≈ homog(AffineMap(SDiagonal(sqrt(5/2), sqrt(5/2), sqrt(5/2)), SVector(0., 0, 0)))
            @test element.geometry isa HyperSphere
            @test radius(element.geometry) ≈ 1
        end

        @testset "heavy sphere" begin
            frame = CartesianFrame3D("com")
            inertia = SpatialInertia(frame, SDiagonal(1., 1, 1), SVector(0., 0, 0), 5.0)
            # I = diag([1, 1, 1]), M = 5
            # so the equivalent ellipsoid must have the same moment of inertia
            # Ixx = 1/5 * M * (b^2 + c^2)
            # so 1 = 1/5 * 5 * (b^2 + c^2) and b = c
            # so 1 = 2 b^2
            # and b^2 = 1/2 -> b = sqrt(1/2)
            element = MechanismGeometries.inertial_ellipsoid(inertia)
            @test element.frame === frame
            @test homog(element.transform) ≈ homog(LinearMap(SDiagonal(sqrt(1/2), sqrt(1/2), sqrt(1/2))))
            @test element.geometry isa HyperSphere
            @test radius(element.geometry) ≈ 1
        end

        @testset "translated sphere" begin
            frame = CartesianFrame3D("com")
            inertia = SpatialInertia(frame, SDiagonal(1., 1, 1), SVector(1., 0, 0), 5.0)
            @test center_of_mass(inertia).v ≈ [0.2, 0, 0]
            element = MechanismGeometries.inertial_ellipsoid(inertia)
            @test element.frame === frame
            @test element.transform(SVector(0, 0, 0)) ≈ center_of_mass(inertia).v
            @test element.geometry isa HyperSphere
            @test radius(element.geometry) ≈ 1
        end
    end

    @testset "urdf" begin
        @testset "acrobot" begin
            urdf = "urdf/Acrobot.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            elements = visual_elements(robot, URDFVisuals(urdf))
            @test length(elements) == 3

            element = elements[1]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "base_link"
            @test element.geometry isa Rect
            @test element.geometry.origin ≈ [-0.1, -0.1, -0.1]
            @test element.geometry.widths ≈ [0.2, 0.2, 0.2]
            @test element.color == RGBA(0, 1, 0, 1)
            @test homog(element.transform) ≈ homog(IdentityTransformation())

            element = elements[2]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "upper_link"
            @test element.geometry isa Cylinder
            @test element.geometry.origin ≈ [0, 0, -1.1 / 2]
            @test element.geometry.extremity ≈ [0, 0, 1.1 / 2]
            @test radius(element.geometry) ≈ 0.05
            @test element.color == RGBA(1, 0, 0, 1)
            @test homog(element.transform) ≈ homog(Translation(0, 0, -0.5))

            element = elements[3]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "lower_link"
            @test element.geometry isa Cylinder
            @test element.geometry.origin ≈ [0, 0, -2.1 / 2]
            @test element.geometry.extremity ≈ [0, 0, 2.1 / 2]
            @test radius(element.geometry) ≈ 0.05
            @test element.color == RGBA(0, 0, 1, 1)
            @test homog(element.transform) ≈ homog(Translation(0, 0, -1))
        end

        @testset "acrobot with fixed elbow" begin
            urdf = "urdf/Acrobot_fixed.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            rbd.remove_fixed_tree_joints!(robot)
            elements = visual_elements(robot, URDFVisuals(urdf))
            @test length(elements) == 3

            element = elements[1]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "world"
            @test element.geometry isa Rect
            @test element.geometry.origin ≈ [-0.1, -0.1, -0.1]
            @test element.geometry.widths ≈ [0.2, 0.2, 0.2]
            @test element.color == RGBA(0, 1, 0, 1)
            @test homog(element.transform) ≈ homog(IdentityTransformation())

            element = elements[2]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "upper_link"
            @test element.geometry isa Cylinder
            @test element.geometry.origin ≈ [0, 0, -1.1 / 2]
            @test element.geometry.extremity ≈ [0, 0, 1.1 / 2]
            @test radius(element.geometry) ≈ 0.05
            @test element.color == RGBA(1, 0, 0, 1)
            @test homog(element.transform) ≈ homog(Translation(0, 0, -0.5))

            element = elements[3]
            @test string(rbd.body_fixed_frame_to_body(robot, element.frame)) == "upper_link"
            @test element.geometry isa Cylinder
            @test element.geometry.origin ≈ [0, 0, -2.1 / 2]
            @test element.geometry.extremity ≈ [0, 0, 2.1 / 2]
            @test radius(element.geometry) ≈ 0.05
            @test element.color == RGBA(0, 0, 1, 1)
            @test homog(element.transform) ≈ homog(Translation(0, 0, -1))
        end

        @testset "acrobot submechanism" begin
            urdf = "urdf/Acrobot.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            s = submechanism(robot, bodies(robot)[3])
            elements = visual_elements(s, URDFVisuals(urdf))
            @test length(elements) == 2
            @test string(rbd.body_fixed_frame_to_body(robot, elements[1].frame)) == "upper_link"
            @test string(rbd.body_fixed_frame_to_body(robot, elements[2].frame)) == "lower_link"
        end

        @testset "ground plane" begin
            urdf = "urdf/ground_plane.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            elements = visual_elements(robot, URDFVisuals(urdf; tag="collision"))
            @test length(elements) == 2
            @test elements[1].frame === elements[2].frame
            @test elements[1].geometry isa Rect
            @test elements[2].geometry isa HyperPlane
            @test elements[2].geometry.normal == Vec(0., 0, 1)
            @test homog(elements[2].transform) ≈ homog(Translation(0, 0, 0.025))
        end

        @testset "file with missing mesh" begin
            urdf = "urdf/missing_meshfile.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            elements = visual_elements(robot, URDFVisuals(urdf; tag="collision"))
        end

        @testset "link_colors keyword argument" begin
            urdf = "urdf/Acrobot.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            link = last(bodies(robot))
            for override_color in [RGBA(0.1f0, 0.2f0, 0.3f0, 0.4f0), BGR(0.1, 0.2, 0.3)]
                link_colors = Dict(string(link) => override_color)
                elements_override = visual_elements(robot, URDFVisuals(urdf; link_colors=link_colors))
                elements_base = visual_elements(robot, URDFVisuals(urdf))
                for (element_override, element_base) in zip(elements_override, elements_base)
                    if element_base == last(elements_base)
                        @test RGBA{Float32}(element_override.color) == RGBA{Float32}(override_color)
                    else
                        @test element_override.color == element_base.color
                    end
                end
            end
        end

        @testset "anymal" begin
            # The main reason for testing this URDF is that it has `scale` tags for the meshes.
            urdf = "urdf/anymal_dummy_meshes.urdf"
            robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
            elements = visual_elements(robot, URDFVisuals(urdf; package_path=["urdf"]))
            @test length(elements) == 17
            for element in elements
                _, s, _ = svd(element.transform.linear)
                expected_scale = 1e-3
                @test all(x -> isapprox(x, expected_scale; atol=1e-10), s)
            end
        end
    end

    @testset "valkyrie" begin
        robot = ValkyrieRobot.Valkyrie().mechanism
        visual_elements(robot, Skeleton())
        visual_elements(robot, Skeleton(inertias=false, randomize_colors=true))
        visual_elements(robot, URDFVisuals(ValkyrieRobot.urdfpath(), package_path=[ValkyrieRobot.packagepath()]))
    end

end


text:
module LinLogQuantization

    export  LinQuantArray, LogQuantArray,
    LinQuant8Array, LinQuant16Array, LinQuant24Array, LinQuant32Array,
    LogQuant8Array, LogQuant16Array, LogQuant24Array, LogQuant32Array,
    minpos

    # enable UInt24 support
    import BitIntegers
    BitIntegers.@define_integers 24

    export UInt24

    include("linquantarrays.jl")
    include("logquantarrays.jl")

end


text:
# Author:
# Jürgen Vorloeper, Hochschule Ruhr West, 2018

# code partially based on Julia package "Gaston"
# see: https://github.com/mbaz/Gaston.jl

# verify gnuplot is present on this system
if !success(`gnuplot --version`)
    error("*** Error: gnuplot is not available on this system.")
end

type GnuplotState

running::Bool  # true when gnuplot is already running
fid
figs::Array  # list of identifiers to all figures

function GnuplotState(running, fid, figs)

  new(running,fid,figs)

end

end

function GnuPlotOutput_Init()

global gnuplot_state

pin = Base.Pipe()
r = spawn(`gnuplot`, pin)

gnuplot_state.running = true
gnuplot_state.fid = (pin, r)

w = write(pin, string("set terminal wxt","\n"))
if !(w > 0)
  println("Something went wrong writing to gnuplot STDIN.")
  return
end

return (pin, r)

end

# close gnuplot pipe
function GnuPlotOutput_Close(x...)
    global gnuplot_state

    if gnuplot_state.running
        close(gnuplot_state.fid[1])
    end

    # reset gnuplot_state
    gnuplot_state.running = false
    gnuplot_state.figs = Any[]

    return 0
end

function GnuPlotOutput_Direct(x::Vector{Float64}, y::Vector{Float64})

global gnuplot_state

if gnuplot_state.running == false
  GnuPlotOutput_Init()
end

@assert(length(x) == length(y), "Vectors x and y must have the same number of elements")

s = string("set xlabel \"x\" ");
write(gnuplot_state.fid[1], string(s,"\n"))

s = string("set ylabel \"y\" ");
write(gnuplot_state.fid[1], string(s,"\n"))

s = string("plot [",x[1],":",x[end],"] '-' notitle with lines")
write(gnuplot_state.fid[1], string(s,"\n"))

for i=1:length(x)
  s = string(i, "  ", y[i])
  write(gnuplot_state.fid[1], string(s,"\n"))
end

write(gnuplot_state.fid[1],string("e","\n"))
flush(gnuplot_state.fid[1])

end

function GnuPlotOutput(x::Vector{Float64}, y::Vector{Float64}, fname::String, title::String)

if length(x)!=length(y)
			println("*** Error")
			return -1
end

fid = open(fname*".gnu","w")

if fid==0
  error("*** GnuPlotOutput: Unable to open file")
end

@printf(fid, "set title \"%s\" \n", title)
@printf(fid, "set xlabel \"x\"\n")
@printf(fid, "set ylabel \"y\"\n")
@printf(fid, "plot \"%s.dat\" notitle with lines\n", fname)
@printf(fid, "reset\n")
close(fid)

fid = fopen(fname*".dat", "w")
if fid==0
    error("*** GnuPlotOutput: Unable to open file")
end

for k=1:length(x)
			@printf(fid, "%10.8e %10.8e\n", x[k], y[k])
end

close(fid)

println("GnuPlotOutput: output in files '%s.dat + .gnu'\n", fname)
end

max_stars_repo_name: j-vorloeper/HRW

max_issues_repo_name: j-vorloeper/HRW

max_forks_repo_name: j-vorloeper/HRW


text:
# Files with testcode contain ordinary Symata code, except that line
# that start a new expression and begin with "T " are tested The
# expression should return true or false

## FIXME: remove most of all of these imported symbols. Use fully qualified symbols
import Symata
import Symata: @sym, @ex, @exsimple, SymataPlainTest, print_test_results

## Standard test type
## TODO: set this in test.jl so printing test results is an option
symatatest = SymataPlainTest()

# "context_test.sj",
function runalltests()
    for f in (
        "symbol_test.sj",
        "block_test.sj",
        "ntable_test.sj",
        "flowcontrol_test.sj",
        "input_test.sj",
        "output_test.sj",
        "numeric_test.sj",
        "sympy_test.sj",
        "pattern_test.sj",
        "integral_derivative_test.sj",
        "evaluation_test.sj",
        "list_test.sj",
        "math_function_test.sj",
        "numcomb_test.sj",
        "ast_translation_test.sj",
        "system_test.sj",
        "stack_exchange_test.sj",
        "part_test.sj",
        "shifrin_test.sj",
        "function_test.sj",
        "up_downvalue_test.sj",
        "linalg_test.sj",
        "string_test.sj",
        "julia_interface_test.sj",
        "latex_test.sj",
        "measurement_test.sj",
        "io_test.sj",
        # "old_comparison_test.sj"
        "newcomparison_test.sj",
        "predicate_test.sj",
        "trig_test.sj",
        ## FIXME! reimplement the following after change in comparison code
        # "context_test.sj",
        "mxpr_test.sj",
        "algebraic_transformation_test.sj",
        "attribute_test.sj",
        "module_test.sj",
        "orderless_test.sj",
        "arithmetic_test.sj",
        "expression_test.sj",
        "spec_fun_test.sj"
    )
        Symata.runtest(symatatest, f)
        println(stdout, "Done testing $f")
        Symata.symataevaluate(Meta.parse("Apply(ClearAll, UserSyms())"),
                              Symata.EvaluateJuliaSyntaxSimple())
    end
end

# Test once with Int as default integer type
save_biginput_state = Symata.setkerneloptions(:bigint_input, false)

try
    runalltests()
catch e
    warn("Failed running Symata tests")
    rethrow(e)
finally
    Symata.setkerneloptions(:bigint_input, save_biginput_state)
end

# Test once with BigInt as default integer type

println()

setkerneloptions(:bigint_input, true)

try
    runalltests()
catch
    warn("Failed running Symata tests")
finally
    setkerneloptions(:bigint_input, save_biginput_state)
end

@exsimple ClearAll(testUserSyms)
print_test_results(symatatest)


text:
"""
    rna(R::AbstractRecurrenceMatrix)
    rna(args...; kwargs...)

Calculate a set of Recurrence Network parameters.
The input `R` can be a symmetric recurrence matrix that is
interpreted as the adjacency matrix of an undirected complex network,
such that linked vertices are neighboring points in the phase space.

Alternatively, the inputs can be a graph object or any
valid inputs to the `SimpleGraph` constructor of the
[LightGraphs](https://github.com/JuliaGraphs/LightGraphs.jl) package.

## Return

The returned value is a dictionary that contains the following entries,
with the corresponding global network properties[1, 2]:
* `:density`: edge density, approximately equivalent to the global recurrence rate in the phase space.
* `:transitivity`: network transitivity, which describes the
global clustering of points following Barrat's and Weigt's formulation [3].
* `:averagepath`: mean value of the shortest path lengths taken over
all pairs of connected vertices, related to the average separation
between points in the phase.
* `:diameter`: maximum value of the shortest path lengths between
pairs of connected vertices, related to the phase space diameter.

## References

[1]: R.V. Donner *et al.* "Recurrence networks — a novel paradigm for nonlinear time series analysis",
*New Journal of Physics* 12, 033025 (2010)
[DOI:10.1088/1367-2630/12/3/033025](https://doi.org/10.1088/1367-2630/12/3/033025)

[2]: R.V. Donner *et al.*, The geometry of chaotic dynamics — a complex network perspective,
*Eur. Phys. J.* B 84, 653–672 (2011)
[DOI:10.1140/epjb/e2011-10899-1](https://doi.org/10.1140/epjb/e2011-10899-1)

[3]: A. Barrat & M. Weight, "On the properties of small-world network models",
*The European Physical Journal B* 13, 547–560 (2000)
[DOI:10.1007/s100510050067](https://doi.org/10.1007/s100510050067)
"""
function rna(args...; kwargs...)
    graph = SimpleGraph(args...; kwargs...)
    return Dict{Symbol, Float64}(
        :density => density(graph),
        :transitivity => global_clustering_coefficient(graph),
        :averagepath => averagepath(graph),
        :diameter => diameter(graph)
    )
end

"""
    averagepath(graph)

Calculates average minimum path length for a SimpleGraph.
Uses Donner, 2010, Eq. 26 on pg. 18.
"""
function averagepath(graph::SimpleGraph)
    num_verts = nv(graph)

    dist_mat = floyd_warshall_shortest_paths(graph).dists
    @. dist_mat[dist_mat == typemax(dist_mat)] = 0.

    return sum(dist_mat) / (num_verts * (num_verts - 1))
end


text:
##################################################

## promote up to symbolic so that math ops work
promote_rule(::Type{T}, ::Type{S})  where {T<:SymbolicObject, S<:Number}= T
Base.promote_type(::Type{Irrational{T}}, ::Type{Sym}) where {T} = Sym
Base.promote_op(::T, ::Type{S}, ::Type{Sym}) where {T, S <: Number} = Sym
Base.promote_op(::T, ::Type{Sym}, ::Type{S}) where {T, S <: Number} = Sym
Base.promote_op(::T, ::Type{Sym}, ::Type{Sym}) where {T} = Sym # This helps out linear algebra conversions

## Conversion
convert(::Type{T}, o::PyCall.PyObject) where {T <: SymbolicObject} = T(o)
convert(::Type{PyObject}, s::Sym) = s.__pyobject__

## rational
convert(::Type{T}, x::Rational) where {T<:SymbolicObject} = sympy.Rational(x.num, x.den)::T

## big. Need mpmath installed separately -- not as a SymPy module as that is how it is called in PyCall
convert(::Type{T}, x::BigFloat) where {T<:SymbolicObject} = Sym(PyCall.PyObject(x))::T
convert(::Type{Sym}, x::Complex{BigFloat}) = Sym(PyCall.PyObject(x))::Sym

## real
convert(::Type{S}, x::T) where {S<:SymbolicObject, T <: Real}= sympify(x)::S
convert(::Type{T}, x::Sym) where {T <: Real} = convert(T, PyObject(x))

## complex
## cf math.jl for `complex` of a value
## IM is SymPy's "i" (synpy.I, not Python's
## Sym(PyCall.PyObject(im)) which gives 1j.
function convert(::Type{Sym}, x::Complex)
    y = ifelse(isa(x, Complex{Bool}), real(x) + imag(x) * im, x)
    real(y) + imag(y) * IM
end
convert(::Type{Complex{T}}, x::Sym) where {T} = complex(map(x -> convert(T, x), x.as_real_imag())...)

## Irrationals
Base.convert(::Type{Sym}, x::Irrational{:π}) = PI
Base.convert(::Type{Sym}, x::Irrational{:ℯ}) = sympy.exp(1)
Base.convert(::Type{Sym}, x::Irrational{:γ}) = Sym(sympy.EulerGamma)
Base.convert(::Type{Sym}, x::Irrational{:catalan}) = Sym(sympy.Catalan)
Base.convert(::Type{Sym}, x::Irrational{:φ}) = (1 + Sym(5)^(1//2))/2

## utility functions to test type of value
function is_integer(x::Sym)

    pycall_hasproperty(x, :is_integer) && return x.is_integer
    x.__class__.__name__ == "int" && return true

    return false
end

function is_rational(x::Sym)
    pycall_hasproperty(x, :is_rational) && return x.is_rational

    return false
end

function is_real(x::Sym)
    pycall_hasproperty(x, :is_real) && return x.is_real
    x.__class__.__name__ == "real" && return true
    x.__class__.__name__ == "mpg" && return true

    return false
end

function is_complex(x::Sym)
    pycall_hasproperty(x, :is_complex) && return x.is_complex
    x.__class__.__name__ == "complex" && return true
    x.__class__.__name__ == "mpc" && return true

    return false
end

"""

Convert a `Sym` value to a numeric Julian value.

The `N` function of SymPy is an alias for `evalf`. Within SymPy, either may be used to
find numeric values from symbolic values.

For example, symbolic roots can be computed numerically, even if not
available symbolically, by calling `N` on the values.

Using `SymPy` within `Julia` makes having two such functions useful:

* one to do the equivalent of SymPy's `evalf` function
* one to *also* convert these expressions back into `Julia` objects.

We use `N` to return a `Julia` object and `evalf` to return a symbolic
object.

Examples:
```
x = Sym("x")
p = subs(x, x, pi)
N(p)                            # float version of pi
evalf(p, 60)                    # 60 digits of pi, as a symbolic value
N(p, 60)                        # when a precision is given, "Big" values are returned
r = subs(x,x,1.2)
N(r)                            # float
q = subs(x, x, 1//2)
N(q)                            # 1//2
z = solve(x^2 + 1)[1]           # -ⅈ
N(z)                            # 0 - 1im
evalf(z)

rts = solve(x^5 - x + 1)
[N(r) for r in rts]             # numeric solutions to quintic
```

The `evalf` function is similar, though it leaves the expression as a symbolic object.
This breaks the similarity of N and evalf for sympy users.

Returns the value unchanged when it has free symbols.

`N` is type unstable.

"""
function N(x::Sym)

    length(free_symbols(x)) > 0 && return x
    # many different possible types, and not all have some nice property
    # python int
    # mpmath int
    # SymPy big int
    # python float
    # Pi, Half, Rational

    for  (u,v) in sympy_core_numbers
        if pycall_hasproperty(x, :__class__)
            if x.__class__.__name__ == string(u)
                return v
            end
        end
    end
    if is_(:real, x)
        if is_(:zero, x)
            return 0
        elseif is_(:infinite, x)
            return (x.is_negative ? -1 : 1) * Inf
        elseif is_(:integer, x)
            u = abs(x)
            if sympy.Le(u, typemax(Int)) == Sym(true)
                return convert(Int, x)
#            elseif sympy.Le(u, typemax(Int128)) == Sym(true)
#                # no PyCall support for this conversion?
#                return convert(Int128, x)
            else
                return convert(BigInt, x)
            end
        elseif x.__class__.__name__ == "Float"
            if x._prec <= 64
                return convert(Float64, x)
            else
                return convert(BigFloat, x)
            end
        elseif is_(:rational, x)
            return N(numer(x)) // N(denom(x))
        elseif pycall_hasproperty(x, :args) && length(x.args) > 1
            def_precision_decimal = ceil(Int, log10(big"2"^Base.MPFR.DEFAULT_PRECISION.x))
            convert(BigFloat, x.evalf(def_precision_decimal))
        elseif length(x.args) > 0
            return N(x.evalf())
        else
            return convert(BigFloat, x)
        end
    elseif x.__class__.__name__ == "ComplexRootOf"
        u = x.evalf(16)
        return N(u)
    elseif is_(:complex, x)
        return complex(N(sympy.re(x)), N(sympy.im(x)))
    elseif is_(:imaginary, x)
        return complex(0, N(sympy.im(x)))
    elseif x.__class__.__name__ == "int"
        convert(Int, x)
    elseif x.__class__.__name__ == "float"
        convert(Float64, x)
    elseif x.__class__.__name__ == "complex"
        return complex(N(sympy.re(x)), N(sympy.im(x)))
    elseif x.__class__.__name__ == "mpf"
        convert(Float64, x)
    elseif x.__class__.__name__ == "mpc"
        return complex(N(sympy.re(x)), N(sympy.im(x)))
    elseif Eq(x,Sym(true)) == Sym(true)
        return true
    elseif Eq(x, Sym(false)) == Sym(true)
        return false
    else
        @info "FAILED to find type for $x. Please report"
        x
    end
end
N(x::Number) = x  # implies N(x::Sym) = x if ...
N(m::AbstractArray{Sym}) = map(N, m)

# special case numbers in sympy.core.numbers
sympy_core_numbers = ((:Zero, 0),
                      (:One, 1),
                      (:NegativeOne, -1),
                      (:Half, 1//2),
                      (:NaN, NaN),
                      (:Exp1, ℯ),
                      (:ImaginaryUnit, im),
                      (:Pi, pi),
                      (:EulerGamma, Base.MathConstants.eulergamma),
                      (:Catalan, Base.MathConstants.catalan),
                      (:GoldenRation, Base.MathConstants.golden),
                      (:TribonacciConstant, big(1)/3 + (-big(3)*sqrt(big(33)) + 19)^(1//3)/3 + (3*sqrt(big(33)) + 19)^(1//3)/3))

# fix me XXX
"""
`N` can take a precision argument.

When given as an integer greater than 16, we try to match the digits of accuracy using `BigFloat` precision on conversions to floating point.

"""
function N(x::Sym, digits::Int)

    ## check
    digits <= 16 && return(N(x))
    if is_integer(x) == nothing
        out = x.evalf(digits)
        return( N(out, digits) )
    end

    ex = x.evalf(digits)
    if is_integer(x)
        return(convert(BigInt, x))
    elseif is_rational(x)
        return N(numer(x)) / N(denom(x))
    elseif is_real(x) == true
        p = round(Int,log2(10)*digits)

        out = setprecision(p) do
            convert(BigFloat, ex)
        end
        return(out)
    elseif is_complex(x) == true
        r, i = ex.as_real_imag()
        u, v = promote(N(r, digits), N(i, digits))
        return(Complex(u, v))
    end

    throw(DomainError())
end

##################################################
##
## infix logical operators

## XXX Experimental! Not sure these are such a good idea ...
## but used with piecewise
import Base: &, |, !
Base.:&(x::Sym, y::Sym) = PyCall.pycall(PyObject(x).__and__, Sym, y)
Base.:|(x::Sym, y::Sym) =  PyCall.pycall(PyObject(x).__or__, Sym, y)
!(x::Sym)         =       PyCall.pycall(PyObject(x).__invert__, Sym)::Sym

## use ∨, ∧, ¬ for |,&,! (\vee<tab>, \wedge<tab>, \neg<tab>)
∨(x::Sym, y::Sym) = x | y
∧(x::Sym, y::Sym) = x & y
¬(x::Sym) = !x
export ∨, ∧, ¬

## In SymPy, symbolic equations are not represented by `=` or `==`
## rather ther function `Eq` is used. Here we use the unicode
## `\Equal<tab>` for an infix operator. There are also unicode values to represent analogs of `<`, `<=`, `>=`. `>`. These are

## * `<`  is `\ll<tab>`
## * `<=` is `\leqq<tab>
## * `==` is `\Equal<tab>`
## * `>=` is `\geqq<tab>`
## * `>`  is `\gg<tab>`

## Instead we have:
## We use unicode for visual appeal of infix operators, but the Lt, Le, Eq, Ge, Gt are the proper way:

"This is `\\ll<tab>` mapped as an infix operator to `Lt`"
(≪)(a::Sym, b::Sym) = Lt(a,b)  # \ll<tab>
(≪)(a::Sym, b::Number) = Lt(a,Sym(b)) # \ll<tab>
(≪)(a::Number, b::Sym) = Lt(Sym(a),b) # \ll<tab>

## could just do this, but it would interfere with other uses outside of SymPy
## (≪)(a::Number, b::Number) = Lt(promote(a,b)...)  # \ll<tab>

"This is `\\leqq<tab>` mapped as an infix operator to `Le`"
(≦)(a::Sym, b::Sym) = Le(a,b)   # \ll<tab>
(≦)(a::Sym, b::Number) = Le(a,Sym(b))  # \ll<tab>
(≦)(a::Number, b::Sym) = Le(Sym(a),b)   # \ll<tab>

"This is `\\gg<tab>` mapped as an infix operator to `Gt`"
(≫)(a::Sym, b::Sym) = Gt(a,b) |> asBool
(≫)(a::Sym, b::Number) = Gt(a,Sym(b))
(≫)(a::Number, b::Sym) = Gt(Sym(a),b)

"This is `\\geqq<tab>` mapped as an infix operator to `Ge`"
(≧)(a::Sym, b::Sym) = Ge(a,b) |> asBool
(≧)(a::Sym, b::Number) = Ge(a,Sym(b))
(≧)(a::Number, b::Sym) = Ge(Sym(a),b)

"For infix `Eq` one can use \\Equal<tab> unicode operator"
(⩵)(a::Sym, b::Sym) = Eq(a,b)  # \Equal<tab>
(⩵)(a::Sym, b::Number) = Eq(a,Sym(b))  # \Equal<tab>
(⩵)(a::Number, b::Sym) = Eq(Sym(a),b)  # \Equal<tab>

export ≪,≦,⩵,≧,≫


text:
const jennrich_sampson = let res_init=zeros(10), jac_init=zeros(10,2), x_init=[0.3, 0.4]

    function res(x, r)
        for i = 1:10
            ti = Float64(i)
            r[i] = 2.0e0 + 2.0e0*ti - (exp(ti*x[1]) + exp(ti*x[2]))
        end
        return r
    end

    function jac(x, j)
        for i = 1:10
            ti = Float64(i)
            j[i,1] = -ti * exp(ti*x[1])
            j[i,2] = -ti * exp(ti*x[2])
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("jennrich_sampson", res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end