text:
## Load Packages
using Omega
using StatsBase
using Random
using Plots
using Distributions

Random.seed!(123)

"""
Run basic Lotka-Volterra Simulation
"""

function gillespie_(rng, transitions, N, t0=0.0)

    """
    Implementation of the Gillespie algorithm for Lotka-Volterra model

    Args:
    transitions(dict): A dictionary with keys "prey" and "pred", and each element
        is a list containing the row of the transition matrix.
    N(int): Number of desired simulation iterations.
    t0(float): initial time point.

    Returns:
    A tuple of lists of the predator, prey and time trajectories
    """

    initial = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
    prey_list = [initial["prey"]]
    pred_list = [initial["pred"]]
    times = [t0]
    ecology = initial

    theta = Dict("spawn_prey" => spawn_prey(rng),
                 "prey2pred" => prey2pred(rng),
                 "pred_dies" => pred_dies(rng))
    t = times[1]

    for i = 1:N
        hazards = get_hazards(ecology, theta)
        transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
        t = t + sum(values(hazards))

        ecology["prey"] += transition[1]
        ecology["pred"] += transition[2]
        # Enforce only positive integers
        ecology["prey"] = max(1, ecology["prey"])
        ecology["pred"] = max(1, ecology["pred"])

        append!(prey_list, ecology["prey"])
        append!(pred_list, ecology["pred"])
        append!(times, t)
    end

    return (prey_list, pred_list, times)
end

function get_hazards(ecology, theta)
    """
    Compute the hazard function given the current states.  Note that in this
        model the function depends only on state, not time, though this is not
        the case in general. "spawn_prey" represents the event of a prey being born,
        "prey2pred" represents a predator consuming a new prey and consequently spawning
        a new predator, "pred_dies" represents the death of a predator.

    args:
        ecology(dict): a dictionary containing species counts.
        theta(dict): A dictionary where keys represent events (reactions) and
            values are the location hyperparameters for rate parameter constants
            corresponding to those events.
    """

    return Dict(
        "spawn_prey" => theta["spawn_prey"] * ecology["prey"],
        "prey2pred" => theta["prey2pred"] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"] * ecology["pred"]
        )
end

# Define gillespie parameters
Pre = [[1, 0], [1, 1], [0, 1]]
Post = [[2, 0], [0, 2], [0, 0]]
transition_mat = Post - Pre
transitions = Dict("spawn_prey" => transition_mat[1,],
                    "prey2pred" => transition_mat[2,],
                    "pred_dies" => transition_mat[3,])
t0=0.0
N = 1000

# Initiate starting values randomely
prey_init = normal(1., .3)
pred_init = normal(1., .3)

# Random variables for rates
spawn_prey = normal(1.0, .01)
prey2pred = normal(0.1, .0001)
pred_dies = normal(1.0, .0075)

# Make gillespie into random variable
gillespie = ciid(gillespie_, transitions, N, 0.0)

# Sample from random variables
samples = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie),
                5, alg = RejectionSample)

plot_vals = hcat(values(samples[2][6][1]), values(samples[2][6][2]))
basic_plot = plot(1:N+1, plot_vals[1:N+1,1:2],
        title = "Basic Simulation",
        label = ["Prey" "Pred"],
        lw = 1.25)
#savefig(basic_plot,"./plots/basic_simulation.png")
__________________________________________________________________________________

"""
Run a conditioned model of LV using Gillespie

The conditional statement reads as:
Condition the simulation as the mean number of prey over cond_time_range is over
the mean_limit
"""

#Define helper functions for our conditional statement
function get_prey_list(g_output)
    return g_output[1]
end

function cond_statement_func(list, N)
    mean = 0
    for n = 1:N
        mean += list[end-n]
    end
    return mean/N
end

#define initial parameters for our conditional statement
cond_time_range = 250
mean_limit = 80

#Turn gillespie into a random variable
gillespie = ciid(gillespie_, transitions, N, 0.0)

#Apply our helper functions to get a conditional statement
prey_list = lift(get_prey_list)(gillespie)
cond_statement = lift(cond_statement_func)(prey_list, cond_time_range)

#run the simulation get the conditioned trajectories
samples_cond = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie, cond_statement),
                cond_statement > mean_limit, 5, alg = RejectionSample)
samples_cond[2][7]

#plot the conditioned predator and prey trajectories
plot_vals = hcat(values(samples_cond[2][6][1]), values(samples_cond[2][6][2]))
cond_plot = plot(1:N+1, plot_vals[1:N+1,1:2],
        title = "Conditioned Model",
        label = ["Prey" "Pred"],
        lw = 1.25)
#savefig(cond_plot,"./plots/conditioned_model.png")

__________________________________________________________________________
"""
Run an Interventional model of LV using Gillespie

The interventional statement reads as:
Intervene the simulation at n_int by setting the (predator, prey) value to
(prey_int, pred_int)
"""
function gillespie_int_(rng, transitions, N, n_int, prey_int, pred_int, t0=0.0)
    """
    Implementation of the Gillespie algorithm with intervention for Lotka-Volterra model

    Args:
    transitions(dict): A dictionary with keys "prey" and "pred", and each element
        is a list containing the row of the transition matrix.
    N(int): Number of desired simulation iterations.
    n_int: time point of intervention
    prey_int: value prey is intervened to at time of intervention
    pred_int: value pred is intervened to at time of intervention
    t0(float): initial time point.

    Returns:
    A tuple of lists of the intervened predator, prey and time trajectories
    """

    initial = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
    prey_list = [initial["prey"]]
    pred_list = [initial["pred"]]
    times = [t0]
    ecology = initial

    theta = Dict("spawn_prey" => spawn_prey(rng),
                 "prey2pred" => prey2pred(rng),
                 "pred_dies" => pred_dies(rng))
    t = times[1]

    for i = 1:N
        hazards = get_hazards(ecology, theta)
        transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
        t = t + sum(values(hazards))

        ecology["prey"] += transition[1]
        ecology["pred"] += transition[2]
        # Enforce only positive integers

        if i == n_int
            if !isnothing(prey_int)
                ecology["prey"] = prey_int
            end
            if !isnothing(pred_int)
                ecology["pred"] = pred_int
            end
        end

        ecology["prey"] = max(1, ecology["prey"])
        ecology["pred"] = max(1, ecology["pred"])

        append!(prey_list, ecology["prey"])
        append!(pred_list, ecology["pred"])
        append!(times, t)
    end

    return (prey_list, pred_list, times)
end

# Define parameters for our interventional statement
t_int = 250
prey_int = 10
pred_int= nothing

# Sample from random variables
gillespie_int = ciid(gillespie_int_, transitions, N, t_int, prey_int, pred_int, 0.0)
samples_int = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie_int),
                5, alg = RejectionSample)

#plot the trajectories of the intervened model
plot_vals = hcat(values(samples_int[2][6][1]), values(samples_int[2][6][2]))
cull_prey_plot = plot(1:N+1, plot_vals[1:N+1,1:2],
        title = "Action: Cull Prey",
        label = ["Prey" "Pred"],
        lw = 1.25)
#savefig(cull_prey_plot,"./plots/cull_prey_plot.png")

#generate histograms for treatment effect
t_int = 250
prey_int = 10
pred_int= nothing
gillespie_int = ciid(gillespie_int_, transitions, N, t_int, prey_int, pred_int, 0.0)
samples_prey_cull_treat = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie, gillespie_int),
                5, alg = RejectionSample)

prey_cull_treat_effect_hist = Plots.histogram(samples_prey_cull_treat[2][7][1] - samples_prey_cull_treat[2][6][1], bins = 10,
        title = "Prey Cull Treatment Effect")
#savefig(prey_cull_treat_effect_hist,"./plots/prey_cull_treat_effect_hist.png")
t_int = 250
prey_int = nothing
pred_int= 60
gillespie_int = ciid(gillespie_int_, transitions, N, t_int, prey_int, pred_int, 0.0)
samples_pred_cull_treat = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie, gillespie_int),
                5, alg = RejectionSample)

pred_inc_treat_effect_hist = Plots.histogram(samples_prey_cull_treat[2][7][1] - samples_prey_cull_treat[2][6][1], bins = 10,
        title = "Pred Inc Treatment Effect")
#savefig(pred_inc_treat_effect_hist,"./plots/pred_inc_treat_effect_hist.png")

__________________________________________________________________________
"""
Run a Counterfactual model of LV using Gillespie

The interventional statement reads as:
Given some simulation what would have happened if at n_int the (predator, prey)
values were (prey_int, pred_int)
"""
function gillespie_count_(prey_inc, pred_inc, n_int)
    """
    Implementation of a Counterfactual Gillespie algorithm for Lotka-Volterra model

    Args:
    n_int: time point of intervention
    prey_int: value prey is intervened to at time of intervention
    pred_int: value pred is intervened to at time of intervention
    t0(float): initial time point.

    Returns:
    A gillespie function with the args defined as constants
    """
    function gillespie_(rng, transitions, N, t0=0.0)
        """
        Implementation of the Gillespie algorithm for Lotka-Volterra model

        Args:
        transitions(dict): A dictionary with keys "prey" and "pred", and each element
            is a list containing the row of the transition matrix.
        N(int): Number of desired simulation iterations.
        t0(float): initial time point.

        Returns:
        A tuple of lists of the intervened predator, prey and time trajectories
        """
        initial = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
        prey_list = [initial["prey"]]
        pred_list = [initial["pred"]]
        times = [t0]
        ecology = initial

        theta = Dict("spawn_prey" => spawn_prey(rng),
                     "prey2pred" => prey2pred(rng),
                     "pred_dies" => pred_dies(rng))
        t = times[1]
        for i = 1:N
            hazards = get_hazards(ecology, theta)
            transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
            t = t + sum(values(hazards))

            ecology["prey"] += transition[1]
            ecology["pred"] += transition[2]
            # Enforce only positive integers

            if i == n_int
                ecology["prey"] = prey_inc
                ecology["pred"] = pred_inc
            end

            ecology["prey"] = max(1, ecology["prey"])
            ecology["pred"] = max(1, ecology["pred"])

            append!(prey_list, ecology["prey"])
            append!(pred_list, ecology["pred"])
            append!(times, t)
        end
        return (prey_list, pred_list, times)
    end
end

#Define the baseline model
gillespie = ciid(gillespie_, transitions, N, 0.0)
gillispie_count_func_ = gillespie_count_(0, 20, 250)

#Apply the counterfactual step
gillispie_count = replace(gillespie, gillespie_ => gillispie_count_func_)

#Sample and run the simulation of the base and counterfactual models
samples_count = rand((prey_init, pred_init, spawn_prey, prey2pred, pred_dies, gillespie, gillispie_count),
                5, alg = RejectionSample)

#Plot the trajectories of the baseline model
plot_vals = hcat(values(samples_count[2][6][1]), values(samples_count[2][6][2]))
count_none_plot = plot(1:N+1, plot_vals[1:N+1,1:2],
        title = "No Counterfactual",
        label = ["Prey" "Pred"],
        lw = 1.25)
#Plot the trajectories of the counterfactual applied to the base model
plot_vals = hcat(values(samples_count[2][7][1]), values(samples_count[2][7][2]))
count_ind_pred_plot = plot(1:N+1, plot_vals[1:N+1,1:2],
        title = "Counterfactual: Inc Predators",
        label = ["Prey" "Pred"],
        lw = 1.25)

#savefig(count_ind_pred_plot,"./plots/count_ind_pred_plot.png")

max_stars_repo_name: neelgeek/causalML

max_issues_repo_name: neelgeek/causalML

max_forks_repo_name: neelgeek/causalML


text:
module StringDistances

using Distances
import Distances: evaluate, result_type

isnormalized(dist::SemiMetric) = false

include("utils.jl")
include("edit.jl")
include("qgram.jl")
include("modifier.jl")

const StringDistance = Union{Jaro, Levenshtein, DamerauLevenshtein, RatcliffObershelp, QGramDistance, Winkler, Partial, TokenSort, TokenSet, TokenMax}
function result_type(dist::StringDistance, s1, s2)
    typeof(evaluate(dist, "", ""))
end

"""
    compare(s1, s2, dist)

return a similarity score between 0 and 1 for the strings `s1` and
`s2` based on the distance `dist`.

### Examples
```julia-repl
julia> compare("martha", "marhta", Levenshtein())
0.6666666666666667
```
"""
compare(s1, s2, dist::StringDistance; min_score = 0.0) = 1 - evaluate(normalize(dist), s1, s2, 1 - min_score)

const metrics = [Cosine, DamerauLevenshtein, Jaccard, Jaro, Levenshtein, Overlap, Partial, QGram,
    RatcliffObershelp, SorensenDice, StringDistance,  TokenMax,
    TokenSet, TokenSort, Winkler]

for M in metrics
    @eval @inline (dist::$M)(a::AbstractString, b::AbstractString) = evaluate(dist, a, b)
end

include("find.jl")

##############################################################################
##
## Export
##
##############################################################################

export
StringDistance,
Levenshtein,
DamerauLevenshtein,
Jaro,
RatcliffObershelp,
QGram,
Cosine,
Jaccard,
SorensenDice,
Overlap,
Winkler,
Partial,
TokenSort,
TokenSet,
TokenMax,
evaluate,
compare,
result_type,
qgrams,
normalize
end

##############################################################################
##
## Some things about Strings

# length: number of characters
# ncodeunits: Return the number of code units in a string (aking to index of vector).
# Not all such indices  are valid – they may not be the start of a character,.
# sizeof:  Size, in bytes, of the string str. Equal to the number of code units in str
# multiplied by the size, in bytes, of one code unit in str.

# lastindex: Return the last index of a collection
# nextinds(s, i):  return the index of the start of the character whose encoding starts after index i
# nextind(s, 0, N): return the index of the Nth character of s (or, if there are
# less than N characters, return ncodeunits(str) + (N - length(s))

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

ext: jl


text:
include("TestInstance.jl")
include("TestQueueRequest.jl")
include("TestDevice.jl")
include("TestSampler.jl")
include("TestFramebuffer.jl")
include("TestGraphicsPipeline.jl")
include("TestComputePipeline.jl")
include("TestImage.jl")
include("TestBuffer.jl")
include("TestDescriptorSet.jl")

instanceRef = Ref{lava.Instance}()
glfwRef = Ref{features.GlfwOutput}()
@test testInstance(instanceRef, glfwRef)

@test testQueueRequest()
queues = [lava.createGraphics(lava.QueueRequest, "graphics")]

deviceRef = Ref{lava.Device}()
@test testDevice(instanceRef[], queues, deviceRef)
@test testDevice(instanceRef[], queues, lava.NthGroupStrategy(0))

pipelineRef = Ref{lava.GraphicsPipeline}()
passRef = Ref{lava.RenderPass}()
@test testGraphicsPipelineFlag()
@test testCreateGraphicsPipeline(deviceRef[], glfwRef[], pipelineRef, passRef)
@test testComputePipeline(deviceRef[])

@test testSampler(deviceRef[].getLogicalDevice())

@test testImageCreateInfo(deviceRef[].getPhysicalDevice())
@test testFramebuffer(deviceRef[], passRef[])

@test testCreateBufferCreateInfo()
@test testBufferFunction(deviceRef[])

@test testDescriptorSet(deviceRef[])


text:
"""
given appropriate constants related to slice like rate of true pasitives, fals positive, false negative
and false negative
it calculates wanted  segmentation metrics
in order to reduce thread divergence and enable easy omitting this step - in case we want only global statistics
we will  calculate it on  in a seperate kernel than the data about true positives ... is calculated
"""
function calculateBasicMetricsPerSlice(tp,tn,fp,fn
                                        ,intermediateResTp
                                        ,intermediateResFp
                                        ,intermediateResFn
                                        ,)

    dice, jaccard, gce
end

MainOverlap.dice(tp,fp, fn)
MainOverlap.jaccard(tp,fp, fn)
MainOverlap.gce(tp,fp, fn)
RandIndex.calculateAdjustedRandIndex(tn,tp,fp, fn)
ProbabilisticMetrics.calculateCohenCappa(tp,fp, fn )
VolumeMetric.getVolumMetric(tp,fp, fn )
InformationTheorhetic.mutualInformationMetr(tn,tp,fp, fn)
InformationTheorhetic.variationOfInformation(tn,tp,fp, fn)


text:
using NonconvexMMA

abstract type AbstractOptimizer end

@params mutable struct Optimizer{Tdev} <: AbstractOptimizer
    model::AbstractModel
    alg::Any
    workspace::Any
    dev::Tdev
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::Optimizer) = println("TopOpt optimizer")

multol!(o::Optimizer, m::Real) = o.options.tol *= m
function setbounds!(o::Optimizer, x, w)
    o.model.box_min .= max.(0, x .- w)
    return o.model.box_max .= min.(1, x .+ w)
end

function Optimizer(
    obj,
    constr,
    vars,
    opt=MMA87(),
    device::Tdev=CPU();
    options=MMAOptions(),
    convcriteria=KKTCriteria(),
) where {Tdev<:AbstractDevice}
    T = eltype(vars)
    nvars = length(vars)
    x0 = copy(vars)
    model = Nonconvex.Model(obj)
    addvar!(model, zeros(T, nvars), ones(T, nvars))
    add_ineq_constraint!(model, Nonconvex.FunctionWrapper(constr, length(constr(x0))))
    workspace = NonconvexMMA.Workspace(
        NonconvexCore.tovecmodel(model)[1],
        opt,
        x0;
        options=options,
        convcriteria=convcriteria,
    )
    return Optimizer(model, opt, workspace, device)
end

Utilities.getpenalty(o::Optimizer) = getpenalty(o.solver)
function Utilities.setpenalty!(o::Optimizer, p)
    setpenalty!(o.solver, p)
    return o
end

function (o::Optimizer)(x0::AbstractVector)
    @unpack workspace = o
    @unpack options, model = workspace
    Nonconvex.NonconvexCore.reset!(workspace, x0)
    setoptions!(workspace, options)
    mma_results = Nonconvex.NonconvexCore.optimize!(workspace)
    return mma_results
end
function setoptions!(workspace::Nonconvex.Workspace, options)
    workspace.options = options
    return workspace
end

# For adaptive SIMP
function (o::Optimizer)(workspace::Nonconvex.Workspace)
    mma_results = Nonconvex.NonconvexCore.optimize!(workspace)
    return mma_results
end

@params struct MMAOptionsGen
    maxiter::Any
    outer_maxiter::Any
    tol::Any
    s_init::Any
    s_incr::Any
    s_decr::Any
    store_trace::Any
    show_trace::Any
    auto_scale::Any
    dual_options::Any
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::MMAOptionsGen)
    return println("TopOpt MMA options generator")
end
function (g::MMAOptionsGen)(i)
    return MMAOptions(;
        maxiter=g.maxiter(i),
        outer_maxiter=g.outer_maxiter(i),
        tol=g.tol(i),
        s_init=g.s_init(i),
        s_incr=g.s_incr(i),
        s_decr=g.s_decr(i),
        store_trace=g.store_trace(i),
        show_trace=g.show_trace(i),
        auto_scale=g.auto_scale(i),
        dual_options=g.dual_options(i),
    )
end

function (g::MMAOptionsGen)(options, i)
    return MMAOptions(;
        maxiter=optionalcall(g, :maxiter, options, i),
        outer_maxiter=optionalcall(g, :outer_maxiter, options, i),
        tol=optionalcall(g, :tol, options, i),
        s_init=optionalcall(g, :s_init, options, i),
        s_incr=optionalcall(g, :s_incr, options, i),
        s_decr=optionalcall(g, :s_decr, options, i),
        store_trace=optionalcall(g, :store_trace, options, i),
        show_trace=optionalcall(g, :show_trace, options, i),
        auto_scale=optionalcall(g, :auto_scale, options, i),
        dual_options=optionalcall(g, :dual_options, options, i),
    )
end
function optionalcall(g, s, options, i)
    if getproperty(g, s) isa Nothing
        return getproperty(options, s)
    else
        return getproperty(g, s)(i)
    end
end


text:
using DataStructures

function main()
    K = parse(Int, readline())
    INF = 10^9

    dists = INF * ones(Int, K)
    dists[2] = 1
    que = Deque{Int}()
    pushfirst!(que, 1)

    while !isempty(que)
        v = popfirst!(que)

        v′ = 10v % K
        if dists[v′+1] > dists[v+1]
            dists[v′+1] = dists[v+1]
            pushfirst!(que, v′)
        end

        v′ = (v + 1) % K
        if dists[v′+1] > dists[v+1] + 1
            dists[v′+1] = dists[v+1] + 1
            push!(que, v′)
        end
    end

    println(dists[1])
end

main()


text:
using MPI

juliapath = joinpath(Sys.BINDIR, "julia")

for np in 1:4
	@info "tests with np = $np"
	cmd = `$(MPI.mpiexec_path) -np $np $juliapath --startup=no --project tests.jl`
	run(cmd)
end


text:
export find_treewidth_from_order
export find_similar_groups
export graph_from_gr, graph_to_gr, load_example_graph

###
### IO functions for reading and writing graphs in gr format.
###

"""
    graph_from_gr(filename::String)

Read a graph from the provided gr file.
"""
function graph_from_gr(filename::String)
    lines = readlines(filename)

    # Create a Graph with the correct number of vertices.
    num_vertices = parse.(Int, split(lines[1], ' ')[3])
    G = lg.SimpleGraph(num_vertices)

    # Add an edge to the graph for every other line in the file.
    for line in lines[2:end]
        src, dst = parse.(Int, split(line, ' '))
        lg.add_edge!(G, src, dst)
    end

    G
end

"""
    graph_to_gr(G::AbstractGraph, filename::String)

Write the provided graph to a file in gr format.
"""
function graph_to_gr(G::lg.AbstractGraph, filename::String)
    open(filename, "w") do io
        write(io, "p tw $(nv(G)) $(ne(G))\n")
        for e in lg.edges(G)
            write(io, "$(e.src) $(e.dst)\n")
        end
    end
end

"""
    load_example_graph(graph_name::String="test")

Return one of the example graphs from the `VertexEliminationOrders` package.

# Possible arguments:
- bristlecone\\_48\\_1-16-1\\_0
- rectangular\\_2x2\\_1-2-1\\_0
- rectangular\\_4x4\\_1-16-1\\_0
- sycamore\\_53\\_8\\_0
- sycamore\\_53\\_20\\_0
- test
"""
function load_example_graph(graph_name::String="test")
    filename = joinpath(dirname(dirname(@__FILE__)), "examples", graph_name * ".gr")
    graph_from_gr(filename)
end

###
### functions to analyse elimination orders.
###

"""
    find_treewidth_from_order(G::AbstractGraph, order::Array{Symbol, 1})

Return the treewidth of `G` with respect to the elimination order in `order`.
"""
function find_treewidth_from_order(G::lg.AbstractGraph, order::Array{<:Integer, 1})
    G = deepcopy(G)
    labels = collect(1:lg.nv(G))
    τ = 0
    for v_label in order
        v = findfirst(vl -> vl == v_label, labels)
        τ = max(τ, lg.degree(G, v))
        lg_eliminate!(G, labels, v)
    end
    τ
end

###
### Functions to analyse a graph.
###

"""Return a vector containing the similar groups of G as defined by Yaun_2011"""
function find_similar_groups(G::lg.SimpleGraph{Int})
    groups = Vector{Int}[]
    for v in lg.vertices(G)
        v_in_a_group = false
        for group in groups
            if are_similar(G, group[1], v)
                push!(group, v)
                v_in_a_group = true
                break
            end
        end
        v_in_a_group || push!(groups, [v])
    end

    groups
end

function are_similar(G::lg.SimpleGraph, u::Int, v::Int)
    return Set(setdiff(lg.all_neighbors(G, u), [v])) == Set(setdiff(lg.all_neighbors(G, v), [u]))
end

"""
    cliqueness(G::AbstractGraph, v::Integer)
Return the number of edges that need to be added to `G` in order to make the neighborhood of
vertex `v` a clique.
"""
function cliqueness(g::lg.AbstractGraph, v::Integer)::Int
    neighborhood = lg.all_neighbors(g, v)::Array{Int64, 1}
    count = 0
    for i in 1:length(neighborhood)-1
        for j in i+1:length(neighborhood)
            vi = neighborhood[i]
            ui = neighborhood[j]
            if !lg.has_edge(g, ui, vi)::Bool
                count += 1
            end
        end
    end
    count
end

###
### Functions to modify a graph.
###

"""
    eliminate!(graph, labels, vertex)

Connects the neighbours of the given vertex into a clique before removing
it from the graph.

The array of vertex labels are also updated to reflect the
reording of vertex indices when the graph is updated.
"""
function lg_eliminate!(g::lg.AbstractGraph, labels, v)
    Nᵥ = lg.all_neighbors(g, v)::Array{Int64, 1}
    for i = 1:length(Nᵥ)-1
        vi = Nᵥ[i]
        for j = i+1:length(Nᵥ)
            vj = Nᵥ[j]
            lg.add_edge!(g, vi, vj)
        end
    end

    lg.rem_vertex!(g, v)
    labels[v] = labels[end]
    pop!(labels)
    g
end

"""
    eliminate!(graph, labels, vertex)
Connects the neighbours of the given vertex into a clique before removing
it from the graph.
The arrays of vertex labels and cliqueness are also updated to reflect the
reording of vertex indices when the graph is updated.
"""
function lg_eliminate!(g::lg.AbstractGraph, labels, c_map, v)
    Nᵥ = lg.all_neighbors(g, v)::Array{Int64, 1}
    for i = 1:length(Nᵥ)-1
        vi = Nᵥ[i]
        for j = i+1:length(Nᵥ)
            vj = Nᵥ[j]

            # Try add an edge connecting vi and ui. If successful, update `c_map`.
            edge_added = lg.add_edge!(g, vi, vj)
            if edge_added
                # Common neighbours of vi and vj have one less edge to add
                # when being eliminated after vi and ui are connected.
                Nvi = lg.all_neighbors(g, vi)::Array{Int64, 1}
                Nvj = lg.all_neighbors(g, vj)::Array{Int64, 1}
                for n in Nvi
                    if n in Nvj
                        c_map[n] -= 1
                    end
                end

                # ui and vi are now neighbours so their cliqueness may increase.
                for n in Nvi
                    if !(n == vj) && !(lg.has_edge(g, n, vj)::Bool)
                        c_map[vi] += 1
                    end
                end
                for n in Nvj
                    if !(n == vi) && !(lg.has_edge(g, n, vi)::Bool)
                        c_map[vj] += 1
                    end
                end
            end
        end
    end

    # Removing v from G means it's also removed from its neighbour's neighbourhood, so their
    # cliqueness may be reduced.
    for n in Nᵥ
        Nₙ = lg.all_neighbors(g, n)::Array{Int64, 1}
        for u in Nₙ
            if !(u == v)
                if !lg.has_edge(g, v, u)
                    c_map[n] -= 1
                end
            end
        end
    end

    lg.rem_vertex!(g, v)
    c_map[v] = c_map[end]
    labels[v] = labels[end]
    g
end


text: Dict{String,Any}("a b" => Dict{String,Any}("value" => "1","type" => "integer"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
export proxL2!, normL2

"""
    proxL2!(x::Vector{T}, λ::Float64; kargs...) where T

proximal map for Tikhonov regularization.
"""
function proxL2!(x::T, λ::Float64; kargs...) where T<:AbstractArray
  x[:] .*= 1. / (1. + 2. *λ)#*x
end

"""
    normL2(x::Vector{T}, λ::Float64, kargs...)

returns the value of the L2-regularization term
"""
normL2(x::T, λ::Float64; kargs...) where T<:AbstractArray = λ*norm(x,2)^2


text:
import LinearAlgebra: norm

function get_errors(BIEsol, RCWAsol)
    x, y₁, y₂ = BIEsol.x, BIEsol.y₁, BIEsol.y₂
    err = [
        getproperty(BIEsol, sym) - getproperty(RCWAsol, sym)
        for sym in (:I₁, :O₁, :I₂, :O₂)
    ]
    (x=x, y₁=y₁, y₂=y₂, I₁=err[1], O₁=err[2], I₂=err[3], O₂=err[4])
end

function L2cubature(sol; dims)
    margL2 = [
        mapslices(norm, e; dims=dims) ./ sqrt(prod(size(e)[dims]))
        for e in (sol.I₁, sol.O₁, sol.I₂, sol.O₂)
    ]
    (I₁=margL2[1], O₁=margL2[2], I₂=margL2[3], O₂=margL2[4])
end

function errorL2(solerrors)
    err = L2cubature(solerrors; dims=1:2)
    (I₁=err.I₁[1], O₁=err.O₁[1], I₂=err.I₂[1], O₂=err.O₂[1])
end
function errorbydistance(solerrors)
    err = L2cubature(solerrors; dims=1)
    (y₁=solerrors.y₁, y₂=solerrors.y₂, I₁=vec(err.I₁), O₁=vec(err.O₁), I₂=vec(err.I₂), O₂=vec(err.O₂))
end


text:
using BenchmarkTools, StaticArrays, Parameters, NamedTuples, LinearAlgebra
#Generic functions used by all variations
f(p) = p.a + 1.0
g(f,p) = f(p)
h(x) = x^2
i(f, p) = f(p) * sum(h.(p.c))
function j(f, p)
    @unpack c, b, d = p #Benchmarked that @unpack is identical
    g(f,p) .+ sum(h.(c)) .+ b .+ c * d
end
k(p) = sum(p.c) + p.a * p.b

#Higher-order function function
function l(p)
    return x-> x^(p.b) #Returns a function
end

#Uses the higher-order function
function m(p, l2)
    @unpack a, b, e = p
    f2 = l2(p) #Uses the higher order function
    return f2(a) + h(b) + e(b)
end

#Checking if there is a higher order function penalty
n1(p) = 1.1^(p.b)
function n2(p)
    n3 = l(p)
    return n3(1.1)
end

#Named Tuple
np = (a=1.0, b=2, c=[1 2 3], e=x -> x^3, d=@SVector [1, 2, 3])

#v0.6 named tuple
np2 = @NT(a=1.0, b=2, c=[1 2 3], e=x -> x^3, d=@SVector [1, 2, 3])

#Typed Structure
struct X{F}
    a::Float64
    b::Int64
    c::Array{Int64,2}
    e::F
    d::SArray{Tuple{3},Int64,1,3}
end
np3 = X(1.0, 2, [1 2 3], x -> x^3, @SVector [1, 2, 3])
#Untyped structure
struct X2
    a
    b
    c
    e
    d
end
np4 = X2(1.0, 2, [1 2 3], x -> x^3, @SVector [1, 2, 3])

#Trying the variations
@btime i($f,$np)
@btime i($f,$np2)
@btime i($f,$np3)
@btime i($f,$np4)

@btime j($f,$np)
@btime j($f,$np2)
@btime j($f,$np3)
@btime j($f,$np4)

@btime k($np)
@btime k($np2)
@btime k($np3)
@btime k($np4)

@btime m($np, $l)
@btime m($np2, $l)
@btime m($np3, $l)
@btime m($np4, $l)

#Checking if overhead from higher order functions
@btime n1($np)
@btime n2($np)
@btime n1($np2)
@btime n2($np2)
@btime n1($np3)
@btime n2($np3)
@btime n1($np4)
@btime n2($np4)


text:
using POMDPModels
using POMDPToolbox

let
    rng = MersenneTwister(7)

    mdp = GridWorld()
    s = initial_state(mdp, rng)
    a = rand(rng, actions(mdp))
    @inferred generate_sri(mdp, s, a, rng)

    pomdp = TigerPOMDP()
    s = initial_state(pomdp, rng)
    a = rand(rng, actions(pomdp))
    @inferred generate_sori(pomdp, s, a, rng)

    policy = RandomPolicy(pomdp, rng=rng)
    @inferred action_info(policy, s)

    solver = RandomSolver(rng=rng)
    policy, sinfo = solve_info(solver, pomdp)
    @test isa(sinfo, Void)

    up = updater(policy)
    d = initial_state_distribution(pomdp)
    b = initialize_belief(up, d)
    a = action(policy, b)
    sp, o = generate_so(pomdp, rand(rng, d), a, rng)
    @inferred update_info(up, b, a, o)
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()
config_with_basis = MOIT.TestConfig(basis = true)

@testset "Scalar functionize" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.Constraint.ScalarFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    sx = MOI.SingleVariable(x)
    sy = MOI.SingleVariable(y)
    ci = MOI.add_constraint(bridged_mock, sx, MOI.GreaterThan(0.0))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ sx
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, sy)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ sy
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(0.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.GreaterThan(1.0))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(1.0)
    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),),
    )

    MOIT.basic_constraint_tests(
        bridged_mock,
        config,
        include = [
            (F, S) for F in [MOI.SingleVariable],
            S in [MOI.GreaterThan{Float64}, MOI.LessThan{Float64}]
        ],
    )

    for T in [Int, Float64], S in [MOI.GreaterThan{T}, MOI.LessThan{T}]
        @test MOIB.added_constraint_types(
            MOIB.Constraint.ScalarFunctionizeBridge{T,S},
        ) == [(MOI.ScalarAffineFunction{T}, S)]
    end

    @testset "linear2" begin
        MOI.empty!(bridged_mock)
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                [1, 0],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1],
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0, 1],
                con_basis = [
                    (
                        MOI.ScalarAffineFunction{Float64},
                        MOI.LessThan{Float64},
                    ) => [MOI.NONBASIC],
                    (
                        MOI.ScalarAffineFunction{Float64},
                        MOI.GreaterThan{Float64},
                    ) => [MOI.BASIC, MOI.NONBASIC],
                ],
            ),
        )
        MOIT.linear2test(bridged_mock, config_with_basis)

        cis = MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.SingleVariable,
                MOI.GreaterThan{Float64},
            }(),
        )

        @testset "$attr" for attr in [
            MOI.ConstraintPrimalStart(),
            MOI.ConstraintDualStart(),
        ]
            for ci in cis
                @test MOI.supports(bridged_mock, attr, typeof(ci))
                MOI.set(bridged_mock, attr, ci, 2.0)
                @test MOI.get(bridged_mock, attr, ci) == 2.0
            end
        end

        attr = MOIT.UnknownConstraintAttribute()
        err = ArgumentError(
            "Bridge of type `$(MathOptInterface.Bridges.Constraint.ScalarFunctionizeBridge{Float64,MathOptInterface.GreaterThan{Float64}})` does not support setting the attribute `$attr` because `MOIB.Constraint.invariant_under_function_conversion($attr)` returns `false`.",
        )
        @test_throws err MOI.set(bridged_mock, attr, ci, 1.0)
        err = ArgumentError(
            "Bridge of type `$(MathOptInterface.Bridges.Constraint.ScalarFunctionizeBridge{Float64,MathOptInterface.GreaterThan{Float64}})` does not support accessing the attribute `$attr` because `MOIB.Constraint.invariant_under_function_conversion($attr)` returns `false`.",
        )
        @test_throws err MOI.get(bridged_mock, attr, ci)

        @testset "delete" begin
            for (i, ci) in enumerate(cis)
                test_delete_bridge(
                    bridged_mock,
                    ci,
                    2,
                    ((
                        MOI.ScalarAffineFunction{Float64},
                        MOI.GreaterThan{Float64},
                        0,
                    ),),
                    num_bridged = 3 - i,
                )
            end
        end
    end
end

@testset "Vector functionize" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.Constraint.VectorFunctionize{Float64}(mock)
    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    z = MOI.add_variable(bridged_mock)
    v1 = MOI.VectorOfVariables([x, y, z])
    v2 = MOI.VectorOfVariables([x, y, y])
    ci = MOI.add_constraint(bridged_mock, v1, MOI.PowerCone(0.1))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v1
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, v2)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ v2
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.1)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.PowerCone(0.2))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.PowerCone(0.2)
    test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((MOI.VectorAffineFunction{Float64}, MOI.PowerCone{Float64}, 0),),
    )

    MOIT.basic_constraint_tests(
        bridged_mock,
        config,
        include = [
            (F, S) for F in [MOI.VectorOfVariables],
            S in [MOI.Nonnegatives, MOI.Nonpositives]
        ],
    )

    for T in [Int, Float64], S in [MOI.Nonnegatives, MOI.Nonpositives]
        @test MOIB.added_constraint_types(
            MOIB.Constraint.VectorFunctionizeBridge{T,S},
        ) == [(MOI.VectorAffineFunction{T}, S)]
    end

    @testset "lin1v" begin
        mock.optimize! =
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                [1.0, 0.0, 2.0],
                (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                    [[0, 2, 0]],
                (MOI.VectorAffineFunction{Float64}, MOI.Zeros) =>
                    [[-3, -1]],
            )
        MOIT.lin1vtest(bridged_mock, config)
        ci = first(
            MOI.get(
                bridged_mock,
                MOI.ListOfConstraintIndices{
                    MOI.VectorOfVariables,
                    MOI.Nonnegatives,
                }(),
            ),
        )
        func = MOI.get(bridged_mock, MOI.ConstraintFunction(), ci)
        MOI.delete(bridged_mock, func.variables[2])
        new_func = MOI.VectorOfVariables(func.variables[[1, 3]])
        @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) == new_func
        @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) ==
              MOI.Nonnegatives(2)

        @testset "$attr" for attr in [
            MOI.ConstraintPrimalStart(),
            MOI.ConstraintDualStart(),
        ]
            @test MOI.supports(bridged_mock, attr, typeof(ci))
            MOI.set(bridged_mock, attr, ci, [1.0, 2.0])
            @test MOI.get(bridged_mock, attr, ci) == [1.0, 2.0]
        end

        test_delete_bridge(
            bridged_mock,
            ci,
            2,
            ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),),
        )
    end
end


text:
using JLD2
using StatsPlots
using DataFrames

@load "paper/box_all.jld2" box_all_df
@load "paper/box_at.jld2" box_at_df
@load "paper/lg_all.jld2" lg_all_df

@df box_at_df scatter(:n, :mean)
p = @df box_all_df plot(:n, :mean, group=:type, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks=200:200:1600, size=(450,250),
  legend=:topleft)

savefig(p, "paper/simpleComparisonATASBS.pdf")

p = @df lg_all_df plot(:n, :mean, group=:type, marker=4,
  xlabel="T", ylabel="mean coupling time", xticks=200:200:1600, size=(450,250),
  legend=:topleft)

savefig(p, "paper/simpleComparison_lg_ATASBS.pdf")


text:
using CleanCMB
using Healpix
using Optim, LinearAlgebra
using Printf, CSV, DataFrames
using Mmap
using Statistics
using Tables
dir = "../data"
# %% Simulation parameters
nrz = 1500
ℓmin, ℓmax = 30, 260 # ℓ range for fitting
ℓko = 60 # remove information from ℓ<=ℓko
Alens = 1
Δℓ = 10 # multipole binning size
# %% Foreground cleaning parameters
showβ = false # Show the fitted foreground parameters (for the smoothed covariance)?
β0 = βs0, βd0, Td0 = [-3.0, 1.6, 19.6] # starting foreground parameters for minimisation of -log(likelihood) by `res = optimize(func, [βs0, βd0, Td0])`
ℓswitch = 30 # the multipole below which the foreground parameters are fitted for each band-power
smooth_FWHM = 0.5 # smoothing for the covariance matrix in units of degrees
σβ = σβs, σβd, σTd = [0.5, 0.5, 5] # Gaussian priors: -2*log(likelihood) = (x[1] - βs0)^2 / σβs^2 + (x[2] - βd0)^2 / σβd^2 + (x[3] - Td0)^2 / σTd^2
Σβ = Diagonal(σβ .^ 2)
lnlike_fgprior(x) = -0.5 * (x .- β0)' * (Σβ \ (x .- β0))
# %% Specification of the experiments
ν = [27, 39, 93, 145, 225, 280, 350, 410, 850] # in GHz
νused = [true, true, true, true, true, true, true, true, true] # Which frequencies to use for fitting?
#νused = [true, true, true, true, true, true, false, false, false]
nν = length(ν)
kν = findall(x -> x == true, νused)
FWHM = [91, 63, 30, 17, 11, 9, 0.58, 0.5, 0.23] # in arcmin
σ = FWHM * π / 10800 / √(8 * log(2)) # in radians
uKarcmin = [35, 21, 2.6, 3.3, 6.3, 16, 105, 372, 5.7e5] # in μK arcmin (for temperature; x√2 for pol)
lknee = [30, 30, 50, 50, 70, 100, 700, 700, 700]
αknee = [-2.4, -2.4, -2.5, -3, -3, -3, -1.4, -1.4, -1.4]
# polarisation beam, Eq.(5.8) of Ng & Liu, Int.J.Mod.Phys.D, 8, 61 (1999)
bPl(ℓ, σb) = ifelse(ℓ >= 2, exp(-(ℓ * (ℓ + 1) - 4) * σb^2 / 2), 0)
# %% Read in the hits map and calculate fsky need for the likelihood
nhitsfile = joinpath(dir, "nhits_SAT_r7.FITS")
nhits = readMapFromFITS(nhitsfile, 1, Float64)
nside = nhits.resolution.nside
nhits /= maximum(nhits)
fsky = mean(nhits)^2 / mean(nhits .^ 2)
# %% Read in binned theory power spectra
clt_th = CSV.read(joinpath(dir, "tensor_eebb_binned.csv"), DataFrame)
cls_th = CSV.read(joinpath(dir, "scalar_eebb_binned.csv"), DataFrame)
ell_eff = clt_th.leff
nbands = length(ell_eff)
rclass = 0.01
# %% Read in covariance matrices and perform parametric maxiimum likelihood fitting
io = open(joinpath(dir, "cov_bb_foreg.dat"), "r")
cov3 = Mmap.mmap(io, Array{Float64,3}, (nν, nν, nbands))
close(io)
βout = zeros(3, nrz)
for irz = 1:nrz
    ioin = open(joinpath(dir, @sprintf("cov_bb_total_irz%04d.dat", irz)), "r")
    cov1 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    ioin = open(joinpath(dir, @sprintf("cov_bb_noise_irz%04d.dat", irz)), "r")
    cov2 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    nij = zeros(nν, nν, nbands) # Noise covariance matrix for the ML analysis
    for iν = 1:nν
        ell_fact = ell_eff .* (ell_eff .+ 1) / 2π
        nij[iν, iν, :] =
            2 *
            (uKarcmin[iν] * π / 10800)^2 *
            (1 .+ (ell_eff / lknee[iν]) .^ αknee[iν]) .* ell_fact .*
            bPl.(ell_eff, σ[iν]) .^ -2
    end
    ## Knock out low multipole data of FYST
    iko = findall(x -> x <= ℓko, ell_eff)
    for iν = 7:9
        nij[iν, iν, iko] .= Inf
    end
    ## Use Parametric Maximum Likelihood method to obtain power spectra of clean maps of the CMB
    A(x) = [ones(length(kν)) synch.(ν[kν], βs = x[1]) dust1.(
        ν[kν],
        βd = x[2],
        Td = x[3],
    )] # Frequency response matrix
    # For ℓ > ℓswitch: Apply parametric maximum likelihood method using a smoothed covariance matrix.
    func_sum(x, cij) =
        -lnlike_fgprior(x) -
        fsky *
        Δℓ *
        sum(
            (2 * ell_eff[jb] + 1) *
            loglike_beta(nij[kν, kν, jb], A(x), cij[kν, kν, jb]) for
            jb = 1:nbands
        ) # -log(likelihood) to minimise by `optimize`
    # Smooth covariance matrices to `smooth_FWHM` resolution
    cij = zeros(nν, nν, nbands)
    cov1ko = copy(cov1)
    cov2ko = copy(cov2)
    σsmo = smooth_FWHM * π / 180 / sqrt(8 * log(2))
    ## Knock out low multipole data of FYST
    for iν = 7:9
        cov1ko[iν, iν, iko] =
            (cov1[iν, iν, iko] .- cov2[iν, iν, iko]) .+ 1e4 * cov2[iν, iν, iko]
        cov2ko[iν, iν, iko] .*= 1e4
    end
    for ib = 1:nbands
        cij[:, :, ib] = cov1ko[:, :, ib] * bPl(ell_eff[ib], σsmo)^2
    end
    res = optimize(x -> func_sum(x, cij), β0)
    if showβ
        println("Fitted parameters: B-mode")
        @show res
    end
    β = Optim.minimizer(res)
    βout[1:3, irz] = β[1:3]
end
# %% Save to the file
t = Tables.table([1:nrz βout[1,:] βout[2,:] βout[3,:]])
CSV.write("beta_results_removelowell.csv", t, header = ["irz", "beta_s", "beta_d", "T_d"])


text:
# Font handling when pango and fontconfig are not available.

# Define this even if we're not calling pango, since cairo needs it.
const PANGO_SCALE = 1024.0

# Serialized glyph sizes for commont fonts.
const glyphsizes = open(fd -> JSON.parse(readstring(fd)),
                        joinpath(dirname(@__FILE__), "..", "data",
                                 "glyphsize.json"))

# It's better to overestimate text extents than to underestimes, since the later
# leads to overlaping where the former just results in some extra space. So, we
# scale estimated text extends by this number. Width on the other hand tends be
# overestimated since it doesn't take kerning into account.
const text_extents_scale_x = 1.0
const text_extents_scale_y = 1.0

# Normalized Levenshtein distance between two strings.
function levenshtein(a::AbstractString, b::AbstractString)
    a = replace(lowercase(a), r"\s+", "")
    b = replace(lowercase(b), r"\s+", "")
    n = length(a)
    m = length(b)
    D = zeros(UInt, n + 1, m + 1)

    D[:,1] = 0:n
    D[1,:] = 0:m

    for i in 2:n, j in 2:m
        if a[i - 1] == b[j - 1]
            D[i, j] = D[i - 1, j - 1]
        else
            D[i, j] = 1 +  min(D[i - 1, j - 1], D[i, j - 1], D[i - 1, j])
        end
    end

    return D[n,m] / max(n, m)
end

# Find the nearst typeface from the glyph size table.
let
    matched_font_cache = Dict{AbstractString, AbstractString}()
    global match_font

    function match_font(families::AbstractString)
        if haskey(matched_font_cache, families)
            return matched_font_cache[families]
        end

        smallest_dist = Inf
        best_match = "Helvetica"
        for family in [lowercase(strip(family, [' ', '"', '\''])) for family in split(families, ',')]
            for available_family in keys(glyphsizes)
                d = levenshtein(family, available_family)
                if d < smallest_dist
                    smallest_dist = d
                    best_match = available_family
                end
            end
        end

        matched_font_cache[families] = best_match
        return best_match
    end
end

# Approximate width of a text in millimeters.
#
# Args:
#   widths: A glyph width table from glyphsizes.
#   text: Any string.
#   size: Font size in points.
#
# Returns:
#   Approximate text width in millimeters.
#
function text_width(widths::Dict, text::AbstractString, size::Float64)
    stripped_text = replace(text, r"<[^>]*>", "")
    width = 0
    for c in stripped_text
        width += get(widths, string(c), widths["w"])
    end
    width
end

function max_text_extents(font_family::AbstractString, size::Measure,
                          texts::AbstractString...)
    if !isa(size, AbsoluteLength)
        error("text_extents requries font size be in absolute units")
    end
    scale = size / 12pt
    font_family = match_font(font_family)
    glyphheight = glyphsizes[font_family]["height"]
    widths = glyphsizes[font_family]["widths"]

    fontsize = size/pt
    chunkwidths = Float64[]
    textheights = Float64[]
    for text in texts
        textheight = 0.0
        for chunk in split(text, '\n')
            chunkheight = glyphheight
            if match(r"<su(p|b)>", chunk) != nothing
                chunkheight *= 1.5
            end
            textheight += chunkheight
            push!(chunkwidths, text_width(widths, chunk, fontsize))
        end
        push!(textheights, textheight)
    end

    width = maximum(chunkwidths)
    height = maximum(textheights)
    (text_extents_scale_x * scale * width * mm,
     text_extents_scale_y * scale * height * mm)
end

function text_extents(font_family::AbstractString, size::Measure, texts::AbstractString...)
    scale = size / 12pt
    font_family = match_font(font_family)
    glyphheight = glyphsizes[font_family]["height"]
    glyphwidths = glyphsizes[font_family]["widths"]
    fontsize = size/pt

    extents = Array{@compat Tuple{Measure, Measure}}(length(texts))
    for (i, text) in enumerate(texts)
        chunkwidths = Float64[]
        textheight = 0.0
        for chunk in split(text, "\n")
            chunkheight = glyphheight
            if match(r"<su(p|b)>", text) != nothing
                chunkheight *= 1.5
            end
            textheight += chunkheight
            push!(chunkwidths, text_width(glyphwidths, chunk, fontsize))

        end
        width = maximum(chunkwidths)
        extents[i] = (text_extents_scale_x * scale * width * mm,
                      text_extents_scale_y * scale * textheight * mm)
    end

    return extents
end

# Amazingly crude fallback to parse pango markup into svg.
function pango_to_svg(text::AbstractString)
    pat = r"<(/?)\s*([^>]*)\s*>"
    input = convert(Array{UInt8}, text)
    output = IOBuffer()
    lastpos = 1

    baseline_shift = 0.0
    open_tag = false

    for mat in eachmatch(pat, text)
        write(output, input[lastpos:mat.offset-1])

        closing_tag = mat.captures[1] == "/"

        if open_tag && !closing_tag
            write(output, "</tspan>")
        end

        if mat.captures[2] == "sup"
            if mat.captures[1] == "/"
                write(output, "</tspan>")
            else
                # write(output, "<tspan style=\"dominant-baseline:inherit\" baseline-shift=\"super\">")
                write(output, "<tspan style=\"dominant-baseline:inherit\" dy=\"-0.6em\" font-size=\"83%\">")
                baseline_shift = -0.6 * 0.83
            end
        elseif mat.captures[2] == "sub"
            if mat.captures[1] == "/"
                write(output, "</tspan>")
            else
                # write(output, "<tspan style=\"dominant-baseline:inherit\" baseline-shift=\"sub\">")
                write(output, "<tspan style=\"dominant-baseline:inherit\" dy=\"0.6em\" font-size=\"83%\">")
                baseline_shift = 0.6 * 0.83
            end
        elseif mat.captures[2] == "i"
            if mat.captures[1] == "/"
                write(output, "</tspan>")
            else
                write(output, "<tspan style=\"dominant-baseline:inherit\" font-style=\"italic\">")
            end
        elseif mat.captures[2] == "b"
            if mat.captures[1] == "/"
                write(output, "</tspan>")
            else
                write(output, "<tspan style=\"dominant-baseline:inherit\" font-weight=\"bold\">")
            end
        end

        if closing_tag && baseline_shift != 0.0
            @printf(output, "<tspan dy=\"%fem\">", -baseline_shift)
            baseline_shift = 0.0
            open_tag = true
        end

        lastpos = mat.offset + length(mat.match)
    end
    write(output, input[lastpos:end])
    if open_tag
        write(output, "</tspan>")
    end
    String(take!(output))
end


text:
# Basic Model Interface (BMI) implementation based on
# https://github.com/Deltares/BasicModelInterface.jl

# At the current state this should be considered as proof of concept.
# To actually use it, it should not only be cleaned up,
# but also made general enough to support different scenarios.

"""
    BMI.initialize(::Type{<:Athesis.Simulation}, config_file)

Initialize the simulation.
Will return a Simulation that is ready to run.
"""
function BMI.initialize(::Type{<:Simulation}, config_file)
    myFloat = Float64
    input = getDefaultInput(myFloat, config_file)
    useGPU = false
    simulation = initSimulation(input, useGPU, myFloat)
    return simulation
end

"""
    BMI.update(simulation::Simulation)

Advance model state by one time step.

Perform all tasks that take place within one pass through the model's
time loop.
"""
function BMI.update(simulation::Simulation, dt)
    time = simulation.timeData.time + dt
    @info("Δt is $dt")
    @info("Updating until $time")
    while simulation.timeData.time < time
        doTimeStep!(simulation)
    end
end

"""
    BMI.update_until(simulation::Simulation)

Advance model state until the given time.

The given `time` must be a model time later than the current model time.
"""
function BMI.update_until(simulation::Simulation, time)
    while simulation.timeData.time < time
        doTimeStep!(simulation)
    end
end

"""
    BMI.get_start_time(simulation::Simulation)

Start time of the simulation.
"""
function BMI.get_start_time(simulation::Simulation)
    return 0.0
end

"""
    BMI.get_end_time(simulation::Simulation)

End time of the simulation.
"""
function BMI.get_end_time(simulation::Simulation)
    return simulation.timeData.tend
end

"""
    BMI.get_current_time(simulation::Simulation)

Current time of the simulation.
"""
function BMI.get_current_time(simulation::Simulation)
    return simulation.timeData.time
end

"""
    BMI.get_time_step(simulation::Simulation)

Time step of the simulation.
"""
function BMI.get_time_step(simulation::Simulation)
    return simulation.timeData.Δt
end

function BMI.finalize(simulation::Simulation)
end


text:
"""
credit account (A/c)
"""

"""General Accounting formula

Assets = Liabilities + Capital

(long-term + short-term )Assets = accounts payable (creditor) + losses_incurred + fees(one_time + recurring)
"""

"""
Credit Accounts

Assets:
-Short-term (within 1 year, or less)
-Cash
-Accounts Recievable

-long-term assets (some books specify 5 years )
(on average: a life-time of 3 years)
-House
-Car
-PC

----------
Liabilities:

-Adding a creditor(accounts payable): i.e. pay you later  A/c
creditoraccount is opened, if existing, the amount is accumulated

----------
Capital:

-Bank account
#######
-shares (needs planning )  <-------------------------------------#TODO #issues: gains, loss handling, share splits (number is doubled, while value is halved )
--stock <-----#TODO #issues: gains, loss handling, share splits (number is doubled, while value is halved )
--ETFs
--Bonds <---what type of bonds , class?AA, AAA, return, payment cycle (default: bi-annual )
#######

-Recieving gains: #Paying off a service
--One time: gifts, Inheritance,
--recurring: Alimony,
---other recurring: passive income , youtube earnings ,

-Revenue gain(s) (realized) #from (running small) busineess

----------

"""

struct crAccount <: Account

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# 2d 2scale mesh
# element easily described by x and y coordinates
# element boundary easily described by a vector of boundary location (ordering: x increases first y second) [write functions]
# oversampling element easily specified by x and y coordinates; oversampling boundary obtained by a vector of locations [write functions]
struct ExpMsFEM_2d2ScaleUnifQuadMesh{Ti,Tf}
    Nce::Ti # number of coarse elements in each dimension
    Nfe::Ti # number of fine elements in each coarse element
    Ne::Ti # number of total fine elements in each dimension
    CGrid_x::Vector{Tf} # coarse uniform grid in x axis, boundary included
    CGrid_y::Vector{Tf} # coarse uniform grid in y axis, boundary included
    ElemNode_loc2glo::Function # local index (node of an element) to global index (global indexed node); one runs over i first then j
    # all but the first attribute are automatically generated
end

function ExponentialMsFEM_GlobalAssembly(MsFEMparam, PDEparam, l)
    # N_c num of coarse intervals each dimension
    # N_f num of fine intervals in each coarse interval
    # l: num of edge basis functions per edge
    N_c = MsFEMparam.Nc
    N_f = MsFEMparam.Nf
    coarse_x = MsFEMparam.coarse_x
    coarse_y = MsFEMparam.coarse_y

    Nfinc = N_f+1;
    N_basis = l+1;
    nNodes = (N_c+1)*(N_c+1)+2*N_c*(N_c+1)*N_basis;

    # sparse assembling
    numb = (4+4*N_basis)^2*(N_c)^2;
    I = zeros(numb)
    J = copy(I)
    K = copy(I)
    F = zeros(nNodes);
    count = 4 + 4*N_basis;
    val = zeros(Nfinc^2,count,N_c^2);

    # there might be inner patches, patches on the edge, and nodal patches
    for i = 1:N_c
        for j = 1:N_c
            [value,k, f] = ExponentialMsFEM_LocalAssembly(PDEparam, coarse_x, coarse_y, i, j,N_f, N_basis);
            val[:,:,(i-1)*N_c+j] = value;
            for p = 1:count
                if p<5 # nodal basis
                    global_p = PDEparam.loc2glo(N_c, i, j, p);
                elseif p<5+N_basis # first edge
                    global_p = (i-1+(j-1)*N_c)*N_basis+(N_c+1)^2+p-4;
                elseif p<5+2*N_basis # second edge
                    global_p = (i-1+(j)*N_c)*N_basis+(N_c+1)^2+p-4-N_basis;
                elseif p<5+3*N_basis # third edge
                    global_p = (i-1+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+p-4-2*N_basis;
                else # last edge
                    global_p = (i+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+p-4-3*N_basis;
                end
                for q = 1:count
                    if q<5
                        global_q = MsFEMparam.loc2glo(N_c, i, j, q);
                    elseif q<5+N_basis
                        global_q = (i-1+(j-1)*N_c)*N_basis+(N_c+1)^2+q-4;
                    elseif q<5+2*N_basis
                        global_q = (i-1+(j)*N_c)*N_basis+(N_c+1)^2+q-4-N_basis;
                    elseif q<5+3*N_basis
                        global_q = (i-1+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+q-4-2*N_basis;
                    else
                        global_q = (i+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+q-4-3*N_basis;
                    end
                    index = (4+4*N_basis)^2*(N_c)*(i-1)+(4+4*N_basis)^2*(j-1)+(4+4*N_basis)*(p-1)+q;
                    I[index] = global_p;
                    J[index] = global_q;
                    K[index] = k[p, q];
                end
                F[global_p] = F[global_p] + f[p];
            end
        end
    end

    # location
    b = zeros(4*(N_c)+4*N_c*N_basis);
    b[1:4*(N_c)+2*N_c*N_basis]=reduce(vcat,collect.([
        1:N_c+1,N_c+2:N_c+1:(N_c+1)*(N_c+1),2*(N_c+1):N_c+1:(N_c+1)*(N_c+1),
        N_c*N_c+N_c+2:(N_c+1)*(N_c+1)-1,(N_c+1)^2+1:(N_c+1)^2+N_basis*N_c,
        (N_c+1)^2+1+N_c^2*N_basis:(N_c+1)^2+N_basis*N_c+N_c^2*N_basis
        ]));

    for i=1:N_basis
        b[4*N_c+2*N_c*N_basis+1+(i-1)*2*N_c:4*N_c+2*N_c*N_basis+(i)*2*N_c]=
        reduce(vcat,collect.([i+N_c*(N_c+1)*N_basis+(N_c+1)^2:N_basis*(N_c+1):i+N_basis*(N_c+1)*(N_c-1)+N_c*(N_c+1)*N_basis+(N_c+1)^2,
            i+N_c*N_basis+N_c*(N_c+1)*N_basis+(N_c+1)^2:N_basis*(N_c+1):i+N_c*N_basis+N_basis*(N_c+1)*(N_c-1)+N_c*(N_c+1)*N_basis+(N_c+1)^2]));
    end

    # assemby
    A = sparse(I,J,K,nNodes,nNodes);
    A[b,:] .= 0; A[:,b] .= 0; F[b] .= 0;
    A[b,b] .= sparse(LinearAlgebra.I, length(b),length(b));

    return A, F
end

function ExponentialMsFEM_LocalAssembly(PDEparam, X, Y, m, n, N_f, N_basis)
    N_c = length(X)-1;

    xlow = coarse_x[i];
    xhigh = coarse_x[i+1];
    ylow = coarse_y[j];
    yhigh = coarse_y[j+1];
    x = collect(LinRange(xlow, xhigh, N_f+1))
    y = collect(LinRange(ylow, yhigh, N_f+1))

    # location
    b=reduce(vcat,collect.([1:N_f+1,N_f+2:N_f+1:(N_f+1)*(N_f+1),2*(N_f+1):N_f+1:(N_f+1)*(N_f+1),N_f*N_f+N_f+2:(N_f+1)*(N_f+1)-1]))

    # edge space (matrix)
    f=collect.([LinRange(1, 0, N_f+1), LinRange(1-1/N_f, 0, N_f),zeros(1,N_f),zeros(1,N_f-1);
        LinRange(0, 1, N_f+1),zeros(1,N_f),LinRange(1-1/N_f, 0, N_f),zeros(1,N_f-1);
        zeros(1,N_f+1),zeros(1,N_f),LinRange(1/N_f, 1, N_f),LinRange(1/N_f, 1-1/N_f, N_f-1);
        zeros(1,N_f+1),LinRange(1/N_f, 1, N_f),zeros(1,N_f),LinRange(1-1/N_f, 1/N_f, N_f-1)]); # need to update

    A = basefun(MsFEMparam, PDEparam, X, Y, m, n, N_f); # harmonic extension matrix
    B = A;
    F = -A[:,b] * f';
    A[b,:] .= 0; A[:,b] .= 0; F[b,:] .= f';
    A[b,b] .= sparse(I, length(b),length(b));
    value = A\F; # get all harmonic extended solutions

    # edge basis
    count = 4;
    if n>1 # (m,n) location
        [L1,L2,N] = harmext(MsFEMparam, PDEparam, X, Y, m, n-1, N_f,1);
        [R,P,bub] = restrict(X, Y, m, n-1,N_f,1);
        [V,D] = eigs(R'*N*R,P,N_basis);
        value[:,count+1:count+N_basis] = L2*R*V;
        value[:,count+N_basis+1] = L2*bub;
    else
        value[:,count+1:count+N_basis+1] = 0;
    end
    count += N_basis +1;
    if n<N_c
        [L1,L2,N] = harmext(MsFEMparam, PDEparam, X, Y, m, n, N_f, 1);
        [R,P,bub] = restrict(X, Y, m, n, N_f, 1);
        [V,D] = eigs(R'*N*R,P,N_basis);
        value[:,count+1:count+N_basis] = L1*R*V;
        value[:,count+N_basis+1] = L1*bub;
    else
        value[:,count+1:count+N_basis+1] .= 0;
    end

    count += N_basis + 1;
    if m>1
        [L1,L2,N] = harmext(MsFEMparam, PDEparam, X, Y, m-1, n, N_f, 2);
        [R,P,bub] = restrict(X, Y, m-1, n, N_f, 2);
        [V,D] = eigs(R'*N*R,P,N_basis);
        value[:,count+1:count+N_basis] = L2*R*V;
        value[:,count+N_basis+1] = L2*bub;
    else
        value[:,count+1:count+N_basis+1] .= 0;
    end
    count += N_basis + 1;
    if m<N_c
        [L1,L2,N] = harmext(MsFEMparam, PDEparam, X, Y, m, n, N_f, 2);
        [R,P,bub] = restrict(X, Y, m, n, N_f, 2);
        [V,D] = eigs(R'*N*R,P,N_basis);
        value[:,count+1:count+N_basis] = L1*R*V;
        value[:,count+N_basis+1] = L1*bub;
    else
        value[:,count+1:count+N_basis+1] .= 0;
    end
    count += N_basis + 1;

    K=value'*B*value;
    # k and f separate using local stiffness matrix%%%%%%

    F = zeros(N_f^2,count);
    for u = 1:N_f
        for v = 1:N_f
            for x1=0:1
                for y1=0:1
                    xf=(x(u+1)+x(u))/2;
                    yf=(y(v+1)+y(v))/2;
                    F[v*N_f+u-N_f,:] = F[v*N_f+u-N_f,:] + PDEparam.rhs(xf,yf)*value[(v+y1-1)*(N_f+1)+u+x1,:]*(x(2)-x(1))^2/4;
                end
            end
        end
    end
    f = sum(F);
    return value, K, f
end

function ExponentialMsFEM_Solver(MsFEMparam,PDEparam,l)

    A, F = ExponentialMsFEM_GlobalAssembly(MsFEMparam, PDEparam, l)
    @info "[Assembly] finish assembly of stiffness matrice"
    u = A\F;

    @info "[Linear solver] linear system solved"
    result = zeros((N_c*N_f+1));

    for i = 1:N_c
        for j = 1:N_c

            xlow = coarse_x[i];
            xhigh = coarse_x[i+1];
            ylow = coarse_y[j];
            yhigh = coarse_y[j+1];

            xs = LinRange(xlow, xhigh, Nfinc);
            ys = LinRange(ylow, yhigh, Nfinc);
            zs = zeros(Nfinc, Nfinc);
            value = val[:,:,(i-1)*N_c+j];
            for p = 1:count
                if p<5
                    global_p = MsFEMparam.loc2glo(N_c, i, j, p);
                elseif p<5+N_basis
                    global_p = (i-1+(j-1)*N_c)*N_basis+(N_c+1)^2+p-4;
                elseif p<5+2*N_basis
                    global_p = (i-1+(j)*N_c)*N_basis+(N_c+1)^2+p-4-N_basis;
                elseif p<5+3*N_basis
                    global_p = (i-1+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+p-4-2*N_basis;
                else
                    global_p = (i+(j-1)*(N_c+1)+N_c*(N_c+1))*N_basis+(N_c+1)^2+p-4-3*N_basis;
                end
                nodevalue = u[global_p];
                valuep = reshape(value[:,p],nSamples,nSamples);
                zs = zs + nodevalue* valuep;
            end
            # bubble part
            zs = zs + bubble(x,y,i,j,N_f);
            result[(i-1)*N_f+1:(i-1)*N_f+nSamples,(j-1)*N_f+1:(j-1)*N_f+Nfinc] = zs;
            # surf(xs, ys, zs'); # for plot the solution
        end
    end
    # result=reshape(result,(N_c*N_f+1)^2,1);
    return result
end

function basefun(MsFEMparam, PDEparam, X, Y, m, n, N_f)
    # % solve a combination of fine basis functions, assembling the stiffness
    # % matrix

    # % the function serves as to solve local boundary value problem on the
    # % coarse mesh, based on the stiffness matrix in "elementstiff1"

    # % for the sake of simplicity, we temporarily use nodes of the fine mesh
    # % as the same nodes that we use when we invoke numerical quadrature to
    # % assemble the global stiffness matrix, the rationale being that we
    # % can either increase the scale of fine mesh or increase the quadrature
    # % accuracy if we want to do higher order.... (because locally on fine
    # % meshes, the basis function is simply a linear function

    # % setting parameters H is the size of the coarse mesh, N is the number of
    # %nodes on each coarse rectangular, i.e. fine mesh size is H/N

    # %N needs to be inputted !!!!

    x = collect(LinRange(X(m), X(m+1), N_f+1));
    y = collect(LinRange(Y(n), Y(n+1), N_f+1));
    nNodes = (N_f+1)*(N_f+1);

    # sparse assembling
    I = zeros(16*N_f^2,1);
    J = copy(I)
    K = copy(I)

    for i = 1:N_f
        for j = 1:N_f
            [k] = elementstiff1(PDEparam, x, y, i, j);
            for p = 1:4
                global_p = MsFEMparam.loc2glo(N_f, i, j, p);
                for q = 1:4
                    index=16*N_f*(i-1)+16*(j-1)+4*(p-1)+q;
                    global_q = MsFEMparam.loc2glo(N_f, i, j, q);
                    I[index] = global_p;
                    J[index] = global_q;
                    K[index] = k[p, q];
                end
            end
        end
    end
    # assemble boundary
    A = sparse(I,J,K,nNodes,nNodes);
    return A
end

function elementstiff1(PDEparam, X, Y, m, n)
    # %the stiffness matrix required locally to compute local nodal basis
    # %for the bilinear boundary value with 1 at node "node"

    # %X Y are the fine mesh, we shall use local bilinear basis on the
    # %fine mesh as basis to compute the harmonic problem for bilinear
    # %boundary value problem on the coarse mesh

    # %H is the size of the coarse mesh

    # %H, nSamples and coefficients a needs to be inputted  !!!!!!!!!

    # %nSamples = 3; % sample points in trapz in each direction
    K = zeros(4, 4);
    xlow = X[m];
    xhigh = X[m+1];
    ylow = Y[n];
    yhigh = Y[n+1];

    x = (xlow+xhigh)/2;
    y = (ylow+yhigh)/2;

    for i = 1:4
        for j = 1:4
            if i==j
                K[i, j] = 2/3*PDEparam.a(x,y);
            elseif i==j+2 || i==j-2
                    K[i,j] = -1/3*PDEparam.a(x,y);
            else
                K[i,j] = -1/6*PDEparam.a(x,y);
            end
        end
    end
    return K
end

function harmext(MsFEMparam, PDEparam, X, Y, m, n, N_f, i)
    # %compute the local harmonic extension corresponding to each coarse edge
    # %i=1, corresponds to the horizontal edges, i=2, corresponds to the vertical
    # %ones. m, n are indices

    # %use dirichlet solver for adajacent patches, obtaining L1.L2 as a linear
    # %combination of the corresponding patches, N is matrix of the inner product

    b = reduce(vcat,collect.([1:N_f+1,N_f+2:N_f+1:(N_f+1)*(N_f+1),2*(N_f+1):N_f+1:(N_f+1)*(N_f+1),N_f*N_f+N_f+2:(N_f+1)*(N_f+1)-1]));
    K1 = basefun(MsFEMparam, PDEparam, X, Y, m, n, N_f);
    f1= zeros(4*N_f,N_f-1);
    f2= zeros(4*N_f,N_f-1);
    if i==1
        K2 = basefun(MsFEMparam, PDEparam, X, Y, m, n+1, N_f);
        f1[3*N_f+2:4*N_f,:] .= sparse(LinearAlgebra/I,N_f-1,N_f-1)
        f2[2:N_f,:] = sparse(LinearAlgebra.I,N_f-1,N_f-1)
    else
        K2 = basefun(MsFEMparam, PDEparam, X, Y, m+1, n, N_f);
        f1[2*N_f+2:3*N_f,:] .= sparse(LinearAlgebra.I,N_f-1,N_f-1)
        f2[N_f+2:2*N_f,:] .= sparse(LinearAlgebra.I,N_f-1,N_f-1)
    end
    M1 = K1;
    M2 = K2;
    F1 = -K1[:,b]*f1;
    K1[b,:] .=0; K1[:,b] .= 0; F1[b,:] = f1;
    K1[b,b] = sparse(LinearAlgebra.I,length(b),length(b))
    L1 = K1\F1;
    F2 = -K2[:,b]*f2;
    K2[b,:] .= 0; K2[:,b] .= 0; F2[b,:] = f2;
    K2[b,b] .= sparse(LinearAlgebra.I,length(b),length(b))
    L2=K2\F2;

    N1 = L1'*M1*L1;
    N2 = L2'*M2*L2;
    N = N1 + N2;
    return L1, L2, N
end

function basefun1(MsFEMparam, PDEparam, X, Y, m, n, N_f,t)
    # % solve a combination of fine basis functions, assembling the oversampled
    # % stiffness matrix
    # %t=1, corresponds to the horizontal edges, t=2, corresponds to the vertical
    # %ones. m, n are indices

    # % the function serves as to solve local boundary value problem on the
    # % coarse mesh, based on the stiffness matrix in "elementstiff1"

    # % for the sake of simplicity, we temporarily use nodes of the fine mesh
    # % as the same nodes that we use when we invoke numerical quadrature to
    # % assemble the global stiffness matrix, the rationale being that we
    # % can either increase the scale of fine mesh or increase the quadrature
    # % accuracy if we want to do higher order.... (because locally on fine
    # % meshes, the basis function is simply a linear function

    # % setting parameters H is the size of the coarse mesh, N is the number of
    # %nodes on each coarse rectangular, i.e. fine mesh size is H/N

    # %N needs to be inputted !!!!
    N_c = length(X)-1;
    if t==1
        if m==1
            x = linspace(X(m), X(m+2), 2*N_f+1);
        elseif m==N_c
            x = linspace(X(m-1), X(m+1), 2*N_f+1);
        else
            x = linspace(X(m-1), X(m+2), 3*N_f+1);
        end
        y = linspace(Y(n), Y(n+2), 2*N_f+1);
    else
        if n==1
            y = linspace(Y(n), Y(n+2), 2*N_f+1);
        elseif n==N_c
            y = linspace(Y(n-1), Y(n+1), 2*N_f+1);
        else
            y = linspace(Y(n-1), Y(n+2), 3*N_f+1);
        end
        x = linspace(X(m), X(m+2), 2*N_f+1);
    end

    nNodes = length(x)*length(y);
    N_x=length(x)-1;
    N_y=length(y)-1;
    # sparse assembling
    I=zeros(16*N_x*N_y,1);J=I;K=I;

    for i = 1:N_x
        for j = 1:N_y
            [k] = elementstiff1(PDEparam, x, y, i, j);
            for p = 1:4
                global_p = MsFEMparam.loc2glo(N_x, i, j, p);
                for q = 1:4
                    index=16*N_y*(i-1)+16*(j-1)+4*(p-1)+q;
                    global_q = MsFEMparam.loc2glo(N_x, i, j, q);
                    I[index] = global_p;
                    J[index] = global_q;
                    K[index] = k[p, q];
                end
            end
        end
    end

    # assemble boundary

    A=sparse(I,J,K,nNodes,nNodes);

    return A, N_x, N_y, x, y
end

function restrict(X, Y, m, n, N_f, t)
    # %compute the matrix corresponding to the oversampled patch and its
    # %restriction operator on the edge, based on the stiffness matrix assembled
    # %in basefun1, and use local subproblem to solve for the bubble and harmonic
    # %parts.

    # %P is the inner product corresponding to the oversampled norm, R is the
    # %restriction operator
    fine = (X(2)-X(1))/N_f;
    N_c = length(X)-1;

    [A,N_x,N_y,x,y]=basefun1(MsFEMparam, PDEparam, X, Y, m, n,N_f,t);
    b=[1:N_x+1,N_x+2:N_x+1:(N_x+1)*(N_y+1),2*(N_x+1):N_x+1:(N_x+1)*(N_y+1),N_x*N_y+N_y+2:(N_x+1)*(N_y+1)-1];
    # %assemble boundary b

    # %solve the basis function corresponding to the harmonic part, whose dof
    # %corresponds to the boundary values

    # %pay attention to the case of boundary patches
    if t==1 && n==1 && 2<m<N_c-1 || t==2 && n<3 && 1<m<N_c-1
        f = zeros(2*(N_x+N_y),N_x+2*N_y-1);
        f[N_x+2:2*(N_x+N_y),:] .= sparse(LinearAlgebra.I, N_x+2*N_y-1, N_x+2*N_y-1)
    elseif t==1 && n==N_c-1 && 2<m<N_c-1 || t==2 && n>N_c-2 && 1<m<N_c-1
        f = zeros(2*(N_x+N_y),N_x+2*N_y-1);
        f[1:N_x+N_y,1:N_x+N_y] = sparse(LinearAlgebra.I, N_x+N_y, N_x+N_y);
        f[N_x+N_y+2:N_x+2*N_y,N_x+N_y+1:N_x+2*N_y-1] .= sparse(LinearAlgebra.I, N_y-1, N_y-1);
    elseif t==1 && m<3 && 1<n<N_c-1 || t==2 && m==1 && 2<n<N_c-1
        f = zeros(2*(N_x+N_y),2*N_x+N_y-1);
        f[2:N_x+1,1:N_x] .= sparse(LinearAlgebra.I, N_x, N_x);
        f[N_x+N_y+2:2*N_x+2*N_y,N_x+1:2*N_x+N_y-1] .= sparse(LinearAlgebra.I, N_x+N_y-1,N_x+N_y-1);
    elseif t==1 && m>N_c-2 && 1<n<N_c-1 || t==2 && m==N_c-1 && 2<n<N_c-1
        f = zeros(2*(N_x+N_y),2*N_x+N_y-1);
        f[1:N_x,1:N_x] .= sparse(LinearAlgebra.I, N_x, N_x);
        f[N_x+2:N_x+N_y+1,N_x+1:N_x+N_y] .= sparse(LinearAlgebra.I, N_y, N_y);
        f[N_x+2*N_y+2:2*(N_x+N_y),N_x+N_y+1:2*N_x+N_y-1] .= sparse(LinearAlgebra.I, N_x-1, N_x-1);
    elseif t==1 && n==1 && m<3 || t==2 && n<3 && m==1
        f = zeros(2*(N_x+N_y),N_x+N_y-1);
        f[N_x+N_y+2:2*(N_x+N_y),:] .= sparse(LinearAlgebra.I, N_x+N_y-1, N_x+N_y-1);
    elseif t==1 && n==1 && m>N_c-2 || t==2 && n<3 && m==N_c-1
        f = zeros(2*(N_x+N_y),N_x+N_y-1);
        f[N_x+2:N_x+N_y+1,1:N_y] .= sparse(LinearAlgebra.I, N_y, N_y);
        f[N_x+2*N_y+2:2*N_x+2*N_y,N_y+1:N_x+N_y-1] .= sparse(LinearAlgebra.I, N_x-1, N_x-1);
    elseif t==1 && n==N_c-1 && m<3 || t==2 && n>N_c-2 && m==1
        f = zeros(2*(N_x+N_y),N_x+N_y-1);
        f[2:N_x+1,1:N_x] .= sparse(LinearAlgebra.I, N_x, N_x);
        f[N_x+N_y+2:N_x+2*N_y,N_x+1:N_x+N_y-1] .= sparse(LinearAlgebra.I, N_y-1, N_y-1);
    elseif t==1 && n==N_c-1 && m>N_c-2 || t==2 && n>N_c-2 && m==N_c-1
        f = zeros(2*(N_x+N_y),N_x+N_y-1);
        f[1:N_x,1:N_x] .= sparse(LinearAlgebra.I, N_x, N_x);
        f[N_x+2:N_x+N_y,N_x+1:N_x+N_y-1] .= sparse(LinearAlgebra.I, N_y-1, N_y-1);
    else
        f = zeros(2*(N_x+N_y),2*(N_x+N_y)-1);
        f[1:2*(N_x+N_y)-1,:] .= sparse(LinearAlgebra.I, 2*(N_x+N_y)-1, 2*(N_x+N_y)-1);
    end

    B = A;
    F = -A(:,b)*f;
    B[b,:] .= 0; B[:,b] .= 0; F[b,:] = f;
    B[b,b] = sparse(LinearAlgebra.I, length(b),length(b))
    harm = B\F;

    # %solve the basis function corresponding to the bubble part, whose dof
    # %corresponds to the right hand side values

    # %assemble rhs
    G = zeros((N_x+1)*(N_y+1), 1);
    for i = 1:N_x
        for j = 1:N_y
            for p=1:4
                xf=(x(i+1)+x(i))/2;
                yf=(y(j+1)+y(j))/2;
                G[MsFEMparam.loc2glo(N_x, i, j, p)]= G[MsFEMparam.loc2glo(N_x, i, j, p)] + PDEparam.rhs(xf,yf)*fine^2/4;
            end
        end
    end
    G[b] .= 0;
    bub = B\G;

    # %compute inner product of energy
    # %P=fine^2*speyes(N_x*N_y+2*(N_x+N_y));
    # %P(2*(N_x+N_y),2*(N_x+N_y))=harm'*A*harm;
    P=harm'*A*harm;
    # %indentifying the edge
        leng=N_f+1;
    if t==1
        if m==1
            c=[(N_x+1)*N_f+1:(N_x+1)*N_f+N_f+1];
        else
            c=[(N_x+1)*N_f+N_f+1:(N_x+1)*N_f+2*N_f+1];
        end
    else
        if n==1
            c=[N_f+1:N_x+1:(N_x+1)*N_f+N_f+1];
        else
            c=[(N_x+1)*N_f+N_f+1:N_x+1:(N_x+1)*N_f*2+N_f+1];
        end
    end

    # %intepolate
    R = harm[c,:];
    R = R-linspace(1,0,leng)'*R(1,:)-linspace(0,1,leng)'*R(leng,:);
    R = R[2:leng-1,:]
    bub = bub[c];
    bub = bub-linspace(1,0,leng)'*bub(1)-linspace(0,1,leng)'*bub(leng);
    bub = bub[2:leng-1];
    return R, P, bub
end

function bubble(X, Y, m, n,N_f)
    # %solve the bubble part
    x = linspace(X(m), X(m+1), N_f+1);
    y = linspace(Y(n), Y(n+1), N_f+1);
    nNodes = (N_f+1)*(N_f+1);
    # %sparse assembling
    I = zeros(16*N_f^2,1);J=I;K=I;
    F = zeros(nNodes, 1);

    for i = 1:N_f
        for j = 1:N_f
            [k,f] = elementstiff2(x, y, i, j);
            for p = 1:4
                global_p = loc2glo(N_f, i, j, p);
                for q = 1:4
                    index=16*N_f*(i-1)+16*(j-1)+4*(p-1)+q;
                    global_q = loc2glo(N_f, i, j, q);
                    I[index] = global_p;
                    J[index] = global_q;
                    K[index] = k(p, q);
                end
                F[global_p] = F[global_p] + f[p];
            end
        end
    end

    # %assemble boundary

    A=sparse(I,J,K,nNodes,nNodes);
    # %assemble boundary
    b=[1:N_f+1,N_f+2:N_f+1:(N_f+1)*(N_f+1),2*(N_f+1):N_f+1:(N_f+1)*(N_f+1),N_f*N_f+N_f+2:(N_f+1)*(N_f+1)-1,];

    A=sparse(I,J,K,nNodes,nNodes);

    A[b,:] .= 0; A[:,b] .= 0; F[b] .= 0;
    A[b,b] = sparse(LinearAlgebra.I,length(b),length(b))
    u = A\F;
    re = reshape(u,N_f+1,N_f+1);

    return re
end


text:
push!(LOAD_PATH, string(pwd(), "/src/"))
using WRIdual
using LinearAlgebra, SparseArrays
using Optim, LineSearches
using PyPlot
using JLD2
using Dates

## Load slowness model and data

@load "./data/Gausslens/Gausslens_data.jld"

## Defining objective functional

# Background model
mb = 1e0./2000e0.^2e0.*ones(Float64, size(m))

# Setting model preprocessing and gradient preconditioning
mask = BitArray(undef, size(m)); mask .= false
mask[3:end-2, 3:end-2] .= true

# Source encoder
src_enc = C.(Matrix(sparse(I, length(x_src), length(x_src))))

# Pre-/post-processing for objective functionals
preproc(x) = contr2abs(x, mask, mb)
postproc(g) = gradprec_contr2abs(g, mask, mb)

# WRI dual functional setup
ε = 0.01
δ = 1e-1*compute_srcfcsδ(mb[1, 1], freqs[1])
weight_array = src_focusing_array(size(mb), nabc, h, x_src, z_src; δ = δ)
fun!(F, G, x) = objWRIdual_ε_yres!(F, G, preproc(x), h, [freqs[1]], x_rcv, z_rcv, x_src, z_src, src_enc, [dat[1]], ε; nabc = nabc, weight_array = weight_array, gradmprec_fun = postproc)

# # FWI functional setup
# fun!(F, G, x) = objFWI!(F, G, preproc(x), h, [freqs[1]], x_rcv, z_rcv, x_src, z_src, src_enc, [dat[1]]; nabc = nabc, gradprec_fun = postproc)

## Optimization parameters

# Starting guess
x0 = zeros(Float64, length(findall(mask)))

# Options
# method = LBFGS(; alphaguess = InitialStatic(; alpha = 1f4), linesearch = BackTracking())
method = LBFGS()
niter = 20
opt = Optim.Options(iterations = niter, store_trace = true, show_trace  = true, show_every  = 1)

## Run optimization

# Minimization
result = optimize(Optim.only_fg!(fun!), x0, method, opt)
x0 = Optim.minimizer(result)
minv = preproc(x0)
loss_log = Optim.f_trace(result)

# Plot result
imshow(minv', aspect = "auto", cmap = "jet")


text:
using TMLEEpistasis
using ArgParse

function parse_commandline()
    s = ArgParseSettings(description="Generate a summary file of the full pipeline")

    @add_arg_table s begin
        "prefix"
            arg_type = String
            help = "Prefix to both estimates/sieve .hdf5 files"
            required = true
        "out"
            arg_type = String
            help = "Path where the summary file will be saved"
            required = true
        "--sieve"
            arg_type = Bool
            help = "If a sieve file is actually present in the hdf5 files"
            action = :store_true
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

build_summary(parsed_args)


text:
## @lc code=start
using LeetCode

function minimum_removal(beans::Vector{Int})::Int
    beans = sort(beans)
    n, same = length(beans), 1
    min_cost = cost = sum(beans) - n * first(beans)
    for k in 2 : n
        if beans[k] > beans[k - 1]
            cost = cost + same * beans[k - 1] - (beans[k] - beans[k-1]) * (n - k + 1)
            min_cost, same = min(cost, min_cost) , 1
        else
            same += 1
        end
    end
    min_cost
end
minimum_removal_oneline(beans::Vector{Int}) = sum(beans) - maximum((length(beans) - i + 1) * v for (i, v) in enumerate(sort(beans)))
## @lc code=end

## @lc test=start
@testset "2171.removing-minimum-number-of-magic-beans.jl" begin
    @test minimum_removal([4,1,6,5]) == 4
    @test minimum_removal([2,10,3,2]) == 7
    @test minimum_removal([25,27,1,10,8,35,17,5,4,16]) == 68
    @test minimum_removal_oneline([4,1,6,5]) == 4
    @test minimum_removal_oneline([2,10,3,2]) == 7
    @test minimum_removal_oneline([25,27,1,10,8,35,17,5,4,16]) == 68
end
## @lc test=end

## @lc info=start
# link: [solution to 2171](https://leetcode-cn.com/problems/removing-minimum-number-of-magic-beans/solution/6006-na-chu-zui-shao-shu-mu-de-mo-fa-dou-bsey/)
## @lc info=end