text:
module Connectors

using Logging
using ..LogCompose
import ..LogCompose: logcompose, log_min_level

#--------------------------------------------------------------------
# The connectors for individual logger implementations should ideally
# be in their own packages. Here we implement connectors only for the
# loggers in the stdlib Logging package.
#--------------------------------------------------------------------

logcompose(::Type{Logging.NullLogger}, config::Dict{String,Any}, logger_config::Dict{String,Any}) = Logging.NullLogger()

function logcompose(::Type{Logging.SimpleLogger}, config::Dict{String,Any}, logger_config::Dict{String,Any})
    level = log_min_level(logger_config, "Info")

    streamname = strip(get(logger_config, "stream", "stdout"))
    @assert !isempty(streamname)

    stream = streamname == "stdout" ? stdout :
             streamname == "stderr" ? stderr :
             open(streamname, "a+")

    Logging.SimpleLogger(stream, level)
end

function logcompose(::Type{Logging.ConsoleLogger}, config::Dict{String,Any}, logger_config::Dict{String,Any})
    level = log_min_level(logger_config, "Info")

    streamname = strip(get(logger_config, "stream", "stdout"))
    @assert !isempty(streamname)

    stream = streamname == "stdout" ? stdout :
             streamname == "stderr" ? stderr :
             open(streamname, "a+")

    color = get(logger_config, "color", nothing)
    color === nothing || (stream = IOContext(stream, :color=>color))

    displaysize = get(logger_config, "displaysize", nothing)
    if displaysize !== nothing
        if !(displaysize isa AbstractVector) || length(displaysize) != 2
            error("Expected [height,width] but got displaysize=$displaysize")
        end
        stream = IOContext(stream, :displaysize=>Tuple(convert(Vector{Int},displaysize)))
    end

    show_limited = get(logger_config, "show_limited", true)
    color isa Bool || error("Expected boolean but got show_limited=$show_limited")

    Logging.ConsoleLogger(stream, level; show_limited=show_limited)
end

end # module Connectors


text:
"""
    SeisRadon_tx(In,Par,itype);
RADON_FX: Operators for f-x, tau-p forward and adjoint Radon transforms
          to compute linear or  parabolic Radon transform.
IN    Radon coefficients if itype =  1 (Forward transform)  with In(np,nt)
      CMP gather         if itype = -1 (Adjoint transform)  with In(nh,nt)

OUT   CMP gather         if itype =  1 (Forward transform)  with Out(nh,nt)
      Radon coeffcients  if itype = -1 (Adjoint transform)  with Out(np,nt)
# Parameters
- `Par.h`:vector containing the nh offsets
- `Par.p`:vector containing the np curvatures normalized by far offset (s) if itype='parab'
- `Par.p`:vector containing the np dips in s/m if itype='linear'
- `Par.dt`:sampling interval
- `Par.f`:frequency corners of BP operator - this acts like a zero phase wavelet.
- `Par.transf `:'parab', 'linear', 'hyperb'
itype =  1  :  forward
itype = -1  :  adjoint
"""
function SeisRadon_tx(In,Par,itype);
    h = Par[1][:h];
    p = Par[2][:p];
    dt = Par[3][:dt];
    f  = Par[4][:f];
    f1 = f[1];
    f2 = f[2];
    f3 = f[3];
    f4 = f[4];
    np = length(p);
    nh = length(h);
    if itype == 1;
       m = In;
       (nt, np) = size(m);
       d = zeros(nt,nh);
       m = SeisBPFilter(m,dt,f1,f2,f3,f4);
    end;
    if itype == -1;
       d = In;
       (nt, nh) = size(d);
       m = zeros(nt,np);
    end;
   tmax = (nt-1)*dt;
   tau = collect(0:1:nt-1)*dt;

     for itau = 1:1:nt
       for ih = 1:nh
        for ip = 1:np
           if cmp(Par[5][:transf],"parab")==0 ;
            hmax = maximum(abs.(h));
            t=tau[itau]+ (h[ih]/hmax)^2*p[ip];
         elseif cmp(Par[5][:transf],"linear")==0;
            t=tau[itau] + (h[ih]*p[ip]);
         elseif cmp(Par[5][:transf],"hyperb")==0;
            t=sqrt(tau[itau]^2 + (h[ih]^2/p[ip]^2));
         end
          if t>0 && t<tmax
            it    = t/dt+1;
            it1   = floor(Int,it);
            it2   = it1 + 1;
            alpha = it  - it1;
               if itype == 1;
                    d[it1,ih] = d[it1,ih]+ (1-alpha) * m[itau,ip];
                    d[it2,ih] = d[it2,ih] +     alpha * m[itau,ip];
                 else
                    m[itau,ip] = m[itau,ip] + (1-alpha) * d[it1,ih] + alpha*d[it2,ih];

          end
        end
       end
      end
   end
   if itype == 1;
      Out = d;
   end;

   if itype == -1;
      m = SeisBPFilter(m,dt,f1,f2,f3,f4);
      Out = m;
   end;
   return Out
end

ext: jl


text:
module TestCxxType
using Test
using CxxInterface

@testset "C++ types" begin
    ctypesymbols = [:Cstring, :Cuchar, :Cuint, :Cchar, :Cdouble, :Cfloat, :Cvoid, :Cwchar_t, :Cint, :Cptrdiff_t, :Clong, :Clonglong,
                    :Cssize_t, :Culong, :Csize_t, :Cshort, :Cwstring, :Culonglong, :Cushort]
    ctypes = Dict(sym => eval(sym) for sym in ctypesymbols)
    for T in values(ctypes)
        @test haskey(cxxtype, T)
    end

    libRoundTrip = joinpath(pwd(), "libRoundTrip")

    eval(cxxsetup())
    eval(cxxnewfile("RoundTrip.cxx", ""))

    make_arg(T::Type) = FnArg(:x, T, "x", cxxtype[T])
    make_body(T::Type) = "return x;"
    # void is not a valid argument type
    make_arg(T::Type{Cvoid}) = FnArg(:x, Cint, "x", "int", Cvoid, _ -> 1)
    make_body(::Type{Cvoid}) = "return;"
    for (sym, T) in pairs(ctypes)
        ex = cxxfunction(FnName(:roundtrip, "roundtrip_$sym", libRoundTrip), FnResult(T, cxxtype[T]), [make_arg(T)], make_body(T))
        eval(ex)
    end

    cxx_write_code!()
    using Libdl: dlext
    run(`c++ -fPIC -shared -o libRoundTrip.$dlext RoundTrip.cxx`)

    example(T) = one(T)
    example(::Type{Cvoid}) = Cvoid()
    example(T::Type{Cstring}) = Base.cconvert(Cstring, C_NULL)
    example(T::Type{Cwstring}) = Base.cconvert(Cwstring, C_NULL)
    if Sys.ARCH ≡ :x86_64
        for T in values(ctypes)
            x = example(T)::T
            @test roundtrip(x) === x
        end
    end
end

end#module


text:
using Test;

function format_re(s)
    out = ""
    indent = ""
    for c in s
        out *= c
        if c == '('
            indent *= "  "
            out *= ("\n" * indent)
        elseif c == ')'
            indent = indent[1:length(indent) - 2]
            out *= ("\n" * indent)
        end
    end
    return out
end

function expand(rules, num, version=1)::String
    # base case
    if typeof(rules[num]) == String
        # println(offset, rules[num])
        return rules[num]
    else
    # recursive case
        options = []
        if version == 2
            rules[11] = [[42],[31]]
        end
        for choice in rules[num] # things separated by '|'
            val = ""
            for r in choice # numbers in the rule
                val *= expand(rules, r, version)
            end
            push!(options, val)
        end
        out = join(options, "|")
        if version == 2
            if num == 8
                return "($out)+"
            elseif num == 11
                a = options[1]
                b = options[2]
                # ab,aabb,aaabbb,aaaabbbb,...
                # 4 seems to be enough to get all answers
                return "((($a)($b))|(($a){2}($b){2})|(($a){3}($b){3})|(($a){4}($b){4}))"
            else
                return "($out)"
            end
        end
        return "($out)"
    end
end

function parse_input(lines, version=2)
    messages = []
    rules = Dict()
    for l in lines
        if l == ""
            continue
        elseif ':' in l
            if version == 2
                # part 2: add repeats
                # if startswith(l, "8: ")
                #     # same number repeated over and over
                #     l = "8: 42 | 42 42 | 42 42 42 | 42 42 42 42 | 42 42 42 42 42 "
                # if startswith(l, "11: ")
                # #     # (a^x)(b^x) where x is the same number... is hard to represent in regex, so just keep extending until I find all the hits
                #     l = "11: 42 31 | 42 42 31 31 | 42 42 42 31 31 31"
                #     # l = "11: 42 31 | 42 42 31 31 | 42 42 42 31 31 31 | 42 42 42 42 31 31 31 31"
                # #     # l = "11: 42 31 | 42 42 31 31 | 42 42 42 31 31 31 | 42 42 42 42 31 31 31 31| 42 42 42 42 42 31 31 31 31 31"
                # end
            end
            num, rule = parse_rule(l)
            rules[num] = rule
        else
            push!(messages, l)
        end
    end
    return rules, messages
end

function parse_rule(l)
    parts = split(l, ": ")
    num = parse(Int, parts[1])
    if in('"', parts[2])
        # it's just a char
        return num, replace(parts[2], "\"" => "")
    end

    out = []
    for choice in split(parts[2], "|")
        items = split(choice)
        nums = map((x) -> parse(Int, x), items)
        push!(out, nums)
    end
    return num, out
end

function run(fname, version=1)
    println("[part ", version, "][START] Running file = ", fname)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    # Parse input
    rules, messages = parse_input(lines, version)
    rule_regex = expand(rules, 0, version)
    println("Rule:")
    println("")
    println(rule_regex)
    # println(format_re(rule_regex))
    println("")
    re = Regex("^$(rule_regex)\$")
    result = 0
    for m in messages
        if occursin(re, m)
            # println("match:    ", m)
            result += 1
        else
            # println("no match: ", m)
        end
    end

    println("[part ", version, "][END] Running file = ", fname, " ... got result = ", result)
    return result
end

# example_rules = """
# 0: 4 1 5
# 1: 2 3 | 3 2
# 2: 4 4 | 5 5
# 3: 4 5 | 5 4
# 4: "a"
# 5: "b"
# """

# expected_rules = Dict()
# expected_rules[0] = [[4,1,5]]
# expected_rules[1] = [[2,3], [3,2]]
# expected_rules[2] = [[4,4], [5,5]]
# expected_rules[3] = [[4,5], [5,4]]
# expected_rules[4] = "a"
# expected_rules[5] = "b"

# actual_rules, _ = parse_input(split(example_rules, "\n"))
# @test actual_rules == expected_rules

# println("== START ==")
# println(expand(actual_rules, 0))
# println("== END ==")

# @test run("19ex.txt") == 2

# run("19.txt")

# @test run("19ex.part2.txt") == 3
# @test run("19ex.part2.txt", 2) == 12

run("19.txt", 2)


text:
"""
    setCrossoverOperator(method::Function, varArgs...;
                        genj::GeneticJulia = GenJ,
                        probability::Number = 0.7,
                        nChildren::Int = -1)

Set the crossover operator of parents for `genj`, receiving a function that cross multiples individual representations and return one child or multiple children.
    This function constructs an structure of type `CrossoverOperator` and adds it to the main structure.

# Arguments
- `method::Function`: function that receives, first, a list of parents as arguments which are individual representation and `"varArgs"`, that is a list of arguments, as own arguments, and return a number between 1 and nParents of children, where each one is an individual representation. (For information about representation look [`IndividualRep`](@ref)). See `"Examples"`.
- `varArgs`: List of arguments for crossover method. See `"Examples"`.

# Keyword Arguments
- `genj::GeneticJulia = GenJ`: the main structure. For `code user`: don't modify the default value unless multiple experiments are goint to be run.
- `probability::Number = 0.7`: A floating point number, between `0` and `1`, that represents the probability of doing each single crossover operation.
- `nChildren::Int = -1`: The number of children as result of the crossover method. It must be either a number between `0` and `nParents` or `-1` which is equal to `nParents`.

## Note
The `CrossoverOperator` uses the `nParents` variable that define the number of parents per cross and it can be set through `setNParentsPerCross` function. The default value is `2`.

# Examples
```jldoctest
#Function with 2 parents, no more arguments and 1 child
function ExampleCrossParents1(Parent1::IndividualRep, Parent2::IndividualRep)
    #···Code···
    return children::IndividualRep # Type expressed here for user awareness
end

julia> setCrossoverOperator(ExampleCrossParents1, probability = 0.4, nChildren = 1)

_______________________________

#Function with more than 2 parents, no more arguments and the same number of children as parents.
function ExampleCrossParents2(Parent1::IndividualRep, Parent2::IndividualRep, Parent3::IndividualRep)
    #···Code···
    return children::Array{IndividualRep}(undef, 3) # Type expressed here for user awareness
end

julia> setCrossoverOperator(ExampleCrossParents2, probability = 0.6, nChildren = -1)

________________________________

#Function with more arguments and the same number of children as parents
function ExampleCrossParents3(Parent1::IndividualRep, Parent2::IndividualRep, arg1, arg2)
    #···Code···
    return children::Array{IndividualRep}(undef, 2) # Type expressed here for user awareness
end

julia> setCrossoverOperator(ExampleCrossParents3, arg1, arg2, probability = 0.6, nChildren = -1)
```

See also: [`CrossoverOperator`](@ref)
"""
function setCrossoverOperator(method::Function, varArgs...; genj::GeneticJulia = GenJ, probability::Number = 0.7, nChildren::Int = -1)

    if !(0 <= probability <= 1)
        error("Crossover probability must be within 0 and 1")
    end

    if nChildren != -1 && nChildren < 1
        error("Number of children per cross must be either -1 or greater than 0")
    end

    varArgs = [i for i in varArgs]
    genj._crossoverOp = CrossoverOperator(method, probability, nChildren, varArgs)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(BiPerplex(0, 0.0))
end

@test begin
    !iszerodivisor(BiPerplex(0,1,2,0))
end

@test begin
    iszerodivisor(BiPerplex(1,2,3,4))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(BiPerplex(a))
end

@test begin
    a = 1
    !isreal(BiPerplex(a, a))
end

@test begin
    z = random(BiPerplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(BiPerplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, BiPerplex(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, R: 2, Q: 3, QR: 4]"
    l == r
end

@test begin
    x = random(BiPerplex{BigFloat})
    y = random(BiPerplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(BiPerplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(BiPerplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Prplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    barR(barR(x)) == x
end

@test begin
    x = random(BiPerplex{BigInt})
    barQ(barQ(x)) == x
end

@test begin
    x = random(BiPerplex{BigInt})
    barQ(barR(x)) == barR(barQ(x))
end

@test begin
    x = random(BiPerplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barR(x * y)
    r = barR(x) * barR(y)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barR(x + y)
    r = barR(x) + barR(y)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barR(x - y)
    r = barR(x) - barR(y)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barQ(x * y)
    r = barQ(x) * barQ(y)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barQ(x + y)
    r = barQ(x) + barQ(y)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = barQ(x - y)
    r = barQ(x) - barQ(y)
    l == r
end

@test begin
    x = random(BiPerplex{Rational{BigInt}})
    y = random(BiPerplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(BiPerplex{Rational{BigInt}})
    y = random(BiPerplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(BiPerplex{Rational{BigInt}})
    y = random(BiPerplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(BiPerplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    n = norm(x)
    n > big(0) || n < big(0) || n == big(0)
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(BiPerplex(0, 0.0))
end

@test_throws ErrorException begin
    random(BiPerplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(BiPerplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(BiPerplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(BiPerplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(BiPerplex{BigInt})
    z == z + zero(BiPerplex{BigInt})
end

@test begin
    z = random(BiPerplex{BigInt})
    z == z * one(BiPerplex{BigInt})
end

@test begin
    z = random(BiPerplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(BiPerplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = commutator(x, y)
    r = zero(BiPerplex{BigInt})
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = alternatorL(x, y)
    r = zero(BiPerplex{BigInt})
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = alternatorR(x, y)
    r = zero(BiPerplex{BigInt})
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    l = flexator(x, y)
    r = zero(BiPerplex{BigInt})
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = associator(x, y, z)
    r = zero(BiPerplex{BigInt})
    l == r
end

@test begin
    x = random(BiPerplex{BigInt})
    y = random(BiPerplex{BigInt})
    z = random(BiPerplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = BiPerplex(Prplex(1,2))
    r = BiPerplex(1,2,0)
    l == r
end

@test begin
    l = BiPerplex(1, 2, 3.0)
    r = BiPerplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = BiPerplex(1, 2, 3, 4.0)
    r = BiPerplex(1.0, 2.0, 3.0, 4.0)
    l == r
end


text:
mutable struct DynamicSimulation
    dyn_system::DynamicSystem
    problem::DiffEqBase.DAEProblem
    callbacks::DiffEqBase.DiscreteCallback
    tstops::Vector{Float64}
    x0_init::Vector{Float64}
    initialized::Bool
    solution::Union{Nothing, DiffEqBase.DAESolution}
end

function DynamicSimulation(dyn_system::DynamicSystem,
                           tspan,
                           control,
                           callback,
                           x0_init)

    if !is_indexed(dyn_system)
        _index_dynamic_system!(dyn_system)
    end

    dx0 = zeros(get_total_rows(dyn_system))
    prob = DiffEqBase.DAEProblem(system_model!,
                              dx0,
                              x0_init,
                              tspan,
                              (control, dyn_system),
                              differential_vars = dyn_system.DAE_vector)

    return DynamicSimulation(dyn_system,
                             prob,
                             callback,
                             [1.0],
                             x0_init,
                             true,
                             nothing)

end

function run_simulation!(sim::DynamicSimulation, solver; kwargs...)
    sim.solution = DiffEqBase.solve(sim.problem,
                           solver;
                           callback = sim.callbacks,
                           tstops = sim.tstops, kwargs...)

    return
end


text:
### A Pluto.jl notebook ###
# v0.17.7

using Markdown
using InteractiveUtils

# ╔═╡ 4df3170c-83a7-11ec-0eb7-9fe0a20bfc16
md"""
# Introduction to Wordle

[Wordle](https://www.powerlanguage.co.uk/wordle/) is a daily word guessing game where the objective is to guess the daily, five-letter "wordle" within six guesses. Each time you make a guess, you are given the following feedback:
- Letters that are in the correct position will be highlighted in green;
- Letters that are in the wordle but not in the correct position are highlighted in yellow;
- And, the rest are highlighted in gray.

Here are the instructions from the [game](https://www.powerlanguage.co.uk/wordle/):
![instructions-screenshot](https://images.squarespace-cdn.com/content/v1/50eca855e4b0939ae8bb12d9/80a523b2-edff-41ac-8a9e-d515285a0b74/How+to+play+original+Wordle.png?format=1000w)

The game went viral as friends would compete with each other to see who could guess the wordle with the fewest number of guesses. [Try it out for yourself!](https://www.powerlanguage.co.uk/wordle/).
"""

# ╔═╡ a97ec29b-bdc6-4015-940e-a84fe8001e01
md"""
## Guessing Strategy

When starting, one might expect that any given word is as good as any other as a first guess. However, this is not the case as certain letters are more frequent in English than others. For example, most words have at least one vowel and consonants like `s` and `r` are more common than others (see [this Wikipedia page](https://www.wikiwand.com/en/Letter_frequency) for more information about this).

Believe it or not, there is a significant amount of nerdiness 🤓 and snobbery 😎 surrounding which words are the best ones to guess first. For example, many have suggested that you should prioritize words that have the most frequent letters, including the author of [this Medium article](https://towardsdatascience.com/a-frequency-analysis-on-wordle-9c5778283363) and many others on Twitter. But is this really the best strategy? Let's find out.
"""

# ╔═╡ 9b5ff53f-0d0f-4aa8-8484-25b5d4ecdc3f
md"""
## Fetching the Word Lists

Before we begin to try different guessing strategies, we need to first find some words to try. It turns out that the full lists of valid words and wordles can be found in the game source code. Once extracted, we will be able to use these lists of words to validate strategies.

However, we should take care in the way we do this extraction to avoid spoilers. Ideally, the code we write to do this extraction requires no explicit knowledge of past or future wordles. In this way, no wordles are specifically "hard-coded" into the extraction routines. Today's wordle (whenever you are reading this) is also not considered to be "known". Additionally, in this spirit, we should only use historical wordles to do our validation.

I won't include this fetching code in this notebook as it is not super interesting. If you are interested, you can find it in the `src/wordle.jl` file or [on GitHub](https://github.com/jowch/wordle-entropy/blob/main/src/wordle.jl).

With these caveats, let's get started! 🚀
"""

# ╔═╡ Cell order:
# ╟─4df3170c-83a7-11ec-0eb7-9fe0a20bfc16
# ╟─a97ec29b-bdc6-4015-940e-a84fe8001e01
# ╟─9b5ff53f-0d0f-4aa8-8484-25b5d4ecdc3f

max_stars_repo_name: jowch/wordle-entropy

max_issues_repo_name: jowch/wordle-entropy

max_forks_repo_name: jowch/wordle-entropy


text:
function extract_equation_data(obj::EquationsMapper, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "extractEquationData", Vector{jdouble}, (jint, Vector{jdouble}), arg0, arg1)
end

function get_number_of_equations(obj::EquationsMapper)
    return jcall(obj, "getNumberOfEquations", jint, ())
end

function get_total_dimension(obj::EquationsMapper)
    return jcall(obj, "getTotalDimension", jint, ())
end

function insert_equation_data(obj::EquationsMapper, arg0::jint, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return jcall(obj, "insertEquationData", void, (jint, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function map_state_and_derivative(obj::EquationsMapper, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return jcall(obj, "mapStateAndDerivative", ODEStateAndDerivative, (jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end


text:
function SPICE2Matrix(FileName::String)
    SPICE_DF = LTSpiceLoad(FileName)
    NodeList,InputList,NumVSources = ProcessSPICE_DF(SPICE_DF)
    return SPICE_DF,NodeList,InputList,NumVSources
end

function ProcessSPICE_DF(DF::DataFrame)

    NodeList = []
    NodeList = reverse(sort(unique(vcat(DF[:,3],DF[:,4])))) #picking out unique nodes, and making 0 at the end

    InputList = NodeList;
    InputList = convert(Array{String,1}, InputList)
    InputList = "Iin_".*InputList

    NumVSources=0
    for k in 1:length(DF.Type)
           if DF.Type[k]=='V'
               NumVSources = NumVSources+1
               push!(InputList, "Vin_".*DF.Name[k])
           end
    end
    for j in 1:length(NodeList)

            DF[:,3] = replace(DF[:,3],NodeList[j]=>j)
            DF[:,4] = replace(DF[:,4],NodeList[j]=>j)

    end
    DF[:,3] = replace(DF[:,3],length(NodeList)=>999) #Making GND node 999
    DF[:,4] = replace(DF[:,4],length(NodeList)=>999)

    filter!(x-> x!="0",NodeList) #Removing GND node
    filter!(x-> x!="Iin_0",InputList)

return NodeList,InputList,NumVSources
end

function addG!(G,Value,N1,N2)

    Z=Value
    if N1<999
        G[N1,N1] = G[N1,N1]+ 1 ./Z
    end
    if N2<999
        G[N2,N2] = G[N2,N2]+ 1 ./Z
    end
    if ((N1<999)&(N2<999))
        G[N1,N2] = G[N1,N2]- 1 ./Z
        G[N2,N1] = G[N2,N1]- 1 ./Z
    end

    return G
end

function SPICE_DF2Matrix_ω(DF,ω,InputList)

    Y = complex(zeros(length(InputList),length(InputList)))

    SrcMat = zeros(size(Y))

 for i in 1:size(DF,1)
     N1 = DF.Node1[i]
     N2 = DF.Node2[i]
     if (DF.Type[i]=='R')|(DF.Type[i]=='L')|(DF.Type[i]=='C')
         #println(i)

         if DF.Type[i]=='R'
             Z = DF.Value[i]

             addG!(Y,Z,N1,N2)
         elseif DF.Type[i]=='L'
             Z = im*ω*DF.Value[i]
             if DF.ESR[i]>0
                 Z += DF.ESR[i]
             end
             addG!(Y,Z,N1,N2)
         elseif DF.Type[i]=='C'
            Z = 1 /(im*ω*DF.Value[i])
            if DF.ESR[i]>0
                Z += DF.ESR[i]
            end
            addG!(Y,Z,N1,N2)
         end

     elseif DF.Type[i]=='V'

         VecLoc = findfirst( x -> x=="Vin_"*DF.Name[i], InputList)

         if N1<999
             SrcMat[VecLoc,N1] = 1
             SrcMat[N1,VecLoc] = 1

         end
         if N2<999

             SrcMat[VecLoc,N2] = -1
             SrcMat[N2,VecLoc] = -1
         end

     end

 end
 FullMat = Y .+ SrcMat

 return FullMat

end


text:
"""
    gauss_lobatto_quad(α, β, N)

Computes Legendre-Gauss-Lobatto quadrature points and weights with Jacobi weights α,β.
"""
function gauss_lobatto_quad(α::Int, β::Int, N)

    if !(iszero(α) && iszero(β))
        error("alpha/beta must be zero for Lobatto")
    end

    x = zeros(N+1)
    w = zeros(N+1)
    if N == 0
        x[1] = 0.
        w[1] = 2.
    elseif N == 1
        x[1] = -1.0
        x[2] = 1.0
        w[1] = 1.0
        w[2] = 1.0
    elseif N==2
        x[1] = -1.
        x[2] = 0.
        x[3] = 1.
        w[1] = 0.333333333333333
        w[2] = 1.333333333333333
        w[3] = 0.333333333333333
    else
        xint, wint = gauss_quad(α+1, β+1, N-2)
        x = vcat(-1,xint,1)
        wend = 2/(N*(N+1))
        w = vcat(wend, (@. wint / (1-xint^2)) ,wend)
    end
    return x[:],w[:]
end

"""
    gauss_quad(α, β, N)

Compute nodes and weights for Gaussian quadrature with Jacobi weights (α,β).
"""
function gauss_quad(α, β, N)
    if N == 0
        x = [-(α-β)/(α+β+2)]
        w = [2.]
        return x, w
    end

    J = zeros(N+1, N+1)
    h₁ = @. 2*(0:N)+α+β
    J = diagm(0 => @. -1/2*(α^2-β^2)/(h₁+2)/h₁) +
        diagm(1 => @. 2/(h₁[1:N]+2) * sqrt((1:N) * ((1:N)+α+β) * ((1:N)+α) * ((1:N)+β) / (h₁[1:N]+1) / (h₁[1:N]+3)))
    if α + β < 10*eps()
        J[1,1] = 0.0
    end

    x, V = eigen(Symmetric(J+J'))
    w = @. transpose(V[1,:])^2 * 2^(α+β+1) / (α+β+1) * gamma(α+1) * gamma(β+1) / gamma(α+β+1)
    return x[:], w[:]
end

"""
    jacobiP(x, α, β, N, tmp_array=zeros(N+1))
    jacobiP(x::AbstractArray, α, β, N)

Evaluate the Jacobi Polynomial (α, β) of order N at points x.
Can optionally specify `tmp_array` of length > N+1 to avoid allocations.
"""
function jacobiP(x, α, β, N, tmp_array=zeros(N+1))

    PL = tmp_array
    γ₀ = 2^(α+β+1) / (α+β+1) * gamma(α+1) * gamma(β+1) / gamma(α+β+1)
    PL[1] = 1.0 / sqrt(γ₀)
    if N == 0
        return PL[1]
    end

    γ₁ = (α+1) * (β+1) / (α+β+3) * γ₀
    PL[2] = ((α+β+2) * x/2 + (α-β)/2) / sqrt(γ₁)
    if N == 1
        return PL[N+1]
    end

    aold = 2 / (2+α+β) * sqrt((α+1) * (β+1) / (α + β + 3))

    for i = 1:N-1
        h₁ = 2i + α + β
        anew = 2 / (h₁+2) * sqrt((i+1) * (i+1+α+β) * (i+1+α) * (i+1+β) / (h₁+1) / (h₁+3))
        bnew = -(α^2-β^2) / h₁ / (h₁+2)
        PL[i+2] = 1 / anew * (-aold * PL[i] + (x - bnew) * PL[i+1])
        aold = anew
    end

    return PL[N+1]
end

# vectorized version.
# TODO: deprecate?
function jacobiP(x::AbstractArray, α, β, N)

    PL = zeros(N+1, length(x))
    γ₀ = 2^(α+β+1) / (α+β+1) * gamma(α+1) * gamma(β+1) / gamma(α+β+1)
    @. PL[1,:] = 1.0/sqrt(γ₀)
    if N == 0
        return vec(PL)
    end

    γ₁ = (α+1) * (β+1) / (α+β+3) * γ₀
    @. PL[2, :] = ((α+β+2) * x/2 + (α-β)/2) / sqrt(γ₁)
    if N == 1
        return PL[N+1,:]
    end

    aold = 2 / (2+α+β) * sqrt((α+1) * (β+1) / (α + β + 3))

    for i = 1:N-1
        h₁ = 2i+α+β
        anew = 2 / (h₁+2) * sqrt((i+1) * (i+1+α+β) * (i+1+α) * (i+1+β) / (h₁+1) / (h₁+3))
        bnew = -(α^2-β^2) / h₁ / (h₁+2)
        @. PL[i+2,:] = 1 / anew * (-aold * PL[i,:] .+ (x .- bnew) .* PL[i+1,:])
        aold = anew
    end

    return PL[N+1,:]
end

"""
    grad_jacobiP(r, α, β, N, tmp_array=())

Evaluate derivative of Jacobi polynomial (α, β) of order N at points r.
Pass in `tmp_array = tuple(zeros(N+1))` to avoid allocations.
"""

function grad_jacobiP(r, α, β, N, tmp_array=())
    if N != 0
        return sqrt(N * (N+α+β+1)) * jacobiP(r, α+1, β+1, N-1, tmp_array...)
    else
        return zero(r)
    end
end


text: include("html/HMTLFormat.jl")


text:
#################################
# Imagery API                   #
#################################
# Get radar and satellite images.

# Contains a centralized function to make calls to the API.
include("request.jl")

function get_radar_and_satellite_imagery(apikey::String, locationKey::String, resolution::String, language::String)
    # Returns 1/5/10/15 day(s) of weather alarms for a specific location.

    url = "http://dataservice.accuweather.com/imagery/v1/maps/radsat/"

    # Add required parameters to the URL
    url = url * resolution * "/"
    # Allowed resolutions: '480x480', '640x480' and '1024x1024'.

    url = url * locationKey
    url = url * "?apikey=" * apikey

    # Add optional parameters to the URL
    if language != Nothing url = url * "&language=" * language end

    return make_API_call(url)
end


text:
using GeneSets
using Test

@testset "GeneSets.jl" begin
    # Write your tests here.
end


text:
export HTTPRequest,
    HTTPResponse,
    HTTP_1_1

function preparehttpresponseheaders(r::Response)
    h_chunk, h_chunk_value = r.protocol.header_msg_chunk

    if chunked(r)
        if haskey(r.headers, r.protocol.header_msg_size)
            throw(ResponseException("invalid_header",
                str("The '", r.protocol.header_msg_size,
                    "' header doesn't valid on chunked data.");
                header=r.protocol.header_msg_size)
            )
        end
    end

    if r.status in r.protocol.status_wo_msg
        h = chunked(r) ? h_chunk : r.protocol.header_msg_size
        if haskey(r.headers, h)
            throw(ResponseException("status_not_accept_msg_body",
                str("The Response status ", r.status,
                    "' doesn't accept response body.");
                header=h, status=r.status))
        end
    end
end

const HTTP_1_1 = Protocol(
    "HTTP",
    "1.1",
    STATUS_CODES,
    sort(filter(x -> (x > 90 && x < 200) || x in (204, 304),
        [_ for _ in keys(STATUS_CODES)]
    )),
    String[],
    preparehttpresponseheaders,
    DefaultStatus(200, 400, 404, 500),
    ("Transfer-Encoding", "chunked"),
    "Content-Type",
    "Content-Length",
    (req) -> N,
    (res) -> N
)

protocol!(HTTP_1_1)
protocol!(("HTTP", "1.0"), HTTP_1_1)

HTTPRequest(io::AbstractIOSocket) = Request(io, HTTP_1_1)
HTTPResponse(io::AbstractIOSocket) = Response(io, HTTP_1_1)


text:
module Metropolis

using Random
using RandomNumbers.Xorshifts
using ProgressMeter
using StaticArrays
using LinearAlgebra: norm
using FileIO
using CellListMap
using Parameters
using InteractiveUtils
using Base.Threads

# Includes
include(joinpath("ensembles", "ensemble.jl"))
include("types.jl")
include("caches.jl")
include(joinpath("potentials", "potential.jl"))
include(joinpath("potentials", "lennard_jones.jl"))
include(joinpath("ensembles", "nvt.jl"))
include("utils.jl")
include("simulate.jl")
include("packing.jl")

# * Exports
# Types
export System, Simulation, Discrete, Continuous, LennardJones, NVT
# Interfaces
export potential_energy, forces, simulate!

end


text:
using GlyphSheets
using Base.Test

library = create_library()
typeface = create_typeface(library, "Butler_Regular.otf")

glyph = create_glyph(typeface, 'A', 12pt)
outlines = outline_decompose(glyph)

@test outlines isa Vector{GlyphSheets.DSL.PathOperation}
@test Any[GlyphSheets.DSL.Move((134, 173)), GlyphSheets.DSL.Line((227, 452)), GlyphSheets.DSL.Line((327, 173)), GlyphSheets.DSL.Line((134, 173)), GlyphSheets.DSL.Move((516, 9)), GlyphSheets.DSL.Cubic((495, 11), (478, 10), (453, 72)), GlyphSheets.DSL.Cubic((453, 72), (358, 328), (264, 576)), GlyphSheets.DSL.Line((256, 576)), GlyphSheets.DSL.Line((128, 192)), GlyphSheets.DSL.Cubic((61, 42), (42, 39), (13, 38)), GlyphSheets.DSL.Line((13, 28)), GlyphSheets.DSL.Line((184, 28)), GlyphSheets.DSL.Line((184, 38)), GlyphSheets.DSL.Cubic((77, 6), (111, 102), (132, 163)), GlyphSheets.DSL.Line((330, 163)), GlyphSheets.DSL.Line((357, 96)), GlyphSheets.DSL.Cubic((375, 53), (383, 9), (341, 9)), GlyphSheets.DSL.Line((341, 0)), GlyphSheets.DSL.Line((516, 0)), GlyphSheets.DSL.Line((516, 9))] == outlines

close_typeface(typeface)
close_library(library)


text:
module TestInner
using LinearAlgebra
using SARA
using Test
using GaussianDistributions: mean, std

@testset "inner loop" begin
    n = 8
    x = randn(n)
    sort!(x)
    σ = .05
    f(x) = cos(π*x)
    y = f.(x)
    @. y += σ*randn()

    domain = (-2, 2)
    next, posteriors, uncertainty = inner(x, [y], domain, σ)

    for i in eachindex(posteriors)
        p = posteriors[i]
        p_mean, p_std = evaluate(p, x)
        @test isapprox(p_mean, y, atol = 3σ)
        @test all(<(3σ), p_std)
    end

    # using Plots
    # xs = range(-3, 3, length = 128)
    # plot(xs, f.(xs))
    # scatter!(x, y)
    # for i in eachindex(posteriors)
    #     p = posteriors[i]
    #     μ_xs, σ_xs = evaluate(p, xs, normalization[i])
    #     plot!(xs, μ_xs, ribbon = σ_xs)
    # end
    # pnext = evaluate(posteriors[1], next, normalization[1])
    # scatter!(next, pnext) # label = "next point")
    # gui()
end

end


text:
# TODO: find/replace all [PkgName] with CamelCase, all [pkgname] with lowercase

# [WEBSITE]

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

function _create_plot(pkg::[PkgName]Package; kw...)
  d = Dict(kw)
  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc
  Plot(nothing, pkg, 0, d, Dict[])
end

function _add_series(::[PkgName]Package, plt::Plot; kw...)
  d = Dict(kw)
  # TODO: add one series to the underlying package
  push!(plt.seriesargs, d)
  plt
end

function _add_annotations{X,Y,V}(plt::Plot{[PkgName]Package}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    # TODO: add the annotation to the plot
  end
end

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

function _before_update_plot(plt::Plot{[PkgName]Package})
end

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{[PkgName]Package}, d::Dict)
end

function _update_plot_pos_size(plt::PlottingObject{[PkgName]Package}, d::Dict)
end

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

# accessors for x/y data

function Base.getindex(plt::Plot{[PkgName]Package}, i::Int)
  series = plt.o.lines[i]
  series.x, series.y
end

function Base.setindex!(plt::Plot{[PkgName]Package}, xy::Tuple, i::Integer)
  series = plt.o.lines[i]
  series.x, series.y = xy
  plt
end

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

function _create_subplot(subplt::Subplot{[PkgName]Package}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
end

function _expand_limits(lims, plt::Plot{[PkgName]Package}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{[PkgName]Package}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

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

function Base.writemime(io::IO, ::MIME"image/png", plt::PlottingObject{[PkgName]Package})
  # TODO: write a png to io
end

function Base.display(::PlotsDisplay, plt::Plot{[PkgName]Package})
  # TODO: display/show the plot
end

function Base.display(::PlotsDisplay, plt::Subplot{[PkgName]Package})
  # TODO: display/show the subplot
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

const TEMP_REF = 20.0
const PRESSURE_REF = 1.0
const TEMP_REF_UNITFUL = TEMP_REF * u"°C"
export TEMP_REF, PRESSURE_REF

const DISPFORM_NAMES = ["Schott", "Sellmeier1", "Herzberger", "Sellmeier2", "Conrady", "Sellmeier3", "HandbookOfOptics1", "HandbookOfOptics2", "Sellmeier4", "Extended1", "Sellmeier5", "Extended2", "Extended3"]
const STATUS = ["Standard", "Preferred", "Obsolete", "Special", "Melt"]


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "FLTK"
version = v"1.3.8"

# Collection of sources required to build FLTK
sources = [
    ArchiveSource("https://www.fltk.org/pub/fltk/$(version)/fltk-$(version)-source.tar.gz",
                  "f3c1102b07eb0e7a50538f9fc9037c18387165bc70d4b626e94ab725b9d4d1bf")
]

# Bash recipe for building across all platforms
script = raw"""
apk add fltk fltk-dev fltk-fluid # Host FLUID install required for cross-compilation
if [[ ${target} == x86_64-linux-musl* ]]; then
    # Remove host system libraries otherwise get "undefined reference to getrandom" error message on x86_64
    rm /usr/lib/libexpat*
    rm /lib/libuuid*
    rm /usr/lib/libfontconfig*
fi
cd ${WORKSPACE}/srcdir/fltk-*
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DOPTION_BUILD_SHARED_LIBS=ON \
    -DOPTION_USE_SYSTEM_LIBJPEG=ON \
    -DOPTION_USE_SYSTEM_LIBPNG=ON \
    -DOPTION_USE_SYSTEM_ZLIB=ON \
    -DOPTION_USE_THREADS=ON
make -j${nproc}
make install
install_license ../COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libfltk", :libfltk),
    LibraryProduct("libfltk_forms", :libfltk_forms),
    LibraryProduct("libfltk_gl", :libfltk_gl),
    LibraryProduct("libfltk_images", :libfltk_images)
]

# Some dependencies are needed only on Linux and FreeBSD
x11_platforms = filter(p ->Sys.islinux(p) || Sys.isfreebsd(p), platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"; platforms=x11_platforms),
    Dependency("Fontconfig_jll"),
    Dependency("FreeType2_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("Libglvnd_jll"; platforms=x11_platforms),
    Dependency("libpng_jll"),
    Dependency("Xorg_libX11_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXext_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXfixes_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXft_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXinerama_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXrender_jll"; platforms=x11_platforms),
    Dependency("Zlib_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"9")


text:
# parameters
mass_bin_name = ARGS[1] # "1540_1560"
tslice = "t1"
path_wavelist = "src"
path_to_working_folder = "data"
const M3pi = Meta.parse(mass_bin_name[1:4])/1000
const Nwave_extended = 88
#########################################################

@show ARGS
push!(LOAD_PATH,"src")
using SDMHelper
using PWAHelper
using amplitudes_compass
using DelimitedFiles
using Plots
using LinearAlgebra
using PlotHelper

# MAIN FIT
##### build the model:

# set names
for app in ["rd", "mc", "fu"]
    pwf = path_to_working_folder
    @eval $(Symbol("basisfunc_"*app)) = joinpath($pwf,"variables_$(mass_bin_name)_$(tslice)_"*$app*".txt")
end
# calculate normfact by getting number of lines
nlines_rd = Meta.parse(split(read(`wc -l $(basisfunc_rd)`,String)," ")[1])
nlines_mc = Meta.parse(split(read(`wc -l $(basisfunc_mc)`,String)," ")[1])
nlines_fu = Meta.parse(split(read(`wc -l $(basisfunc_fu)`,String)," ")[1])
normfact = nlines_fu/nlines_mc*nlines_rd
# read matrix of integrals
BmatFU = read_cmatrix("data/integrmat_$(mass_bin_name)_$(tslice)_fu.txt");
size(readdlm(joinpath(path_wavelist,"wavelist_formated.txt")),1)
# Model description
wavelist = get_wavelist(joinpath(path_wavelist,"wavelist_formated.txt");
     path_to_thresholds=joinpath(path_wavelist,"thresholds_formated.txt"),
     M3pi=M3pi)
const noϵ =  [1] # flat wave
const posϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="+"]
const negϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="-"]
const ModelBlocks = [noϵ, posϵ, negϵ, negϵ]
const Nwaves = size(wavelist,1)

println("Start calculating SDMs")
SDM = []
list_of_files = split(readline("data/llhfit/$(mass_bin_name)/llh_attmpts_$(mass_bin_name)_$(tslice).txt"),"\t")[1]#ARGS[3:end]
path_and_filename = list_of_files
minpars = readdlm(path_and_filename);
SDM = normfact*pars_to_SDM(minpars, BmatFU, ModelBlocks)

threshold_mask = get_threshold_mask(joinpath(path_wavelist,"thresholds_formated.txt"), M3pi, 88)
SDM_enlarged=enlarge_with_zeros!(SDM,threshold_mask)

# PUBLISHED DATA

SDM_mass_bin_name =string(trunc(Int,100- ((2480 - Meta.parse(split(mass_bin_name,"_")[1])) ./ 20)))                #ARGS[2]
#Compass SDM
path_to_SDM = "/mnt/data/compass/2008/pwa_results/SDMs/0.100000-0.112853/";
SDM_RD = read_compass_SDM(joinpath(path_to_SDM,"sdm$(SDM_mass_bin_name)."),
    path_to_wavelist   = joinpath("src", "wavelist_formated.txt");
    M3pi = M3pi)

SDM_RD_err = read_compass_SDM(joinpath(path_to_SDM,"sdm$(SDM_mass_bin_name)-err."),
    path_to_wavelist   = joinpath("src", "wavelist_formated.txt");
    M3pi = M3pi)

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

# SDM MAIN FIT PLOTTING

bar(real.(diag(SDM_enlarged)), lab="S.S",color = "red" ,size=(800,500))
scatter!(real.(diag(SDM_RD)),yerror=real.(diag(SDM_RD_err)), lab="F.H.-D.R.",color = "green",
                        xlab="# wave",ylab ="Magnitude", title = "Diagonal of the SDM_$(mass_bin_name)_$(tslice)", size=(800,500))
savefig(joinpath("plots/sdm_results/final","sdm_results_$(mass_bin_name)_$(tslice)_test.png"))

# BOOTSTRAP PLOTTING

bootstrap_file = "data/llhfit_bootstrap_$(mass_bin_name)_$(tslice).txt"
BootstrapResults = readdlm(bootstrap_file)
Nbstrap_attempts =vcat(size(BootstrapResults[1,:])...)[1]
SDMs = [Array{Complex{Float64}}(undef, Nwave_extended,Nwave_extended) for idx in 1:Nbstrap_attempts]
@time for b in 1:Nbstrap_attempts
    _br = BootstrapResults[2:end,b]
    SDMs[b] .= enlarge_with_zeros!(normfact*pars_to_SDM(_br, BmatFU, ModelBlocks),threshold_mask)
end

plotBSTsummary(SDMs, SDM_enlarged, SDM_RD, SDM_RD_err) #  tosort=true, toannotate=true
savefig(joinpath("plots/bootstrap_combined","bootstrap_combined_$(mass_bin_name)_$(tslice)_test.pdf"))
#saveBSTSplot(joinpath("plots/bootstrap_summary","bootstrap_summary_$(mass_bin_name)_$(tslice)_test.pdf"),
            # SDMs, SDM, SDM_RD, SDM_RD_err)


text:
function print_indent(n)
  i = 0
  while i < n
    print("  ")
    i += 1
  end
end

function json_print_node(val, lv, pretty)
  if pretty
    print_indent(lv)
  end

  if val isa String || val isa SubString
    print('"' * val * '"')
  elseif val isa Int
    print(val)
  elseif val isa Vector
    json_print_list(val, lv, pretty)
  else
    error("type: $(typeof(val))")
  end
end

function json_print_list(list, lv, pretty)
  print("[")
  if pretty
    print("\n")
  end

  i = 0
  for x in list
    json_print_node(x, lv + 1, pretty)

    if i < length(list) - 1
      if pretty
        print(",")
      else
        print(", ")
      end
    end
    if pretty
      print("\n")
    end

    i += 1
  end

  if pretty
    print_indent(lv)
  end
  print("]")
end

function json_print(list, pretty)
  json_print_list(list, 0, pretty)
end

function json_parse_node(rest)
  if 0 < match_int(rest)
    # byte size
    bsize = match_int(rest)
    s = rest[1:bsize]
    n = parse(Int, s)
    [n, bsize]
  elseif 0 <= match_str(rest)
    bsize = match_str(rest)
    s = rest[2:prevind(rest, bsize + 2)]
    [s, bsize + 2]
  elseif rest[1] == '['
    json_parse_list(rest)
  else
    error("unexpected pattern")
  end
end

function json_parse_list(rest)
  list = Any[]
  li = 1

  bi = 2 # byte index
  while bi <= lastindex(rest)
    _rest = rest[bi:end]

    if _rest[1] == ']'
      bi += 1
      break
    elseif _rest[1] == '\n'
      bi += 1
    elseif _rest[1] == ' '
      bi += 1
    elseif _rest[1] == ','
      bi += 1
    else
      val_bi = json_parse_node(_rest)
      push!(list, val_bi[1])
      bi += val_bi[2]
    end
  end

  [list, bi]
end

function json_parse(json)
  list_bi = json_parse_list(json)
  list_bi[1]
end


text:
type TridiagonalOperator
  diagonal::Vector{Float64}
  lowerDiagonal::Vector{Float64}
  upperDiagonal::Vector{Float64}
  temp::Vector{Float64}
  n::Int
end

TridiagonalOperator(n::Int) = TridiagonalOperator(zeros(n), zeros(n - 1), zeros(n - 1), zeros(n), n)
TridiagonalOperator() = TridiagonalOperator(Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), 0)
TridiagIdentity(n::Int) = TridiagonalOperator(ones(n), zeros(n - 1), zeros(n - 1), zeros(n), n)

function set_first_row!(L::TridiagonalOperator, valB::Float64, valC::Float64)
  L.diagonal[1] = valB
  L.upperDiagonal[1] = valC

  return L
end

function set_mid_row!(L::TridiagonalOperator, i::Int, valA::Float64, valB::Float64, valC::Float64)
  i >= 2 && i <= L.n - 1 || error("out of range for Tridiagonal System: set mid row")

  L.lowerDiagonal[i - 1] = valA
  L.diagonal[i] = valB
  L.upperDiagonal[i] = valC

  return L
end

function set_last_row!(L::TridiagonalOperator, valA::Float64, valB::Float64)
  L.lowerDiagonal[L.n - 1] = valA
  L.diagonal[L.n] = valB

  return L
end

function solve_for(L::TridiagonalOperator, rhs::Vector{Float64})
  # create tridiagonal
  tri = Tridiagonal(L.lowerDiagonal, L.diagonal, L.upperDiagonal)

  # lu factorization
  Alu = lufact(tri)

  return Alu \ rhs
end

function solve_for!(L::TridiagonalOperator, rhs::Vector{Float64}, result::Vector{Float64})
  bet = L.diagonal[1]
  result[1] = rhs[1] / bet

  for j = 2:L.n
    L.temp[j] = L.upperDiagonal[j - 1] / bet
    bet = L.diagonal[j] - L.lowerDiagonal[j - 1] * L.temp[j]
    result[j] = (rhs[j] - L.lowerDiagonal[j - 1] * result[j - 1]) / bet
  end

  for j = L.n-1:-1:2
    result[j] -= L.temp[j + 1] * result[j + 1]
  end

  result[1] -= L.temp[2] * result[2]

  return L, result
end

function apply_to(L::TridiagonalOperator, v::Vector{Float64})
  L.n == length(v) || error("v of incorrect length")

  res = L.diagonal .* v

  res[1] = L.upperDiagonal[1] * v[2]
  for j = 2:L.n - 1
    res[j] += L.lowerDiagonal[j-1] * v[j-1] + L.upperDiagonal[j] * v[j+1]
  end

  res[end] += L.lowerDiagonal[end-1] * v[end-1]

  return res
end

# basic operators
import Base.+, Base.-, Base.*, Base./, Base.copy

function -(D::TridiagonalOperator)
  low = -D.lowerDiagonal
  mid = -D.diagonal
  high = -D.upperDiagonal

  return TridiagonalOperator(mid, low, high, D.temp, D.n)
end

function +(D1::TridiagonalOperator, D2::TridiagonalOperator)
  D1.n == D2.n || error("Dimension mismatch")
  low = D1.lowerDiagonal + D2.lowerDiagonal
  mid = D1.diagonal + D2.diagonal
  high = D1.upperDiagonal + D2.upperDiagonal

  return TridiagonalOperator(mid, low, high, D1.temp, D1.n)
end

function -(D1::TridiagonalOperator, D2::TridiagonalOperator)
  D1.n == D2.n || error("Dimension mismatch")
  low = D1.lowerDiagonal - D2.lowerDiagonal
  mid = D1.diagonal - D2.diagonal
  high = D1.upperDiagonal - D2.upperDiagonal

  return TridiagonalOperator(mid, low, high, D1.temp, D1.n)
end

function *(a::Number, D::TridiagonalOperator)
  low = D.lowerDiagonal * a
  mid = D.diagonal * a
  high = D.upperDiagonal * a

  return TridiagonalOperator(mid, low, high, D.temp, D.n)
end

function *(D::TridiagonalOperator, a::Number)
  low = D.lowerDiagonal * a
  mid = D.diagonal * a
  high = D.upperDiagonal * a

  return TridiagonalOperator(mid, low, high, D.temp, D.n)
end

function /(D::TridiagonalOperator, a::Number)
  low = D.lowerDiagonal / a
  mid = D.diagonal / a
  high = D.upperDiagonal / a

  return TridiagonalOperator(mid, low, high, D.temp, D.n)
end

copy(L::TridiagonalOperator) = TridiagonalOperator(L.diagonal, L.lowerDiagonal, L.upperDiagonal, L.temp, L.n)


text:
###########################################
# Unconstrained optimization
###########################################

f(x) = sin(x[1] + x[2]) + cos(x[1])^2
g(x) = [cos(x[1] + x[2]) - 2*cos(x[1])*sin(x[1]); cos(x[1] + x[2])]

f(x1,x2) = f([x1;x2])

f([0; 0])
f(0, 0)

#### Exercise

using Plots

xs = range(-3, 1, length = 40)
ys = range(-2, 1, length = 40)

contourf(xs, ys, f, color = :jet)

#### Exercise

finite_difference(f, x::Real; h=1e-8) = (f(x+h) - f(x)) / h

#### Exercise

x = [-2; -1]
fin_diff(h) = finite_difference(y -> f(x[1], y), x[2]; h=h)

true_grad = g(x)[2]

hs = 10. .^ (-15:0.01:-1)

plot(hs, fin_diff,
    xlabel = "h",
    ylabel = "Partial gradient wrt y",
    label = ["Approximation" "True gradient"],
    xscale = :log10,
)

hline!([true_grad]; label =  "True gradient")

#### Numerical errors

x = 1
h = 1e-13

(x+h)^2 - x^2
2*x*h + h^2

#### Exercise

x = [-2; -1]
α = 0.25
x_grad = [x x.+α.*g(x)]

contourf(xs, ys, f; color = :jet)
plot!(x_grad[1, :], x_grad[2, :];
    line = (:arrow, 4, :black),
    label = "",
)

#### Exercise

function optim(f, g, x, α; max_iter=100)
    xs = zeros(length(x), max_iter+1)
    xs[:,1] = x
    for i in 1:max_iter
        x -= α*g(x)
        xs[:,i+1] = x
    end
    return xs
end

#### Animation

using Random

function create_anim(
    f,
    path,
    xlims,
    ylims,
    file_name = joinpath(pwd(), randstring(12) * ".gif");
    xbounds = xlims,
    ybounds = ylims,
    fps = 15,
)
    xs = range(xlims...; length = 100)
    ys = range(ylims...; length = 100)
    plt = contourf(xs, ys, f; color = :jet)

    # add constraints if provided
    if !(xbounds == xlims && ybounds == ylims)
        x_rect = [xbounds[1]; xbounds[2]; xbounds[2]; xbounds[1]; xbounds[1]]
        y_rect = [ybounds[1]; ybounds[1]; ybounds[2]; ybounds[2]; ybounds[1]]

        plot!(x_rect, y_rect; line = (2, :dash, :red), label="")
    end

    # add an empty plot
    plot!(Float64[], Float64[]; line = (4, :arrow, :black), label = "")

    # extract the last plot series
    plt_path = plt.series_list[end]

    # create the animation and save it
    anim = Animation()
    for x in eachcol(path)
        push!(plt_path, x[1], x[2]) # add a new point
        frame(anim)
    end
    gif(anim, file_name; fps = fps, show_msg = false)
    return nothing
end

#### Exercise

x_gd = optim([], g, [0; -1], 0.1)

xlims = (-3, 1)
ylims = (-2, 1)
create_anim(f, x_gd, xlims, ylims, "anim1.gif")

f_gd = [f(x) for x in eachcol(x_gd)]

plot(f_gd, label="", xlabel="Iteration", ylabel="Function value")

#### Different stepsizes

x_gd = optim([], g, [0; -1], 0.01)

create_anim(f, x_gd, xlims, ylims, "anim2.gif")

x_gd = optim([], g, [0; -1], 1)

create_anim(f, x_gd, xlims, ylims, "anim3.gif")

#### Stepsize selection

abstract type Step end

struct GD <: Step
    α::Float64
end

#### Gradient descent

optim_step(s::GD, f, g, x) = -s.α*g(x)

function optim(f, g, x, s::Step; max_iter=100)
    xs = zeros(length(x), max_iter+1)
    xs[:,1] = x
    for i in 1:max_iter
        x += optim_step(s, f, g, x)
        xs[:,i+1] = x
    end
    return xs
end

gd = GD(0.1)
x_opt = optim(f, g, [0;-1], gd)

create_anim(f, x_opt, xlims, ylims, "anim4.gif")

#### Exercise

struct Armijo <: Step
    c::Float64
    α_max::Float64
end

function optim_step(s::Armijo, f, g, x)
    fun = f(x)
    grad = g(x)
    α = s.α_max
    while f(x .- α*grad) > fun - s.c*α*(grad'*grad)
        α /= 2
        if α <= 1e-6
            warning("Armijo line search failed.")
            break
        end
    end
    return -α*grad
end

gd = Armijo(1e-4, 1)
x_opt = optim(f, g, [0;-1], gd)

create_anim(f, x_opt, xlims, ylims, "anim5.gif")

###########################################
# Constrained optimization
###########################################

f(x) = sin(x[1] + x[2]) + cos(x[1])^2
g(x) = [cos(x[1] + x[2]) - 2*cos(x[1])*sin(x[1]); cos(x[1] + x[2])]

f(x1,x2) = f([x1;x2])

#### Projected gradients

function optim(f, g, P, x, α; max_iter=100)
    xs = zeros(length(x), max_iter+1)
    ys = zeros(length(x), max_iter)
    xs[:,1] = x
    for i in 1:max_iter
        ys[:,i] = xs[:,i] - α*g(xs[:,i])
        xs[:,i+1] = P(ys[:,i])
    end
    return xs, ys
end

P(x, x_min, x_max) = min.(max.(x, x_min), x_max)

x_min = [-1; -1]
x_max = [0; 0]

xs, ys = optim(f, g, x -> P(x,x_min,x_max), [0;-1], 0.1)

#### Plot 1

xlims = (-3, 1)
ylims = (-2, 1)

create_anim(f, xs, xlims, ylims, "anim6.gif";
    xbounds=(x_min[1], x_max[1]),
    ybounds=(x_min[2], x_max[2]),
)

#### Plot 2

xys = hcat(reshape([xs[:,1:end-1]; ys][:], 2, :), xs[:,end])

create_anim(f, xys, xlims, ylims, "anim7.gif";
    xbounds=(x_min[1], x_max[1]),
    ybounds=(x_min[2], x_max[2]),
)

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 2a0680b2-a0f5-11eb-2658-d794c2f5d8a6
begin
	using LightGraphs
	using LinearAlgebra
	using Plots
	using PlutoUI
	using GraphRecipes
	using Statistics
	using LsqFit
	using Random
end

# ╔═╡ 2ea044d4-0ef4-4777-be1d-120b6e33b853
html"""
Term Paper 4

This work is a submission for the course Modelling Complex Systems(IDC-631) for the year 2020-21 Spring (Even) Semester by Dhruva Sambrani, MS18163
"""

# ╔═╡ 521ee215-35b3-432f-b584-1cc98daa4c5d
md"""
# Abstract

This paper looks at the robustness of the mean path length of many network topologies to removal of edges. A definition of "fitness" of a graph is also provided which attempts to capture real-world constraints on the networks we are considering.

## Structure of the paper

- **§1** Graphs and Networks - A brief introduction
- **§2** The Storm Story - The problem statement
- **§3** Gyms and Strength - A definition of network robustness
- **§4** Results and Conclusions
- **§5** References

!!! note "Note to reader"
    The author believes that papers which involve a lot of simulation and code should also be accompanied with the code they use so that readers can gain deeper insight into the processes described in the paper. Hence, the author has attempted to create a paper which interlaces the information with the corresponding code. However, in the interest of keeping the code aligned to the topic at hand, certain auxillary functions have been hidden. The content of these hidden cells can be found in the executable version of this notebook [here]()

The coding language used is `Julia`[^julia] and following document is made with `Pluto.jl` [^pluto].
"""

# ╔═╡ 5d01718d-aa67-467d-9a9c-1f65ac6d6dca
TableOfContents(depth=6)

# ╔═╡ 6b4e7cd4-d8ab-4d1b-a0df-a299595e4739
md"
# §1 Graphs and Networks - A brief introduction

Most things in our world do not exist in isolation, and are related to each other in some way or another. Take all people. They can be connected by, well, relationships. Houses are connected over a powergrid, and computers over a network. These networks are modelled mathematically by a structure known as Graphs.

Graphs are completely definied by a set of vertices(nodes) $\{v_i\}$ and edges $\{e_{i\rightarrow j}\}$. Often, the vertices are considered to be individuals in the system and the edges are considered to be the existance of a relationship between the edges.
"

# ╔═╡ 7fd589b6-e188-4f67-a6eb-e04d4c886c51
md"""## Example Graphs
$(plot(
	graphplot(barbell_graph(5, 5), title="Barbell"),
	graphplot(complete_graph(10), title="Complete"),
	graphplot(watts_strogatz(10, 4, 0.3), title="Watts-Strogtaz"),
	graphplot(erdos_renyi(10, 0.4), title="Erdos-Renyi"),
))"""

# ╔═╡ 1d8f1d47-d0f7-4901-b91f-4308f79681f5
md"""## Graph Properties

Graphs, like all other mathematical objects have certain properties associated with them.
"""

# ╔═╡ 63aed5e9-f17a-4468-839b-e5f01d39e334
md"""
### Adjacency Matrix

An adjacency matrix is a bit matrix which is defined as $A_{ij} = e_{i\rightarrow j}$. For undirected, unweighted graphs, this is a symmetric matrix. The adjacency matrix has all the information about the topology of the graph.

For example, the adjacency matrix of a Barbell graph looks like this-
"""

# ╔═╡ 5e5460e6-180f-4474-895a-8eee1e7c181a
adjacency_matrix(barbell_graph(5,5))

# ╔═╡ 92b62f06-c57d-43f5-87d7-a5ae409e1402
md"### Connectedness

A graph is connected if there exists a path between any two nodes on the graph.

As an example, an Erdos-Renyi graph may not be connected
"

# ╔═╡ 8aec44e0-3aa9-4228-87a4-d605d23e2ded
g = erdos_renyi(10,20)

# ╔═╡ 48a7f6cd-a457-4c3a-b47a-2ccc612b0878
graphplot(g, size=(300,225))

# ╔═╡ 069d1c41-f378-4364-ad1e-cd80c3dcd8e4
is_connected(g) ? md"**Graph is connected**" : md"**Graph is not connected**"

# ╔═╡ ff945901-8c13-4404-b71c-33616a478aef
md"""# §2 The Storm Story - The problem statement

Consider this. It is a stormy night when your workplace network stops working and you are cutoff from the world. You are only scrolling Social Media, and you imagine remember the good times when you could "like" your friends' pictures. Living in the 21st century, you'd expect the network to have multiple paths to the internet.

This paper asks exactly that. How many connections can be dropped from the network before some devices are completely disconnected from the rest of the world. While you could wait for another 10 minutes before your site loading, there are some situations which demand that the network stay up for as long as possible before complete failure. These include banking and security systems, emergency response networks, remote control systems and food ordering systems. Ok, the last one was a stretch.

Similar models can also be used to model divergence of species in a population where some geo-physical/biological barrier prevents genes from spreading between the sub-populations.

Along the same lines, power grids too can be tested for robustness. However, powergrids add an extra layer of complexity, as
1. They are directed graphs
2. Source and sink nodes are not equivalent
Due to these complexities, we do not look at this particular problem and leave it for future work.

However, we can add certain further complexity to the problem. While the network may remain connected, the mean path length can increase to such high values, that it becomes _effectively_ disconnected. This is especially important in the case where there are competing processes at play.
- In the case of genes, it is the rate of mutations and the rate of change of the environment.
- In the case of emergency response, it is the priority of the case
- In remote control, too high ping rates can lead to poor control
- In the case of internet, the competing process is often the patience of the user.

Hence, we need a way to measure the **robustness** of the network to connection failures.
"""

# ╔═╡ 89fdc6e6-a16a-4c53-b64f-409798d1be7f
md"""
# §3 Gyms and Strength - A definition of network robustness

## The Graph theory model

Graph theory is an elegant way to translate the problem into a mathematical and computational struture that can be manipulated.

In this scenario, the nodes are the devices, organisms, or some equivalent entity. The edges are the "connectivity" between two nodes. In the case of an internet network, it is the channel for information to pass through."""

# ╔═╡ 926e2dcd-5ee7-4507-b41d-879b58f42662
md"""## Network failure

As described in the previous section, the network "fails" when it becomes disconnected or mean path length increases beyond a particular cutoff.

!!! note "Network Failure"
	- A network **fails** when,
	  - It becomes disconnected OR
	  - Its mean path length increases beyond some cutoff, as a multiple of initial mean path length or an absolute value."""

# ╔═╡ ffc21b99-e466-424c-918a-d2cd30844278
md"## Topology dependance

Naively, it seems like that the network which has more edges will fare better than one which has few. But this is not the case. Consider the Barbell graph.
"

# ╔═╡ 5a5e98c6-d983-4fbb-b073-c347968a6415
graphplot(barbell_graph(5,5), size=(300,225))

# ╔═╡ e4c5c26b-4976-44af-8aa9-2d7d22f84af0
md"In this graph, if the edge which connects the two subgraphs(cliques) is broken, the graph immediatly fails"

# ╔═╡ 01455f58-0492-41a6-ae44-3cd39750aaa1
let
	g = barbell_graph(5,5)
	rem_edge!(g, 5, 6)
	isc = is_connected(g) ? "connected" : "not connected"
	md"$(graphplot(g, size=(300,225))) Graph is $isc"
end

# ╔═╡ b19b1856-c46b-4812-bc3e-9e7dce4b5b69
md"Hence, the robustness of the graph is most definitely dependant on the topology of the network."

# ╔═╡ 8fe47f6d-1a41-4e58-b7dd-eb365065404f
md"""## Robustness

The simulation then could involve creating a graph, checking connectedness, breaking edges one by one and checking connectedness after each edge removal. Then, a **measure of robustness is the amount of edges we can remove before the network fails**. However, it would be unfair to a graph of 100 edges which fails on removing 50 edges to be considered worse than a network of 1000 edges that fails on removing 100 edges. So a better measure of robustness is the **_fraction_ of edges which can fail before the network fails**. As we saw in the previous subsection, robustness clearly depends _which_ edge we remove. Hence, we again amend our definition of robustness to read as **_the average_ fraction of edges which can fail before the network fails**

!!! note "Definition of Robustness"
	Robustness of a network is defined as the the _average fraction_ of edges which can fail before the network _fails_.
"""

# ╔═╡ 07953ab8-6449-4968-8d3d-75548706170e
md"""### Small scale Robustness

Sometimes, it may be more useful to measure how the mean path length increases if few edges are removed, especially when edge failure is a rare event, or faults can be fixed quickly, and when responsiveness to the network is of higher importance than long scale robustness. For this, we define small scale robustness of a network.

!!! note "Small scale robustness"
	It is defined as the slope of the mean path length vs the fraction of nodes removed for small fractions.

There are some evident disadvantages to such a definition
- Definition of "small" - this can be solved by choosing a range of interest
- Difficult to define for small networks - Since the definition involves fitting, small networks suffer having too few data points which can not be increased.
- Computationally expensive - Compared to our definition of robustness, SSR is more computationally expensive.
- Linear nature of plot is not guaranteed - Consider the Barbell network.
- Linear nature may not remain over the range of interest - More generalized fits may be more useful
"""

# ╔═╡ 0a470cac-b6a0-45bf-9b67-b1ad22458011
md"# §4 Results and Conclusion

In the following section, we implement the setup above for multiple topologies, and look at which topologies are more robust according to our definition.
"

# ╔═╡ 611f678e-48a7-4a0b-8f8d-9abfea658eba
md"## Implementation"

# ╔═╡ 894715a6-9ca5-49fb-be98-ad9131efbbfc
begin
	# Remove a random edge
	breakrandom!(g::AbstractGraph) = rem_edge!(g, rand(collect(edges(g))))
	mean_path_length(g) =
		sum(floyd_warshall_shortest_paths(g).dists)/(nv(g)^2-nv(g))
end

# ╔═╡ 2902dc16-2a73-42bd-b2bb-7e43e3803c39
md"### Mean Path Length

This is the mean of the shortest path lengths between all nodes in the graph.

It is given by the formula

$$\mu(g) = \frac{\sum_{i\in\text{v}(g), j\in\text{v}(g),\ j\ne i} d(i,j)}{n(n-1)}$$

where $d(i,j)$ is the shortest possible path $i\rightarrow j$. If no such path exists, some definitions make it infinity, others make it 0. The appropriate definition should be chosen as per requirements.

The Floyd Warshall shortest paths algorithm provides a matrix of the shortest paths between all the nodes in a graph.

For example, the mean path length of a completely connected graph is $(mean_path_length(complete_graph(10))), as expected.
"

# ╔═╡ bef640f9-95c9-4e74-b7ad-91ac793f3a97
function ssr(pathlengths, frac)
	datarange = floor(Int, length(pathlengths)*frac)
	datarange > 5 || return -1 # network robustness is too low to measure SSR
	@. model(x, (m, c)) = m*x+c
	fit = curve_fit(
		model,
		collect(range(0, stop=frac, length=datarange)),
		pathlengths[1:datarange],
		[0.,0.]
	)
	return fit.param[1]
end

# ╔═╡ 8e9b6375-bfc4-4e55-9afd-a6e8c810b900
begin
	mpl_increase(mpls) = last(mpls)/first(mpls)
	function simulate!(g, cutoff_mpl_scale; seed=-1)
		if seed >= 0
			Random.seed!(seed)
		end
		pathlengths = [mean_path_length(g)]
		nedges = ne(g)
		while is_connected(g) || (mpl_increase(pathlengths) > cutoff_mpl_scale)
			breakrandom!(g)
			push!(pathlengths, mean_path_length(g))
		end
		pathlengths = length(pathlengths)>1 ? pathlengths[1:end-1] : pathlengths
		return (
			length(pathlengths)/nedges,
			pathlengths,
			length(pathlengths)>1 ? range(0, stop=length(pathlengths)/nedges, length=length(pathlengths)) : [0]
		)
	end
end

# ╔═╡ c71ed1c4-a13b-432c-93a3-c9d27cc6ce94
function makeplot((pls, x), name)
	plot(
		x,
		map(pls) do i;
			i/first(pls);
		end,
		xlims=(0,1),
		ylims=(1,5),
		xlabel="Fraction of edges removed",
		ylabel="Ratio of MPL",
		title=name
	)
end

# ╔═╡ c4e14fd7-d0de-4370-ad7f-06563cff1009
md"## Test Runs"

# ╔═╡ 22c6f858-d074-49ae-bf41-b3502c08ed15
seed = 5176

# ╔═╡ a1d6a627-fc7c-4e9c-9a64-14c6059dbb3f
begin
	p1 = makeplot(simulate!(watts_strogatz(50, 6, 0.1, seed=seed), Inf, seed=seed)[2:3], "Watts-Strogatz")
	p2 = makeplot(simulate!(complete_graph(50), Inf, seed=seed)[2:3], "Complete")
	p3 = makeplot(simulate!(barabasi_albert(50, 6, seed=seed), Inf, seed=seed)[2:3], "Barabasi-Albert")
	plot(p1, p2, p3, size=(600, 600))
end

# ╔═╡ fda41875-dced-4ea3-953a-b0691fc3aa63
function get_robustness_measures(generator, steps; cutoff_mpl_scale=Inf, frac=0.1)
	rbms = map(1:steps) do i
		g = generator()
		data = simulate!(g, cutoff_mpl_scale)
		println(length(data[2]))
		(data[1], ssr(data[2], frac))
	end
	return first.(rbms), filter(i->i>0, last.(rbms))
end

# ╔═╡ 57b518a2-be25-47e6-b556-f1707c35405b
begin
	function getstats(rbs, ssrs, name)
		ssrs = length(ssrs)==0 ? [-1] : ssrs
		md"""
	### $(name)

	| Property | Value for Robustness             |Value for SSR                      |
	|:--------:|:--------------------------------:|:---------------------------------:|
	| Median   |$(round(median(rbs), sigdigits=3))|$(round(median(ssrs), sigdigits=3))|
	| Mean     |$(round(mean(rbs), sigdigits=3))  |$(round(mean(ssrs), sigdigits=3))  |
	| StdDev   |$(round(std(rbs), sigdigits=3))   |$(round(std(ssrs), sigdigits=3))   |
	| Length   |$(length(rbs))                    |$(length(ssrs))                    |
	"""
	end
	getstats(name) = (dat) -> getstats(dat[1], dat[2], name)
end

# ╔═╡ 6370353a-fbc4-4f1b-be96-1942b486e617
md"## Runs over different topologies"

# ╔═╡ df107cca-e86f-48a8-9c7a-6c9c6036bc94
get_robustness_measures(()->watts_strogatz(50, 6, 0.01), 500) |> getstats("Watts-Strogatz low β(0.01)")

# ╔═╡ b1468ecd-c3e8-4009-9290-e944bb433ef5
get_robustness_measures(()->watts_strogatz(50, 6, 0.1), 500) |> getstats("Watts-Strogatz mid β(0.1)")

# ╔═╡ 8b0af9d6-e60b-483e-a848-0a1d8963e351
get_robustness_measures(()->watts_strogatz(50, 6, 0.5), 500) |> getstats("Watts-Strogatz high β(0.5)")

# ╔═╡ 576d7b10-8fbf-4aa8-893e-5dbf7b567812
get_robustness_measures(()->watts_strogatz(50, 20, 0.1), 500) |> getstats("Watts-Strogatz high k(20)")

# ╔═╡ 3ac2fc49-b17e-48d1-8991-d7f970331c31
get_robustness_measures(()->barabasi_albert(50, 6), 500) |> getstats("Barabasi-Albert mid k(6)")

# ╔═╡ 64b84265-b49d-4869-adcb-b9c68d505c4f
get_robustness_measures(()->barabasi_albert(50, 20), 500) |> getstats("Barabasi-Albert high k(20)")

# ╔═╡ 58be2740-f3ab-46af-9d64-b3897a7deeb3
get_robustness_measures(()->complete_graph(50), 100) |> getstats("Complete Graph")

# ╔═╡ 51a0b307-cb71-4280-af8f-8ad5afafa8c9
get_robustness_measures(()->barbell_graph(25,25), 100) |> getstats("Barbell Graph")

# ╔═╡ 4c13a4fe-f642-469b-8203-0d7796457581
get_robustness_measures(()->erdos_renyi(50, 100), 100) |> getstats("Erdos-Renyi Graph low k(10%)")

# ╔═╡ f35f13d3-03d3-4ed1-9246-4ae0f06835ac
get_robustness_measures(()->erdos_renyi(50, 625), 100) |> getstats("Erdos-Renyi Graph high k(50%)")

# ╔═╡ 67c579be-f039-4352-bb9b-1791ec1ec858
md"""## Conclusions

1. As a comparision between kinds of topology, a Complete graph has the highest robustness. This is not surprising, because after all, there are so many connections that random edges being removed don't affect the topology of the graph, and no particular edge increases leads to network failure.
2. For the Erdos-Renyi graph[^erdos_renyi], a very low k can lead to multiple nodes outside the main graph. For a high k results are similar to that of the Complete graph, which is not surprising, as a complete graph with random nodes removed is a random graph.
3. For the Watts-Strogatz[^watts_strogatz] graph, while increasing β decreases SSR, it also decreases robustness. If k is increased, it increases robustness and decreases SSR. This is in contrast to simply adding more nodes in the Erdos-Renyi, as it decreases SSR, which is preferable. This strongly shows that Watts-Strogatz model is resistant to small scale perturbations, leading to an evolutionary benefit, as predicted before by [^wiki].
4. The Barabasi-Albert[^barabasi_albert] graph is interesting in context of the other graphs, especially when compared to the Watts-Strogatz model. While Barabasi-Albert is more robust and has lower SSR for low k, for higher ks, the Barabasi-Albert has lower robustness and comparable SSR. This suggests that for low density edge graphs, a Barabasi-Albert topology is "better" than the Watts-Strogatz, but as density increases, the Watts-Strogatz model scales better.
5. The Barbell Graph reacts surprisingly well, but it has low robustness and high SSR.
"""

# ╔═╡ f7393bf5-ad86-4619-a69a-9ff7108b2870
md"# References
[^robustness]:[wiki: Small World Network#Network_robustness](https://en.wikipedia.org/wiki/Small-world_network#Network_robustness)
[^erdos_renyi]: [wiki: Erdos-Renyi_model](https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model)
[^watts_strogatz]:[Watts-Strogatz_model](https://en.wikipedia.org/wiki/Watts%E2%80%93Strogatz_model)
[^barabasi_albert]: [wiki: Barabasi-Albert_model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model)

## Packages used

[^lightgraphs]: JuliaGraphs/LightGraphs.jl: an optimized graphs package for the Julia programming language. Seth Bromberger, James Fairbanks, and other contributors. [10.5281/zenodo.889971](https://doi.org/10.5281/zenodo.889971)
[^graphrecepies]: GraphRecipes. Thomas Breloff, and other contributors. [github](https://github.com/JuliaPlots/GraphRecipes.jl)
[^julia]: Julia: A Fresh Approach to Numerical Computing. Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah. (2017) SIAM Review, 59: 65–98. doi: [10.1137/141000671](https://dx.doi.org/10.1137/141000671)
[^plots]: Plots - powerful convenience for visualization in Julia. Thomas Breloff. [docs.juliaplots.org](docs.juliaplots.org)
[^pluto]: Pluto - Reactive Notebooks for Julia. Fons van der Plas and Mikołaj Bochenski. [https://github.com/fonsp/Pluto.jl](https://github.com/fonsp/Pluto.jl)
[^plutoui]: PlutoUI. Fons van der Plas. [https://github.com/fonsp/PlutoUI.jl](https://github.com/fonsp/PlutoUI.jl)
[^stats]: Statistics - Julia stdlib for Statistics. [https://docs.julialang.org/en/v1/stdlib/Statistics/](https://docs.julialang.org/en/v1/stdlib/Statistics/)
[^lsqfit]: LsqFit.jl: Basic least-squares fitting in pure Julia under an MIT license. [https://github.com/JuliaNLSolvers/LsqFit.jl](https://github.com/JuliaNLSolvers/LsqFit.jl)
"

# ╔═╡ Cell order:
# ╟─2ea044d4-0ef4-4777-be1d-120b6e33b853
# ╟─521ee215-35b3-432f-b584-1cc98daa4c5d
# ╟─5d01718d-aa67-467d-9a9c-1f65ac6d6dca
# ╟─2a0680b2-a0f5-11eb-2658-d794c2f5d8a6
# ╟─6b4e7cd4-d8ab-4d1b-a0df-a299595e4739
# ╟─7fd589b6-e188-4f67-a6eb-e04d4c886c51
# ╟─1d8f1d47-d0f7-4901-b91f-4308f79681f5
# ╟─63aed5e9-f17a-4468-839b-e5f01d39e334
# ╟─5e5460e6-180f-4474-895a-8eee1e7c181a
# ╟─2902dc16-2a73-42bd-b2bb-7e43e3803c39
# ╟─92b62f06-c57d-43f5-87d7-a5ae409e1402
# ╟─8aec44e0-3aa9-4228-87a4-d605d23e2ded
# ╟─48a7f6cd-a457-4c3a-b47a-2ccc612b0878
# ╟─069d1c41-f378-4364-ad1e-cd80c3dcd8e4
# ╟─ff945901-8c13-4404-b71c-33616a478aef
# ╟─89fdc6e6-a16a-4c53-b64f-409798d1be7f
# ╟─926e2dcd-5ee7-4507-b41d-879b58f42662
# ╟─ffc21b99-e466-424c-918a-d2cd30844278
# ╟─5a5e98c6-d983-4fbb-b073-c347968a6415
# ╟─e4c5c26b-4976-44af-8aa9-2d7d22f84af0
# ╟─01455f58-0492-41a6-ae44-3cd39750aaa1
# ╟─b19b1856-c46b-4812-bc3e-9e7dce4b5b69
# ╟─8fe47f6d-1a41-4e58-b7dd-eb365065404f
# ╟─07953ab8-6449-4968-8d3d-75548706170e
# ╟─0a470cac-b6a0-45bf-9b67-b1ad22458011
# ╟─611f678e-48a7-4a0b-8f8d-9abfea658eba
# ╠═894715a6-9ca5-49fb-be98-ad9131efbbfc
# ╠═bef640f9-95c9-4e74-b7ad-91ac793f3a97
# ╠═8e9b6375-bfc4-4e55-9afd-a6e8c810b900
# ╠═c71ed1c4-a13b-432c-93a3-c9d27cc6ce94
# ╟─c4e14fd7-d0de-4370-ad7f-06563cff1009
# ╠═22c6f858-d074-49ae-bf41-b3502c08ed15
# ╠═a1d6a627-fc7c-4e9c-9a64-14c6059dbb3f
# ╠═fda41875-dced-4ea3-953a-b0691fc3aa63
# ╟─57b518a2-be25-47e6-b556-f1707c35405b
# ╟─6370353a-fbc4-4f1b-be96-1942b486e617
# ╟─df107cca-e86f-48a8-9c7a-6c9c6036bc94
# ╟─b1468ecd-c3e8-4009-9290-e944bb433ef5
# ╟─8b0af9d6-e60b-483e-a848-0a1d8963e351
# ╟─576d7b10-8fbf-4aa8-893e-5dbf7b567812
# ╟─3ac2fc49-b17e-48d1-8991-d7f970331c31
# ╟─64b84265-b49d-4869-adcb-b9c68d505c4f
# ╟─58be2740-f3ab-46af-9d64-b3897a7deeb3
# ╟─51a0b307-cb71-4280-af8f-8ad5afafa8c9
# ╟─4c13a4fe-f642-469b-8203-0d7796457581
# ╟─f35f13d3-03d3-4ed1-9246-4ae0f06835ac
# ╟─67c579be-f039-4352-bb9b-1791ec1ec858
# ╟─f7393bf5-ad86-4619-a69a-9ff7108b2870


text:
# This file investigates mean square error (MSE) with respect to free variables in 1D interpolation

using FractalTools
using Makie

# Construct interpolation data
f = FractalTools.wen            # For irregular data
# f = FractalTools.sinusoid     # For regular data
line = [[0.], [1.]]
pts = getdata(f, line, 11)

# Construct test data
tpts = getdata(line, 101)
npts = length(tpts)

# Compute errors
fvals = map(pnt -> f(pnt[1]), tpts)
freevars = 0.001 : 0.01 : 0.999
mse = map(freevars) do freevar
    interp = interpolate(pts, Interp1D(freevar))
    ivals = map(pnt -> interp(pnt[1]), tpts)
    sum((fvals - ivals).^2) / npts
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="MSE", title="1D Interpolation MSE")
stem!(ax, freevars, mse, color=:black)
save(joinpath(@__DIR__, "interp1d_mse.png"), fig)
display(fig)


text:
mutable struct CurrentPlot
    nullableplot::Union{AbstractPlot, Nothing}
end
const CURRENT_PLOT = CurrentPlot(nothing)

isplotnull() = CURRENT_PLOT.nullableplot === nothing

"""
    current()
Returns the Plot object for the current plot
"""
function current()
    if isplotnull()
        error("No current plot/subplot")
    end
    CURRENT_PLOT.nullableplot
end
current(plot::AbstractPlot) = (CURRENT_PLOT.nullableplot = plot)

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

Base.string(plt::Plot) = "Plot{$(plt.backend) n=$(plt.n)}"
Base.print(io::IO, plt::Plot) = print(io, string(plt))
Base.show(io::IO, plt::Plot) = print(io, string(plt))

getplot(plt::Plot) = plt
getattr(plt::Plot, idx::Int = 1) = plt.attr
convertSeriesIndex(plt::Plot, n::Int) = n

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

"""
The main plot command. Use `plot` to create a new plot object, and `plot!` to add to an existing one:

```
    plot(args...; kw...)                  # creates a new plot window, and sets it to be the current
    plot!(args...; kw...)                 # adds to the `current`
    plot!(plotobj, args...; kw...)        # adds to the plot `plotobj`
```

There are lots of ways to pass in data, and lots of keyword arguments... just try it and it will likely work as expected.
When you pass in matrices, it splits by columns. To see the list of available attributes, use the `plotattr([attr])`
function, where `attr` is the symbol `:Series:`, `:Subplot:`, `:Plot` or `:Axis`. Pass any attribute to `plotattr`
as a String to look up its docstring; e.g. `plotattr("seriestype")`.
"""
function plot(args...; kw...)
    # this creates a new plot with args/kw and sets it to be the current plot
    plotattributes = KW(kw)
    preprocessArgs!(plotattributes)

    # create an empty Plot then process
    plt = Plot()
    # plt.user_attr = plotattributes
    _plot!(plt, plotattributes, args)
end

# build a new plot from existing plots
# note: we split into plt1 and plts_tail so we can dispatch correctly
function plot(plt1::Plot, plts_tail::Plot...; kw...)
    plotattributes = KW(kw)
    preprocessArgs!(plotattributes)

    # build our plot vector from the args
    n = length(plts_tail) + 1
    plts = Array{Plot}(undef, n)
    plts[1] = plt1
    for (i,plt) in enumerate(plts_tail)
        plts[i+1] = plt
    end

    # compute the layout
    layout = layout_args(plotattributes, n)[1]
    num_sp = sum([length(p.subplots) for p in plts])

    # create a new plot object, with subplot list/map made of existing subplots.
    # note: we create a new backend figure for this new plot object
    # note: all subplots and series "belong" to this new plot...
    plt = Plot()

    # TODO: build the user_attr dict by creating "Any matrices" for the args of each subplot

    # TODO: replace this with proper processing from a merged user_attr KW
    # update plot args, first with existing plots, then override with plotattributes
    for p in plts
        _update_plot_args(plt, copy(p.attr))
        plt.n += p.n
    end
    _update_plot_args(plt, plotattributes)

    # pass new plot to the backend
    plt.o = _create_backend_figure(plt)
    plt.init = true

    series_attr = KW()
    for (k,v) in plotattributes
        if haskey(_series_defaults, k)
            series_attr[k] = pop!(plotattributes,k)
        end
    end

    # create the layout
    plt.layout, plt.subplots, plt.spmap = build_layout(layout, num_sp, copy(plts))

    # do we need to link any axes together?
    link_axes!(plt.layout, plt[:link])

    # initialize the subplots
    cmdidx = 1
    for (idx, sp) in enumerate(plt.subplots)
        _initialize_subplot(plt, sp)
        serieslist = series_list(sp)
        if sp in sp.plt.inset_subplots
            push!(plt.inset_subplots, sp)
        end
        sp.plt = plt
        sp.attr[:subplot_index] = idx
        for series in serieslist
            merge!(series.plotattributes, series_attr)
            _add_defaults!(series.plotattributes, plt, sp, cmdidx)
            push!(plt.series_list, series)
            _series_added(plt, series)
            cmdidx += 1
        end
    end

    # first apply any args for the subplots
    for (idx,sp) in enumerate(plt.subplots)
        _update_subplot_args(plt, sp, plotattributes, idx, false)
    end

    # finish up
    current(plt)
    _do_plot_show(plt, get(plotattributes, :show, default(:show)))
    plt
end

# this adds to the current plot, or creates a new plot if none are current
function  plot!(args...; kw...)
    local plt
    try
        plt = current()
    catch
        return plot(args...; kw...)
    end
    plot!(current(), args...; kw...)
end

# this adds to a specific plot... most plot commands will flow through here
function plot!(plt::Plot, args...; kw...)
    plotattributes = KW(kw)
    preprocessArgs!(plotattributes)
    # merge!(plt.user_attr, plotattributes)
    _plot!(plt, plotattributes, args)
end

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

# this is the core plotting function.  recursively apply recipes to build
# a list of series KW dicts.
# note: at entry, we only have those preprocessed args which were passed in... no default values yet
function _plot!(plt::Plot, plotattributes::KW, args::Tuple)
    plotattributes[:plot_object] = plt

    if !isempty(args) && !isdefined(Main, :StatsPlots) &&
            first(split(string(typeof(args[1])), ".")) == "DataFrames"
        @warn("You're trying to plot a DataFrame, but this functionality is provided by StatsPlots")
    end

    # --------------------------------
    # "USER RECIPES"
    # --------------------------------

    kw_list = _process_userrecipes(plt, plotattributes, args)

    # @info(1)
    # map(DD, kw_list)

    # --------------------------------
    # "PLOT RECIPES"
    # --------------------------------

    # "plot recipe", which acts like a series type, and is processed before
    # the plot layout is created, which allows for setting layouts and other plot-wide attributes.
    # we get inputs which have been fully processed by "user recipes" and "type recipes",
    # so we can expect standard vectors, surfaces, etc.  No defaults have been set yet.
    still_to_process = kw_list
    kw_list = KW[]
    while !isempty(still_to_process)
        next_kw = popfirst!(still_to_process)
        _process_plotrecipe(plt, next_kw, kw_list, still_to_process)
    end

    # @info(2)
    # map(DD, kw_list)

    # --------------------------------
    # Plot/Subplot/Layout setup
    # --------------------------------
    _plot_setup(plt, plotattributes, kw_list)
    _subplot_setup(plt, plotattributes, kw_list)

    # !!! note: At this point, kw_list is fully decomposed into individual series... one KW per series.          !!!
    # !!!       The next step is to recursively apply series recipes until the backend supports that series type !!!

    # --------------------------------
    # "SERIES RECIPES"
    # --------------------------------

    # @info(3)
    # map(DD, kw_list)

    for kw in kw_list
        sp::Subplot = kw[:subplot]
        # idx = get_subplot_index(plt, sp)

        # # we update subplot args in case something like the color palatte is part of the recipe
        # _update_subplot_args(plt, sp, kw, idx, true)

        # set default values, select from attribute cycles, and generally set the final attributes
        _add_defaults!(kw, plt, sp, command_idx(kw_list,kw))

        # now we have a fully specified series, with colors chosen.   we must recursively handle
        # series recipes, which dispatch on seriestype.  If a backend does not natively support a seriestype,
        # we check for a recipe that will convert that series type into one made up of lower-level components.
        # For example, a histogram is just a bar plot with binned data, a bar plot is really a filled step plot,
        # and a step plot is really just a path.  So any backend that supports drawing a path will implicitly
        # be able to support step, bar, and histogram plots (and any recipes that use those components).
        _process_seriesrecipe(plt, kw)
    end

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

    current(plt)

    # do we want to force display?
    # if plt[:show]
    #     gui(plt)
    # end
    _do_plot_show(plt, plt[:show])

    plt
end

# we're getting ready to display/output.  prep for layout calcs, then update
# the plot object after
function prepare_output(plt::Plot)
    _before_layout_calcs(plt)

    w, h = plt.attr[:size]
    plt.layout.bbox = BoundingBox(0mm, 0mm, w*px, h*px)

    # One pass down and back up the tree to compute the minimum padding
    # of the children on the perimeter.  This is an backend callback.
    _update_min_padding!(plt.layout)
    for sp in plt.inset_subplots
        _update_min_padding!(sp)
    end

    # now another pass down, to update the bounding boxes
    update_child_bboxes!(plt.layout)

    # update those bounding boxes of inset subplots
    update_inset_bboxes!(plt)

    # the backend callback, to reposition subplots, etc
    _update_plot_object(plt)
end

function backend_object(plt::Plot)
    prepare_output(plt)
    plt.o
end

# --------------------------------------------------------------------
# plot to a Subplot

function plot(sp::Subplot, args...; kw...)
    plt = sp.plt
    plot(plt, args...; kw..., subplot = findfirst(isequal(sp), plt.subplots))
end
function plot!(sp::Subplot, args...; kw...)
    plt = sp.plt
    plot!(plt, args...; kw..., subplot = findfirst(isequal(sp), plt.subplots))
end

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


text:
immutable PoolingLayerState <: AbstractScatteredLayerState
    blobs::Vector{Mocha.Blob}
    layer::PoolingLayer
end

function Mocha.setup(
        backend::Mocha.Backend,
        layer::PoolingLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    blobs = Vector{ScatteredBlob}(length(inputs))
    for idblob in eachindex(inputs)
        innodes = inputs[idblob].nodes
        outnodes = DataStructures.SortedDict{Path,AbstractNode,
            Base.Order.ForwardOrdering}()
        for inpath in keys(innodes)
            innode = innodes[inpath]
            inranges = innode.ranges
            inpathkeys = [ pair.first for pair in collect(inranges) ]
            subscripts = vcat(
                [find(pathkey.==inpathkeys) for pathkey in layer.pathkeys]...)
            indata = innode.data
            outdata = layer.pooling(indata, subscripts)
            outranges = collect(inranges)
            for subscript in subscripts
                outranges[subscript] =
                    Pair(inranges[subscript].first,
                        StepRange(
                        inranges[subscript].second.start,
                        inranges[subscript].second.stop -
                        inranges[subscript].second.start + 1,
                        inranges[subscript].second.stop))
            end
            outranges = tuple(outranges...)
            outnode = Node(outdata, outranges)
            outpath = inpath
            outnodes[outpath] = outnode
        end
        blobs[idblob] = ScatteredBlob(outnodes)
    end
    PoolingLayerState(blobs, layer)
end


text:
using LinearAlgebra
using FermiCG
using Printf
using Test
using LinearMaps
using Arpack
using Random
using Profile

atoms = []
push!(atoms,Atom(1,"H",[0,0,0]))
push!(atoms,Atom(2,"H",[0,0,1]))
push!(atoms,Atom(3,"H",[1,0,2]))
push!(atoms,Atom(4,"H",[1,0,3]))
push!(atoms,Atom(5,"H",[2,0,4]))
push!(atoms,Atom(6,"H",[2,0,5]))
push!(atoms,Atom(7,"H",[3,0,6]))
push!(atoms,Atom(8,"H",[3,0,7]))
#push!(atoms,Atom(9,"H",[0,0,8]))
#push!(atoms,Atom(10,"H",[0,0,9]))
#push!(atoms,Atom(11,"H",[0,0,10]))
#push!(atoms,Atom(12,"H",[0,0,11]))
basis = "6-31g"
#basis = "sto-3g"

mol     = Molecule(0,1,atoms, basis)
mf = FermiCG.pyscf_do_scf(mol)
e_mf = mf.e_tot - mf.energy_nuc()
@printf(" Mean-field energy %12.8f", e_mf)
ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
@time e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,4,4)
# @printf(" FCI Energy: %12.8f\n", e_fci)

norbs = size(ints.h1)[1]

problem = StringCI.FCIProblem(norbs, 4, 4)

display(problem)

nr = 1
v0 = rand(problem.dim,nr)
v0[:,1] .= 0
v0[1,1] = 1
v0 = v0 * inv(sqrt(v0'*v0))

Hmap = StringCI.get_map(ints, problem)
Random.seed!(3);
A = Diagonal(rand(20)) + .0001*rand(20,20)
A = A'+A

#davidson = FermiCG.Davidson(A,max_iter=400, nroots=nr, tol=1e-5)
davidson = FermiCG.Davidson(Hmap,v0=v0,max_iter=100, nroots=nr, tol=1e-5)
Adiag = StringCI.compute_fock_diagonal(problem,mf.mo_energy, e_mf)
#FermiCG.solve(davidson)
FermiCG.solve(davidson, Adiag=Adiag)
#FermiCG.solve(davidson, Adiag=Diagonal(A))

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using MaxMinFilters
using Test
using Random
# using Serialization

@test [] == detect_ambiguities(Base, Core)

tests = [
         "movmaxmin",
         "movmaxmin_stateful",
        ]

if length(ARGS) > 0
    tests = ARGS
end

@testset "MaxMinFilters" begin

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end

end # @testset


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

using Convex
using ECOS
using PiecewiseQuadratics
using LinearAlgebra

@testset "Only quadratic" begin
    P = ones(1, 1)
    q = -2 * ones(1)
    A = zeros(0, 1)
    b = zeros(0)
    g = zeros(PiecewiseQuadratic, 1)

    ρ = zeros(0)
    σ = ones(1)

    params = AdmmParams(P, q, A, b, g)
    settings = Settings(; ρ=ρ, σ=σ)
    vars, _ = optimize(params, settings)

    @test abs(vars.x[1] - 2) < 1e-3
    @test abs(vars.xt[1] - 2) < 1e-3
end

@testset "Only separable" begin
    P = zeros(1, 1)
    q = zeros(1)
    A = zeros(0, 1)
    b = zeros(0)
    g = [PiecewiseQuadratic([BoundedQuadratic(2.0, 3.0, 0.0, -1.0, 10.0),
                             BoundedQuadratic(3.0, 5.0, 0.0, 1.0, 4.0)])]

    ρ = zeros(0)
    σ = ones(1)

    params = AdmmParams(P, q, A, b, g)
    settings = Settings(; ρ=ρ, σ=σ)
    vars, _ = optimize(params, settings)

    display(vars)
    @test abs(vars.x[1] - 3) < 1e-3
    @test abs(vars.xt[1] - 3) < 1e-3
end

@testset "Random quadratic program" begin
    # Original generating code:
    # P = randn(5, 5)
    # P = P'*P
    # A = randn(2, 5)

    P = [1.29585 2.90775 -0.421317 1.52657 0.209435
         2.90775 10.5604 0.298516 4.64058 -2.92845
         -0.421317 0.298516 1.94608 -1.61207 -0.198397
         1.52657 4.64058 -1.61207 4.06981 -2.02592
         0.209435 -2.92845 -0.198397 -2.02592 4.66317]
    q = -10 * ones(5)
    A = [1.33694 -0.745464 -0.0531773 -2.11537 1.22246
         -0.0731486 -1.22006 -0.165136 -0.066768 0.567695]
    b = A * ones(5)
    g = repeat([indicator(0, Inf)], 5)

    x = Variable(5)
    problem = Convex.minimize(0.5 * quadform(x, P) + q' * x, [x >= 0, A * x == b])
    solve!(problem, ECOS.Optimizer)

    ρ = ones(2)
    σ = ones(5)
    params = AdmmParams(P, q, A, b, g)
    settings = Settings(; ρ=ρ, σ=σ, term_cond_freq=1000, max_iters=1000)
    vars, stats = optimize(params, settings)

    @test norm(vars.x - x.value) < 1e-3
    @test norm(vars.xt - x.value) < 1e-3
end

@testset "Random Lasso problem" begin
    n = 2  # num features
    m = 0  # num constraints

    X = [1.33694 -0.745464
         -0.0531773 -2.11537
         1.22246 -0.0731486
         -1.22006 -0.165136
         -0.066768 0.567695]
    y = X * ones(n)
    λ = 1.0

    β = Variable(n)
    problem = Convex.minimize(0.5 * sumsquares(X * β - y) + λ * sum(abs(β)))
    solve!(problem, ECOS.Optimizer)

    ρ = ones(m)
    σ = ones(n)
    P = X' * X
    q = -X' * y
    A = zeros(m, n)
    b = zeros(m)
    g = [PiecewiseQuadratic([BoundedQuadratic(-Inf, 0.0, 0.0, -λ, 0.0),
                             BoundedQuadratic(0.0, Inf, 0.0, λ, 0.0)]) for i in 1:n]
    params = AdmmParams(P, q, A, b, g)
    settings = Settings(; ρ=ρ, σ=σ, term_cond_freq=1000, max_iters=1000)
    vars, _ = optimize(params, settings)

    @test norm(vars.x - β.value) < 1e-3
    @test norm(vars.xt - β.value) < 1e-3

    # first vars term cache
    settings = Settings(; ρ=ρ, σ=σ, term_cond_freq=1000, max_iters=1000,
                        term_cond_type=FIRST_VARS_TERM_COND_FLAG, compute_stats=true)
    vars, stats = optimize(params, settings)

    @test norm(vars.x - β.value) < 1e-3
    @test norm(vars.xt - β.value) < 1e-3
end

@testset "small example" begin
    # NOTE: example from README.md

    # construct problem data (ensuring the problem is feasible)
    x0 = [0.0036165677461501566, 0.9751164348173793, 0.4825907494313493,
          0.47578842808561417]
    A = [0.477388 0.796388 0.509418 0.426954;
         0.772124 0.381946 0.412414 0.815745]
    b = A * x0
    X = [0.444226 0.331317 0.870612 0.576343;
         0.236767 0.392398 0.000545117 0.0703828;
         0.413522 0.893505 0.366298 0.654319;
         0.149137 0.0738117 0.737134 0.403253]
    P = X'X  # ensure P is positive definite
    @assert isposdef(P)
    q = [0.007228400155782522, 0.9979590358695292, 0.6284406861482683, 0.08716574428367818]

    # x1 has to be in union([-1, 2], [2.5, 3.5]) and has a quadratic penalty if
    # it lies in [-1, 2] and a linear penalty if it lies in [2.5, 3.5]
    g1 = PiecewiseQuadratic([BoundedQuadratic(-1, 2, 1, 0, 0),
                             BoundedQuadratic(2.5, 3.5, 0, 1, 0)])
    # x2 has to be between -20 and 10
    g2 = indicator(-20, 10)

    # x3 has to be between -5 and 10
    g3 = indicator(-5, 10)

    # x4 has to be exactly 1.2318
    g4 = indicator(1.2318, 1.2318)

    g = [g1, g2, g3, g4]

    # solve
    params = AdmmParams(P, q, A, b, g)

    # display
    display(params)

    vars, _ = optimize(params)

    expected = [-0.39769, 1.73815, -0.96784, 1.2318]
    @test norm(vars.x - expected) < 1e-4
end