text:
using Printf

# Total circles area: https://rosettacode.org/wiki/Total_circles_area

xc = [1.6417233788, -1.4944608174, 0.6110294452, 0.3844862411, -0.2495892950, 1.7813504266,
     -0.1985249206, -1.7011985145, -0.4319462812, 0.2178372997, -0.6294854565, 1.7952608455,
      1.4168575317, 1.4637371396, -0.5263668798, -1.2197352481, -0.1389358881, 1.5293954595,
     -0.5258728625, -0.1403562064, 0.8055826339, -0.6311979224, 1.4685857879, -0.6855727502,
      0.0152957411]
yc = [1.6121789534, 1.2077959613, -0.6907087527, 0.2923344616, -0.3832854473, 1.6178237031,
     -0.8343333301, -0.1263820964, 1.4104420482, -0.9499557344, -1.3078893852, 0.6281269104,
      1.0683357171, 0.9463877418, 1.7315156631, 0.9144146579, 0.1092805780, 0.0030278255,
      1.3782633069, 0.2437382535, -0.0482092025, 0.7184578971, -0.8347049536, 1.6465021616,
      0.0638919221]
r  = [0.0848270516, 1.1039549836, 0.9089162485, 0.2375743054, 1.0845181219, 0.8162655711,
      0.0538864941, 0.4776976918, 0.7886291537, 0.0357871187, 0.7653357688, 0.2727652452,
      1.1016025378, 1.1846214562, 1.4428514068, 1.0727263474, 0.7350208828, 1.2472867347,
      1.3495508831, 1.3804956588, 0.3327165165, 0.2491045282, 1.3670667538, 1.0593087096,
      0.9771215985]

# Size of my grid -- higher values => higher accuracy.
function main(xc::Vector{<:Real}, yc::Vector{<:Real}, r::Vector{<:Real}, ngrid::Integer=10000)
    r2 = r .* r
    ncircles = length(xc)

    # Compute the bounding box of the circles.
    xmin = minimum(xc .- r)
    xmax = maximum(xc .+ r)
    ymin = minimum(yc .- r)
    ymax = maximum(yc .+ r)
    # Keep a counter.
    inside = 0
    # For every point in my grid.
    for x in linspace(xmin, xmax, ngrid), y = linspace(ymin, ymax, ngrid)
        inside += any(r2 .> (x - xc) .^ 2 + (y - yc) .^ 2)
    end
    boxarea = (xmax - xmin) * (ymax - ymin)
    return boxarea * inside / ngrid ^ 2
end

println(@time main(xc, yc, r, 1000))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module NationalInsuranceCalculations
#
# This module calculates National Insurance for an individual.
# It's based mainly on the descriptions and examples in the excellent Melville's Taxation:
# Melville, Alan. 2019. Melville’s Taxation: Finance Act 2019. 25th ed. London. (and 2020 edition).
#
using BudgetConstraints #: BudgetConstraint, get_x_from_y
using Dates
using Dates: Date, now, TimeType, Year
using Parameters: @with_kw

using ScottishTaxBenefitModel
using .STBIncomes
using .Definitions
using .ModelHousehold: Person
using .STBParameters: NationalInsuranceSys
using .GeneralTaxComponents:
    RateBands,
    TaxResult,
    *,
    calctaxdue;

using .Utils:
    BC_SETTINGS,
    eq_nearest_p,
    get_if_set

using .Results:
    IndividualResult,
    NIResult

export
    calc_class1_secondary,
    calculate_national_insurance!

function calc_class1_secondary( gross :: Real, pers::Person, sys :: NationalInsuranceSys ) :: Real
    rates = copy( sys.secondary_class_1_rates )
    # FIXME parameterise this
    if pers.age <= 21 # or  age <= 25 and apprentice
        rates[2] = 0.0
    end
    tres = calctaxdue(
        taxable = gross, # get(pers.income,wages, 0.0)
        rates = rates,
        thresholds = sys.secondary_class_1_bands )
    tres.due
    ## TODO apprentiships
end

function make_one_net( data :: Dict, gross :: Real ) :: Real
    pers = data[:pers]
    sys  = data[:sys]
    # pers.income[wage] = gross
    ni = calc_class1_secondary( gross, pers, sys )
    return gross - ni
end

function make_gross_wage_bc( pers :: Person, sys :: NationalInsuranceSys ) :: BudgetConstraint
    data = Dict(
        :pers=>pers,
        :sys=>sys
    )
    # hourly wage is irrelevant
    return makebc( data, make_one_net, Utils.BC_SETTINGS)
end

"""
FIXME pass in AgeLimits record
"""
function calculate_national_insurance!(
    pres :: IndividualResult,
    pers :: Person,
    sys  :: NationalInsuranceSys )
    # employer's NI on any wages
    if pers.age < 16 || sys.abolished # must be 16+
        return
    end
    bc = make_gross_wage_bc( pers, sys )
    wage = isum(
        pres.income,
        sys.class_1_income )

    wage = max(0.0, wage)
    gross = gross_from_net( bc, wage )
    pres.ni.class_1_secondary = calc_class1_secondary( gross, pers, sys )
    @assert isapprox(gross - wage, pres.ni.class_1_secondary, atol=3 ) "gross $gross wage $wage pres.ni.class_1_secondary $(pres.ni.class_1_secondary)"
    pres.ni.assumed_gross_wage = gross

    # class 1 on any wages, se only on main ..
    if pers.age < sys.state_pension_age # FIXME pass in the age limit thing
        tres = calctaxdue(
            taxable = wage,
            rates = sys.primary_class_1_rates,
            thresholds = sys.primary_class_1_bands )
        pres.ni.class_1_primary = tres.due
        pres.ni.above_lower_earnings_limit = tres.end_band > 1

        if( pers.employment_status in [Full_time_Self_Employed, Part_time_Self_Employed])
            seinc = isum(
                pres.income,
                sys.class_4_income )
            # maybe? pers.principal_employment_type != An_Employee
            # FIXME do I need *any* check on whether someone is classed as SE, & not just se income present?
            if seinc > sys.class_2_threshold
                pres.ni.class_2 = sys.class_2_rate
            end
            pres.ni.class_4 = calctaxdue(
                taxable = seinc,
                rates = sys.class_4_rates,
                thresholds = sys.class_4_bands ).due
        end # self emp
    end

    # do something random for class 3

    # don't count employers NI here
    pres.income[NATIONAL_INSURANCE] =
        pres.ni.class_1_primary +
        pres.ni.class_2 +
        pres.ni.class_4

 end

function gross_from_net( bc :: BudgetConstraint, net :: Real )::Real
    return get_x_from_y( bc, net )
end

end # module

ext: jl

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
@testset "946.validate-stack-sequences.jl" begin
    @test validate_stack_sequences([1, 2, 3, 4, 5], [4, 5, 3, 2, 1])
    @test validate_stack_sequences([1, 2, 3, 4, 5], [4, 3, 5, 1, 2]) == false
    @test validate_stack_sequences(Int[], Int[])
end


text:
#####
##### Writing EDFs
#####

function _write(io::IO, value::T) where T<:Union{PatientID,RecordingID}
    nb = 0
    for f in fieldnames(T)
        x = getfield(value, f)
        if T <: RecordingID && f === :startdate
            nb += _write(io, "Startdate ")
        end
        nb += _write(io, x) + Base.write(io, 0x20)
    end
    return nb
end

_write(io::IO, value) = Base.write(io, string(value))
_write(io::IO, date::Date) = Base.write(io, uppercase(Dates.format(date, dateformat"dd-u-yyyy")))
_write(io::IO, ::Missing) = Base.write(io, "X")
_write(io::IO, x::AbstractFloat) = Base.write(io, string(isinteger(x) ? trunc(Int, x) : x))

function _write(io::IO, tal::AnnotationsList)
    nb = _write(io, tal.offset >= 0 ? '+' : '-') + _write(io, tal.offset)
    if tal.duration !== nothing
        nb += Base.write(io, 0x15) + _write(io, tal.duration)
    end
    nb += Base.write(io, 0x14)
    mark = position(io)
    join(io, tal.event, '\x14')
    nb += position(io) - mark
    nb += Base.write(io, 0x14, 0x0)
    return nb
end

function _write(io::IO, anno::RecordAnnotation)
    nb = _write(io, anno.offset >= 0 ? '+' : '-') +
         _write(io, anno.offset) +
         Base.write(io, 0x14, 0x14)
    mark = position(io)
    join(io, anno.event, '\x14')
    nb += position(io) - mark
    nb += Base.write(io, 0x0)
    for tal in anno.annotations
        nb += _write(io, tal)
    end
    while nb < anno.n_bytes
        nb += Base.write(io, 0x0)
    end
    @assert nb == anno.n_bytes
    return nb
end

function write_padded(io::IO, value, n::Integer)
    b = _write(io, value)
    @assert b <= n
    while b < n
        b += Base.write(io, 0x20)
    end
    return b
end

function write_header(io::IO, file::File)
    h = file.header
    has_anno = file.annotations !== nothing
    b = write_padded(io, h.version, 8) +
        write_padded(io, h.patient, 80) +
        write_padded(io, h.recording, 80) +
        _write(io, Dates.format(h.start, dateformat"dd\.mm\.yyHH\.MM\.SS")) +
        write_padded(io, h.nb_header, 8) +
        write_padded(io, h.continuous ? "EDF+C" : "EDF+D", 44) +
        write_padded(io, h.n_records, 8) +
        write_padded(io, h.duration, 8) +
        write_padded(io, h.n_signals + has_anno, 4)
    pads = [16, 80, 8, 8, 8, 8, 8, 80, 8]
    av = Any["EDF Annotations", "", "", -1, 1, -32768, 32767, ""]
    has_anno && push!(av, div(first(file.annotations).n_bytes, 2))
    for (i, w) in zip(1:fieldcount(Signal)-1, pads)
        for s in file.signals
            b += write_padded(io, getfield(s, i), w)
        end
        if has_anno
            b += write_padded(io, av[i], w)
        end
    end
    ns = 32 * (length(file.signals) + has_anno)
    for _ = 1:ns
        b += Base.write(io, 0x20)
    end
    return b
end

function write_data(io::IO, file::File)
    b = 0
    for i = 1:file.header.n_records
        for signal in file.signals
            n = signal.n_samples
            s = (i - 1) * n
            stop = min(s + n, length(signal.samples))
            b += Base.write(io, view(signal.samples, s+1:stop))
        end
        if file.annotations !== nothing
            b += _write(io, file.annotations[i])
        end
    end
    return b
end

"""
    EDF.write(io::IO, edf::EDF.File)
    EDF.write(path::AbstractString, edf::EDF.File)

Write the given `EDF.File` object to the given stream or file and return the number of
bytes written.
"""
write(io::IO, file::File) = write_header(io, file) + write_data(io, file)
write(path::AbstractString, edf::File) = open(io->write(io, edf), path, "w")


text:
function assemble_traction_force_component_linear_form!(
    systemrhs,
    tractionfunc,
    basis,
    facequads,
    mesh,
    onboundary,
    component
)

    dim = dimension(mesh)
    facemidpoints = reference_face_midpoints()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    ncells = number_of_cells(mesh)

    for cellid in 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid = 1:nfaces
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )
                if cellsign == 0 || cellsign == +1
                    assemble_face_traction_force_component_linear_form!(
                        systemrhs,
                        tractionfunc,
                        basis,
                        facequads,
                        component,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                    )
                end
                if cellsign == 0 || cellsign == -1
                    assemble_face_traction_force_component_linear_form!(
                        systemrhs,
                        tractionfunc,
                        basis,
                        facequads,
                        component,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                    )
                end
            end
        end
    end
end

function assemble_face_traction_force_component_linear_form!(
    systemrhs,
    tractionfunc,
    basis,
    facequads,
    component,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
)
    dim = dimension(mesh)
    cellmap = cell_map(mesh, cellsign, cellid)
    rhs = component_linear_form(
        tractionfunc,
        basis,
        facequads[cellsign, faceid, cellid],
        component,
        cellmap,
        dim,
        facedetjac[faceid],
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test
using JuliaFEM.Preprocess
using SparseArrays, Test

# test tet4 + volume load

X = Dict(1 => [2.0, 3.0, 4.0], 2 => [6.0, 3.0, 2.0],
         3 => [2.0, 5.0, 1.0], 4 => [4.0, 3.0, 6.0])
element1 = Element(Tet4, (1, 2, 3, 4))
element2 = Element(Tri3, (1, 2, 3))
update!((element1, element2), "geometry", X)
update!(element1, "youngs modulus", 96.0)
update!(element1, "poissons ratio", 1/3)
update!(element1, "displacement load 3", 784.0/110.0)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
update!(element2, "displacement 3", 0.0)
problem1 = Problem(Elasticity, "tetra", 3)
problem2 = Problem(Dirichlet, "bc", 3, "displacement")
add_element!(problem1, element1)
add_element!(problem2, element2)
analysis = Analysis(Linear)
add_problems!(analysis, problem1, problem2)
run!(analysis)
u4 = problem1("displacement", 0.0)[4]
u4_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
@test isapprox(u4, u4_expected)

# tet4 + surface load

element3 = Element(Tri3, (4, 3, 2))
update!(element3, "geometry", X)
update!(element3, "surface pressure", -96.0)
update!(element1, "displacement load 3", 0.0)
add_element!(problem1, element3)
run!(analysis)
u4 = problem1("displacement", 0.0)[4]
u4_expected = [-17.0/14.0, -27.0/14.0, 1.0]
@debug("displacement at node 4", u4, u4_expected)

#= TODO: Fix test. Make linear perturbation solver.
@testset "test tet4 + buckling" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [2.0, 3.0, 4.0],
        2 => [6.0, 3.0, 2.0],
        3 => [2.0, 5.0, 1.0],
        4 => [4.0, 3.0, 6.0])
    u = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0, 0.0],
        2 => [0.0, 0.0, 0.0],
        3 => [0.0, 0.0, 0.0],
        4 => [-0.25, -0.25, -0.25])
    e1 = Element(Tet4, [1, 2, 3, 4])
    update!(e1, "geometry", X)
    update!(e1, "displacement", u)
    update!(e1, "youngs modulus", 96.0)
    update!(e1, "poissons ratio", 1/3)
    p1 = Problem(Elasticity, "tetra", 3)
    p1.properties.formulation = :continuum_buckling
    push!(p1, e1)
    assemble!(p1, 0.0)
    free_dofs = [10, 11, 12]
    Km = sparse(p1.assembly.K)[free_dofs, free_dofs]
    Kg = sparse(p1.assembly.Kg)[free_dofs, free_dofs]
    dump(full(Km))
    dump(full(Kg))
    la = sort(eigs(Km, -Kg)[1])
    la_expected = [1.0, 4.0]
    @info("la = $la")
    @info("la_expected = $(la_expected)")
    @test isapprox(la, la_expected)
end
=#


text:
@testset "Duals: Test Basic" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))

    @variable(model, x[1:2] in POI.Parameter.(ones(2) .* 4.0))
    @variable(model, y[1:6])

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

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

    JuMP.optimize!(model)

    @test 5 / 3 ≈ JuMP.dual(ctr1) atol = 1e-3
    @test [-35 / 3, 0.0] ≈ MOI.get.(model, POI.ParameterDual(), 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(model) atol = 1e-3
end

@testset "Duals: Test Multiple Parameters" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))

    @variable(model, x[1:6] in POI.Parameter.(ones(6) .* 4.0))
    @variable(model, y[1:6])

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

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

    JuMP.optimize!(model)

    @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] ≈
          MOI.get.(model, POI.ParameterDual(), 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(model) atol = 1e-3
end

@testset "Duals: Test LessThan" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(-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 MOI.get(model, POI.ParameterDual(), α) == -1.0

    MOI.set(model, POI.ParameterValue(), α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == -1.0
    @test MOI.get(model, POI.ParameterDual(), α) == -1.0
end

@testset "Duals: Test EqualTo" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(-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 MOI.get(model, POI.ParameterDual(), α) == -1.0

    MOI.set(model, POI.ParameterValue(), α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == -1.0
    @test MOI.get(model, POI.ParameterDual(), α) == -1.0
end

@testset "Duals: Test GreaterThan" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(1.0))
    MOI.set(model, POI.ParameterValue(), α, -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 MOI.get(model, POI.ParameterDual(), α) == 1.0

    MOI.set(model, POI.ParameterValue(), α, 2.0)
    JuMP.optimize!(model)
    @test JuMP.value(x) == 2.0
    @test JuMP.dual(cref) == 1.0
    @test MOI.get(model, POI.ParameterDual(), α) == 1.0
end

@testset "Duals: Test Multiple Parameters 2" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α[1:10] in POI.Parameter.(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 MOI.get(model, POI.ParameterDual(), α[3]) == 2.0
end

@testset "Duals: Test Mixing Params and vars 1" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α[1:5] in POI.Parameter.(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 MOI.get(model, POI.ParameterDual(), α[3]) == 2 / 5
end

@testset "Duals: Test mixing Params and vars 2" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α[1:5] in POI.Parameter.(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 MOI.get(model, POI.ParameterDual(), α[3]) == 2 / 5
end

@testset "Duals: Test mixing Params and vars 3" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α[1:5] in POI.Parameter.(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 MOI.get(model, POI.ParameterDual(), α[3]) == 2 / 5
end

@testset "Duals: Test add after solve" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(1.0))
    MOI.set(model, POI.ParameterValue(), α, -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 MOI.get(model, POI.ParameterDual(), α) == -1.0

    @variable(model, b in POI.Parameter(-2.0))
    cref = @constraint(model, x <= b)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -2.0
    @test JuMP.dual(cref) == -1.0
    @test MOI.get(model, POI.ParameterDual(), α) == 0.0
    @test MOI.get(model, POI.ParameterDual(), b) == -1.0
end

@testset "Duals: Test add ctr alaternative" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(-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 MOI.get(model, POI.ParameterDual(), α) == -1.0
end

@testset "Duals: Test deletion constraint" begin
    model = Model(() -> POI.Optimizer(GLPK.Optimizer()))
    @variable(model, α in POI.Parameter(-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 MOI.get(model, POI.ParameterDual(), α) == -1.0

    JuMP.delete(model, cref2)
    JuMP.optimize!(model)
    @test JuMP.value(x) == -0.5
    @test JuMP.dual(cref1) == -1.0
    @test MOI.get(model, POI.ParameterDual(), α) == -0.5
end


text:
"This abstract type needs to be implemented by all equations."
abstract type Equation end

"This abstract type needs to be implementred by all scenarios"
abstract type Scenario end

"""
This function initializes the discrete solution at t = 0 for an equation `eq` and a scenario `scenario`.
It evaluates get\\_initial\\_values for `eq` and `scenario` at the global coordinates of our nodes.
It then interpolates this on the `basis`.
"""
function interpolate_initial_dofs(eq::Equation, scenario::Scenario, celldofs, cell, basis)
    celldofs .= project_to_reference_basis(basis, get_ndofs(eq)) do px, py
        pxg, pyg = globalposition(cell, (px, py))
        get_initial_values(eq, scenario, (pxg, pyg))
    end
end

"This function returns the number of variables for equation `eq::Equation`"
function get_nvars(eq::Equation)
    get_ndofs(eq) - get_nparams(eq)
end

"This function returns the number of material parameters for equation `eq::Equation`."
function get_nparams(eq::Equation)
    0
end

"""

    get_variable_name(eq::Equation, variable_index)

Return the name of the variable with index `variable_index` for equation `eq::Equation`.
# Implementation
Custom equations need to define `get_variable_names(eq)` for this to work.
"""
function get_variable_name(eq::Equation, variable_index)
    get_variable_names(eq)[variable_index]
end

"""
    is_periodic_boundary(equation::Equation, scenario::Scenario)

Returns true if all boundary conditions should be periodic for equation `eq::Equation` and
scenario `scenario::Scenario`.
"""
function is_periodic_boundary(equation::Equation, scenario::Scenario)
    false
end

"""
    evaluate_boundary(eq::Equation, scenario::Equation, face, normalidx, dofsface, dofsfaceneigh)

Evaluate the boundary condition, defined for equation `eq` and `scenario`, on face `face` with `normalidx`.
Store dofs of ghost cell in `dofsfaceneigh`.
You can use inner (face)-dofs `dofsface` to construct to boundary.
"""
function evaluate_boundary(eq::Equation, scenario::Equation, face, normalidx, dofsface, dofsfaceneigh)
    error("evaluate_boundary called but not define.")
end

"""
    get_initial_values(eq::Equation, scenario::Scenario, global_position; t=0.0)

Return a vector of initial values for equation `eq` and scenario `scenario` at
`global_positon` and at time `t`.
Values other than `t = 0` only have to be supported if
`is_analytical_solution` returns true.
"""
function get_initial_values(eq::Equation, scenario::Scenario, global_position; t=0.0)
    error("Function get_initial_values not defined for equation/scenario")
end

"""
    is_analytical_solution(equation::Equation, scenario::Scenario)

Returns true if the initial condition is also an analytical solution
for equation `equation` and scenario `scenario`.
"""
function is_analytical_solution(equation::Equation, scenario::Scenario)
    false
end

"""
    evaluate_flux(eq::Equation, celldofs, cellflux)

Evaluate the flux for equation `eq` and dofs `celldofs`,
Stores them in `cellflux`.
"""
function evaluate_flux(eq::Equation, celldofs, cellflux)
    error("Function evaluate_flux not defined for equation/scenario")
end

#function evaluate_flux_with_source(eq::Equation, celldofs, cellflux, cell, basis)
#    error("Function evaluate_flux not defined for equation/scenario")
#end

function evaluate_flux_with_stiffinclusion(eq::Equation, scenario::Scenario, celldofs, cellflux, global_position)
    error("Function evaluate_flux_with_stiffinclusion not defined for equation/scenario")
end
"""
    max_eigenval(eq::Equation, celldata, normalidx)

Return the maximum eigenvalue in direction of `normalidx` for
`celldata`.

"""
function max_eigenval(eq::Equation, celldata, normalidx)
    error("Function max_eigenval not defined for equation/scenario")
end

"""
    declare_dofs(eq, dof_names, num_params=0)

Generate boilerplate code for equation type `eq` which has
degrees of freedom named `dof_names` (should be iterable of symbols)
and `num_params` parameters.
Generates

    get_variable_names(eq)
    get_ndofs(eq)
    get_nvars(eq)
    get_nparams(eq)

Also generates a struct called [typename]Shortcuts.

# Example

    struct Advection <: Equation end
    @declare_dofs Advection [:ρ_1, :ρ_2, :ρ_3]

Generates:

    get_variable_names(eq::Advection) = [:ρ_1, :ρ_2, :ρ_3]
    get_ndofs(eq::Advection) = 3
    get_nvars(eq::Advection) = 3
    get_nparams(eq::Advection) = 0

You can use the struct in this way:

    s = AdvectionShortcuts()
    @assert s.ρ_1 == 1
    @assert s.ρ_2 == 2
    @assert s.ρ_3 == 3
"""
macro declare_dofs(eq, dof_names, num_params=0)
    local expressions::Array{Expr, 1} = []
    local shortcut_struct = Symbol(eq, "Shortcuts")
    push!(expressions, quote
        function $(esc(Symbol("get_variable_names")))(eq::$eq)
            $dof_names
        end

        function $(esc(Symbol("get_ndofs")))(eq::$eq)
            $(length(dof_names.args))
        end

        function $(esc(Symbol("get_nvars")))(eq::$eq)
            $(length(dof_names.args) - num_params)
        end

        function $(esc(Symbol("get_nparams")))(eq::$eq)
            $num_params
        end

        struct $(esc(shortcut_struct))

        end

        @inline function Base.getproperty(x::$(esc(shortcut_struct)), name::Symbol)
            get_variable_shortcut(x, Val(name))
        end
    end
    )

    for (i, name) in enumerate(dof_names.args)
        push!(expressions, quote
            @inline function $(esc(Symbol("get_variable_shortcut")))(eq::$(esc(shortcut_struct)), ::Val{$name})
                return $i
            end

            @inline function $(esc(Symbol("get_variable_shortcut")))(eq::$eq, ::Val{$name})
                return $i
            end
        end)
        end

    quote $(expressions...) end
end

include("equations/advection.jl")
include("equations/acoustic.jl")
include("equations/euler.jl")
include("equations/elastic.jl")

"""
    make_equation(config::Configuration)

Returns equation object for configuration `config`.
"""
function make_equation(config::Configuration)
    if config.equation_name == "advection"
        return Advection()
    elseif config.equation_name == "acoustic"
        return Acoustic()
    elseif config.equation_name == "euler"
        return Euler()
    elseif config.equation_name == "elastic"
        return Elastic()
    else
        error(string("Unknown equation name: ", config.equation_name))
    end
end

"""
    make_scenario(config::Configuration)

Returns scenario object for configuration `config`.
"""
function make_scenario(config::Configuration)
    if config.scenario_name == "planar_waves"
        return PlanarWaves()
    elseif config.scenario_name == "sibson"
        return Sibson()
    elseif config.scenario_name == "gaussian_wave"
        return GaussianWave()
    elseif config.scenario_name == "gaussian_wave_euler"
        return GaussianWaveeuler()
    elseif config.scenario_name == "sod_shock_tube"
        return Sod_shock_tube()
    elseif config.scenario_name == "elastic_waves"
        return Elastic_waves()
    elseif config.scenario_name == "gaussian_point"
        return GaussianPoint()
    elseif config.scenario_name == "elastic_source"
        return Elastic_source()
    elseif config.scenario_name == "stiff_inclusion"
        return StiffInclusion()
    else
        error(string("Unknown scenario name: ", scenario_name))
    end
end


text:
function lagrange_interpolant(R, xs, i)
    n = length(xs)
    x = AbstractAlgebra.gen(R)

    f = R(prod(
        j == i ? 1 : (x - xs[j])
        for j in 1:n
    ))

    c = prod(
        j == i ? 1 : (xs[i] - xs[j])
        for j in 1:n
    )

    AbstractAlgebra.map_coefficients(t -> t // c, f)
end

# univariate
# returns :: AbstractAlgebra.Generic.Poly
function interpolate_polynomial(R, xs, ys)
    f = sum(
        ys[i] * lagrange_interpolant(R, xs, i)
        for i in 1:length(xs)
    )
end

"""
    For given polynomial g and m computes the rational fraction r//t,
    s.t. r//t approximates g modulo m

        r == t*g (mod m)

    and
        deg(r) + deg(t) <= deg(m)

    Throws if ??? TODO
"""
function polynomial_reconstruction(g, m)
    # The function implements Algorithm 12.9 from
    #   "Modern Computer Algebra", second edition,
    #   Joachim von zur Gathen, Jürgen Gerhard

    n = AbstractAlgebra.degree(m)
    k = div(n, 2)
    FF = AbstractAlgebra.parent(g)

    U = (FF(1), FF(0), m)
    V = (FF(0), FF(1), g)
    while AbstractAlgebra.degree(V[3]) >= k
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T
    end

    r, t = V[3], V[2]

    if gcd(r, t) == FF(1) && AbstractAlgebra.degree(t) <= n - k
        return r // t
    end

    throw(DomainError(
        g, "reconstruction of $g (mod $m) does not exist"
    ))
end

"""
    Interpolates rational function in the given ring R with
    the given interpolation points xs and ys by
    interpolating as a polynomial first
    and reconstructing into a polynomial fraction after
"""
function interpolate_rational_function(R, xs, ys)
    # The function implements algorithm derived from Corollary 5.20 of
    #   "Modern Computer Algebra", second edition,
    #   Joachim von zur Gathen, Jürgen Gerhard

    x = AbstractAlgebra.gen(R)
    g = Nemo.interpolate(R, xs, ys)

    # TODO: this computation can be done faster
    m = prod(x - xs[i] for i in 1:length(xs))

    polynomial_reconstruction(g, m)
end

# TODO: why is it still here?
function interpolate_multivariate_rational_function(MR, R, xs, ys)

    f = interpolate_rational_function(R, xs, ys)

    f
end

"""
    Returns reversed vector of base-adic representation of n
"""
function decompose_by_degrees(n, base, len)
    powers = [base^(i-1) for i in 1:len]
    exps = zeros(Int, len)
    for i in 1:len
        exps[len - i + 1] = div(n, powers[len - i + 1])
        n = n - exps[len - i + 1] * powers[len - i + 1]
    end
    exps
end

# TODO : not used
function random_linear_shift(ground_ring, n)
    functors = []
    inverses = []

    for i in 1:n
        a, b = ground_ring(rand(1:10000)), ground_ring(rand(1:10000))
        f  = ( x -> x + b )
        fi = ( x -> x - b )

        push!(functors, f)
        push!(inverses, fi)
    end

    return functors, inverses
end

"""
    Given a callable functor f tries to find degrees of univariate polynomials
    A, B approximating f as a rational function f = A // B
"""
# TODO: not used
function predict_degrees(f)
    # staring from small amount of points..
    n = 8

    # = deg(A), deg(B) initially
    predicted_degrees = (n, n)

    polyring, = AbstractAlgebra.PolynomialRing(ground, "u")

    # "safety condition" to assume f is interpolated correctly indeed
    # failes on first iteration
    while 4*(sum(predicted_degrees) + 2) > n
        xs = [ rand(f) for _ in 1:n  ]
        ys = f.(xs)

        interpolated = interpolate_rational_function(polyring, xs, ys)

        predicted_degrees = map(degree, (numerator(interpolated), denominator(interpolated)))

        @debug "for $n points degrees are " predicted_degrees

        n = n * 2
    end

    predicted_degrees
end

"""
    Same as the function above, key difference is that here
    f returns an array of results on each evaluation
"""
# TODO: not used
function simultaneous_predict_degrees(f)
    # staring from small amount of points..
    n = 8

    polyring = base_ring(f)
    nfrom, nto = arity(f)

    predicted_degrees = [ (n, n) for _ in 1:nto ]
    already_interpolated = 0

    while already_interpolated < nto

        4*(sum(values(predicted_degrees)) + 2) > n
        xs = [ rand(f) for _ in 1:n  ]
        ys = f.(xs)

        for (j, deg) in enumerate(predicted_degrees)
            if 4*(sum(deg) + 2) < n
                continue
            end
            yparticular = [y[j] for y in ys]
            interpolated = interpolate_rational_function(polyring, xs, yparticular)

            predicted_degrees[j] = map(degree, (numerator(interpolated), denominator(interpolated)))

        end

        @debug "for $n points degrees are " predicted_degrees

        n = n * 2
    end

    predicted_degrees
end

"""

"""
# TODO: not used
function backward_kronecker(
                     f,
                     target_ring,
                     maxexp) where {T}

    polybuilder = MPolyBuildCtx(target_ring)
    nvariables = length(gens(target_ring))

    # TODO: fix 0:degree(f)
    for (e, c) in zip(0:degree(f), coefficients(f))
        if iszero(c)
            continue
        end
        push_term!(
            polybuilder,
            c,
            decompose_by_degrees(e, maxexp, nvariables)
        )
    end

    finish(polybuilder)
end

function decompose_by_good_degrees(e, prods, n)
    ans = zeros(Int, n)

    for i in 1:n
        ans[n - i + 1] = div(e, prods[n - i + 1])
        e = mod(e, prods[n - i + 1])
    end

    ans
end

function backward_good_kronecker(
                     f,
                     target_ring,
                     exps) where {T}

    prods = deepcopy(exps)
    prods[1] = 1
    for i in 2:length(prods)
        prods[i] = prods[i - 1] * (exps[i - 1] + 1)
    end

    polybuilder = MPolyBuildCtx(target_ring)
    nvariables = length(gens(target_ring))

    for (e, c) in zip(0:degree(f), coefficients(f))
        if iszero(c)
            continue
        end
        push_term!(
            polybuilder,
            c,
            decompose_by_good_degrees(e, prods, nvariables)
        )
    end

    ans = finish(polybuilder)
    ans
end

function backward_kronecker(
                  f::AbstractAlgebra.Generic.Frac{T},
                  target_ring,
                  maxexp) where {T}

    backward_kronecker(numerator(f), target_ring, maxexp) //
        backward_kronecker(denominator(f), target_ring, maxexp)
end

function backward_good_kronecker(
                  f::AbstractAlgebra.Generic.Frac{T},
                  target_ring,
                  maxexp) where {T}

    backward_good_kronecker(numerator(f), target_ring, maxexp) //
        backward_good_kronecker(denominator(f), target_ring, maxexp)
end

function generate_kronecker_points(ground, maxexp, nvariables)
    npoints = (maxexp^nvariables + 2) * 2
    xs = [rand(ground) for _ in 1:npoints]
    points = [
        [ ground(j)^i for i in [ (maxexp)^k for k in 0:(nvariables - 1) ] ]
          for j in xs
    ]
    xs, points
end

function generate_good_kronecker_points(ground, exps, nvariables)
    npoints = prod( map(x -> x + 1, exps) ) * 2

    prods = deepcopy(exps)
    # univariate case? =(
    prods[1] = 1
    for i in 2:length(prods)
        prods[i] = prods[i - 1] * (exps[i - 1] + 1)
    end

    xs = Set([rand(ground) for _ in 1:npoints])
    while length(xs) < npoints
        xs = union!(xs, Set([rand(ground) for _ in 1:(npoints-length(xs))]))
    end
    xs = collect(xs)

    points = [
        [ p^prods[k] for k in 1:nvariables ]
        for p in xs
    ]

    xs, points
end


text:
@testset "Randgraphs" begin
    r1 = SimpleGraph(10, 20)
    r2 = SimpleDiGraph(5, 10)

    @test nv(r1) == 10
    @test ne(r1) == 20
    @test nv(r2) == 5
    @test ne(r2) == 10
    @test eltype(r1) == Int
    @test eltype(r2) == Int

    @test eltype(Graph(0x5, 0x2)) == eltype(Graph(0x5, 2)) == UInt8
    for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt, Int]
        @test eltype(Graph{T}(5, 2)) == T
        @test eltype(DiGraph{T}(5, 2)) == T
    end

    @test SimpleGraph(10, 20, seed=3) == SimpleGraph(10, 20, seed=3)
    @test SimpleDiGraph(10, 20, seed=3) == SimpleDiGraph(10, 20, seed=3)
    @test SimpleGraph(10, 20, seed=3) == erdos_renyi(10, 20, seed=3)
    @test ne(Graph(10, 40, seed=3)) == 40
    @test ne(DiGraph(10, 80, seed=3)) == 80

    er = erdos_renyi(10, 0.5)
    @test nv(er) == 10
    @test is_directed(er) == false
    er = erdos_renyi(10, 0.5, is_directed=true)
    @test nv(er) == 10
    @test is_directed(er) == true

    er = erdos_renyi(10, 0.5, seed=17)
    @test nv(er) == 10
    @test is_directed(er) == false

    cl = expected_degree_graph(zeros(10), seed=17)
    @test nv(cl) == 10
    @test ne(cl) == 0
    @test is_directed(cl) == false

    cl = expected_degree_graph([3, 2, 1, 2], seed=17)
    @test nv(cl) == 4
    @test is_directed(cl) == false

    cl = expected_degree_graph(fill(99, 100), seed=17)
    @test nv(cl) == 100
    @test all(degree(cl) .> 90)

    ws = watts_strogatz(10, 4, 0.2)
    @test nv(ws) == 10
    @test ne(ws) == 20
    @test is_directed(ws) == false

    ws = watts_strogatz(10, 4, 0.2, is_directed=true)
    @test nv(ws) == 10
    @test ne(ws) == 20
    @test is_directed(ws) == true

    ba = barabasi_albert(10, 2)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, 2)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 4, 2)
    @test nv(ba) == 10
    @test ne(ba) == 12
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 17
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 17
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 4, 2, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == false

    ba = barabasi_albert(10, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 16
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 4, 2, is_directed=true)
    @test nv(ba) == 10
    @test ne(ba) == 12
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 2, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 18
    @test is_directed(ba) == true

    ba = barabasi_albert(10, 4, 2, is_directed=true, complete=true)
    @test nv(ba) == 10
    @test ne(ba) == 24
    @test is_directed(ba) == true

    fm = static_fitness_model(20, rand(10))
    @test nv(fm) == 10
    @test ne(fm) == 20
    @test is_directed(fm) == false

    fm = static_fitness_model(20, rand(10), rand(10))
    @test nv(fm) == 10
    @test ne(fm) == 20
    @test is_directed(fm) == true

    sf = static_scale_free(10, 20, 2.0)
    @test nv(sf) == 10
    @test ne(sf) == 20
    @test is_directed(sf) == false

    sf = static_scale_free(10, 20, 2.0, 2.0)
    @test nv(sf) == 10
    @test ne(sf) == 20
    @test is_directed(sf) == true

    rr = random_regular_graph(5, 0)
    @test nv(rr) == 5
    @test ne(rr) == 0
    @test is_directed(rr) == false

    rd = random_regular_digraph(10, 0)
    @test nv(rd) == 10
    @test ne(rd) == 0
    @test is_directed(rd)

    rr = random_regular_graph(6, 3, seed=1)
    @test nv(rr) == 6
    @test ne(rr) == 9
    @test is_directed(rr) == false

    rr = random_regular_graph(1000, 50)
    @test nv(rr) == 1000
    @test ne(rr) == 25000
    @test is_directed(rr) == false
    for v in vertices(rr)
        @test degree(rr, v) == 50
    end

    rr = random_configuration_model(10, repeat([2,4], 5), seed=3)
    @test nv(rr) == 10
    @test ne(rr) == 15
    @test is_directed(rr) == false
    num2 = 0; num4 = 0
    for v in vertices(rr)
        d = degree(rr, v)
        @test  d == 2 || d == 4
        d == 2 ? num2 += 1 : num4 += 1
    end
    @test num4 == 5
    @test num2 == 5

    rr = random_configuration_model(1000, zeros(Int, 1000))
    @test nv(rr) == 1000
    @test ne(rr) == 0
    @test is_directed(rr) == false

    rr = random_configuration_model(3, [2,2,2], check_graphical=true)
    @test nv(rr) == 3
    @test ne(rr) == 3
    @test is_directed(rr) == false

    rd = random_regular_digraph(1000, 4)
    @test nv(rd) == 1000
    @test ne(rd) == 4000
    @test is_directed(rd)
    outdegree_rd = @inferred(outdegree(rd))
    @test all(outdegree_rd .== outdegree_rd[1])

    rd = random_regular_digraph(1000, 4, dir=:in)
    @test nv(rd) == 1000
    @test ne(rd) == 4000
    @test is_directed(rd)
    indegree_rd = @inferred(indegree(rd))
    @test all(indegree_rd .== indegree_rd[1])

    rr = random_regular_graph(10, 8, seed=4)
    @test nv(rr) == 10
    @test ne(rr) == 40
    @test is_directed(rr) == false
    for v in vertices(rr)
        @test degree(rr, v) == 8
    end

    rd = random_regular_digraph(10, 8, dir=:out, seed=4)
    @test nv(rd) == 10
    @test ne(rd) == 80
    @test is_directed(rd)

    rt = random_tournament_digraph(10)
    @test nv(rt) == 10
    @test ne(rt) == 45
    @test is_directed(rt)
    @test all(degree(rt) .== 9)
    Edges = edges(rt)
    for i = 1:10, j = 1:10
        if i != j
            edge = Edge(i, j)
            @test xor(edge ∈ Edges, reverse(edge) ∈ Edges)
        end
    end

    g = stochastic_block_model(2., 3., [100,100])
    @test  4.0 < mean(degree(g)) < 6.0
    g = stochastic_block_model(3., 4., [100,100,100])
    @test  10.0 < mean(degree(g)) < 12.0

    function generate_nbp_sbm(numedges, sizes)
        density = 1
        between = density * 0.90
        intra = density * -0.005
        noise = density * 0.00501
        sbm = nearbipartiteSBM(sizes, between, intra, noise)
        edgestream = make_edgestream(sbm)
        g = SimpleGraph(sum(sizes), numedges, edgestream)
        return sbm, g
    end

    function test_sbm(sbm, bp)
        @test sum(sbm.affinities) != NaN
        @test all(sbm.affinities .> 0)
        @test sum(sbm.affinities) != 0
        @test all(bp .>= 0)
        @test all(bp .!= NaN)
    end

    numedges = 100
    sizes = [10, 10, 10, 10]

    n = sum(sizes)
    sbm, g = generate_nbp_sbm(numedges, sizes)
    @test ne(g) >= 0.9numedges
    bc = blockcounts(sbm, g)
    bp = blockfractions(sbm, g) ./ (sizes * sizes')
    ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
    test_sbm(sbm, bp)
    @test norm(collect(ratios)) < 0.25

    sizes = [200, 200, 100]
    internaldeg = 15
    externaldeg = 6
    internalp = Float64[internaldeg / i for i in sizes]
    externalp = externaldeg / sum(sizes)
    numedges = internaldeg + externaldeg #+ sum(externaldeg.*sizes[2:end])
    numedges *= div(sum(sizes), 2)
    sbm = StochasticBlockModel(internalp, externalp, sizes)
    g = SimpleGraph(sum(sizes), numedges, sbm)
    @test ne(g) >= 0.9numedges
    @test ne(g) <= numedges
    @test nv(g) == sum(sizes)
    bc = blockcounts(sbm, g)
    bp = blockfractions(sbm, g) ./ (sizes * sizes')
    test_sbm(sbm, bp)
    ratios = bp ./ (sbm.affinities ./ sum(sbm.affinities))
    @test norm(collect(ratios)) < 0.25

    # check that average degree is not too high
    # factor of two is cushion for random process
    @test mean(degree(g)) <= 4 // 2 * numedges / sum(sizes)
    # check that the internal degrees are higher than the external degrees
    # 5//4 is cushion for random process.
    @test all(sum(bc - diagm(0 => diag(bc)), dims=1) .<= 5 // 4 .* diag(bc))

    sbm2 = StochasticBlockModel(0.5 * ones(4), 0.3, 10 * ones(Int, 4))
    sbm  = StochasticBlockModel(0.5, 0.3, 10, 4)
    @test sbm == sbm2
    sbm.affinities[1,1] = 0
    @test sbm != sbm2

    kg = @inferred kronecker(5, 5)
    @test nv(kg) == 32
    @test is_directed(kg)

    g = @inferred(dorogovtsev_mendes(10))
    @test nv(g) == 10 && ne(g) == 17
    g = dorogovtsev_mendes(11)
    @test nv(g) == 11 && ne(g) == 19
    @test δ(g) == 2
    g = dorogovtsev_mendes(3)
    @test nv(g) == 3 && ne(g) == 3
    # testing domain errors
    @test_throws DomainError dorogovtsev_mendes(2)
    @test_throws DomainError dorogovtsev_mendes(-1)
end


text:
# Implement code for proposal evotech/complexity/cce/'Robustness evolvability compleity proposal 8_25_20.docx'.
# Includes evolving multiple chromosomes that correspond to every possible 3-input 1-output phenotype.
# In other words, evolves 3-input 1-output chromosomes (genotypes) whose output is every possible 8-bit UInt.
# There are
export run_geno_pheno_evolution, gp_result, gp_result_to_tuple, run_gp_evolve!
using DataFrames
using CSV
using Statistics
using Distributions
geno_pheno_result_type = Main.CGP.geno_pheno_result_type
#context = construct_contexts(numinputs)[numinputs]
#p = Parameters(numinputs, numoutputs, nodearity, numinteriors, numlevelsback)
#iterations = 20
#print_parameters(p)
using CSV
using Dates

function run_geno_pheno_evolution( iterations::Int64, numinputs::IntRange, numoutputs::IntRange,
    numinteriors::IntRange, goallistlength::IntRange, maxsteps::IntRange,
    levelsback::IntRange, csvfile::String;
    base::Float64=2.0, allgoals::Bool=false, active_only::Bool=false, gl_repetitions::IntRange=1)
  println(default_funcs(2))
  (df,ttime) =
    @timed run_geno_pheno( iterations, numinputs, numoutputs, numinteriors, goallistlength, maxsteps, levelsback,
        base=base, allgoals=allgoals, active_only=active_only, gl_repetitions=gl_repetitions )
  println("run time in minutes: ",ttime/60)
  println("cor( df.logsteps, df.complex): ",cor( df.logsteps, df.complex))
  println("cor( df.logsteps, df.gb_complex): ",cor( df.logsteps, df.gb_complex))
  println("cor( df.logsteps, df.degen): ",cor( df.logsteps, df.degen))
  println("cor( df.logsteps, df.sdegen): ",cor( df.logsteps, df.sdegen))
  println("cor( df.logsteps, df.f_mutinf): ",cor( df.logsteps, df.f_mutinf))
  println("cor( df.complex, df.f_mutinf): ",cor( df.complex, df.f_mutinf))
  println("cor( df.gb_complex, df.f_mutinf): ",cor( df.gb_complex, df.f_mutinf))
  println("cor( df.logsteps, df.mutrobust): ",cor( df.logsteps, df.mutrobust))
  println("cor( df.mutrobust, df.complex: ",cor( df.mutrobust, df.complex))
  println("cor( df.mutrobust, df.degen: ",cor( df.mutrobust, df.degen))
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  open( csvfile, "w" ) do f
    println(f,"# date and time: ",Dates.now())
    println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
    println(f,"# run time in minutes: ",ttime/60)
    println(f,"# funcs: ", Main.CGP.default_funcs(numinputs[end]))
    println(f,"# cor( df.logsteps, df.complex): ",cor( df.logsteps, df.complex))
    println(f,"# cor( df.logsteps, df.gb_complex): ",cor( df.logsteps, df.gb_complex))
    println(f,"# cor( df.logsteps, df.degen): ",cor( df.logsteps, df.degen))
    println(f,"# cor( df.logsteps, df.sdegen): ",cor( df.logsteps, df.sdegen))
    println(f,"# cor( df.logsteps, df.f_mutinf): ",cor( df.logsteps, df.f_mutinf))
    println(f,"# cor( df.complex, df.f_mutinf): ",cor( df.complex, df.f_mutinf))
    println(f,"# cor( df.gb_complex, df.f_mutinf): ",cor( df.gb_complex, df.f_mutinf))
    println(f,"# cor( df.logsteps, df.mutrobust): ",cor( df.logsteps, df.mutrobust))
    println(f,"# cor( df.mutrobust, df.complex: ",cor( df.mutrobust, df.complex))
    println(f,"# cor( df.mutrobust, df.degen: ",cor( df.mutrobust, df.degen))
    #println(f,"# nodearity: ",nodearity)
    #println(f,"# active_only: ",active_only)
    CSV.write( f, df, append=true, writeheader=true )
  end
  df
end

# Tentative conclusion:  10/6/20:  The robustness and evolvability computed are genotypic rather than phenotypic.
# Not what we really want.  The computations in geno_robustness() and geno_complexity() in Evolvability.jl should be correct.
function run_geno_pheno( iterations::Int64, numinputs::IntRange, numoutputs::IntRange,
    numinteriors::IntRange, goallistlength::IntRange, maxsteps::IntRange, levelsback::IntRange ;
    base::Float64=2.0, allgoals::Bool=false, active_only::Bool=false, gl_repetitions::IntRange=1)
  max_numinteriors = collect(numinteriors)[end]
  maxints_for_degen = Main.CGP.maxints_for_degen
  test_MyInt(max_numinteriors)
  nodearity = 2
  hamming_sel = true
  all_gl = allgoals ? [ [[x]] for x = MyInt(0):MyInt(2^2^numinputs-1) ] : [Vector{MyInt}[]]
  repeat_limit=10
  gp_result_list = geno_pheno_result_type[]
  df = DataFrame()
  df.goallist=Vector{Vector{MyInt}}[]
  df.numinputs=Int64[]
  df.numoutputs=Int64[]
  df.numints=Int64[]
  df.levsback=Int64[]
  df.ngoals=Int64[]
  df.maxsteps=Int64[]
  df.gl_reps=Int64[]
  df.steps=Int64[]
  df.logsteps=Float64[]
  df.avgfit=Float64[]
  df.nactive=Int64[]
  df.redund=Float64[]
  df.complex=Float64[]
  df.gb_complex=Float64[]
  df.degen=Float64[]
  df.sdegen=Float64[]
  df.f_mutinf=Float64[]
  df.mutrobust=Float64[]
  df.evolvability=Float64[]
  println("size(df): ",size(df))
  for num_inputs = numinputs
    for num_outputs = numoutputs
      #fit_limit = Float64(num_outputs)
      funcs = default_funcs(num_inputs)
      for num_goals = goallistlength
        for gl_reps = gl_repetitions
          #gl = Vector{MyInt}[]
          for num_interiors = numinteriors
            for levsback = levelsback
              println("numinputs: ",num_inputs,"  numoutputs: ",num_outputs,"  numints: ",num_interiors,
                "  levsback: ", levsback, "  numgoals: ", num_goals)
              for max_steps = maxsteps
                for gl = (allgoals ? all_gl : [Vector{MyInt}[]])
                  for _ = 1:iterations
                    p = Parameters( num_inputs, num_outputs, num_interiors, levsback )
                    rr = gp_result( gl, p, num_goals, hamming_sel, active_only, max_steps, gl_reps )
                    push!(gp_result_list,rr)
                  end
                end
              end
            end
          end
        end
      end
    end
  end
  new_gp_result_list = pmap(r->run_gp_evolve!(r,maxints_for_degen=maxints_for_degen,gl_repetitions=gl_repetitions,base=base,repeat_limit=repeat_limit),gp_result_list)
  #new_gp_result_list = map(r->run_gp_evolve!(r,maxints_for_degen=Main.CGP.maxints_for_degen,gl_repetitions=gl_repetitions,base=base,repeat_limit=repeat_limit),gp_result_list)
  for r = new_gp_result_list
    new_row = gp_result_to_tuple(r)
    #println("length new row: ",length(new_row))
    Base.push!( df, new_row )
  end
  df
end

# Tentative conclusion:  10/6/20:  The robustness and evolvability computed are genotypic rather than phenotypic.
function run_gp_evolve!( rr::geno_pheno_result_type; maxints_for_degen::Int64, gl_repetitions::Int64=1, base::Float64=2.0,
      repeat_limit::Int64=10 )
  #repeat_limit = 10
  #println("run_gp_evolve! fault_tol: ",rr.fault_tol)
  nodearity = 2   # built-in default
  p = Parameters( numinputs=rr.numinputs, numoutputs=rr.numoutputs, numinteriors=rr.numints, numlevelsback=rr.levelsback )
  #print_parameters( p )
  #=
  if length(rr.goallist) > 0
    println("rr.goallist: ",rr.goallist)
  end
  =#
  if length(rr.goallist) == 0
    rr.goallist = randgoallist(rr.ngoals,rr.numinputs,rr.numoutputs,repetitions=gl_repetitions)
  end
  funcs = default_funcs(rr.numinputs)
  c = random_chromosome( p, funcs )
  repeat = 0
  (c,rr.steps,output,matched_goals,matched_goals_list) = mut_evolve(c,rr.goallist,funcs,rr.maxsteps)
  if rr.steps == rr.maxsteps
    println("goal: ",rr.goallist[1])
  end
  while repeat < repeat_limit && rr.steps == rr.maxsteps
    rr.numints += 1
    rr.levelsback += 1
    p_repeat = Parameters( numinputs=rr.numinputs, numoutputs=rr.numoutputs, numinteriors=rr.numints, numlevelsback=rr.levelsback )
    println("repeating function mut_evolve with numints: ",rr.numints,"  and with levsback: ",rr.levelsback)
    c = random_chromosome( p_repeat, funcs )
    (c,rr.steps,output,matched_goals,matched_goals_list) =
        mut_evolve(c,rr.goallist,funcs,rr.maxsteps)
    #print_build_chromosome( c )
    repeat += 1
  end
  if repeat == repeat_limit
    println("run_gp_evolve!  repeat == repeat_limit ")
  end
  #if rr.numints > p.numinteriors
  #  println("rr.numints - p.numinteriors: ",rr.numints - p.numinteriors)
  #end
  rr.logsteps = log( 1.0 + rr.steps + rr.maxsteps*(rr.numints - p.numinteriors) )
  rr.avgfit = c.fitness
  rr.nactive = number_active( c )
  rr.redund = redundancy( c, base=base )
  rr.complex = rr.numints <= maxints_for_degen ? complexity5( c, base=base ) : 0.0
  #rr.gb_complex = rr.numints <= maxints_for_degen ? gb_complexity_chrome( c, base=base ) : 0.0 # meaningless 9/27
  rr.degen = rr.numints <= maxints_for_degen ? degeneracy( c, base=base ) : 0.0
  rr.sdegen = rr.numints <= maxints_for_degen ? degeneracy( c, base=base, mutinf=mutinf2 ) : 0.0
  rr.f_mutinf = fmi_chrome( c )
  mutall = mutate_all( c, funcs, robustness_only=true )
  rr.mutrobust = mutall[1]
  rr.evolvability = mutall[2]
  rr
end

function gp_result( gl::Vector{Vector{MyInt}}, p::Parameters, num_goals::Int64, hamming_sel::Bool, active_only::Bool, max_steps::Int64,
      gl_reps::Int64 )
  geno_pheno_result_type(
    #Vector{MyInt}[],
    gl,
    p.numinputs,
    p.numoutputs,
    p.numinteriors,
    p.numlevelsback,
    num_goals,
    hamming_sel,
    active_only,
    max_steps,
    gl_reps,
    0,      # steps
    0.0,    # logsteps
    0.0,    # avgfit
    0,      # nactive
    0.0,    # redund
    0.0,    # complex
    0.0,    # gb_complex
    0.0,    # degen
    0.0,    # sdegen
    0.0,    # f_mutinf
    0.0,    # mutrobust
    0.0     # evolvability
  )
end

function gp_result_to_tuple( rr::geno_pheno_result_type )
 (
  rr.goallist,
  rr.numinputs,
  rr.numoutputs,
  rr.numints,
  rr.levelsback,
  rr.ngoals,
  rr.maxsteps,
  rr.gl_reps,
  rr.steps,
  rr.logsteps,
  rr.avgfit,
  rr.nactive,
  rr.redund,
  rr.complex,
  rr.gb_complex,
  rr.degen,
  rr.sdegen,
  rr.f_mutinf,
  rr.mutrobust,
  rr.evolvability
 )
end
#=
# Analyze a dataframe that was created by run_geno_pheno_evolution()
function analyze_dataframe( df_filename::String )
  df = CSV.read(df_filename,header=true,comment="#")
  println("cor(logsteps,complex): ",cor(df.logsteps,df.complex))
  df
end
=#


text:
const c=299792458.
const mu0=4*pi*1e-7
const mur=1
const eps0=1./(mu0*c^2)
const epsr=1

using PyPlot

function urc(Q,f,f0)
    w=f/f0
    I=1/sqrt(1+Q^2*(w-1/w)^2)
    return I
end

const a=2.9#3.11
const b=3.7#3.62
const d=8.7#4.84

const df=2e6
const nf=1
const f=[274e6]#[224e6]#[208e6]#[382e6]#linspace(f0,f1,nf)#array([20e6,40e6,45e6,56e6,65e6])#linspace(20e6,.3e9,281)#linspace(20e6,3e8,1121)
const V=a*b*d
const Q4=16*pi^2*V*f.^3/(1*c^3)
const Q1=3*V./(2*1*(2./sqrt(2*pi*f.*mur*mu0*1*1e5)))/500
const Q=1./(1./Q1+1./Q4)

const x=1
dy=0.04
dz=0.04
ny=int(b/dy)+1
nz=int(d/dz)+1
const y=linspace(0,b,ny)
const z=linspace(0,d,nz)

n_steps=180
phi=linspace(2*pi/n_steps,2*pi,n_steps)
x_b=ones(n_steps)
y_b=cos(phi)+1.5
z_b=sin(phi)+1.5

#z_b=[2,2.8,4.5]
#x_b=[1,1,1]
#y_b=[0.5,3,1]

n_steps=length(x_b)

E0=1
H0=1/(c*mu0*mur)

Ex=zeros(Complex128,nf,n_steps,ny,nz)
#Ey=zeros(Complex128,nf,n_steps,Npt)
#Ez=zeros(Complex128,nf,n_steps,Npt)
for u=1:nf
  n_x=int(f[u]*(1+100/Q[u])*2*a/c)
  n_y=int(f[u]*(1+100/Q[u])*2*b/c)
  n_z=int(f[u]*(1+100/Q[u])*2*d/c)
  M=[]
  N=[]
  P=[]
  F=[]
  AMP=[]
  for mm=0:n_x
    for nn=0:n_y
      for pp=0:n_z
        fmnp=c/2*sqrt((mm/a)^2+(nn/b)^2+(pp/d)^2)
        #if abs(fmnp-f[u])<100*f[u]/Q[u]
        Amp=urc(Q[u],f[u],fmnp)
        if Amp>0.01
          M=vcat(M,mm)
          N=vcat(N,nn)
          P=vcat(P,pp)
          F=vcat(F,fmnp)
          AMP=vcat(AMP,Amp)
        end
      end
    end
  end
  nm=length(M)
  println("$nm resonant freq.")
  f_coeffTMx_x=zeros(n_steps,nm)
  f_coeffTMx_y=zeros(n_steps,nm)
  f_coeffTMx_z=zeros(n_steps,nm)
  f_coeffTEx_x=zeros(Complex128,n_steps,nm)
  f_coeffTEx_y=zeros(Complex128,n_steps,nm)
  f_coeffTEx_z=zeros(Complex128,n_steps,nm)
  for pos=1:n_steps
    ETMx_x=zeros(ny,nz)
    ETMx_y=zeros(ny,nz)
    ETMx_z=zeros(ny,nz)
    ETEx_y=zeros(Complex128,ny,nz)
    ETEx_z=zeros(Complex128,ny,nz)
    xe=x_b[pos]
    ye=y_b[pos]
    ze=z_b[pos]
    println("Step $pos/$n_steps")
    for k=1:nm
      fmnp=F[k]
      Amp=AMP[k]
      kx=M[k]*pi/a
      ky=N[k]*pi/b
      kz=P[k]*pi/d
      k2mnp=kx^2+ky^2+kz^2
      ssc_ze=sum(sin(kx*xe).*sin(ky*ye).*cos(kz*ze))
      scs_ze=sum(sin(kx*xe).*cos(ky*ye).*sin(kz*ze))
      css_ze=sum(cos(kx*xe).*sin(ky*ye).*sin(kz*ze))
      ssc_xe=sum(sin(ky*ye).*sin(kz*ze).*cos(kx*xe))
      scs_xe=sum(sin(ky*ye).*cos(kz*ze).*sin(kx*xe))
      css_xe=sum(cos(ky*ye).*sin(kz*ze).*sin(kx*xe))
      ssc_ye=sum(sin(kz*ze).*sin(kx*xe).*cos(ky*ye))
      scs_ye=sum(sin(kz*ze).*cos(kx*xe).*sin(ky*ye))
      css_ye=sum(cos(kz*ze).*sin(kx*xe).*sin(ky*ye))

      if M[k]>0 && N[k]>0 && P[k]>0
        #TM modes z
        f_coeffTMx_z[pos,k] = -kx*kz/(k2mnp-kz^2)*css_ze
        #f_coeffTMy_z[pos,k] = ky*kz/(k2mnp-kz^2)*scs_ze
        #f_coeffTMz_z[pos,k] = ssc_ze
        #TE modes z
        f_coeffTEx_z[pos,k] = -1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kz^2)*css_ze
        #f_coeffTEy_z[pos,k] = 1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-kz^2)*scs_ze

        #TM modes x
        #f_coeffTMy_x[pos,k] = -ky*kx/(k2mnp-kx^2)*css_xe
        #f_coeffTMz_x[pos,k] = kz*kx/(k2mnp-kx^2)*scs_xe
        f_coeffTMx_x[pos,k] = ssc_xe
        #TE modes x
        #f_coeffTEy_x[pos,k] = -1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-kx^2)*css_xe
        #f_coeffTEz_x[pos,k] = 1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kx^2)*scs_xe

        ##TM modes y
        #f_coeffTMz_y[pos,k] += -kz*ky/(k2mnp-ky^2)*css_ye
        f_coeffTMx_y[pos,k] = kx*ky/(k2mnp-ky^2)*scs_ye
        #f_coeffTMy_y[pos,k] = ssc_ye
        ##TE modes y
        #f_coeffTEz_y[pos,k] = -1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-ky^2)*css_ye
        f_coeffTEx_y[pos,k] = 1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-ky^2)*scs_ye
      end

      #if M[k]>0 && N[k]>0 && P[k]==0
        #f_coeffTMz_z[pos,k] = ssc_ze #TM mode z
        #f_coeffTEz_y[pos,k] = -1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-ky^2)*css_ye #TE mode y
        #f_coeffTEz_x[pos,k] = 1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kx^2)*scs_xe #TE mode x
      #end

      if M[k]==0 && N[k]>0 && P[k]>0
        f_coeffTEx_z[pos,k] += -1im*2*pi*fmnp*mur*mu0*ky/(k2mnp-kz^2)*css_ze #TE mode z
        f_coeffTMx_x[pos,k] += ssc_xe #TM mode x
        f_coeffTEx_y[pos,k] += 1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-ky^2)*scs_ye #TE mode y
      end

      #if M[k]>0 && N[k]==0 && P[k]>0
        #f_coeffTEy_z[pos,k] = 1im*2*pi*fmnp*mur*mu0*kx/(k2mnp-kz^2)*scs_ze #TE mode z
        #f_coeffTMy_y[pos,k] = ssc_ye #TM mode y
        #f_coeffTEy_x[pos,k] = -1im*2*pi*fmnp*mur*mu0*kz/(k2mnp-kx^2)*css_xe #TE mode x
      #end

      for i=1:ny
        for j=1:nz
          ssc_z=sin(kx*x)*sin(ky*y[i])*cos(kz*z[j])
          scs_z=sin(kx*x)*cos(ky*y[i])*sin(kz*z[j])
          css_z=cos(kx*x)*sin(ky*y[i])*sin(kz*z[j])
          ssc_x=sin(ky*y[i])*sin(kz*z[j])*cos(kx*x)
          scs_x=sin(ky*y[i])*cos(kz*z[j])*sin(kx*x)
          css_x=cos(ky*y[i])*sin(kz*z[j])*sin(kx*x)
          ssc_y=sin(kz*z[j])*sin(kx*x)*cos(ky*y[i])
          scs_y=sin(kz*z[j])*cos(kx*x)*sin(ky*y[i])
          css_y=cos(kz*z[j])*sin(kx*x)*sin(ky*y[i])

          if M[k]>0 && N[k]>0 && P[k]>0#TM modes z
            #TM modes z
            ETMx_z[i,j] += E0*Amp*f_coeffTMx_z[pos,k]*css_z
            #ETMy_z[i,j] += E0*Amp*f_coeffTMy_z[pos,k]*scs_z
            #ETMz_z[i,j] += E0*Amp*f_coeffTMz_z[pos,k]*ssc_z

            #TE modes z
            ETEx_z[i,j] += H0*Amp*f_coeffTEx_z[pos,k]*css_z
            #ETEy_z[i,j]+=H0*Amp*f_coeffTEy_z[pos,k]*scs_z

            #TM modes x
            #ETMy_x[i,j]+=E0*Amp*f_coeffTMy_x[pos,k]*css_x
            #ETMz_x[i,j]+=E0*Amp*f_coeffTMz_x[pos,k]*scs_x
            ETMx_x[i,j]+=E0*Amp*f_coeffTMx_x[pos,k]*ssc_x

            #TE modes x
            #ETEy_x[i,j]+=H0*Amp*f_coeffTEy_x[pos,k]*css_x
            #ETEz_x[i,j]+=H0*Amp*f_coeffTEz_x[pos,k]*scs_x

            ##TM modes y
            #ETMz_y[i,j]+=E0*Amp*f_coeffTMz_y[pos,k]*css_y
            ETMx_y[i,j]+=E0*Amp*f_coeffTMx_y[pos,k]*scs_y
            #ETMy_y[i,j]+=E0*Amp*f_coeffTMy_y[pos,k]*ssc_y

            ##TE modes y
            #ETEz_y[i,j]+=H0*Amp*f_coeffTEz_y[pos,k]*css_y
            ETEx_y[i,j]+=H0*Amp*f_coeffTEx_y[pos,k]*scs_y
          end

          #if M[k]>0 && N[k]>0 && P[k]==0
            #ETMz_z[i,j] += E0*Amp*f_coeffTMz_z[pos,k]*ssc_z #TM mode z
            #ETEz_y[i,j]+=H0*Amp*f_coeffTEz_y[pos,k]*css_y #TE mode y
            #ETEz_x[i,j]+=H0*Amp*f_coeffTEz_x[pos,k]*scs_x #TE mode x
          #end

          if M[k]==0 && N[k]>0 && P[k]>0
            ETEx_z[i,j] += H0*Amp*f_coeffTEx_z[pos,k]*css_z #TE mode z
            ETMx_x[i,j]+=E0*Amp*f_coeffTMx_x[pos,k]*ssc_x #TM mode x
            ETEx_y[i,j]+=H0*Amp*f_coeffTEx_y[pos,k]*scs_y #TE mode y
          end

          #if M[k]>0 && N[k]==0 && P[k]>0
            #ETEy_z[i,j]+=H0*Amp*f_coeffTEy_z[pos,k]*scs_z #TE mode z
            #ETMy_y[i,j]+=E0*Amp*f_coeffTMy_y[pos,k]*ssc_y #TM mode y
            #ETEy_x[i,j]+=H0*Amp*f_coeffTEy_x[pos,k]*css_x #TE mode x
          #end
        end
      end
    end
    Ex[u,pos,:]=ETMx_x+ETMx_y+ETMx_z+ETEx_y+ETEx_z
    #Ey[u,pos,:]=ETMy_x+ETMy_y+ETMy_z+ETEy_x+ETEy_z
    #Ez[u,pos,:]=ETMz_x+ETMz_y+ETMz_z+ETEz_y+ETEz_x

  end
  M=[]
  N=[]
  P=[]
end
u=1
cmax=maximum(abs(Ex[u,:,:,:]))

using PyCall
@pyimport matplotlib.pyplot as plt

for o=1:n_steps
  plt.figure(num=2,figsize=(8,4))
  ftext=round(f[u]/1e6*100)/100
  phitext=round(phi[o]*180/pi)
  plt.title("\$\\vert E_x\\vert\$, in plane \$x=1\$ m, \$ \\phi= $phitext^{\\circ}\$",fontsize=15)
  plt.pcolor(z,y,squeeze(squeeze(abs(Ex[u,o,:,:]),1),1),cmap="RdBu_r")
  plt.plot(z_b[o],y_b[o],marker="o",markerfacecolor="w", markersize=10)
  plt.clim(0,0.15)#cmax)
  plt.colorbar(format="%.2f",orientation="horizontal")
  plt.xlabel("\$ z\$/m",fontsize=15)
  plt.ylabel("\$ y\$/m",fontsize=15)
  plt.axis("image")
  plt.grid()
  plt.savefig("Ex_$o.png",format="png",bbox_inches="tight",dpi=96)# , transparent="True")
  plt.clf()
end


text:
module Jasmin

using LinearAlgebra, Combinatorics, Test

import Base.println

#=
export AbstractStandard, BaseEnumeration, StandardSimplexe, solve!,  Simplexe,
     AbstractStatus, Optimal, Infeasible, Unknown, Unbounded,
     xstar, vstar
=#
include("utils.jl")
include("Status.jl")

abstract type AbstractSolver{T} end
abstract type AbstractLP{T} end

include("LinearProblemRepresentation/main.jl")
include("Solvers/main.jl")

include("test/main.jl")
end # module

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
"""
spin13.jl - vanilla
"""

WDIR = get(ENV, "ROBUST_QOC_PATH", "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

using Altro
using HDF5
using LinearAlgebra
using RobotDynamics
using StaticArrays
using TrajectoryOptimization
const RD = RobotDynamics
const TO = TrajectoryOptimization

# paths
const EXPERIMENT_META = "spin"
const EXPERIMENT_NAME = "spin13"
const SAVE_PATH = abspath(joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME))

# problem
const CONTROL_COUNT = 1
const STATE_COUNT = 2
const ASTATE_SIZE = STATE_COUNT * HDIM_ISO + 3 * CONTROL_COUNT
const ACONTROL_SIZE = CONTROL_COUNT
# state indices
const STATE1_IDX = 1:HDIM_ISO
const STATE2_IDX = STATE1_IDX[end] + 1:STATE1_IDX[end] + HDIM_ISO
const INTCONTROLS_IDX = STATE2_IDX[end] + 1:STATE2_IDX[end] + CONTROL_COUNT
const CONTROLS_IDX = INTCONTROLS_IDX[end] + 1:INTCONTROLS_IDX[end] + CONTROL_COUNT
const DCONTROLS_IDX = CONTROLS_IDX[end] + 1:CONTROLS_IDX[end] + CONTROL_COUNT
# control indices
const D2CONTROLS_IDX = 1:CONTROL_COUNT

# model
struct Model <: AbstractModel
end
@inline RD.state_dim(::Model) = ASTATE_SIZE
@inline RD.control_dim(::Model) = ACONTROL_SIZE

# dynamics
function RD.discrete_dynamics(::Type{RK3}, model::Model, astate::SVector,
                              acontrol::SVector, time::Real, dt::Real)
    h_prop = exp(dt * (FQ_NEGI_H0_ISO + astate[CONTROLS_IDX[1]] * NEGI_H1_ISO))
    state1 = h_prop * astate[STATE1_IDX]
    state2 = h_prop * astate[STATE2_IDX]
    intcontrols = astate[INTCONTROLS_IDX[1]] + astate[CONTROLS_IDX[1]] * dt
    controls = astate[CONTROLS_IDX[1]] + astate[DCONTROLS_IDX[1]] * dt
    dcontrols = astate[DCONTROLS_IDX[1]] + acontrol[D2CONTROLS_IDX[1]] * dt

    astate_ = [
        state1; state2; intcontrols; controls; dcontrols;
    ]

    return astate_
end

# main
function run_traj(;gate_type=zpiby2, evolution_time=30., solver_type=altro,
                  sqrtbp=false, integrator_type=rk3, smoke_test=false,
                  dt_inv=Int64(1e1), constraint_tol=1e-8, al_tol=1e-4,
                  pn_steps=2, max_penalty=1e11, verbose=true, save=true, max_iterations=Int64(2e5),
                  max_cost_value=1e8, qs=[1e0, 1e0, 1e0, 1e-1, 1e-1],
                  benchmark=false,)
    # model configuration
    model = Model()
    n = state_dim(model)
    m = control_dim(model)
    t0 = 0.

    # initial state
    x0 = zeros(n)
    x0[STATE1_IDX] = IS1_ISO_
    x0[STATE2_IDX] = IS2_ISO_
    x0 = SVector{n}(x0)

    # final state
    if gate_type == xpiby2
        target_state1 = Array(XPIBY2_ISO_1)
        target_state2 = Array(XPIBY2_ISO_2)
    elseif gate_type == ypiby2
        target_state1 = Array(YPIBY2_ISO_1)
        target_state2 = Array(YPIBY2_ISO_2)
    elseif gate_type == zpiby2
        target_state1 = Array(ZPIBY2_ISO_1)
        target_state2 = Array(ZPIBY2_ISO_2)
    end
    xf = zeros(n)
    xf[STATE1_IDX] = target_state1
    xf[STATE2_IDX] = target_state2
    xf = SVector{n}(xf)

    # control amplitude constraint
    x_max = fill(Inf, n)
    x_max[CONTROLS_IDX] .= MAX_CONTROL_NORM_0
    x_max = SVector{n}(x_max)
    x_min = fill(-Inf, n)
    x_min[CONTROLS_IDX] .= -MAX_CONTROL_NORM_0
    x_min = SVector{n}(x_min)

    # control amplitude constraint at boundary
    x_max_boundary = fill(Inf, n)
    x_max_boundary[CONTROLS_IDX] .= 0
    x_max_boundary = SVector{n}(x_max_boundary)
    x_min_boundary = fill(-Inf, n)
    x_min_boundary[CONTROLS_IDX] .= 0
    x_min_boundary = SVector{n}(x_min_boundary)

    # initial trajectory
    dt = dt_inv^(-1)
    N = Int(floor(evolution_time * dt_inv)) + 1
    U0 = [SVector{m}(
        fill(1e-4, CONTROL_COUNT)
    ) for k = 1:N-1]
    X0 = [SVector{n}([
        fill(NaN, n);
    ]) for k = 1:N]
    Z = Traj(X0, U0, dt * ones(N))

    # cost function
    Q = Diagonal(SVector{n}([
        fill(qs[1], STATE_COUNT * HDIM_ISO); # ψ1, ψ2
        fill(qs[2], CONTROL_COUNT); # ∫a
        fill(qs[3], CONTROL_COUNT); # a
        fill(qs[4], CONTROL_COUNT); # ∂a
    ]))
    Qf = Q * N
    R = Diagonal(SVector{m}([
        fill(qs[5], CONTROL_COUNT); # ∂2a
    ]))
    objective = LQRObjective(Q, R, Qf, xf, N)

    # constraints
    # must satisfy control amplitude bound
    control_bnd = BoundConstraint(n, m, x_max=x_max, x_min=x_min)
    # must statisfy conrols start and end at 0
    control_bnd_boundary = BoundConstraint(n, m, x_max=x_max_boundary, x_min=x_min_boundary)
    # must reach target state, must have zero net flux
    target_astate_constraint = GoalConstraint(xf, [STATE1_IDX; STATE2_IDX; INTCONTROLS_IDX])
    # must obey unit norm
    norm_constraints = [NormConstraint(n, m, 1, TO.Equality(), idx)
                        for idx in [STATE1_IDX, STATE2_IDX]]

    constraints = ConstraintList(n, m, N)
    add_constraint!(constraints, control_bnd, 2:N-2)
    add_constraint!(constraints, control_bnd_boundary, N-1:N-1)
    add_constraint!(constraints, target_astate_constraint, N:N);
    for norm_constraint in norm_constraints
        add_constraint!(constraints, norm_constraint, 2:N-1)
    end

    # solve problem
    prob = Problem{IT_RDI[integrator_type]}(model, objective, constraints,
                                            x0, xf, Z, N, t0, evolution_time)
    solver = ALTROSolver(prob)
    verbose_pn = verbose ? true : false
    verbose_ = verbose ? 2 : 0
    projected_newton = solver_type == altro ? true : false
    constraint_tolerance = solver_type == altro ? constraint_tol : al_tol
    iterations_inner = smoke_test ? 1 : 300
    iterations_outer = smoke_test ? 1 : 30
    n_steps = smoke_test ? 1 : pn_steps
    set_options!(solver, square_root=sqrtbp, constraint_tolerance=constraint_tolerance,
                 projected_newton_tolerance=al_tol, n_steps=n_steps,
                 penalty_max=max_penalty, verbose_pn=verbose_pn, verbose=verbose_,
                 projected_newton=projected_newton, iterations_inner=iterations_inner,
                 iterations_outer=iterations_outer, iterations=max_iterations,
                 max_cost_value=max_cost_value)
    if benchmark
        benchmark_result = Altro.benchmark_solve!(solver)
    else
        benchmark_result = 0
        Altro.solve!(solver)
    end

    # post-process
    acontrols_raw = TO.controls(solver)
    acontrols_arr = permutedims(reduce(hcat, map(Array, acontrols_raw)), [2, 1])
    astates_raw = TO.states(solver)
    astates_arr = permutedims(reduce(hcat, map(Array, astates_raw)), [2, 1])
    Q_raw = Array(Q)
    Q_arr = [Q_raw[i, i] for i in 1:size(Q_raw)[1]]
    Qf_raw = Array(Qf)
    Qf_arr = [Qf_raw[i, i] for i in 1:size(Qf_raw)[1]]
    R_raw = Array(R)
    R_arr = [R_raw[i, i] for i in 1:size(R_raw)[1]]
    cidx_arr = Array(CONTROLS_IDX)
    d2cidx_arr = Array(D2CONTROLS_IDX)
    cmax = TrajectoryOptimization.max_violation(solver)
    cmax_info = TrajectoryOptimization.findmax_violation(TO.get_constraints(solver))
    iterations_ = Altro.iterations(solver)

    result = Dict(
        "acontrols" => acontrols_arr,
        "controls_idx" => cidx_arr,
        "d2controls_dt2_idx" => d2cidx_arr,
        "evolution_time" => evolution_time,
        "astates" => astates_arr,
        "Q" => Q_arr,
        "Qf" => Qf_arr,
        "R" => R_arr,
        "cmax" => cmax,
        "cmax_info" => cmax_info,
        "dt" => dt,
        "solver_type" => Integer(solver_type),
        "sqrtbp" => Integer(sqrtbp),
        "max_penalty" => max_penalty,
        "constraint_tol" => constraint_tol,
        "al_tol" => al_tol,
        "gate_type" => Integer(gate_type),
        "save_type" => Integer(jl),
        "integrator_type" => Integer(integrator_type),
        "iterations" => iterations_,
        "max_iterations" => max_iterations,
        "pn_steps" => pn_steps,
        "max_cost_value" => max_cost_value,
        "benchmark_result" => benchmark_result,
    )

    # save
    if save
        save_file_path = generate_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        println("Saving this optimization to $(save_file_path)")
        h5open(save_file_path, "cw") do save_file
            for key in keys(result)
                write(save_file, key, result[key])
            end
        end
        result["save_file_path"] = save_file_path
    end

    result = benchmark ? benchmark_result : result

    return result
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
using StatsFuns
using StatsBase
using Statistics
using Distributed
using Printf
using Crayons
using Base.Iterators
using Dates

# Probably will not be used in favor of the weighted mean/cov functions provided by
# StatsBase
function pmcmc_mean_cov(param_matrix::Array{Float64,2}, log_w::Vector{Float64})
    mean_vec = fill(0.0, size(param_matrix)[1])
    cov_mat = fill(0.0, size(param_matrix)[1], size(param_matrix)[1])
    total_weight = 0.0
    w_max = maximum(log_w)
    w_n = Vector{Float64}(undef, size(param_matrix)[2])
    for n = 1:size(param_matrix)[2]
        w_n[n] = exp(log_w[n] - w_max)
        total_weight += w_n[n]
        mean_vec .+= w_n[n] * (@view param_matrix[:, n])
    end
    mean_vec ./= total_weight
    p_vec = Vector{Float64}(undef, size(param_matrix)[1])
    for n = 1:size(param_matrix)[2]
        p_vec = param_matrix[:, n] .- mean_vec
        cov_mat += w_n[n] * (p_vec * p_vec')
    end
    cov_mat ./= total_weight
    return (mean_vec, cov_mat)
end

function pmcmc_mean_cov(params::Array{Array{Q,1},1}, log_w::Vector{Q}) where {Q<:Any}
    pmcmc_mean_cov(hcat(params...), log_w)
end

function pmcmc_mean_cov(params::Vector{T}, log_w::Vector{Q}) where {T<:NamedTuple,Q<:Any}
    pmcmc_mean_cov(hcat(map(x -> vcat(x...), params)...), log_w)
end

# Given a parameter (for dispatch), sample from a truncated MVT normal distribution
function trunc_normal_jump(
    x::Vector,
    mu::Vector,
    sigma::Matrix;
    lower_bound = (x, n) -> -Inf,
    upper_bound = (x, n) -> Inf,
    perturb = 1e-4,
)
    sigma = sigma + I * perturb
    A = Array{Float64,2}(undef, length(x) - 1, length(x))
    sigma_n = Array{Float64,1}(undef, length(x))
    n_x = length(x)
    x_current = copy(x)
    for n = 1:length(x)
        A = sigma[1:n_x.!=n, 1:n_x.!=n] \ sigma[1:n_x.!=n, n]
        sigma_n = sqrt(sigma[n, n] - (sigma[n, 1:n_x.!=n]' * A))
        mean_n = mu[n] .+ (A'*(x_current[1:n_x.!=n].-mu[1:n_x.!=n]))[1, 1]
        l = lower_bound(x_current, n)
        u = upper_bound(x_current, n)
        x_current[n] = rand(Truncated(Normal(mean_n, sigma_n), l, u))
    end
    return x_current
end

function trunc_normal_dens(
    x::Vector,
    x2::Vector,
    mu::Vector,
    sigma::Matrix;
    lower_bound = (x, n) -> -Inf,
    upper_bound = (x, n) -> Inf,
    perturb = 1e-4,
)
    sigma = sigma + I * perturb
    A = Array{Float64,2}(undef, length(x) - 1, length(x))
    sigma_n = Array{Float64,1}(undef, length(x))
    n_x = length(x)
    x_current = copy(x)
    ldens = 0.0
    for n = 1:length(x)
        A[:, n] = sigma[1:n_x.!=n, 1:n_x.!=n] \ sigma[1:n_x.!=n, n]
        sigma_n[n] = sqrt(sigma[n, n] - (sigma[n, 1:n_x.!=n]' * A[:, n]))
        mean_n = mu[n] .+ (A[:, n]'*(x_current[1:n_x.!=n].-mu[1:n_x.!=n]))[1, 1]
        l = lower_bound(x_current, n)
        u = upper_bound(x_current, n)
        ldens += logpdf(Truncated(Normal(mean_n, sigma_n[n]), l, u), x2[n])
        x_current[n] = x2[n]
    end
    return ldens
end

function smc_pmcmc_proposal(
    theta_current::Vector,
    pmcmc_mean::Vector,
    pmcmc_cov::Array;
    C = (2.38^2) / length(theta_current),
    kwargs...,
)
    return trunc_normal_jump(theta_current, theta_current, pmcmc_cov * C; kwargs...)
end

function smc_pmcmc_proposal_logdens(
    theta_current::Vector,
    theta_proposed::Vector,
    pmcmc_mean::Vector,
    pmcmc_cov::Matrix;
    C = (2.38^2) / length(theta_current),
    kwargs...,
)
    return trunc_normal_dens(
        theta_current,
        theta_proposed,
        theta_current,
        pmcmc_cov * C;
        kwargs...,
    )
end

# Utilities

function calculate_ess(xs)
    x_max = maximum(xs)
    x_sum = 0.0
    x_squared_sum = 0.0
    for n = 1:length(xs)
        exp_x_n = exp(xs[n] - x_max)
        x_sum += exp_x_n
        x_squared_sum += exp_x_n^2
    end

    return x_sum^2 / x_squared_sum
end

# This is a batch function which will run a filter from 1:T for a given vector of parameters in parallel

function dt_smc2_estimation(
    thetas,
    loglik_fun,
    prior_fun;
    grid_steps = 1e-8,
    ess_threshold = 0.5,
    pmcmc_theta_steps = 10,
    parallelize = true,
    kwargs...,
)
    N_θ = length(thetas)
    if parallelize
        ## Currently this pmap is type unstable, but this should
        ## not be too detrimental to performance
        logliks = pmap(loglik_fun, thetas)
    else
        logliks = map(loglik_fun, thetas)
    end
    ξ = 0.0

    acceptances = Array{Bool,3}(undef, N_θ, pmcmc_theta_steps, 0)

    while ξ < 1.0
        print(crayon"blue", "Picking new value for ξ\n")
        ξ_diff = grid_steps
        ess = calculate_ess(logliks * ξ_diff)

        ξ_diff, ess = xi_grid_search(logliks, ess_threshold)

        ξ = min(ξ_diff + ξ, 1.0)

        print(crayon"blue", @sprintf("Xi set to %.3e\nESS = %.3f\n", ξ, ess))

        thetas, logliks, acceptances_i = density_tempered_pmcmc(
            thetas,
            logliks,
            loglik_fun,
            prior_fun,
            ξ,
            ξ_diff;
            pmcmc_theta_steps = pmcmc_theta_steps,
            kwargs...,
        )
        acceptances = cat3(acceptances, acceptances_i)
    end

    return thetas, logliks, acceptances
end

function density_tempered_pmcmc(
    thetas,
    logliks,
    loglik_fun,
    prior_fun,
    ξ,
    ξ_diff;
    pmcmc_theta_steps = 10,
    rprop_pmh = smc_pmcmc_proposal,
    dprop_pmh = smc_pmcmc_proposal_logdens,
    kwargs...,
)
    N_θ = length(thetas)
    log_ws = logliks * ξ_diff

    ws = exp.(log_ws .- maximum(log_ws))
    pmcmc_mean, pmcmc_cov = pmcmc_mean_cov(thetas, log_ws)
    println("Param" * repeat(" ", 11) * "\tMean\t\tStd\t\t")
    for i = 1:length(pmcmc_mean)
        name = first(string(keys(thetas[1])[i]), 16)
        name *= repeat(" ", 17 - length(name))
        p_str = @sprintf("%s\t%.3e\t%.3e", name, pmcmc_mean[i], sqrt(pmcmc_cov[i, i]))
        println(p_str)
    end

    # A                 = wsample(1:N_θ, ws, N_θ)
    A = residual_resample(Random.GLOBAL_RNG, ws, N_θ)
    thetas = thetas[A]
    logliks = logliks[A]

    acceptances = Array{Bool,2}(undef, N_θ, pmcmc_theta_steps)
    for step_idx = 1:pmcmc_theta_steps
        starttime = Dates.now()
        res = pmap(
            pmcmc_propose_accept_reject,
            thetas,
            logliks,
            cycle([loglik_fun]),
            cycle([prior_fun]),
            cycle([ξ]),
            cycle([pmcmc_mean]),
            cycle([pmcmc_cov]),
            cycle([rprop_pmh]),
            cycle([dprop_pmh]),
        )
        for (i, r) in enumerate(res)
            thetas[i] = r[1]
            logliks[i] = r[2]
            acceptances[i, step_idx] = r[3]
        end
        timetaken = Dates.now() - starttime
        print(
            crayon"green",
            @sprintf(
                "Step %02d: Accepted %0.2f; Filter time: %f secs \n",
                step_idx,
                mean(acceptances[:, step_idx]),
                timetaken.value / 1000
            )
        )
    end
    print(
        crayon"green",
        @sprintf("Total θ moved at least once: %0.2f \n", mean(any(acceptances, dims = 2)))
    )
    return thetas, logliks, acceptances
end

function pmcmc_propose_accept_reject(
    theta,
    loglik,
    loglik_fun,
    prior_fun,
    ξ,
    pmcmc_mean,
    pmcmc_cov,
    rprop_pmh,
    dprop_pmh,
)
    # Propose new parameter
    theta_prop = rprop_pmh(theta, pmcmc_mean, pmcmc_cov)
    d_prop = dprop_pmh(theta, theta_prop, pmcmc_mean, pmcmc_cov)
    d_current = dprop_pmh(theta_prop, theta, pmcmc_mean, pmcmc_cov)

    # Evaluate likelihood of parameter
    loglik_prop = loglik_fun(theta_prop)

    # Evaluate prior of parameter
    prior = prior_fun(theta)
    prior_prop = prior_fun(theta_prop)

    # Calculate log acceptance probability
    log_prior_ratio = prior_prop - prior
    log_lik_ratio = ξ * (loglik_prop - loglik)
    log_prop_ratio = d_current - d_prop

    log_posterior_proposed = loglik_prop + prior_prop
    log_acceptance = log_prior_ratio + log_lik_ratio + log_prop_ratio
    accepted = (log_posterior_proposed > -Inf && log(rand()) < log_acceptance)

    if accepted
        theta = theta_prop
        loglik = loglik_prop
    end

    return (theta, loglik, accepted)
end

function xi_grid_search(w::AbstractArray{Float64}, ess_threshold = 0.5, depth = 100)
    ξ = (2.0)^(-1.0)
    for d = 2:depth
        ess = calculate_ess(w * ξ) / length(w)
        if ess < ess_threshold
            ξ -= (2.0)^(-d)
        else
            ξ += (2.0)^(-d)
        end
    end

    ess = calculate_ess(w * ξ) / length(w)
    if (ess < ess_threshold)
        ξ -= (2.0)^(-depth)
    end

    if ξ == 0
        println("Less than ess_threshold of observations are zero")
        ξ = 1e-12
    end

    return ξ, ess
end


text:
function title!(simple_plot::SimplePlot, cur_string::AbstractString)
  simple_plot.layout["title"] = cur_string
  validate!(simple_plot)
end

function xlabel!(simple_plot::SimplePlot, cur_string::AbstractString)
  simple_plot.layout["xaxis"]["title"] = cur_string
  validate!(simple_plot)
end

function ylabel!(simple_plot::SimplePlot, cur_string::AbstractString)
  simple_plot.layout["yaxis"]["title"] = cur_string
  validate!(simple_plot)
end

title!(cur_string::AbstractString) = title!(_plot, cur_string)

xlabel!(cur_string::AbstractString) = xlabel!(_plot, cur_string)
ylabel!(cur_string::AbstractString) = ylabel!(_plot, cur_string)

export title!

export xlabel!
export ylabel!


text:
module ClimateMARGO

if VERSION < v"1.5"
    @error "ClimateMARGO requires Julia v1.5 or newer."
end

include("Models/Models.jl")
include("Utils/Utils.jl")
include("Diagnostics/Diagnostics.jl")
include("Optimization/Optimization.jl")
include("IO/IO.jl")
include("PolicyResponse/PolicyResponse.jl")
if get(ENV, "JULIA_MARGO_LOAD_PYPLOT", "1") == "1"
    # default
    include("Plotting/Plotting.jl")
else
    # (in our API, we don't load the plotting functions.)
    @info "Not loading plotting functions"
end
end


text:
using PowerModelsReliability

using Logging
# suppress warnings during testing
Logging.configure(level=ERROR)

using JuMP
using PowerModels
PMs = PowerModels

using Ipopt

using Base.Test

# default setup for solvers
ipopt_solver = IpoptSolver(tol=1e-6, print_level=0)

# this will work because PowerModels is a dependency
case_files = [
    "../../PowerModelsReliability/test/data/case5_tf.m",
]

#include("prob/tfopf.jl")

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
struct VariableReference
    kind::Symbol
    value
end

mutable struct DebuggerState
    last_exception
    top_level_expressions::Vector{Any}
    current_top_level_expression::Int
    frame
    not_yet_set_function_breakpoints::Set{Any}
    debug_mode::Symbol
    compile_mode
    sources::Dict{Int,String}
    next_source_id::Int
    varrefs::Vector{VariableReference}

    function DebuggerState()
        return new(nothing, [], 0, nothing, Set{String}(), :unknown, JuliaInterpreter.finish_and_return!, Dict{Int,String}(), 1, VariableReference[])
    end
end

is_toplevel_return(frame) = frame.framecode.scope isa Module && JuliaInterpreter.isexpr(JuliaInterpreter.pc_expr(frame), :return)

function attempt_to_set_f_breakpoints!(bps)
    for bp in bps
        @debug "Trying to set function breakpoint for '$(bp.name)'."
        try
            f = Core.eval(bp.mod, bp.name)

            signat = if bp.signature !== nothing
                Tuple{(Core.eval(Main, i) for i in bp.signature)...}
            else
                nothing
            end

            JuliaInterpreter.breakpoint(f, signat, bp.condition)
            delete!(bps, bp)

            @debug "Setting function breakpoint for '$(bp.name)' succeeded."
        catch err
            @debug "Setting function breakpoint for '$(bp.name)' failed."
        end
    end
end

function get_next_top_level_frame(state)
    state.current_top_level_expression += 1

    if state.current_top_level_expression > length(state.top_level_expressions)
        return nothing
    else
        next_top_level = state.top_level_expressions[state.current_top_level_expression]
        next_frame = JuliaInterpreter.prepare_thunk(next_top_level)
        return next_frame
    end
end

function our_debug_command(cmd, state)
    empty!(state.varrefs)
    while true
        @debug "Running a new frame." state.frame state.compile_mode

        ret = Base.invokelatest(JuliaInterpreter.debug_command, state.compile_mode, state.frame, cmd, true)

        attempt_to_set_f_breakpoints!(state.not_yet_set_function_breakpoints)

        @debug "Finished running frame." ret

        if ret !== nothing && is_toplevel_return(ret[1])
            ret = nothing
        end

        if ret !== nothing
            state.frame = ret[1]
            return ret[2]
        end

        state.frame = get_next_top_level_frame(state)

        if state.frame === nothing
            return nothing
        end

        ret !== nothing && error("Invalid state.")

        if ret === nothing && (cmd == :n || cmd == :s || cmd == :finish || JuliaInterpreter.shouldbreak(state.frame, state.frame.pc))
            return state.frame.pc
        end
    end
end

function send_stopped_msg(conn, ret_val, state)
    if ret_val isa JuliaInterpreter.BreakpointRef
        if ret_val.err === nothing
            send_notification(conn, "STOPPEDBP")
        else
            state.last_exception = ret_val.err
            error_msg = Base.invokelatest(string, ret_val.err)
            send_notification(conn, "STOPPEDEXCEPTION", error_msg)
        end
    elseif ret_val isa Number
        send_notification(conn, "STOPPEDSTEP")
    elseif ret_val === nothing
        send_notification(conn, "STOPPEDSTEP")
    end
end

max_stars_repo_name: cloudgiant/julia-vscode

max_issues_repo_name: cloudgiant/julia-vscode

max_forks_repo_name: cloudgiant/julia-vscode


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

module Utils

import ACE.RPI: BasicPSH1pBasis, SparsePSHDegree, RPIBasis, get_maxn
import ACE: PolyTransform, transformed_jacobi
import ACE.PairPotentials: PolyPairBasis

# - simple ways to construct a radial basis
# - construct a descriptor
# - simple wrappers to generate RPI basis functions (ACE + relatives)

export rpi_basis, descriptor, pair_basis, ace_basis

function rpi_basis(; species = :X, N = 3,
      # transform parameters
      r0 = 2.5,
      trans = PolyTransform(2, r0),
      # degree parameters
      D = SparsePSHDegree(),
      maxdeg = 8,
      # radial basis parameters
      rcut = 5.0,
      rin = 0.5 * r0,
      pcut = 2,
      pin = 0,
      constants = false,
      rbasis = transformed_jacobi(get_maxn(D, maxdeg, species), trans, rcut, rin;
                                  pcut=pcut, pin=pin),
      # one-particle basis
      Basis1p = BasicPSH1pBasis,
      basis1p = Basis1p(rbasis; species = species, D = D) )

   return RPIBasis(basis1p, N, D, maxdeg, constants)
end

descriptor = rpi_basis
ace_basis = rpi_basis

function pair_basis(; species = :X,
      # transform parameters
      r0 = 2.5,
      trans = PolyTransform(2, r0),
      # degree parameters
      maxdeg = 8,
      # radial basis parameters
      rcut = 5.0,
      rin = 0.5 * r0,
      pcut = 2,
      pin = 0,
      rbasis = transformed_jacobi(maxdeg, trans, rcut, rin; pcut=pcut, pin=pin))

   return PolyPairBasis(rbasis, species)
end

end


text:
@testset "idpad" begin
    const idpad = Manifolds.idpad

    @test idpad([1, 2], 4) == [1, 2, 0, 0]
    @test idpad([1 2; 3 4], 3) == [1 2 0;
                                   3 4 0;
                                   0 0 1]
    @test idpad([1 2 3; 4 5 6], (3, 4)) == [1 2 3 0;
                                            4 5 6 0;
                                            0 0 0 1]
    @test idpad([1 2 3; 4 5 6], (4, 3)) == [1 2 3;
                                            4 5 6;
                                            0 0 0;
                                            0 0 0]

    sizes = 2.^(1:round(Int, sqrt(nruns)))
    for s1 in sizes, s2 in sizes
        v = rand(s1)
        m = rand(s1, s1)
        padv = idpad(v, s2)
        padm = idpad(m, s2)
        if s1 > s2
            @test padv == v[1:s2]
            @test padm == m[1:s2, 1:s2]
        elseif s1 < s2
            @test padv == [v; zeros(s2-s1)]
            @test padm == [m zeros(s1, s2-s1);
                           zeros(s2-s1, s1) I]
        else
            @test padv == v
            @test padm == m
        end

        @test size(padv) == (s2,)
        @test size(padm) == (s2, s2)

        @test IndexStyle(padv) == IndexStyle(typeof(padv)) == IndexLinear()
        @test IndexStyle(padm) == IndexStyle(typeof(padm)) == IndexCartesian()
    end
end

@testset "Frame" begin
    # dodaj test: tnbframe(Interval(), 0.) * [1,2,3,4,5]
    const Frame = Manifolds.Frame

    v3 = rand(3)
    v5 = rand(5)

    # Identities:
    @test Frame(eye(1), [0.0])      * v3 == v3
    @test Frame(eye(100), [0.0])    * v3 == v3
    @test Frame(eye(1), zeros(100)) * v3 == v3

    fr  = Frame(eye(3)[randcycle(3), :], zeros(3))
    id1 = Frame(eye(1), [0.0])
    id2 = Frame(eye(10), [0.0])
    id3 = Frame(eye(10), zeros(10))
    @test id1 * fr == fr * id1 == fr
    @test id2 * fr == fr * id2 == fr
    @test id3 * fr == fr * id3 == fr

    t = rand()
    tnbframe(Interval(), t) * v3 == [t; v3]
    tnbframe(Interval(), t) * v5 == [t; v5]
end


text:
import Dagger.Sch: SchedulerOptions, ThunkOptions

function inc(x)
    x+1
end
@everywhere begin
function checkwid(x...)
    @assert myid() == 1
    return 1
end
function checktid(x...)
    @assert Threads.threadid() != 1 || Threads.nthreads() == 1
    return 1
end
end

@testset "Scheduler" begin
    @testset "Scheduler options: single worker" begin
        options = SchedulerOptions(;single=1)
        a = delayed(checkwid)(1)
        b = delayed(checkwid)(2)
        c = delayed(checkwid)(a,b)

        @test collect(Context(), c; options=options) == 1
    end
    @testset "Thunk options: single worker" begin
        options = ThunkOptions(;single=1)
        a = delayed(checkwid; options=options)(1)

        @test collect(Context(), a) == 1
    end
    @static if VERSION >= v"1.3.0-DEV.573"
        if Threads.nthreads() == 1
            @warn "Threading tests running in serial"
        end
        @testset "Scheduler options: threads" begin
            options = SchedulerOptions(;threads=true)
            a = delayed(checktid)(1)
            b = delayed(checktid)(2)
            c = delayed(checktid)(a,b)

            @test collect(Context(), c; options=options) == 1
        end
        @testset "Thunk options: threads" begin
            options = ThunkOptions(;threads=true)
            a = delayed(checktid; options=options)(1)

            @test collect(Context(), a) == 1
        end
    end
end


text:
module geometry_msgs

using Cxx
using CxxStd: WrappedCppPrimArray

cxx"""
#include <geometry_msgs/Vector3.h>
#include <geometry_msgs/Point.h>
#include <geometry_msgs/Quaternion.h>
"""

TVals = Union{}
TRefs = Union{}
TPtrs = Union{}

for dtype in [ :Vector3, :Point, :Quaternion ]
    cppname  = string("geometry_msgs::", dtype)
    cxxtval  = Expr(:macrocall, Symbol("@cxxt_str"), cppname)
    cxxtref  = Expr(:macrocall, Symbol("@cxxt_str"), string(cppname, "&"))
    cxxtptr  = Expr(:macrocall, Symbol("@cxxt_str"), string(cppname, "*"))

    symtval  = Symbol(dtype, :Msg)
    symtref  = Symbol(dtype, :Msg, :R)
    symtptr  = Symbol(dtype, :Msg, :P)
    symtall  = Symbol(dtype, :Msg, :U)

    @eval begin
        global const $(symtval) = $(cxxtval)
        global const $(symtref) = $(cxxtref)
        global const $(symtptr) = $(cxxtptr)
        global const $(symtall) = Union{$(cxxtval), $(cxxtref)}
        #export $(symtval), $(symtref), $(symtall)
    end

    TVals = Union{TVals, @eval $(symtval)}
    TRefs = Union{TRefs, @eval $(symtref)}
    TPtrs = Union{TPtrs, @eval $(symtptr)}

    body = Expr(:macrocall, Symbol("@icxx_str"), string(cppname, "();"))
    @eval (::Type{$(symtval)})() = $(body)

    body = Expr(:macrocall, Symbol("@icxx_str"), string("(", cppname, "&) *\$v;"))
    @eval Base.convert(::Type{$(cxxtref)}, v::$(cxxtptr)) = $(body)
end

function Base.convert(::Type{CT}, v::T) where {CT<:AbstractArray,T<:Union{TVals,TRefs,TPtrs}}
    convert(CT, unsafe_wrap(DenseArray, v))
end

Base.unsafe_wrap(::Type{DenseArray}, v::T) where {T<:TPtrs} = unsafe_wrap(DenseArray, icxx"*$v;")

function Base.unsafe_wrap(::Type{DenseArray}, v::T) where {T<:Union{TVals,TRefs}}
    if typeof(icxx"$v;") <: QuaternionMsgU
        WrappedCppPrimArray(icxx"&($v.x);",4)
    else
        WrappedCppPrimArray(icxx"&($v.x);",3)
    end
end

end #module

eval(geometry_msgs, :(Vector3     = Vector3Msg))
eval(geometry_msgs, :(Vector3R    = Vector3MsgR))
eval(geometry_msgs, :(Vector3P    = Vector3MsgP))
eval(geometry_msgs, :(Vector3U    = Vector3MsgU))
eval(geometry_msgs, :(Quaternion  = QuaternionMsg))
eval(geometry_msgs, :(QuaternionR = QuaternionMsgR))
eval(geometry_msgs, :(QuaternionP = QuaternionMsgP))
eval(geometry_msgs, :(QuaternionU = QuaternionMsgU))

export geometry_msgs

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS