text:
using Keldysh, Test

@testset "time_grid" begin
  let c = FullContour(tmax=2.0, β=5.0), nt = 21, ntau=51
    grid = FullTimeGrid(c, nt, ntau)
    @test_throws MethodError KeldyshTimeGrid(c, nt)
    @test_throws MethodError ImaginaryTimeGrid(c, ntau)

    @test length(grid) == 2*nt + ntau

    x = [grid[i] for i in 1:length(grid)]
    @test length(x) == 2*nt + ntau

    @test length(grid[forward_branch]) == nt
    @test length(grid[backward_branch]) == nt
    @test length(grid[imaginary_branch]) == ntau

    @test step(grid, forward_branch) == 0.1
    @test step(grid, backward_branch) == -0.1
    @test step(grid, imaginary_branch) == -0.1im

    @test integrate(t -> 1, grid) ≈ -1.0im * c.β

    for b in c.branches
      for t in b.(range(0.0, 1.0, length=3))
        tl = Keldysh.find_lower(grid, t)
        @test tl.bpoint.domain == t.domain
        @test tl.cidx != length(grid)
        tu = grid[tl.cidx + 1]

        @test heaviside(c, t, tl.bpoint)
        @test heaviside(c, tu.bpoint, t)
        @test heaviside(c, t, grid[1].bpoint)
        @test heaviside(c, grid[end].bpoint, t)
      end
    end
  end

  let c = KeldyshContour(tmax=2.0), nt = 21, ntau=51
    grid = KeldyshTimeGrid(c, nt)
    @test_throws MethodError FullTimeGrid(c, nt, ntau)
    @test_throws MethodError ImaginaryTimeGrid(c, ntau)

    @test length(grid) == 2*nt

    x = [grid[i] for i in 1:length(grid)]
    @test length(x) == 2*nt

    @test length(grid[forward_branch]) == nt
    @test length(grid[backward_branch]) == nt
    @test_throws AssertionError grid[imaginary_branch]

    @test step(grid, forward_branch) == 0.1
    @test step(grid, backward_branch) == -0.1
    @test_throws AssertionError step(grid, imaginary_branch)

    N = length(grid)
    A = zeros(ComplexF64, N, N)
    B = zeros(ComplexF64, N, N)
    for t1 in grid
      for t2 in grid
        A[t1.cidx, t2.cidx] = integrate(t -> 1.0, grid, t1, t2)
        B[t1.cidx, t2.cidx] = t1.bpoint.val - t2.bpoint.val
      end
    end
    @test maximum(abs.(A .- B)) < 1e-12

  end

  let c = ImaginaryContour(β=5.0), nt = 21, ntau=51
    grid = ImaginaryTimeGrid(c, ntau)
    @test_throws MethodError FullTimeGrid(c, nt, ntau)
    @test_throws MethodError KeldyshTimeGrid(c, nt)

    @test length(grid) == ntau

    x = [grid[i] for i in 1:length(grid)]
    @test length(x) == ntau

    @test length(grid[imaginary_branch]) == ntau
    @test_throws AssertionError grid[forward_branch]
    @test_throws AssertionError grid[backward_branch]

    @test_throws AssertionError step(grid, forward_branch)
    @test_throws AssertionError step(grid, backward_branch)
    @test step(grid, imaginary_branch) == -0.1im
  end

end


text:
module Grids
#    using DataTypes

    export StaggeredGrid2D, getExtrapolatedValue

    abstract Grid

    type StaggeredGrid2D <: Grid
        nx::Int
        ny::Int
        nu::Int
        nghost::Int
        np::Int
        dx::Real
        dy::Real
        sizes::Dict

        u_left::Vector{Real}
        u_bottom::Vector{Real}
        u_right::Vector{Real}
        u_top::Vector{Real}
        v_left::Vector{Real}
        v_bottom::Vector{Real}
        v_right::Vector{Real}
        v_top::Vector{Real}

        x::Vector{Real}
        y::Vector{Real}

        configure::Function
        init::Function
        createGrid::Function
        getLengths::Function
        getSizes::Function

        function StaggeredGrid2D()
            this = new()
            thisnghost = 1     # using 1 element for the ghost region in each direction

            this.configure = function(config_options_grid)
                this.nx = config_options_grid["nx"]
                this.ny = config_options_grid["ny"]
                this.x = linspace(0,config_options_grid["length"],this.nx+1)
                this.y = linspace(0,config_options_grid["height"],this.ny+1)
                this.dx = 1/this.nx
                this.dy = 1/this.ny
                this.sizes = Dict()
                this.sizes = ["u" => [this.nx-1,this.ny], "v" => [this.nx,this.ny-1], "p" => [this.nx,this.ny]]
            end #fun config

            this.init = function()
                this.nu = this.sizes["u"][1] * this.sizes["u"][2] + this.sizes["v"][1] * this.sizes["v"][2]
                this.np = this.sizes["p"][1] * this.sizes["p"][2]
                this.u_left     = ones(1,this.ny)
                this.u_bottom   = ones(1,this.nx-1)
                this.u_right    = ones(1,this.ny)
                this.u_top      = ones(1,this.nx-1)
                this.v_left     = ones(1,this.ny-1)
                this.v_bottom   = ones(1,this.nx)
                this.v_right    = ones(1,this.ny-1)
                this.v_top      = ones(1,this.nx)
            end #fun init

            this.getLengths = function()
                return this.nu, this.np
            end #fun getLengths

            this.getSizes = function()
                return this.sizes
            end #fun getSizes
            return this
        end
    end #type

    # linear extrapolation of value inside grid to next point in ghost region
    function getExtrapolatedValue(ubc,u)
        return 2.0*ubc-u
    end #fun

end #module

ext: jl


text:
@testset "Chains" begin
    a = Chain(Int8)
    @test dim(a) == 0
    @test length(a) == 0
    @test keytype(a) == UInt
    @test valtype(a) == Int8

    a = Chain(Int, Int)
    @test keytype(a) == Int
    @test valtype(a) == Int

    b = Chain(Int[],Int[])
    @test keytype(b) == Int
    @test valtype(b) == Int
    push!(b, (1, 2))
    @test first(keys(b)) == 1
    @test b[1] == 2
    push!(b, 1, -1)
    @test b[1] == 1
    @test (b + (1, -1))[1] == 0
    @test b[1] == 1

    c = Chain(Int,Int)
    push!(c, 0=>1)
    @test length(c) == 1
    @test first(keys(c)) == 0
    @test first(values(c)) == 1

    d = Chain(Int,UInt8)
    d+=(1,0x01)
    @test d[1] == 1
    @test (d+(1,0x01))[1] == 2
    @test d[1] == 1

    a+=c
    a+=b
    a+=(2,1)
    @testset "loop thrugh chain" for i in keys(a)
        @test a[i] == 1
    end

    b+=(2,1)
    @test b[2] == 1

    append!(c,b)
    @testset "two chain compare" for (aa,cc) in zip(a,c)
        @test cc == aa
    end

    @testset "sum" for (i,v) in a+c
        @test v == 2
    end

    @test (2+b+2)[1] == 5
    @test (2*b*2)[1] == 4
    @test (1+2*b*3+4)[1] == 11

    @test iszero(simplify(a-c))
end


text:
export Domain, ODEtime, DEsteady, solveSS

"""
    DEmeta

Meta inofmration of differential equation. The family of ['DEmeta'](@ref)
"""
abstract type DEmeta end

"""
    ODEtime(func, u0, p, tspan)

Struct for solveing time-series of differential-equations. The data type of ode function is referenced to DifferentialEquations.jl.

Argument
--------
- `func`: ODE function.
- `u0`: initial values
- `p`: parameters
- `tspan`: time span

References
----------

1. [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/ode_example/)

Reset method
------------

Reset a field of `ODEtime` struct with broadcast method.

Usage
-----
When using the `LabelledArray.jl` to define function and subtypes of `DEmeta`. Use this method to update the values of named arrays without changing the type of the field.

Purpose
-------
This feature is to solve the LabelledArray problem. when using the field vector to define function, the initial values or parameters need to be Named vectors. However, the grid search iterator returns vector which gets error when apply with the funtion of name vector.

!!! compat ModelingToolkit
    When using ['jacobian'](@ref), should not use `LabelledArray` or name tuples for model function.

Arguement
---------
- `u0`: should be vector or types that can be broadcast. The length of `u0` should be same as `ODEtime.u0`
- `key`: field name of the stuct (default: `:u0`).

Example
-------
```julia-repl
julia> using LabelledArrays, DifferentialEqiations
julia> u = LVector(s1=1.0,s2=0.2)
2-element LArray{Float64,1,Array{Float64,1},(:s1, :s2)}:
 :s1 => 1.0
 :s2 => 0.2

julia> de = ODEtime(func=x->x, u0=u, p=1.0, tspan=(0.0,1.0))
ODEtime
  func: #9 (function of type var"#9#10")
  u0: LArray{Float64,1,Array{Float64,1},(:s1, :s2)}
  p: Float64 1.0
  tspan: Tuple{Float64,Float64}
  method: CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType},Rational{Int64},Int64}}

julia> de_new = de([1.3,1.4];key=:u0)
ODEtime
  func: #9 (function of type var"#9#10")
  u0: LArray{Float64,1,Array{Float64,1},(:s1, :s2)}
  p: Float64 1.0
  tspan: Tuple{Float64,Float64}
  method: CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType},Rational{Int64},Int64}}

julia> de_new.u0 # The updated u0 is the struct of LArray
2-element LArray{Float64,1,Array{Float64,1},(:s1, :s2)}:
 :s1 => 1.3
 :s2 => 1.4

julia> de.u0
2-element LArray{Float64,1,Array{Float64,1},(:s1, :s2)}:
 :s1 => 1.0
 :s2 => 0.2
```
"""
@with_kw struct ODEtime <: DEmeta
    func
    u0
    p
    tspan
    method= AutoTsit5(Rosenbrock23()) # For autonomaous detecion for stiff and non-stiff problem
end

function (self::ODEtime)(u0; key=:u0)
    return update_struct_broadcast(self, u0, key)
end

"""
    DEsteady(func, p u0, method)
    (self::DEsteady)(u0; key=:u0)

- Struct for solving steady state of an differential equation model.

- Update steady state meta, and return another DEsteady object.

Argument
--------
- `func`: ODE function.
- `u0`: initial values
- `p`: parameters
- `method`: Method for solving steady-states. (i.e. `DifferentialEqiations.Tsit5()`, `DifferentialEquations.AutoTsit5(Rosenbrock23())`)

Reference
---------
1. [ODE solvers of DifferentialEquations.jl](https://diffeq.sciml.ai/stable/solvers/split_ode_solve/)

Example
-------

```julia-repl
julia> using LabelledArrays, DifferentialEqiations
julia> deS = DEsteady(func=x->x, u0= LVector(s1=1.0,s2=2.0), p=1.0)
DEsteady
  func: #17 (function of type var"#17#18")
  p: Float64 1.0
  u0: LArray{Float64,1,Array{Float64,1},(:s1, :s2)}
  method: Tsit5 Tsit5()

julia> deS_new = deS([1000.0,200.0];key=:u0)
DEsteady
  func: #17 (function of type var"#17#18")
  p: Float64 1.0
  u0: LArray{Float64,1,Array{Float64,1},(:s1, :s2)}
  method: Tsit5 Tsit5()

julia> deS_new.u0
2-element LArray{Float64,1,Array{Float64,1},(:s1, :s2)}:
 :s1 => 1000.0
 :s2 => 200.0
```
"""
@with_kw struct DEsteady <: DEmeta
    func
    p
    u0
    method = SSRootfind()
end

function (self::DEsteady)(u0; key=:u0)
    return update_struct_broadcast(self, u0, key)
end

@with_kw struct Domain
    low::Number
    high::Number
    function Domain(low, high)
        if low > high
            @warn("Domain.low should <= Domain.high. Reset to Domain($low, $high)")
            new(high, low)
        else
            new(low, high)
        end
    end
end

"""
    update_struct_broadcast(var_struc, vals, key)

Update a field of struct varibale by broadcasting.

Argument
--------
-  `var_struc`: A struct variable.
- `vals`: A `Vector` which has the same length with `Base.getfield(var_struc, key)` for update the new struct variable.
- `key` : A `Symbol` which can apply `Base.getfield(var_struc, key)`

Return
------
- Struct variable similar to `var_struc` but with renew `vals` on `key` site.
"""
function update_struct_broadcast(var_struc, vals, key)
    vals_ = deepcopy(getfield(var_struc, key)) # derive the structure of `var_struc.key`
    vals_ .= vals #broadcast method
    struc = typeof(var_struc) # type of the struct variable

    fnames = map(fieldnames(struc)) do n
        if n == key
            v = vals_
        else
            v = getfield(var_struc, n)
        end

        return v
    end
    return typeof(var_struc)(fnames...) # use type of `var_struc` to create new struct variable
end


text:
## Make a simple enough GUI to plot variables in a data frame
using GtkInteract
using RDatasets, DataFrames
using Plots
backend(:immerse)

cars = dataset("MASS", "Cars93")
nms = names(cars)

## vartypes
nums = filter(nm -> eltype(cars[nm]) <: Number, nms)
bools = filter(nm -> eltype(cars[nm]) <: Bool, nms)
facs = filter(nm -> isa(cars[nm], PooledDataArray), nms)
facs = setdiff(facs, [:Manufacturer])
others = setdiff(nms, vcat(nums, bools, facs))

nums = map(string, nums)
xvar = dropdown(nums, label="X variable", value_label=nums[1])
yvar = dropdown(nums, label="Y variable", value_label=nums[2])
vb = formlayout()
append!(vb, [xvar, yvar])

fb = formlayout()
factors = Dict()
for fac in facs
    levs = levels(cars[fac])
    widget = buttongroup(levs, label=string(fac), value=levs)
    factors[fac] = widget
    push!(fb, widget)
end

## A graphics device
cg = cairographic()

## show toolbar
## * a togglebutton to hold state on whether titles should be displayed
## * a button to display an "about" message
do_titles = togglebutton(value=true, label="titles?")

about = button("about")

## layout toolbar with space between the items.
tb = toolbar(do_titles, vskip(),  about |> tooltip("Simple GUI"))

w = window(title="simple GUI");
b = vbox(tb);
push!(w, b)
push!(b, hbox(vbox(halign(:start, GtkInteract.bold("Select variables: ")),
                 vb,
                 halign(:start, GtkInteract.bold("Filter by: ")),
                   fb),
              cg |> grow |> padding(5)
              ))

## How to generate the graphic: first subset data, then plot
function make_graphic(args...)
    df = copy(cars)
    for fac in facs
        df = df[ Bool[val in value(factors[fac]) for val in df[fac]], :]
    end

    if size(df)[1] > 0
        xs = df[symbol(value(xvar))]
        ys = df[symbol(value(yvar))]
        p = scatter(xs, ys)
        if value(do_titles)
            plot!(p, title="Cars93", xlabel = value(xvar), ylabel=value(yvar))
        end
        push!(cg, p)
    else
        info("No cases left after selection")
    end
end

## actually show GUI
display(w)

## Hook up signals

f(args...) = info("A simple GUI") #messagebox("a simple GUI")
## ## XXX THIS IS FLAKY XXX
map(f, about)

## make update happen when a widget is updated
map(vcat(do_titles, xvar, yvar, values(factors)...)...) do args...
    make_graphic()
end


text:
"""
    waveguide_domains = defect_waveguide(domain::Domain, x0, y0, direction, width)

domain should be of the :pc type.
"""
function pc_waveguide_domains(pc_domain; width, direction, x0=0, y0=0, waveguide_number)

        xb0, yb0 = bravais_coordinates(x0, y0, pc_domain)
        if direction ∈ [:u, :up, :n, :north, :t, :top]
            axis = :b
            N1_bulk = floor(Int,yb0/pc_domain.lattice.b)
            N2_bulk = +Int(1e4)
            Nt = floor(Int,xb0/pc_domain.lattice.a)
            N_width = ceil(Int,width/pc_domain.lattice.a)
            bulk_domain_type = :bulk_pc_waveguide_y
            which_asymptote = :top
        elseif direction ∈ [:b, :bottom, :s, :south, :d, :down]
            axis = :b
            N1_bulk = -Int(1e4)
            N2_bulk = floor(Int,yb0/pc_domain.lattice.b)
            Nt = floor(Int,xb0/pc_domain.lattice.a)
            N_width = ceil(Int,width/pc_domain.lattice.a)
            bulk_domain_type = :bulk_pc_waveguide_y
            which_asymptote = :bottom
        elseif direction ∈ [:r, :right, :e, :east]
            axis = :a
            N1_bulk = floor(Int,xb0/pc_domain.lattice.a)
            N2_bulk = +Int(1e4)
            Nt = floor(Int,yb0/pc_domain.lattice.b)
            N_width = ceil(Int,width/pc_domain.lattice.b)
            bulk_domain_type = :bulk_pc_waveguide_x
            which_asymptote = :right
        elseif direction ∈ [:l, :left, :w, :west]
            axis = :a
            N1_bulk = -Int(1e4)
            N2_bulk = floor(Int,xb0/pc_domain.lattice.a)
            Nt = floor(Int,yb0/pc_domain.lattice.b)
            N_width = ceil(Int,width/pc_domain.lattice.b)
            bulk_domain_type = :bulk_pc_waveguide_x
            which_asymptote = :left
        elseif direction ∈ [:h, :horizontal, :x]
            axis = :a
            N1_bulk = -Int(1e4)
            N2_bulk = +Int(1e4)
            Nt = floor(Int,yb0/pc_domain.lattice.b)
            N_width = ceil(Int,width/pc_domain.lattice.b)
            bulk_domain_type = :bulk_pc_waveguide_x
            which_asymptote1 = :left
            which_asymptote2 = :right
        elseif direction ∈ [:v, :vertical, :y]
            axis = :b
            N1_bulk = -Int(1e4)
            N2_bulk = +Int(1e4)
            Nt = floor(Int,xb0/pc_domain.lattice.a)
            N_width = ceil(Int,width/pc_domain.lattice.a)
            bulk_domain_type = :bulk_pc_waveguide_y
            which_asymptote1 = :bottom
            which_asymptote2 = :top
        else
            throw(ArgumentError("invalid waveguide direction $(direction)"))
        end

        defect_bulk = line_defect(pc_domain, axis, N1_bulk, N2_bulk, Nt; N_width=N_width)
        wvg_bulk = vcat(Domain(defect_bulk; :domain_type => bulk_domain_type))

        defect_bnd = line_defect(pc_domain, axis, -Int(1e4), Int(1e4), Nt; N_width=N_width)

        if direction ∈ [:v, :vertical, :y, :x, :horizontal, :h]
            wvg = vcat(
                    Domain(defect_bnd; :which_waveguide => waveguide_number, :which_asymptote => which_asymptote1, :domain_type => :pc_waveguide),
                    Domain(pc_domain; :which_waveguide => waveguide_number, :which_asymptote => which_asymptote1, :domain_type => :pc_waveguide_background, :is_in_domain => Universe()),
                    Domain(defect_bnd; :which_waveguide => waveguide_number+1, :which_asymptote => which_asymptote2, :domain_type => :pc_waveguide),
                    Domain(pc_domain; :which_waveguide => waveguide_number+1, :which_asymptote => which_asymptote2, :domain_type => :pc_waveguide_background, :is_in_domain => Universe())
                    )
        else
            wvg = vcat(
                    Domain(defect_bnd; :which_waveguide => waveguide_number, :which_asymptote => which_asymptote, :domain_type => :pc_waveguide),
                    Domain(pc_domain; :which_waveguide => waveguide_number, :which_asymptote => which_asymptote, :domain_type => :pc_waveguide_background, :is_in_domain => Universe())
                    )

        end

    return vcat(wvg, wvg_bulk)
end

"""
    sim = add_defect_waveguide(sim::Simulation; x0, y0, direction, width)
    sys = add_defect_waveguide(sys::System; x0, y0, direction, width)
"""
function add_pc_waveguide(sys::System; width::Number, direction::Symbol, x0::Number, y0::Number)
    waveguide_number = length(sys.waveguides)+1
    temp_sys = System(sys.domains[.!isWaveguide.([sys.domains[i] for i ∈ eachindex(sys.domains)]) .& .!isBulkWaveguide.([sys.domains[i] for i ∈ eachindex(sys.domains)]) .& .!isDefect.([sys.domains[i] for i ∈ eachindex(sys.domains)]) ])
    domain = which_domain(x0,y0, Boundary(∂Ω=[-Inf Inf;-Inf Inf],bc=DirichletBC), temp_sys)
    waveguide_domains = pc_waveguide_domains(temp_sys.domains[domain]; x0=x0, y0=y0, direction=direction, width=width, waveguide_number=waveguide_number)
    return System(vcat(waveguide_domains,sys.domains))
end
function add_pc_waveguide(sim::Simulation; x0=0, y0=0, direction, width)
    sys = add_defect_waveguide(sim.sys; x0=x0, y0=y0, direction=direction, width=width)
    return Simulation(sim; :sys => sys)
end


text:
abstract type Cycle end
struct V <: Cycle end
struct W <: Cycle end
struct F <: Cycle end

struct Grid{TA, Pre, Post}
    A::TA
    presmoother::Pre
    postsmoother::Post
end

Base.size(grid::Grid) = size(grid.A)
Base.eltype(grid::Grid) = eltype(grid.A)

struct MultigridWorkspace{Tx}
    x::Vector{Tx}
    b::Vector{Tx}
    r::Vector{Tx}
end

struct Multigrid{TA, Tx, Coarse}
    grids::Vector{TA}
    coarsesolver::Coarse
    workspace::MultigridWorkspace{Tx}
end

Multigrid(A::Poisson2{T}; kwargs...) where {T} = Multigrid(T, A; kwargs...)

function Multigrid(
    T,
    A::Poisson2;
    presmoother = RedBlackGaussSeidel(ForwardSweep(), iter=2),
    coarsesolver = EarlyStopping(RedBlackGaussSeidel(ForwardSweep(), iter=100), sqrt(eps(T)), 1000),
    postsmoother = RedBlackGaussSeidel(ForwardSweep(), iter=2),
    maxlevel = max(2, floor(Int, log2(minimum(size(A)))-2)),
)
    maxlevel > 1 || throw(ArgumentError(""))

    sz = size(A)
    Tx = Array{T}

    presmoother = expand(presmoother, maxlevel)
    postsmoother = expand(postsmoother, maxlevel)
    grids = Vector{Grid}([Grid(A, presmoother[1], postsmoother[1])])

    x = Tx(undef, sz .* 0) # placeholder
    b = Tx(undef, sz .* 0) # placeholder
    r = tfill!(Tx(undef, sz), zero(T))
    workspace = MultigridWorkspace([x], [b], [r])

    for l in 2:maxlevel
        Ac = restrict(grids[l-1].A)
        push!(grids, Grid(Ac, presmoother[l], postsmoother[l]))

        szc = size(Ac)
        x = tfill!(Tx(undef, szc), zero(T))
        b = tfill!(Tx(undef, szc), zero(T))
        r = tfill!(Tx(undef, szc), zero(T))
        push!(workspace.x, x)
        push!(workspace.b, b)
        push!(workspace.r, r)
    end

    return Multigrid(grids, coarsesolver, workspace)
end

function solve(M::Multigrid, b::AbstractArray, args...; kwargs...)
    return solve!(tzero(b), M, b, args...; kwargs...)
end

function solve!(
    x::AbstractArray{<:AbstractFloat, N},
    M::Multigrid,
    b::AbstractArray{T, N},
    cycle::Cycle = V();
    atol::Real = zero(T),
    rtol::Real = sqrt(eps(T)),
    maxit::Integer = 100,
    residual::Bool = true,
    log::Bool = false,
    verbose::Bool = false,
) where {T<:AbstractFloat, N}
    size(x) == size(b) || throw(DimensionMismatch())
    size(x) == size(M.grids[1].A) || throw(DimensionMismatch())

    M.workspace.x[1] = x
    M.workspace.b[1] = b

    A = M.grids[1].A

    if residual
        nr = norm(b)
        ϵ = max(rtol * nr, atol)

        if log
            ch = [nr]
        end

        if verbose
            @printf("==== Multigrid ====\n")
            @printf("Depth: %d\n", length(M.grids))
            @printf("Tolerance: %1.4e\n", ϵ)
            @printf("%6s %3s %9s %7s\n", "iter", " ", "resnorm", "ratio")
            @printf("[%3d/%3d]  %1.4e\n", 0, maxit, nr)
        end
    end

    for i in 1:maxit
        cycle!(M, cycle)

        if residual
            nr0 = nr
            nr = __norm_residual(b, A, x)

            if log
                push!(ch, nr)
            end

            if verbose
                @printf("[%3d/%3d]  %1.4e  %0.2f\n", i, maxit, nr, nr / nr0)
            end

            if residual && nr ≤ ϵ
                break
            end
        end
    end

    if residual && verbose
        println()
    end

    return log ? (x, ch) : x
end

@noinline function cycle!(M::Multigrid, cycle::Cycle, l::Integer = 1)
    A = M.grids[l].A
    Ac = M.grids[l+1].A

    presmooth! = M.grids[l].presmoother
    postsmooth! = M.grids[l].postsmoother

    x = M.workspace.x[l]
    b = M.workspace.b[l]
    r = M.workspace.r[l]

    xc = M.workspace.x[l+1]
    bc = M.workspace.b[l+1]

    presmooth!(x, A, b)
    residual!(r, b, A, x)
    restrict!(bc, Ac, r)

    xc = tfill!(xc, 0, Ac)

    if l+1 == length(M.grids)
        M.coarsesolver(xc, Ac, bc)
    else
        _cycle!(M, cycle, l+1)
    end

    correct_prolong!(x, A, xc)
    postsmooth!(x, A, b)

    return M
end

function _cycle!(M, cycle::V, l)
    cycle!(M, cycle, l)
end

function _cycle!(M, cycle::W, l)
    cycle!(M, cycle, l)
    cycle!(M, cycle, l)
end

function _cycle!(M, cycle::F, l)
    cycle!(M, cycle, l)
    cycle!(M, V(), l)
end


text:
push!(LOAD_PATH, "./")

using ExpSim

simulate(ARGS[1], 800, 480)


text:
#=
Response surfaces implementantion, following:
"A Taxonomy of Global Optimization Methods Based on Response Surfaces"
by DONALD R. JONES
=#
mutable struct RadialBasis{F,X,Y,B,C} <: AbstractSurrogate
    phi::F
    dim_poly::Int
    x::X
    y::Y
    bounds::B
    coeff::C
end

"""
Calculates current estimate of array value 'val' with respect to RadialBasis object.
"""
function (rad::RadialBasis)(val)
    n = length(rad.x)
    d = length(rad.x[1])
    q = rad.dim_poly
    lb, ub = rad.bounds
    sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:d)
    mean_half_diameter = sum_half_diameter/d
    central_point = _center_bounds(first(rad.x), lb, ub)
    approx = sum(rad.coeff[i]*rad.phi(val .- rad.x[i]) for i = 1:n)
    for k = 1:q
        approx += rad.coeff[n+k]*centralized_monomial(val, k-1, mean_half_diameter, central_point)
    end
    return approx
end

function _scaled_chebyshev(x, k, lb, ub)
    return cos(k*acos(-1 + 2*(x-lb)/(ub-lb)))
end

_center_bounds(x::Tuple, lb, ub) = ntuple(i -> (ub[i] - lb[i])/2, length(x))
_center_bounds(x, lb, ub) = (ub .- lb) ./ 2

"""
Calculates current estimate of value 'val' with respect to RadialBasis object.
"""
function (rad::RadialBasis)(val::Number)
    approx = zero(eltype(rad.x[1]))
    n = length(rad.x)
    q = rad.dim_poly
    for i = 1:n
        approx += rad.coeff[i]*rad.phi(val .- rad.x[i])
    end
    for k = 1:q
        approx += rad.coeff[n+k] * _scaled_chebyshev(val, k-1, rad.bounds[1], rad.bounds[2])
    end
    return approx
end

"""
    RadialBasis(x,y,a::Number,b::Number,phi::Function,q::Int)

Constructor for RadialBasis surrogate
- (x,y): sampled points
- '(a,b)': interval of interest
-'phi': radial basis of choice
-'q': number of polynomial elements
"""
function RadialBasis(x, y, lb::Number, ub::Number, phi::Function, q::Int)
    coeff = _calc_coeffs(x, y, lb, ub, phi, q)
    return RadialBasis(phi, q, x, y, (lb, ub), coeff)
end

function _calc_coeffs(x, y, lb, ub, phi, q)
    n = length(x)
    m = n + q
    D = zeros(eltype(x[1]), m, m)
    d = zeros(eltype(x[1]), m)

    @inbounds for i = 1:n
        d[i] = y[i]
        for j = 1:n
            D[i,j] = phi(x[i] .- x[j])
        end
        if i <= n
            for k = 1:q
                D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    coeff = D_sym\d
    return coeff
end

"""
RadialBasis(x,y,bounds,phi::Function,q::Int)

Constructor for RadialBasis surrogate

- (x,y): sampled points
- bounds: region of interest of the form [[a,b],[c,d],...,[w,z]]
- phi: radial basis of choice
- q: number of polynomial elements
"""
function RadialBasis(x, y, bounds, phi::Function, q::Int)
    coeff = _calc_coeffs(x, y, bounds, phi, q)
    return RadialBasis(phi, q, x, y, bounds, coeff)
end

function _calc_coeffs(x, y, bounds, phi, q)
    n = length(x)
    d = length(x[1])
    lb, ub = bounds
    central_point = _center_bounds(first(x), lb, ub)
    sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:d)
    mean_half_diameter = sum_half_diameter / d
    m = n+q
    D = zeros(eltype(x[1]), m, m)
    d = zeros(eltype(x[1]), m)

    @inbounds for i = 1:n
        d[i] = y[i]
        for j = 1:n
            D[i,j] = phi(x[i] .- x[j])
        end
        if i < n + 1
            for k = 1:q
                D[i,n+k] = centralized_monomial(x[i], k-1, mean_half_diameter, central_point)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    return coeff = D_sym\d
end

"""
    centralized_monomial(vect,alpha,mean_half_diameter,central_point)

Returns the value at point vect[] of the alpha degree monomial centralized.

#Arguments:
-'vect': vector of points i.e [x,y,...,w]
-'alpha': degree
-'mean_half_diameter': half diameter of the domain
-'central_point': central point in the domain
"""
function centralized_monomial(vect,alpha,mean_half_diameter,central_point)
    if iszero(alpha) return one(eltype(vect)) end
    centralized_product = prod(vect .- central_point)
    return (centralized_product / mean_half_diameter)^alpha
end

"""
    add_point!(rad::RadialBasis,new_x,new_y)

Add new samples x and y and updates the coefficients. Return the new object radial.
"""
function add_point!(rad::RadialBasis,new_x,new_y)
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(rad.bounds[1])>1)
        push!(rad.x,new_x)
        push!(rad.y,new_y)
        if length(rad.bounds[1]) == 1
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds[1],rad.bounds[2],rad.phi,rad.dim_poly)
        else
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds,rad.phi,rad.dim_poly)
        end
    else
        append!(rad.x,new_x)
        append!(rad.y,new_y)
        if length(rad.bounds[1]) == 1
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds[1],rad.bounds[2],rad.phi,rad.dim_poly)
        else
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds,rad.phi,rad.dim_poly)
        end
    end
    nothing
end


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "GitLabPages" begin
    @testset "Plugin creation" begin
        p = GitLabPages()
        @test p.gitignore == ["/docs/build/", "/docs/site/"]
        @test isempty(p.assets)
        p = GitLabPages(; assets=[test_file])
        @test p.assets == [test_file]
        @test_throws ArgumentError GitLabPages(; assets=[fake_path])
    end

    @testset "Badge generation" begin
        p = GitLabPages()
        @test badges(p, me, test_pkg) ==  [
            "[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://$me.gitlab.io/$test_pkg.jl/dev)"
        ]
    end

    @testset "File generation" begin
        p = GitLabPages()
        @test gen_plugin(p, t, test_pkg) == ["docs/"]
        @test isdir(joinpath(pkg_dir, "docs"))
        @test isfile(joinpath(pkg_dir, "docs", "make.jl"))
        make = readchomp(joinpath(pkg_dir, "docs", "make.jl"))
        @test occursin("assets=String[]", make)
        @test !occursin("deploydocs", make)
        @test isdir(joinpath(pkg_dir, "docs", "src"))
        @test isfile(joinpath(pkg_dir, "docs", "src", "index.md"))
        index = readchomp(joinpath(pkg_dir, "docs", "src", "index.md"))
        @test occursin("autodocs", index)
        rm(joinpath(pkg_dir, "docs"); recursive=true)
        p = GitLabPages(; assets=[test_file])
        @test gen_plugin(p, t, test_pkg) == ["docs/"]
        make = readchomp(joinpath(pkg_dir, "docs", "make.jl"))
        # Check the formatting of the assets list.
        @test occursin(
            strip("""
            assets=[
                    "assets/$(basename(test_file))",
                ]
            """),
            make,
        )
        @test isfile(joinpath(pkg_dir, "docs", "src", "assets", basename(test_file)))
        rm(joinpath(pkg_dir, "docs"); recursive=true)
    end

    @testset "Package generation with GitLabPages plugin" begin
        temp_dir = mktempdir()
        t = Template(; user=me, dir=temp_dir, plugins=[GitLabCI(), GitLabPages()])
        generate(test_pkg, t; gitconfig=gitconfig)

        gitlab = read(joinpath(t.dir, test_pkg, ".gitlab-ci.yml"), String)
        @test occursin("pages:", gitlab)
    end
end

rm(pkg_dir; recursive=true)


text:
"""
    function threegrouptransmission(du,u,p,t)

In-place creation of vector field for the three group transmission model.
"""
function threegrouptransmission(du,u,p,t)
    #Get state
	S₁,E₁,I₁,R₁,W₁,C₁,C_sero₁,S₂,E₂,I₂,R₂,W₂,C₂,C_sero₂,S₃,E₃,I₃,R₃,W₃,C₃,C_sero₃ = u
    #Get parameters
	R₀,α,γ,ϵ,σ,N₁,N₂,N₃,ω,schooleffect,ct_min1,ct_min2,ct_min3 = p
	#Adjust for schools being open or shut
    Reff = R₀*(t < 24 || t ≥ 319) + R₀*schooleffect*(t≥ 24 && t< 319)
	#Contact rates for each group
    ct1 = KenyaCoVSD.ct_kenya(t,ct_min1)
    ct2 = KenyaCoVSD.ct_kenya(t,ct_min2)
    ct3 = KenyaCoVSD.ct_kenya(t,ct_min3)
    #Force of infection on each group
    λ₁ = γ*Reff*(ct1*ϵ*I₁ + ct2*(1-ϵ)*I₂*(N₁/(N₁ + N₃)) + ct3*(1-ϵ)*I₃*(N₁/(N₁ + N₂)))
    λ₂ = γ*Reff*(ct1*(1-ϵ)*I₁*(N₂/(N₂ + N₃)) + ct2*ϵ*I₂ + ct3*(1-ϵ)*I₃*(N₂/(N₂ + N₁)))
    λ₃ = γ*Reff*(ct1*(1-ϵ)*I₁*(N₃/(N₂ + N₃)) + ct2*(1-ϵ)*I₂*(N₃/(N₃ + N₁)) + ct3*ϵ*I₃)
	#Dynamics
	du[1] = -ct1*(S₁/N₁)*λ₁
	du[2] = ct1*((S₁+σ*W₁)/N₁)*λ₁ - α*E₁
	du[3] = α*E₁ - γ*I₁
	du[4] = γ*I₁ - ω*R₁
	du[5] = ω*R₁ - σ*ct1*(W₁/N₁)*λ₁
	du[6] = ct1*((S₁+σ*W₁)/N₁)*λ₁
    du[7] = ct1*(S₁/N₁)*λ₁

	du[8] = -ct2*(S₂/N₂)*λ₂
	du[9] = ct2*((S₂+σ*W₂)/N₂)*λ₂ - α*E₂
	du[10] = α*E₂ - γ*I₂
	du[11] = γ*I₂ - ω*R₂
	du[12] = ω*R₂ - σ*ct2*(W₂/N₂)*λ₂
	du[13] = ct2*((S₂+σ*W₂)/N₂)*λ₂
    du[14] = ct2*(S₂/N₂)*λ₂

	du[15] = -ct3*(S₃/N₃)*λ₃
	du[16] = ct3*((S₃+σ*W₃)/N₃)*λ₃ - α*E₃
	du[17] = α*E₃ - γ*I₃
	du[18] = γ*I₃ - ω*R₃
	du[19] = ω*R₃ - σ*ct3*(W₃/N₃)*λ₃
	du[20] = ct3*((S₃+σ*W₃)/N₃)*λ₃
    du[21] = ct3*(S₃/N₃)*λ₃
	return nothing
end

u0_threegroup = [4.3e6*0.25 - 100.,100.,0.,0.,0.,0.,0.,
        4.3e6*0.5 - 100.,100.,0.,0.,0.,0.,0.,
        4.3e6*0.25 - 100.,100.,0.,0.,0.,0.,0.]
f_threegroup = ODEFunction(threegrouptransmission;syms = [:S₁,:E₁,:I₁,:R₁,:W₁,:C₁,:C_sero₁,
															:S₂,:E₂,:I₂,:R₂,:W₂,:C₂,:C_sero₂,
															:S₃,:E₃,:I₃,:R₃,:W₃,:C₃,:C_sero₃])
prob_three_group = ODEProblem(f_threegroup,u0_threegroup,(0,365))

"""
    function ll_threegroup_newvariant_infboost(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate,cts;σ = 0.16,ω = 1/180,ct_min2 = 0.445)

Log-likelihood of the case and serology data in `model` assuming: 1) Three group dynamics, and 2) that new variant is more transmissible rather than immune evading.
"""
function ll_threegroup_newvariant_infboost(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate;σ = 0.16,ω = 1/180)
    @unpack PCR_cases,sero_cases,baseline_sero_array,PCR_array,sero_sensitivity,sero_specificity,N,M_BB,prob,α,γ,relative_testing_rate = model
    @unpack R₀,ϵ,schooleffect,ct_min1,ct_min2,ct_min3,χ₁,χ₂,χ₃,p_test1,p_test2,p_test3,z_pch1,z_pch2,z_pch3,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    N₁ = 0.25*N
    N₂ = 0.5*N
    N₃ = 0.25*N
    pch1,pch2,pch3 = softmax([z_pch1,z_pch2,z_pch3])

    #Set variance scalers
    clustering_factor_PCR = 0.5
    M_PCR = 30.
    T = eltype(R₀)
	# if T <: Real
	# 	T = eltype(cts)
	# end
    #Set transmission parameters with daily cts appended on end
    p = convert.(T,[R₀,α,γ,ϵ,σ,N₁,N₂,N₃,ω,schooleffect,ct_min1,ct_min2,ct_min3])
    u0 = convert.(T,[N₁-E₀,E₀,0.,0.,0.,0.,0.,
                    N₂-E₀,E₀,0.,0.,0.,0.,0.,
                    N₃-E₀,E₀,0.,0.,0.,0.,0.])

    #Create new variant introduction
    function new_variant_effect!(integrator)
        integrator.p[1] *= extra_transmissibility
        integrator.u[1] -= influx_exposed_new_variant
        integrator.u[2] += influx_exposed_new_variant
        integrator.u[1 + 7] -= influx_exposed_new_variant
        integrator.u[2 + 7] += influx_exposed_new_variant
        integrator.u[1 + 14] -= influx_exposed_new_variant
        integrator.u[2 + 14] += influx_exposed_new_variant
    end

    janendpoint = (Date(2021,1,30) - Date(2020,2,20)).value
	aprilendpoint = (Date(2021,4,30) - Date(2020,2,20)).value
    variant_cb = PresetTimeCallback([janendpoint],new_variant_effect!)

    #Sero-waning
    sero_array = vcat(baseline_sero_array[1:30],[(1-seroreversionrate)^k for k in 1:length(baseline_sero_array[31:end])])

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    LL = T(0.)
    # try
        sol = solve(prob, BS3();tspan = (0,aprilendpoint),reltol = 1e-3,
					callback = variant_cb,
					u0=u0,
					p=p,
					saveat = 1,
					verbose = false,
                    isoutofdomain=(u,p,t) -> any(x -> x < 0, u))

        ι₁ = diff(sol[:C₁])
        ι₂ = diff(sol[:C₂])
        ι₃ = diff(sol[:C₃])
        ι_sero₁ = diff(sol[:C_sero₁])
        ι_sero₂ = diff(sol[:C_sero₂])
        ι_sero₃ = diff(sol[:C_sero₃])

        PCR₁ = KenyaCoVSD.simple_conv(ι₁,PCR_array)
        PCR₂ = KenyaCoVSD.simple_conv(ι₂,PCR_array)
        PCR₃ = KenyaCoVSD.simple_conv(ι₃,PCR_array)
        weighted_PCR_ptest = (p_test1.*PCR₁ .+ p_test2.*PCR₂ .+ p_test3.*PCR₃).*1e-4
        weighted_prop_PCR = pch1.*(χ₁.*PCR₁./((χ₁-1).*PCR₁ .+ N₁)) .+ pch2.*(χ₂.*PCR₂./((χ₂-1).*PCR₂ .+ N₂)) .+ pch3.*(χ₃.*PCR₃./((χ₃-1).*PCR₃ .+ N₃))
        boosted_weighted_prop_PCR = pch1.*(χ_boost.*χ₁.*PCR₁./((χ_boost.*χ₁-1).*PCR₁ .+ N₁)) .+ pch2.*(χ_boost.*χ₂.*PCR₂./((χ_boost.*χ₂-1).*PCR₂ .+ N₂)) .+ pch3.*(χ_boost.*χ₃.*PCR₃./((χ_boost.*χ₃-1).*PCR₃ .+ N₃))

    	sero = sero_sensitivity.*KenyaCoVSD.simple_conv(ι_sero₁.+ι_sero₂.+ι_sero₃,sero_array)

        #Calculate log-likelihood for PCR testing
        for t in 55:janendpoint
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = relative_testing_rate[t]*weighted_PCR_ptest[t] + 0.001 #Covert approximately from total probability of detection in % to daily probability of detection
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = weighted_prop_PCR[t] + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing --- proportion postive
            end
        end

        #Calculate log-likelihood for PCR testing after new variant introduction
        for t in (janendpoint+1):min(size(PCR_cases,1),size(PCR₁,1))
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = p_test_boost*relative_testing_rate[t]*weighted_PCR_ptest[t] + 0.001 #Covert approximately from total probability of detection in % to daily probability of detection
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = boosted_weighted_prop_PCR[t] + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing that day in that age group
            end
        end
        #Calculate log-likelihood contribution from serological testing
        for t in 1:min(size(sero_cases,1),size(sero,1))
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation
            #Overall serological testing
            if sero_cases[t,1,1] + sero_cases[t,2,1] > 0
                p_sero_pred  = sero[t]/N
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,1] + sero_cases[t,2,1],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,1])#Likelihood contribution from sero testing
            end

        end
    # catch
    #     LL = T(-Inf)
    # end
    return LL::T
end

function threegroup_priors(θ)
	@unpack R₀,ϵ,schooleffect,ct_min1,ct_min2,ct_min3,χ₁,χ₂,χ₃,p_test1,p_test2,p_test3,z_pch1,z_pch2,z_pch3,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ
    LP = 0.
    LP += logpdf(Gamma(3,2.5/3),R₀)
    LP += logpdf(Beta(7,3),ϵ)
    LP += logpdf(Beta(8,2),schooleffect)
    LP += logpdf(Beta(7,3),ct_min1)
    LP += logpdf(Beta(7,3),ct_min2)
    LP += logpdf(Beta(55,45),ct_min3)
    LP += logpdf(Gamma(3,4/3),χ₁)
    LP += logpdf(Gamma(3,4/3),χ₂)
    LP += logpdf(Gamma(3,1.5/3),χ₃)
    LP += logpdf(Gamma(3,1/3),p_test1)#In the region of 1% detection rate
    LP += logpdf(Gamma(3,1/3),p_test2)#In the region of 1% detection rate
    LP += logpdf(Gamma(3,5/3),p_test3)#In the region of 5% detection rate
    LP += logpdf(Normal(0,1.0),z_pch1)
	LP += logpdf(Normal(0,1.0),z_pch2)
	LP += logpdf(Normal(0,1.0),z_pch3)
    LP += logpdf(Gamma(10,1.5/10),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,100/3),E₀)
    return LP
end

## Create three group model
nai_three_group = deepcopy(basic_nai_model)
nai_three_group.prob = prob_three_group
nai_three_group.log_likelihood = ll_threegroup_newvariant_infboost
nai_three_group.log_priors = threegroup_priors

#Define variable transformation for three group model
trans_three_groups = as((R₀ = as(Real, 0.0, 7.0),
        ϵ = as(Real,0.5,0.999),
        schooleffect = as(Real,0.25,1.0),
        ct_min1 = as(Real, 0.1, 0.999),
        ct_min2 = as(Real, 0.1, 0.999),
        ct_min3 = as(Real, 0.1, 0.999),
        χ₁=as(Real, 0.1, 30.0),χ₂=as(Real, 0.1, 30.0),χ₃=as(Real, 0.1, 30.0),
        p_test1=as(Real, 0.0, 50.0),p_test2=as(Real, 0.0, 50.0),p_test3=as(Real, 0.0, 50.0),
        z_pch1 = asℝ,z_pch2 = asℝ,z_pch3 = asℝ,
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 0.5, 3.0),
        χ_boost= as(Real, 0.5, 3.0),
        E₀=as(Real, 0.0, 10e3)))

D = Diagonal(0.1*ones(TransformVariables.dimension(trans_three_groups)))


text:
@testset "Maximum likelihood" begin
    b = 1.0
    M0 = 0.5
    mags = rand(GRclassic(b,M0), 100000)
    filter!(x->rand() < cdf(Normal(0.8,0.2),x), mags)
    ΔM = 0.1

    M0 = 1.25 - ΔM/2
    mags = mags[mags .>= M0]

    @test round.( bvalue(fit_mle(GRclassic,mags,M0)), digits=3) == 0.995
    @test round.( fit_stderr(fit_mle(GRclassic,mags,M0), mags), digits=3) == 0.007

    @test round.( bvalue(fit_binned(GRclassic,mags,M0,ΔM)), digits=3) == 0.999
    @test round.( fit_stderrbinned(fit_binned(GRclassic,mags,M0,ΔM), mags, ΔM), digits=3) == 0.007

end


text:
# ---
# title: 175. Combine Two Tables
# id: problem175
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/combine-two-tables/description/>
# hidden: true
# ---
#
# Table: `Person`
#
#
#
#     +-------------+---------+
#     | Column Name | Type    |
#     +-------------+---------+
#     | PersonId    | int     |
#     | FirstName   | varchar |
#     | LastName    | varchar |
#     +-------------+---------+
#     PersonId is the primary key column for this table.
#
#
# Table: `Address`
#
#
#
#     +-------------+---------+
#     | Column Name | Type    |
#     +-------------+---------+
#     | AddressId   | int     |
#     | PersonId    | int     |
#     | City        | varchar |
#     | State       | varchar |
#     +-------------+---------+
#     AddressId is the primary key column for this table.
#
#
#
#
# Write a SQL query for a report that provides the following information for
# each person in the Person table, regardless if there is an address for each of
# those people:
#
#
#
#     FirstName, LastName, City, State
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
  rep(ind::Tuple)

axiliary for cnorm count how many times a block should be counted
```jldoctest
julia> rep((1,2,3))
6

julia> rep((1,2,2))
3

julia> rep((1,1,1))
1
```
"""
function rep(ind::Tuple)
  @inbounds c = counts([ind...])
  div(factorial(length(ind)), mapreduce(factorial, * , c))
end

"""
norm(bt::SymmetricTensor{T, m}, p::Union{Float64, Int}=2) where {T<:AbstractFloat, m}

Returns Float, a p-norm of bt, supported for p ≠ 0

```jldoctest

julia> te = [-0.112639 0.124715 0.124715 0.268717 0.124715 0.268717 0.268717 0.046154];

julia> st = convert(SymmetricTensor, (reshape(te, (2,2,2))));

julia> norm(st)
0.5273572868359742

julia> norm(st, 2.5)
0.4468668679541424

julia> norm(st, 1)
1.339089
```
"""
function norm(bt::SymmetricTensor{T, m}, p::Union{Float64, Int}=2) where {T<:AbstractFloat, m}
  p != 0 || throw(AssertionError("0-norm not supported"))
  ret = 0.
  for i in pyramidindices(m, bt.bln)
    @inbounds ret += rep(i) * norm(getblockunsafe(bt, i), p)^p
  end
  ret^(1/p)
end

"""
  cnorms(c::Vector{SymmetricTensor{T}})

Returns vector of Floats of norms of cumulants of order 3, ..., k , ..., m, normalsed
by √||C₂||ᵏ
.....

"""
function cnorms(c::Vector{SymmetricTensor{T}}) where T <: AbstractFloat
  n = norm(c[2])
  [norm(c[k])/(n)^(k/2) for k in 3:length(c)]
end


text:
using ClickDetection
using SampledSignals
using DSP
using Statistics
using Dates
using FileIO
using Test

@testset "Intervals" begin
    intervals1 = [1..2, 3..4]
    intervals2 = [1..2, 1.5..4]
    @test length(merge_intervals(intervals1)) == 2
    @test length(merge_intervals(intervals2)) == 1
end

@testset "Detection" begin
    wavfile = joinpath(@__DIR__, "data", "MARS_20180201_131212_clip.wav")
    audio = loadstreaming(wavfile)
    x = read(audio)
    close(audio)

    # preprocess raw audio before click detection
    x = convert.(Float32, x)
    fs = samplerate(x)
    butter = digitalfilter(Highpass(5e3, fs=fs), Butterworth(3))
    x = filtfilt(butter, x)

    # Actual click detection
    thresh = 1e-5
    window = 1e-3
    clicks1 = detect_clicks(x, thresh, window)
    @test series(clicks1[1]) === series(clicks1[2])
    c1 = first(clicks1)
    c2 = Click(c1)
    @test length(samples(c1)) == length(samples(c2))
    @test length(times(c1)) == length(times(c2))
    @test all(samples(c1) .== samples(c2))
    @test all(times(c1) .== times(c2))
    @test left(c1) == left(c2)
    @test right(c1) == right(c2)

    for c in [c1, c2]
        t = times(c)
        @test length(samples(c)) == length(t)
        @test left(c) == first(t)
        @test right(c) == last(t)
    end

    clicks2 = detect_clicks(Array(x), samplerate(x), thresh, window)
    @test series(first(clicks1)) !== series(first(clicks2))
    @test length(clicks1) == length(clicks2)
    @test typeof(samples(first(clicks1))) <: Array
    @test all(samples(first(clicks1)) .== samples(first(clicks2)))
    @test interval(first(clicks1)) == interval(first(clicks2))
    @test left(first(clicks1)) == left(first(clicks2))
    @test right(first(clicks1)) == right(first(clicks2))
    @test samplerate(first(clicks1)) == samplerate(first(clicks2))
    @test samplerate(first(clicks2)) == samplerate(first(clicks1))
    @test times(first(clicks1)) == times(first(clicks2))

    required_colnames = [:start, :stop, :samplerate, :waveform]

    df = as_dataframe(clicks1)
    @test all(in(names(df)).(required_colnames))
    @test size(df,1) == length(clicks1)
    df = as_dataframe(clicks1, DateTime(2019,1,1,0,0,0))
    @test size(df,1) == length(clicks1)
    @test all(in(names(df)).([required_colnames; :start_datetime; :stop_datetime]))
    df = as_dataframe(clicks1, DateTime(2019,1,1,0,0,0), "fakefile.wav")
    @test all(in(names(df)).([required_colnames; :start_datetime; :stop_datetime; :wavfile]))
    @test size(df,1) == length(clicks1)
    fakefiles = string.(1:length(clicks1))
    df = as_dataframe(clicks1, DateTime(2019,1,1,0,0,0), fakefiles)
    @test all(in(names(df)).([required_colnames; :start_datetime; :stop_datetime; :wavfile]))
    @test size(df,1) == length(clicks1)
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: PHE.mso 250 2007-04-27 16:32:02Z bicca $
#*------------------------------------------------------------------
type PHE_Geometry
	PHE_Geometry()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Chemical Components",
				:Hidden=>true
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Default=>3.14159265,
				:Hidden=>true,
				:Symbol=>"\\pi"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Auxiliar Constant",
				:Hidden=>true,
				:Default=>15
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Auxiliar Constant",
				:Hidden=>true,
				:Default=>14
			)),
			fill(constant(Dict{Symbol,Any}(
				:Brief=>"First constant in Kumar calculation for Pressure Drop",
				:Hidden=>true
			)),(N1)),
			fill(constant(Dict{Symbol,Any}(
				:Brief=>"Second constant in Kumar calculation for Pressure Drop",
				:Hidden=>true
			)),(N1)),
			fill(constant(Dict{Symbol,Any}(
				:Brief=>"First constant in Kumar calculation for Heat Transfer",
				:Hidden=>true
			)),(N2)),
			fill(constant(Dict{Symbol,Any}(
				:Brief=>"Second constant Kumar calculation for Heat Transfer",
				:Hidden=>true
			)),(N2)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			length(Dict{Symbol,Any}(
				:Brief=>"Vertical Ports Distance",
				:Lower=>0.1
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Total Number of Plates in The Whole Heat Exchanger",
				:Default=>25,
				:Symbol=>"N_{plates}"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Passes for Hot Side",
				:Symbol=>"Npasshot"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Passes for Cold Side",
				:Symbol=>"Npasscold"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Ports Diameter",
				:Lower=>1e-6,
				:Symbol=>"D_{ports}"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Plate Width",
				:Lower=>0.1
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Plate Pitch",
				:Lower=>0.1
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Plate Thickness",
				:Lower=>0.1
			)),
			conductivity(Dict{Symbol,Any}(
				:Brief=>"Plate Thermal Conductivity",
				:Default=>1.0,
				:Symbol=>"K_{wall}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Hot Side Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cold Side Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Enlargement Factor",
				:Lower=>1e-6,
				:Symbol=>"\\phi"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Total Effective  Area",
				:Lower=>1e-6,
				:Symbol=>"A_{total}",
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Port Opening  Area of Plate",
				:Lower=>1e-6,
				:Symbol=>"A_{ports}",
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Cross-Sectional Area for Channel Flow",
				:Lower=>1e-6,
				:Symbol=>"A_{channel}",
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Equivalent Diameter of Channel",
				:Lower=>1e-6,
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Corrugation Depth",
				:Lower=>1e-6,
				:Protected=>true
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Total Number of Channels in The Whole Heat Exchanger",
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Plate Vertical Distance between Port Centers",
				:Lower=>0.1,
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Compact Plate Pack Length",
				:Lower=>0.1,
				:Protected=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Plate Horizontal Distance between Port Centers",
				:Lower=>0.1,
				:Protected=>true
			)),
			[:PP,:NComp,:Pi,:N1,:N2,:Kp1,:Kp2,:Kc1,:Kc2,:M,:Lv,:Nplates,:NpassHot,:NpassCold,:Dports,:Lw,:pitch,:pt,:Kwall,:Rfh,:Rfc,:PhiFactor,:Atotal,:Aports,:Achannel,:Dh,:Depth,:Nchannels,:Lp,:Lpack,:Lh,],
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Pi::constant
	N1::DanaInteger
	N2::DanaInteger
	Kp1::Array{constant}
	Kp2::Array{constant}
	Kc1::Array{constant}
	Kc2::Array{constant}
	M::Array{molweight}
	Lv::length
	Nplates::DanaInteger
	NpassHot::DanaInteger
	NpassCold::DanaInteger
	Dports::length
	Lw::length
	pitch::length
	pt::length
	Kwall::conductivity
	Rfh::positive
	Rfc::positive
	PhiFactor::DanaReal
	Atotal::area
	Aports::area
	Achannel::area
	Dh::length
	Depth::length
	Nchannels::DanaInteger
	Lp::length
	Lpack::length
	Lh::length
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PHE_Geometry
function set(in::PHE_Geometry)
	#"Vector Length of constants for Kumar's calculating Pressure Drop"
	N1 = 15
	 #"Vector Length of constants for Kumar's calculating Heat Transfer"
	N2 = 14
	 #"First constant for Kumar's calculating Pressure Drop"
	Kp1 = [50,19.40,2.990,47,18.290,1.441,34,11.250,0.772,24,3.240,0.760,24,2.80,0.639]
	 #"Second constant for Kumar's calculating Pressure Drop"
	Kp2 = [1,0.589,0.183,1,0.652,0.206,1,0.631,0.161,1,0.457,0.215,1,0.451,0.213]
	 #"First constant for Kumar's calculating Heat Transfer"
	Kc1 = [0.718,0.348,0.718,0.400,0.300,0.630,0.291,0.130,0.562,0.306,0.108,0.562,0.331,0.087]
	 #"Second constant for Kumar's calculating Heat Transfer"
	Kc2 = [0.349,0.663,0.349,0.598,0.663,0.333,0.591,0.732,0.326,0.529,0.703,0.326,0.503,0.718]
	 #"Component Molecular Weight"
	M = PP.MolecularWeight()
	 #"Pi Number"
	Pi = 3.14159265
	 #"Plate Vertical Distance between Port Centers"
	Lp = Lv - Dports
	 #"Corrugation Depth"
	Depth=pitch-pt
	 #"Plate Horizontal Distance between Port Centers"
	Lh=Lw-Dports
	 #"Hydraulic Diameter"
	Dh=2*Depth/PhiFactor
	 #"Ports Area"
	Aports=0.25*Pi*Dports*Dports
	 #"Channel Area"
	Achannel=Depth*Lw
	 #"Pack Length"
	Lpack=Depth*(Nplates-1)+Nplates*pt
	 #"Total Number of  Channels"
	Nchannels = Nplates -1
	 #"Exchange Surface Area"
	Atotal =(Nplates-2)*Lw*Lp*PhiFactor
end
function atributes(in::PHE_Geometry,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Parameters for a gasketed plate heat exchanger."
	drive!(fields,_)
	return fields
end
PHE_Geometry(_::Dict{Symbol,Any})=begin
	newModel=PHE_Geometry()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#=
The nth term of the sequence of triangle numbers is given by, tn = ½n(n+1); so the first ten triangle numbers are:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

By converting each letter in a word to a number corresponding to its alphabetical position and adding these values we form a word value. For example, the word value for SKY is 19 + 11 + 25 = 55 = t10. If the word value is a triangle number then we shall call the word a triangle word.

Using words.txt (right click and 'Save Link/Target As...'), a 16K text file containing nearly two-thousand common English words, how many are triangle words?
=#

using DelimitedFiles

function calc()
  x = sort!(vec(readdlm(dirname(@__FILE__()) * "/../Resources/p42.txt", ',', AbstractString)))
  tn = [div(n*(n+1),2) for n in 1:100]
  trianglewords = 0
  for element in x
    wordsum = 0
    for character in element
      wordsum += Int(character)-64
    end
    trianglewords += in(wordsum, tn) == true ? 1 : 0
  end
  trianglewords
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module Collocated
using SparseArrays, LinearAlgebra
import sbp
using sbp.Metrics: CovariantBasis
using sbp.Sparse: block_matrix, block_matrix_insert

struct Operators2D
        nx::Int64
        ny::Int64
        Dx::AbstractArray
        Dy::AbstractArray
        Hx::AbstractArray
        Hy::AbstractArray
        Bx::AbstractArray
        By::AbstractArray
end

function boundary_matrix(n, is_sparse=true)
        if is_sparse
                B = spzeros(n, n)
        else
                B = zeros(n, n)
        end

        B[1,1] = -1
        B[end,end] = 1

        return B
end

function build_operators_2d(Dx::AbstractArray,
                            Dy::AbstractArray,
                            Hx::AbstractArray,
                            Hy::AbstractArray,
                            Bx::AbstractArray,
                            By::AbstractArray
                            )
        nx = size(Dx, 1)
        ny = size(Dy, 1)

        Ix = sparse(1.0I, nx, nx)
        Iy = sparse(1.0I, ny, ny)

        Dx2 = kron(Dx, Iy)
        Dy2 = kron(Ix, Dy)

        Hx2 = kron(Hx, Iy)
        Hy2 = kron(Ix, Hy)

        Bx2 = kron(Bx, Iy)
        By2 = kron(Ix, By)

        return Operators2D(nx, ny, Dx2, Dy2, Hx2, Hy2, Bx2, By2)

end

"""
Construct the covariant basis at the cell centers (mm) using SBP operators.
The mapping function f(x,y) = (fx, fy) must be defined at the nodes (pp).
"""
function build_covariant_basis(fx::AbstractArray, fy::AbstractArray,
                               op::Operators2D)
        x_r1 = op.Dx * fx
        x_r2 = op.Dy * fx
        y_r1 = op.Dx * fy
        y_r2 = op.Dy * fy

        return CovariantBasis(x_r1, y_r1, x_r2, y_r2)
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python")
include(joinpath(@__DIR__, "stats_evaluation_multi_target_localization.jl"))

function evaluate()
    data_dir = normpath(joinpath(@__DIR__, "../data",
                                 "multi_target_localization"));
    data = load(joinpath(data_dir, "statistical_data_localization_openloop.jld2"))
    ergodic_data = load(joinpath(data_dir, "statistical_data_localization_openloop_ergodic.jld2"));

    safe_rngnums = ergodic_data["safe_rngnums"];

    cost_tlqg = Array{Float64}[]
    tcalc_tlqg = Array{Float64}[]
    tcalc_tlqg_offline = Array{Float64}[]

    time_idx_array = [];

    num_successful_cases = 0;

    for rngnum in safe_rngnums
        println("------- RNG Seed: $(rngnum) & #Cases Tested: $(num_successful_cases) -------")
        println("4. Simulation with T-LQG Controller")
        y_history,tcalc_true_history,s_history,u_history,U_pool =
        simulate_main("tlqg",rng_num=rngnum, animate=false, plotfig=false, verbose=false);
        push!(cost_tlqg,getMetricHistory(s_history));
        push!(tcalc_tlqg,tcalc_true_history);
        if isempty(time_idx_array)
            push!(time_idx_array, unique([s.t for s in s_history]))
        end

        num_successful_cases += 1
    end

    println("5. Simulation with T-LQG Controller (Offline, for Computation Time Statistics Only)")
    rngnum = safe_rngnums[1]
    ~,tcalc_true_history,~,~,~ =
    simulate_main("tlqg_offline",rng_num=rngnum, animate=false, plotfig=false, verbose=false);
    push!(tcalc_tlqg_offline,tcalc_true_history);

    tcalc_tlqg = vcat(tcalc_tlqg...);
    tcalc_tlqg_offline = vcat(tcalc_tlqg_offline...);

    FileIO.save(joinpath(data_dir, "statistical_data_tlqg_localization_openloop.jld2"),
            "cost_tlqg", cost_tlqg, "tcalc_tlqg", tcalc_tlqg,
            "tcalc_tlqg_offline", tcalc_tlqg_offline,
            "safe_rngnums", safe_rngnums, "time_idx_array", time_idx_array);

    cost_greedy = data["cost_greedy"]
    cost_sacbp = data["cost_sacbp"]
    cost_mcts = data["cost_mcts"]
    cost_ergodic = ergodic_data["cost_ergodic"]
    tcalc_greedy = data["tcalc_greedy"]
    tcalc_sacbp = data["tcalc_sacbp"]
    tcalc_mcts = data["tcalc_mcts"]
    tcalc_ergodic = ergodic_data["tcalc_ergodic"];

    # Plot Stats Data
    cost_std_mcts = [std([cost[ii] for cost in cost_mcts]) for ii = 1:length(cost_mcts[1])]
    cost_std_sacbp = [std([cost[ii] for cost in cost_sacbp]) for ii = 1:length(cost_sacbp[1])]
    cost_std_greedy = [std([cost[ii] for cost in cost_greedy]) for ii = 1:length(cost_greedy[1])];
    cost_std_ergodic = [std([cost[ii] for cost in cost_ergodic]) for ii = 1:length(cost_ergodic[1])];
    cost_std_tlqg = [std([cost[ii] for cost in cost_tlqg]) for ii = 1:length(cost_tlqg[1])];

    MetricHistoryMCTS = mean(cost_mcts)
    plot(time_idx_array,MetricHistoryMCTS[1:end-1],label="MCTS-DPW",ribbon=cost_std_mcts,color=:lightcoral,fillalpha=0.2,linewidth=2.,linestyle=:dash)
    MetricHistoryTLQGOnline = mean(cost_tlqg);
    plot!(time_idx_array,MetricHistoryTLQGOnline[1:end-1],label="T-LQG",color=:darkgreen,ribbon=cost_std_tlqg,fillalpha=0.1,linewidth=2.,linestyle=:dot)
    MetricHistoryErgodic = mean(cost_ergodic)
    plot!(time_idx_array,MetricHistoryErgodic[1:end-1],label="Ergodic",ribbon=cost_std_ergodic,color=:darkorange,fillalpha=0.2,linewidth=2.,linestyle=:dashdot)
    MetricHistoryGreedy = mean(cost_greedy)
    plot!(time_idx_array,MetricHistoryGreedy[1:end-1],xlabel="Time [s]",ribbon=cost_std_greedy,fillalpha=0.2,ylabel="Worst Entropy Value",label="Greedy",color=:dimgrey,linewidth=2.,linestyle=:dot)
    MetricHistorySACBP = mean(cost_sacbp)
    plot!(time_idx_array,MetricHistorySACBP[1:end-1],label="SACBP (Ours)",ribbon=cost_std_sacbp,fillalpha=0.2,color=:darkblue,linewidth=2.,size=(400,400))
    savefig(joinpath(data_dir, "tracking_costs_with_tlqg.pdf"))

    # Plot Computation Time Comparison
    tcalcTrueArray = [tcalc_greedy,tcalc_mcts,tcalc_tlqg,tcalc_ergodic,tcalc_sacbp]
    meanTcalcTrueArray = mean.(tcalcTrueArray);
    stdTcalcTrueArray = std.(tcalcTrueArray);
    bar([1,2,3,4,5],meanTcalcTrueArray,yerror=stdTcalcTrueArray,linecolor=:black,marker=stroke(1.5,:black, :dash),c=[:dimgrey,:lightcoral,:palegreen3,:darkorange,:royalblue3],label="",size=(500,400));
    xticks!([1,2,3,4,5],["Greedy","MCTS-DPW","T-LQG","Ergodic","SACBP (Ours)"]);
    plot!([0.5,5.5],[3/4*0.2,3/4*0.2],width=2,label="Targeted",color=:darkred,ylim=(-0.03,0.8))
    ylabel!("Online Replanning Time [s]")
    savefig(joinpath(data_dir, "tracking_time_with_tlqg.pdf"))
end

if abspath(PROGRAM_FILE) == @__FILE__
    evaluate();
end


text:
using Gaugefields
using Test
using Random
import Wilsonloop:loops_staple

const eps = 1e-1

@testset "gradientflow_general" begin
    println("gradientflow with general action")
    include("gradientflow_general.jl")
end

@testset "gradientflow nowing" begin
    println("gradientflow nowing")
    include("gradientflow_test_nowing.jl")
end

@testset "gradientflow" begin
    println("gradientflow")
    include("gradientflow_test.jl")
end

@testset "HMC nowing" begin
    println("HMC nowing")
    include("HMC_test_nowing.jl")
end

@testset "HMC" begin
    println("HMC")
    include("HMC_test.jl")
end

@testset "heatbath" begin
    println("heatbath")
    include("heatbathtest.jl")
end

@testset "heatbath nowing" begin
    println("heatbath nowing")
    include("heatbathtest_nowing.jl")
end

@testset "heatbath with plaq and rect actions" begin
    println("heatbath with plaq and rect actions")
    include("heatbathtest_general.jl")
end

@testset "ScalarNN" begin
    println("Scalar neural networks")
    include("scalarnn.jl")
end

@testset "Initialization" begin
    println("Initialization")
    include("init.jl")
end

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


text:
include("centrality.jl")

# Generate 3-uniform, r-petal sunflower
function sunflower3(r::Int64)
    I1, I2, I3 = Int64[], Int64[], Int64[]
    core = 1
    curr_ind = 2
    for j in 1:r  # each petal
        push!(I1, core)
        push!(I2, curr_ind); curr_ind += 1
        push!(I3, curr_ind); curr_ind += 1
    end
    return SymTensor3(I1, I2, I3, ones(Float64, length(I1)), curr_ind - 1)
end

# Generate 4-uniform, r-petal sunflower
function sunflower4(r::Int64)
    I1, I2, I3, I4 = Int64[], Int64[], Int64[], Int64[]
    core = 1
    curr_ind = 2
    for j in 1:r  # each petal
        push!(I1, core)
        push!(I2, curr_ind); curr_ind += 1
        push!(I3, curr_ind); curr_ind += 1
        push!(I4, curr_ind); curr_ind += 1
    end
    return SymTensor4(I1, I2, I3, I4, ones(Float64, length(I1)), curr_ind - 1)
end

# Generate 5-uniform, r-petal sunflower
function sunflower5(r::Int64)
    I1, I2, I3, I4, I5 = Int64[], Int64[], Int64[], Int64[], Int64[]
    core = 1
    curr_ind = 2
    for j in 1:r  # each petal
        push!(I1, core)
        push!(I2, curr_ind); curr_ind += 1
        push!(I3, curr_ind); curr_ind += 1
        push!(I4, curr_ind); curr_ind += 1
        push!(I5, curr_ind); curr_ind += 1
    end
    return SymTensor5(I1, I2, I3, I4, I5, ones(Float64, length(I1)), curr_ind - 1)
end

# Test of ratios
function sf_test1()
    function ratios(m::Int64, r::Int64)
        cec_ratio = 2 * r * (m - 1) / (sqrt(m^2 + 4 * (m - 1) * (r - 1)) + m - 2)
        zec_ratio = sqrt(r)
        hec_ratio = r^(1/m)
        return (cec_ratio, zec_ratio, hec_ratio)
    end

    function get_sunflower(m::Int64, r::Int64)
        if m == 3; return sunflower3(r); end
        if m == 4; return sunflower4(r); end
        if m == 5; return sunflower5(r); end
    end

    for m in 3:5, r in 2:8
        println("m = $m; r = $r")
        S = get_sunflower(m, r)
        cec_ratio, zec_ratio, hec_ratio = ratios(m, r)
        cec_v = CEC(S)[1]
        zec_v = ZEC(S)[1]
        hec_v = HEC(S)[1]
        cec_normalized = cec_v[1] ./ cec_v
        zec_normalized = zec_v[1] ./ zec_v
        hec_normalized = hec_v[1] ./ hec_v
        println("cec ratio = $(cec_ratio)")
        println("cec vec = $(cec_normalized)")
        println("zec ratio = $(zec_ratio)")
        println("zec vec = $(zec_normalized)")
        println("hec ratio = $(hec_ratio)")
        println("hec vec = $(hec_normalized)")
        println("-----")
    end
end

function sf_test2()
    r = 4
    S = sunflower3(r)
    c = zeros(2 * r + 1)
    cP = [1.0, 2.0, 3.0, 4.0]
    c[2:3] .= cP[1]
    c[4:5] .= cP[2]
    c[6:7] .= cP[3]
    c[8:9] .= cP[4]
    c[1] = sqrt(sum(cP .^ 2))
    rat = apply(S, c) ./ c
    println("Elements should be same: $(rat)")

    S = sunflower4(r)
    c = zeros(3 * r + 1)
    cP = [1.0, 2.0, 3.0, 4.0]
    c[2:4]   .= cP[1]
    c[5:7]   .= cP[2]
    c[8:10]  .= cP[3]
    c[11:13] .= cP[4]
    c[1] = sqrt(sum(cP .^ 2))
    rat = apply(S, c) ./ c
    println("Elements should be different: $(rat)")
end
;

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
export kcore

function kcore!(s::Stinger, labels::Array{Int64}, counts::Array{Int64})
    kout = Int64(0)
    assert(length(labels)==length(counts))
    ccall(
        dlsym(stinger_alg_lib, "kcore_find"),
        Void,
        (
            Ptr{Void},
            Ref{Int64},
            Ref{Int64},
            Int64,
            Ref{Int64}
        ),
        s,
        labels,
        counts,
        length(counts),
        kout
    )
end

"""
    kcore(s::Stinger, nv::Int64)

Finds the kcore of the graph. Uses the `kcore_find` function exposed by the C
STINGER library.
Returns the labels and the counts as two `Vector{Int64}`s.
"""
function kcore(s::Stinger, nv::Int64)
    labels = zeros(Int64, nv)
    counts = zeros(Int64, nv)
    kcore!(s, labels, counts)
    return labels, counts
end

"""
    kcore(s::Stinger)

This version is slower as it makes a call to `stinger_max_active_vertex`,
which is sequential and runs through every vertex in the graph. If you know
the maximum number of active vertices, call `kcore(s, nv)` which is faster.
"""
function kcore(s::Stinger)
    nv = get_nv(s)
    kcore(s, nv)
end


text:
__precompile__()

module Currencies

# Exports
export AbstractMonetary, Monetary
export currency, decimals, majorunit, @usingcurrencies
export currencyinfo, iso4217num, iso4217alpha, shortsymbol, longsymbol
export valuate, ExchangeRateTable, ecbrates
export Basket, StaticBasket, DynamicBasket
export simplefv, compoundfv
export newcurrency!, @usingcustomcurrency
export format

# DeclarativeFormatting (not specific to Currencies; under development)
include("DeclarativeFormatting/DeclarativeFormatting.jl")

# Currency data
include("CurrencyData/CurrencyData.jl")

# Core features
include("CurrenciesBase/CurrenciesBase.jl")
using .CurrenciesBase

include("Baskets/Basket.jl")
using .Baskets

include("Valuation/Valuation.jl")
using .Valuation

# Interface (display/formatting, convenience macro)
include("CurrencyFormatting/CurrencyFormatting.jl")
using .CurrencyFormatting

# Deprecations
include("deprecated.jl")

end # module


text:
using RunTests
using Base.Test

@testmodule TestModuleWithOneXPassingTest begin

@xfail function test_xpassing()
  @test true
end

end


text:
module TestChaining

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta

##############################################################################
##
## @as macro for expression chaining
## from James Porter
##
##############################################################################

function asexpand(expr)
    if isa(expr,Expr) && expr.head == :call && expr.args[1] == :|>
        return [asexpand(expr.args[2]) expr.args[3]]
    else
        return expr
    end
end

macro as(name, bindings)
    if !isa(bindings, Expr)
        error("malformed @as bindings")
    end

    if bindings.head == :block
        exprs = filter(x-> !isa(x,Expr) || x.head != :line, bindings.args)
    elseif bindings.head == :call
        exprs = asexpand(bindings)
    end

    quote
        let $([Expr(:(=),name,expr) for expr in exprs]...)
            $name
        end
    end
end

##############################################################################
##
## @> macro for expression chaining
## adapted from @> by Mike Innes in his Lazy.jl package
## https://github.com/one-more-minute/Lazy.jlJ
##
##############################################################################

thread_left(x) = thread_left(filter(x-> !isa(x,Expr) || x.head != :line, x.args)...)

function thread_left(x, expr, exprs...)
  if typeof(expr) == Symbol
    callexpr = Expr(:call, expr, x)
  elseif typeof(expr) == Expr && expr.head in [:call, :macrocall]
    callexpr = Expr(expr.head, expr.args[1], x, expr.args[2:end]...)
  elseif typeof(expr) == Expr && expr.head == :->
    callexpr = Expr(:call, expr, x)
  else
    error("Unsupported expression $expr in @>")
  end
  isempty(exprs) ? callexpr : :(@> $callexpr $(exprs...))
end

macro >(exprs...)
    esc(thread_left(exprs...))
end

##############################################################################
##
## The test...
##
##############################################################################

srand(1)
n = 100
df = DataFrame(a = rand(1:3, n),
               b = ["a","b","c","d"][rand(1:4, n)],
               x = rand(n))

x = @where(df, :a .> 2)
x = @transform(x, y = 10 * :x)
x = @by(x, :b, meanX = mean(:x), meanY = mean(:y))
x = @orderby(x, :b, -:meanX)
x = @select(x, var = :b, :meanX, :meanY)

x_as = @as _ begin
    df
    @where(_, :a .> 2)
    @transform(_, y = 10 * :x)
    @by(_, :b, meanX = mean(:x), meanY = mean(:y))
    @orderby(_, :b, -:meanX)
    @select(_, var = :b, :meanX, :meanY)
end

x_thread = @> begin
    df
    @where(:a .> 2)
    @transform(y = 10 * :x)
    @by(:b, meanX = mean(:x), meanY = mean(:y))
    @orderby(:b, -:meanX)
    @select(var = :b, :meanX, :meanY)
end

@test x == x_as
@test x == x_thread

end # module


text:
module IterUtils

include("throw.jl")

end # module


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# diet.jl
#
# Solve the classic "diet problem".
# Based on
#  http://www.gurobi.com/documentation/5.6/example-tour/diet_cpp_cpp
#############################################################################

using JuMP

function PrintSolution(status, foods, buy)
    println("RESULTS:")
    if status == :Optimal
        for i = 1:length(foods)
            println("  $(foods[i]) = $(getvalue(buy[i]))")
        end
    else
        println("  No solution")
    end
    println("")
end

function SolveDiet()

    # Nutrition guidelines
    numCategories = 4
    categories = ["calories", "protein", "fat", "sodium"]
    minNutrition = [1800, 91, 0, 0]
    maxNutrition = [2200, Inf, 65, 1779]

    # Foods
    numFoods = 9
    foods = ["hamburger", "chicken", "hot dog", "fries",
                     "macaroni", "pizza", "salad", "milk", "ice cream"]
    cost = [2.49, 2.89, 1.50, 1.89, 2.09, 1.99, 2.49, 0.89, 1.59]
    nutritionValues = [410 24 26 730;
                       420 32 10 1190;
                       560 20 32 1800;
                       380  4 19 270;
                       320 12 10 930;
                       320 15 12 820;
                       320 31 12 1230;
                       100  8 2.5 125;
                       330  8 10 180]

    # Build model
    m = Model()

    # Variables for nutrition info
    @variable(m, minNutrition[i] <= nutrition[i=1:numCategories] <= maxNutrition[i])
    # Variables for which foods to buy
    @variable(m, buy[i=1:numFoods] >= 0)

    # Objective - minimize cost
    @objective(m, Min, dot(cost, buy))

    # Nutrition constraints
    for j = 1:numCategories
        @constraint(m, sum{nutritionValues[i,j]*buy[i], i=1:numFoods} == nutrition[j])
    end

    # Solve
    println("Solving original problem...")
    status = solve(m)
    PrintSolution(status, foods, buy)

    # Limit dairy
    @constraint(m, buy[8] + buy[9] <= 6)
    println("Solving dairy-limited problem...")
    status = solve(m)
    PrintSolution(status, foods, buy)

end

SolveDiet()


text:
function add_operation!(ls_new::LoopSet, included::Vector{Int}, ls::LoopSet, op::Operation, ids::Vector{Int}, issecond::Bool)
  newid = included[identifier(op)]
  iszero(newid) || return operations(ls_new)[newid]
  vparents = Operation[]
  for opp ∈ parents(op)
    # TODO: get it so that
    # a[i] = f(a[i]) will split into one loop computing and storing f(a[i]), and the other loading from that storage if it needs it.
    if issecond && (iscompute(opp) & (!isstore(op)))
      found = false
      for opc ∈ children(opp)
        if isstore(opc) && identifier(opc) ∉ ids
          # @show opp opc op
          # replace opp with a load from opc
          parentsopc = parents(opc)
          parentsnew = length(parentsopc) > 1 ? Operation[] : NOPARENTS
          opnew = Operation(
            length(operations(ls_new)), name(opp), opc.elementbytes, instruction(:getindex), memload,
            loopdependencies(opc), reduceddependencies(opc), parentsnew, opc.ref, reducedchildren(opc)
          )
          addsetv!(ls_new.includedactualarrays, vptr(opc.ref))
          push!(operations(ls_new), opnew)
          push!(vparents, opnew)
          for i ∈ 2:length(parentsopc)
            push!(parentsnew, add_operation!(ls_new, included, ls, parentsopc[i], ids, issecond))
          end
          included[identifier(opp)] = identifier(opnew)
          found = true
          break
        end
      end
      found && continue
    end
    push!(vparents, add_operation!(ls_new, included, ls, opp, ids, issecond))
  end
  opnew = Operation(
    length(operations(ls_new)), name(op), op.elementbytes, instruction(op), op.node_type,
    loopdependencies(op), reduceddependencies(op), vparents, op.ref, reducedchildren(op)
  )
  accesses_memory(op) && addsetv!(ls_new.includedactualarrays, vptr(op.ref))
  push!(operations(ls_new), opnew)
  included[identifier(op)] = identifier(opnew)
  opnew
end

function append_if_included!(vnew, vold, included)
  for (i, v) ∈ vold
    id = included[i]
    iszero(id) || push!(vnew, (id, v))
  end
end

function split_loopset(ls::LoopSet, ids::Vector{Int}, issecond::Bool)
  ls_new = LoopSet(:LoopVectorization)
  included = zeros(Int, length(operations(ls)))
  for i ∈ ids
    add_operation!(ls_new, included, ls, operations(ls)[i], ids, issecond)
  end
  for op ∈ operations(ls_new)
    for l ∈ loopdependencies(op)
      if l ∉ ls_new.loopsymbols
        add_loop!(ls_new, getloop(ls, l))
      end
    end
    length(ls_new.loopsymbols) == length(ls.loopsymbols) && break
  end
  append_if_included!(ls_new.preamble_symsym, ls.preamble_symsym, included)
  append_if_included!(ls_new.preamble_symint, ls.preamble_symint, included)
  append_if_included!(ls_new.preamble_symfloat, ls.preamble_symfloat, included)
  append_if_included!(ls_new.preamble_zeros, ls.preamble_zeros, included)
  append_if_included!(ls_new.preamble_funcofeltypes, ls.preamble_funcofeltypes, included)
  for i ∈ ls.outer_reductions
    id = included[i]
    iszero(id) || push!(ls_new.outer_reductions, id)
  end
  # TODO: allow them to differ. E.g., non-AVX2 x86 cpus don't have efficient integer calculations
  # Therefore, it would be profitable to split for this reason.
  # However, currently the default assumption in vector width will be wrong, so we should calculate
  # it correctly (like ls.vector_width); wrong (too high) value will encourage splitting when
  # it shouldn't.
  # Current behavior is incorrect when VECWIDTH chosen does actually differ between
  # split loops and the loops are statically sized, because code gen will then assume it is correct...
  l1, l2, l3 = cache_sze(ls)
  set_hw!(ls_new, reg_size(ls), reg_count(ls), cache_lnsze(ls), l1, l2, l3)
  ls_new.vector_width = ls.vector_width
  fill_offset_memop_collection!(ls)
  # println("ls_new operations:")
  # display(ls_new.operations)
  # println()
  ls_new
end

function returned_ops(ls::LoopSet)
    ops = operations(ls)
    retops = Int[]
    for op ∈ ops
        isstore(op) && push!(retops, identifier(op))
    end
    for i ∈ ls.outer_reductions
        push!(retops, i)
    end
    retops
end

function lower_and_split_loops(ls::LoopSet, inline::Int)
  split_candidates = returned_ops(ls)
  length(split_candidates) > 1 || return lower(ls, inline)
  order_fused, unrolled_fused, tiled_fused, vectorized_fused, U_fused, T_fused, cost_fused, shouldinline_fused = choose_order_cost(ls)
  remaining_ops = Vector{Int}(undef, length(split_candidates) - 1); split_1 = Int[0];
  # for (ind,i) ∈ enumerate(split_candidates)
  for (ind,i) ∈ enumerate(split_candidates)
    split_1[1] = i
    ls_1 = split_loopset(ls, split_1, false)
    order_1, unrolled_1, tiled_1, vectorized_1, U_1, T_1, cost_1, shouldinline_1 = choose_order_cost(ls_1)
    remaining_ops[1:ind-1] .= @view(split_candidates[1:ind-1]); remaining_ops[ind:end] .= @view(split_candidates[ind+1:end])
    ls_2 = split_loopset(ls, remaining_ops, true)
    order_2, unrolled_2, tiled_2, vectorized_2, U_2, T_2, cost_2, shouldinline_2 = choose_order_cost(ls_2)
    # U_1 = T_1 = U_2 = T_2 = 2
    # return ls_1, ls_2
    # @show cost_1 + cost_2 ≤ 0.9cost_fused, (cost_1 + cost_2) / cost_fused, cost_1, cost_2, cost_fused
    if cost_1 + cost_2 ≤ 0.9cost_fused
      ls_2_lowered = if length(remaining_ops) > 1
        inline = iszero(inline) ? (shouldinline_1 % Int) : inline
        lower_and_split_loops(ls_2, inline)
      else
        doinline = inlinedecision(inline, shouldinline_1 | shouldinline_2)
        lower(ls_2, order_2, unrolled_2, tiled_2, vectorized_2, U_2, T_2, doinline)
      end
      return Expr(
        :block,
        ls.preamble,
        lower(ls_1, order_1, unrolled_1, tiled_1, vectorized_1, U_1, T_1, false),
        ls_2_lowered,
        nothing
      )
    end
    length(split_candidates) == 2 && break
  end
  doinline = inlinedecision(inline, shouldinline_fused)
  lower(ls, order_fused, unrolled_fused, tiled_fused, vectorized_fused, U_fused, T_fused, doinline)
end


text:
struct AllocationReplay{A}
    allocations::A
    step::Base.RefValue{Int}
end

AllocationReplay(record) = AllocationReplay(record.allocations, Ref(1))

Cassette.@context ReplayCtx
new_replay_ctx(record) = new_replay_ctx(AllocationReplay(record))
function new_replay_ctx(replay::AllocationReplay)
    #replay.step[] = 1
    return Cassette.disablehooks(ReplayCtx(metadata=replay))
end

@inline function next_scheduled_alloc!(replay::AllocationReplay)
    step = replay.step[] :: Int
    alloc = replay.allocations[step] :: Array
    replay.step[] = step + 1
    return alloc
end
@inline next_scheduled_alloc!(ctx::ReplayCtx) = next_scheduled_alloc!(ctx.metadata)

@inline function Cassette.overdub(
    ctx::ReplayCtx, ::Type{Array{T,N}}, ::UndefInitializer, dims
)::Array{T,N} where {T,N}
    scheduled = next_scheduled_alloc!(ctx) :: Array{T,N}

    # Commented out until we can workout how to do this without allocations on the happy path
    # It seems like having any branch here makes it allocate
    # TODO: reenable this
    #==
    if  typeof(scheduled) !== Array{T,N} || size(scheduled) !== dims
        @warn "Allocation reuse failed. Indicates value dependent allocations." step=ctx.metadata.step[] expected_T=eltype(scheduled) actual_T=T expected_size=size actual_size=dims
        # Fallback to just doing the allocation
        return Array{T,N}(undef, dims)
    end
    ==#

    return scheduled
end

function avoid_allocations(record, f, args...; kwargs...)
    ctx = new_replay_ctx(record)
    return Cassette.overdub(ctx, f, args...; kwargs...)
end


text:
using TestImages, ImageMagick

is2d_img(s) = ! any(occursin.(["_3d_", "series","stack"], lowercase(s)))

function get2d_images()
    rf = TestImages.remotefiles
    rf = filter(is2d_img, rf)
    return (;img_list=rf,)
end

function send_testimage(;img_name)
    img = testimage(img_name)
    img_type = string(typeof(img))
    return (;img_type, bigarrs=(;img))
end

fns = (;get2d_images, send_testimage, )


text:
y = [1, 2, 3, 4, 5, 6]
plot(y, mcol="red")
y[3] = 0


text:
type ComplementarityType
    lb::Float64
    var::JuMP.Variable
    ub::Float64
    F::JuMP.NonlinearExpression
    lin_idx::Int
end

function MCPModel()
    m = Model()
    m.ext[:MCP] = Array(ComplementarityType,0)
    return m
end

function getMCPData(m::Model)
    if haskey(m.ext, :MCP)
        return m.ext[:MCP]::Array
    else
        error("The 'getMCPData' function is only for MCP models as in ComplementarityType.jl")
    end
end

# The most basic one.
function complements(m::Model, F::JuMP.NonlinearExpression, var::JuMP.Variable)
    lb = getlowerbound(var)
    ub = getupperbound(var)
    new_dimension = ComplementarityType(lb, var, ub, F, linearindex(var))
    mcp_data = getMCPData(m)
    push!(mcp_data, new_dimension)
end
function complements(m::Model, var::JuMP.Variable, F::JuMP.NonlinearExpression)
    complements(m, F, var)
end

function complements(m::Model, F::Array{JuMP.NonlinearExpression}, var::Array{JuMP.Variable})
    vars = collect(var)
    Fs = collect(F)

    @assert length(vars) == length(Fs)

    for i in 1:length(vars)
        complements(m, Fs[i], vars[i])
    end
end
function complements(m::Model, var::Array{JuMP.Variable}, F::Array{JuMP.NonlinearExpression}, )
    complements(m::Model, F, var)
end

function complements(m::Model,  arr1::JuMP.JuMPArray, arr2::JuMP.JuMPArray)
    arr1s = collect(arr1.innerArray)
    arr2s = collect(arr2.innerArray)

    @assert length(arr1s) == length(arr2s)

    for i in 1:length(arr1s)
        complements(m, arr1s[i], arr2s[i])
    end
end

function getBoundsLinearIndex(mcp_data)
    n = length(mcp_data)
    lb = zeros(n)
    ub = ones(n)
    for i in 1:n
        lb[linearindex(mcp_data[i].var)] = mcp_data[i].lb
        ub[linearindex(mcp_data[i].var)] = mcp_data[i].ub
    end
    return lb, ub
end

# Placeholder for multiple methods in the future
function solveMCP(m::Model; solver=:PATH, method=:trust_region)
    if solver == :PATH
        return _solve_path(m)
    elseif solver == :NLsolve
        return _solve_nlsolve(m, method=method)
    end
end

function sortMCPDataperm(obj::Array{ComplementarityType,1})
    n = length(obj)
    ref = Array{Int}(n)
    for i in 1:n
        ref[i] = obj[i].lin_idx
    end

    return sortperm(ref)
end

# Using PATHSolver
function _solve_path(m::Model)

    function myfunc(z)
        # z is in LindexIndex, passed from PATHSolver

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        F_val = zeros(n)
        MathProgBase.eval_g(d, F_val, z)

        # F_val also should be in LindexIndex
        # since it is the order in which constraints are added

        return F_val
    end

    function myjac(z)
        # z is in LindexIndex, passed from PATHSolver

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        I,J = MathProgBase.jac_structure(d)
        jac_val = zeros(size(J))
        MathProgBase.eval_jac_g(d, jac_val, z)

        # return matrix also should be in LindexIndex
        # since it is the order in which constraints are added

        return sparse(I, J, jac_val)
    end

    mcp_data = getMCPData(m)
    n = length(mcp_data)

    # Two Indices
    # MCP_Index: the order stored in MCPModel = array index of Array{ComplementarityType}
    # LinearIndex: the order used in JuMP / MathProgBase

    # Declaring MCP operator F as constraints
    # in order to query Jacobian using AutoDiff thru MathProgBase
    # i = LinearIndex
    # Add constraint in the order of LinearIndex
    p = sortMCPDataperm(mcp_data)
    @NLconstraint(m, dummy[i=1:n], mcp_data[p[i]].F == 0)

    # lb and ub in LinearIndex
    lb, ub = getBoundsLinearIndex(mcp_data)

    # Solve the MCP using PATHSolver
    # ALL inputs to PATHSolver must be in LinearIndex
    status, z, f = PATHSolver.solveMCP(myfunc, myjac, lb, ub)
    # z, f are in LinearIndex

    # After solving set the values in m::JuMP.Model to the solution obtained.
    for i in 1:n
        setvalue(mcp_data[i].var, z[mcp_data[i].lin_idx])
    end

    # This function has changed the content of m already.
    return status
end

function _solve_nlsolve(m::Model; method=:trust_region)

    function myfunc!(z, fvec)
        # z is in LindexIndex, passed from PATHSolver

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        # F_val = zeros(n)
        MathProgBase.eval_g(d, fvec, z)

        # F_val also should be in LindexIndex
        # since it is the order in which constraints are added

        # return fvec
    end

    function myjac!(z, fjac)
        # z is in LindexIndex, passed from PATHSolver

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Grad])
        I,J = MathProgBase.jac_structure(d)
        jac_val = zeros(size(J))
        MathProgBase.eval_jac_g(d, jac_val, z)

        # return matrix also should be in LindexIndex
        # since it is the order in which constraints are added

        sparse_fjac = sparse(I, J, jac_val)
        fjac = full(sparse_fjac)
        # return fjac
    end

    mcp_data = getMCPData(m)
    n = length(mcp_data)

    # Two Indices
    # MCP_Index: the order stored in MCPModel = array index of Array{ComplementarityType}
    # LinearIndex: the order used in JuMP / MathProgBase

    # Declaring MCP operator F as constraints
    # in order to query Jacobian using AutoDiff thru MathProgBase
    # i = LinearIndex
    # Add constraint in the order of LinearIndex
    p = sortMCPDataperm(mcp_data)
    @NLconstraint(m, dummy[i=1:n], mcp_data[p[i]].F == 0)

    # lb and ub in LinearIndex
    lb, ub = getBoundsLinearIndex(mcp_data)
    initial_z = (lb+ub) / 2
    for i = 1:length(initial_z)
        if initial_z[i] == Inf
            initial_z[i] = lb[i]
        end
    end

    # Solve the MCP using PATHSolver
    # ALL inputs to PATHSolver must be in LinearIndex
    # status, z, f = PATHSolver.solveMCP(myfunc, myjac, lb, ub)

    r = NLsolve.mcpsolve(myfunc!, myjac!, lb, ub, initial_z, method = method,
        iterations = 10_000)
    # function mcpsolve{T}(f!::Function,
    #                   g!::Function,
    #                   lower::Vector,
    #                   upper::Vector,
    #                   initial_x::Vector{T};
    #                   method::Symbol = :trust_region,
    #                   reformulation::Symbol = :smooth,
    #                   xtol::Real = zero(T),
    #                   ftol::Real = convert(T,1e-8),
    #                   iterations::Integer = 1_000,
    #                   store_trace::Bool = false,
    #                   show_trace::Bool = false,
    #                   extended_trace::Bool = false,
    #                   linesearch!::Function = Optim.backtracking_linesearch!,
    #                   factor::Real = one(T),
    #                   autoscale::Bool = true)

    # julia> fieldnames(r)
    # 12-element Array{Symbol,1}:
    #  :method
    #  :initial_x
    #  :zero
    #  :residual_norm
    #  :iterations
    #  :x_converged
    #  :xtol
    #  :f_converged
    #  :ftol
    #  :trace
    #  :f_calls
    #  :g_calls

    # r.zero is in LinearIndex

    # After solving set the values in m::JuMP.Model to the solution obtained.
    for i in 1:n
        setvalue(mcp_data[i].var, r.zero[mcp_data[i].lin_idx])
    end

    # This function has changed the content of m already.
    return r
end

# The below will be useful for creating wrapper macros
# so that users don't need to do 'using Complementarity, JuMP'
# but, just 'using Complementarity'

# Modification of deprecate_macro from JuMP.jl
macro macro_wrapper(old, new)
    oldmac = symbol(string("@",old))
    newmac = symbol(string("@",new))
    s = string(oldmac," is just a wrapper of Jump.", newmac, ".")
    if VERSION > v"0.5-"
        # backtraces are ok on 0.5
        # depwarn = :(Base.depwarn($s,$(Meta.quot(oldmac))))
        depwarn = :()
    else
        # backtraces are junk on 0.4
        # depwarn = :(Base.warn_once($s))
        depwarn = :()
    end
    @eval macro $old(args...)
        return Expr(:block, $depwarn, Expr(:macrocall, $(Meta.quot(newmac)), [esc(x) for x in args]...))
    end
    eval(Expr(:export,oldmac))
    return
end

@macro_wrapper mapping defNLExpr


text:
# # "Cosine advection-diffusion" Spatial resolution convergence test

using PyPlot

using Oceananigans.Grids

# Define a few utilities for running tests and unpacking and plotting results

include("ConvergenceTests/ConvergenceTests.jl")

run_test = ConvergenceTests.OneDimensionalCosineAdvectionDiffusion.run_test

""" Run advection-diffusion test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions...)

    # Determine save time-step
             Lx = 2π
      stop_time = 0.01
              h = Lx / maximum(resolutions)
    proposal_Δt = 1e-3 * min(h / U, h^2 / κ)

    # Adjust time-step
    stop_iteration = round(Int, stop_time / proposal_Δt)
                Δt = stop_time / stop_iteration

    # Run the tests
    results = [run_test(Nx=Nx, Δt=Δt, stop_iteration=stop_iteration, U=U, κ=κ) for Nx in resolutions]

    return results
end

#####
##### Run test
#####

Nx = 2 .^ (3:7) # N = 64 through N = 256
diffusion_results = run_convergence_test(1e-1, 0, Nx...)
advection_results = run_convergence_test(1e-6, 3, Nx...)
advection_diffusion_results = run_convergence_test(1e-2, 1, Nx...)

#####
##### Plot solution and error profile
#####

all_results = (diffusion_results, advection_results, advection_diffusion_results)
names = ("diffusion only",  "advection only",  "advection-diffusion")
linestyles = ("-", "--", ":")
specialcolors = ("xkcd:black", "xkcd:indigo", "xkcd:wine red")

# Solutions
close("all")

fig, axs = subplots(nrows=2, figsize=(12, 6), sharex=true)

ConvergenceTests.OneDimensionalUtils.plot_solutions!(axs, all_results, names, linestyles, specialcolors)

savefig("figs/cosine_advection_diffusion_solutions.png", dpi=480)

# Error profile
fig, axs = subplots()

ConvergenceTests.OneDimensionalUtils.plot_error_convergence!(axs, Nx, all_results, names)

savefig("figs/cosine_advection_diffusion_error_convergence.png", dpi=480)


text:
module AminoacidTranslations

using BioSequences

include("codon_usage.jl")

ReverseCodonTable = Dict(AA_A => (RNACodon("GCU"),RNACodon("GCC"),RNACodon("GCA"),RNACodon("GCG")), AA_R => (RNACodon("CGU"),RNACodon("GCC"),RNACodon("GCA"),RNACodon("GCG")), AA_N =>(RNACodon("AAU"),RNACodon("AAC")), AA_D => (RNACodon("GAU"),RNACodon("GAC")),AA_C =>(RNACodon("UGU"),RNACodon("UGC")), AA_Q=>(RNACodon("CAA"),RNACodon("CAG")), AA_E=>(RNACodon("GAA"),RNACodon("GAG")), AA_G=>(RNACodon("GGU"),RNACodon("GGC"),RNACodon("GGA"),RNACodon("GGG")),AA_H=>(RNACodon("CAU"),RNACodon("CAC")),AA_I=>(RNACodon("AUU"),RNACodon("AUC"),RNACodon("AUA")),AA_L=>(RNACodon("UUA"),RNACodon("UUG"),RNACodon("CUU"),RNACodon("CUC"),RNACodon("CUA"),RNACodon("CUG")), AA_K=>(RNACodon("AAA"),RNACodon("AAG")), AA_M=>(RNACodon("AUG")),AA_F=>(RNACodon("UUU"),RNACodon("UUC")),AA_P=>(RNACodon("CCU"),RNACodon("CCC"),RNACodon("CCA"),RNACodon("CCG")),AA_S=>(RNACodon("UCU"),RNACodon("UCC"),RNACodon("UCA"),RNACodon("UCG"),RNACodon("AGU"),RNACodon("AGC")),AA_T=>(RNACodon("ACU"),RNACodon("ACC"),RNACodon("ACA"),RNACodon("ACG")), AA_W=>(RNACodon("UGG")), AA_Y=>(RNACodon("UAU"),RNACodon("UAC")), AA_V=>(RNACodon("GUU"),RNACodon("GUC"),RNACodon("GUA"),RNACodon("GUG")))

ReverseCodonTableDNA = Dict(AA_A => (DNACodon("GCT"),DNACodon("GCC"),DNACodon("GCA"),DNACodon("GCG")), AA_R => (DNACodon("CGT"),DNACodon("GCC"),DNACodon("GCA"),DNACodon("GCG")), AA_N =>(DNACodon("AAT"),DNACodon("AAC")), AA_D => (DNACodon("GAT"),DNACodon("GAC")),AA_C =>(DNACodon("TGT"),DNACodon("TGC")), AA_Q=>(DNACodon("CAA"),DNACodon("CAG")), AA_E=>(DNACodon("GAA"),DNACodon("GAG")), AA_G=>(DNACodon("GGT"),DNACodon("GGC"),DNACodon("GGA"),DNACodon("GGG")),AA_H=>(DNACodon("CAT"),DNACodon("CAC")),AA_I=>(DNACodon("ATT"),DNACodon("ATC"),DNACodon("ATA")),AA_L=>(DNACodon("TTA"),DNACodon("TTG"),DNACodon("CTT"),DNACodon("CTC"),DNACodon("CTA"),DNACodon("CTG")), AA_K=>(DNACodon("AAA"),DNACodon("AAG")), AA_M=>(DNACodon("ATG")),AA_F=>(DNACodon("TTT"),DNACodon("TTC")),AA_P=>(DNACodon("CCT"),DNACodon("CCC"),DNACodon("CCA"),DNACodon("CCG")),AA_S=>(DNACodon("TCT"),DNACodon("TCC"),DNACodon("TCA"),DNACodon("TCG"),DNACodon("AGT"),DNACodon("AGC")),AA_T=>(DNACodon("ACT"),DNACodon("ACC"),DNACodon("ACA"),DNACodon("ACG")), AA_W=>(DNACodon("TGG")), AA_Y=>(DNACodon("TAT"),DNACodon("TAC")), AA_V=>(DNACodon("GTT"),DNACodon("GTC"),DNACodon("GTA"),DNACodon("GTG")))

function translate(aminoacid; to=RNA)
	if to==DNA
		ReverseCodonTableDNA[aminoacid]
       else
		ReverseCodonTable[aminoacid]
       end
end

end # module


text:
using Graphs
using SimpleValueGraphs
using Test

using Base.Iterators: product

include("testutils.jl")

tests = Dict{String, String}()

function register_test(path::String)

    tests[path] = path
end

register_test("AbstractTuples.jl")
register_test("abstractvaluegraph.jl")
register_test("graphwrappers.jl")
register_test("integrations/SimpleGraphs.jl")
register_test("interface/constructors.jl")
register_test("interface/edges.jl")
register_test("interface/edgeval_accessors.jl")
register_test("interface/iterators.jl")
register_test("interface/modifiers.jl")
register_test("interface/vertexval_accessors.jl")
# register_test("interface/valuegraph.jl")
register_test("matrices.jl")
register_test("operators.jl")

@testset "ValueGraphs" begin

    regex = r""
    if !isempty(ARGS)
        length(ARGS) > 1 && error("At most one regex can be specified for test_args")
        regex = Regex(ARGS[1])
    end

    filtered_test_names = filter(test_name -> occursin(regex, test_name), keys(tests)) |> collect |> sort

    for (i, test_name) in enumerate(filtered_test_names)
        println("$i / $(length(filtered_test_names)) $test_name")
        @time include(tests[test_name])
    end
end # testset


text:
function ObjectiveFunction(arg0::MultivariateFunction)
    return ObjectiveFunction((MultivariateFunction,), arg0)
end

function get_objective_function(obj::ObjectiveFunction)
    return jcall(obj, "getObjectiveFunction", MultivariateFunction, ())
end


text:
@testset "frohlich_α" begin

    # NaCl ∼ 5 (Feynman 1955)
    @test frohlich_α(2.3, 5.6, 4.9e13 / (2 * π), 1.0) ≈ 5.0 atol = 0.3

    # CdTe ∼ 0.39 (Stone), 0.29 (Devreese)
    @test frohlich_α(7.1, 10.4, 5.08e12, 0.095) ≈ 0.3 atol = 0.1

    # GaAs ∼ 0.068 (Devreese)
    @test frohlich_α(10.89, 12.9, 8.46e12, 0.063) ≈ 0.068 atol = 0.01

end


text:
export ExcelSheet
mutable struct ExcelSheet
	name::String
	data::DataFrame
end

export ExcelData
mutable struct ExcelData
	sheets::Array{ExcelSheet,1}
    #charts::Array{Chart,1}
    function ExcelData()
        return new(Array{ExcelSheet}(undef,0))
    end
    function ExcelData(a)
        return  new(a)
    end
end

function create_custom_dict(df::DataFrame)
	header=propertynames(df)
	d=Dict{AbstractString,Array{Any,1}}()
	for i=1:length(header)
		d[string(header[i])]=df[!,i]
	end
	return d
end

export writeExcel
function writeExcel(table,xlFile::T,write_header::Bool,write_index::Bool, printfilename::Bool = false) where {T <: AbstractString}
    #convert table to ExcelData
    table
    local df
    if isa(table,DataFrame)
        df = table
    else
        df=DataFrame(table,:auto)
    end
    xlsheet1=ExcelSheet("sheet1",df)
    xlData=ExcelData([xlsheet1])
    return writeExcel(xlData,xlFile,write_header,write_index, printfilename = printfilename)
end

function writeExcel(excelData::ExcelData,xlFile::T,write_header::Bool,write_index::Bool; printfilename::Bool = false) where {T <: AbstractString}
	#writing an Excel file seems very slow if the file already exists!
	isfile(xlFile)&&rm(xlFile)

	writer = writeDFtoExcel(excelData,xlFile,0,0,write_header,write_index)
	#workbook = writer.book
    #=
		#Plot charts
		for c in excelData.charts
			sheetWhereChartIsLocated = c.sheet
			worksheet = writer.sheets[sheetWhereChartIsLocated]
			addChartToWorkbook!(workbook,worksheet,c.chartDict,c.location);
		end
    =#

	writer.close()
	if printfilename
        println(xlFile)
    end
	return nothing
end

function writeDFtoExcel(excelData::ExcelData,existingFile::T,row::Int,col::Int,write_header::Bool,write_index::Bool) where {T <: AbstractString}
#http://search.cpan.org/~jmcnamara/Excel-Writer-XLSX/lib/Excel/Writer/XLSX.pm
	@assert min(row,col)>=0
    writer=pyModPandas.ExcelWriter(existingFile, engine = "xlsxwriter")

	for xlSheet in excelData.sheets
		df=xlSheet.data
	    sheet=xlSheet.name
		#create python dataframe
		    dataDict = create_custom_dict(df)
		pyDF=PyCall.pycall(pyModPandas.DataFrame, PyCall.PyObject, dataDict,columns=propertynames(df))
		PyCall.pycall(pyDF."to_excel",PyCall.PyAny,writer, header=write_header,index=write_index, sheet_name = sheet,startrow=row, startcol=col, encoding="utf-8")  #index=false suppress the rowcount
	end
	return writer
end


text:
module test_stdlib_profile

using Profile

f() = 0

@profile f()

# statprofilehtml()

end # module test_stdlib_profile