text:
"""
    constants
"""

const GSL_SQRT_DBL_EPSILON = 1.4901161193847656e-08
const GSL_LOG_DBL_MIN = -7.0839641853226408e02

const synch_c0 = π / sqrt(3.0)

const c01 = 0.2257913526447274323630976
const cond1 = 2 * sqrt(2.0) * GSL_SQRT_DBL_EPSILON
const cond3 = -8.0 * GSL_LOG_DBL_MIN / 7.0

const synchrotron1_data = [ 30.364682982501076273,
                            17.079395277408394574,
                            4.560132133545072889,
                            0.549281246730419979,
                            0.372976075069301172e-01,
                            0.161362430201041242e-02,
                            0.481916772120371e-04,
                            0.10512425288938e-05,
                            0.174638504670e-07,
                            0.22815486544e-09,
                            0.240443082e-11,
                            0.2086588e-13,
                            0.15167e-15
                            ]

const synchrotron1a_data = [ 2.1329305161355000985,
                             0.741352864954200240e-01,
                             0.86968099909964198e-02,
                             0.11703826248775692e-02,
                             0.1645105798619192e-03,
                             0.240201021420640e-04,
                             0.35827756389389e-05,
                             0.5447747626984e-06,
                             0.838802856196e-07,
                             0.13069882684e-07,
                             0.2053099071e-08,
                             0.325187537e-09,
                             0.517914041e-10,
                             0.83002988e-11,
                             0.13352728e-11,
                             0.2159150e-12,
                             0.349967e-13,
                             0.56994e-14,
                             0.9291e-15,
                             0.152e-15,
                             0.249e-16,
                             0.41e-17,
                             0.7e-18
                            ]

const synchrotron2_data = [  0.4490721623532660844,
                             0.898353677994187218e-01,
                             0.81044573772151290e-02,
                             0.4261716991089162e-03,
                             0.147609631270746e-04,
                             0.3628633615300e-06,
                             0.66634807498e-08,
                             0.949077166e-10,
                             0.1079125e-11,
                             0.10022e-13,
                             0.77e-16,
                             0.5e-18
                            ]

@inline @fastmath function cheb_eval(coeff::Vector{<:Real}, order::Integer, x::Real)

    d = 0.0
    dd = 0.0

    x2 = 2x

    @inbounds for j = order:-1:2
        temp = d
        d = x2 * d - dd + coeff[j]
        dd = temp
    end

    d = x * d - dd + 0.5 * coeff[1]

    return d
end

"""
    F(x::Real)

Computes the first synchrotron function at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.

``F(x) = x \\int_x^\\infty K_{\\frac{5}{3}}(t) dt``

"""
@inline function F(x::Real)

    if x < cond1

        z = x^(1/3)
        cf = 1 - 8.43812762813205e-01 * z * z
        return 2.14952824153447863671 * z * cf

    elseif x <= 4.0

        px = ∛(x)
        px11 = px*px*px*px*px*px*px*px*px*px*px
        t = x^2 / 8.0 - 1.0
        result_c1 = cheb_eval( synchrotron1_data, 13, t )

        result_c2 = cheb_eval( synchrotron2_data, 12, t )

        return px * result_c1 - px11 * result_c2 - synch_c0 * x

    elseif x < cond3

        t = (12.0 - x) / (x + 4.0)

        result_c1 = cheb_eval( synchrotron1a_data, 23, t )

        return sqrt(x) * result_c1 * exp(c01 - x)

    else
        return 0.0
    end
end

"""
    ℱ(x::Real)

Computes the first synchrotron function at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.

``F(x) = x \\int_x^\\infty K_{\\frac{5}{3}}(t) dt``

Wrapper for [`F`](@ref).
"""
ℱ(x::Real) = F(x)


text:
using Pkg

using BinDeps
import CondaBinDeps

if lowercase(get(ENV, "CI", "false")) == "true"

    let basepython = get(ENV, "PYTHON", "python2")
        envpath = joinpath(@__DIR__, "env")
        run(`pip install --user virtualenv`)
        run(`virtualenv --python=$basepython $envpath`)

        if Sys.iswindows()
            python = joinpath(envpath, "Scripts", "python.exe")
        else
            python = joinpath(envpath, "bin", "python2")
        end
        run(`$python -m pip install numpy`)
        run(`$python -m pip install scipy`)
        run(`$python -m pip install matplotlib`)
        run(`$python -m pip install https://github.com/matplotlib/basemap/archive/v1.0.7rel.tar.gz`)
        run(`$python -m pip install git+https://github.com/matplotlib/cmocean`)

        ENV["PYTHON"] = python
        Pkg.build("PyCall")

        if VERSION >= v"0.7.0"
            using Libdl
        end

        function validate_netcdf_version(name,handle)
            f = Libdl.dlsym_e(handle, "nc_inq_libvers")
            #
            # Example
            # banner = "4.5.6.7 of Apr 1 2000 00:00:00"
            banner = unsafe_string(ccall(f,Ptr{UInt8},()))

            # remove the date
            # verstr = "4.5.6.7"
            verstr = split(banner)[1]

            # vernumbers = ["4","5","6","7"]
            vernumbers = split(verstr,r"[.-]")

            # major_minor_patch = ["4","5","6"]
            major_minor_patch = vernumbers[1:min(3,length(vernumbers))]

            # ver = v"4.5.6"
            ver = VersionNumber([parse(Int,s) for s in major_minor_patch]...)

            return ver > v"4.2"
        end

        @BinDeps.setup
        libnetcdf = library_dependency("libnetcdf", aliases = ["libnetcdf4","libnetcdf-7","netcdf"], validate = validate_netcdf_version)

        CondaBinDeps.Conda.add_channel("anaconda")
        provides(CondaBinDeps.Manager, "libnetcdf", libnetcdf)
        provides(AptGet, "libnetcdf-dev", libnetcdf, os = :Linux)
        # provides(Yum, "netcdf-devel", libnetcdf, os = :Linux)

        @BinDeps.install Dict(:libnetcdf => :libnetcdf)

    end
end

ext: jl


text:
function timestep(fishes)
    no_of_fish_to_add = 0
    for (i, fish) in enumerate(fishes)
        if fish == 0
            no_of_fish_to_add += 1
            fishes[i] = 6
        else
            fishes[i] = fish - 1
        end
    end
    for i in 1:no_of_fish_to_add
        fishes = push!(fishes, 8)
    end
    #println(fishes)
    return fishes
end

function timestep2(fishes)
    toadd = zeros(Int, 9)
    for i in 1:9
        if i == 1
            toadd[9] = fishes[i]
            toadd[7] = fishes[i]
        else
            toadd[i-1] += fishes[i]
        end
    end
    return toadd
end

function solve(filename, timesteps = 80)
    #load txt file
    values = readlines(filename*".txt")
    values = split(values[1], ',')
    fishes = parse.(Int, values)
    #println(fishes)
    for i in 1:timesteps
        #print("Day "*string(i))
        fishes = timestep(fishes)
    end
    println("Number of list elements after "*string(timesteps)*" days: "*string(length(fishes)))
end

function solve2(filename, timesteps = 80)
    #load txt file
    values = readlines(filename*".txt")
    values = split(values[1], ',')
    values = parse.(Int, values)
    fishes = zeros(Int, 9)
    for v in values
        fishes[v+1] += 1
    end
    for i in 1:timesteps
        fishes = timestep2(fishes)
    end
    res = 0
    for v in fishes
        res += v
    end
    println("Number of list elements after "*string(timesteps)*" days: "*string(res))
end

solve("smalltest")
solve("input")
solve2("smalltest", 256)
solve2("input", 256)


text:
module Yama
export create_surface, create_trajectory, SurfaceArgs

using Base.Iterators: product
using BSON: @save, @load

using CuArrays
using Flux: params, gradient, gpu, cpu, loadparams!, testmode!
using Flux.Data: DataLoader
using LinearAlgebra: norm, dot
using MultivariateStats: PCA, fit, projection

using Parameters: @with_kw
using ProgressMeter: @showprogress
using Zygote: Params

include("utils.jl")
include("surface.jl")
include("trajectory.jl")

end


text:
using Flux.Optimise
using Flux.Tracker
using Test
@testset "Optimise" begin
  w = randn(10, 10)
  @testset for Opt in [SGD, Nesterov, Momentum, ADAM, AdaMax, RMSProp, ps -> ADAGrad(ps, 0.1), ADADelta, AMSGrad, NADAM]
    w′ = param(randn(10, 10))
    loss(x) = Flux.mse(w*x, w′*x)
    opt = Opt([w′])
    for t=1:10^5
      l = loss(rand(10))
      back!(l)
      opt()
    end
    @test Flux.mse(w, w′) < 0.01
  end
end

@testset "Training Loop" begin
  i = 0
  l = param(1)

  Flux.train!(() -> (sleep(0.1); i += 1; l),
              Iterators.repeated((), 100),
              ()->(),
              cb = Flux.throttle(() -> (i > 3 && stop()), 1))

  @test 3 < i < 50
end


text:
@testset "stopping criteria as controls" begin

    # A stopping criterion that ignores training losses:

    m = SquareRooter(4)
    c = NumberLimit(2)

    IC.train!(m, 3)
    state = IC.update!(c, m, 0, 1)
    @test state == 1
    @test !IC.done(c, state)
    IC.train!(m, 3)
    state = IC.update!(c, m, 0, 2, state)
    @test state == 2
    @test IC.done(c, state)
    report = @test_logs (:info, r"NumberLimit\(2\)") IC.takedown(c, 1, state)
    @test report.done
    @test report.log ==
        "Stop triggered by NumberLimit(2) stopping criterion. "

    # A stopping criterion than uses training losses:

    m = SquareRooter(4)
    c = PQ(k=5)

    # Note that `SquareRooter` does not cache training losses. It only
    # makes the most recent losses available, which is the worst-case
    # scenario.

    IC.train!(m, 3)
    state = IC.update!(c, m, 0, 1)
    train_losses = m.training_losses
    @test reverse(state.training_losses) == train_losses # length 3
    @test !IC.done(c, state)

    IC.train!(m, 2)
    state = IC.update!(c, m, 0, 2, state)
    train_losses = vcat(train_losses, m.training_losses) # length 2+3 = 5

    @test reverse(state.training_losses) == train_losses
    @test !IC.done(c, state)
    report = IC.takedown(c, 1, state)
    @test !report.done
    @test report.log == ""
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestNonlinearRead

using Test
import MathOptInterface

const MOI = MathOptInterface
const NL = MOI.FileFormats.NL

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_next_int()
    model = NL._CacheModel()
    for i in (1, 9, 10, 11, 99)
        io = IOBuffer()
        print(io, " $i\n")
        seekstart(io)
        @test NL._next(Int, io, model) == i
        @test peek(io) == UInt8('\n')
    end
    return
end

function test_next_float()
    model = NL._CacheModel()
    for (s, x) in [
        ("-1", -1.0),
        ("1", 1.0),
        ("-2.1e+1", -21.0),
        ("2.1e+1", 21.0),
        ("-2.1e-01", -0.21),
        ("2.1e-01", 0.21),
        ("-2.1E-01", -0.21),
        ("2.1E-01", 0.21),
        ("-20.1e+0", -20.1),
        ("20.1e+00", 20.1),
        ("  20.1e+00", 20.1),
    ]
        io = IOBuffer()
        print(io, "$s\n")
        seekstart(io)
        @test NL._next(Float64, io, model) == x
        @test peek(io) == UInt8('\n')
    end
    return
end

function test_parse_expr()
    model = NL._CacheModel()
    io = IOBuffer()
    write(io, "o2\nv0\no2\nn2\no2\nv3\nv1\n")
    # (* x1 (* 2 (* x4 x2)))
    seekstart(io)
    x = MOI.VariableIndex.(1:4)
    @test NL._parse_expr(io, model) == :(*($(x[1]), *(2, *($(x[4]), $(x[2])))))
    @test eof(io)
    return
end

function test_parse_expr_nary()
    model = NL._CacheModel()
    io = IOBuffer()
    write(io, "o54\n4\no5\nv0\nn2\no5\nv2\nn2\no5\nv3\nn2\no5\nv1\nn2\n")
    seekstart(io)
    x = MOI.VariableIndex.(1:4)
    @test NL._parse_expr(io, model) ==
          :(+($(x[1])^2, $(x[3])^2, $(x[4])^2, $(x[2])^2))
    @test eof(io)
    return
end

function test_parse_expr_minimum()
    model = NL._CacheModel()
    io = IOBuffer()
    write(io, "o11\n3\nv0\nv1\nv2\n")
    seekstart(io)
    x = MOI.VariableIndex.(1:3)
    @test NL._parse_expr(io, model) == :(min($(x[1]), $(x[2]), $(x[3])))
    @test eof(io)
    return
end

function test_parse_expr_maximum()
    model = NL._CacheModel()
    io = IOBuffer()
    write(io, "o12\n3\nv0\nv1\nv2\n")
    seekstart(io)
    x = MOI.VariableIndex.(1:3)
    @test NL._parse_expr(io, model) == :(max($(x[1]), $(x[2]), $(x[3])))
    @test eof(io)
    return
end

function test_parse_header_binary()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(io, "b3 1 1 0\n")
    seekstart(io)
    @test_throws(
        ErrorException("Unable to parse NL file : unsupported mode b"),
        NL._parse_header(io, model),
    )
    return
end

function test_parse_header_assertion_errors()
    model = NL._CacheModel()
    for header in [
        "g4 1 1 0\n3 3 2 0 0 0\n",
        "g4 1 1 0\n3 3 1 0 0 0\n0 2\n",
        "g4 1 1 0\n3 3 1 0 0 0\n0 0\n1 0\n",
        "g4 1 1 0\n3 3 1 0 0 0\n0 0\n0 1\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n1 0 0 1\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 1 0 1\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 1 1\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n-1 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n0 -1\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n1 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 1\n",
    ]
        io = IOBuffer()
        write(io, header)
        seekstart(io)
        @test_throws(AssertionError, NL._parse_header(io, model))
    end
    return
end

function test_parse_header_assertion_errors()
    model = NL._CacheModel()
    err = ErrorException(
        "Unable to parse NL file : we don't support common exprs",
    )
    for header in [
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 0\n1 0 0 0 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 0\n0 1 0 0 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 0\n0 0 1 0 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 0\n0 0 0 1 0\n",
        "g3 1 1 0\n4 2 1 0 1 0\n2 1\n0 0\n4 0 0\n0 0 0 1\n0 0 0 2 0\n8 4\n0 0\n0 0 0 0 1\n",
    ]
        io = IOBuffer()
        write(io, header)
        seekstart(io)
        @test_throws(err, NL._parse_header(io, model))
    end
    return
end

function test_parse_y_error()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(io, "y\n")
    seekstart(io)
    @test_throws(
        ErrorException("Unable to parse NL file: unhandled header y"),
        NL._parse_section(io, model),
    )
    return
end

function test_parse_O()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """
O0 0# can stick a comment anywhere
o2
v0  # can stick a comment anywhere
o2
v2
o2
v3
v1
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.sense == MOI.MIN_SENSE
    x = MOI.VariableIndex.(1:4)
    @test model.objective == :($(x[1]) * ($(x[3]) * ($(x[4]) * $(x[2]))))
    return
end

function test_parse_O_max()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(io, "O0 1\nv0\n")
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.sense == MOI.MAX_SENSE
    x = MOI.VariableIndex(1)
    @test model.objective == :(+$x)
    return
end

function test_parse_x()
    model = NL._CacheModel()
    NL._resize_variables(model, 5)
    io = IOBuffer()
    write(
        io,
        """
x3
0 1.1# can stick a comment anywhere
3 2.2
2 3.3  # can stick a comment anywhere
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.variable_primal == [1.1, 0.0, 3.3, 2.2, 0.0]
    return
end

function test_parse_d()
    model = NL._CacheModel()
    io = IOBuffer()
    write(
        io,
        """
d3
0 1.1# can stick a comment anywhere
3 2.2
2 3.3  # can stick a comment anywhere
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    return
end

function test_parse_r()
    model = NL._CacheModel()
    NL._resize_constraints(model, 5)
    io = IOBuffer()
    write(
        io,
        """
r# can stick a comment anywhere
1 3.3
3
0 1.1 2.2  # can stick a comment anywhere
4 5.5
2 4.4# can stick a comment anywhere
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.constraint_lower == [-Inf, -Inf, 1.1, 5.5, 4.4]
    @test model.constraint_upper == [3.3, Inf, 2.2, 5.5, Inf]
    return
end

function test_parse_b()
    model = NL._CacheModel()
    NL._resize_variables(model, 5)
    io = IOBuffer()
    write(
        io,
        """
b# can stick a comment anywhere
1 3.3
3# can stick a comment anywhere
0 1.1 2.2
4 5.5
2 4.4 # can stick a comment anywhere
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.variable_lower == [-Inf, -Inf, 1.1, 5.5, 4.4]
    @test model.variable_upper == [3.3, Inf, 2.2, 5.5, Inf]
    return
end

function test_parse_k()
    model = NL._CacheModel()
    NL._resize_variables(model, 3)
    io = IOBuffer()
    write(
        io,
        """
k
2 # can stick a comment anywhere
4
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    return
end

function test_parse_J()
    model = NL._CacheModel()
    NL._resize_constraints(model, 3)
    io = IOBuffer()
    write(
        io,
        """
J1 2  # can stick a comment anywhere
2 1.1
3 2.2  # can stick a comment anywhere
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.constraints[1] == :()
    x, y = MOI.VariableIndex(3), MOI.VariableIndex(4)
    @test model.constraints[2] == :(1.1 * $x + 2.2 * $y)
    @test model.constraints[3] == :()
    return
end

function test_parse_J_zeros()
    model = NL._CacheModel()
    NL._resize_constraints(model, 1)
    io = IOBuffer()
    write(
        io,
        """
J0 2
0 0
1 0
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.constraints[1] == :()
    return
end

function test_parse_C_J()
    model = NL._CacheModel()
    NL._resize_constraints(model, 1)
    io = IOBuffer()
    write(
        io,
        """
C0
o2
v0
v1
J0 2
0 1.1
1 2.2
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    NL._parse_section(io, model)
    @test eof(io)
    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
    @test model.constraints[1] == :((1.1 * $x + 2.2 * $y) + $x * $y)
    return
end

function test_parse_J_C()
    model = NL._CacheModel()
    NL._resize_constraints(model, 1)
    io = IOBuffer()
    write(
        io,
        """
J0 2
0 1.1
1 2.2
C0
o2
v0
v1
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    NL._parse_section(io, model)
    @test eof(io)
    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
    @test model.constraints[1] == :((1.1 * $x + 2.2 * $y) + $x * $y)
    return
end

function test_parse_G()
    model = NL._CacheModel()
    NL._resize_constraints(model, 3)
    io = IOBuffer()
    write(
        io,
        """
G0 2 # can stick a comment anywhere
2 1.1  # can stick a comment anywhere
3 2.2
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    x, y = MOI.VariableIndex(3), MOI.VariableIndex(4)
    @test model.objective == :(1.1 * $x + 2.2 * $y)
    return
end

function test_parse_G_zeros()
    model = NL._CacheModel()
    NL._resize_constraints(model, 3)
    io = IOBuffer()
    write(
        io,
        """
G0 3
0 0
1 0
2 0
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    @test eof(io)
    @test model.objective == :()
    return
end

function test_parse_O_G()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """
O0 0
o2
v0
v1
G0 2
0 1.1
1 2.2
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    NL._parse_section(io, model)
    @test eof(io)
    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
    @test model.objective == :((1.1 * $x + 2.2 * $y) + $x * $y)
    return
end

function test_parse_G_O()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """
G0 2
0 -1.1
1 2.2e+01
O0 0
o2
v0
v1
""",
    )
    seekstart(io)
    NL._parse_section(io, model)
    NL._parse_section(io, model)
    @test eof(io)
    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
    @test model.objective == :((-1.1 * $x + 22.0 * $y) + $x * $y)
    return
end

function test_hs071()
    model = NL.Model()
    open(joinpath(@__DIR__, "data", "hs071.nl"), "r") do io
        return read!(io, model)
    end
    model_print = sprint(print, model)
    _in_ = @static Sys.iswindows() ? "in" : "∈"
    for line in [
        "(1.1 * v[1] + 1.4 * v[2] + 1.2 * v[3] + 1.3 * v[4]) + +2.0",
        "v[1] >= 1.1",
        "v[2] >= 1.4",
        "v[3] >= 1.2",
        "v[4] >= 1.3",
        "v[1] <= 5.1",
        "v[2] <= 5.4",
        "v[3] <= 1.0",
        "v[4] <= 5.3",
        "v[3] $_in_ ℤ",
        "v[4] $_in_ ℤ",
        "v[1] * (v[3] * (v[4] * v[2])) >= 25.0",
        "v[1] ^ 2.0 + v[3] ^ 2.0 + v[4] ^ 2.0 + v[2] ^ 2.0 == 40.0",
    ]
        @test occursin(line, model_print)
    end
    return
end

function test_parse_header_integrality_obj()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """g3 0 1 0	# problem test_simple
2 1 1 0 0	# vars, constraints, objectives, ranges, eqns
1 1	# nonlinear constraints, objectives
0 0	# network constraints: nonlinear, linear
1 2 0	# nonlinear vars in constraints, objectives, both
0 0 0 1	# linear network variables; functions; arith, flags
0 0 0 0 0	# discrete variables: binary, integer, nonlinear (b,c,o)
1 1	# nonzeros in Jacobian, gradients
0 0	# max name lengths: constraints, variables
0 0 0 0 0	# common exprs: b,c,o,c1,o1\n""",
    )
    seekstart(io)
    NL._parse_header(io, model)
    @test model.variable_type[1] == NL._CONTINUOUS
    @test model.variable_type[2] == NL._CONTINUOUS
    return
end

function test_parse_header_integrality_obj_int()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """g3 0 1 0	# problem test_simple
2 1 1 0 0	# vars, constraints, objectives, ranges, eqns
1 1	# nonlinear constraints, objectives
0 0	# network constraints: nonlinear, linear
1 2 0	# nonlinear vars in constraints, objectives, both
0 0 0 1	# linear network variables; functions; arith, flags
0 0 0 0 1	# discrete variables: binary, integer, nonlinear (b,c,o)
1 1	# nonzeros in Jacobian, gradients
0 0	# max name lengths: constraints, variables
0 0 0 0 0	# common exprs: b,c,o,c1,o1\n""",
    )
    seekstart(io)
    NL._parse_header(io, model)
    @test model.variable_type[1] == NL._CONTINUOUS
    @test model.variable_type[2] == NL._INTEGER
    return
end

function test_parse_header_integrality_both()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """g3 0 1 0	# problem test_simple
2 1 1 0 0	# vars, constraints, objectives, ranges, eqns
1 1	# nonlinear constraints, objectives
0 0	# network constraints: nonlinear, linear
2 1 1	# nonlinear vars in constraints, objectives, both
0 0 0 1	# linear network variables; functions; arith, flags
0 0 0 0 0	# discrete variables: binary, integer, nonlinear (b,c,o)
1 1	# nonzeros in Jacobian, gradients
0 0	# max name lengths: constraints, variables
0 0 0 0 0	# common exprs: b,c,o,c1,o1\n""",
    )
    seekstart(io)
    NL._parse_header(io, model)
    @test model.variable_type[1] == NL._CONTINUOUS
    @test model.variable_type[2] == NL._CONTINUOUS
    return
end

function test_parse_header_integrality_both_int()
    model = NL._CacheModel()
    NL._resize_variables(model, 4)
    io = IOBuffer()
    write(
        io,
        """g3 0 1 0	# problem test_simple
2 1 1 0 0	# vars, constraints, objectives, ranges, eqns
1 1	# nonlinear constraints, objectives
0 0	# network constraints: nonlinear, linear
2 1 1	# nonlinear vars in constraints, objectives, both
0 0 0 1	# linear network variables; functions; arith, flags
0 0 1 0 0	# discrete variables: binary, integer, nonlinear (b,c,o)
1 1	# nonzeros in Jacobian, gradients
0 0	# max name lengths: constraints, variables
0 0 0 0 0	# common exprs: b,c,o,c1,o1\n""",
    )
    seekstart(io)
    NL._parse_header(io, model)
    @test model.variable_type[1] == NL._INTEGER
    @test model.variable_type[2] == NL._CONTINUOUS
    return
end

"""
    test_mac_minlp()

This function tests reading the files in Sven Leyffer's MINLP testset, which is
available at: https://wiki.mcs.anl.gov/leyffer/index.php/MacMINLP.

If the folder isn't available locally, the test is skipped.
"""
function test_mac_minlp()
    dir = joinpath(@__DIR__, "MacMINLP")
    if !isdir(dir)
        return
    end
    exclude = [
        # Has subexpressions
        "space-25-r.nl",
        "space-960-ir.nl",
        "space-960-r.nl",
        "top1-15x05.nl",
    ]
    for file in filter(f -> !(f in exclude) && endswith(f, ".nl"), readdir(dir))
        model = NL.Model()
        open(joinpath(dir, file), "r") do io
            return read!(io, model)
        end
    end
    return
end

end

TestNonlinearRead.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
module Uploaders

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export uploader

function __init__()
  Genie.Renderer.Html.register_normal_element("q__uploader", context = Genie.Renderer.Html)
end

function uploader(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  wrap() do
    Genie.Renderer.Html.q__uploader(args...; kwargs...)
  end
end

end


text:
"""
    LD(A::Diagonal)
    LD(A::HBlikDiag)
    LD(A::DenseMatrix)

Return the value of `log(det(triu(A)))` calculated in place.
"""
LD(d::Diagonal) = sum(log, d.diag)

function LD(d::HBlkDiag)
    aa = d.arr
    p, q, k = size(aa)
    pq = p * q
    dd = diagind(p, q)
    r = sum(i -> log(aa[i]), dd)
    for j in 1:(k - 1)
        r += sum(i -> log(aa[i]), dd + j * pq)
    end
    r
end

LD(d::DenseMatrix) = sum(i -> log(d[i]), diagind(d))

"""
    logdet(m::LinearMixedModel)

Return the value of `log(det(Λ'Z'ZΛ + I))` calculated in place.
"""
function Base.logdet{T}(m::LinearMixedModel{T})
    k, R = length(m.Λ), m.R
    res = zero(T)
    for i in 1 : length(m.Λ)
        res += T(LD(R[i, i]))
    end
    2 * res
end


text:
"""
`Basics.modify()`  ... modifies some data, either interactively or due to some given rule.

  + `("level energies: interactive", multiplet::Multiplet)`
    ... to shift the total energies of the atomic levels in multiplet interactively; a (new) multiplet::Multiplet is
        returned in which the total energies are just modified accordingly. **Not yet implemented !**
"""
function Basics.modify(sa::String, multiplet::Multiplet)
    !(sa == "level energies: interactive")   &&   error("Unsupported keystring = $sa")
    error("Not yet implemented !")

    return( nothing )
end


text:
using Distributions
using Plots

function rt_dickey(n, μ, ν, σ)
    # sample y firstly
    dist_y = InverseGamma(ν/2, ν/2)
    ys = rand(dist_y, n)
    # sample x
    xs = σ * randn(n) .* sqrt.(ys) .+ μ
    return xs, ys
end

function cmp_res(n, μ, ν, σ)
    xs, ys = rt_dickey(n, μ, ν, σ)
    cum_δm = cumsum( exp.(- xs .^2) )
    cum_δm_star = cumsum( 1 ./ sqrt.(2 * σ^2 .* ys .+ 1) .* exp.(-μ^2 ./ (1 .+ 2*σ^2 .* ys)) )
    δm = cum_δm ./ (1:n)
    δm_star = cum_δm_star ./ (1:n)
    p = plot(δm, label = "MC")
    plot!(p, δm_star, label = "RB")
    return p
end

using Random
Random.seed!(123)
p1 = cmp_res(10000, 0, 4.6, 1)
p2 = cmp_res(10000, 3, 5, 0.5)
plot(p1, p2)
savefig("two-situations.svg")

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
struct LocalVEMOperators{T}
  order::Int
  D::Matrix{T}
  G::Matrix{T}
  B::Matrix{T}
  H::Matrix{T}
  b::Vector{T}
end

get_order(lvem::LocalVEMOperators) = lvem.order

struct VEMOperators
    dofs::DofHandler
    elements::Vector{LocalVEMOperators}
end

struct CellCache{ct}
  cell::ct
  ci::Int
  tess
  quad
  f
end

function VEMOperators(dof::DofHandler, u::TrialFunction; load = (x->0.0)) where {dim,T}
    elements = get_local_elements(dof, getelement(getfunctionspace(u)), load)
    VEMOperators(dof,elements)
end

function get_local_elements(dof::DofHandler{2}, element::PoissonVirtualElement{2}, load)
    local_elements = Vector{LocalVEMOperators}()
    for (ci, cell) in enumerate(getcells(dof.mesh))
        vertices = getverticescoords(dof.mesh, ci)
        tess = get_2Dtesselation(vertices) #tesselate cell
        quad = QuadratureRule{Triangle}(Strang(),get_degree(element))
        cache = CellCache(cell, ci, tess, quad,load)
        centroid = cell_centroid(dof.mesh, ci)
        diameter = cell_diameter(dof.mesh, ci)
        lelement = LocalPoissonVirtualElement(2,get_degree(element),centroid,diameter)
        D = _compute_local_D(lelement, dof, cache)
        G = _compute_local_G(lelement, dof, cache)
        B = _compute_local_B(lelement, dof, cache)
        H = _compute_local_H(lelement, dof, cache)
        b = _compute_local_b(lelement, dof, cache)
        push!(local_elements, LocalVEMOperators(get_degree(element),D,G,B,H,b))
    end
    local_elements
end

# NOTE: From here the methods work only for 2D
""" Compute local D matrix """
function _compute_local_D(element::LocalPoissonVirtualElement, dof::DofHandler{2}, cache::CellCache{CT}) where {CT}
    cell = cache.cell; ci = cache.ci
    degree = get_degree(element)
    ndofs = ndofs_per_cell(dof, ci)
    nk = Int((degree+1)*(degree+2)/2)
    D = zeros(ndofs,nk)
    nv = getnvertices(dof.mesh,cell)
    vertices = getverticescoords(dof.mesh, ci)
    # compute value in vertices
    for i in 1:nv, j in 1:nk
            D[i,j] = value(element.Pk_basis, j, vertices[i])
    end
    if degree >= 2
        ne = getnedges(dof.mesh,cell)
        # Compute value for edges
        for i in (nv+1):(nv+ne), j in 1:nk
            D[i,j] = value(element.Pk_basis, j,
                map_unit_to_segment(element.edge_quad.points[2],
                        getverticescoords(dof.mesh, EdgeIndex(ci,i-nv))))
        end
        #Compute the Internal dofs
        nf = Int((degree-1)*degree/2)   #dimension of Pk_basis with order k-2
        lap_basis = Monomials(2,degree-2,element.Pk_basis.centroid, element.Pk_basis.diameter)
        for i in (nv+ne*(degree-1)+1):ndofs, j in 1:nk
             for k in cache.tess
               pt = mapPointsFromReference(Triangle,k,cache.quad.points);
               for g in 1:size(pt,1)
                 D[i,j] = D[i,j] + 2*simplex_area(k)*
                   cache.quad.weights[g]*
                   (value(element.Pk_basis, j, pt[g])*
                    value(lap_basis, i-(nv+ne*(degree-1)),
                   pt[g]))/cell_volume(dof.mesh, ci)
               end
             end
           end
         end
    return D
end

function p0m(element::LocalPoissonVirtualElement, dof::DofHandler{2},cache::CellCache{CT}, i::Int) where {CT}
  cell = cache.cell; ci = cache.ci
  degree = get_degree(element)
  nv = getnvertices(dof.mesh,cell)
  P0m = 0.0
   if degree == 1
     vertices = getverticescoords(dof.mesh, ci)
     for j in 1:nv
       P0m = P0m + value(element.Pk_basis, i, vertices[j])
     end
     P0m = P0m/nv
   else
     for k in cache.tess
       pt = mapPointsFromReference(Triangle,k,cache.quad.points);
       for g in 1:size(pt,1)
         P0m = P0m + 2*simplex_area(k)*
           cache.quad.weights[g]*value(element.Pk_basis, i, pt[g])
       end
     end
     P0m = P0m/cell_volume(dof.mesh, ci)
   end
   return P0m
end

""" Compute local G matrix """
function _compute_local_G(element::LocalPoissonVirtualElement, dof::DofHandler{2}, cache::CellCache{CT}) where {CT}
  cell = cache.cell; ci = cache.ci
  degree = get_degree(element)
  nk = Int((degree+1)*(degree+2)/2)
  ne = getnedges(dof.mesh,cell)
  G = zeros(nk, nk );
  # Gradient terms
  # First row
  for j in 1:nk
    G[1,j] = p0m(element, dof, cache, j);
  end
  if degree == 1
     # No need to integrate over the domain in case of k=1
     for i in 2:nk, j in 2:nk
         for k in 1:ne
           v = getverticescoords(dof.mesh, EdgeIndex(ci, k))
           normal = get_Normal(dof.mesh, EdgeIndex(ci, k))
           distance = norm(v[2]-v[1])

           # Use trapezium rule
           grad_m_i_v1 = gradient_value(element.Pk_basis, i, v[1])
           m_j_v1 = value(element.Pk_basis, j, v[1])

           grad_m_i_v2 = gradient_value(element.Pk_basis, i, v[2])
           m_j_v2 = value(element.Pk_basis, j, v[2])

           G[i,j] = G[i,j] + distance/2 * (
             dot(grad_m_i_v1,normal)*m_j_v1 +
             dot(grad_m_i_v2,normal)*m_j_v2 )
         end
     end
   else
     for i in 2:nk, j in 2:nk
         for k in cache.tess
           pt = mapPointsFromReference(Triangle,k,cache.quad.points);
           for g in 1:size(pt,1)
             G[i,j] = G[i,j] + 2*simplex_area(k)*
                cache.quad.weights[g]*
                dot(gradient_value(element.Pk_basis, i, pt[g]),
                    gradient_value(element.Pk_basis, j, pt[g]),)
           end
         end
       end
   end
   return G
end

function compute_d_αβ(element::LocalPoissonVirtualElement, dof::DofHandler{2},cache::CellCache{CT}, k::Int) where {CT}
  cell = cache.cell; ci = cache.ci
  degree = get_degree(element)
  nv = getnvertices(dof.mesh,cell)
  nkm2 = Int((degree-1)*degree/2)
  # Dimension of M_{k-2}: n_{k-2}
  d = zeros(nkm2,1);
  h = cell_diameter(dof.mesh, ci)
  # Build dx^2 monomial
  dxx_monomial_index = element.Pk_basis.indices[k][1]-2
  laplacian_indices = get_monomial2DIndices(degree - 2)
  if (dxx_monomial_index >= 0)
    for i in 1:nkm2
      if (laplacian_indices[i] == (dxx_monomial_index, element.Pk_basis.indices[k][2]))
        d[i] = d[i] + element.Pk_basis.indices[k][1]*(element.Pk_basis.indices[k][1]-1)/h^2;
      end
    end
  end
  # Build dy^2 monomial
  dyy_monomial_index = element.Pk_basis.indices[k][2]-2
  if (dyy_monomial_index >= 0)
    for i in 1:nkm2
      if (laplacian_indices[i] == (element.Pk_basis.indices[k][1],dyy_monomial_index))
        d[i] = d[i] + element.Pk_basis.indices[k][2]*(element.Pk_basis.indices[k][2]-1)/h^2;
      end
    end
  end
  return d
end

""" Compute the local B matrix """
function _compute_local_B(element::LocalPoissonVirtualElement, dof::DofHandler{2}, cache::CellCache{CT}) where {CT}
  cell = cache.cell; ci = cache.ci
  degree = get_degree(element)
  nk = Int((degree+1)*(degree+2)/2)
  ne = getnedges(dof.mesh,cell)
  ndofs = ndofs_per_cell(dof, ci)
  nv = getnvertices(dof.mesh,cell)
  B = zeros(nk, ndofs);

   if (degree == 1)
     # First row
     B[1,:] .= 1/nv
     for i in 2:nk, j in 1:ndofs
         for k in 1:ne
           v = getverticescoords(dof.mesh, EdgeIndex(ci, k))
           normal = get_Normal(dof.mesh, EdgeIndex(ci, k))
           distance = norm(v[2]-v[1])

           # Use trapezium rule
           grad_m_i_v1 = gradient_value(element.Pk_basis, i, v[1])
           φ_j_v1 = δ(j, getverticesindices(dof.mesh, EdgeIndex(ci,k))[1])

           grad_m_i_v2 = gradient_value(element.Pk_basis, i, v[2])
           φ_j_v2 = δ(j, getverticesindices(dof.mesh, EdgeIndex(ci,k))[2])

           B[i,j] = B[i,j] + distance/2 * (
             dot( grad_m_i_v1,normal ) * φ_j_v1 +
             dot( grad_m_i_v2,normal ) * φ_j_v2 )

         end
     end
   else
     # First row
     for j in (nv*degree+1):ndofs
       B[1,j] = p0m(element,dof,cache,j-(nv*degree))
     end
     for i in 2:nk
       d_αβ = compute_d_αβ(element,dof,cache,i)
       for j in 1:ndofs
         # First compute the components owning the Laplacian
         if j > (nv*degree)
           B[i,j] = B[i,j] -
                cell_volume(dof.mesh, ci)*d_αβ[j - (nv*degree)]
         end
         # Now the boundary part
         for k in 1:ne
           v = getverticescoords(dof.mesh, EdgeIndex(ci, k))
           normal = get_Normal(dof.mesh, EdgeIndex(ci, k))
           distance = norm(v[2]-v[1])
           for g=1:size(element.edge_quad.points,1)
             grad_m_i = gradient_value(element.Pk_basis, i,
                  map_unit_to_segment(element.edge_quad.points[g],v))
             edge_dof = -1
             if g == 1
               edge_dof = reference_edge_vertices(dof.mesh,cell)[k][1]
             elseif g == size(element.edge_quad.points,1)
               edge_dof = reference_edge_vertices(dof.mesh,cell)[k][2]
             else
               edge_dof = nv + (degree-1)*(k-1)+(g-1)
             end
             φ_j = δ(j,edge_dof)

             B[i,j] = B[i,j] +
               distance *(element.edge_quad.weights[g]*dot(grad_m_i,normal) * φ_j )
           end
        end
       end
     end
   end
   return B
end

""" Compute the local H matrix """
function _compute_local_H(element::LocalPoissonVirtualElement, dof::DofHandler{2}, cache::CellCache{CT}) where {CT}
  cell = cache.cell; ci = cache.ci
  degree = get_degree(element)
  nk = Int((degree+1)*(degree+2)/2)
  ne = getnedges(dof.mesh,cell)
  ndofs = ndofs_per_cell(dof, ci)
  nv = getnvertices(dof.mesh,cell)
  H = zeros(nk, nk)

  for i in 1:nk, j in 1:nk
      for k in cache.tess
        pt = mapPointsFromReference(Triangle,k,cache.quad.points);
        for g=1:size(pt,1)
          H[i,j] = H[i,j] + 2*simplex_area(k)*
            cache.quad.weights[g]*
            ( value(element.Pk_basis, i, pt[g])* value(element.Pk_basis, j, pt[g]))
        end
      end
  end
  return H
end

""" Compute the local load vector """
function _compute_local_b(element::LocalPoissonVirtualElement, dof::DofHandler{2}, cache::CellCache{CT}) where {CT}
  cell = cache.cell; ci = cache.ci; load_func = cache.f
  degree = get_degree(element)
  if degree == 1
    ndofs = ndofs_per_cell(dof, ci)
    b = load_func(element.Pk_basis.centroid)*ones(ndofs)*cell_volume(dof.mesh, ci)/getnedges(cell)
  else
    #TODO: Should we use k-1?
    nk = Int((degree+1)*(degree+2)/2)
    b = zeros(nk)
    for i in 1:nk
        for k in cache.tess
          pt = mapPointsFromReference(Triangle,k,cache.quad.points);
          for g=1:size(pt,1)
            b[i] = b[i] + 2*simplex_area(k)*
              cache.quad.weights[g]*
              ( value(element.Pk_basis, i, pt[g])* load_func(pt[g]))
          end
        end
    end
  end
  return b
end

function assemble_load(op::VEMOperators)
      b = zeros(ndofs(op.dofs))
      for k = 1:getncells(op.dofs.mesh)
        cell_dofs = Vector{Int}(undef, ndofs_per_cell(op.dofs, k))
        if get_order(op.elements[k]) == 1
          b_local = op.elements[k].b
        elseif get_order(op.elements[k])  < 3
          # In this case: Π^∇ = Π^0
          Π∇s = op.elements[k].G\op.elements[k].B
          b_local = Π∇s'*op.elements[k].b
        else
          throw("The L^2 operator is not yet supported for k>2.")
        end
        # Assembly
        celldofs!(cell_dofs, op.dofs, k)
        assemble!(b,cell_dofs,b_local)
      end
      return b
end

""" Compute stiffness Matrix """
function assemble_stiffnessMat(op::VEMOperators)
    K = create_sparsity_pattern(op.dofs);
    assembler = start_assemble(K)
    for k = 1:getncells(op.dofs.mesh)
      cell_dofs = Vector{Int}(undef, ndofs_per_cell(op.dofs, k))
      # Build local stiffness matrix
      Π∇s = op.elements[k].G\op.elements[k].B
      Π∇ = op.elements[k].D*Π∇s
      nk = size(op.elements[k].G,2)
      G_tilde = [ zeros(1, nk) ;
                  op.elements[k].G[2:end,:] ]
      α = 1;
      K_local = Π∇s'*G_tilde*Π∇s + α*(I-Π∇)'*(I-Π∇)

      # Assembly
      celldofs!(cell_dofs, op.dofs, k)
      assemble!(assembler, cell_dofs, K_local)
    end
    return K
end

""" Compute Mass Matrix """
function assemble_massMat(op::VEMOperators)
  M = create_sparsity_pattern(op.dofs);
  assembler = start_assemble(M)

  for k = 1:getncells(op.dofs.mesh)
    cell_dofs = Vector{Int}(undef, ndofs_per_cell(op.dofs, k))
    if get_order(op.elements[k])  < 3
      # In this case: Π^∇ = Π^0
      Π∇s = op.elements[k].G\op.elements[k].B
      Π∇ = op.elements[k].D*Π∇s

      M_local = Π∇s'*op.elements[k].H*Π∇s
    else
      throw("The L^2 operator is not yet supported for k>2.")
    end
    # Assembly
    celldofs!(cell_dofs, op.dofs, k)
    assemble!(assembler, cell_dofs, M_local)
  end
  return M
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
export get_initial_compat

"""
    get_initial_compat(filename_manifest::S, filename_project::S) where {S<:AbstractString}

A string containing a "first guess" for the `[compat]`
section in a Julia `Project.toml` file, based on existing
packages in the Manifest.toml.
"""
function get_initial_compat(
  filename_manifest::S,
  filename_project::S,
  excludes =
  ["Test",
  "Pkg",
  "InteractiveUtils",
  "Profile"]
  ) where {S<:AbstractString}
  manifest = open(f->read(f, String), filename_manifest)
  project = open(f->read(f, String), filename_project)

  deps = last(split(project, "[deps]"))
  deps = first(split(deps, "[compat]"))
  deps = split(deps, "\n")
  filter!(x->!isempty(x), deps)
  deps = [first(split(x, " = ")) for x in deps]
  filter!(x->!any([x==y for y in excludes]), deps)

  man = split(manifest, "\n\n")
  filter!(x->occursin("[", x), man)
  man = [split(x, "\n") for x in man]
  s = "version ="
  man = Dict([replace(replace(first(x), "[" =>""), "]"=>"") => last(split(last(x),s)) for x in man if occursin(s, last(x))])
  for (k,v) in man
    man[k] = strip(replace(v, "\"" => ""))
  end

  compat = ["$x = \"$(man[x])\"" for x in deps]
  pushfirst!(compat, "[compat]")
  return join(compat, "\n")
end


text:
using DataFrames, Dates

# This contains code which makes up the pipeline of csv -> interaction seqs

export add_page_column!, get_path_sequences, my_filter

# Some intial helper functions
# Take row (vector) and output a label which merges the entries with seperator
# '-', e.g. [1,2,1] -> '1-2-1'

function make_label(X::AbstractVector; nan="NaN")
    return join(X[findall(X.!=nan)], "-")
end

function make_labels(X::AbstractArray; nan="NaN")
    Y = Vector{String}()
    for row in eachrow(X)
        push!(Y, make_label(row))
    end
    return Y
end

# tmp = rand(1:10, 100,10)
# make_labels(tmp)
# tmp1

# Now a function for dataframes

function make_labels(X::AbstractDataFrame, cols::Vector{Int}; nan="NaN")
    X = Matrix(X[:,cols])
    return make_labels(X, nan=nan)
end

function make_labels(X::AbstractDataFrame, cols::Int; nan="NaN")
    return make_labels(X, [cols], nan=nan)
end

function add_page_column!(
    X::AbstractDataFrame,
    cols::Vector{Int};
    level::Int=length(cols), nan="NaN")
    # @show level
    X.page = make_labels(X, cols[1:level], nan=nan)
end

## Splitting into sessions

function get_paths(
    X::AbstractDataFrame;
    T::TimePeriod=Minute(15),
    vertex_col_name::String="vertex",
    time_col_name::String="time_gmt")

    @assert(issorted(X, time_col_name), "Must be sorted according to time.")

    function recurse_split!(
        X::AbstractDataFrame,
        Y::AbstractVector,
        ind::BitArray{1},
        i::Int;
        vertex_col_name::String)

        if findnext(ind, i) == nothing
            push!(Y, X[i:end, vertex_col_name])
            return
        else
            j = findnext(ind, i)
            # println("Path from $(i) to $(j)")
            # println(X[i:j])
            push!(Y, X[i:j, vertex_col_name])
            recurse_split!(X, Y, ind, j+1, vertex_col_name=vertex_col_name)
        end
    end
    Y = Vector{Vector{eltype(X[:,vertex_col_name])}}() # Initialise storage
    ind = diff(X[:,time_col_name]).>T
    recurse_split!(X, Y, ind, 1, vertex_col_name=vertex_col_name)
    return Y
end

## Now do on a grouped data frame (will return a vector of vector of vectors)
"""
Given a `DataFrames.GroupedDataFrame` object and some time threshold of type `Dates.TimePeriod` this will convert the data to interaction sequences.

Column names for vertex labels and timestamps are passed names arguments.

**Note** the time column must be of type `Dates.DateTime`.
"""
function get_path_sequences(
    X::GroupedDataFrame;
    T::TimePeriod=Minute(15),
    vertex_col_name::String="vertex",
    time_col_name::String="time_gmt")

    data = Vector{Vector{Vector{String}}}(undef, length(X))
    for (i, key) in enumerate(keys(X))
        data[i] = get_paths(X[key], T=T, vertex_col_name=vertex_col_name, time_col_name=time_col_name)
    end
    return data
end

## Filtering functions
# On raw dataframe

function max_length(
    X::AbstractDataFrame;
    time_col_name::String="time_gmt",
    T::TimePeriod=Minute(15))

    tmp = diff(X[:,time_col_name]) .> T
    push!(tmp, true)
    pushfirst!(tmp, true)
    return maximum(diff(findall(tmp)))
end

function get_num_paths(
    X::AbstractDataFrame;
    time_col_name::String="time_gmt",
    T::TimePeriod=Minute(15)
    )
    return sum(diff(X[:, time_col_name]) .> T) + 1
end

function my_filter(
    X::GroupedDataFrame;
    T::TimePeriod=Minute(15),
    time_col_name::String="time_gmt",
    max_clicks=7000,
    max_path_length=750,
    max_num_paths=450
    )
    clicks = [size(x)[1] for x in X]
    path_lengths = [max_length(group, time_col_name=time_col_name, T=T) for group in X]
    num_paths = [get_num_paths(group, time_col_name=time_col_name, T=T) for group in X]

    ind = (clicks .< max_clicks) .& (path_lengths .< max_path_length) .& (num_paths .< max_num_paths)

    X[ind]
end

# On session data

function max_length(X::Vector{Vector{T}}) where T
    return maximum(map(length, X))
end

function collapsed_length(X::Vector{Vector{T}}) where T
    z = 0
    for x in X
        z += length(x)
    end
    return z
end

function my_filter(
    X::Vector{Vector{Vector{String}}};
    max_clicks=7000,
    max_path_length=750,
    max_num_paths=450
    )
    ind = (collapsed_length.(X) .< max_clicks) .& (max_length.(X) .< max_path_length) .& (length.(X) .< max_num_paths)
    return X[ind]
end

max_stars_repo_path: src/Data Processing/PathSequences.jl

max_issues_repo_path: src/Data Processing/PathSequences.jl

max_forks_repo_path: src/Data Processing/PathSequences.jl


text:
#### Multiperiod DC Optimal Power Shutoff ####

# This file provides a basic implentation of the multiperiod DC Optimal Power Shutoff
# problem for experimentation on algorithm and problem variations.

# This file can be run by calling `include("mp-dc-ops.jl")` from the Julia REPL

# Developed by Noah Rhodes(@noahrhodes) and Line Roald (@lroald)

###############################################################################
# 0. Initialization
###############################################################################
using Pkg; Pkg.activate(".")
Pkg.instantiate()
Pkg.precompile()

# Load Julia Packages
#--------------------
using PowerModels
using PowerModelsWildfire
using Gurobi
using JuMP
using CSV
using DataFrames

# Load System Data
# ----------------

# load the data file
data = PowerModels.parse_file(joinpath(@__DIR__, "data", "case5_risk_sys1.m"))
load_data = CSV.read(joinpath(@__DIR__, "data", "case5_load.csv"),DataFrame)
risk_data = CSV.read(joinpath(@__DIR__, "data", "case5_risk.csv"),DataFrame)

# Modify risk weighting
# data["risk_weight"] = 0.2

# Add zeros to turn linear objective functions into quadratic ones
# so that additional parameter checks are not required
PowerModels.standardize_cost_terms!(data, order=2)

# Adds reasonable rate_a values to branches without them
PowerModels.calc_thermal_limits!(data)

# Create multiperiod data model
mp_data = PowerModels.replicate(data, size(load_data,1))

for row in eachrow(load_data)
    for col in names(row[2:end])
        col
        nwid = string(row[1])
        comp_type, comp_id = split(col,"_")
        load = row[col]

        mp_data["nw"][nwid][comp_type][comp_id]["pd"] = load
    end
end

for row in eachrow(risk_data)
    for col in names(row[2:end])
        col
        nwid = string(row[1])
        comp_type, comp_id = split(col,"_")
        risk = row[col]

        mp_data["nw"][nwid][comp_type][comp_id]["power_risk"] = risk
    end
end

# use build_ref to filter out inactive components
ref = PowerModels.build_ref(mp_data)
ref_add_on_off_va_bounds!(ref, mp_data)
ref = ref[:it][:pm]

# Note: ref contains all the relevant system parameters needed to build the OPS model
# When we introduce constraints and variable bounds below, we use the parameters in ref.

###############################################################################
# 1. Building the Optimal Power Shutoff Model
###############################################################################

# Initialize a JuMP Optimization Model
#-------------------------------------
# model = Model(GLPK.Optimizer)
model = Model(Gurobi.Optimizer)

# Add Optimization and State Variables
# ------------------------------------

# Add voltage angles va for each bus
@variable(model, va[t in keys(ref[:nw]),i in keys(ref[:nw][1][:bus])])
# note: [i in keys(ref[:nw][1][:bus])] adds one `va` variable for each bus in the network

# Add active power generation variable pg for each generator
@variable(model, pg[t in keys(ref[:nw]),i in keys(ref[:nw][1][:gen])])

# Add power flow variables p to represent the active power flow for each branch
@variable(model,  p[t in keys(ref[:nw]), (l,i,j) in ref[:nw][1][:arcs_from]])

# Build JuMP expressions for the value of p[(l,i,j)] and p[(l,j,i)] on the branches
p_expr = Dict([((t,(l,i,j)), 1.0*p[t,(l,i,j)]) for (l,i,j) in ref[:nw][1][:arcs_from], t in keys(ref[:nw])])
p_expr = merge(p_expr, Dict([((t,(l,j,i)), -1.0*p[t,(l,i,j)]) for (l,i,j) in ref[:nw][1][:arcs_from], t in keys(ref[:nw])]))
# note: this is used to make the definition of nodal power balance simpler

# Add binary variables to represent that state of all components
@variable(model, z_bus[t in keys(ref[:nw]), i in keys(ref[:nw][1][:bus])], Bin)
@variable(model, z_gen[t in keys(ref[:nw]), i in keys(ref[:nw][1][:gen])], Bin)
@variable(model, z_branch[t in keys(ref[:nw]), (l,i,j) in ref[:nw][1][:arcs_from]], Bin)

# Add continous load shed variable
@variable(model, 0.0 <= x_load[t in keys(ref[:nw]), l in keys(ref[:nw][1][:load])] <= 1.0)
@variable(model, 0.0 <= x_shunt[t in keys(ref[:nw]), s in keys(ref[:nw][1][:shunt])] <= 1.0)

# Add Objective Function
# ----------------------
# Maximize power delivery while minimizing wildfire risk

if haskey(ref[:nw][1], :risk_weight)
    alpha = ref[:nw][1][:risk_weight]
else
    @warn "network data should specify risk_weight, using 0.5 as a default"
    alpha = 0.5
end

for comp_type in [:gen, :load, :bus, :branch]
    for (nwid,nw) in ref[:nw]
        for (id,comp) in  nw[comp_type]
            if ~haskey(comp, "power_risk")
                @warn "$(comp_type) $(id) does not have a power_risk value, using 0.1 as a default"
                comp["power_risk"] = 0.1
            end
        end
    end
end

JuMP.@objective(model, Max,
    sum(
        (1-alpha)*(
                sum(x_load[t,i]*load["pd"] for (i,load) in nw[:load])
        )
        - alpha*(
            sum(z_gen[t,i]*gen["power_risk"] for (i,gen) in nw[:gen])
            + sum(z_bus[t,i]*bus["power_risk"] for (i,bus) in nw[:bus])
            + sum(z_branch[t,(l,i,j)]*nw[:branch][l]["power_risk"] for (l,i,j) in nw[:arcs_from])
            + sum(x_load[t,i]*load["power_risk"] for (i,load) in nw[:load])
        )
        for (t,nw) in ref[:nw]
    )
)

# Add Constraints
# ---------------

#  constraints for each period of OPS
for (t,nw) in ref[:nw]

    # Fix the voltage angle to zero at the reference bus
    for (i,bus) in nw[:ref_buses]
        @constraint(model, va[t,i] == 0)
    end

    # Nodal power balance constraints
    for (i,bus) in nw[:bus]
        # Build a list of the loads and shunt elements connected to the bus i
        bus_loads = [(l,nw[:load][l]) for l in nw[:bus_loads][i]]
        bus_shunts = [(s,nw[:shunt][s]) for s in nw[:bus_shunts][i]]

        # Active power balance at node i
        @constraint(model,
            sum(p_expr[(t,a)] for a in nw[:bus_arcs][i]) ==                    # sum of active power flow on lines from bus i +
            sum(pg[t,g] for g in nw[:bus_gens][i]) -                         # sum of active power generation at bus i -
            sum(x_load[t,l]*load["pd"] for (l,load) in bus_loads) -           # sum of active load * load shed  at bus i -
            sum(x_shunt[t,s]*shunt["gs"] for (s,shunt) in bus_shunts)*1.0^2   # sum of active shunt element injections * shed at bus i
        )
    end

    # Branch power flow physics and limit constraints
    for (i,branch) in nw[:branch]
        # Build the from variable id of the i-th branch, which is a tuple given by (branch id, from bus, to bus)
        f_idx = (i, branch["f_bus"], branch["t_bus"])
        p_fr = p[t,f_idx]                     # p_fr is a reference to the optimization variable p[f_idx]
        z_br = z_branch[t,f_idx]
        va_fr = va[t,branch["f_bus"]]         # va_fr is a reference to the optimization variable va on the from side of the branch
        va_to = va[t,branch["t_bus"]]         # va_fr is a reference to the optimization variable va on the to side of the branch
        # Compute the branch parameters and transformer ratios from the data
        g, b = PowerModels.calc_branch_y(branch)

        # Voltage angle difference limit
        JuMP.@constraint(model, va_fr - va_to <= branch["angmax"]*z_br + nw[:off_angmax]*(1-z_br))
        JuMP.@constraint(model, va_fr - va_to >= branch["angmin"]*z_br + nw[:off_angmin]*(1-z_br))

        # DC Power Flow Constraint
        if b <= 0
            JuMP.@constraint(model, p_fr <= -b*(va_fr - va_to + nw[:off_angmax]*(1-z_br)) )
            JuMP.@constraint(model, p_fr >= -b*(va_fr - va_to + nw[:off_angmin]*(1-z_br)) )
        else # account for bound reversal when b is positive
            JuMP.@constraint(model, p_fr >= -b*(va_fr - va_to + nw[:off_angmax]*(1-z_br)) )
            JuMP.@constraint(model, p_fr <= -b*(va_fr - va_to + nw[:off_angmin]*(1-z_br)) )
        end

        # Thermal limit
        JuMP.@constraint(model, p_fr <=  branch["rate_a"]*z_br)
        JuMP.@constraint(model, p_fr >= -branch["rate_a"]*z_br)

        # Connectivity constraint
        JuMP.@constraint(model, z_br <= z_bus[t,branch["f_bus"]])
        JuMP.@constraint(model, z_br <= z_bus[t,branch["t_bus"]])
    end

    # Generator constraints
    for (i,gen) in nw[:gen]
        # Power limit
        JuMP.@constraint(model, pg[t,i] <= gen["pmax"]*z_gen[t,i])
        JuMP.@constraint(model, pg[t,i] >= gen["pmin"]*z_gen[t,i])

        # Connectivity constraint
        JuMP.@constraint(model, z_gen[t,i] <= z_bus[t,gen["gen_bus"]])
    end

    # Load constraints
    for (i, load) in nw[:load]
        # Connectivity constraints
        JuMP.@constraint(model, x_load[t,i] <= z_bus[t,load["load_bus"]])
    end
end

# Inter-period constraints
for t in keys(ref[:nw])
    if t != 1
        # if component is disabled, it must stay disabled for future periods
        for (i,branch) in ref[:nw][t][:branch]
            f_idx = (i, branch["f_bus"], branch["t_bus"])
            JuMP.@constraint(model, z_branch[t-1,f_idx] >= z_branch[t,f_idx])
        end

        for (i,bus) in ref[:nw][t][:bus]
            JuMP.@constraint(model, z_bus[t-1,i] >= z_bus[t,i])
        end
        for (i,gen) in ref[:nw][t][:gen]
            JuMP.@constraint(model, z_gen[t-1,i] >= z_gen[t,i])
        end
    end
end

###############################################################################
# 3. Solve the Optimal Power Flow Model and Review the Results
###############################################################################

# Solve the optimization problem
optimize!(model)

# Check that the solver terminated without an error
println("The solver termination status is $(termination_status(model))")

# Check the value of the objective function
println("The objective value is $(objective_value(model))")
total_load = sum(sum(value(x_load[t,i])*load["pd"] for (i,load) in nw[:load]) for (t,nw) in ref[:nw])
println("The load served is $total_load ")
total_risk = sum(
    sum(value(z_branch[t,(l,i,j)])*nw[:branch][l]["power_risk"] for (l,i,j) in nw[:arcs_from]) +
    sum(value(z_gen[t,g])*gen["power_risk"] for (g,gen) in nw[:gen]) +
    sum(value(z_bus[t,i])*bus["power_risk"] for (i,bus) in nw[:bus]) +
    sum(value(x_load[t,i])*load["power_risk"] for (i,load) in nw[:load])
    for (t,nw) in ref[:nw]
)
println("The system risk is $total_risk")

# Get inactive devices
for t in sort(collect(keys(ref[:nw])))
    nw = ref[:nw][t]
    for l in sort(collect(keys(nw[:branch])))
        f_idx = (l, nw[:branch][l]["f_bus"], nw[:branch][l]["t_bus"])
        value(z_branch[t,f_idx])==0 ? println("Branch $l is disabled in period $t.") : false
    end
    for g in sort(collect(keys(nw[:gen])))
        value(z_gen[t,g])==0 ? println("Generator $g is disabled in period $t.") : false
    end
    for i in sort(collect(keys(nw[:bus])))
        value(z_bus[t,i])==0 ? println("Bus $i is disabled in period $t.") : false
    end
    for i in sort(collect(keys(nw[:load])))
        value(x_load[t,i])!=1 ? println("Load $i is serving $(value(x_load[t,i])*100)% of demand in period $t.") : false
    end
end

max_stars_repo_name: WISPO-POP/PowerModelsWildfireAnnex

max_issues_repo_name: WISPO-POP/PowerModelsWildfireAnnex

max_forks_repo_name: WISPO-POP/PowerModelsWildfireAnnex


text:
using EngEconomics

# Given
i = 0.08
N = 20
g = 0.05
A1 = 1000
A_13_20 = 2000
m = 12

# Find.
# How much should be invested now so that the interest pays after all this can
# pay off everything?

# First Find ic
ic = i / m
Nc = 20

# First figure out how much in total this guy owes in present money
owe = (A1 * presentWorthConversionFactorGeometric(g, ic, 11) + A_13_20 * seriesPresentAmountFactor(ic, 8) * presentWorthFactor(ic, 11))
payNow = owe * presentWorthFactor(ic, 1)


text:
##################### Interfaces for Lattices #####################
"""
    AbstractLattice{D, Z}

Abstract lattice type, where
    * D: the dimension.
    * Z: the coordination number.
"""
abstract type AbstractLattice{D, Z} end

"""
    nsite(c) -> Int

Returns number of site.
"""
function nsite end

"""
    neighbors(c::AbstractLattice, i::Int) -> Vector

Returns indices of neighbor sites.
"""
function neighbors end

"""
    nunitcells(c::AbstractLattice) -> Int

Returns the number of sites in a unit cell.
"""
function nunitcells end

##################### Interfaces for Bins #####################
"""
    AbstractBin

Abstract bin for storing data, with support to `push!`, `length`, `fitted` and `print_lastbin`.
"""
abstract type AbstractBin end

"""
    fitted(bin::AbstractBin) -> Bool

Return true if the data size just fit integer numbers of bins.
"""
function fitted(bin::AbstractBin) end

"""
    print_lastbin(io::IO, bin::AbstractBin)

Show means of data in last bin.
"""
function print_lastbin(io::IO, bin::AbstractBin) end

################## Interfaces for Models #################
"""
    AbstractModel{LT<:AbstractLattice}

Physical model that defines physical parameters.
"""
abstract type AbstractModel{LT<:AbstractLattice} end


text:
@inline function assert_package_found(package_was_found::Bool, name, uuid)
    if !package_was_found
        if uuid === nothing
            msg = "Package with name `$(name)` was not found."
        else
            msg = "Package with name `$(name)` and UUID `$(uuid)` was not found."
        end
        throw(PackageNotFoundException(msg))
    end
    return nothing
end

@inline function assert_no_ambiguity(ambiguity_exists::Bool, name, uuid)
    if ambiguity_exists
        if uuid === nothing
            msg = "Package with name `$(name)` is ambiguous. Disambiguate by providing the UUID."
        else
            msg = "There are multiple different packages with name `$(name)` and UUID `$(uuid)`."
        end
        throw(PackageAmbiguousException(msg))
    end
    return nothing
end

@inline function assert_no_disagreement(a::VersionAndTree, b::VersionAndTree)
    if a != b
        msg = "Two registries disagree. One registry has `$(a)`, and another registry has `$(b)`."
        throw(RegistriesDisagreeException(msg))
    end
    return nothing
end


text:
module BangBang

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) BangBang

export @!,
       @set!!,
       Empty,
       add!!,
       append!!,
       broadcast!!,
       collector,
       delete!!,
       deleteat!!,
       empty!!,
       finish!,
       lmul!!,
       materialize!!,
       merge!!,
       mergewith!!,
       mul!!,
       pop!!,
       popfirst!!,
       prefermutation,
       push!!,
       pushfirst!!,
       rmul!!,
       setdiff!!,
       setindex!!,
       setproperties!!,
       setproperty!!,
       singletonof,
       splice!!,
       union!!,
       unique!!

using Base.Broadcast:
    Broadcasted,
    broadcasted,
    combine_eltypes,
    copyto_nonleaf!,
    instantiate,
    materialize!,
    preprocess
using Base: HasEltype, IteratorEltype, promote_typejoin
using Compat: hasproperty
using ConstructionBase: constructorof
using InitialValues
using LinearAlgebra
using Requires
using Tables: Tables
using ZygoteRules: @adjoint

include("utils.jl")

# Used in NoBang:
function implements end
function push!! end
function unique!! end

include("NoBang/NoBang.jl")
using .NoBang: Empty, SingletonVector, singletonof

include("core.jl")
include("base.jl")
include("linearalgebra.jl")
include("extras.jl")
include("broadcast.jl")
include("collectors.jl")
include("initials.jl")
include("macro.jl")
include("dataframes_impl.jl")
include("zygote.jl")

include("setfield.jl")
using .SetfieldImpl: @set!!, prefermutation

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin
        include("staticarrays.jl")
    end
    @require StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" begin
        include("structarrays.jl")
    end
    @require TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" begin
        include("typedtables.jl")
    end
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("dataframes.jl")
    end
end

end # module


text:
include("analyse_closed_order.jl")
include("check_for_new_closed.jl")
include("check_grid_buy.jl")
include("check_grid_sell.jl")
include("create_trade_grid.jl")
include("finish_grid_bot.jl")
include("freeze_orders.jl")
include("get_buy_value.jl")
include("get_current_volume.jl")
include("get_grid_bots.jl")
include("get_next_gridpoints.jl")
include("init_grid_bot.jl")
include("place_grid_order.jl")
include("place_init_grid_order.jl")
include("update_status.jl")


text:
module Utils

export stack, ambiguity

"""
    check_length_consistency(inputs::AbstractVector...)

Raise a `DimensionMismatch` exception if the input vectors are not all the same
length. Does nothing otherwise.

# Examples
```jldoctest
julia> check_length_consistency([1,2,3], [4,5])
ERROR: DimensionMismatch("Input vectors should have the same length!")
[...]
```
"""
function check_length_consistency(inputs::AbstractVector...)
    if any(!=(length(first(inputs))), map(length, inputs))
        throw(DimensionMismatch("Input vectors should have the same length!"))
    end
end

"""
    check_length_consistency(inputs::AbstractMatrix...)

Raise a `DimensionMismatch` exception if the input matrices don't have the same
number of columns. Does nothing otherwise.

# Examples
```jldoctest
julia> check_length_consistency([1 2 3; 4 5 6], [7 8; 9 10])
ERROR: DimensionMismatch("Input vectors should have the same length!")
[...]
```
"""
function check_length_consistency(inputs::AbstractMatrix...)
    if any(!=(size(first(inputs))[2]), map(i -> size(i)[2], inputs))
        throw(DimensionMismatch("input matrices should have the same number of columns"))
    end
end

# TODO: Consider renaming it to stack_rows
"""
    stack(inputs::AbstractVector...)

Stack the input vectors as rows of a matrix.

# Examples
```jldoctest
julia> stack([1,2,3], [4,5,6])
2×3 Array{Int64,2}:
 1  2  3
 4  5  6
```
"""
function stack(inputs::AbstractVector...)
    check_length_consistency(inputs...)

    mapreduce(permutedims, vcat, inputs)
end

"""
    stack(inputs::AbstractMatrix...)

Stack the input matrices vertically.
"""
function stack(inputs::AbstractMatrix...)
    check_length_consistency(inputs...)

    reduce(vcat, inputs)
end

# TODO: Untested
function accuracy(y_predict::AbstractVector{T}, y_actual::AbstractVector{T}) where {T}
    check_length_consistency(y_predict, y_actual)

    sum(y_predict .== y_actual) / length(y_predict)
end

"""
    hasdraw(X::Tuple)

Return `true` if there is draw in the tuple `X` and `false` otherwise. A draw
occurs when the largest value in the tuple is non-unique.
"""
function hasdraw(X::Tuple)
    length(X) == 0 && throw(ArgumentError("X must have at least one element."))

    repeats = 0
    largest = first(X)

    for x in X
        if x > largest
            largest = x
            repeats = 1
        elseif x == largest
            repeats += 1
        end
    end

    return repeats != 1
end

# TODO: Add example to docstring
"""
    ambiguity(responses...)

Measure the ambiguity in a multidiscriminator classifier prediction.

The ambiguity is the average number of ties in the responses of the
discriminators that compose a multidiscriminator classifier when predicting the
classes of a test set.
"""
function ambiguity(responses::AbstractVector{T}...) where {T <: Integer}
    check_length_consistency(responses...)
    count = 0

    for res in Iterators.zip(responses...)
        if hasdraw(res)
            count += 1
        end
    end

    return count / (length ∘ first)(responses)
end

end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# abstract types
mutable struct MarkovProcess <: AbstractMarkovProcess end

# properties
const markov_process_props = [:renewal, :markovian, :dynamic]

# parameters
function set_parameters!(std::AbstractSTD, cls::AbstractMarkovProcess)
    set_rates!(std, cls)
end

# stochastic process
function homogeneous_markov_process(du, u, p, t)
    G, ρ, info = p
    for ns in 1:_LG.nv(G)
        if info[ns].trapping
            du[ns] = sum(ρ[_LG.Edge(nt,ns)]*u[nt]
                            for nt in _LG.inneighbors(G,ns)
                            if ns ≠ nt)
        else
            du[ns] = sum(ρ[_LG.Edge(nt,ns)]*u[nt]
                            for nt in _LG.inneighbors(G,ns)
                            if ns ≠ nt) -
                     sum(ρ[_LG.Edge(ns,nt)]*u[ns]
                            for nt in _LG.outneighbors(G,ns)
                            if ns ≠ nt)
    end end
end
function inhomogeneous_markov_process(du, u, p, t)
    G, ρ, info = p
    for ns in 1:_LG.nv(G)
        if info[ns].trapping
            du[ns] = sum(ρ[_LG.Edge(nt,ns)](t)*u[nt]
                            for nt in _LG.inneighbors(G,ns)
                            if ns ≠ nt)
        else
            du[ns] = sum(ρ[_LG.Edge(nt,ns)](t)*u[nt]
                            for nt in _LG.inneighbors(G,ns)
                            if ns ≠ nt) -
                     sum(ρ[_LG.Edge(ns,nt)](t)*u[ns]
                            for nt in _LG.outneighbors(G,ns)
                            if ns ≠ nt)
    end end
end
function solve!(std::AbstractSTD, cls::AbstractMarkovProcess; tsim::Number=1.0u"yr", tol::Real=1e-8)
    set_parameters!(std, cls)
    p   = [std.graph,get_tprop(std,:rate),get_sprop(std,:info)]
    u₀  = get_sprop(std,:init)
    ts  = (0.0_UF.unit(tsim),tsim)

    if get_info(std,:time_homogeneous)
        prob = _ODE.ODEProblem(homogeneous_markov_process, u₀, ts, p)
    else
        prob = _ODE.ODEProblem(inhomogeneous_markov_process, u₀, ts, p)
    end
    sol = _ODE.solve(prob, _ODE.Tsit5(), reltol = tol, abstol = tol)

    set_prop!(std, :cls, cls)
    set_prop!(std, :time, sol.t)
    set_prop!(std, states(std), :prob, [sol[ns,:] for ns in states(std)])

    set_info!(std, :solved, true)
end


text:
"""
Solve the constrained least-squares problem

    minimize  ||Ax-b||₂  st  0 ≤ x ≤ u.

Rephrases the problem as linear SOC program:

    minimize  [0]'[x]  st  [A 0 I][x] = [b] m
              [0] [r]             [t]
                                  [r]
                            n 1 m

                               x, (t,r) ∈  Rn₊ x (Qn₊)

with dimensions  len(x) = n  len(b) = m
                 len(r) = m

Solve using ECOS, which formulates linear SOC programs as

    minimize  c'x  st  A*x=b, G*x + s = h, s in cone
"""
function lsbox(Ain::Matrix{Float64}, b::Vector{Float64},
               u::Vector{Float64}=zeros(0); kwargs...)
    m, n = size(Ain)
    nu = length(u)       # number of vars x with upper bound
    n_ecos = n + 1 + m   # number of vars == len(x)
    p = m                # number of equality constraints == len(b)
    l = n + nu           # number of vars in positive orthant
    ncones = 1           # number of SOC cones
    q = [1+m]            # number of vars in the SOC
    e = 0                # number of exponential cones

    GG = [ sparse(-1.0I,n ,n)  spzeros(n ,1)       spzeros(n ,m)
           sparse(+1.0I,nu,n)  spzeros(nu,1)       spzeros(nu,m)
                spzeros(1 ,n)      -1.0            spzeros( 1,m)
                spzeros(m ,n)  spzeros(m ,1)  sparse(-1.0I,m ,m)]
    G = ECOS.ECOSMatrix(GG)
    h = vcat(zeros(n), u,   0,  zeros(m))
    m_ecos =       n + nu + 1 + m
    A = ECOS.ECOSMatrix([Ain spzeros(m,1) sparse(1.0I,m,m)])
    c = zeros(n_ecos); c[n+1] = 1.0

    ecos_ptr = ECOS.setup(n_ecos, m_ecos, p, l, ncones, q, e,
                          G, A, c, h, b; kwargs...)
    ECOS.solve(ecos_ptr)
    # finalizer(cleanup!, problem)

    # Extract the primal and dual solutions of the ECOS problem
    ecos_prob = unsafe_wrap(Array, ecos_ptr, 1)[1]
    xtr = unsafe_wrap(Array, ecos_prob.x, n_ecos)[:]

    ECOS.cleanup(ecos_ptr, 0)
    x = xtr[1:n]
    t = xtr[n+1] # t = norm(r)
    r = xtr[n+2:end]
    return x, r, t
end