text:
function test0(args...)
    @testset "basic test" begin
        m_slave = ModelWithParams(args...)

        x = add_parameters(m_slave, 4.0*ones(2))
        @test x == all_parameters(m_slave)
        @variable(m_slave, y[1:6])

        @constraint(m_slave, ctr1, 3*y[1] >= 2 - 7*x[1])

        @objective(m_slave, Min, 5*y[1])

        JuMP.optimize!(m_slave)

        @test 5/3 ≈ JuMP.dual(ctr1) atol=1e-3
        @test [-35/3, 0.0] ≈ JuMP.dual.(x) atol=1e-3
        @test [-26/3, 0.0, 0.0, 0.0, 0.0, 0.0] ≈ JuMP.value.(y) atol=1e-3
        @test -130/3 ≈ JuMP.objective_value(m_slave) atol=1e-3
        @test parametrized_dual_objective_value(m_slave) ≈ -35/3 * x[1] + 10/3 atol=1e-3
    end
end

function test1(args...)
    @testset "multiple parameters" begin
        m_slave = ModelWithParams(args...)

        x = add_parameters(m_slave, 4.0*ones(6))
        @variable(m_slave, y[1:6])

        @constraint(m_slave, ctr1, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr2, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr3, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr4, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr5, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr6, 3*y[1] >= 2 - 7*x[3])
        @constraint(m_slave, ctr7, sum(3*y[i]+x[i] for i in 2:4) >= 2 - 7*x[3])
        @constraint(m_slave, ctr8, sum(3*y[i]+7.0*x[i]-x[i] for i in 2:4) >= 2 - 7*x[3])

        @objective(m_slave, Min, 5*y[1])

        JuMP.optimize!(m_slave)

        @test 5/3 ≈ JuMP.dual(ctr1) + JuMP.dual(ctr2) + JuMP.dual(ctr3) + JuMP.dual(ctr4) + JuMP.dual(ctr5) + JuMP.dual(ctr6) atol=1e-3
        @test 0.0 ≈ JuMP.dual(ctr7) atol=1e-3
        @test 0.0 ≈ JuMP.dual(ctr8) atol=1e-3
        @test [0.0, 0.0, -35/3, 0.0, 0.0, 0.0] ≈ JuMP.JuMP.dual.(x) atol=1e-3
        @test [-26/3, 0.0, 0.0, 0.0, 0.0, 0.0] ≈ JuMP.JuMP.value.(y) atol=1e-3
        @test -130/3 ≈ JuMP.objective_value(m_slave) atol=1e-3
        @test parametrized_dual_objective_value(m_slave) ≈ -35/3 * x[3] + 10/3 atol=1e-3
    end
end

function test2(args...)
    @testset "LessThan modification" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x ≤ α)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α - 2 atol=1e-3

        fix(α, 2.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 2.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α + 4 atol=1e-3
    end
end

function test3(args...)
    @testset "EqualTO modification" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x == α)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α - 2 atol=1e-3

        fix(α, 2.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 2.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α + 4 atol=1e-3
    end
end

function test4(args...)
    @testset "GreaterThan modification" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x >= α)
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ 1α

        fix(α, 2.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 2.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ 1α
    end
end

function test5(args...)
    @testset "Multiple parameters" begin
        model = ModelWithParams(args...)
        α = add_parameters(model, ones(10))
        @variable(model, x)
        cref = @constraint(model, x == sum(2 * α[i] for i in 1:10))
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 20.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α[3]) == 2.0
        @test parametrized_dual_objective_value(model) ≈ 2sum(α)
    end
end

function test6(args...)
    @testset "Mixing parameters and vars 1" begin
        model = ModelWithParams(args...)
        α = add_parameters(model, ones(5))
        @variable(model, x)
        cref = @constraint(model, sum(x for i in 1:5) == sum(2 * α[i] for i in 1:5))
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 2.0
        @test JuMP.dual(cref) == 1/5
        @test JuMP.dual(α[3]) == 2/5
        @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α)
    end
end

function test7(args...)
    @testset "Mixing parameters and vars 2" begin
        model = ModelWithParams(args...)
        α = add_parameters(model, ones(5))
        @variable(model, x)
        cref = @constraint(model, 0.0 == sum(-x + 2 * α[i] for i in 1:5))
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 2.0
        @test JuMP.dual(cref) == 1/5
        @test JuMP.dual(α[3]) == 2/5
        @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α)
    end
end

function test8(args...)
    @testset "Mixing parameters and vars 3" begin
        model = ModelWithParams(args...)
        α = add_parameters(model, ones(5))
        @variable(model, x)
        cref = @constraint(model, 0.0 == sum(-x + 2.0 + 2 * α[i] for i in 1:5))
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 4.0
        @test JuMP.dual(cref) == 1/5
        @test JuMP.dual(α[3]) == 2/5
        @test parametrized_dual_objective_value(model) ≈ 2/5 * sum(α) + 2
    end
end

function test9(args...)
    @testset "Test bad model error" begin
        model_1 = Model(args...)
        @test_throws ErrorException x = add_parameters(model_1, ones(5))
        @test_throws ErrorException y = add_parameter(model_1, 1.0)
    end
    @testset "Test lazy duals errors" begin
        model_2 = ModelWithParams(args...)
        ParameterJuMP.set_lazy_duals(model_2)
        ParameterJuMP.set_lazy_duals(model_2) # warn
        ParameterJuMP.set_not_lazy_duals(model_2)
        ParameterJuMP.set_not_lazy_duals(model_2) # warn
        y = add_parameter(model_2, 1.0)
        @test_throws ErrorException ParameterJuMP.set_lazy_duals(model_2)

        model_3 = ModelWithParams(args...)
        ParameterJuMP.set_lazy_duals(model_3)
        y = add_parameter(model_3, 1.0)
        @test_throws ErrorException ParameterJuMP.set_not_lazy_duals(model_3)
        @test !ParameterJuMP.is_sync(model_3)
    end
end

function test10(args...)
    @testset "Add after solve" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x <= α)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α - 2

        b = add_parameter(model, -2.0)
        @test all_parameters(model) == [α, b]
        cref = @constraint(model, x <= b)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -2.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 0.0
        @test JuMP.dual(b) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -b - 4
    end
end

function test11(args...)
    @testset "Set coefficient" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        @variable(model, x)
        cref = @constraint(model, x <= 0.0)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 0.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 0.0
        @test parametrized_dual_objective_value(model) ≈ 0α

        set_coefficient(cref, α, 1.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ α - 2
    end
    @testset "Change coefficient" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x >= α)
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ 1α

        ParameterJuMP.set_coefficient(cref, α, -2.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -2.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 2.0
        @test parametrized_dual_objective_value(model) ≈ 2α
    end
    @testset "Set coefficient with lazy" begin
        model = ModelWithParams(args...)
        ParameterJuMP.set_lazy_duals(model)
        α = add_parameter(model, 1.0)
        @variable(model, x)
        cref = @constraint(model, x <= 0.0)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 0.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 0.0
        @test parametrized_dual_objective_value(model) ≈ 0α

        set_coefficient(cref, α, 1.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ α - 2
    end
    @testset "Set coefficient with lazy 2" begin
        model = ModelWithParams(args...)
        ParameterJuMP.set_lazy_duals(model)
        α = add_parameter(model, 1.0)
        b = add_parameter(model, 0.0)
        @variable(model, x)
        cref = @constraint(model, x <= b)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 0.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 0.0
        @test JuMP.dual(b) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -b

        set_coefficient(cref, α, 1.0)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == 1.0
        @test JuMP.dual(b) == -1.0
        @test parametrized_dual_objective_value(model) ≈ α - b - 2
    end
end

function test12(args...)
    @testset "Remove parameter from constraint" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x >= α)
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ 1α

        ParameterJuMP.delete_from_constraint(cref, α)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 0.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 0.0
        @test parametrized_dual_objective_value(model) ≈ 0α
    end
    @testset "Remove parameter from all constraint" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x >= α)
        @objective(model, Min, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 1.0
        @test parametrized_dual_objective_value(model) ≈ 1α

        ParameterJuMP.delete_from_constraints(α)
        JuMP.optimize!(model)
        @test JuMP.value(x) == 0.0
        @test JuMP.dual(cref) == 1.0
        @test JuMP.dual(α) == 0.0
        @test parametrized_dual_objective_value(model) ≈ 0α
    end
end

function test13(args...)
    @testset "Test no duals errors" begin
        model = ModelWithParams(args...)
        ParameterJuMP.set_no_duals(model)
        α = add_parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        cref = @constraint(model, x == α)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test isnan(JuMP.dual(α))

        model_2 = ModelWithParams(args...)
        y = add_parameter(model_2, 1.0)
        @test_throws ErrorException ParameterJuMP.set_no_duals(model_2)
    end
end

using SparseArrays

# macros from jump
macro test_expression(expr)
    esc(quote
            @test JuMP.isequal_canonical(@expression(m, $expr), $expr)
    end)
end
macro test_expression_with_string(expr, str)
    esc(quote
            @test string(@inferred $expr) == $str
            @test_expression $expr
    end)
end

function test14(args...)
    @testset "Test ParametrizedAffExpr" begin
        m = ModelWithParams()
        @variable(m, x)
        @variable(m, y)
        a = add_parameter(m)
        @test name(a) == ""
        set_name(a, "a")
        @test name(a) == "a"
        b = add_parameter(m)
        set_name(b, "b")
        @test all_parameters(m) == [a, b]

        exp1 = x + y + a
        @test typeof(exp1) == ParameterJuMP.PAE{Float64}
        @test length(exp1.v.terms) == 2
        exp1 = exp1 + y
        @test length(exp1.v.terms) == 2

        @test iszero(ParameterJuMP.PAE{Float64}())
        @test iszero(zero(exp1))
        @test iszero(one(exp1) - one(ParameterJuMP.PAE{Float64}))
        @test iszero(SparseArrays.dropzeros((exp1 - copy(exp1))))

        empty_func(empty_arg) = 0
        exp2 = map_coefficients(Base.zero, exp1)
        @test iszero(SparseArrays.dropzeros(exp2))
        @test iszero(constant(exp2))

        @test isequal_canonical(exp1, copy(exp1))

        exp4 = exp1 - copy(exp1)
        @test iszero(SparseArrays.dropzeros(exp4))

        # var + num
        @test_expression_with_string 7.1 * x + 2.5 "7.1 x + 2.5"
        @test_expression_with_string 1.2 * y + 1.2 "1.2 y + 1.2"
        # par + num
        @test_expression_with_string 7.1 * a + 2.5 "7.1 a + 2.5"
        @test_expression_with_string 1.2 * b + 1.2 "1.2 b + 1.2"

        # par + par + num
        @test_expression_with_string b + a + 1.2 "b + a + 1.2"
        # par - par + num
        @test_expression_with_string b - a + 1.2 "b - a + 1.2"
        # par - (par - num)
        @test_expression_with_string b - (a - 1.2) "b - a + 1.2"
        # par - par - num
        @test_expression_with_string b - a - 1.2 "b - a - 1.2"
        # var + par + num
        @test_expression_with_string x + a + 1.2 "x + a + 1.2"
        # var + var + par + num
        @test_expression_with_string x + y + a + 1.2 "x + y + a + 1.2"
        # var + var - par + num
        @test_expression_with_string x + y - a + 1.2 "x + y - a + 1.2"
        # var + var - par + par + num
        @test_expression_with_string x + y - a + b + 1.2 "x + y + b - a + 1.2"
        # par + par - var + var + num
        @test_expression_with_string a + b - x + y + 1.2 "y - x + a + b + 1.2"

        exp5 = x + y
        @test_expression_with_string x + y "x + y"
        @test_expression_with_string convert(ParameterJuMP.PAE{Float64}, exp5) "x + y"

    end
end

function test15(args...)
    @testset "Test add ctr - alternative" begin
        model = ModelWithParams(args...)
        α = Parameter(model, 1.0)
        fix(α, -1.0)
        @variable(model, x)
        exp = x - α
        cref = @constraint(model, exp ≤ 0)
        @objective(model, Max, x)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α - 2
    end
end

function test15(args...)
    @testset "Test macro" begin
        m = ModelWithParams()
        @variable(m, x)
        @variable(m, y)
        @variable(m, a, ParameterJuMP.Param())
        @variable(m, b, ParameterJuMP.Param())
        # TODO - escaping
        # @variable(m, b, Param())
        # TODO - anonimous
        # b = @variable(m, Param())
        # set_name(b, "b")
        @test_expression_with_string a + b - x + y + 1.2 "y - x + a + b + 1.2"
    end
end

function test16(args...)
    @testset "Test deletion of constraint" begin
        model = ModelWithParams(args...)
        α = add_parameter(model, -1.0)
        @variable(model, x)
        cref1 = @constraint(model, x ≤ α/2)
        cref2 = @constraint(model, x ≤ α)
        cref3 = @constraint(model, x ≤ 2α)
        @objective(model, Max, x)
        JuMP.delete(model, cref3)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -1.0
        @test JuMP.dual(cref1) == 0.0
        @test JuMP.dual(cref2) == -1.0
        @test JuMP.dual(α) == -1.0
        @test parametrized_dual_objective_value(model) ≈ -α - 2

        JuMP.delete(model, cref2)
        JuMP.optimize!(model)
        @test JuMP.value(x) == -0.5
        @test JuMP.dual(cref1) == -1.0
        @test JuMP.dual(α) == -0.5
        @test parametrized_dual_objective_value(model) ≈ -0.5α - 1
    end
end


text:
#=
test-io.jl
test AVSfldIO. through FileIO
=#

using FileIO: File, @format_str
using Test: @test, @testset, @test_throws
import AVSfldIO

@testset "loadsave" begin
     data1 = reshape(Int32.(1:20), 4, 5)

     file = tempname() * ".fld"

     ff = File{format"AVSfld"}(file)
     AVSfldIO.save(ff, data1)
     data2 = AVSfldIO.load(ff)
     @test data2 == data1

     rm(file, force=true)
end

@testset "save type" begin
    file = tempname() * ".fld"
    ff = File{format"AVSfld"}(file)
    @test_throws ArgumentError AVSfldIO.save(ff, "string")
end

ext: jl

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO


text:
using DashBootstrapComponents, DashHtmlComponents

buttons = html_div([
    dbc_button("Primary", color="primary", className="mr-1"),
    dbc_button("Secondary", color="secondary", className="mr-1"),
    dbc_button("Success", color="success", className="mr-1"),
    dbc_button("Warning", color="warning", className="mr-1"),
    dbc_button("Danger", color="danger", className="mr-1"),
    dbc_button("Info", color="info", className="mr-1"),
    dbc_button("Light", color="light", className="mr-1"),
    dbc_button("Dark", color="dark", className="mr-1"),
    dbc_button("Link", color="link"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""

# Download stock market data from Yahoo Finance

    get_symbols(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols(symbol::String, from::String, to::String)

    from = string(Integer(datetime2unix(DateTime(from * "T12:00:00"))))
    to = string(Integer(datetime2unix(DateTime(to * "T12:00:00"))))
    parse(Int, from) > parse(Int, to) ? error("in get_symbols: to must be older than from") : nothing

    url = "https://query1.finance.yahoo.com/v7/finance/chart/$symbol?&interval=1d&period1=$from&period2=$to"

    response = HTTP.get(url, cookies = true)
    body = JSON.parse(String(response.body))["chart"]["result"][1]
    values = body["indicators"]["quote"][1]

    x = DataFrame(
        Open = values["open"],
        High =  values["high"],
        Low = values["low"],
        Close = values["close"],
        Volume = values["volume"],
        Adjusted = body["indicators"]["adjclose"][1]["adjclose"],
        Time = Dates.Date.(unix2datetime.(body["timestamp"]))
    )

    deleterows!(x, isnothing.(x).Close)
    x = TimeArray(x, timestamp = :Time)

    return(x)
end


text:
include("EllCrv/EllCrv.jl")
include("EllCrv/Misc.jl")
include("EllCrv/QQ.jl")
include("EllCrv/Finite.jl")


text:
"""
MurmurHash3 was written by Austin Appleby, and is placed in the public
domain. The author hereby disclaims copyright to this source code.

This version was translated into Julia by Scott P. Jones
It is licensed under the MIT license

Note - The x86 and x64 versions do _not_ produce the same results, as the
algorithms are optimized for their respective platforms. You can still
compile and run any of them on any platform, but your performance with the
non-native version will be less than optimal.
"""
module MurmurHash3
export mmhash128_a, mmhash128_u, mmhash128_c, mmhash32

u8(val)   = val%UInt8
u32(val)  = val%UInt32
u64(val)  = val%UInt64
u128(val) = val%UInt128

@inline rotl(x::Unsigned, r) = (x << r) | (x >>> (sizeof(typeof(x))*8 - r))

@inline xor33(k::UInt64) = xor(k, k >>> 33)

@inline rotl27(k) = rotl(k, 27)
@inline rotl31(k) = rotl(k, 31)
@inline rotl33(k) = rotl(k, 33)

#-----------------------------------------------------------------------------
# Finalization mix - force all bits of a hash block to avalanche

export dbf
const dbf = Ref(false)

@inline fmix(k::UInt64) = xor33(xor33(xor33(k) * 0xff51afd7ed558ccd) * 0xc4ceb9fe1a85ec53)

const c1 = 0x87c37b91114253d5
const c2 = 0x4cf5ad432745937f

@inline mhtail1(h1, k1) = xor(h1, rotl31(k1 * c1) * c2)
@inline mhtail2(h2, k2) = xor(h2, rotl33(k2 * c2) * c1)

@inline function mhblock(h1, h2, k1, k2)
    # dbf[] && print("mhblock($(repr(h1)), $(repr(h2)), $(repr(k1)), $(repr(k2))) => ")
    h1 = (rotl27(mhtail1(h1, k1)) + h2) * 5 + 0x52dce729
    h2 = (rotl31(mhtail2(h2, k2)) + h1) * 5 + 0x38495ab5
    # dbf[] && println(repr(h1), ", ", repr(h2))
    h1, h2
end

@inline mhblock(h1, h2, k1) = mhblock(h1, h2, u64(k1), u64(k1 >>> 64))

@inline function mhbody(nblocks, pnt, h1, h2)
    for i = 1:nblocks
        h1, h2 = mhblock(h1, h2, unsafe_load(pnt), unsafe_load(pnt + 8))
        pnt += 16
    end
    pnt, h1, h2
end

@inline function mhfin(len, h1, h2)
    h1 = xor(h1, u64(len))
    h2 = xor(h2, u64(len))

    h1 += h2
    h2 += h1

    h1 = fmix(h1)
    h2 = fmix(h2)

    h1 += h2
    h1, h1 + h2
end

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

# mmhash128_8 implementation for strings stored little-ending packed in an Unsigned value

# simplified case for when length == 0
@inline function mhfin(h1)
    h2 = h1
    h1 += h2
    h2 += h1
    h1 = fmix(h1)
    h2 = fmix(h2)
    h1 += h2
    h1, h1 + h2
end

mask_val(val, left) = u64(val) & ((UInt64(1) << ((left & 7) << 3)) - 0x1)

# For val that is stored in up to 128 bits, we handle it without any loops
# len can only be 1-7 in this case
_mmhash128(len::Integer, val::Union{UInt32,UInt64}, seed::UInt32) =
    mhfin(len, mhtail1(u64(seed), mask_val(val, len)), u64(seed))

# For val that is stored in up to 128 bits, we handle it without any loops
# len can only be 1-15 in this case
_mmhash128(len::Integer, val::UInt128, seed::UInt32) =
    mhfin(len, mhtail1(u64(seed), len < 8 ? mask_val(val, len) : u64(val)),
          len > 8 ? mhtail2(u64(seed), mask_val(val >> 64, len)) : u64(seed))

@inline function mhbody(nblocks, val::Unsigned, h1, h2)
    for i = 1:nblocks
        h1, h2 = mhblock(h1, h2, u64(val), u64(val>>>64))
        val >>>= 128
    end
    val, h1, h2
end

# Handle values that are more than 128 bits long
function _mmhash128(len::Integer, val::Unsigned, seed::UInt32)
    val, h1, h2 = mhbody(len >>> 4, val, u64(seed), u64(seed))
    if (left = len & 15) > 0
        h1 = mhtail1(h1, left < 8 ? mask_val(val, left) : u64(val))
        left > 8 && (h2 = mhtail2(h2, mask_val(val >> 64, left)))
    end
    mhfin(len, h1, h2)
end

mmhash128_8_a(len::Integer, val::Unsigned, seed::UInt32) =
    len === 0 ? mhfin(seed) : _mmhash128(len, val, seed)

mmhash128_8_c(len::Integer, val::Unsigned, seed::UInt32) = mmhash128_8_a(len, val, seed)

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

up8(val)  = u32(val) << 8
up16(val) = u32(val) << 16
up24(val) = u32(val) << 24
up32(val) = u64(val) << 32
up40(val) = u64(val) << 40
up48(val) = u64(val) << 48
up56(val) = u64(val) << 56

up13b(val) = u128(val) << 104
up14b(val) = u128(val) << 112
up15b(val) = u128(val) << 120

dn6(val) = u8(val >>> 6)
dn12(val) = u8(val >>> 12)
dn18(val) = u8(val >>> 18)

msk6(val) = u8(val & 0x3f)

# Support functions for UTF-8 handling
@inline get_utf8_2(ch) = 0x000080c0 | dn6(ch) | up8(msk6(ch))
@inline get_utf8_3(ch) = 0x008080e0 | dn12(ch) | up8(msk6(dn6(ch))) | up16(msk6(ch))
@inline get_utf8_4(ch) =
    0x808080f0 | dn18(ch) | up8(msk6(dn12(ch))) | up16(msk6(dn6(ch))) | up24(msk6(ch))

# Optimized in-place conversion to UTF-8 for hashing compatibly with isequal / String
@inline shift_n(v, n) = u128(v) << (n<<3)

# if cnt == 0 - 4, bytes must fit in k1
# cnt between 5 - 8, may overflow into k2
# if l == 8 - 12,  bytes must fit in k2
# cnt between 12 - 15, may overflow into k3

@inline function add_utf8(cnt, ch::Char, k1::UInt128)
    v = bswap(reinterpret(UInt32, ch))
    (cnt + ifelse(v < 0x80, 1, ifelse((v & 0xff) < 0xe0, 2, 3 + ((v & 0xff) >= 0xf0))),
     k1 | shift_n(v, cnt))
end

@inline function add_utf8_split(cnt, ch::Char, k1::UInt128)
    v = bswap(reinterpret(UInt32, ch))
    v < 0x80 && return cnt + 1, k1 | shift_n(ch, cnt), u64(0)
    if (v & 0xff) < 0xe0
        nc = cnt + 2
        v1, v2 = cnt == 15 ? (up15b(v), u64(v) >>> 8) : (shift_n(v, cnt), u64(0))
    else
        nc = cnt + 3 + ((v & 0xff) >= 0xf0)
        v1, v2 = cnt == 13 ? (up13b(v), u64(v) >>> 24) :
            cnt == 14 ? (up14b(v), u64(v) >>> 16) :
            (up15b(v), u64(v) >>> 8)
    end
    return (nc, k1 | v1, v2)
end

@inline function add_utf8(cnt, chr, k1::UInt128)
    ch = u32(chr)
    if ch <= 0x7f
        cnt + 1, k1 | shift_n(ch, cnt)
    elseif ch <= 0x7ff
        cnt + 2, k1 | shift_n(get_utf8_2(ch), cnt)
    elseif ch <= 0xffff
        cnt + 3, k1 | shift_n(get_utf8_3(ch), cnt)
    else
        cnt + 4, k1 | shift_n(get_utf8_4(ch), cnt)
    end
end

@inline function add_utf8_split(cnt, chr, k1::UInt128)
    ch = u32(chr)
    ch <= 0x7f && return (cnt + 1, k1 | shift_n(ch, cnt), u64(0))
    # dbf[] && print("add_utf_split($cnt, $(repr(ch)), $(repr(k1))")
    if ch <= 0x7ff
        nc = cnt + 2
        v = get_utf8_2(ch)
        v1, v2 = cnt == 15 ? (up15b(v), u64(v) >>> 8) : (shift_n(v, cnt), u64(0))
    elseif ch <= 0xffff
        nc = cnt + 3
        v = get_utf8_3(ch)
        v1, v2 = cnt == 13 ? (up13b(v), u64(0)) :
                 cnt == 14 ? (up14b(v), u64(v >>> 16)) :
                 (up15b(v), u64(v) >>> 8)
    else
        # This will always go over, may be 1, 2, 3 bytes in second word
        nc = cnt + 4
        v = get_utf8_4(ch)
        # dbf[] && println(" : cnt=$cnt, v=$(repr(v))")
        v1, v2 = cnt == 13 ? (up13b(v), u64(v) >>> 24) :
                 cnt == 14 ? (up14b(v), u64(v) >>> 16) :
                 (up15b(v), u64(v) >>> 8)
    end
    # dbf[] && println(" -> ($nc, $(repr(v1)) => $(repr(k1|v1)), $(repr(v2)))")
    return (nc, k1 | v1, v2)
end

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

# AbstractString MurmurHash3, converts to UTF-8 on the fly
function mmhash128_8_c(str::AbstractString, seed::UInt32)
    k1 = UInt128(0)
    h1 = h2 = u64(seed)
    cnt = len = 0
    @inbounds for ch in str
        if cnt < 13
            cnt, k1 = add_utf8(cnt, ch, k1)
            if cnt == 16
                h1, h2 = mhblock(h1, h2, k1)
                k1 = 0%UInt128
                len += 16
                cnt = 0
            end
        else
            cnt, k1, k2 = add_utf8_split(cnt, ch, k1)
            # When k1 and k2 are full, then hash another block
            if cnt > 15
                h1, h2 = mhblock(h1, h2, k1)
                k1 = k2%UInt128
                len += 16
                cnt &= 15
            end
        end
    end
    # We should now have characters in k1 and k2, and total length in len
    if cnt != 0
        h1 = mhtail1(h1, u64(k1))
        cnt > 8 && (h2 = mhtail2(h2, u64(k1>>>64)))
    end
    mhfin(len + cnt, h1, h2)
end

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

# Note: this is designed to work on the Str/String types, where I know in advance that
# the start of the strings are 8-byte aligned, and it is safe to access a full
# 8-byte chunk always at the end (simply masking off the remaining 1-7 bytes)

@inline mask_load(pnt, left) = unsafe_load(pnt) & ((UInt64(1) << ((left & 7) << 3)) - 0x1)

function mmhash128_8_a(len::Integer, pnt::Ptr, seed::UInt32)
    pnt8, h1, h2 = mhbody(len >>> 4, reinterpret(Ptr{UInt64}, pnt), u64(seed), u64(seed))
    if (left = len & 15) > 0
        h1 = mhtail1(h1, left < 8 ? mask_load(pnt8, left) : unsafe_load(pnt8))
        left > 8 && (h2 = mhtail2(h2, mask_load(pnt8 + 8, left)))
    end
    mhfin(len, h1, h2)
end

function mmhash128_8_a(seed::Integer)
    h1 = fmix(2 * u64(seed))
    h2 = fmix(3 * u64(seed))
    h1 + h2, h1 + 2 * h2
end

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

# Combine bits from k1, k2, k3
@inline shift_mix(shft, k1::T, k2::T, k3::T) where {T<:Union{UInt32,UInt64}} =
    k1 >>> shft | (k2 << (sizeof(T)*8 - shft)),
    k2 >>> shft | (k3 << (sizeof(T)*8 - shft)),
    k3 >>> shft

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

function mmhash128_8_u(len::Integer, unaligned_pnt::Ptr, seed::UInt32)
    # Should optimize handling of short (< 16 byte) unaligned strings
    ulp = reinterpret(UInt, unaligned_pnt)
    pnt = reinterpret(Ptr{UInt64}, ulp & ~UInt(7))
    fin = reinterpret(Ptr{UInt64}, (ulp + len + 0x7) & ~UInt(7)) - 8
    shft = (ulp & UInt(7))<<3
    h1 = h2 = u64(seed)
    k1 = unsafe_load(pnt) # Pick up first 1-7 bytes
    k2 = u64(0)
    while pnt < fin
        k1, k2, k3 = shift_mix(shft, k1, unsafe_load(pnt += 8), unsafe_load(pnt += 8))
        h1, h2 = mhblock(h1, h2, k1, k2)
        k1 = k3
    end
    # We should now have characters in k1 and k2, and total length in len
    if (len & 15) != 0
        h1 = mhtail1(h1, k1)
        (len & 15) > 8 && (h2 = mhtail2(h2, k2))
    end
    mhfin(len, h1, h2)
end

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

# 32-bit MurmurHash3 (see MurmurHash3_x86_32)

@inline xor16(k::UInt32) = xor(k, k >>> 16)
@inline xor13(k::UInt32) = xor(k, k >>> 13)

# I don't think these help the generated code anymore (but they make the code easier to read)
@inline rotl13(k) = rotl(k, 13)
@inline rotl15(k) = rotl(k, 15)
@inline rotl16(k) = rotl(k, 16)
@inline rotl17(k) = rotl(k, 17)
@inline rotl18(k) = rotl(k, 18)
@inline rotl19(k) = rotl(k, 19)

# Constants for mmhash_32
const d1 = 0xcc9e2d51
const d2 = 0x1b873593

@inline fmix(h::UInt32) = xor16(xor13(xor16(h) * 0x85ebca6b) * 0xc2b2ae35)

@inline mhblock(h1, k1) = u32(rotl13(xor(h1, rotl15(k1 * d1) * d2))*0x00005) + 0xe6546b64

@inline function mhbody(nblocks, pnt, h1)
    for i = 1:nblocks
        h1 = mhblock(h1, unsafe_load(pnt))
        pnt += 4
    end
    pnt, h1
end

@inline mhtail32(h, v) = xor(h, rotl15(v * d1) * d2)
@inline mask32(v, res) = v & ifelse(res==1, 0x000ff, ifelse(res==2, 0x0ffff, 0xffffff))

@inline function calc32(len, pnt::Ptr, seed)
    res = len & 3
    res != 0 && (seed = mhtail32(seed, mask32(unsafe_load(pnt), res)))
    fmix(xor(seed, u32(len)))
end

@inline function calc32(len, val::UInt32, seed)
    res = len & 3
    res != 0 && (seed = mhtail32(seed, mask32(val, res)))
    fmix(xor(seed, u32(len)))
end

mmhash32(len, pnt::Ptr, seed::UInt32) =
    calc32(len, mhbody(len >>> 2, reinterpret(Ptr{UInt32}, pnt), seed)...)

# length must be 0-3
mmhash32(len, val::UInt32, seed::UInt32) = calc32(len, val, seed)

# length must be 0-7
mmhash32(len, val::UInt64, seed::UInt32) =
    (len > 3
     ? calc32(len, u32(val>>>32), mhblock(seed, u32(val)))
     : calc32(len, u32(val), seed))

function mmhash32(len, val::Unsigned, seed::UInt32)
    for i = 1:(len>>>2)
        seed = mhblock(seed, u32(val))
        val >>>= 32
    end
    calc32(len, u32(val), seed)
end

@inline function mhfin(len, h1, h2, h3, h4)
    h1 = xor(h1, u32(len))
    h2 = xor(h2, u32(len))
    h3 = xor(h3, u32(len))
    h4 = xor(h4, u32(len))

    h1 += h2; h1 += h3; h1 += h4; h2 += h1; h3 += h1; h4 += h1

    h1 = fmix(h1)
    h2 = fmix(h2)
    h3 = fmix(h3)
    h4 = fmix(h4)

    h1 += h2; h1 += h3; h1 += h4; h2 += h1; h3 += h1; h4 += h1

    up32(h2) | h1, up32(h4) | h3
end

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

# Calculate MurmurHash for 32-bit platforms

# Constants for mmhash128_4
const e1 = 0x239b961b
const e2 = 0xab0e9789
const e3 = 0x38b34ae5
const e4 = 0xa1e38b93

mhtail4_1(h, v) = xor(h, rotl15(v * e1) * e2)
mhtail4_2(h, v) = xor(h, rotl16(v * e2) * e3)
mhtail4_3(h, v) = xor(h, rotl17(v * e3) * e4)
mhtail4_4(h, v) = xor(h, rotl18(v * e4) * e1)

@inline mask_v32(val, left) = u32(val) & ((UInt32(1) << ((left & 3) << 3)) - 0x1)

@inline function mhblock(h1, h2, h3, h4, k1, k2, k3, k4)
    h1 = (rotl19(mhtail4_1(h1, k1)) + h2)*5 + 0x561ccd1b
    h2 = (rotl17(mhtail4_2(h2, k2)) + h3)*5 + 0x0bcaa747
    h3 = (rotl15(mhtail4_3(h3, k3)) + h4)*5 + 0x96cd1c35
    h4 = (rotl13(mhtail4_4(h4, k4)) + h1)*5 + 0x32ac3b17
    h1, h2, h3, h4
end

@inline function mhbody(nblocks, pnt, h1, h2, h3, h4)
    for i = 1:nblocks
        h1, h2, h3, h4 =
            mhblock(h1, h2, h3, h4,
                    unsafe_load(pnt), unsafe_load(pnt+4), unsafe_load(pnt+8), unsafe_load(pnt+12))
        pnt += 16
    end
    pnt, h1, h2, h3, h4
end

# degenerate case, hash for 0 length strings, based entirely on seed
function mmhash128_4(seed::UInt32)
    h = fmix(5*seed)*5
    up32(h) | fmix(4*seed)*4, up32(h) | h
end

# For val that is stored in up to 32 bits, we handle it without any loops
# len can only be 1-3 in this case
function mmhash128_4(len::Integer, val::UInt32, seed::UInt32)
    len == 0 && return mmhash128_4(seed)
     mhfin(len, mhtail4_1(seed, mask_v32(val, len)), seed, seed, seed)
end

# For val that is stored in up to 64 bits, we handle it without any loops
# len can only be 1-7 in this case
function mmhash128_4(len::Integer, val::UInt64, seed::UInt32)
    len == 0 && return mmhash128_4(seed)
    mhfin(len, mhtail4_1(seed, len < 4 ? mask_v32(val, len) : u32(val)),
          len > 4 ? mhtail4_2(seed, mask_v32(val>>>32, len)) : seed,
          seed, seed)
end

# For val that is stored in up to 128 bits, we handle it without any loops
# len can only be 1-15 in this case
function mmhash128_4(len::Integer, val::UInt128, seed::UInt32)
    len == 0 && return mmhash128_4(seed)
    h2 = h3 = h4 = seed
    if len > 4
        val >>>= 32
        h2  = mhtail4_2(h2, len < 8 ? mask_v32(val, len) : u32(val))
        if len > 8
            val >>>= 32
            h3  = mhtail4_3(h3, len < 12 ? mask_v32(val, len) : u32(val))
            len > 12 && (h4 = mhtail4_4(h4, mask_v32(val>>>32, len)))
        end
    end
    mhfin(len, mhtail4_1(seed, len < 4 ? mask_v32(val, len) : u32(val)), h2, h3, h4)
end

function mmhash128_4(len::Integer, val::Unsigned, seed::UInt32)
    h1 = h2 = h3 = h4 = seed
    for i = 1:(len>>>4)
        h1, h2, h3, h4 =
            mhblock(h1, h2, h3, h4, u32(val), u32(val>>>32), u32(val>>>64), u32(val>>>96))
        val >>>= 128
    end
    if (left = len & 15) != 0
        # Pick up 32-bit
        h1  = mhtail4_1(h1, left < 4 ? mask_v32(val, left) : u32(val))
        if left > 4
            val >>>= 32
            h2  = mhtail4_2(h2, left < 8 ? mask_v32(val, left) : u32(val))
            if left > 8
                val >>>= 32
                h3  = mhtail4_3(h3, left < 12 ? mask_v32(val, left) : u32(val))
                left > 12 && (h4 = mhtail4_4(h4, mask_v32(val>>>32, left)))
            end
        end
    end
    mhfin(len, h1, h2, h3, h4)
end

function mmhash128_4(len::Integer, pnt::Ptr{UInt32}, seed::UInt32)
    pnt, h1, h2, h3, h4 = mhbody(len >>> 4, pnt, seed, seed, seed, seed)
    if (left = len & 15) != 0
        h1  = mhtail4_1(h1, unsafe_load(pnt))
        if left > 4
            h2  = mhtail4_2(h2, unsafe_load(pnt+4))
            if left > 8
                h3  = mhtail4_3(h3, unsafe_load(pnt+8))
                left > 12 && (h4  = mhtail4_4(h4, unsafe_load(pnt+12)))
            end
        end
    end
    mhfin(len, h1, h2, h3, h4)
end

mmhash128_4(len::Integer, pnt::Ptr, seed::UInt32) =
    mmhash128_4(len, reinterpret(Ptr{UInt32}, pnt), seed)

# Handle value stored in an unsigned value
import Base.GC: @preserve

# AbstractString MurmurHash3, converts to UTF-8 on the fly (not optimized yet!)
function mmhash128_4(s::AbstractString, seed::UInt32)
    str = string(s)
    @preserve str mmhash128_4(sizeof(str), pointer(str), seed)
end

@inline function get_utf8(cnt, ch)
    if ch <= 0x7f
        cnt + 1, u32(ch)
    elseif ch <= 0x7ff
        cnt + 2, get_utf8_2(ch)
    elseif ch <= 0xffff
        cnt + 3, get_utf8_3(ch)
    else
        cnt + 4, get_utf8_4(ch)
    end
end

const mmhash128_c = @static sizeof(Int) == 8 ? mmhash128_8_c : mmhash128_4
const mmhash128_a = @static sizeof(Int) == 8 ? mmhash128_8_a : mmhash128_4
const mmhash128_u = @static sizeof(Int) == 8 ? mmhash128_8_c : mmhash128_4 # Todo: fix unaligned

end # module MurmurHash3


text:
# Script to showcase how to use the model scf

push!(LOAD_PATH, "../src/")
using modelscf
using HDF5
using LinearAlgebra
using PyPlot

# flag to save the data

# number of MD simulations and
Nit = 3

# getting all the parameters
dx = 0.5;
Nunit = 8;   # number of units
Lat = 10;     # size of the lattice
Ls = Nunit*Lat;
Ns = round(Integer, Ls / dx); # number of discretization points

# using the default values in Lin's code
YukawaK = 0.0100
n_extra = 10;
epsil0 = 10.0;
T_elec = 100.0;

kb = 3.1668e-6;
au2K = 315774.67;
Tbeta = au2K / T_elec;

betamix = 0.5;
mixdim = 10;

Ndist = 1;
Natoms = round(Integer, Nunit / Ndist);

sigma  = ones(Natoms,1)*(1.0);  # insulator
omega  = ones(Natoms,1)*0.03;
Eqdist = ones(Natoms,1)*10.0;
mass   = ones(Natoms,1)*42000.0;
nocc   = ones(Natoms,1)*2;     # number of electrons per atom
Z      = nocc;

# defining the position of the nuclei
R = reshape([ (j-0.5)*Lat*Ndist+dx for j=1:Natoms ], Natoms, 1)

# creating an atom structure
atoms = Atoms(Natoms, R, sigma,  omega,  Eqdist, mass, Z, nocc);
# allocating a Hamiltonian
ham = Ham(Lat, Nunit, n_extra, dx, atoms,YukawaK, epsil0, Tbeta)
# total number of occupied orbitals
Nocc = round(Integer, sum(atoms.nocc) / ham.nspin);

mixOpts = andersonMixOptions(ham.Ns, betamix, mixdim )
eigOpts = eigOptions(1.e-10, 1000, "eigs");
#eigOpts = eigOptions(1.e-10, 1000, "lobpcg_sep");
scfOpts = scfOptions(1.e-8, 300, eigOpts, mixOpts)

# initialize the potentials within the Hemiltonian, setting H[\rho_0]
init_pot!(ham, Nocc)

Profile.clear()
@profile VtoterrHist = scf!(ham, scfOpts)
Juno.profiletree()
Juno.profiler()

@profiler scf!(ham, scfOpts)
# running the scf iteration

if VtoterrHist[end] > scfOpts.SCFtol
    println("convergence not achieved!! ")
end

println(length(VtoterrHist))

# we compute the forces
get_force!(ham)
# computing the energy
Vhar = hartree_pot_bc(ham.rho+ham.rhoa,ham);

# NOTE: ham.Fband is only the band energy here.  The real total energy
# is calculated using the formula below:
Etot = ham.Eband + 1/2*sum((ham.rhoa-ham.rho).*Vhar)*dx;

# we need to add a proper title
plot(ham.ev)

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
module Depository
using SolveBio

function list()
    return SolveBio.request("GET", "v1/depositories")
end

function get(id)
    path = "v1/depositories/$id"
    return SolveBio.request("GET", path)
end

function versions(id)
    path = "v1/depositories/$id/versions"
    return SolveBio.request("GET", path)
end

end


text:
"""
    AbstractBridge

Subtype of [`MathOptInterface.Bridges.AbstractBridge`](@ref) for constraint
bridges.
"""
abstract type AbstractBridge <: MOIB.AbstractBridge end

# TODO [breaking] merge with `Bridges.Variable.IndexInVector` into `Bridges.IndexInVector`
"""
    IndexInVector

Index of variable in vector of variables.
"""
struct IndexInVector
    value::Int
end

"""
    bridge_constraint(BT::Type{<:AbstractBridge}, model::MOI.ModelLike,
                      func::AbstractFunction, set::MOI.AbstractSet)

Bridge the constraint `func`-in-`set` using bridge `BT` to `model` and returns
a bridge object of type `BT`. The bridge type `BT` should be a concrete type,
that is, all the type parameters of the bridge should be set. Use
[`concrete_bridge_type`](@ref) to obtain a concrete type for given function
and set types.
"""
function bridge_constraint end

"""
    MOI.get(b::AbstractBridge, ::MOI.NumberOfVariables)

The number of variables created by the bridge `b` in the model.
"""
MOI.get(::AbstractBridge, ::MOI.NumberOfVariables) = 0

"""
    MOI.get(b::AbstractBridge, ::MOI.NumberOfVariables)

The list of variables created by the bridge `b` in the model.
"""
MOI.get(::AbstractBridge, ::MOI.ListOfVariableIndices) = MOI.VariableIndex[]

"""
    MOI.supports_constraint(BT::Type{<:AbstractBridge}, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging `F`-in-`S` constraints.
"""
function MOI.supports_constraint(
    ::Type{<:AbstractBridge},
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

"""
    added_constrained_variable_types(BT::Type{<:MOI.Bridges.Constraint.AbstractBridge},
                                     F::Type{<:MOI.AbstractFunction},
                                     S::Type{<:MOI.AbstractSet})

Return a list of the types of constrained variables that bridges of type `BT`
add for bridging `F`-in-`S` constraints. This falls back to
`added_constrained_variable_types(concrete_bridge_type(BT, F, S))`
so bridges should not implement this method.
"""
function MOIB.added_constrained_variable_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.added_constrained_variable_types(concrete_bridge_type(BT, F, S))
end

"""
    added_constraint_types(BT::Type{<:MOI.Bridges.Constraint.AbstractBridge},
                           F::Type{<:MOI.AbstractFunction},
                           S::Type{<:MOI.AbstractSet})

Return a list of the types of constraints that bridges of type `BT` add for
bridging `F`-in-`S` constraints. This falls back to
`added_constraint_types(concrete_bridge_type(BT, F, S))`
so bridges should not implement this method.
"""
function MOIB.added_constraint_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.added_constraint_types(concrete_bridge_type(BT, F, S))
end

"""
    concrete_bridge_type(
        BT::Type{<:AbstractBridge},
        F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet}
    )::DataType

Return the concrete type of the bridge supporting `F`-in-`S` constraints. This
function can only be called if `MOI.supports_constraint(BT, F, S)` is `true`.

## Examples

As a [`MathOptInterface.SingleVariable`](@ref)-in-[`MathOptInterface.Interval`](@ref)
constraint is bridged into a
[`MathOptInterface.SingleVariable`](@ref)-in-[`MathOptInterface.GreaterThan`](@ref)
and a
[`MathOptInterface.SingleVariable`](@ref)-in-[`MathOptInterface.LessThan`](@ref)
by the [`SplitIntervalBridge`](@ref):

```jldoctest; setup=:(using MathOptInterface; const MOI = MathOptInterface)
MOI.Bridges.Constraint.concrete_bridge_type(
    MOI.Bridges.Constraint.SplitIntervalBridge{Float64},
    MOI.SingleVariable,
    MOI.Interval{Float64},
)

# output

MathOptInterface.Bridges.Constraint.SplitIntervalBridge{Float64,MathOptInterface.SingleVariable,MathOptInterface.Interval{Float64},MathOptInterface.GreaterThan{Float64},MathOptInterface.LessThan{Float64}}
```
"""
function concrete_bridge_type(
    bridge_type::DataType,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return bridge_type
end

function concrete_bridge_type(
    b::MOIB.AbstractBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return concrete_bridge_type(MOIB.bridge_type(b, F, S), F, S)
end


text:
module Setters

using Test

using Setfield: set
using Unitful
using UnitfulAtomic

using QuantumESPRESSOBase.Inputs.PWscf
using QuantumESPRESSOBase.Setters

@testset "Test `VerbositySetter` on `ControlNamelist`" begin
    control = ControlNamelist()
    @test (  # Default values
        control.verbosity,
        control.wf_collect,
        control.tstress,
        control.tprnfor,
        control.disk_io,
    ) == ("low", true, false, false, "low")
    control = set(control, VerbositySetter(:high))
    @test (
        control.verbosity,
        control.wf_collect,
        control.tstress,
        control.tprnfor,
        control.disk_io,
    ) == ("high", true, true, true, "high")
    control = set(control, VerbositySetter(:low))
    @test (
        control.verbosity,
        control.wf_collect,
        control.tstress,
        control.tprnfor,
        control.disk_io,
    ) == ("low", false, false, false, "low")
end # testset

@testset "Test `FiniteTemperatureSetter` on `SystemNamelist`" begin
    system = SystemNamelist()
    @test (system.occupations, system.degauss, system.smearing) ==
          ("fixed", 0.0, "gaussian")
    system = set(system, FiniteTemperatureSetter(0.0019000869380733452))
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.0019000869380733452, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(3e8 * u"μK"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.0019000869380733452, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(2.5851999786e-8 * u"MeV"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.001900086938034304, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(0.00095004346903668 * u"hartree"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.00190008693807336, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(6.250985736 * u"THz"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.0019000869377698247, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(20851.044012 * u"m^-1"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.001900086937837879, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(208.51044012 * u"cm^-1"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.001900086937837879, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(4.141947e-24 * u"kJ"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.0019000869380663146, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(4.6085375610000005e-35 * u"g"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.0019000869377760491, "fermi-dirac")
    system = set(system, FiniteTemperatureSetter(4.6085375610000005e-32 * u"mg"))  # 300K
    @test (system.occupations, system.degauss, system.smearing) ==
          ("smearing", 0.001900086937776049, "fermi-dirac")
end # testset

end


text:
struct PoissonProcess end
generate(::Type{PoissonProcess}, np::Int) = rand(np, 2)

function pointstogrid(coords; gridsize = 100)
    intcoords = map(x -> Int(ceil(x * gridsize)), coords)
    z = zeros(gridsize, gridsize)
    for c in eachrow(intcoords)
        z[c...] = 1
    end
    return z
end


text:
module SingleStackUtils

export get_vars_from_nodal_stack, get_vars_from_element_stack

using OrderedCollections
using StaticArrays

using ..DGmethods
using ..DGmethods.Grids
using ..MPIStateArrays
using ..VariableTemplates

"""
    get_vars_from_nodal_stack(
        grid::DiscontinuousSpectralElementGrid{T, dim, N},
        Q::MPIStateArray,
        vars;
        vrange::UnitRange = 1:size(Q, 3),
        i::Int = 1,
        j::Int = 1,
        exclude::Vector{String} = String[],
    ) where {T, dim, N}

Return a dictionary whose keys are the `flattenednames()` of the variables
specified in `vars` (as returned by e.g. `vars_state_conservative()`), and
whose values are arrays of the values for that variable along the vertical
dimension in `Q`. Only a single element is expected in the horizontal as
this is intended for the single stack configuration and `i` and `j` identify
the horizontal nodal coordinates.

Variables listed in `exclude` are skipped.
"""
function get_vars_from_nodal_stack(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
    Q::MPIStateArray,
    vars;
    vrange::UnitRange = 1:size(Q, 3),
    i::Int = 1,
    j::Int = 1,
    exclude::Vector{String} = String[],
) where {T, dim, N}

    # extract grid information and bring `Q` to the host if needed
    FT = eltype(Q)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    nvertelem = grid.topology.stacksize
    host_Q = Array ∈ typeof(Q).parameters ? Q.realdata : Array(Q.realdata)

    # set up the dictionary to be returned
    var_names = flattenednames(vars)
    stack_vals = OrderedDict()
    num_vars = varsize(vars)
    vars_wanted = Int[]
    for vi in 1:num_vars
        if !(var_names[vi] in exclude)
            stack_vals[var_names[vi]] = FT[]
            push!(vars_wanted, vi)
        end
    end

    # extract values from `host_Q`
    for ev in vrange
        for k in 1:Nqk
            ijk = i + Nq * ((j - 1) + Nq * (k - 1))
            for v in vars_wanted
                push!(stack_vals[var_names[v]], host_Q[ijk, v, ev])
            end
        end
    end

    return stack_vals
end

"""
    get_vars_from_element_stack(
        grid::DiscontinuousSpectralElementGrid{T, dim, N},
        Q::MPIStateArray,
        vars;
        vrange::UnitRange = 1:size(Q, 3),
        exclude::Vector{String} = String[],
    ) where {T, dim, N}
        exclude = [],
    )

Return an array of [`get_vars_from_nodal_stack()`](@ref)s whose dimensions
are the number of nodal points per element in the horizontal plane.

Variables listed in `exclude` are skipped.
"""
function get_vars_from_element_stack(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
    Q::MPIStateArray,
    vars;
    vrange::UnitRange = 1:size(Q, 3),
    exclude::Vector{String} = String[],
) where {T, dim, N}
    Nq = N + 1
    return [
        get_vars_from_nodal_stack(
            grid,
            Q,
            vars,
            vrange = vrange,
            i = i,
            j = j,
            exclude = exclude,
        ) for i in 1:Nq, j in 1:Nq
    ]
end

end # module


text:
function OrekitStepNormalizer(arg0::jdouble, arg1::OrekitFixedStepHandler)
    return OrekitStepNormalizer((jdouble, OrekitFixedStepHandler), arg0, arg1)
end

function handle_step(obj::OrekitStepNormalizer, arg0::OrekitStepInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (OrekitStepInterpolator, jboolean), arg0, arg1)
end

function init(obj::OrekitStepNormalizer, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function requires_dense_output(obj::OrekitStepNormalizer)
    return jcall(obj, "requiresDenseOutput", jboolean, ())
end


text:
# test AprilTag detections

using Images, ImageDraw, ImageMagick
using AprilTags

# ENV["PYTHON"] = "python"
# using PyCall
#
# @pyimport numpy as np
# @pyimport cv2

# deprecated, use AprilTags.drawTags instead
function showImage(image, tags, K)
    # Convert image to RGB
    imageCol = RGB.(image)
    #traw color box on tag corners
    foreach(tag->drawTagBox!(imageCol, tag, width = 2, drawReticle = false), tags)
    foreach(tag->drawTagAxes!(imageCol,tag, K), tags)
    imageCol
end

global detector = AprilTagDetector()

# include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","cameraUtils.jl"))

cx = 327.986
cy = 198.066
fx = 349.982
fy = 349.982
K = [-fx 0  cx;
0 fy cy]

img = load(dirname(@__FILE__)*"/example_tag_image.jpeg");
tags = detector(img)
showImage(img, tags, K)

tags

cTt = homographytopose(tags[1].H, fx,fy, cx, cy, taglength=0.172)

# getAprilTagImage(1, AprilTags.tag36h11)

using CoordinateTransformations

#


text:
#####################################
#ScalarExpr##########################
#####################################

immutable ScalarExpr <: AbstractScalar
	ex::Expr
end

isequal(a::ScalarExpr, b::ScalarExpr) = isequal(a.ex, b.ex)
==(a::ScalarExpr, b::ScalarExpr) = ==(a.ex, b.ex)

convert{N<:Union(InnerProduct, Number)}(::Type{ScalarExpr}, i::N) = ScalarExpr(:(1*$i))
promote_rule{N<:Union(InnerProduct, Number)}(::Type{ScalarExpr}, ::Type{N}) = ScalarExpr
promote_rule{N<:Union(InnerProduct, Number), T}(::Type{Array{ScalarExpr,T}}, ::Type{Array{N,T}}) = Array{ScalarExpr,T}
promote_rule{N<:Union(InnerProduct, Number), T}(::Type{SparseMatrixCSC{ScalarExpr,T}}, ::Type{SparseMatrixCSC{N,T}}) = SparseMatrixCSC{ScalarExpr,T}

scalar(ex::Expr) = ScalarExpr(ex)
scalar(x) = convert(ScalarExpr, x)

#####################################
#Misc Functions######################
#####################################
qexpr(s::ScalarExpr) = s.ex
qexpr(x) = x
qeval(f::Function, s::ScalarExpr) = eval(qreduce(f, s.ex))
qeval(f::Function, i::InnerProduct) = eval(f(i.bra, i.ket))
qeval(f::Function, n::Number) = n

function qreduce(f::Function, ex::Expr)
	ex = copy(ex)
	for i=1:length(ex.args)
		if typeof(ex.args[i])<:InnerProduct
			ex.args[i] = f(ex.args[i].bra, ex.args[i].ket)
		elseif typeof(ex.args[i])==Expr
			ex.args[i] = qreduce(f, ex.args[i])
		end
	end
	return ex
end

#####################################
#Arithmetic Operations###############
#####################################
function ^(d::DiracCoeff, n::Integer)
	if n==1
		return scalar(d)
	elseif n==0
		return scalar(1)
	else
		return scalar(:($(qexpr(d))^$(qexpr(n))))
	end
end

function ^(a::DiracCoeff, b::DiracCoeff)
	if b==1
		return scalar(a)
	elseif b==0
		return scalar(1)
	else
		return scalar(:($(qexpr(a))^$(qexpr(b))))
	end
end

function *(a::DiracCoeff, b::DiracCoeff)
	if a==1
		return scalar(b)
	elseif b==1
		return scalar(a)
	elseif a==0 || b==0
		return scalar(0)
	else
		return scalar(:($(qexpr(a))*$(qexpr(b))))
	end
end

function +(a::DiracCoeff, b::DiracCoeff)
	if a==0
		return scalar(b)
	elseif b==0
		return scalar(a)
	else
		scalar(:($(qexpr(a))+$(qexpr(b))))
	end
end

function -(a::DiracCoeff, b::DiracCoeff)
	if a==0
		return scalar(-b)
	elseif b==0
		return scalar(a)
	elseif a==b
		return scalar(0)
	else
		return scalar(:($(qexpr(a))-$(qexpr(b))))
	end
end

-(d::DiracCoeff) = ScalarExpr(:(-$(qexpr(d))))
-(s::ScalarExpr) = length(s.ex.args)==2 && s.ex.args[1]==:- ? ScalarExpr(s.ex.args[2]) :  ScalarExpr(:(-$(qexpr(s))))

abs(i::InnerProduct) = ScalarExpr(:(abs($i)))
abs(s::ScalarExpr) = length(s.ex.args)==2 && s.ex.args[1]==:abs ? s :  ScalarExpr(:(abs($(qexpr(s)))))

exp(a::AbstractScalar) = ScalarExpr(:(exp($(qexpr(a)))))

function /(a::DiracCoeff, b::DiracCoeff)
	if a==0
		return b
	elseif b==0
		error("Zero in denominator of ScalarExpr")
	elseif a==b
		return 1
	else
		ScalarExpr(:($(qexpr(a))/$(qexpr(b))))
	end
end

conj(s::ScalarExpr)	= length(s.ex.args)==2 && s.ex.args[1]==:conj ? ScalarExpr(s.ex.args[2]) :  ScalarExpr(:(conj($(qexpr(s)))))
ctranspose(s::ScalarExpr) = conj(s)

for op=(:*,:-,:+,:/,:^) #define for elementwise operators
	elop = symbol(string(:.) * string(op))
	@eval ($elop)(a::DiracCoeff, b::DiracCoeff) = ($op)(a,b)
	@eval ($elop)(a::AbstractScalar, arr::AbstractArray) = broadcast(($op), [a], arr)
	@eval ($elop)(arr::AbstractArray, a::AbstractScalar) = broadcast(($op), arr, [a])
end

*(a::AbstractScalar, arr::AbstractArray) = a.*arr
*(arr::AbstractArray, a::AbstractScalar) = arr.*a
/(arr::AbstractArray, a::AbstractScalar) = arr./a


text:
function ueberpruefeerror(Fehler,PhantomVektor,RekonstruktionVektor)

tol = 1e-12;

if norm(Fehler) - norm(PhantomVektor-RekonstruktionVektor) < tol
    println("Beachtet, dass wir Vektoren betrachten. Dafür müsst ihr die Länge der Vektoren betrachten.")
    println("Zudem habt ihr den absolten Fehler berechnet. Berechnet den Fehler relativ zum originalen Wert.")
elseif Fehler - norm(PhantomVektor-RekonstruktionVektor) < tol
    println("Ihr habt den absolten Fehler berechnet. Berechnet den Fehler relativ zum originalen Wert.")
elseif Fehler - norm(PhantomVektor-RekonstruktionVektor)/norm(PhantomVektor) < tol
    println("Fast richtig. Lasst euch das Ergebnis in Prozent ausgeben.")
elseif Fehler - norm(PhantomVektor-RekonstruktionVektor)/norm(PhantomVektor)*100 < tol
    Fehler2 = round(Fehler;digits=4);
    println("Eure Formel für den relativen Fehler ist korrekt! Der Fehler beträgt: ", Fehler2)
else
    println("Googlt noch einmal nach Fehlerarten, die ihr kennt und schaut in Workspace nach welche Variablen ihr habt, die euch hier weiterhelfen können.")
end
      # return Fehler2;

end


text:
# Values from modbus.h and unit-test.h files of libmodbus
#
if !((@isdefined _MODBUS_CONSTANTS_) && _MODBUS_CONSTANTS_)
    _MODBUS_CONSTANTS_ = true

    LM = LibModbus

    const MODBUS_BROADCAST_ADDRESS = 0

    const MODBUS_MAX_READ_BITS = 2000
    const MODBUS_MAX_WRITE_BITS = 1968
    const MODBUS_MAX_READ_REGISTERS = 125
    const MODBUS_MAX_WRITE_REGISTERS = 123

    # Protocol exceptions

    const MODBUS_ENOBASE =  112345678

    const EMBXILFUN = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_ILLEGAL_FUNCTION))
    const EMBXILADD = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_ILLEGAL_DATA_ADDRESS))
    const EMBXILVAL = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_ILLEGAL_DATA_VALUE))
    const EMBXSFAIL = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_SLAVE_OR_SERVER_FAILURE))
    const EMBXACK =   (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_ACKNOWLEDGE))
    const EMBXSBUSY = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_SLAVE_OR_SERVER_BUSY))
    const EMBXNACK =  (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_NEGATIVE_ACKNOWLEDGE))
    const EMBXMEMPAR = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_MEMORY_PARITY))
    const EMBXGPATH = (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_GATEWAY_PATH))
    const EMBXGTAR =  (MODBUS_ENOBASE + UInt8(LM.MODBUS_EXCEPTION_GATEWAY_TARGET))

    # Native libmodbus error codes
    const EMBBADCRC =   (EMBXGTAR + 1)
    const EMBBADDATA =  (EMBXGTAR + 2)
    const EMBBADEXC =   (EMBXGTAR + 3)
    const EMBUNKEXC =   (EMBXGTAR + 4)
    const EMBMDATA =    (EMBXGTAR + 5)
    const EMBBADSLAVE =  (EMBXGTAR + 6)

    # Values from tests/unit-test-server.h file of libmodbus
    const UT_BITS_ADDRESS = 0x130
    const UT_BITS_NB = 0x25
    const UT_BITS_TAB = [ 0xCD, 0x6B, 0xB2, 0x0E, 0x1B ]

    const UT_INPUT_BITS_ADDRESS = 0x1C4
    const UT_INPUT_BITS_NB = 0x16
    const UT_INPUT_BITS_TAB = [ 0xAC, 0xDB, 0x35 ]

    const UT_REGISTERS_ADDRESS = 0x160
    const UT_REGISTERS_NB = 0x3
    const UT_REGISTERS_NB_MAX = 0x20
    const UT_REGISTERS_TAB = [ 0x022B, 0x0001, 0x0064 ]

    const UT_INPUT_REGISTERS_ADDRESS = 0x108
    const UT_INPUT_REGISTERS_NB = 0x1
    const UT_INPUT_REGISTERS_TAB = [ 0x000A ]

    const INVALID_SERVER_ID = 18
    const UT_REGISTERS_ADDRESS_SPECIAL = 0x0170
    const UT_REGISTERS_ADDRESS_INVALID_TID_OR_SLAVE = 0x0171
    const UT_REGISTERS_ADDRESS_SLEEP_500_MS = 0x0172
    const UT_REGISTERS_ADDRESS_BYTE_SLEEP_5_MS = 0x0173
    const UT_REGISTERS_NB_SPECIAL = 0x0002

    const NB_REPORT_SLAVE_ID = 10
end
nothing


text:
struct SimulationResults
    global_index::MAPPING_DICT
    bus_lookup::Dict{Int, Int}
    system::PSY.System
    time_log::Dict{Symbol, Any}
    solution::SciMLBase.AbstractODESolution
    function SimulationResults(
        inputs::SimulationInputs,
        system::PSY.System,
        time_log,
        solution,
    )
        new(make_global_state_map(inputs), get_lookup(inputs), system, time_log, solution)
    end
end

get_global_index(res::SimulationResults) = res.global_index
get_bus_count(res::SimulationResults) = get_n_buses(res.system)
get_bus_lookup(res::SimulationResults) = res.bus_lookup
get_system(res::SimulationResults) = res.system
get_solution(res::SimulationResults) = res.solution

"""
Internal function to obtain as a Vector of Float64 of a specific state. It receives the solution and the
global index for a state.

"""
function _post_proc_state_series(solution, ix::Int, dt::Union{Nothing, Float64})
    if dt === nothing
        ix_t = unique(i -> solution.t[i], eachindex(solution.t))
        ts = solution.t[ix_t]
        state = solution[ix, ix_t]
    else
        ts = range(0, stop = solution.t[end], step = dt)
        state = solution(collect(ts); idxs = ix)
    end
    return ts, state
end

"""
Function to obtain the state time series of a specific state. It receives the simulation, and a tuple
containing the name of the Dynamic Device and the symbol of the state.
"""
function post_proc_state_series(
    res::SimulationResults,
    ref::Tuple{String, Symbol},
    dt::Union{Nothing, Float64},
)
    global_state_index = get_global_index(res)
    ix = get(global_state_index[ref[1]], ref[2], 0)
    return _post_proc_state_series(get_solution(res), ix, dt)
end

"""
Function to obtain voltage and output currents for a dynamic device. It receives the simulation, and the name
of the Dynamic Device.
"""
function post_proc_voltage_current_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)::NTuple{5, Vector{Float64}}
    #Note: Type annotation since get_dynamic_injector is type unstable and solution is Union{Nothing, DAESol}
    system = get_system(res)
    bus_lookup = get_bus_lookup(res)
    n_buses = length(bus_lookup)
    solution = res.solution
    device = PSY.get_component(PSY.StaticInjection, system, name)
    if isnothing(device)
        error("Device $(name) not found in the system")
    end
    bus_ix = get(bus_lookup, PSY.get_number(PSY.get_bus(device)), -1)
    ts, V_R, V_I = post_proc_voltage_series(solution, bus_ix, n_buses, dt)
    dyn_device = PSY.get_dynamic_injector(device)
    _, I_R, I_I = compute_output_current(res, dyn_device, V_R, V_I, dt)
    return ts, V_R, V_I, I_R, I_I
end

"""
Function to obtain voltage using the bus index (and not the bus number). It receives the solution, the bus index and
the total number of buses.

"""
function post_proc_voltage_series(
    solution,
    bus_ix::Int,
    n_buses::Int,
    dt::Union{Nothing, Float64},
)
    bus_ix < 0 && error("Bus number $(bus_number) not found.")
    ts, V_R = _post_proc_state_series(solution, bus_ix, dt)
    _, V_I = _post_proc_state_series(solution, bus_ix + n_buses, dt)
    return collect(ts), collect(V_R), collect(V_I)
end

"""
Function to compute the real current output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_real_current_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    ts, _, _, I_R, _ = post_proc_voltage_current_series(res, name, dt)
    return ts, I_R
end
"""
Function to compute the imaginary current output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_imaginary_current_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    ts, _, _, _, I_I = post_proc_voltage_current_series(res, name, dt)
    return ts, I_I
end

"""
Function to compute the active power output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_activepower_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    ts, V_R, V_I, I_R, I_I = post_proc_voltage_current_series(res, name, dt)
    return ts, V_R .* I_R + V_I .* I_I
end

"""
Function to compute the active power output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_reactivepower_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    ts, V_R, V_I, I_R, I_I = post_proc_voltage_current_series(res, name, dt)
    return ts, V_I .* I_R - V_R .* I_I
end

"""
Function to compute the field current output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_field_current_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    system = get_system(res)
    bus_lookup = get_bus_lookup(res)
    n_buses = length(bus_lookup)
    solution = res.solution
    device = PSY.get_component(PSY.StaticInjection, system, name)
    bus_ix = get(bus_lookup, PSY.get_number(PSY.get_bus(device)), -1)
    ts, V_R, V_I = post_proc_voltage_series(solution, bus_ix, n_buses, dt)
    dyn_device = PSY.get_dynamic_injector(device)
    _, I_fd = compute_field_current(res, dyn_device, V_R, V_I, dt)
    return ts, I_fd
end

"""
Function to compute the field voltage output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
string name of the Dynamic Injection device.

"""
function post_proc_field_voltage_series(
    res::SimulationResults,
    name::String,
    dt::Union{Nothing, Float64},
)
    system = get_system(res)
    device = PSY.get_component(PSY.StaticInjection, system, name)
    dyn_device = PSY.get_dynamic_injector(device)
    ts, Vf = compute_field_voltage(res, dyn_device, dt)
    return ts, Vf
end

"""
    get_state_series(
        res::SimulationResults,
        ref::Tuple{String, Symbol};
        dt::Union{Nothing, Float64} = nothing
    )
    end

Function to obtain series of states out of DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `ref:Tuple{String, Symbol}` : Tuple used to identify the dynamic device, via its name, as a `String`, and the associated state as a `Symbol`.
"""
function get_state_series(res::SimulationResults, ref::Tuple{String, Symbol}; dt = nothing)
    return post_proc_state_series(res, ref, dt)
end

"""
    get_voltage_magnitude_series(
        res::SimulationResults,
        bus_number::Int
    )

Function to obtain the voltage magnitude series out of the DAE Solution.

# Arguments:

- `res::SimulationResults` : Simulation Results object that contains the solution
- `bus_number::Int` : Bus number identifier
"""
function get_voltage_magnitude_series(res::SimulationResults, bus_number::Int; dt = nothing)
    n_buses = get_bus_count(res)
    bus_ix = get(get_bus_lookup(res), bus_number, 0)
    ts, V_R, V_I = post_proc_voltage_series(res.solution, bus_ix, n_buses, dt)
    return ts, sqrt.(V_R .^ 2 .+ V_I .^ 2)
end

"""
    get_voltage_angle_series(
        res::SimulationResults,
        bus_number::Int
    )

Function to obtain the voltage angle series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `bus_number::Int` : Bus number identifier
"""
function get_voltage_angle_series(res::SimulationResults, bus_number::Int; dt = nothing)
    n_buses = get_bus_count(res)
    bus_ix = get(get_bus_lookup(res), bus_number, 0)
    ts, V_R, V_I = post_proc_voltage_series(res.solution, bus_ix, n_buses, dt)
    return ts, atan.(V_I, V_R)
end

"""
    get_real_current_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the real current time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_real_current_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_real_current_series(res, name, dt)
end

"""
    get_imaginary_current_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the imaginary current time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_imaginary_current_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_imaginary_current_series(res, name, dt)
end

"""
    get_activepower_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the active power output time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_activepower_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_activepower_series(res, name, dt)
end

"""
    get_reactivepower_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the reactive power output time series of a Dynamic Injection series out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_reactivepower_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_reactivepower_series(res, name, dt)
end

"""
    get_field_current_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the field current time series of a Dynamic Generator out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_field_current_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_field_current_series(res, name, dt)
end

"""
    get_field_voltage_series(
            res::SimulationResults,
            name::String,
    )

Function to obtain the field voltage time series of a Dynamic Generator out of the DAE Solution.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
- `name::String` : Name to identify the specified device
"""
function get_field_voltage_series(res::SimulationResults, name::String; dt = nothing)
    return post_proc_field_voltage_series(res, name, dt)
end

"""
    show_states_initial_value(res::SimulationResults)

Function to print initial states.

# Arguments

- `res::SimulationResults` : Simulation Results object that contains the solution
"""
function show_states_initial_value(res::SimulationResults)
    bus_size = get_bus_count(res)
    system = get_system(res)
    x0_init = res.solution.u[1]
    println("Voltage Variables")
    println("====================")
    buses_sorted =
        sort(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
    for bus in buses_sorted
        name = PSY.get_name(bus)
        println(name)
        println("====================")
        bus_n = PSY.get_number(bus)
        bus_ix = get_bus_lookup(res)[bus_n]
        V_R = x0_init[bus_ix]
        V_I = x0_init[bus_ix + bus_size]
        Vm = sqrt(V_R^2 + V_I^2)
        θ = angle(V_R + V_I * 1im)
        print("Vm ", round(Vm, digits = 4), "\n")
        print("θ ", round(θ, digits = 4), "\n")
        println("====================")
    end
    println("====================")
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        println("Differential States")
        println(name)
        println("====================")
        global_index = get_global_index(res)[name]
        for s in states
            print(s, " ", round(x0_init[global_index[s]], digits = 4), "\n")
        end
        println("====================")
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        println("====================")
        println("Line Current States")
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            printed_name = "Line " * name
            println("====================")
            println(printed_name)
            global_index = get_global_index(res)[name]
            for (i, s) in enumerate(states)
                print(s, " ", round(x0_init[global_index[s]], digits = 5), "\n")
            end
            println("====================")
        end
    end
    return
end


text:
using GLAbstraction, ModernGL, GLWindow, GLFW

window 	= createwindow("GLPlot", 1920,1280, windowhints=[(GLFW.SAMPLES, 0)])
ocamera = OrthographicPixelCamera(window.inputs)

points = Vec3[Vec3(sin(i), i, 0) * 100f0 for i=0:0.1:2pi]

data = [
:vertex         	=> GLBuffer(Float32[1:4], 1), # simple quad
:index          	=> indexbuffer(GLuint[0, 1, 2, 2, 3, 0]),

:points 			=> Texture(points),
:linewidth 			=> 20f0,
:shadow 			=> 0.0f0,
:projection			=> ocamera.projection,
:view 				=> ocamera.view,
]

program = TemplateProgram(
	Pkg.dir("GLDraw","src","line.vert"), Pkg.dir("GLDraw","src","line.frag"),
	attributes=data, fragdatalocation=[(0, "fragment_color"),(1, "fragment_objectid")]
)
obj = instancedobject(data, length(points)-3, program)
prerender!(obj, enabletransparency)

glClearColor(1, 1, 1, 1.0)

while !GLFW.WindowShouldClose(window.nativewindow)
    yield()
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0,0,window.inputs[:framebuffer_size].value...)
    render(obj)

    GLFW.SwapBuffers(window.nativewindow)
	GLFW.PollEvents()
end
GLFW.Terminate()


text:
using Primes

function getmersenneprimes(n)
    t1 = time()
    count = 0
    i = 2
    while(n > count)
        if(isprime(i) && ismersenneprime(2^BigInt(i) - 1))
            println("M$i, cumulative time elapsed: $(time() - t1) seconds")
            count += 1
        end
        i += 1
    end
end

getmersenneprimes(50)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("individual.jl")

########################################
########################################
#########      POPULATION      #########
########################################
########################################
mutable struct Population
    pop::Array{Individual}
    Population() = new()
end

max_stars_repo_name: mahmoud-al-najar/neuroevolution

max_issues_repo_name: mahmoud-al-najar/neuroevolution

max_forks_repo_name: mahmoud-al-najar/neuroevolution


text:
module TabularTMazeUtils
using Curiosity
using GVFHordes
using SparseArrays
import ..TMazeCumulantSchedules
import ..GVFSRHordes
import ..FeatureCreator
import ..SRCreationUtils

using Distributions
using ...StatsBase
const TTMCS = TMazeCumulantSchedules

const NUM_DEMONS = 4
const NUM_ACTIONS = 4

DrifterDistractor(parsed) = begin
    c_dist = Uniform(parsed["constant_target"][1],parsed["constant_target"][2])
    c1,c2 = rand(c_dist,2)

    if "drifter" ∈ keys(parsed)
        TTMCS.DrifterDistractor(
            c1,
            c2,
            parsed["drifter"][1],
            parsed["drifter"][2],
            parsed["distractor"][1],
            parsed["distractor"][2])
    else

        TTMCS.DrifterDistractor(
            c1,
            c2,
            parsed["drifter_init"],
            parsed["drifter_std"],
            parsed["distractor_mean"],
            parsed["distractor_std"])
    end
end

function get_cumulant_schedule(parsed)
    sched = parsed["cumulant_schedule"]
    if parsed["cumulant_schedule"] == "DrifterDistractor"
        DrifterDistractor(parsed)
    elseif parsed["cumulant_schedule"] == "Constant"
        if parsed["cumulant"] isa Number
            TTMCS.Constant(parsed["cumulant"])
        else
            TTMCS.Constant(parsed["cumulant"]...)
        end

    else
        throw("$(sched) Not Implemented")
    end
end

function pseudoterm(;kwargs...)
    obs = kwargs[:state_tp1]
    term = false
    if obs[1] == 1 || obs[1] == 5 || obs[1] == 17 || obs[1] == 21
        term = true
    end
    return term
end

function valid_state_mask()
    """
    return a mask of valid states that is 9x7
    """
    world = [["G1", "0", "0", "0", "0", "0", "G3"],
             ["1", "0", "0", "0", "0", "0", "1"],
             ["1", "1", "1", "1", "1", "1", "1"],
             ["1", "0", "0", "1", "0", "0", "1"],
             ["G2", "0", "0", "1", "0", "0", "G4"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"]]

    world = permutedims(hcat(world...))
    valid_states = findall(x-> x!="0", world)
    return valid_states
end

struct TTMazeStateActionCumulant <: GVFParamFuncs.AbstractCumulant
    state_num::Int
    action::Int
end

function Base.get(cumulant::TTMazeStateActionCumulant; kwargs...)
    state = kwargs[:constructed_state_t]
    action = kwargs[:action_t]
    if state.nzind[1] == cumulant.state_num && action == cumulant.action
        return 1
    else
        return 0
    end
end

## Create Reward feature Projector
struct StateActionFeatures <: FeatureCreator
    num_states::Int
    num_actions::Int
end

function project_features(fc::StateActionFeatures, state, action, state_tp1)
    new_s = spzeros(fc.num_actions * fc.num_states)
    ind = (state[1] -1)* fc.num_actions + action
    new_s[convert(Int,ind)] = 1
    return new_s
end

(FP::StateActionFeatures)(state,action,state_tp1) = project_features(FP, state,action,state_tp1)
Base.size(FP::StateActionFeatures) = FP.num_states * FP.num_actions

function make_behaviour_gvf(discount, state_constructor_func, learner, exploration_strategy)
    function b_π(state_constructor_func, learner, exploration_strategy; kwargs...)
        s = state_constructor_func(kwargs[:state_t])
        preds = learner(s)
        return exploration_strategy(preds)[kwargs[:action_t]]
    end
    GVF_policy = GVFParamFuncs.FunctionalPolicy((;kwargs...) -> b_π(state_constructor_func, learner, exploration_strategy; kwargs...))
    BehaviourGVF = GVF(GVFParamFuncs.RewardCumulant(), GVFParamFuncs.StateTerminationDiscount(discount, pseudoterm), GVF_policy)
end
function demon_target_policy(gvf_i; kwargs...)

    # state = convert(Int,observation[1])
    state = convert(Int,kwargs[:state_t][1])
    action = kwargs[:action_t]

    policy_1 =  [1 0 0 0 0 0 3;
                 1 0 0 0 0 0 3;
                 1 4 4 4 4 4 4;
                 1 0 0 1 0 0 1;
                 1 0 0 1 0 0 1;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0]

    policy_2  = [3 0 0 0 0 0 3;
                 3 0 0 0 0 0 3;
                 3 4 4 4 4 4 4;
                 3 0 0 1 0 0 1;
                 3 0 0 1 0 0 1;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0]
    policy_3 = [3 0 0 0 0 0 1;
                 3 0 0 0 0 0 1;
                 2 2 2 2 2 2 1;
                 1 0 0 1 0 0 1;
                 1 0 0 1 0 0 1;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0]
    policy_4 = [3 0 0 0 0 0 3;
                 3 0 0 0 0 0 3;
                 2 2 2 2 2 2 3;
                 1 0 0 1 0 0 3;
                 1 0 0 1 0 0 3;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0;
                 0 0 0 1 0 0 0]
    gvfs = [policy_1,policy_2,policy_3, policy_4]
    mask = valid_state_mask()

    action_prob = if gvfs[gvf_i][mask][state] == action
        1.0
    else
        0.0
    end
    return action_prob
end

struct GoalPolicy <: GVFHordes.GVFParamFuncs.AbstractPolicy
    goal::Int
end
(π::GoalPolicy)(s) = begin
    sample(Weights([demon_target_policy(π.goal, state_t=s, action_t=a) for a ∈ 1:4]))
end

function Base.get(π::GoalPolicy; state_t, action_t, kwargs...)
    s = state_t[1]
    demon_target_policy(π.goal; state_t=s, action_t)
end

function get_true_values(env::Curiosity.TabularTMaze, eval_set)
    copy_eval_est = deepcopy(eval_set)
    num_gvfs = 4
    goal_cumulants = TTMCS.get_cumulant_eval_values(env.cumulant_schedule)
    for i in 1:num_gvfs
        copy_eval_est[i, :] .*= goal_cumulants[i]
    end
    return copy_eval_est
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""Bitwise XOR delta. Elements include A are XORed with the previous one. The
first element is left unchanged.
E.g. [0b0011,0b0010] -> [0b0011,0b0001] """
function xor_delta!(A::Array{T,1}) where {T<:Unsigned}
    a = A[1]
    @inbounds for i in 2:length(A)
        b = A[i]
        A[i] = a ⊻ b
        a = b
   end
end

"""Undo bitwise XOR delta. Elements include A are XORed again to reverse xor_delta.
E.g. [0b0011,0b0001] -> [0b0011,0b0010] """
function unxor_delta!(A::Array{T,1}) where {T<:Unsigned}
    a = A[1]
    @inbounds for i in 2:length(A)
        b = A[i]
        a = a ⊻ b
        A[i] = a
   end
end

"""Bitwise XOR delta. Elements include A are XORed with the previous one. The
first element is left unchanged.
E.g. [0b0011,0b0010] -> [0b0011,0b0001]. """
function xor_delta(A::Array{T,1}) where {T<:Unsigned}
   B = copy(A)
   xor_delta!(B)
   return B
end

"""Undo bitwise XOR delta. Elements include A are XORed again to reverse xor_delta.
E.g. [0b0011,0b0001] -> [0b0011,0b0010] """
function unxor_delta(A::Array{T,1}) where {T<:Unsigned}
    B = copy(A)
    unxor_delta!(B)
    return B
end

"""Bitwise XOR delta. Elements include A are XORed with the previous one. The
first element is left unchanged.
E.g. [0b0011,0b0010] -> [0b0011,0b0001]. """
function xor_delta(::Type{UIntT},A::Array{T,1}) where {UIntT<:Unsigned,T<:AbstractFloat}
   B = reinterpret.(UIntT,A)
   xor_delta!(B)
   return reinterpret.(T,B)
end

"""Undo bitwise XOR delta. Elements include A are XORed again to reverse xor_delta.
E.g. [0b0011,0b0001] -> [0b0011,0b0010] """
function unxor_delta(::Type{UIntT},A::Array{T,1}) where {UIntT<:Unsigned,T<:AbstractFloat}
   B = reinterpret.(UIntT,A)
   unxor_delta!(B)
   return reinterpret.(T,B)
end

"""Bitwise XOR delta. Elements include A are XORed with the previous one. The
first element is left unchanged.
E.g. [0b0011,0b0010] -> [0b0011,0b0001]. """
xor_delta(A::Array{T,1}) where {T<:AbstractFloat} = xor_delta(whichUInt(T),A)

"""Undo bitwise XOR delta. Elements include A are XORed again to reverse xor_delta.
E.g. [0b0011,0b0001] -> [0b0011,0b0010] """
unxor_delta(A::Array{T,1}) where {T<:AbstractFloat} = unxor_delta(whichUInt(T),A)