text:
# In this tutorial, we will learn
#    - How to solve a simple nonlinear PDE in Gridap
#    - How to define the weak residual and its Jacobian
#    - How to setup and use a nonlinear solver
#    - How to define new boundaries from a given discrete model
#
# ## Problem statement
#
# The goal of this tutorial is to solve a nonlinear PDE in Gridap. For the sake of simplicity, we consider the $p$-Laplacian equation as the model problem. Specifically, the PDE  we want to solve is: find the scalar-field $u$ such that
#
# ```math
# \left\lbrace
# \begin{aligned}
# -\nabla \cdot \left( |\nabla u|^{p-2} \ \nabla u \right) = f\ &\text{in}\ \Omega,\\
# u = 0 \ &\text{on} \ \Gamma_0,\\
# u = g \ &\text{on} \ \Gamma_g,\\
# \left( |\nabla u|^{p-2}\ \nabla u \right)\cdot n = 0 \ &\text{on} \ \Gamma_{\rm N},
# \end{aligned}
# \right.
# ```
# with $p>2$.
# The computational domain $\Omega$ is the one depicted in next figure, which is the same as in the first tutorial. However, we slightly change the boundary conditions here. We impose homogeneous Dirichlet and homogeneous Neumann boundary conditions on $\Gamma_0$ and $\Gamma_{\rm N}$  respectively, and in-homogeneous Dirichlet conditions on $\Gamma_g$. The Dirichlet boundaries $\Gamma_0$ and $\Gamma_g$ are defined as the closure of the green and blue surfaces in next figure respectively, whereas the Neumann boundary is $\Gamma_{\rm N}\doteq\partial\Omega \setminus (\Gamma_0\cup\Gamma_g)$. In this example, we consider the values $p=3$, $f=1$, and $g=2$.
#
# ![](../assets/p_laplacian/model.png)
#
# ## Numerical scheme
#
# We discretize the problem with conforming Lagrangian FE spaces. For this formulation, the nonlinear weak form reads: find $u\in U_g$ such that $[r(u)](v) = 0$ for all $v\in V_0$. As in previous tutorials, the space $U_g$ is the set of functions in $H^1(\Omega)$ that fulfill the Dirichlet boundary conditions, whereas $V_0$ is composed by functions in $H^1(\Omega)$ that vanish at the Dirichlet boundary. The weak residual $r(u)$ evaluated at a function  $u\in U_g$ is the linear form defined as
#
# ```math
# [r(u)](v) \doteq \int_\Omega \nabla v \cdot \left( |\nabla u|^{p-2}\ \nabla u \right) \ {\rm d}\Omega - \int_\Omega v\ f \ {\rm d}\Omega.
# ```
#
# In order to solve this nonlinear weak equation, we consider a Newton-Raphson method, which is associated with a linearization of the problem in an arbitrary direction $\delta u\in V_0$, namely $[r(u+\delta u)](v)\approx [r(u)](v) + [j(u)](\delta u,v)$. In previous formula,  $j(u)$ is the Jacobian evaluated at $u\in U_g$, which is the bilinear form
#
# ```math
# [j(u)](\delta u,v) = \int_\Omega \nabla v \cdot \left( |\nabla u|^{p-2}\ \nabla \delta u \right) \ {\rm d}\Omega + (p-2) \int_\Omega \nabla v \cdot \left(  |\nabla u|^{p-4} (\nabla u \cdot \nabla \delta u) \nabla u  \right) \ {\rm d}\Omega.
# ```
#
# Note that the solution of this nonlinear PDE with a Newton-Raphson method, will require to discretize both the residual $r$ and the Jacobian $j$. In Gridap, this is done by following an approach similar to the one already shown in previous tutorials for discretizing the bilinear and linear forms associated with a linear FE problem. The specific details are discussed now.
#
# ## Discrete model
#
# As in previous tutorials, the first step to solve the PDE is to load a discretization of the computational domain. In this case, we load the model from the same file as in the first tutorial

using Gridap
model = DiscreteModelFromFile("../models/model.json")

# As stated before, we want to impose Dirichlet boundary conditions on $\Gamma_0$ and $\Gamma_g$,  but none of these boundaries is identified in the model. E.g., you can easily see by writing the model in vtk format

writevtk(model,"model")

# and by opening the file `"model_0"` in Paraview that the boundary identified as `"sides"` only includes the vertices in the interior of $\Gamma_0$, but here we want to impose Dirichlet boundary conditions in the closure of $\Gamma_0$, i.e., also on the vertices on the contour of $\Gamma_0$. Fortunately, the objects on the contour of $\Gamma_0$ are identified  with the tag `"sides_c"` (see next figure). Thus, the Dirichlet boundary $\Gamma_0$ can be built as the union of the objects identified as `"sides"` and `"sides_c"`.
#
# ![](../assets/p_laplacian/sides_c.png)
#
# Gridap provides a convenient way to create new object identifiers (referred to as "tags") from existing ones. First, we need to extract from the model, the object that holds the information about the boundary identifiers (referred to as `FaceLabeling`):

labels = get_face_labeling(model)

# Then, we can add new identifiers (aka "tags") to it. In the next line, we create a new tag called `"diri0"` as the union of the objects identified as `"sides"` and `"sides_c"`, which is precisely what we need to represent the closure of the Dirichlet boundary $\Gamma_0$.

add_tag_from_tags!(labels,"diri0",["sides", "sides_c"])

# We follow the same approach to build a new identifier for the closure of the Dirichlet boundary $\Gamma_g$. In this case, the boundary is expressed as the union of the objects identified with the tags `"circle"`, `"circle_c"`, `"triangle"`, `"triangle_c"`, `"square"`, `"square_c"`. Thus, we create a new tag for  $\Gamma_g$, called `"dirig"` simply as follows:

add_tag_from_tags!(labels,"dirig",
  ["circle","circle_c", "triangle", "triangle_c", "square", "square_c"])

# ## FE Space
#
# Now, we can build the FE space by using the newly defined boundary tags.

reffe = ReferenceFE(lagrangian,Float64,1)
V0 = TestFESpace(model,reffe,conformity=:H1,labels=labels,dirichlet_tags=["diri0", "dirig"])

# The construction of this space is essentially the same as in the first tutorial (we build a continuous scalar-valued Lagrangian interpolation of first order). However, we also pass here the `labels` object (that contains the newly created boundary tags). From this FE space, we define the trial FE spaces

g = 1
Ug = TrialFESpace(V0,[0,g])

# ## Nonlinear FE problem
#
# At this point, we are ready to build the nonlinear FE problem. To this end, we need to define the weak residual and also its corresponding Jacobian. This is done following a similar procedure to the one considered in previous tutorials to define the bilinear and linear forms associated with linear FE problems. We first need to define the usual objects for numerical integration:

degree=2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

#
# On the one hand, the weak residual is built as follows

using LinearAlgebra: norm
const p = 3
flux(∇u) = norm(∇u)^(p-2) * ∇u
f(x) = 1
res(u,v) = ∫( ∇(v)⊙(flux∘∇(u)) - v*f )*dΩ

# Function `res` is the one representing the integrand of the weak residual $[r(u)](v)$. The first argument of function `res` stands for the function $u\in U_g$, where the residual is evaluated, and the second argument stands for a generic test function $v\in V_0$.
#
# On the other hand,  we (optionally) implement a function `jac` representing the Jacobian.
dflux(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*(∇u⊙∇du)*∇u+norm(∇u)^(p-2)*∇du
jac(u,du,v) = ∫( ∇(v)⊙(dflux∘(∇(du),∇(u))) )*dΩ

# The first argument of function `jac` stands for function $u\in U_g$, where the Jacobian is evaluated. The second argument is a test function $v\in V_0$, and the third argument represents an arbitrary direction $\delta u \in V_0$.
#
# We finally construct the nonlinear FE problem

op = FEOperator(res,jac,Ug,V0)

# Here, we have constructed an instance of `FEOperator`, which is the type that represents a general nonlinear FE problem in Gridap. The constructor takes the functions representing the weak residual and Jacobian, and the test and trial spaces. If only the function for the residual is provided, the Jacobian is computed internally with automatic differentiation.

#
# ## Nonlinear solver phase
#
# We have already built the nonlinear FE problem. Now, the remaining step is to solve it. In Gridap, nonlinear (and also linear) FE problems can be solved with instances of the type `FESolver`.
#
# We construct an instance of `FESolver` as follows:

using LineSearches: BackTracking
nls = NLSolver(
  show_trace=true, method=:newton, linesearch=BackTracking())
solver = FESolver(nls)

# Note that the `NLSolver` function used above internally calls the `nlsolve` function of the [NLsolve](https://github.com/JuliaNLSolvers/NLsolve.jl) package with the provided key-word arguments. Thus, one can use any of the nonlinear methods available via the function `nlsolve` to solve the nonlinear FE problem. Here, we have selected a Newton-Raphson method with a back-tracking line-search from the [LineSearches](https://github.com/JuliaNLSolvers/LineSearches.jl) package.
#
# We are finally in place to solve the nonlinear FE problem. The initial guess is a `FEFunction`, which we build from a vector of random (free) nodal values:

import Random
Random.seed!(1234)
x = rand(Float64,num_free_dofs(Ug))
uh0 = FEFunction(Ug,x)
uh, = solve!(uh0,solver,op)

# We finish this tutorial by writing the computed solution for visualization (see next figure).

writevtk(Ω,"results",cellfields=["uh"=>uh])

# ![](../assets/p_laplacian/sol-plap.png)
#


text:
include("../examples/vrptw_example.jl")
include("../examples/vrptw_example_dists.jl")

ext: jl


text:
module BetterMGF
using ..Concepts
import LinearAlgebra

struct MGF{T<:Any}  <: AbstractMGF
    OPTION_LOG_SCALE::Optional{Bool}

    function MGF{T}() where T<:Any
        return new{T}()
    end

    function MGF{T}(object::T;logscale::Optional{Bool} = nothing) where T<:Any
        if isnothing(logscale)
            return new{T}()
        end
        return new{T}(logscale)
    end
end

const MGF() = MGF{Any}()

const MGF(object::Union{T, Symbol}; logscale::Optional{Bool} = nothing) where T<:Any =
    !isa(object,Symbol) ? MGF{T}(object, logscale = logscale) : MGF(type_conversion(object), logscale=logscale)

struct SampleMGF <: AbstractMGF
    OPTION_LOG_SCALE::Optional{Bool}

    function SampleMGF(;logscale::Optional{Bool} = nothing)
        if isnothing(logscale)
            return new()
        end
        return new(logscale)
    end
end

export MGF,
    SampleMGF

# @overload
# function Base.convert(::Type{FrequencyDomainObjects},x::Symbol)
#     if x == :MGF
#         return MGF
#     end
# end

@overload
function Concepts.type_conversion(::Type{FrequencyDomainObjects}, x::Symbol)
    if x == :MGF
        return MGF
    end
end

@overload
function Concepts.evaluate(object::MGF{AbstractPoisson},t::AutoboxedArray{Real};
                           λ::S) where S<:Real
    if λ <= 0
        throw(DomainError("λ should be a positive number"))
    end
    if !isnothing(object.OPTION_LOG_SCALE) && object.OPTION_LOG_SCALE == true
        @warn "[evaluate(MGF{Poisson}): using log-scale]"
        return λ .* (exp.(t) .- 1)
    end
    return exp.(λ .* (exp.(t) .- 1))
end

@overload
function Concepts.evaluate(object::MGF{AbstractBernoulli},t::AutoboxedArray{Real};
                           p::S) where S<:Real
    if p < 0 || p >1
        throw(DomainError("p should be a positive number in [0,1]"))
    end
    if !isnothing(object.OPTION_LOG_SCALE) && object.OPTION_LOG_SCALE == true
        @warn "[evaluate(MGF{Bernoulli}): using log-scale]"
        return log.((1 .- p) .+ p .* exp.(t))
    end
    return (1 .- p) .+ p .* exp.(t)
end

@overload
function Concepts.evaluate(object::MGF{AbstractGamma},t::AutoboxedArray{Real};
                           α::S,θ::S) where S<:Real
    if α <= 1e-5 || θ <= 1e-5
        throw(DomainError("α and θ should be a positive number"))
    end
    if !isnothing(object.OPTION_LOG_SCALE) && object.OPTION_LOG_SCALE == true
        @warn "[evaluate(MGF{Gamma}): using log-scale]"
        return (-α) .* log.(1 .- t .* θ)
    end
    return (1 .- t .* θ) .^ (-α)
end

@overload
function Concepts.evaluate(object::MGF{AbstractGaussian},t::AutoboxedArray{Real};
                           μ::S = 0,σ::S=1) where S<:Real
    if abs(σ) <= 1e-4
        throw(DomainError("σ should be non-zero. Otherwise, it will be degenerate."))
    end
    if !isnothing(object.OPTION_LOG_SCALE) && object.OPTION_LOG_SCALE == true
        @warn "[evaluate(MGF{Gaussian}): using log-scale]"
        return t .* μ .+ σ^2 / 2 .* t.^2
    end
    return exp.(t .* μ .+ σ^2 / 2 .* t.^2)
end

import Distributions
@overload
function Concepts.evaluate(object::MGF{AbstractNegativeBinomial},t::AutoboxedArray{Real};
                           r::S = 0,p::S=1) where S<:Real
    if r <= 1e-5 || p <= 1e-5 || p >=1
        throw(DomainError("[Negative Binomial]: r should be > 0, p should be in (0,1)"))
    end
    if !isnothing(object.OPTION_LOG_SCALE) && object.OPTION_LOG_SCALE == true
        @warn "[evaluate(MGF{NegativeBinomial})]: using log-scale"
        return (r .* (log.(1 .- p) .+ t)) .- (r .* log.(1 .- p .* exp.(t)))
    end
    return ((1 .- p) .* exp.(t)).^r ./ (1 .- p .* exp.(t)).^r
end

@overload
function Concepts.evaluate(object::SampleMGF,t::VecOrMat{T};data::VecOrMatOf{Real},order::Integer = 8,data_layout::Symbol=:flatten) where T <: Real
    if isa(t,Matrix) || isa(data,Matrix)
        #        throw(DomainError("Unimplemented"))
        if data_layout == :flatten
            @info "MGF of the matrix is evaluated by flattening the data matrix into a vector"
            return dp_get_sample_basis(t,order) * dp_get_sample_mean(data[:],order)
        end
        if data_layout == :bycol
            throw(DomainError("Unimplemented"))
        end
        if data_layout == :byrow
            throw(DomainError("Unimplemented"))
        end
    end
    return dp_get_sample_basis(t,order) * dp_get_sample_mean(data,order)
end

@private
function dp_get_denominator(order::Integer = 10)
    denom = ones(order+1)
    for i = 2:order+1
        denom[i] = denom[i-1]*(i-1)
    end
    return 1 ./ denom
end

@private
function dp_get_numerator(t::Array{S,1},order::Integer = 10) where S<:Real
    n = length(t)
    T = ones(n,order+1)
    for i = 2: order +1
        T[:,i] = T[:,i-1] .* t
    end
    return T
end

@private
function dp_get_sample_basis(t::Array{S,1}, order::Integer = 10) where S<:Real
    num = dp_get_numerator(t,order)
    denom = dp_get_denominator(order)
    for i = 1:size(num)[1]
       num[i,:] = num[i,:] .* denom
    end
    return num
end

@private
function dp_get_sample_mean(data::Array{S,1}, order::Integer = 10) where S<: Real
    ret = ones(order + 1)
    copy_of_data = deepcopy(data)
    n = length(copy_of_data)
    for i = 2:order + 1
        ret[i] = sum(copy_of_data)/n
        copy_of_data = copy_of_data .* data
    end
    return ret
end

end


text:
const K_CART_TYPE{T} = Quantity{T,Unitful.𝐋^-1,Unitful.FreeUnits{(Ang^-1,),Unitful.𝐋^-1,nothing}}

phases(kpoints::Vector{<:Vec3}, R::Vec3) = exp.(-2im * π .* dot.(kpoints, (R,)))

uniform_shifted_kgrid(::Type{T}, nkx::Integer, nky::Integer, nkz::Integer) where {T} =
	reshape([Vec3{T}(kx, ky, kz) for kx = 0.5/nkx:1/nkx:1, ky = 0.5/nky:1/nky:1, kz = 0.5/nkz:1/nkz:1], nkx*nky*nkz)

uniform_shifted_kgrid(nkx::Integer, nky::Integer, nkz::Integer) = uniform_shifted_kgrid(Float64, nkx, nky, nkz)

uniform_kgrid(nkx::Integer, nky::Integer, nkz::Integer) =
    reshape([Vec3{Float64}(kx, ky, kz) for kz in range(0, 1-1/nkz, length=nkz), ky in range(0, 1-1/nky, length=nky), kx in range(0, 1-1/nkx, length=nkx)], nkx*nky*nkz)

abstract type AbstractKGrid{T} end

core_kgrid(x::AbstractKGrid) = x.core

k_cryst(x::AbstractKGrid)    = core_kgrid(x).k_cryst
k_cryst(x::Vec3) = x
# k_cart(x::AbstractKGrid)     = core_kgrid(x).k_cart
# phase(x::AbstractKGrid)      = core_kgrid(x).phase
Base.length(kgrid::AbstractKGrid) = length(core_kgrid(kgrid))

struct CoreKGrid{T} <: AbstractKGrid{T}
    k_cryst::Vector{Vec3{T}}
    # k_cart ::Vector{Vec3{T}} #not needed for now
    # phase  ::Vector{Complex{T}}
end

core_kgrid(x::CoreKGrid) = x

Base.length(x::CoreKGrid) = length(x.k_cryst)

abstract type KPoint{T<:AbstractFloat} end

struct KBond{T<:AbstractFloat} # All matrices/operators are in wannier gauge, i.e. block-like gauge
    k_id1   ::Int #TODO: Optimize, can we drop these?
    k_id2   ::Int
    vr      ::Vec3{K_CART_TYPE{T}}
end

@with_kw mutable struct AbInitioKPoint{T} <: KPoint{T}
    k_cryst  ::Vec3{T}
    k_cart   ::Vec3{K_CART_TYPE{T}}
    eigvals  ::Vector{T}         #original eigenvalues, in hamiltonian gauge
    H        ::Matrix{Complex{T}} # Hk in wannier gauge
    neighbors::Vector{KBond{T}} = KBond{T}[]
    overlaps ::Vector{Matrix{Complex{T}}} = Matrix{Complex{T}}[] #already in wannier gauge
    hamis    ::Vector{Matrix{Complex{T}}} = Matrix{Complex{T}}[] #Hamiltonian element between the block-like states in wannier gauge
    uHu      ::Matrix{Matrix{Complex{T}}} = Matrix{Matrix{Complex{T}}}(undef, 0, 0)
end

k_cryst(k::AbInitioKPoint) = k.k_cryst

struct AbInitioKGrid{T, SA} <: AbstractKGrid{T}
    kpoints::SA
    neighbor_weights::Vector{T} #ordered in the same way as neighbors in kpoints
end

num_states(chk, ik) = chk.ndimwin[ik]

function disentanglement_range(chk, ik)
    first_band_id = findfirst(view(chk.lwindow, :, ik))
    return first_band_id : first_band_id + num_states(chk, ik) - 1
end

function AbInitioKGrid(::Type{T},
                       eig_filename::AbstractString,
                       chk_filename::AbstractString,
                       nnkp_filename::AbstractString,
                       mmn_filename::AbstractString,
                       uHu_filename::AbstractString) where {T}
    eigenvalues = read_eig(eig_filename)
    wannier_chk_params = read_chk(chk_filename)
    v_mat = wannier_chk_params.V_matrix
    nwann = wannier_chk_params.n_wann
    kpoints = map(enumerate(wannier_chk_params.kpoints)) do (ik, k)
        @views v = v_mat[1:num_states(wannier_chk_params, ik), 1:nwann, ik]
        eigvals = eigenvalues[:, ik]
        return AbInitioKPoint{T}(k_cryst = k,
                                 k_cart  = wannier_chk_params.recip_cell * k,
                                 eigvals = eigenvalues[:, ik],
                                 H       = v' * diagm(eigvals[disentanglement_range(wannier_chk_params, ik)]) * v)

    end
    fill_k_neighbors!(kpoints, nnkp_filename, wannier_chk_params.recip_cell)
    fill_overlaps!(kpoints, mmn_filename, uHu_filename, wannier_chk_params)

    return AbInitioKGrid(StructArray(kpoints), wannier_chk_params.neighbor_weights)
end

AbInitioKGrid(eig_filename::AbstractString, chk_filename::AbstractString, nnkp_filename::AbstractString, mmn_filename::AbstractString, uHu_filename::AbstractString) =
    AbInitioKGrid(Float64, eig_filename, chk_filename, nnkp_filename, mmn_filename, uHu_filename)

function AbInitioKGrid(job::Job)
    wancalc = getfirst(x->x isa Calculation{Wannier90}, job.calculations)
    if wancalc === nothing
        error("Coulnd't find a wannier calculation in job $job.")
    end
    wname = name(wancalc)

    wan_file = ext -> begin
        files = DFC.searchdir(job, "$wname.$ext")
        if isempty(files)
            error("Couldn't find $wname.$ext in job directory $(job.local_dir).")
        end
        return files[1]
    end

    return AbInitioKGrid(wan_file("eig"), wan_file("chk"), wan_file("nnkp"), wan_file("mmn"), wan_file("uHu"))
end

n_wannier_functions(grid::AbInitioKGrid) = size(grid.kpoints.overlaps[1], 1)
n_nearest_neighbors(grid::AbInitioKGrid) = length(grid.kpoints.neighbors[1])

Base.length(grid::AbInitioKGrid) = length(grid.kpoints)


text:
"""
secant(f,x0; rtol=1.e-6, atol=1.e-12, maxit=10,
        armmax=10, armfix=false, pdata=nothing,
        printerr=true, keepsolhist=true, stagnationok=false)

C. T. Kelley, 2021

The secant method for scalar equations.

Input:\n
f: function\n
x0: initial iterate

Keyword Arguments (kwargs):\n
rtol, atol: real and absolute error tolerances\n

maxit: upper bound on number of nonlinear iterations\n

If you use secant and your initial iterate is poor, you have made
a mistake. You will get an error message.

armmax: upper bound on stepsize reductions in linesearch

armfix:\n
The default is a parabolic line search (ie false). Set to true and
the stepsize will be fixed at .5. Don't do this unless you are doing
experiments for research.

printerr:\n
I print a helpful message when the solver fails. To supress that
message set printerr to false.

keepsolhist:\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

stagnationok:\n
Set this to true if you want to disable the line search and either
observe divergence or stagnation. This is only useful for research
or writing a book.

Output:\n
A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

solution = converged result
functionval = F(solution)
history = the vector of residual norms (||F(x)||) for the iteration
stats = named tuple of the history of (ifun, ijac, iarm), the number
of functions/derivatives/steplength reductions at each iteration.
For the secant method, ijac = 0.

idid=true if the iteration succeeded and false if not.

errcode = 0 if if the iteration succeeded
        = -1 if the initial iterate satisifies the termination criteria
        = 10 if no convergence after maxit iterations
        = 1  if the line search failed

solhist:\n
This is the entire history of the iteration if you've set
keepsolhist=true\n
secant builds solhist with a function from the Tools directory. For
systems, solhist is an N x K array where N is the length of x and K
is the number of iteration + 1. So, for scalar equations (N=1), solhist
is a row vector. Hence the use of solhist' in the example below.

### Example for secant.jl

```jldoctest

julia> secout=secant(atan,1.0;maxit=6,atol=1.e-12,rtol=1.e-12);

julia> secout.history
7-element Array{Float64,1}:
 7.85398e-01
 5.18729e-01
 5.39030e-02
 4.86125e-03
 4.28860e-06
 3.37529e-11
 2.06924e-22
```
"""
function secant(
    f,
    x0;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 10,
    solver = "secant",
    armmax = 5,
    armfix = false,
    dx = 1.e-7,
    pdata = nothing,
    printerr = true,
    keepsolhist = true,
    stagnationok = false,
)
    itc = 0
    idid = true
    errcode = 0
    iline = false
    #=
     The theory does not support convergence of the secant-Armijo iteration
     and you assume a risk when you use it. The same is true for Broyden
     and any other quasi-Newton method.
    =#
    fc = 0.0
    fc = EvalF!(f, fc, x0, pdata)
    fm = fc
    xm = copy(x0)
    xm = x0 * 1.0001
    if xm == 0
        xm = 0.0001
    end
    fm = fc
    fm = EvalF!(f, fm, xm, pdata)
    newfun0 = 1
    derivative_is_old = false
    resnorm = abs(fc)
    jfact = nothing
    stagflag = stagnationok && (armmax == 0)
    (ItRules, x, n) =
        Secantinit(x0, dx, f, solver, armmax, armfix, maxit, printerr, pdata, jfact)
    #
    # Initialize the iteration statistics
    #
    newiarm = -1
    ItData = ItStats(resnorm, 2)
    newfun = 0
    newjac = 0
    newsol = x
    xt = x
    keepsolhist ? (solhist = solhistinit(n, maxit, x)) : (solhist = [])
    #
    # Fix the tolerances for convergence and define the derivative df
    # outside of the main loop for scoping.
    #
    tol = rtol * resnorm + atol
    residratio = 1
    df = 0.0
    armstop = true
    #
    # If the initial iterate satisfies the termination criteria, tell me.
    #
    toosoon = (resnorm <= tol)
    #
    # The main loop stops on convergence, too many iterations, or a
    # line search failure after a derivative evaluation.
    #
    while (resnorm > tol) && (itc < maxit) && (armstop || stagnationok)
        newfun = 0
        #
        # Extra function call at the start.
        #
        newjac = 0
        newfun = 0
        #
        df = (fc - fm) / (x - xm)
        derivative_is_old = (newjac == 0) && (solver == "newton")
        #
        # Compute the Newton direction and call the line search.
        #
        xm = x
        fm = fc
        ft = fc
        d = -fc / df
        AOUT = armijosc(xt, x, ft, fc, d, resnorm, ItRules, derivative_is_old)
        #
        # update solution/function value
        #
        xm = x
        x = AOUT.ax
        fm = fc
        fc = AOUT.afc
        #
        # If the line search fails and the derivative is current,
        # stop the iteration.
        #
        armstop = AOUT.idid || derivative_is_old
        iline = ~armstop && ~stagflag
        newiarm = AOUT.aiarm
        #
        # Keep the books.
        #
        residm = resnorm
        resnorm = AOUT.resnorm
        residratio = resnorm / residm
        updateStats!(ItData, newfun, newjac, AOUT)
        #
        itc += 1
        ~keepsolhist || (@views solhist[:, itc+1] .= x)
    end
    (idid, errcode) = NewtonOK(resnorm, iline, tol, toosoon, itc, ItRules)
    newtonout = CloseIteration(x, fc, ItData, idid, errcode, keepsolhist, solhist)
    return newtonout
end


text:
#
# Types supporting Parameters and their connections
#

abstract type ModelParameter <: MimiStruct end

mutable struct ScalarModelParameter{T} <: ModelParameter
    value::T
    is_shared::Bool

    function ScalarModelParameter{T}(value::T; is_shared::Bool = false) where T
        new(value, is_shared)
    end

    function ScalarModelParameter{T}(value::T, is_shared::Bool) where T
        new(value, is_shared)
    end

    function ScalarModelParameter{T1}(value::T2; is_shared::Bool = false) where {T1, T2}
        try
            new(T1(value), is_shared)
        catch err
            error("Failed to convert $value::$T2 to $T1")
        end
    end
end

mutable struct ArrayModelParameter{T} <: ModelParameter
    values::T
    dim_names::Vector{Symbol} # if empty, we don't have the dimensions' name information
    is_shared::Bool

    function ArrayModelParameter{T}(values::T, dims::Vector{Symbol}; is_shared::Bool = false) where T
        new(values, dims, is_shared)
    end

    function ArrayModelParameter{T}(values::T, dims::Vector{Symbol}, is_shared::Bool) where T
        new(values, dims, is_shared)
    end
end

ScalarModelParameter(value) = ScalarModelParameter{typeof(value)}(value)
ScalarModelParameter(value, is_shared) = ScalarModelParameter{typeof(value)}(value, is_shared)

Base.convert(::Type{ScalarModelParameter{T}}, value::Number) where {T} = ScalarModelParameter{T}(T(value))
Base.convert(::Type{T}, s::ScalarModelParameter{T}) where {T} = T(s.value)

ArrayModelParameter(value, dims::Vector{Symbol}) = ArrayModelParameter{typeof(value)}(value, dims)
ArrayModelParameter(value, dims::Vector{Symbol}, is_shared::Bool) = ArrayModelParameter{typeof(value)}(value, dims, is_shared)

# Allow values to be obtained from either parameter type using one method name.
value(param::ArrayModelParameter)  = param.values
value(param::ScalarModelParameter) = param.value

Base.copy(obj::ScalarModelParameter{T}) where T = ScalarModelParameter(obj.value, obj.is_shared)
Base.copy(obj::ArrayModelParameter{T}) where T = ArrayModelParameter(obj.values, obj.dim_names, obj.is_shared)

dim_names(obj::ArrayModelParameter) = obj.dim_names
dim_names(obj::ScalarModelParameter) = []

is_shared(obj::ArrayModelParameter) = obj.is_shared
is_shared(obj::ScalarModelParameter) = obj.is_shared

abstract type AbstractConnection <: MimiStruct end

struct InternalParameterConnection <: AbstractConnection
    src_comp_path::ComponentPath
    src_var_name::Symbol
    dst_comp_path::ComponentPath
    dst_par_name::Symbol
    ignoreunits::Bool
    backup::Union{Symbol, Nothing} # a Symbol identifying the model param providing backup data, or nothing
    backup_offset::Union{Int, Nothing}

    function InternalParameterConnection(src_path::ComponentPath, src_var::Symbol,
                                         dst_path::ComponentPath, dst_par::Symbol,
                                         ignoreunits::Bool, backup::Union{Symbol, Nothing}=nothing;
                                         backup_offset::Union{Int, Nothing}=nothing)
        self = new(src_path, src_var, dst_path, dst_par, ignoreunits, backup, backup_offset)
        return self
    end
end

struct ExternalParameterConnection  <: AbstractConnection
    comp_path::ComponentPath
    param_name::Symbol      # name of the parameter in the component
    model_param_name::Symbol  # name of the parameter stored in model_params
end

# Converts symbol to component path
function ExternalParameterConnection(comp_name::Symbol, param_name::Symbol, model_param_name::Symbol)
    return ExternalParameterConnection(ComponentPath(comp_name), param_name, model_param_name)
end

Base.pathof(obj::ExternalParameterConnection) = obj.comp_path
Base.nameof(obj::ExternalParameterConnection) = obj.param_name

##
## DEPRECATIONS - Should move from warning --> error --> removal
##

function Base.getproperty(epc::ExternalParameterConnection, field::Symbol)
    if field == :external_param
        @warn "ExternalParameterConnection's `external_param` field is renamed to `model_param_name`, please change code accordingly."
        field = :model_param_name
    end
    return getfield(epc, field)
end


text:
using KernelFunctions, Test, Distances
#TODO Test metric weights for ARD, test equivalency for different constructors,
# test type conversion
l = 2.0
vl = [l,l]
s = ScaleTransform(l)

## Add tests for Transformed Kernel and Scaled Kernel

## SqExponentialKernel
@testset "Exponential" begin
    @test KernelFunctions.metric(ExponentialKernel()) == Euclidean()
    @test KernelFunctions.metric(SqExponentialKernel()) == SqEuclidean()
    @test KernelFunctions.metric(GammaExponentialKernel()) == SqEuclidean()
    @test KernelFunctions.metric(GammaExponentialKernel(γ=2.0)) == SqEuclidean()
end

## MaternKernel
@testset "MaternKernel" begin
    @test KernelFunctions.metric(MaternKernel()) == Euclidean()
    @test KernelFunctions.metric(MaternKernel(ν=2.0)) == Euclidean()
    @test KernelFunctions.metric(Matern32Kernel()) == Euclidean()
    @test KernelFunctions.metric(Matern52Kernel()) == Euclidean()
end

@testset "Exponentiated" begin
    @test KernelFunctions.metric(ExponentiatedKernel()) == KernelFunctions.DotProduct()
end

@testset "Constant" begin
    @test KernelFunctions.metric(ConstantKernel()) == KernelFunctions.Delta()
    @test KernelFunctions.metric(ConstantKernel(c=2.0)) == KernelFunctions.Delta()
    @test KernelFunctions.metric(WhiteKernel()) == KernelFunctions.Delta()
    @test KernelFunctions.metric(ZeroKernel()) == KernelFunctions.Delta()
end

@testset "Polynomial" begin
    @test KernelFunctions.metric(LinearKernel()) == KernelFunctions.DotProduct()
    @test KernelFunctions.metric(LinearKernel(c=2.0)) == KernelFunctions.DotProduct()
    @test KernelFunctions.metric(PolynomialKernel()) == KernelFunctions.DotProduct()
    @test KernelFunctions.metric(PolynomialKernel(d=3.0)) == KernelFunctions.DotProduct()
    @test KernelFunctions.metric(PolynomialKernel(d=3.0,c=2.0)) == KernelFunctions.DotProduct()
end

@testset "RationalQuadratic" begin
    @test KernelFunctions.metric(RationalQuadraticKernel()) == SqEuclidean()
    @test KernelFunctions.metric(RationalQuadraticKernel(α=2.0)) == SqEuclidean()
    @test KernelFunctions.metric(GammaRationalQuadraticKernel()) == SqEuclidean()
    @test KernelFunctions.metric(GammaRationalQuadraticKernel(γ=2.0)) == SqEuclidean()
    @test KernelFunctions.metric(GammaRationalQuadraticKernel(γ=2.0,α=3.0)) == SqEuclidean()
end


text:
@doc """
  Computes the error-free addition `a + b = x + y`, where `x = fl(a + b)`, for
  general floating-point numbers `a` and `b`. The algorithm requires 6 flops;
  twice as many as `err_fast_add`. However, failed branch prediction because of
  the required comparison `|a| ≥ |b|` in the general case, cause this algorithm
  to be faster on modern pipelined architectures.

  References:

    * D.E. Knuth, The Art of Computer Programming, 3rd ed., vol. 2,
      Addison-Wesley, 1998.
    * D. Møller, Quasi double-precision in floating-point addition, BIT
      Numerical Mathematics 5, pp. 37–50, 1965.
""" ->
# TODO: Extend the definition to floating-point arrays
function err_add{T<:IEEE754}(a::T, b::T)
  x = a + b
  c = x - a
  y = (a - (x-c)) + (b-c)
  x, y
end

@doc """
  Computes the error-free addition `a + b = x + y`, where `x = fl(a + b)`, for
  floating-point numbers `a` and `b` such that the exponent of `a` is larger
  than or equal to that of `b`. If `|a| ≥ |b|`, this condition is always
  satisfied. The algorithm requires 3 flops.

  References:

    * T.J. Dekker, A floating-point technique for extending the available
      precision, Numerische Mathematik 18, pp. 224–242, 1971.
  """ ->
# TODO: Extend the definition to floating-point arrays
function err_fast_add{T<:IEEE754}(a::T, b::T)
  x = a + b
  y = b - (x-a)
  x, y
end


text:
"""
`module JAC.Math`
    ... a submodel of JAC that contains various (adopted) procedures from different libraries and sources.
"""
module Math

  using GSL
  using JAC, ..Radial

  """
  `Math.finiteDifferenceWeights(position, npoints; order=1)`
    ... computes the weights for a n-point finite-difference approximation. position=0 corresponds to the centered stencil,
        while positive and negative values refer to the non-symmetric stencils, respecitvely. The order specifies for which
        derivative the coefficients are computed, the default ist the first derivative.
  """
  function finiteDifferenceWeights(position, npoints; order=1)

      if !isodd(npoints)
          error("Number of points must be odd")
      end

      x = collect(-Int64((npoints-1)/2):Int64((npoints-1)/2))
      c = Array{Float64}(undef, order+1, npoints)

      c1 = 1
      c4 = x[1] - position
      c[1, 1] = 1

      for i = 2:npoints
          mn = min(i, order+1)
          c2 = 1
          c5 = c4
          c4 = x[i] - position
          for j = 1:i-1
              c3 = x[i] - x[j]
              c2 = c2 * c3
              if j == i - 1
                  c[2:mn, i] = c1 / c2 * ( collect(1:mn -1 ) .* c[1:mn-1, i-1] - c5 * c[2:mn, i-1])
                  c[1, i] = -c1 * c5 * c[1, i-1]/c2
              end
              c[2:mn, j] = 1 / c3 * (c4 * c[2:mn, j] - collect(1:mn-1) .* c[1:mn-1, j])
              c[1, j] = c4 * c[1, j] / c3
          end
          c1 = c2
      end

      return( c )
   end

  """
  `Math.derivative(input::Array{Float64}, i::Int64)`
    ... computes the derivative of a given input array at position i; a Float64 is returned. An equally spaced grid is
        assumed, the result must be multiplied by the inverse lattice spacing h manually. If an exponential grid is
        used, the transformation must be applied too.
  """
  function derivative(input::Array{Float64}, i::Int64) :: Float64

      n = JAC.Defaults.FINITE_DIFFERENCE_NPOINTS
      result ::Float64 = 0.

      if i <= n

        for j = 1:2*n + 1
          result += JAC.Defaults.weights[i, j] * input[j]
        end
        return result
      end

      if size(input, 1) - i < n

        for j = 1:2*n + 1
          result += JAC.Defaults.weights[2*n + 1 + i - size(input, 1), j] * input[size(input, 1) - 2*n - 1 + j]
        end

        return result
      end

      for j = i - n:i+n
        result += JAC.Defaults.weights[n + 1, j - i + n + 1] * input[j]
      end

      return result
  end

  """
  `Math.integrateFitTransform(f::Function, mtp::Int64, grid::Radial.Grid)`
    ... to integrate a (point-wise) given function over mtp grid points on an exponential grid. The transformation of the
        exponential grid is automatically performed such that no modifications to the integrand are necessary; a value::Float64
        is returned.
  """
  function integrateFitTransform(f::Function, mtp::Int64, grid::Radial.Grid)
      n      = size(JAC.Defaults.newtonCotesCoefficients, 1)
      newtonCotesCoefficientsOnGrid = copy(JAC.Defaults.newtonCotesCoefficients) * grid.h
      i0     = 2
      gamma  = 0.
      result = 0.

      while abs(f(i0)) == 0 || f(i0)/f(i0+1) < 0
          i0 += 1
      end
      gamma = log(f(i0)/f(i0+1)) / log(grid.r[i0]/grid.r[i0+1])
      result = 1/(gamma + 1) * grid.r[i0] * f(i0)

      result += f(i0) * grid.rp[i0] * newtonCotesCoefficientsOnGrid[n]
      for i = i0+1:mtp
          result += newtonCotesCoefficientsOnGrid[(i - i0) % (n - 1) + 1] * f(i) * grid.rp[i]
      end

      if (mtp - i0) % (n - 1) + 1 == 1
          result -= newtonCotesCoefficientsOnGrid[n] * f(mtp) * grid.rp[mtp]
      end

    return( result )
  end

  """
  `Math.integrateFit(f::Function, mtp::Int64, grid::Radial.Grid)`  ... to integrate a (point-wise) given function over mtp grid points
       on an exponential grid. This function assumes that the ingegrand was already transformed to the exponential grid by multiplying it with
       the derivative of the transformation function; a value::Float64 is returned.
  """
  function integrateFit(f::Function, mtp::Int64, grid::Radial.Grid)
      n      = size(JAC.Defaults.newtonCotesCoefficients, 1)
      newtonCotesCoefficientsOnGrid = copy(JAC.Defaults.newtonCotesCoefficients) * grid.h

      i0     = 2
      gamma  = 0.
      result = 0.

      while abs(f(i0)) == 0 || f(i0)/f(i0+1) < 0
          i0 += 1
      end
      gamma = log(f(i0)/f(i0+1) * grid.rp[i0+1] / grid.rp[i0]) / log(grid.r[i0]/grid.r[i0+1])
      result = 1/(gamma + 1) * grid.r[i0] * f(i0) / grid.rp[i0]

      result += f(i0) * newtonCotesCoefficientsOnGrid[n]
      for i = i0+1:mtp
          result += newtonCotesCoefficientsOnGrid[(i - i0) % (n - 1) + 1] * f(i)
      end

      if (mtp - i0) % (n - 1) + 1 == 1
          result -= newtonCotesCoefficientsOnGrid[n] * f(mtp)
      end

      return result
  end

  """
  `Math.integrateTransform(f::Function, n0::Int64, mtp::Int64, grid::Radial.Grid)`  ... to integrate a (point-wise) given function over
       mtp grid points ...; a value::Float64 is returned.
  """
  function integrateTransform(f::Function, n0::Int64, mtp::Int64, grid::Radial.Grid)
      newtonCotesCoefficientsOnGrid             = Array{Float64}([2 * 7, 32, 12, 32, 7]) * 2 / 45  # * grid.h
      ## newtonCotesCoefficientsOnGrid            = Array{Float64}([2 * 41, 216, 27, 272, 27, 216, 41]) * 1 / 140
      ## newtonCotesCoefficientsOnGrid             = Array{Float64}([2 * 989, 5888, -928, 10496, -4540, 10496, -928, 5888, 989]) * 4 / 14175
      newtonCotesCoefficientsOnGrid *= grid.h
      n      = size(newtonCotesCoefficientsOnGrid, 1)

      i0     = n0
      gamma  = 0.
      result = 0.

      result += f(i0) * grid.rp[i0] * newtonCotesCoefficientsOnGrid[n]
      for i = i0+1:mtp
          result += newtonCotesCoefficientsOnGrid[(i - i0) % (n - 1) + 1] * f(i) * grid.rp[i]
      end

      if (mtp - i0) % (n - 1) + 1 == 1
          ##x result -= newtonCotesCoefficientsOnGrid[n] * f(mtp)
          result -= newtonCotesCoefficientsOnGrid[n] * f(mtp) * grid.rp[mtp]
      end

      return( result )
  end

  """
  `Math.polylogExp(x::Float64, s::Int64)`  ... to evaluate the polylog(s, -exp(x)) function.
  """
  function polylogExp(x::Float64, s::Int64)
      return -sf_fermi_dirac_int(s-1, x)
  end

end # module


text:
# Starting custom types example file for Learning Julia

# Create a new custom type using the struct keyword
struct MyType
    field1::Int
    field2::String
end

# instantiate the type
x = MyType(10, "ABC")
println(x)
println("field1 of x is ", x.field1, "\n")

# by default, fields in a type is immutable
# ex. x.field1 = 20 will result in an error

# use "mutable" to make a type that can be altered
mutable struct MyMutableType
    field1::Int
    field2::String
end

y = MyMutableType(10, "ABC")
y.field2 = "CDE"
println(y, "\n")

# use the isa() method to see if a variable is a given type
println("Is x a MyMutableType: ", isa(x, MyMutableType))
println("Is y a MyMutableType: ", isa(y, MyMutableType))

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
using Pkg
# change pwd
cd(joinpath(dirname(@__FILE__),".."))
# activate environment
Pkg.activate(pwd())
@info pwd()
using Pluto
Pluto.run()


text:
# === Design Notes ===
#
# `NullableArray{T, N}` is a struct-of-arrays representation of
# `Array{Nullable{T}, N}`. This makes it easy to define complicated operations
# (e.g. matrix multiplication) by reusing the existing definition for
# `Array{T}`.
#
# One complication when defining functions that operate on the internal fields
# of a `NullableArray` is that developers must take care to ensure that they
# do not index into an undefined entry in the `values` field. This is not a
# problem for `isbits` types, which are never `#undef`, but will trigger
# an exception for any other type.
#
# TODO: Ensure that size(values) == size(isnull) using inner constructor.
# TODO: Implement outer constructor required once we add an inner constructor.
@doc """
`NullableArray{T, N}` is an efficient alternative to `Array{Nullable{T}, N}`.
It allows users to easily define operations on arrays with null values by
reusing operations that only work on arrays without any null values.
""" ->
immutable NullableArray{T, N} <: AbstractArray{Nullable{T}, N}
    values::Array{T, N}
    isnull::Array{Bool, N}
    # extra field for potentially holding a reference to a parent memory block
    # (think mmapped file, for example) that `values` is actually derived from
    parent::Vector{UInt8}

    function NullableArray(d::AbstractArray{T, N}, m::Array{Bool, N}, parent::Vector{UInt8}=Vector{UInt8}())
        if size(d) != size(m)
            msg = "values and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        new(d, m, parent)
    end
end
typealias NullableVector{T} NullableArray{T, 1}
typealias NullableMatrix{T} NullableArray{T, 2}


text:
export Pool, MaxPool

"""
$(TYPEDEF)

Represents a pooling operation.

`p(x)` is shorthand for [`pool(x, p)`](@ref) when `p` is an instance of `Pool`.

## Fields:
$(FIELDS)
"""
struct Pool{N} <: Layer
    strides::NTuple{N,Int}
    pooling_function::Function
end

function Base.show(io::IO, p::Pool)
    (_, stride_height, stride_width, _) = p.strides
    function_display_name = Dict(MIPVerify.maximum => "max")
    print(
        io,
        "$(function_display_name[p.pooling_function]) pooling with a $(stride_height)x$(stride_width) filter and a stride of ($stride_height, $stride_width)",
    )
end

"""
$(SIGNATURES)

Convenience function to create a [`Pool`](@ref) struct for max-pooling.
"""
function MaxPool(strides::NTuple{N,Int}) where {N}
    Pool(strides, MIPVerify.maximum)
end

"""
$(SIGNATURES)

For pooling operations on an array where a given element in the output array
corresponds to equal-sized blocks in the input array, returns (for a given
dimension) the index range in the input array corresponding to a particular
index `output_index` in the output array.

Returns an empty array if the `output_index` does not correspond to any input
indices.

# Arguments
* `stride::Integer`: the size of the operating blocks along the active
     dimension.

"""
function getsliceindex(
    input_array_size::Integer,
    stride::Integer,
    output_index::Integer,
)::Array{Int,1}
    parent_start_index = (output_index - 1) * stride + 1
    parent_end_index = min((output_index) * stride, input_array_size)
    if parent_start_index > parent_end_index
        return []
    else
        return parent_start_index:parent_end_index
    end
end

"""
$(SIGNATURES)

For pooling operations on an array, returns a view of the parent array
corresponding to the `output_index` in the output array.
"""
function getpoolview(
    input_array::AbstractArray{T,N},
    strides::NTuple{N,Int},
    output_index::NTuple{N,Int},
)::SubArray{T,N} where {T,N}
    it = zip(size(input_array), strides, output_index)
    input_index_range = map(x -> getsliceindex(x...), it)
    return view(input_array, input_index_range...)
end

"""
$(SIGNATURES)

For pooling operations on an array, returns the expected size of the output
array.
"""
function getoutputsize(
    input_array::AbstractArray{T,N},
    strides::NTuple{N,Int},
)::NTuple{N,Int} where {T,N}
    output_size = ((x, y) -> round(Int, x / y, RoundUp)).(size(input_array), strides)
    return output_size
end

"""
$(SIGNATURES)

Returns output from applying `f` to subarrays of `input_array`, with the windows
determined by the `strides`.
"""
function poolmap(f::Function, input_array::AbstractArray{T,N}, strides::NTuple{N,Int}) where {T,N}
    output_size = getoutputsize(input_array, strides)
    output_indices = collect(CartesianIndices(output_size))
    return ((I) -> f(getpoolview(input_array, strides, I.I))).(output_indices)
end

"""
$(SIGNATURES)

Computes the result of applying the pooling function `params.pooling_function` to
non-overlapping cells of `input` with sizes specified in `params.strides`.
"""
function pool(input::AbstractArray{T,N}, params::Pool{N}) where {T<:JuMPReal,N}
    poolmap(params.pooling_function, input, params.strides)
end

(p::Pool)(x::Array{<:Real}) = MIPVerify.pool(x, p)
(p::Pool)(x::Array{<:JuMPLinearType}) =
    (Memento.info(MIPVerify.LOGGER, "Specifying $p ... "); MIPVerify.pool(x, p))


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

using Markdown
using InteractiveUtils

# ╔═╡ 7cbf25aa-3def-11eb-0c74-af66e15f45ca
s_ex = """mask = XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X
mem[8] = 11
mem[7] = 101
mem[8] = 0"""

# ╔═╡ 70c2c378-3e01-11eb-34c7-d17201765158
lines = readlines("14.txt")

# ╔═╡ 8fe5aaf2-3def-11eb-2708-d395cd38ffde
lines_ex = split(s_ex, "\n")

# ╔═╡ 0a9f638a-3dfc-11eb-3c1a-d3d4f9eab152
function maskedvalue(mask, value)
	bin_arr = lpad(string(value, base=2), 36, "0")
	masked = join([m == 'X' ? v : m for (m, v) in zip(mask, bin_arr)])
	parse(Int64, masked, base=2)
end

# ╔═╡ 40bb7a5e-3e05-11eb-046e-d9cef0009c51
function flipidx(binarr, idx, value="1")
	binarr[1:idx-1] * value * binarr[idx+1:end]
end

# ╔═╡ 7f768d70-3e04-11eb-0a0c-91bb5677385a
function maskedaddresses(mask, value)
	N = length(mask)
	bin_arr = lpad(string(value, base=2), N, "0")
	floating_bits = [i for i in 1:N if mask[i] == 'X']
	masked = join([m == '0' ? v : m for (m, v) in zip(mask, bin_arr)])
	addrs = [masked]
	for idx in floating_bits
		new_addrs = []
		for addr in addrs
			@assert addr[idx] == 'X'
			push!(new_addrs, flipidx(addr, idx, "0"))
			push!(new_addrs, flipidx(addr, idx, "1"))
		end
		addrs = new_addrs
	end
	# floating_bits, masked, addrs
	[parse(Int64, s, base=2) for s in addrs]
end

# ╔═╡ 9d343d2c-3def-11eb-155e-6ff7501e4fb6
function parseline(line, mask, memory::Dict)
	m = match(r"mask = ([X01]+)", line)
	if m !== nothing
		mask = m.captures[1]
		return mask, memory
	end
	m = match(r"mem\[([0-9]+)\] = ([0-9]+)", line)
	if m !== nothing
		mem_idx = parse(Int64, m.captures[1])
		value = parse(Int64, m.captures[2])
		memory[mem_idx] = maskedvalue(mask, value)
		return mask, memory
	end
	println("Unparsed `$(line)`")
end

# ╔═╡ fe31a3d8-3e00-11eb-3b37-cdc824855e51
begin
	mask = join("X" for i=1:36)
	memory = Dict()
	for line in lines
		mask, memory = parseline(line, mask, memory)
	end
	sum(values(memory))
end

# ╔═╡ e665342a-3e09-11eb-2e93-fb0813ae24af
function parseline2(line, mask, memory::Dict)
	m = match(r"mask = ([X01]+)", line)
	if m !== nothing
		mask = m.captures[1]
		return mask, memory
	end
	m = match(r"mem\[([0-9]+)\] = ([0-9]+)", line)
	if m !== nothing
		mem_idx = parse(Int64, m.captures[1])
		value = parse(Int64, m.captures[2])
		for addr in maskedaddresses(mask, mem_idx)
			memory[addr] = value
		end
		return mask, memory
	end
	println("Unparsed `$(line)`")
end

# ╔═╡ c0ee620a-3e09-11eb-2ded-7d1ffe2c0243
begin
	mask2 = join("X" for i=1:36)
	memory2 = Dict()
	for line in lines
		mask2, memory2 = parseline2(line, mask2, memory2)
	end
	sum(values(memory2))
end

# ╔═╡ Cell order:
# ╠═7cbf25aa-3def-11eb-0c74-af66e15f45ca
# ╠═70c2c378-3e01-11eb-34c7-d17201765158
# ╠═8fe5aaf2-3def-11eb-2708-d395cd38ffde
# ╠═0a9f638a-3dfc-11eb-3c1a-d3d4f9eab152
# ╠═40bb7a5e-3e05-11eb-046e-d9cef0009c51
# ╠═7f768d70-3e04-11eb-0a0c-91bb5677385a
# ╠═9d343d2c-3def-11eb-155e-6ff7501e4fb6
# ╠═fe31a3d8-3e00-11eb-3b37-cdc824855e51
# ╠═e665342a-3e09-11eb-2e93-fb0813ae24af
# ╠═c0ee620a-3e09-11eb-2ded-7d1ffe2c0243

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
# This file was generated, do not modify it. # hide
coefs = mtm.fitresult.fitresult.machine.fitresult
round.(coefs, sigdigits=2)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using COSMO
using Distributions
using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils
plotly()

##################### Read Prices #####################
prices = get_test_data();
numD, numA = size(prices) # A: Assets    D: Days

##################### Calculating returns #####################
returns_series = percentchange(prices);

# market_history = VolumeMarketHistory(returns_series)

##################### backtest Parameters #####################
DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 900000
)

date_range = timestamp(returns_series)[100:end];

## Strategies parameters
# Look back days
k_back = 60
# results
backtest_results = Dict()

##################### backtest stochastic strategies #####################

backtest_results["markowitz_limit_R"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = 0.0015 / market_budget(market)

    formulation = PortfolioFormulation(MIN_SENSE,
        ObjectiveTerm(SqrtVariance(d)),
        RiskConstraint(ExpectedReturn(d), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

backtest_results["markowitz_limit_var"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    V = 0.003 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(d)),
        RiskConstraint(SqrtVariance(d), LessThan(V)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

##################### backtest with RO strategies #####################

backtest_results["normal_limit_soyster"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = -0.05 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(d)),
        RiskConstraint(ExpectedReturn(BudgetSet(d, maximum(abs.(returns); dims=1)'[:] .- r̄, numA * 1.0), WorstCase), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

backtest_results["normal_limit_bertsimas_30"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = -0.05 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(d)),
        RiskConstraint(ExpectedReturn(BudgetSet(d, maximum(abs.(returns); dims=1)'[:] .- r̄, numA * 0.3), WorstCase), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

backtest_results["normal_limit_bertsimas_60"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = -0.05 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(d)),
        RiskConstraint(ExpectedReturn(BudgetSet(d, maximum(abs.(returns); dims=1)'[:] .- r̄, numA * 0.6), WorstCase), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

backtest_results["normal_limit_bental"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = -0.05 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(d)),
        RiskConstraint(ExpectedReturn(EllipticalSet(d, 0.1), WorstCase), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

##################### backtest with Data-Driven RO strategies #####################

backtest_results["mixed_signals_limit_betina"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)

    # Parameters
    j_robust = 25
    R = -0.015 / market_budget(market)
    num_train = 20
    kst = 14
    klt = 30
    kmom = 20
    Q = 2

    # Predict return
    r̄ = zeros(numA)
    for i in 1:numA
        r̄[i] = mixed_signals_predict_return(
            returns[:, i], num_train, kst, klt, kmom, Q
        )
    end

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedReturn(Dirac(r̄))),
        RiskConstraint(ConditionalExpectedReturn{WorstCase}(Inf, DeterministicSamples(returns'[:,:]), j_robust), GreaterThan(R)),
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

##################### backtest with Moment DRO strategies #####################

backtest_results["delage_utility"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    # Parameters
    d = MvNormal(r̄, Σ)
    R = -0.015 / market_budget(market)

    formulation = PortfolioFormulation(MAX_SENSE,
        ObjectiveTerm(ExpectedUtility(MomentUncertainty(d, 0.05, 1.3),
            PieceWiseUtility(
                [1.0], [0.0]
            ),
            WorstCase
        ))
    )

    pointers = change_bids!(market, formulation, DEFAULT_SOLVER)
    return pointers
end

##################### extra strategies #####################

backtest_results["max_sharpe"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)
    returns = values(past_returns)
    Σ, r̄ = mean_variance(returns[(end - k_back):end, :])

    change_bids!(market, market_budget(market) * max_sharpe(Σ, r̄, risk_free_rate(market)))

    return nothing
end

backtest_results["equal_weights"], _ = sequential_backtest_market(
    VolumeMarketHistory(returns_series), date_range,
) do market, past_returns, ext
    # Prep
    numD, numA = size(past_returns)

    change_bids!(market, market_budget(market) * equal_weights(numA))

    return nothing
end

##################### plot results  #####################
using Plots
using Plots.PlotMeasures

plt = plot(;title="Culmulative Wealth",
    xlabel="Time",
    ylabel="Wealth",
    legend=:outertopright,
    left_margin=10mm
);
for (strategy_name, recorders) in backtest_results
    plot!(plt,
        axes(get_records(recorders[:wealth]), 1), get_records(recorders[:wealth]).data;
        label=strategy_name,
    )
end
plt

plt2 = scatter(;
    title="Mean Vs Std",
    xlabel="σ",
    ylabel="E[r]",
    legend=:outertopright,
    left_margin=10mm
);
for (strategy_name, recorders) in backtest_results
    scatter!(plt2,
    [std(get_records(recorders[:returns]).data)], [mean(get_records(recorders[:returns]).data)];
        label=strategy_name,
    )
end
plt2

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
using AWS
using AWS: AWSCredentials
using AWS: AWSServices
using AWS.AWSExceptions: AWSException, NoCredentials
using AWS.AWSMetadataUtilities: _clean_documentation, _filter_latest_service_version,
    _generate_low_level_definition, _generate_high_level_definition, _generate_high_level_definitions,
    _get_aws_sdk_js_files, _get_service_and_version, _get_function_parameters, _clean_uri, _format_function_name,
    InvalidFileName, ProtocolNotDefined
using Base64
using Compat: mergewith
using Dates
using GitHub
using HTTP
using JSON
using OrderedCollections: LittleDict, OrderedDict
using MbedTLS: digest, MD_SHA256, MD_MD5
using Mocking
using Pkg
using Retry
using Suppressor
using Test
using UUIDs
using XMLDict

Mocking.activate()

include("patch.jl")

aws = AWSConfig()

function _now_formatted()
    return lowercase(Dates.format(now(Dates.UTC), dateformat"yyyymmdd\THHMMSSsss\Z"))
end

@testset "AWS.jl" begin
    include("AWS.jl")
    include("AWSCredentials.jl")
    include("AWSExceptions.jl")
    include("AWSMetadataUtilities.jl")
    include("issues.jl")
    include("test_pkg.jl")
    include("utilities.jl")
end


text:
mutable struct kvh_raw_imu_batch_t <: LCMType
    utime::Int64
    num_packets::Int32
    raw_imu::Vector{kvh_raw_imu_t}
end

@lcmtypesetup(kvh_raw_imu_batch_t,
    raw_imu => (num_packets,)
)


text:
"""
    Module EntropyStableEuler

Includes general math tools
"""

module EntropyStableEuler
using StaticArrays

const γ = 1.4

export logmean, γ
include("./logmean.jl")

# # submodules
# export Fluxes1D,Fluxes2D,Fluxes3D

# export entropy_scaling,scale_entropy_input,scale_entropy_output
# changes definition of entropy variables by a constant scaling
#const entropy_scaling = 1/γ # implies vE = -1 /(ι * γ) and γι = non-dim temperature if R* = cp
const entropy_scaling = 1
scale_entropy_output(V...) = (x->@. x*entropy_scaling).(V)
scale_entropy_input(V...) = (x->@. x/entropy_scaling).(V)

#####
##### one-dimensional fluxes
#####

module Fluxes1D
import ..γ
import ..entropy_scaling
import ..EntropyStableEuler: logmean
import ..EntropyStableEuler: scale_entropy_output, scale_entropy_input

export wavespeed_1D

export primitive_to_conservative, conservative_to_primitive_beta
export u_vfun, v_ufun
export euler_fluxes
export Sfun,pfun,betafun
include("./euler_fluxes_1D.jl")
include("./entropy_variables.jl")
end

#####
##### two-dimensional fluxes
#####
module Fluxes2D
import ..γ
import ..entropy_scaling
import ..EntropyStableEuler: logmean, entropy_scaling
import ..EntropyStableEuler: scale_entropy_output, scale_entropy_input

export primitive_to_conservative,conservative_to_primitive_beta
export u_vfun, v_ufun
export Sfun,pfun,betafun
export euler_fluxes

include("./entropy_variables.jl")
include("./euler_fluxes_2D.jl")
end

#####
##### three-dimensional fluxes
#####
module Fluxes3D
import ..γ
import ..entropy_scaling
import ..EntropyStableEuler: logmean, entropy_scaling
import ..EntropyStableEuler: scale_entropy_output, scale_entropy_input

export primitive_to_conservative,conservative_to_primitive_beta
export u_vfun, v_ufun
export Sfun,pfun,betafun
export euler_fluxes

include("./entropy_variables.jl")
include("./euler_fluxes_3D.jl")
end

# export u_vfun, v_ufun, betafun, pfun, rhoe_ufun
# export dVdU_explicit, dUdV_explicit
# export wavespeed # c
# export Sfun, sfun # math/physical entropies
# export u_vfun1D, v_ufun1D, betafun1D # specialization to 1D
# export primitive_to_conservative
# include("./euler_variables.jl")

# export euler_fluxes
# export euler_flux_x, euler_flux_y # separate x,y fluxes for faster implicit assembly using ForwardDiff
# include("./euler_fluxes.jl")

# export vortex
# include("./analytic_solutions.jl")

end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
export step!, run!

"""
update e.elites with the best individuals from the current population or
existing elites, based on fitness
"""
function elites_generation!(e::AbstractEvolution)
    pop = sort([e.population; e.elites])
    n_elites = length(e.elites)
    for i in 1:n_elites
        e.elites[i] = deepcopy(pop[length(e.population)+i])
    end
end

"""
generation should be used for any additional record keeping necessary in an
algorithm. It is called after populate and evaluate, so all individuals should
have genes and fitness reflecting their current state. The default function
calls elites_generation if e.elites is defined, otherwise does nothing.
"""
function generation!(e::AbstractEvolution)
    if hasfield(typeof(e), :elites)
        elites_generation!(e)
    end
end

"""
The generic iteration of an evolution. Calls populate!, evaluate!, and generation!.
Also calls log_gen and save_gen based on the provided config values. Subclasses
of AbstractEvolution should override the populate!, evaluate!, or generation!
functions rather than overriding this function.
"""
function step!(e::AbstractEvolution)
    e.gen += 1
    if e.gen > 1
        populate!(e)
    end
    evaluate!(e)
    generation!(e)
    if ((e.config.log_gen > 0) && mod(e.gen, e.config.log_gen) == 0)
        log_gen(e)
    end
    if ((e.config.save_gen > 0) && mod(e.gen, e.config.save_gen) == 0)
        save_gen(e)
    end
end

"Call step!(e) e.config.n_gen times consecutively"
function run!(e::AbstractEvolution)
    for i in (e.gen+1):e.config.n_gen
        step!(e)
    end
end


text:
## A General Structure for D.R.I.P.
"""
# Summary
    A type srtucture for storing the steady state solution of a Drip.

# Fields
    K      : Kalman gain matrix
    Y      : Weight vector for evolution of actions
    Σ_z    : Covariance matrix of the rational inattention error
    Σ_p    : Steady-state posterior covariance matrix under the solution
    Σ_1    : Steady-state prior covariance matrix under the solution
    Ω      : Dynamic benefit matrix
    err    : Convergence error for the solution
"""
struct SteadyState
    K   :: Array{Float64,2}
    Y   :: Array{Float64,2}
    Σ_z :: Array{Float64,2}
    Σ_p :: Array{Float64,2}
    Σ_1 :: Array{Float64,2}
    Ω   :: Array{Float64,2}
    err :: Float64
end

"""
# Summary
  A Type Structure for LQG Dynamic Rational Inattention Problems (DRIPs)

# Fields
## Primitives of the DRIP
    ω      : Cost of information
    β      : Discount factor
    A      : Transition matrix: x=Ax+Qu
    Q      : Std. Dev. matrix: x=Ax+Qu
    H      : Mapping of shocks to actions: v=-0.5(a'-x'H)(a-H'x)

## Solution of the DRIP in the Steady State
    ss     : Steady State Solution as a SteadyState type
    See also: [`SteadyState`](@ref)
"""
struct Drip
# primitives
    ω   :: Float64
    β   :: Float64
    A   :: Array{Float64,2}
    Q   :: Array{Float64,2}
    H   :: Array{Float64,2}
# steady state
    ss  :: SteadyState
end

"""
    Drip(ω,β,A,Q,H; kwargs...) -> Drip

Solves for the steady state of a Dynamic Rational Inattention Problem (DRIP)
    defined by the arguments and stores the solution in a Drip type.
    See [Afrouzi and  Yang (2019)](http://afrouzi.com/dynamic_inattention.pdf)
    for details.
# Arguments
The function takes the primitives of the Drip as arguments:

    * ω      : Cost of information
    * β      : Discount factor
    * A      : Transition matrix: x=Ax+Qu
    * Q      : Std. Dev. matrix: x=Ax+Qu
    * H      : Mapping of shocks to actions: v=-0.5(a'-x'H)(a-H'x)
## Optional Arguments (kwargs...)
Default values are set unless specified otherwise by user.

    * fcap  [= false]    [if `true` then solves the problem with fixed capacity = ω bits]
    * Ω0    [= H*H']     [initial guess for steady state information matrix]
    * Σ0    [= A*A'+Q*Q'][initial guess for steady state prior]
    * w     [= 1]        [updating weight on the new guess in iteration]
    * tol   [= 1e-4]     [tolerance level for convergence]
    * maxit [= 10000]    [maximum number of iterations]
# Outputs
The function returns a `Drip` structure with the primitives and the solution objects:

    * Y      : Weight vector for evolution of actions
    * Σ_z    : Covariance matrix of the rational inattention error
    * K      : Kalman gain matrix
    * Σ_1    : Steady-state prior covariance matrix under the solution
    * Σ_p    : Steady-state posterior covariance matrix under the solution
    * Ω      : Dynamic benefit matrix
# Examples
```julia-repl
julia> P = Drip(ω,β,A,Q,H)
```
"""
function Drip(ω,β,A,Q,H;         # primitives of the D.R.I.P.
              fcap  = false,     # optional: if true then solves the problem with fixed capacity κ = ω.
              Ω0    = H*H',      # optional: initial guess for steady state information matrix
              Σ0    = A*A'+Q*Q', # optional: initial guess for steady state prior
              w     = 1,         # optional: updating weight in iteration
              tol   = 1e-4,      # optional: tolerance level for convergence
              maxit = 10000      # optional: maximum number of iterations
              )

    A, Q, H = collect(A)[:,:], collect(Q)[:, :], collect(H)[:, :];

    ## initialize
    (n,m) = size(H)
    # n: dimension of state, m: number of actions
    eye   = Matrix{Float64}(I,n,n);
    err   = 1
    iter  = 0
    SqRΣ  = sqrt(Σ0);
    Σ1    = Matrix{Float64}(I,n,n)
    Ω1    = Matrix{Float64}(I,n,n)
    Σp    = Matrix{Float64}(I,n,n)
    κ     = ω
    Ω_c   = H*H';
    Σq    = Q*Q';
    # iterate
    while (err > tol) & (iter < maxit)
        D, U    = eigen!(Symmetric(SqRΣ*Ω0*SqRΣ));
        D       = diagm(D);
        U       = real.(U);

        if fcap == true
            ω = (2^(2*κ)/det(max.(ω*eye,D)))^(-1/n);
        end
        Σp      = Symmetric(ω*SqRΣ*U/(max.(D,ω*eye))*U'*SqRΣ);

        SqRΣ    = real.(sqrt(Σ0));
        invSqRΣ = real.(inv(SqRΣ));

        Σ1      = Symmetric(A*Σp*A' + Σq);
        Ω1      = Ω_c .+ β*A'*invSqRΣ*U*(min.(D,ω*eye))*U'*invSqRΣ*A
        err     = norm(Σ1 - Σ0,1)/norm(Σ0,1) + norm(Ω1 - Ω0,1)/norm(Ω0,1);

        Σ0      = w*Σ1 + (1-w)*Σ0;
        Ω0      = w*Ω1 + (1-w)*Ω0;
        Ω0      = (abs.(Ω0).>1e-10).*Ω0;

        iter   += 1
    end

    if iter == maxit
        print("RI Code hit maxit\n")
    end
    inv_Σ1 = inv(Σ1) ;
    Σ_1  = collect(Σ1)[:,:];
    Σ_p  = collect(Σp)[:,:];
    Y    = collect((eye - Σ_p*inv_Σ1)'*H)[:,:];
    Σ_z  = collect(H'*(Σ_p - Σ_p*inv_Σ1*Σ_p)*H)[:,:];
    K    = collect(Σ1*Y*inv(Y'*Σ1*Y .+ Σ_z))[:,:];
    Ω    = collect(Ω1)[:,:];
    err  = err;
    return(Drip(ω,β,A,Q,H,SteadyState(K,Y,Σ_z,Σ_p,Σ_1,Ω,err)))
end


text:
using HypothesisTests, Test, Statistics, StableRNGs

@testset "Permutation" begin
    @testset "ExactPermutationTest" begin
        @test pvalue(ExactPermutationTest([1,2,3], [4,5,6], mean), tail=:both) ≈ 0.1
        @test pvalue(ExactPermutationTest([4,5,6], [1,2,3], mean), tail=:both) ≈ 0.1
        @test pvalue(ExactPermutationTest([1,2,3], [4,5,6], mean), tail=:left) ≈ 0.05
        @test pvalue(ExactPermutationTest([4,5,6], [1,2,3], mean), tail=:right) ≈ 0.05
    end

    @testset "ApproximatePermutationTest" begin
        rng = StableRNG(12345)
        @test pvalue(ApproximatePermutationTest(rng, [1,2,3], [4,5,6], mean, 100), tail=:both) ≈ 0.14
        @test pvalue(ApproximatePermutationTest(rng, [4,5,6], [1,2,3], mean, 100), tail=:both) ≈ 0.08
        @test pvalue(ApproximatePermutationTest(rng, [1,2,3], [4,5,6], mean, 100), tail=:left) ≈ 0.05
        @test pvalue(ApproximatePermutationTest(rng, [4,5,6], [1,2,3], mean, 100), tail=:right) ≈ 0.05
        # Test that the non-rng method works
        @test ApproximatePermutationTest([4,5,6], [1,2,3], mean, 100) isa HypothesisTests.PermutationTest{Float64}
    end
end


text:
module BeurlingLasso

using DynamicPolynomials
using SemialgebraicSets
using SumOfSquares
using MosekTools
using Optim
using GLMNet
using LineSearches

include("blasso.jl")
include("momsos.jl")
include("sfw.jl")
include("utils.jl")

end


text:
using PrettyMacroExpand
using Test

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


text:
using InteractiveUtils

struct Parameter
  p::Float64
  b::Float64
  r::Float64
end

function lorenz(a, p)
  return [
    p.p * (a[2] - a[1]),
    a[1]*(p.r - a[3]) - a[2],
    a[1]*a[2] - p.b * a[3]
  ]
end

function main()
  x0 = [1.0, 0, 0]
  p0 = Parameter(10, 28, 8/3)
  dt = 0.001
  for t in 0:1000
    x0 += dt * lorenz(x0, p0)
    println(x0)
  end
end

code_llvm(main, ())

max_stars_repo_name: termoshtt/julia-kernel-experiment

max_issues_repo_name: termoshtt/julia-kernel-experiment

max_forks_repo_name: termoshtt/julia-kernel-experiment


text:
# * Envelopes
#   The envelopes implemented below are all /vector potential/
#   envelopes, since that is what is being used in
#   calculations. However, they may be specified using
#   intensity-related quantities, e.g. Gaussian pulses are most often
#   specified using the FWHM duration of their /intensity/ envelopes.
#

envelope_types = Dict{Symbol,Any}()

complex_squared(t) = t*t

# ** Gaussian
@doc raw"""
    GaussianEnvelope

A Gaussian pulse is given by

```math
f(t) = \exp\left(-\frac{t^2}{2\sigma^2}\right),
```

where the standard deviation ``σ`` is related to the FWHM duration τ
of the intensity envelope as

```math
\sigma = \frac{\tau}{2\sqrt{2\ln 2}}.
```

Since we define all fields in terms of the vector potential ``A(t)``,
and the [`instantaneous_intensity`](@ref) is given by
``\abs{-\partial_t A(t)}^2``, we have to find a coefficient ``α`` such
that the field amplitude

```math
F(t) \sim -\partial_t \exp(-\alpha t^2) \sin(\omega t + \phi)
```

has an _intensity envelope_ with the desired FWHM; we do this
iteratively in [`gaussian_common!`](@ref). This is mainly important
for ultrashort pulses, since for long pulse durations we can
approximate ``\exp(-\alpha t^2) \sim 1`` such that only the carrier
contributes to the time derivative.

Since a Gaussian never ends, we specify how many ``σ`` we require; the
resulting time window will be rounded up to an integer amount of
cycles of the fundamental.

Required parameters:
- `λ|T|f|ν|ω|ħω`,
- `τ|σ`,
- `σmax|tmax|Tmax`,
- `env=:gauss` (optional, since [`GaussianEnvelope`](@ref) is the default envelope).
"""
struct GaussianEnvelope{T} <: AbstractEnvelope
    τ::T # Intensity FWHM
    σ::T # Intensity std.dev.
    α::T # Vector potential exponential coefficient
    σmax::T
    tmax::T # Maximum time. Time window: [-tmax,tmax]
end
envelope_types[:gauss] = GaussianEnvelope

(env::GaussianEnvelope)(t) = exp(-env.α*complex_squared(t))

show(io::IO, env::GaussianEnvelope) =
    printfmt(io, "Gaussian envelope of duration {1:s} (intensity FWHM; ±{2:0.2f}σ)",
             au2si_round(env.τ, u"s"), env.σmax)

@doc raw"""
    gaussian_common!(field_params, carrier[; Tmax_rounder, verbosity])

Compute parameters common to Gaussian envelopes,
i.e. [`GaussianEnvelope`](@ref) and
[`TruncatedGaussianEnvelope`](@ref); most importantly, given an
intensity FWHM or σ, we need to figure out the coefficient ``\alpha``
for the envelope of the vector potential ``\exp(-\alpha t^2)``, such
that the electric field amplitude and intensity have the desired
durations. The optional function `Tmax_rounder` determines if the
envelop should be extended to encompass an integer amount of cycles of
the `carrier` (default).
"""
function gaussian_common!(field_params, carrier;
                          Tmax_rounder = Base.Fix1(ceil, Int), verbosity=0)
    @namespace!(field_params) do
        if τ
            σ = τ/(2*√(2log(2)))
        else
            τ = 2*√(2log(2))*σ
        end

        if σmax
            Tmax = Tmax_rounder(σmax*σ/T)
            tmax = Tmax*T
        else
            if tmax
                Tmax = Tmax_rounder(tmax/T)
            elseif Tmax
                tmax = Tmax*T
            end
            tmax = Tmax*T
        end
        σmax = tmax/σ
    end
    τ = austrip(field_params[:τ])
    σ = austrip(field_params[:σ])
    α = 2log(2)/τ^2
    period = austrip(field_params[:T])
    ω = austrip(field_params[:ω])

    σmax = austrip(field_params[:σmax])
    tmax = austrip(field_params[:tmax])

    τ < 0.5period &&
        @warn "Pulse durations smaller than half a period not reliably supported"

    # Find the exponential coefficient for the vector potential that
    # will give the desired FWHM for the intensity profile.
    f = α -> begin
        env = GaussianEnvelope(τ, σ, α[1], σmax, tmax)
        field = make_temp_field(carrier, env, field_params)
        abs2(intensity(field, τ/2) - ω^2/2)
    end
    verbosity > 0 && @info "Finding optimal vector potential coefficient"
    res = optimize(f, [α], BFGS())
    verbosity > 0 && display(res)

    field_params[:α] = res.minimizer[1]
end

function GaussianEnvelope(field_params::Dict{Symbol,Any}, carrier)
    test_field_parameters(field_params, [:T]) # Period time required to round time window up
    test_field_parameters(field_params, [:τ, :σ])
    test_field_parameters(field_params, [:σmax, :tmax, :Tmax])

    gaussian_common!(field_params, carrier)

    @unpack τ, σ, α, σmax, tmax = field_params
    GaussianEnvelope(austrip(τ), austrip(σ), α, austrip(σmax), austrip(tmax))
end

continuity(::GaussianEnvelope) = Inf
span(env::GaussianEnvelope) = -env.tmax..env.tmax

time_integral(env::GaussianEnvelope) = env.σ*√(2π)

# *** Spectrum
@doc raw"""
    spectrum(env::GaussianEnvelope)

Gaussians belong to the [Schwartz
class](https://en.wikipedia.org/wiki/Schwartz_space), i.e. functions
who, under Fourier transform, are mapped back to the same space. That
is to say, the Fourier transform of a Gaussian is a Gaussian:

```math
\exp(-\alpha t^2) \leftrightarrow
\frac{1}{\sqrt{2\alpha}}
\exp\left(-\frac{\omega^2}{4\alpha}\right).
```

Comparing with the above, we find that the spectral standard
deviation

```math
\Omega = \sqrt{2\alpha} = \frac{2\sqrt{\ln 2}}{\tau},
```

and the Gaussian function in the spectral domain is thus

```math
E(\omega) =
\frac{E_0\tau}{2\sqrt{\ln 2}}
\exp\left[-\frac{(\omega\tau)^2}{8\ln2}\right].
```
"""
function spectrum(env::GaussianEnvelope)
    N = env.E₀*env.τ/(2*√(log(2)))
    ω -> N*exp(-(ω*env.τ)^2/8log(2))
end

# ** Truncated Gaussian

@doc raw"""
    TruncatedGaussianEnvelope

Since a Gaussian function never ends, this envelope adds a soft truncation over a time interval:

```math
f(t)=
\begin{cases}
\exp(-\alpha t^2), & \abs{t} \leq t_{\textrm{off}},\\
\exp\left\{
-\alpha\left[
t_{\textrm{off}} + \frac{2}{\pi}(t_{\textrm{max}} - t_{\textrm{off}})
\tan\left(
\frac{\pi}{2}
\frac{\abs{t} - t_{\textrm{off}}}{t_{\textrm{max}} - t_{\textrm{off}}}
\right)
\right]^2
\right\},
& t_{\textrm{off}} < \abs{t} \leq t_{\textrm{max}}, \\
0, & t_{\textrm{max}} < \abs{t}.
\end{cases}
```
This is Eq. (72) of

- Patchkovskii, S., & Muller, H. (2016). Simple, accurate, and
  efficient implementation of 1-electron atomic time-dependent
  Schrödinger equation in spherical coordinates. Computer Physics
  Communications, 199,
  153–169. [10.1016/j.cpc.2015.10.014](http://dx.doi.org/10.1016/j.cpc.2015.10.014)

Required parameters:
- `λ|T|f|ν|ω|ħω`,
- `τ|σ`,
- `toff`,
- `σmax|tmax|Tmax`,
`env=:trunc_gauss`.

Beyond this, everything else is the same as for
[`GaussianEnvelope`](@ref).
"""
struct TruncatedGaussianEnvelope{T} <: AbstractEnvelope
    τ::T # Intensity FWHM
    σ::T # Intensity std.dev.
    α::T # Vector potential exponential coefficent
    toff::T # Start of hard turn-off
    tmax::T # Maximum time; end of hard turn-off. Time window: [-tmax,tmax]
    Tmax::T # Maximum time, in cycles of the fundamental.
end
envelope_types[:trunc_gauss] = TruncatedGaussianEnvelope

function (env::TruncatedGaussianEnvelope{T})(t) where T
    at = abs(real(t))
    at > env.tmax && return zero(T)
    @unpack α, toff, tmax = env
    if at ≤ env.toff
        exp(-α*complex_squared(t))
    else
        exp(-α*complex_squared(toff + 2/π*(tmax-toff)*tan(π/2*(at-toff)/(tmax-toff))))
    end
end

show(io::IO, env::TruncatedGaussianEnvelope) =
    printfmt(io, "Truncated Gaussian envelope of duration {1:.4f} jiffies = {2:s} (intensity FWHM; turn-off from {3:s} to {4:s})",
             env.τ, au2si_round(env.τ, u"s"),
             au2si_round(env.toff, u"s"), au2si_round(env.tmax, u"s"))

function TruncatedGaussianEnvelope(field_params::Dict{Symbol,Any}, carrier)
    test_field_parameters(field_params, [:T]) # Period time required to round time window up
    test_field_parameters(field_params, [:τ, :σ])
    test_field_parameters(field_params, [:toff, :σoff])
    test_field_parameters(field_params, [:σmax, :tmax, :Tmax])

    gaussian_common!(field_params, carrier, Tmax_rounder=NoUnits)
    @namespace!(field_params) do
        if σoff
            toff = σ*σoff
        end
    end

    @unpack τ, σ, α, toff, tmax, Tmax = field_params
    toff < tmax || throw(ArgumentError("Hard turn-off must occur before end of pulse"))
    TruncatedGaussianEnvelope(austrip(τ), austrip(σ), α, austrip(toff), austrip(tmax), Tmax)
end

continuity(::TruncatedGaussianEnvelope) = Inf # This is not exactly true
span(env::TruncatedGaussianEnvelope) = -env.tmax..env.tmax

# TODO: Take truncation into account
time_integral(env::TruncatedGaussianEnvelope) = env.σ*√(2π)

# ** Trapezoidal

@doc raw"""
    TrapezoidalEnvelope

This is a very simple piecewise linear function:

```math
f(t)=
\begin{cases}
r/r_{\textrm{up}},
& 0 \leq r < r_{\textrm{up}},\\
1,
& r_{\textrm{up}} \leq r < r_{\textrm{up}} + r_{\textrm{flat}}, \\
1 - \frac{r-r_{\textrm{up}}-r_{\textrm{flat}}}{r_{\textrm{down}}},
& r_{\textrm{up}} + r_{\textrm{flat}} \leq r < r_{\textrm{up}} + r_{\textrm{flat}} + r_{\textrm{down}}, \\
0,
& \textrm{else},
\end{cases}
\quad
r \defd t/T,
```
where ``T`` is the period time of the carrier.

Required parameters:
- `λ|T|f|ν|ω|ħω`,
- `ramp | (ramp_up & ramp_down)`,
- `flat`,
- `env=:trapezoidal | env=:tophat`.

Beware that this envelope can introduce artifacts at the ends of the
pulse, such that the electric field is non-vanishing, depending of
e.g. the phase of the carrier.
"""
struct TrapezoidalEnvelope{T} <: AbstractEnvelope
    ramp_up::T
    flat::T
    ramp_down::T
    period::T

    function TrapezoidalEnvelope(ramp_up::T, flat::T, ramp_down::T, period::T) where T
        ramp_up >= 0 || error("Negative up-ramp not supported")
        flat >= 0 || error("Negative flat region not supported")
        ramp_down >= 0 || error("Negative down-ramp not supported")
        ramp_up + flat + ramp_down > 0 || error("Pulse length must be non-zero")

        new{T}(ramp_up, flat, ramp_down, period)
    end
end
envelope_types[:trapezoidal] = TrapezoidalEnvelope
# Common alias
envelope_types[:tophat] = TrapezoidalEnvelope

function (env::TrapezoidalEnvelope{T})(t) where T
    t = real(t)/env.period
    if t < 0
        zero(T)
    elseif t < env.ramp_up
        t/env.ramp_up
    elseif t < env.ramp_up + env.flat
        one(T)
    elseif t < env.ramp_up + env.flat + env.ramp_down
        one(T) - (t-env.ramp_up-env.flat)/env.ramp_down
    else
        zero(T)
    end
end

show(io::IO, env::TrapezoidalEnvelope) =
    printfmt(io, "/{1:d}‾{2:d}‾{3:d}\\ cycles trapezoidal envelope",
             env.ramp_up, env.flat, env.ramp_down)

function TrapezoidalEnvelope(field_params::Dict{Symbol,Any}, args...)
    test_field_parameters(field_params, [:T]) # Period time required to relate ramps/flat to cycles
    test_field_parameters(field_params, [:ramp, :ramp_up])
    test_field_parameters(field_params, [:ramp, :ramp_down])
    test_field_parameters(field_params, [:flat])

    @namespace!(field_params) do
        if ramp
            ramp_up = ramp
            ramp_down = ramp
        end
    end

    @unpack ramp_up, flat, ramp_down, T = field_params
    TrapezoidalEnvelope(ramp_up, flat, ramp_down, austrip(T))
end

continuity(::TrapezoidalEnvelope) = 0
span(env::TrapezoidalEnvelope{T}) where T =
    zero(T)..((env.ramp_up + env.flat + env.ramp_down)*env.period)

# ** Cos²

@doc raw"""
    Cos²Envelope

```math
f(t) = \begin{cases}
\cos^2 \left(\frac{\pi t}{cT}\right),
& -1 \leq t/(cT) \leq 1,\\
0, & \textrm{else},
\end{cases}
```
where ``c`` is the number of cycles from zero to zero of the
``\cos^2`` envelope and ``T`` the period time.

Required parameters:
- `λ|T|f|ν|ω|ħω`,
- `cycles`,
- `env=:cos² | env=:cos2`.
"""
struct Cos²Envelope{T} <: AbstractEnvelope
    cycles::T
    period::T

    function Cos²Envelope(cycles::T, period::T) where T
        cycles >= 0 || error("Negative duration not supported")
        new{T}(cycles, period)
    end
end
envelope_types[:cos²] = Cos²Envelope
envelope_types[:cos2] = Cos²Envelope

function (env::Cos²Envelope{T})(t) where T
    t /= (env.cycles*env.period)
    if -1 ≤ 2real(t) ≤ 1
        complex_squared(cospi(t))
    else
        zero(T)
    end
end

show(io::IO, env::Cos²Envelope) =
    printfmt(io, "{1:0.2f} cycles cos² envelope", env.cycles)

function Cos²Envelope(field_params::Dict{Symbol,Any}, args...)
    test_field_parameters(field_params, [:T]) # Period time required to relate ramps/flat to cycles
    test_field_parameters(field_params, [:cycles])

    @unpack cycles, T = field_params
    Cos²Envelope(cycles, austrip(T))
end

continuity(::Cos²Envelope) = 0
function span(env::Cos²Envelope)
    s = env.cycles*env.period/2
    -s..s
end

# ** Exports

export continuity


text:
@testset "Internal Methods" begin
    @testset "_support_sum_coeffs" begin
        @test IOMT._support_sum_coeffs([1, 2]) == [1, 1]
    end
end

# TODO test support_sum with JuMP container of infinite parameters

@testset "Support Sum" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in [0, 1])
    @infinite_parameter(m, y[1:2] in [0, 1])
    @infinite_variable(m, inf(x, y))

    @test InfiniteOpt._index_type(support_sum(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(support_sum(inf, y)) == MeasureIndex
    @test name(support_sum(inf, x)) == "support_sum"
end

@testset "Macro" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in [0, 1])
    @infinite_parameter(m, y[1:2] in [0, 1])
    @infinite_variable(m, inf(x, y))

    @test InfiniteOpt._index_type(@support_sum(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(@support_sum(inf, y)) == MeasureIndex
end


text:
module Modifyfield

@generated function copy_and_modify{fieldname}(x, ::Type{Val{fieldname}}, newval)
    F = fieldnames(x)
    mask = F .== fieldname
    if sum([convert(Int,mask[j]) for j = 1 : length(mask)]) != 1
        error("Type $x does not have a field $fieldname")
    end
    D = Any[:(x.$f) for f in F]
    D[mask] = :newval
    Expr(:call, x, D...)
end

@generated function copy_and_modify_mult{tupfieldnames}(x,
                                                        ::Type{Val{tupfieldnames}},
                                                        newvals...)
    F = fieldnames(x)
    D = Any[:(x.$f) for f in F]
    for (pos,fn) in enumerate(tupfieldnames)
        mask = F .== fn
        if sum([convert(Int,mask[j]) for j = 1 : length(mask)]) != 1
            error("Type $x does not have a field $fn")
        end
        D[mask] = :(newvals[$pos])
    end
    Expr(:call, x, D...)
end

macro modify_fields!(ex)
    ex.head == :(.) && ex.args[2].head == :tuple || error("Invalid usage of @modify_fields! macro")
    lhs = ex.args[1]
    arrfieldnames = Any[]
    newvals = Any[]
    for assn in ex.args[2].args
        assn.head == :(=) || error("Invalid usage of @modify_fields! macro")
        push!(arrfieldnames, Symbol(assn.args[1]))
        push!(newvals, assn.args[2])
    end
    # generate an expression like:
    #  lhs = copy_and_modifymult(lhs, Val{(fieldname1,fieldname2,...), newval1,
    #                            newval2, ...)
    # Do it the long way because I couldn't figure out where to put the
    # esc calls and ...'s.
    res = Expr(:(=))
    push!(res.args, lhs)
    push!(res.args, Expr(:call))
    callargs = res.args[2].args
    push!(callargs, :copy_and_modify_mult)
    push!(callargs, lhs)
    push!(callargs, Expr(:curly))
    curlyargs = callargs[3].args
    push!(curlyargs, :Val)
    push!(curlyargs, Expr(:tuple))
    for fn in arrfieldnames
        push!(curlyargs[2].args, QuoteNode(fn))
    end
    for nv in newvals
        push!(callargs, nv)
    end
    esc(res)
end

macro modify_field!(ex)
    ex.head == :(=) && ex.args[1].head == :(.) || error("Invalid usage of @modify_field! macro")
    lhs, newval = ex.args
    x, g = lhs.args
    esc(:($x = copy_and_modify($x, Val{$g}, $newval)))
    # :($(esc(x)) = copy_and_modify($(esc(x)), Val{$g}, $(esc(newval))))
end

@generated function copy_and_modify_tup{N,i}(x::NTuple{N}, ::Type{Val{i}}, newval)
    (i < 1 || i > N) && error("Tuple subscript $i is out of range for length-$N tuple")
    D = Any[:(x[$j]) for j = 1 : N]
    D[i] = :newval
    Expr(:tuple, D...)
end

macro modify_tuple_entry!(ex)
    ex.head == :(=) && ex.args[1].head == :ref || error("Invalid usage of @modify_tuple_entry! macro")
    lhs, newval = ex.args
    x, sub = lhs.args
    #:($(esc(x)) = copy_and_modify_tup($(esc(x)), Val{$sub}, $(esc(newval))))
    esc(:($x = copy_and_modify_tup($x, Val{$sub}, $newval)))
end

export copy_and_modify,
  copy_and_modify_mult,
  copy_and_modify_tup,
  @modify_field!,
  @modify_tuple_entry!

end


text:
using Test
using BetweenFlags.FeaturedFuncs

@testset "Get flat" begin
  s_i1 = "Some text... {GRAB THIS}, some more text {GRAB THIS TOO}..."
  L_o1 = get_flat(s_i1, ["{"], ["}"])
  s_i2 = "Some text... {GRAB THIS}, some more text {GRAB THIS TOO}..."
  L_o2 = get_flat(s_i2, ["{"], ["}"], false)
  s_i3 = "Some text... {GRAB THIS), } some more text {GRAB THIS TOO}..."
  L_o3 = get_flat(s_i3, ["{"], ["}", ")"])
  s_i4 = "Some text... {GRAB THIS), } some more text {GRAB THIS TOO}..."
  L_o4 = get_flat(s_i4, ["{"], ["}", ")"], false)
  s_i5 = "Some text... GRAB NOTHING), some more text GRAB NOTHING..."
  L_o5 = get_flat(s_i5, ["{"], ["}", ")"], false)

  @test L_o1[1]=="{GRAB THIS}"
  @test L_o2[1]=="GRAB THIS"
  @test L_o3[1]=="{GRAB THIS)"
  @test L_o4[1]=="GRAB THIS"
  @test L_o1[2]=="{GRAB THIS TOO}"
  @test L_o2[2]=="GRAB THIS TOO"
  @test L_o3[2]=="{GRAB THIS TOO}"
  @test L_o4[2]=="GRAB THIS TOO"
  @test L_o5[1]==""
end

@testset "Get level flat" begin
  s_i1 = "Some text... {GRAB {THIS}}, some more text {GRAB THIS TOO}..."
  L_o1 = get_level_flat(s_i1, ["{"], ["}"])
  s_i2 = "Some text... {GRAB {THIS}}, some more text {GRAB THIS TOO}..."
  L_o2 = get_level_flat(s_i2, ["{"], ["}"], false)
  s_i3 = "Some text... {GRAB {THIS}), } some more text {GRAB THIS TOO}..."
  L_o3 = get_level_flat(s_i3, ["{"], ["}", ")"])
  s_i4 = "Some text... {GRAB {THIS}), } some more text {GRAB THIS TOO}..."
  L_o4 = get_level_flat(s_i4, ["{"], ["}", ")"], false)
  s_i5 = "Some text... GRAB {NOTHING),  some more text GRAB NOTHING..."
  L_o5 = get_level_flat(s_i5, ["{"], ["}"], false)

  @test L_o1[1]=="{GRAB {THIS}}"
  @test L_o2[1]=="GRAB {THIS}"
  @test L_o3[1]=="{GRAB {THIS})"
  @test L_o4[1]=="GRAB {THIS}"
  @test L_o1[2]=="{GRAB THIS TOO}"
  @test L_o2[2]=="GRAB THIS TOO"
  @test L_o3[2]=="{GRAB THIS TOO}"
  @test L_o4[2]=="GRAB THIS TOO"
  @test L_o5[1]==""
end

@testset "Get level flat practical" begin
  s_i = ""
  s_i = string(s_i, "\n", "Some text")
  s_i = string(s_i, "\n", "function myfunc()")
  s_i = string(s_i, "\n", "  more stuff")
  s_i = string(s_i, "\n", "  if something")
  s_i = string(s_i, "\n", "    print('something')")
  s_i = string(s_i, "\n", "  else")
  s_i = string(s_i, "\n", "    print('not something')")
  s_i = string(s_i, "\n", "  end")
  s_i = string(s_i, "\n", "  more stuff")
  s_i = string(s_i, "\n", "end")
  s_i = string(s_i, "\n", "more text")
  L_o = get_level_flat(s_i, ["function ", "if "], [" end", "\nend"])

  s_o = ""
  s_o = string(s_o,       "function myfunc()")
  s_o = string(s_o, "\n", "  more stuff")
  s_o = string(s_o, "\n", "  if something")
  s_o = string(s_o, "\n", "    print('something')")
  s_o = string(s_o, "\n", "  else")
  s_o = string(s_o, "\n", "    print('not something')")
  s_o = string(s_o, "\n", "  end")
  s_o = string(s_o, "\n", "  more stuff")
  s_o = string(s_o, "\n", "end")

  @test L_o[1]==s_o
end

@testset "Get level complex" begin
  s_i = ""
  s_i = string(s_i, "\n", "Some text")
  s_i = string(s_i, "\n", "if something")
  s_i = string(s_i, "\n", "  function myfunc()")
  s_i = string(s_i, "\n", "    more stuff")
  s_i = string(s_i, "\n", "    if something")
  s_i = string(s_i, "\n", "      print('something')")
  s_i = string(s_i, "\n", "    else")
  s_i = string(s_i, "\n", "      print('not something')")
  s_i = string(s_i, "\n", "    end")
  s_i = string(s_i, "\n", "    for something")
  s_i = string(s_i, "\n", "      print('something')")
  s_i = string(s_i, "\n", "    else")
  s_i = string(s_i, "\n", "      print('not something')")
  s_i = string(s_i, "\n", "    end")
  s_i = string(s_i, "\n", "    more stuff")
  s_i = string(s_i, "\n", "  end")
  s_i = string(s_i, "\n", "end")
  s_i = string(s_i, "\n", "more text")

  word_boundaries_left = ["\n", " ", ";"]
  word_boundaries_right = ["\n", " ", ";"]
  word_boundaries_right_if = [" ", ";"]

  FS_outer = FlagSet(
    Flag("function", word_boundaries_left, word_boundaries_right),
    Flag("end",      word_boundaries_left, word_boundaries_right)
  )

  FS_inner = [
  FlagSet(
    Flag("if",       word_boundaries_left, word_boundaries_right_if),
    Flag("end",      word_boundaries_left, word_boundaries_right)
  ),
  FlagSet(
    Flag("for",      word_boundaries_left, word_boundaries_right),
    Flag("end",      word_boundaries_left, word_boundaries_right)
  )]

  L_o = get_level(s_i, FS_outer, FS_inner)

  s_o = ""
  s_o = string(s_o,       " function myfunc()") # The extra space is due to the left word boundary of the function...
  s_o = string(s_o, "\n", "    more stuff")
  s_o = string(s_o, "\n", "    if something")
  s_o = string(s_o, "\n", "      print('something')")
  s_o = string(s_o, "\n", "    else")
  s_o = string(s_o, "\n", "      print('not something')")
  s_o = string(s_o, "\n", "    end")
  s_o = string(s_o, "\n", "    for something")
  s_o = string(s_o, "\n", "      print('something')")
  s_o = string(s_o, "\n", "    else")
  s_o = string(s_o, "\n", "      print('not something')")
  s_o = string(s_o, "\n", "    end")
  s_o = string(s_o, "\n", "    more stuff")
  s_o = string(s_o, "\n", "  end\n") # The \n is due to the right word boundary of the "end"

  @test L_o[1]==s_o

  # Edge case (no flags in string):
  s_i = ""
  s_i = string(s_i, "\n", "Some text")
  s_i = string(s_i, "\n", "more text")
  L_o = get_level(s_i, FS_outer, FS_inner)
  s_o = ""
  @test L_o[1]==s_o
end

@testset "Remove flat" begin
  s_i1 = "Here is some text, and {THIS SHOULD BE REMOVED}, FeaturedFuncs offers a simple interface..."
  s_o1 = remove_flat(s_i1, ["{"], ["}"])
  s_i2 = "Here is some text, and {THIS SHOULD BE REMOVED}, FeaturedFuncs offers a simple interface..."
  s_o2 = remove_flat(s_i2, ["{"], ["}"], false)
  s_i3 = "Here is some text, and {THIS SHOULD BE REMOVED), FeaturedFuncs} offers a simple interface..."
  s_o3 = remove_flat(s_i3, ["{"], ["}", ")"])
  s_i4 = "Here is some text, and {THIS SHOULD BE REMOVED), FeaturedFuncs} offers a simple interface..."
  s_o4 = remove_flat(s_i4, ["{"], ["}", ")"], false)
  s_i5 = "Here is some text, and THIS SHOULD REMAIN, FeaturedFuncs} offers a simple interface..."
  s_o5 = remove_flat(s_i5, ["{"], ["}"], false)
  s_i6 = "Here is some text, and THIS SHOULD REMAIN, FeaturedFuncs| offers a simple interface..."
  s_o6 = remove_flat(s_i6, ["|"], ["|"], false)
  s_i7 = "Here is some text, and |THIS SHOULD BE REMOVED|, FeaturedFuncs offers a simple interface..."
  s_o7 = remove_flat(s_i7, ["|"], ["|"], false)
  s_doc = "\"\"\""
  s_i8 = "Here is some text-"*s_doc*" docs... "*s_doc*", and-"*s_doc*"more docs"*s_doc*", FeaturedFuncs offers a simple interface..."
  s_o8 = remove_flat(s_i8, [s_doc], [s_doc], true)
  s_doc = "\"\"\""
  s_i9 = "Here is some text-"*s_doc*" docs... "*s_doc*", and-"*s_doc*"more docs"*s_doc*", and "*s_doc*"more docs"*s_doc*", FeaturedFuncs offers a simple interface..."
  s_o9 = remove_flat(s_i9, [s_doc], [s_doc], true)

  @test s_o1=="Here is some text, and , FeaturedFuncs offers a simple interface..."
  @test s_o2=="Here is some text, and {}, FeaturedFuncs offers a simple interface..."
  @test s_o3=="Here is some text, and , FeaturedFuncs} offers a simple interface..."
  @test s_o4=="Here is some text, and {), FeaturedFuncs} offers a simple interface..."
  @test s_o5=="Here is some text, and THIS SHOULD REMAIN, FeaturedFuncs} offers a simple interface..."
  @test s_o6=="Here is some text, and THIS SHOULD REMAIN, FeaturedFuncs| offers a simple interface..."
  @test s_o7=="Here is some text, and ||, FeaturedFuncs offers a simple interface..."
  @test s_o8=="Here is some text-, and-, FeaturedFuncs offers a simple interface..."
  @test s_o9=="Here is some text-, and-, and , FeaturedFuncs offers a simple interface..."
end


text:
"""
    StaticStartRange{T,B,E,L}

StaticRange subtype with parametric parameters for the first and last parts of a range.
"""
abstract type StaticStartRange{T,B,E,L} <: StaticRange{T,L} end

@pure first(r::StaticStartRange{T,B}) where {T,B<:SVal} = B()
@pure first(::Type{<:StaticStartRange{T,B}}) where {T,B<:SVal} = B()

@pure last(::StaticStartRange{T,B,E}) where {T,B,E<:SVal} = E()
@pure last(::Type{<:StaticStartRange{T,B,E}}) where {T,B,E<:SVal} = E()

@inline first(r::StaticStartRange{T,B}) where {T,B<:BaseNumber} = getfield(r, :start)::B
@inline first!(r::StaticStartRange{T,B}, start::B) where {T,B<:BaseNumber} = setfield!(r, :start, start)

@inline last(r::StaticStartRange{T,B,E}) where {T,B,E<:BaseNumber} = getfield(r, :stop)::E
@inline last!(r::StaticStartRange{T,B,E}, stop::E) where {T,B,E<:BaseNumber} = setfield!(r, :stop, stop)


text:
@testset "make.jl: load_config" begin
    configfile = joinpath(@__DIR__, "ABS.yml")
    config = load_config(configfile)
    @test config["System"]["label"] == "AB/S"
end

@testset "make.jl: make" begin
    configfile = joinpath(@__DIR__, "ABS.yml")
    config = load_config(configfile)
    abs = make(config)
    @test systemtype(abs) == PolymerSolution()
    chainAB = abs.components[1].molecule
    # @test chaintype(chainAB) == LinearArchitecture()
end


text:
function __init__()
    RLDatasets.d4rl_init()
    RLDatasets.d4rl_pybullet_init()
end


text:
using SpecialFunctions: gamma
using GSL: sf_gamma_inc_P
incomplete_gamma = sf_gamma_inc_P

"""
The displaced Poisson distribution (a.k.a the hyper-Poisson distribution).

Ref: https://www.jstor.org/stable/2283992
"""
struct DisplacedPoisson{T<:Real} <: Distributions.DiscreteUnivariateDistribution
    λ::T
    r::T
    function DisplacedPoisson{T}(λ, r) where {T<:Real}
        @assert λ > 0 "λ is not positve"
        @assert r >= 0 "r is not non-negative"
        return new(λ, r)
    end
end
DisplacedPoisson(λ::T, r::T) where {T<:Real} = DisplacedPoisson{T}(λ, r)

function pdf(dp::DisplacedPoisson, k::Int)
    if k < dp.r
        return 0
    else
        p = dp.λ^(k - dp.r) * exp(-dp.λ) / gamma(k - dp.r + 1)
        diff = ceil(dp.r) - dp.r
        # incomplete_gamma(a, b) is undefined for a = 0
        if !iszero(diff)
            return p / incomplete_gamma(diff, dp.λ)
        else
            return p
        end
    end
end

function rand(dp::DisplacedPoisson{T}) where {T}
    u = rand(T)
    cm = 0
    k = -1
    while cm < u
        k += 1
        cm += pdf(dp, k)
    end
    return k
end

function mode(dp::DisplacedPoisson)
    k = 0
    p = pdf(dp, k)
    while true
        p_next = pdf(dp, k + 1)
        if p_next < p
            break
        end
        p = p_next
        k += 1
    end
    return k
end

minimum(dp::DisplacedPoisson) = ceil(Integer, dp.r)
maximum(dp::DisplacedPoisson) = Inf

Base.show(io::IO, dp::DisplacedPoisson) =
    print(io, "DisplacedPoisson($(round(dp.λ; sigdigits=3)), $(round(dp.r; sigdigits=3)))")