text:
@params mutable struct AdaptiveSIMP{T} <: AbstractSIMP
    simp::SIMP{T}
    pstart::T
    pfinish::T
    reuse::Bool
    result::SIMPResult{T}
    ftol_cont
    Δp::T
    maxiter::Int
    innerpolynomial::PolynomialFit{T}
    prev_f::T
    prev_dfdp::T
    prev_ftol::T
    toltrace::AbstractVector{T}
    progresstrace::AbstractVector{T}
end
GPUUtils.whichdevice(a::AdaptiveSIMP) = whichdevice(a.simp)

function AdaptiveSIMP(simp::SIMP{T}, ::Type{Val{filtering}}=Val{false};
    pstart = T(1),
    pfinish = T(5),
    ftol_cont = ExponentialContinuation(T(0.01), T(4log(0.7)), T(0), 100, T(1e-4)),
    Δp = T(0.1),
    reuse = true,
    maxiter = 500,
    ratio = 0.7) where {T, filtering}

    ncells = getncells(simp.optimizer.obj.problem)
    result = NewSIMPResult(T, simp.optimizer, ncells)
    innerpolynomial = PolynomialFit{T}(3, ratio)

    return AdaptiveSIMP(simp, pstart, pfinish, reuse, result, ftol_cont, Δp, maxiter, innerpolynomial, T(NaN), T(NaN), T(0), T[], T[])
end

function (asimp::AdaptiveSIMP)(x0=asimp.simp.optimizer.obj.solver.vars)
    @unpack pstart, pfinish, Δp = asimp
    @unpack obj = asimp.simp.optimizer
    T = eltype(x0)

    p = pstart
    workspace = setup_workspace(asimp, x0, p)
    @unpack primal_data, options, convstate = workspace
    steps = round(Int, (pfinish - pstart)/Δp)

    maxiter = options.maxiter
    fevals1, fevals2, fevals = obj.fevals, obj.fevals, obj.fevals
    for i in 1:steps
        obj.fevals >= asimp.maxiter && break
        workspace.outer_iter = workspace.iter = 0
        options.maxiter = 1

        f_x_previous = primal_data.f_x_previous
        fevals2 = fevals1
        fevals1 = fevals

        workspace.converged = false
        if i == steps
            obj.reuse = false
        else
            obj.reuse = asimp.reuse
        end

        setpenalty!(asimp.simp, p + Δp)
        simp_results = _innersolve!(asimp, workspace)

        if workspace.converged
            if obj.reuse
                primal_data.f_x = primal_data.f_x_previous
                primal_data.f_x_previous = f_x_previous
                primal_data.x .= primal_data.x1
                primal_data.x1 .= primal_data.x2
            end
        else
            obj.reuse = false
            while !convstate.converged && obj.fevals < asimp.maxiter && workspace.iter < maxiter
                options.maxiter += 1
                _innersolve!(asimp, workspace)
            end
        end
        fevals = obj.fevals
        x_hist = asimp.simp.optimizer.obj.topopt_trace.x_hist
        x = x_hist[fevals]
        x1 = x_hist[fevals1]
        x2 = x_hist[fevals2]
        same = true
        for _x in (x1, x2)
            for j in 1:length(x)
                if round(x[j]) != round(_x[j])
                    same = false
                    break
                end
            end
        end
    end
    options.maxiter = maxiter
    asimp.result = asimp.simp.result

    return asimp.result
end

function setup_workspace(asimp::AdaptiveSIMP, x0::AbstractArray{T}, p) where T
    @unpack model, s_init, s_decr, s_incr, optimizer, suboptimizer, dual_caps, obj, constr, x_abschange, x_converged, f_abschange, f_converged, g_residual, g_converged = asimp.simp.optimizer

    @unpack innerpolynomial, simp = asimp
    # Set penalty as starting penalty
    setpenalty!(simp, p)
    simp.optimizer.obj.reuse = false

    # Does the first function evaluation
    # Number of function evaluations is the number of iterations plus 1
    workspace = MMA.Workspace(model, x0, optimizer, suboptimizer, options = MMA.Options(s_init=s_init, s_incr=s_incr, s_decr=s_decr, dual_caps=dual_caps))
    # Record the first value in the polynomial fit struct
    newvalue!(innerpolynomial, obj.fevals, workspace.f_x)

    # Initial workspace.iter is 0
    # Tolerance is pretty low
    MMA.ftol!(workspace.model, asimp.ftol_cont(0))
    MMA.grtol!(workspace.model, zero(T))
    MMA.xtol!(workspace.model, zero(T))

    return workspace
end

function _innersolve!(asimp::AdaptiveSIMP{T}, workspace::MMA.Workspace) where T
    @unpack innerpolynomial, pstart, pfinish, simp = asimp
    @unpack optimizer, penalty = asimp.simp
    @unpack p = penalty
    @unpack obj = optimizer
    @unpack xmin = obj.solver
    @unpack model = workspace

    local simp_results

    while !workspace.converged && obj.fevals < asimp.maxiter && workspace.iter < model.maxiter[]
        if order(innerpolynomial) < maxorder(innerpolynomial)
            ftol = asimp.ftol_cont(0)
        else
            ftol = get_tol(asimp)
        end
        push!(asimp.toltrace, ftol)
        MMA.ftol!(model, ftol)
        simp_results = simp(workspace)
        newvalue!(innerpolynomial, obj.fevals, workspace.f_x)
    end

    return simp_results
end

function get_tol(asimp::AdaptiveSIMP{T}) where {T}
    p = getpenalty(asimp.simp).p
    obj = asimp.simp.optimizer.obj
    poly = asimp.innerpolynomial
    # If not converged yet, time to do tolerance adaptation
    # Fit the quadratic
    innerpolycoeffs = solve!(poly)

    # Set the gradient tolerance as the speed/innerpolynomial''
    if p >= asimp.pfinish - sqrt(eps(T))
        progress = one(T)
    else
        # Get the second derivative of the quadratic
        zero_derivative = poly(obj.fevals)
        first_derivative = poly'(obj.fevals)
        second_derivative = poly''(obj.fevals)
        third_derivative = poly'''(obj.fevals)

        f = first_derivative / zero_derivative
        s = second_derivative / zero_derivative
        t = third_derivative / zero_derivative
        ds = abs(f) + abs(s) + abs(t)
        #_progress = log(1 + 1/max(ds, sqrt(eps(T))))
        #maxprogress = log(1 + 1/sqrt(eps(T)))
        #progress = clamp(_progress/maxprogress, (p - asimp.pstart)/(asimp.pfinish - asimp.pstart), 1)
        #progress = _progress/maxprogress
        progress = (p - asimp.pstart)/(asimp.pfinish - asimp.pstart)
    end
    push!(asimp.progresstrace, progress)
    ftol = asimp.ftol_cont(progress)

    return ftol
end

#=
Utilities.getpenalty(obj) = obj.solver.penalty

dFdp(obj::Objective{<:ComplianceFunction}) = dFdp(obj, getpenalty(obj))
function dFdp(obj::Objective{<:ComplianceFunction{T}}, penalty::PowerPenalty{T}) where T
    p = penalty.p
    cell_comp = obj.cell_comp
    vars = obj.solver.vars
    xmin = obj.solver.xmin

    return -mapreduce(e -> power_penalty_dFdpterm(vars[e], xmin, cell_comp[e], p), +, zero(T), 1:length(cell_comp))
end

function power_penalty_dFdpterm(var, xmin, cellcomp, p)
    return cellcomp*log(var)*var^p
end
=#


text:
# This file assumes that the directories are organized as follows:
#
# ./SomeDirectory/RRTMGP.jl/
# ./SomeDirectory/rte-rrtmgp/
#

include("ConfigRRTMGP.jl")

cd(joinpath(root_dir, "rte-rrtmgp", "build")) do
    rm("Makefile.libs"; force = true)
    rm("Makefile.conf"; force = true)
    rm("Makefile.rules"; force = true)
end

cd(joinpath(root_dir, "rte-rrtmgp", "build")) do
    run(`make clean`)
end

cd(joinpath(root_dir, "rte-rrtmgp")) do
    rm(".DS_store"; force = true)
end
cd(joinpath(root_dir, "rte-rrtmgp", "examples")) do
    rm(".DS_store"; force = true)
end

cd(rte_rrtmgp_rfmip_clear_sky) do
    conds_file = joinpath(".", "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc")
    rm(conds_file; force = true)
    rm(".DS_store"; force = true)
    rm("generate-output-file-templates.py"; force = true)
    rm("rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force = true)
    rm("rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force = true)
    rm("rrtmgp_rfmip_lw"; force = true)
    rm("rrtmgp_rfmip_sw"; force = true)
    rm("rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force = true)
    rm("rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force = true)
    run(`make clean`)
end

cd(rte_rrtmgp_rfmip_all_sky) do
    run(`make clean`)
end

ext: jl

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
import Term: Panel, RenderableText, TextBox
import Term: split_lines
import Term: install_stacktrace

install_stacktrace()

println(Panel("\n", width = 120, style = "dim blue", title_style = "bold bright_yellow"))

# ----------------------------- basic and nested ----------------------------- #
print(
    Panel(
        """
        [green]This is a multiline text[/green]
        [on_black]This line is [bold] different![/bold][/on_black]""",
        subtitle = "aa",
        subtitle_justify = :right,
        style = "red",
        title = "From text",
        width = 44,
    ),
)
println("."^44)
print("\n\n")

w = 32
println(Panel("."^w, title = "TITLE", title_justify = :center, style = "green", width = 2w))

print(
    Panel(
        RenderableText("""
        [green]This is a multiline text[/]
        [on_black]This line is [bold]different!"""),
        style = "blue dim",
        title = "From Segment",
        width = 50,
        justify = :center,
        title_style = "red bold",
    ),
)
println("."^50)
print("\n\n")

println(
    Panel(
        Panel("[bold white]Title panel!!", style = "dim", width = 88),
        Panel(
            RenderableText(
                """
[green]This is a multiline text[/]
[on_black]This line is [bold]different!
""",
                width = 100,
            ),
            style = "blue dim",
            title = "From Segment",
            justify = :center,
            title_style = "red bold",
            subtitle = "test",
            subtitle_justify = :right,
            width = 110,
        ),
        justify = :center,
        title = "created with Term",
        title_style = "gray62",
    ),
)
print("\n\n")

# ---------------------------------- textbox --------------------------------- #
print(TextBox(join("TEST"^4, "sdfs"^8, "\n"), title = "test box!"))
print("\n\n")

print(
    TextBox(join("X"^26, "y"^8), title = "test box!", title_style = "bold red", width = 44),
)
println("."^44)
print("\n\n")

# ---------------------------- textbox with colors --------------------------- #

print(TextBox("[red]TEST[/red]aaa"^10, title = "test box!"; width = 44))
println("."^44)
print("\n\n")

print(
    TextBox(
        "[red]red [blue]blue[/blue][green]green[/green]red[/red]"^50,
        title = "test box!",
        title_style = "bold red",
        width = 44,
    ),
)
println("."^44)
print("\n\n")


text:
using Plots
using LLLplus
using Printf
using Random
using LinearAlgebra
using DoubleFloats
using Quadmath

include("lrtest.jl")

# getIntType is used to indicate what type we want the unimodular integer
# matrix to have. We need to add methods of getIntType to handle the
# DoubleFloats and Quadmath, since they are external packages that LLLplus
# doesn't know about. So yes, we do a bit of type piracy.
import LLLplus.getIntType
getIntType(Td::Type{Tr}) where {Tr<:Float128} = Int128
getIntType(Td::Type{Tr}) where {Tr<:Double64} = Int128

lrtest(40,2 .^[7],[100],[Int32,Int64,Int128,Float32,Float64,Float128,Double64,BigInt,BigFloat],"rand")
savefig("perfVsDataType.png")

lrtest(40,2 .^[0:8;],[1],[Float64],"randn")
savefig("perfVsNfloat64.png")


text:
@symbol_func function alpha_I_star(cur_reactor::AbstractReactor)
  cur_scaling = cur_reactor.mode_scaling

  cur_alpha_star = cur_scaling[:I]

  cur_alpha_star += alpha_n_star(cur_reactor)

  cur_alpha_star
end


text:
"""
    DAS
Abstract supertype of all Discrete Abelian Symmetries
"""
abstract type DAS end

"""
    DASCharges
Abstract supertype of all collections of Charges of Discrete Abelian Symmetries
"""
abstract type DASCharges end

"""
    DASCharge
Abstract supertype of all Charges of Discrete Abelian Symmetries
"""
abstract type DASCharge end
Base.hash(a::DASCharge, h::UInt) = hash(a.ch, h)

"""
    ⊕(a::T, bs::T...) where {T<:Union{DASCharge, DASCharges}}
returns the result of fusing one or more `DASCharge` or `DASCharges` together.
Fusing `DASCharges` yields a `DASCharges` that holds all elements that result
from fusing the `DASCharges`.
"""
⊕(a::DASCharges) = a
⊕(a::T, b::T, c::T...) where {T<:DASCharges} where N = ⊕((a ⊕ b),c...)
const oplus = ⊕

function Base.iterate(a::DASCharges, i = 1)
    i <= length(a) && return (a[i], i+1)
    return nothing
end
Base.IteratorSize(::DASCharges) = Base.HasLength()
Base.IteratorEltype(::DASCharges) = Base.HasEltype()

Base.iszero(a::DASCharge) = iszero(a.ch)
Base.:-(a::T) where {T<:DASCharge} = InOut(-1) ⊗ a
Base.:(==)(a::DASCharge, b::DASCharge) = a.ch == b.ch
Base.isless(a::DASCharge, b::DASCharge) = isless(a.ch, b.ch)

"""
    chargeindex(ch::DASCharge, chs::DASCharges)
returns the index i of ch in chs

# Example
```julia-repl
julia> ch = U1Charge(0); chs = U1Charges(-1:1);
julia> chargeindex(ch, chs)
2
julia> chs[2] == ch
true
```
"""
function chargeindex end

include("dastensor/U1.jl")
include("dastensor/ZN.jl")
include("dastensor/NDAS.jl");

#Sector
"""
    DASSector{N,T}
DASSectors are a configuration of charges that are allowed under a given
symmetry and index degeneracy spaces in DASTensors{T,N}.

# Example
```julia-repl
julia> DASSector(U1Charge(1), U1Charge(2))
DASSector(U1Charge(1), U1Charge(2))
```
"""
struct DASSector{N,T}
    chs::NTuple{N,T}
    DASSector(chs::T...) where T = new{length(chs),T}(chs)
    DASSector{N,T}(chs::T...) where {N,T} = new{N,T}(chs)
end

Base.show(io::IO, s::DASSector{N,S}) where {N,S}  = print(io, "DASSector", s.chs)

function Base.iterate(a::DASSector{N}, i = 1) where N
    i <= length(a) && return (a[i], i+1)
    return nothing
end
Base.IteratorSize(::DASSector) = Base.HasLength()
Base.IteratorEltype(::DASSector) = Base.HasEltype()

Base.length(s::DASSector{N}) where N = N
Base.eltype(::DASSector{N,T}) where {N,T} = T
Base.hash(a::DASSector, h::UInt) = hash(a.chs, h)
⊕(s1::DASSector{N,T}, s2::DASSector{N,T}) where {N,T} = DASSector(map(⊕, s1.chs, s2.chs))
Base.getindex(s::DASSector,i::Int) = s.chs[i]
permute(s::DASSector{N,T}, perm) where {N,T}  = DASSector{N,T}(TT.permute(s.chs,perm)...)

"""
    charge(a::DASSector)
returns the charge which is calculated as minus the sum of all charges it contains.
"""
charge(s::DASSector) = -reduce(⊕, s.chs)
Base.isless(a::T, b::T) where {T<:DASSector} = a.chs < b.chs
Base.:(==)(a::T, b::T) where {T<:DASSector} = hash(a) == hash(b)
Base.:(==)(a::DASSector, b::DASSector) = false

Base.vcat(s::DASSector) = s
Base.vcat(s1::DASSector, s2::DASSector, s3::DASSector...) =
    vcat(vcat(s1,s2),s3...)
Base.getindex(s::DASSector{N,T}, i::NTuple{M}) where {N,T,M} =
    DASSector{M,T}(TT.getindices(s.chs,i)...)
Base.vcat(s1::DASSector{N1,T}, s2::DASSector{N2,T}) where {T,N1,N2} =
    DASSector{N1+N2,T}(TT.vcat(s1.chs, s2.chs)...)
deleteat(s1::DASSector{N,T}, i::NTuple{M}) where {N,T,M} =
    DASSector{N-M,T}(TT.deleteat(s1.chs,i)...)

#InOut
"""
    InOut{N}
InOut describes whether representations of the DAS act on an index of a
DASTensor{T,N} directly or via their dual.
InOut(1,1,1,-1) can be read as the first three indices corresponding to incoming,
the last as an outgoing index w.r.t the group action.
"""
struct InOut{N}
    v::NTuple{N,Int}
    function InOut(v::Int...)
        all(in((1,-1)), v) || throw(ArgumentError("elements must either be 1 or -1"))
        new{length(v)}(v)
    end
end
Base.show(io::IO, s::InOut) = print(io, "InOut", s.v)

function Base.iterate(a::InOut{N}, i = 1) where N
    i <= N && return (InOut(a.v[i]), i+1)
    return nothing
end
Base.IteratorSize(::InOut) = Base.HasLength()
Base.IteratorEltype(::InOut) = Base.HasEltype()

Base.length(::InOut{L}) where L = L
Base.eltype(::InOut) = InOut{1}
⊗(a::InOut{N}, b::T) where {T<:DASSector,N} = T(map(⊗, ntuple(i -> a[i] ,Val{N}()) , b.chs)...)
⊗(a::InOut{1}, b::T) where {T<:DASCharge} = T(a.v[1] * b.ch)
⊗(a::InOut{1}, b::T) where {T<:NDASCharge} = NDASCharge(map(⊗(a),b.ch)...)
⊗(a::InOut{M}, b::NTuple{M,DASCharges}) where {M} = ntuple(i -> a[i] ⊗ b[i], Val(M))
⊗(a::InOut{1}, b::DASCharges) = a.v == (1,) ? b : inv(b)

⊗(io::InOut) = x -> io ⊗ x
Base.inv(io::InOut) = InOut((-1 .* io.v)...)
Base.getindex(a::InOut, i::Int) = InOut(a.v[i])
Base.getindex(a::InOut, i::Union{NTuple{L},UnitRange}) where L = InOut(TT.getindices(a.v,i)...)
Base.lastindex(::InOut{L}) where L = L
Base.vcat(a::InOut) = a
Base.vcat(a::InOut, b::InOut) = InOut(a.v..., b.v...)
Base.vcat(a::InOut, b::InOut, c::InOut...) = vcat(vcat(a,b),c...)

#sector generators
"""
    allsectors(chs)
returns a generator that generates all possible combinations of charges in
`chs` wrapped in a `DASSector`.

#Example
```julia-repl
julia> allsectors((U1Charges(-1:1), U1Charges(4:5))) |> collect
3×2 Array{DASSector{2,U1Charge},2}:
 DASSector(U1Charge(-1), U1Charge(4))  DASSector(U1Charge(-1), U1Charge(5))
 DASSector(U1Charge(0), U1Charge(4))   DASSector(U1Charge(0), U1Charge(5))
 DASSector(U1Charge(1), U1Charge(4))   DASSector(U1Charge(1), U1Charge(5))
```
"""
allsectors(a::NTuple{N,DASCharges}) where N = (DASSector(s...) for s in Iterators.product(a...))

"""
    covariantsectors(chs, io[, ch = zero(chs)])
returns  all sectors in allsectors(io ⊗ chs) that have total charge ch
"""
covariantsectors(charges, io::InOut, ch = zero(eltype(first(charges)))) =
    Iterators.filter(==(ch) ∘ charge ∘ ⊗(io), allsectors(charges))

"""
    covariantsectors(chs, io)
returns all sectors in  allsectors(io ⊗ chs) that have total charge zero.
"""
invariantsectors(charges, io::InOut) = covariantsectors(charges, io)

# TENSOR
mutable struct DASTensor{T,N,SYM,CHARGES,SIZES,CHARGE} <: AbstractTensor{T,N}
    chs     ::NTuple{N, CHARGES}
    dims    ::NTuple{N, SIZES}
    io      ::InOut{N}
    ch      ::Union{CHARGE, Missing}
    tensor  ::Dict{DASSector{N,CHARGE},Array{T,N}}
end

function Base.show(io::IO, A:: DASTensor{T,N,SYM}) where {T,N,SYM}
    print(io, "DASTensor{",T,",",N,",",SYM,"}\n")
    print(io, "charges: ", charges(A), "\n")
    print(io, "charge: ",  charge(A), "\n")
    print(io, "sizes: ",   sizes(A), "\n")
    print(io, "in/out: ",  in_out(A), "\n")
    print(io, "Tensors ",  eltype(A), "\n")
end

##Functions
function DASTensor{T,N}(sym::ST, charges, dims, io,
    ch = zero(chargetype(sym)), dict = Dict()) where {T,N,ST <: DAS}
    length.(dims) == length.(charges) || throw(ArgumentError("incompatible dims/charges"))
    CHARGES = chargestype(sym)
    CHARGE  = chargetype(sym)
    SIZES   = Vector{Int}
    DASTensor{T,N,sym,CHARGES,SIZES,CHARGE}(charges,dims,io,ch,dict)
end

"""
    charges(A::DASTensor[,i])
returns the charges of `A`.
If `i` is specified as either `Int` or `Tuple`, returns only the charges of the
indices in `i`.
"""
charges(A::DASTensor)        = A.chs
charges(A::DASTensor,i::Int) = A.chs[i]
charges(A::DASTensor,i)      = TT.getindices(A.chs,i)

"""
    sizes(A::DASTensor[,i])
returns the sizes of `A` as a tuple of vectors `v` such that the degeneracy
space associated with a charge `ch` has size `v[chargeindex(ch, chs)]` where
`chs` is the `DASCharges` associated with the specified leg.
If `i` is specified as either `Int` or `Tuple`, returns only the charges of the
indices in `i`.
"""
sizes(A::DASTensor)        = A.dims
sizes(A::DASTensor,i::Int) = A.dims[i]
sizes(A::DASTensor,i)      = TT.getindices(A.dims,i)

"""
    in_out(A::DASTensor[,i])
returns the `InOut` of `A` which specifies the action of the symmetry group
on the corresponding leg.
If `i` is specified as either `Int` or `Tuple`, returns only the charges of the
indices in `i`.
"""
in_out(A::DASTensor)   = A.io
in_out(A::DASTensor,i) = A.io[i]

"""
    tensor(A::DASTensor[,i])
returns a dictionary of `DASSectors` and their associated degeneracy spaces.
"""
tensor(A::DASTensor)    = A.tensor
Base.getindex(A::DASTensor, i) = A.tensor[i]
Base.setindex!(A::DASTensor{T,N}, t, i::DASSector{N}) where {T,N} =
    (A.tensor[i] = t; A)

Base.values(A::DASTensor) = values(A.tensor)
Base.keys(A::DASTensor) = keys(A.tensor)
Base.haskey(A::DASTensor,k) = haskey(A.tensor,k)
"""
    symmetry(A::DASTensor)
returns the symmetry of a tensor `A`
"""
symmetry(A::DASTensor{<:Any,<:Any,SYM}) where {SYM} = SYM
Base.eltype(A::DASTensor{T,N}) where {T,N} = T
chargetype(A::DASTensor{<:Any,<:Any,<:Any,<:Any,<:Any,CHARGE}) where {CHARGE} = CHARGE
chargestype(A::DASTensor{<:Any,<:Any,<:Any,CHARGES}) where {CHARGES} = CHARGES
dimstype(A::DASTensor{<:Any,<:Any,<:Any,<:Any,SIZES}) where {SIZES} = SIZES
Base.ndims(A::DASTensor{T,N}) where {T,N} = N

"""
    setcharges!(A::DASTensor, chs)
set the charges of `A` to be `chs` where the latter is a tuple of `DASCharges`.
"""
setcharges!(A::DASTensor, charges) = (A.chs = charges; A)
"""
    setsizes!(A::DASTensor, s)
set the sizes of `A` to be `s` where the latter is a tuple of `Vector{Int}`.
"""
setsizes!(A::DASTensor, sizes)     = (A.dims   = sizes;   A)
"""
    setin_out!(A::DASTensor, io)
set the `InOut` of `A` to be `io` where the latter is a `InOut`.
"""
setin_out!(A::DASTensor, in_out)   = (A.io  = in_out;  A)
settensor!(A::DASTensor, tensor)   = (A.tensor  = tensor;  A)

"""
    charge(a::DASTensor)
returns the charge of a tensor.
"""
charge(A::DASTensor) = A.ch
setcharge!(A::DASTensor, ch) = (A.ch = ch; A)

"""
    isinvariant(a::DASTensor)
return true if `charge(a)` is `zero`.
"""
isinvariant(a::DASTensor) = iszero(charge(a))

function chargesize(charge::DASCharge, charges, dims)
    charge in charges || throw(
        ArgumentError(string(charge, " is not in ", charges)))
    dims[chargeindex(charge, charges)]
end
chargesize(A::DASTensor, i, charge)          = chargesize(charge, charges(A,i), sizes(A,i))
degeneracysize(sector, charges, dims)        = map(chargesize, sector, charges, dims)

"""
    initwith!(A::DASTensor{T}, fun [,ch])
modifies `A` such that each sector with charge `ch` (default=zero) is (independently)
set to `fun(T, dims...)` where `dims` is the size of the degeneracy space for the sector.
"""
function initwith!(A::DASTensor{T,N}, fun) where {T,N}
    empty!(A.tensor)
    dims, chs = sizes(A), charges(A)
    for k in covariantsectors(chs, in_out(A), charge(A))
        A.tensor[k] = fun(T, degeneracysize(k, chs, dims)...)
    end
    return A
end

function Base.rand(::Type{DASTensor{T,N}}, sym, charges, dims, io,
                    ch = zero(eltype(first(charges)))) where {T,N}
    A = DASTensor{T,N}(sym, charges, dims, io)
    initwithrand!(A, ch)
end

"""
    initwithzero!(A::DASTensor)
construct all valid sectors  in `A` and initialize their
degeneracy spaces with zeros.
"""
initwithzero!(A::DASTensor) = initwith!(A, zeros)

"""
    initwithrand!(A::DASTensor)
construct all valid sectors in `A` and initialize their
degeneracy spaces with `rand`.
"""
initwithrand!(A::DASTensor) = initwith!(A, rand)

function Base.convert(::Type{Array{S}}, A::DASTensor{T,N,<:Any,CHARGES,SIZES,CHARGE}) where
        {S,T,N,CHARGES,SIZES,CHARGE}
    iszero(N) && return convert(Array{S}, first(values(tensor(A))))
    cumdims = ntuple(i -> prepend!(cumsum(sizes(A,i)),0) , N)
    array = zeros(S, last.(cumdims)...)
    for (sector, degeneracy) in tensor(A)
        inds   = map(chargeindex, sector.chs, charges(A))
        ranges = map((i, cdims) -> cdims[i]+1 : cdims[i+1], inds, cumdims)
        array[ranges...] = degeneracy
    end
    array
end

toarray(a::DASTensor{T}) where {T} = convert(Array{T}, a)

Base.convert(::Type{DTensor{S}}, A::DASTensor) where S = DTensor(convert(Array{S}, A))
todense(A::DASTensor{T,N}) where {T,N} = convert(DTensor{T},A)

function LA.diag(A::DASTensor{T,2}) where T # TODO!!
    ks = sort!(collect(keys(A)))
    return reduce(vcat, LA.diag(A[k]) for k in ks)
end

#= Equalities =#
function Base.:(==)(A::TT, B::TT) where {TT <: DASTensor{T,N}} where {T,N}
    issimilar(A, B) && _tensor_equ(tensor(A), tensor(B))
end

Base.:(==)(A::DASTensor, B::DASTensor) = false

function Base.:(≈)(A::DASTensor{T,N}, B::DASTensor{T,N}) where {T,N}
    issimilar(A, B) && _tensor_approx(A, B)
end

Base.:(≈)(A::DASTensor, B::DASTensor) = false

function _tensor_equ(A::DASTensor, B::DASTensor)
    keys(A) == keys(B) || return false
    all(k -> A[k] == B[k], keys(A))
end

function _tensor_approx(A::DASTensor, B::DASTensor)
    keys(A) == keys(B) || return false
    all(k -> A[k] ≈ B[k], keys(A))
end

_tensor_equ(A, B) = false

function issimilar(A::DASTensor, B::DASTensor)
    return  charge(A) == charge(B) &&
            charges(A) == charges(B) &&
            sizes(A) == sizes(B) &&
            in_out(A) == in_out(B)
end

Base.:(≈)(A::DASTensor, B) = false
Base.:(≈)(A, B::DASTensor) = false

#= Copy =#
function Base.copy!(dest::TT, source::TT) where {TT <: DASTensor{T,N}} where {T,N}
    setcharges!(dest, charges(source))
    setsizes!(dest, deepcopy(sizes(source)))
    setin_out!(dest, in_out(source))
    settensor!(dest, deepcopy(tensor(source)))
    return dest
end

Base.copy(A::DASTensor) = Base.deepcopy(A)

function Base.copyto!(dest::TT, source::TT) where {TT <: DASTensor}
    setcharges!(dest, charges(source))
    setsizes!(dest, deepcopy(sizes(source)))
    setin_out!(dest, in_out(source))
    setcharge!(dest, charge(source))
    empty!(dest.tensor)
    for (k,v) in source.tensor
        dest.tensor[k] = copy(v)
    end
    return dest
end

Base.similar(A::DASTensor{T,N,SYM}, ::Type{S}) where {T,N,SYM,S} =
    DASTensor{S,N}(SYM, charges(A), deepcopy(sizes(A)), in_out(A),charge(A))

Base.similar(A::DASTensor{T}) where T = similar(A,T)


text:
abstract type                AbstractMarker end
         struct Marker{N} <: AbstractMarker end

const Marker0 = Marker{0} ; const ʊ₀ = Marker0
const Marker1 = Marker{1} ; const ʊ₁ = Marker1
const Marker2 = Marker{2} ; const ʊ₂ = Marker2
const Marker3 = Marker{3} ; const ʊ₃ = Marker3

string(::Type{Marker0}) = "ʊ₀" ; string(x::Marker0) = "ʊ₀"
string(::Type{Marker1}) = "ʊ₁" ; string(x::Marker1) = "ʊ₁"
string(::Type{Marker2}) = "ʊ₂" ; string(x::Marker2) = "ʊ₂"
string(::Type{Marker3}) = "ʊ₃" ; string(x::Marker3) = "ʊ₃"

show(io::IO, ::Type{Marker0}) = print(io, string(Marker0))
show(io::IO, ::Type{Marker1}) = print(io, string(Marker1))
show(io::IO, ::Type{Marker2}) = print(io, string(Marker2))
show(io::IO, ::Type{Marker3}) = print(io, string(Marker3))


text:
using PkgSRA
using Plots, Random
using Flux, DiffEqFlux, OrdinaryDiffEq, Statistics
using DSP
pyplot()
Random.seed!(11)

#####
##### Generate true dynamics matrix and data
#####
# Note: continuous version

# Generate the control signal starts
# include("example_forcing_functions.jl")
include("example_falling_ball.jl")
u0 = u1;
# dat_no_wall = solve_ball_system()
# plot(dat_no_wall[1,:], label="Ball position")

# dat = solve_deformed_ball_system()

U_func(X) = U_func_spring(X);
dat = solve_ball_system(U_func_time_trivial,
                        U_func)
# dat = Array(sol)

# Add a row of constants to the data (not actually needed here)
n, m = size(dat)
dat_and_const = vcat(dat, ones(1,m))

#####
##### Mini-SRA: initial guess for dynamics vs. control
#####
dat_grad = numerical_derivative(dat, ts)
# Plot this numerical derivative vs. the truth
plot(dat_grad[2,:], label="Numerical derivative")
# plot!(core_dyn_true(dat)[2,:], label="True derivative")

U_guess = sra_dmd(dat_and_const, dat_grad,
                    quantile_threshold=0.8)[1]
A, B = dmdc(dat_and_const, dat_grad, U_guess)
# Loop a couple of times
for i in 1:20
    pt = plot(U_guess[2,:], label="Residual")
    # plot!(U_func.(ts), label="True")
    title!("Iteration $i")
    display(pt)
    # Update guess for control signal
    global U_guess = calculate_residual(dat_and_const,
                                        dat_grad,
                                        A)
    sparsify_signal!(U_guess, quantile_threshold=0.8)
    # Update guess for dynamics
    global A, B = dmdc(dat_and_const, dat_grad, U_guess)
end

# Filter the signal to be easier for the NN to learn
# myfilter = digitalfilter(Lowpass(0.5),Butterworth(2));
# U_filter = zeros(size(U_guess))
# for i in 1:size(U_guess, 1)
#     U_filter[i,:] = filtfilt(myfilter, U_guess[i,:]);
# end
# U_filter *= maximum(abs.(U_guess))/maximum(abs.(U_filter))
U_true = zeros(size(dat))
for i in 1:size(dat,2)
    U_true[:,i] = U_func(dat[:,i])
end
# U_true = [U_func(d) for d in dat]
plot(U_true[2,:], label="Control Signal")
plot!(U_guess[2,:], label="Raw Guessed Control Signal")
# plot!(U_filter[2,:], label="Filtered Guessed Control Signal")
U_filter = U_guess

#####
##### Learn a NN that can reconstruct the SRA control signal
#####
# Now, a function of phase space
nn_dim = 24
U_dim = 1 # Expect that the true dimensionality is 1
# Note: need to initialize the first layer so that the whole space is used
ctr_dyn = Chain(Dense(n, nn_dim),
                Dense(nn_dim, nn_dim, σ),
                Dense(nn_dim, U_dim),
                Dense(U_dim,2)) # Equivalent to the B matrix
ps = Flux.params(ctr_dyn)
loss_U_filter() = sum(abs2,ctr_dyn(dat) .- U_filter)

loss_to_use = loss_U_filter

cb_ctr = function (i=2)
    pt = plot(Flux.data(ctr_dyn(dat)[i,:]), label="Prediction")
    plot!(U_filter[i,:], label="Guess for control signals")
    display(pt)
    return pt
end
cb_ctr(2)

# Initial fast learning
tol = 200
dat_dummy = Iterators.repeated((), 100)
opt = ADAM(0.1)
while loss_to_use() > 5*tol
    Flux.train!(loss_to_use, ps, dat_dummy, opt, cb = ()->@show loss_to_use())
    cb_ctr(2)
end
opt = ADAM(0.01)
while loss_to_use() > tol
    Flux.train!(loss_to_use, ps, dat_dummy, opt, cb = ()->@show loss_to_use())
    cb_ctr(2)
end
println("Finished learning control signal")

#####
##### As a sanity check, do the NN syntax with the true dynamics
#####
# B_true = Float32.([0.; 1.])
# function dudt_true(u::AbstractArray,p,t)
#     core_dyn_true(u) + B_true*Float32(U_func(t))
# end
# p = param(Float32[0.0]) # Unused
# prob_true = ODEProblem(dudt_true,u0,tspan,p)
# sol = diffeq_rd(p,prob_true,Tsit5(), saveat=ts)
# plot(Flux.data.(Array(sol))[2,:], label="Ideal Simulation")
# plot!(dat[2,:], label="Data")
# title!("True dynamics")
#
# best_loss = sum(abs2,Flux.data.(Array(sol)) .- dat)
# println("The best loss achievable is $best_loss")

#####
##### DEFINE the unified mixed NODE
#####
# Part 1: Linear dynamics; don't predict time
# Note: use actual constructor, not the convenience one
core_dyn = Dense(param(Float32.(A[:, 1:n])),
                 param(Float32.(A[:, end])),
                 identity)
# Part 2: Neural network for control signal; initialized above
p = param(Float32[1.0; 1.0])

function dudt_(u::AbstractArray,p,t)
    core_dyn(u) + p.*ctr_dyn(u)
end

prob = ODEProblem(dudt_,u0,tspan,p)
# Check the initial solution
# sol = diffeq_rd(p,prob,Tsit5(), saveat=ts)
# plot(Flux.data.(Array(sol))[2,:], label="Simulation")
# plot!(dat[2,:], label="Data")

_x = param(u0) # The initial condition can be learned too

function predict_rd()
    Flux.Tracker.collect(diffeq_rd(p,prob,Tsit5(),u0=_x, saveat=ts))
end
function loss_rd(λ1=0.0001, λ2=0.001)
    return sum(abs2,predict_rd() .- dat) +      # Good fit
            λ1*loss_U_filter() +                # Control signal close to guess
            λ2*sum(abs,ctr_dyn(dat))   # L1 norm on ctr
end

cb_plot = function (i=1)
  # display(loss_rd())
  pt = plot(ts, Flux.data(predict_rd()[i,:]), label="Prediction")
  plot!(ts, dat[i,:], label="Data")
  display(pt)
  return pt
end

#####
##### SOLVE the mixed NODE
#####
ps = params(core_dyn,ctr_dyn,p);

# Display the ODE with the current parameter values.
cb_plot(1)
# Initial quick learning
opt = ADAM(1e-6)
dat_dummy = Iterators.repeated((), 10)
Flux.train!(loss_rd, ps,
            dat_dummy, opt, cb = ()->@show loss_rd())
cb_plot(1)
# Run for a long time
opt = ADAM(1e-9)
tol = 50
dat_dummy = Iterators.repeated((), 5)
while loss_rd() > tol
    Flux.train!(loss_rd, ps,
                dat_dummy, opt, cb = ()->@show loss_rd())
    cb_plot(1)
    break
end

#####
##### Instead: SOLVE the mixed NODE via alternating
#####

# Pure loss
loss_dat() = sum(abs2,predict_rd() .- dat)
loss_ctr(λ=0.1) = loss_dat() + λ*sum(abs,ctr_dyn(dat))

function plot_dat_and_ctr()
    pt1 = cb_plot();
    pt2 = cb_ctr();
    pt = plot(pt1, pt2, layout=@layout [a;b]);
    display(pt)
end

ps_core = params(core_dyn);
ps_ctr = params(ctr_dyn,p);

plot_dat_and_ctr()
dat_dummy = Iterators.repeated((), 10)
tol = 150;
opt = ADAM(1e-10)
while loss_dat() > tol
    println("Learning core dynamics; direct L2 loss")
    Flux.train!(loss_dat, ps_core,
                dat_dummy, opt, cb = ()->@show loss_dat())
    plot_dat_and_ctr()

    println("Learning Control signal; direct L2 loss with L1 penalty")
    Flux.train!(loss_ctr, ps_ctr,
                dat_dummy, opt, cb = ()->@show loss_dat())
    plot_dat_and_ctr()

    # break
end


text:
using LAP_julia, TimerOutputs, DataFrames, TableView, JLD2, FileIO
##
sp_lap_win_sum = DataFrame(
    index = Int[],
    reg_fun = Symbol[],
    img_size = Int[],
    whs = Int[],
    timer = TimerOutput[],
    results = Dict{String,Float64}[]);

window_half_sizes = cat(collect(1:51), collect(61:10:101), collect(141:40:381), dims=1);
img_sizes = [50, 100, 200, 400, 800, 1600];

##
df = sp_lap_win_sum;

let index = 0
    for img_size in img_sizes
        img, imgw, flow = gen_init(:chess, chess_args=[25, img_size/25])

        # limit the max window half size to 1/4 the image
        whs_limit = img_size/4
        whs_modified = filter(x -> x <= whs_limit, window_half_sizes)

        for whs in whs_modified
            for reg_fun in [sparse_lap]
                window = [whs * 2 + 1, whs * 2 + 1]
                # window sum 1
                timer = TimerOutput("reg alg: sp lap")
                flow_est, source_reg, timer, results = test_registration_alg(reg_fun, img, imgw, flow, [whs, window], Dict(:timer => timer), timer=timer, display=false)

                # repeat with the same timer
                if whs <= repeat_thresh
                    for _ in 1:repeat_count
                        flow_est, source_reg, timer, results = test_registration_alg(reg_fun, img, imgw, flow, [whs, window], Dict(:timer => timer), timer=timer, display=false)
                    end
                end

                index = index + 1
                println("at index: ", index, " img_size:", img_size, " whs: ", whs)
                push!(df, Dict(:index => index,
                               :reg_fun => Symbol(reg_fun),
                               :img_size => img_size,
                               :whs => whs,
                               :timer => timer,
                               :results => results))
            end
        end
    end
end

# df

# win_sum1_df = df[(df.reg_fun .== :sparse_lap_win_sum1), :];
# win_sum3_df = df[(df.reg_fun .== :sparse_lap), :]; # [1, :timer]
#
# TimerOutputs.ncalls(df[(df.reg_fun .== :sparse_lap), :][1, :timer][sections[1]][sections[2]][sections[3]])
#
# get_timer(win_sum3_df[1, :timer], sections)

import LAP_julia: sparse_lap

show(win_sum3_df[1, :timer])

sections = ["reg alg: sp lap", "sparse lap", "prepare A and b"]
sections = ["reg alg: sp lap", "sparse lap", "filtering"]

fig, axs = subplots(6, sharex=true, figsize=(10,10))
for (img_size, ax) in zip(img_sizes[1:end-1], axs)
    ax.set_title("Image Size: " * string(img_size))

    n = count(k -> k < (img_size/4), window_half_sizes)
    timers3 = df[(df.reg_fun .== :sparse_lap) .& (df.img_size .== img_size), :][:, :timer][1:n]

    y3 = map(x -> get_avg_time(get_timer(x, sections)), timers3)
    x = filter(x -> x < (img_size/4), window_half_sizes)[1:n]

    xscale("log")
    yscale("log")
    subplots_adjust(top=1.5)
    ax.plot(x, y3)
    xlabel("whs [pixels]")
    ax.set_ylabel("time [s]")
end
gcf()

# get the avg time in seconds for the section in timer
function get_avg_time(timer)
    return (TimerOutputs.time(timer)/TimerOutputs.ncalls(timer))/10e8
end

function get_timer(timer, sections)
    private_sections = copy(sections)
    try
        timer[private_sections[1]]
    catch e
        println(e)
        return TimerOutput()
    end
    out = timer
    while private_sections != []
        out = out[popfirst!(private_sections)]
    end
    return out
end

full_flow_estim, source_reg, flow_estim_at_inds, inds = sparse_lap_filt_one_by_one(img, imgw, 12)

inds

using ImageFiltering: KernelFactors.gaussian
using LAP_julia: inds_to_points, filt_onebyone!, window_sum3_at_inds!, multi_mat_div_qr
function sparse_lap_filt_one_by_one(target,
                                    source,
                                    fhs,
                                    window_size=(2fhs+1, 2fhs+1);
                                    spacing::Integer=15,
                                    point_count::Integer=100,
                                    timer::TimerOutput=TimerOutput("sparse_lap"),
                                    display::Bool=false,
                                    flow_interpolation_method::Symbol=:quad,
                                    base_method_kwargs=Dict(:timer => timer, :display => display))

    mask = parent(padarray(trues(size(target).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
    @timeit_debug timer "find edge points" begin
        inds = find_edge_points(target, spacing=spacing, number=point_count, mask=mask)
    end
    @timeit_debug timer "sparse lap" begin
        flow_estim_at_inds, inds = single_lap_at_points_filt_one_by_one(target, source, fhs, window_size, inds; base_method_kwargs...)
    end
    if all(isnan, flow_estim_at_inds)
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = zeros(size(target)) .* im .+ zeros(size(target))
        end
    else
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = interpolate_flow(flow_estim_at_inds, inds, size(target), method=flow_interpolation_method)
        end
    end
    @timeit_debug timer "generate source_reg" begin
        source_reg = warp_img(source, real(full_flow_estim), imag(full_flow_estim))
    end
    return full_flow_estim, source_reg, flow_estim_at_inds, inds
end

function single_lap_at_points_filt_one_by_one(target::Image,
                                              source::Image,
                                              filter_half_size::Integer,
                                              window,
                                              inds::Array{CartesianIndex,1};
                                              filter_count::Integer=3,
                                              timer::TimerOutput=TimerOutput("Sparse LAP"),
                                              display::Bool=false)

    image_size = size(target)
    pixel_count = length(target)
    ind_count = length(inds)
    filter_size = 2 * filter_half_size + 1

    # Prepare filters:
    # Calculate separable filters from basis:
    sigma = (filter_half_size + 2) / 4
    centered_inds = centered(-filter_half_size:filter_half_size)
    gaus = gaussian(sigma, filter_size)
    gaus_der = gaus .* centered_inds .* (-1)/sigma^2
    gaus_der_flip = reverse(gaus_der, dims=1)

    # prepare 2D filter basis for later use
    basis = similar(target, (filter_size, filter_size, filter_count))

    # dif = (points of forward image - points of backward image) for each filter
    dif = similar(target, (image_size..., filter_count))
    points = inds_to_points(inds)

    @timeit_debug timer "filtering" begin
        if filter_count == 3
            # temporary place to store filtered images
            # tmp_filtered_1 = fill(NaN, image_size)
            # tmp_filtered_2 = fill(NaN, image_size)
            tmp_filtered_1 = similar(target)
            tmp_filtered_2 = similar(target)

            # basis 1 - (gaus, gaus) both forward and backward
            kernf_1 = kernelfactors((gaus, gaus))
            basis[:, :, 1] = broadcast(*, kernf_1...)
            filt_onebyone!(tmp_filtered_1, target, kernf_1, filter_half_size, points)
            filt_onebyone!(tmp_filtered_2, source, kernf_1, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_1, "symmetric")
            # imfilter!(tmp_filtered_2, source, kernf_1, "symmetric")
            dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

            # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
            kernf_2f = kernelfactors((gaus, gaus_der_flip))
            basis[:, :, 2] = broadcast(*, kernf_2f...)
            filt_onebyone!(tmp_filtered_1, target, kernf_2f, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_2f, "symmetric")
            kernf_2b = kernelfactors((gaus, gaus_der))
            filt_onebyone!(tmp_filtered_2, source, kernf_2b, filter_half_size, points)
            # imfilter!(tmp_filtered_2, source, kernf_2b, "symmetric")
            dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

            # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
            kernf_3f = kernelfactors((gaus_der_flip, gaus))
            basis[:, :, 3] = broadcast(*, kernf_3f...)
            filt_onebyone!(tmp_filtered_1, target, kernf_3f, filter_half_size, points)
            # imfilter!(tmp_filtered_1, target, kernf_3f, "symmetric")
            kernf_3b = kernelfactors((gaus_der, gaus))
            filt_onebyone!(tmp_filtered_2, source, kernf_3b, filter_half_size, points)
            # imfilter!(tmp_filtered_2, source, kernf_3b, "symmetric")
            dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
        end
    end # "filtering"

    dif = reshape(dif, (:, filter_count))

    # prepare matrices for linear system of equations
    A = similar(target, (filter_count-1, filter_count-1, ind_count))
    b = similar(target, (filter_count-1, ind_count))

    @timeit_debug timer "prepare A and b" begin
        # TODO: check whether window_sum_around_points calculates only around points and not on whole window.
        for k in 1:filter_count-1
            for l in k:filter_count-1
                @timeit_debug timer "window sum part 1" begin
                    # @views window_sum_around_points!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, points)
                    @views window_sum3_at_inds!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, inds)
                    A[l, k, :] = A[k, l, :]
                end
            end
            @timeit_debug timer "window sum part 2" begin
                # @views window_sum_around_points!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, points)
                @views window_sum3_at_inds!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, inds)
            end
        end
    end # "prepare A and b"

    # Perform Gauss elimination on all pixels in parallel:
    # coeffs will be of shape: pixel_count, filter_count-1
    @timeit_debug timer "multi mat div" begin
        @views coeffs = multi_mat_div_qr(A, b)
    end

    # adding ones so that all base filters have their coefficients even the first one
    all_coeffs = [ones(ind_count) coeffs]

    k = (-filter_half_size:filter_half_size)

    # Get the displacement vector field from the filters
    @timeit_debug timer "calculate flow" begin
        u_est_at_inds = similar(target, ComplexF64, ind_count)
        @views u_est_at_inds[:] = 2 .* ((im .* (.-1 .* all_coeffs[:, 3]) ./ all_coeffs[:, 1]) .+ ((.-1 .* all_coeffs[:, 2]) ./ all_coeffs[:, 1]));
    end # "calculate flow"

    # dont use estimations whose displacement is larger than the filter_half_size
    displacement_mask = (real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .<= filter_half_size^2

    return u_est_at_inds[displacement_mask], inds[displacement_mask]

end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

"""
Return the Lar.Struct of octahedron
"""
function octahedron(l)
    a = l * (sqrt(3) / 2)
    h = sqrt(a^2 - ((l / 2)^2))
    #h=0.70710678118654
    V = [0 l 0 l l/2; 0 0 l l l/2; 0 0 0 0 h]
    VV = [[1], [2], [3], [4], [5]]
    EV = [[1, 2], [2, 3], [3, 4], [4, 1], [1, 5], [2, 5], [3, 5], [4, 5]]
    FV = [[1, 2, 5], [1, 3, 5], [3, 4, 5], [2, 4, 5]]
    CV = [[1, 2, 3, 4, 5]]
    tetrahedron = V, EV, FV, CV
    return Lar.Struct([tetrahedron, Lar.t(0, l, 0), Lar.r(π, 0, 0), tetrahedron])
end
"""
Return the Lar.Struct of an irregular octahedron whit h as height
"""
function octahedron(l, h)
    V = [0 l 0 l l / 2; 0 0 l l l / 2; 0 0 0 0 h]
    VV = [[1], [2], [3], [4], [5]]
    EV = [[1, 2], [3, 4], [4, 5], [1, 5], [1, 5], [2, 5], [3, 5], [4, 5]]
    FV = [[1, 2, 5], [1, 3, 5], [3, 4, 5], [2, 4, 5]]
    CV = [[1, 2, 3, 4, 5]]
    tetrahedron = V, EV, FV, CV
    return Lar.Struct([tetrahedron, Lar.t(0, l, 0), Lar.r(π, 0, 0), tetrahedron])
end
"""
Return the Lar.Struct of a rectangular base octahedron whit h as height
"""
function rectangularbase_octahedron(ab, da, h)
    V = [0 ab 0 ab ab / 2; 0 0 da da da / 2; 0 0 0 0 h]
    V2 = [0 ab 0 ab ab / 2; 0 0 da da da / 2; 0 0 0 0 -h]
    VV = [[1], [2], [3], [4], [5]]
    EV = [[1, 2], [3, 4], [4, 5], [1, 5], [1, 5], [2, 5], [3, 5], [4, 5]]
    FV = [[1, 2, 5], [1, 3, 5], [3, 4, 5], [2, 4, 5]]
    CV = [[1, 2, 3, 4, 5]]
    tetrahedron = V, EV, FV, CV
    tetrahedron2 = V2, EV, FV, CV
    return Lar.Struct([tetrahedron, tetrahedron2])
end


text:
# Finding python. This is slightly complicated in order to support using PyCall
# from libjulia. We check if python symbols are present in the current process
# and if so do not use the deps.jl file, getting everything we need from the
# current process instead.

proc_handle = unsafe_load(cglobal(:jl_exe_handle, Ptr{Cvoid}))

struct Dl_info
    dli_fname::Ptr{UInt8}
    dli_fbase::Ptr{Cvoid}
    dli_sname::Ptr{UInt8}
    dli_saddr::Ptr{Cvoid}
end
EnumProcessModules(hProcess, lphModule, cb, lpcbNeeded) =
    ccall(:K32EnumProcessModules, stdcall, Bool,
        (Ptr{Cvoid}, Ptr{Ptr{Cvoid}}, UInt32, Ptr{UInt32}),
        hProcess, lphModule, cb, lpcbNeeded)

symbols_present = false
@static if Sys.iswindows()
    lpcbneeded = Ref{UInt32}()
    proc_handle = ccall(:GetCurrentProcess, stdcall, Ptr{Cvoid}, ())
    handles = Vector{Ptr{Cvoid}}(undef, 20)
    if EnumProcessModules(proc_handle, handles, sizeof(handles), lpcbneeded) == 0
        resize!(handles, div(lpcbneeded[],sizeof(Ptr{Cvoid})))
        EnumProcessModules(proc_handle, handles, sizeof(handles), lpcbneeded)
    end
    # Try to find python if it's in the current process
    for handle in handles
        sym = ccall(:GetProcAddress, stdcall, Ptr{Cvoid},
            (Ptr{Cvoid}, Ptr{UInt8}), handle, "Py_GetVersion")
        sym != C_NULL || continue
        global symbols_present = true
        global libpy_handle = handle
        break
    end
else
    global symbols_present = hassym(proc_handle, :Py_GetVersion)
end

if !symbols_present
    # Python not present. Use deps.jl
    const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    isfile(depfile) || error("PyCall not properly installed. Please run Pkg.build(\"PyCall\")")
    include(depfile) # generated by Pkg.build("PyCall")
    # Only to be used at top-level - pointer will be invalid after reload
    libpy_handle = Libdl.dlopen(libpython, Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
    # need SetPythonHome to avoid warning, #299
    Py_SetPythonHome(libpy_handle, pyversion_build, PYTHONHOME)
else
    @static if Sys.iswindows()
        pathbuf = Vector{UInt16}(undef, 1024)
        ret = ccall(:GetModuleFileNameW, stdcall, UInt32,
            (Ptr{Cvoid}, Ptr{UInt16}, UInt32),
            libpy_handle, pathbuf, length(pathbuf))
        @assert ret != 0
        pathlen = something(findfirst(iszero, pathbuf)) - 1
        libname = String(Base.transcode(UInt8, pathbuf[1:pathlen]))
        if (Libdl.dlopen_e(libname) != C_NULL)
            const libpython = libname
        else
            const libpython = nothing
        end
    else
        libpy_handle = proc_handle
        # Now determine the name of the python library that these symbols are from
        some_address_in_libpython = Libdl.dlsym(libpy_handle, :Py_GetVersion)
        some_address_in_main_exe = Libdl.dlsym(proc_handle, Sys.isapple() ? :_mh_execute_header : :main)
        dlinfo1 = Ref{Dl_info}()
        dlinfo2 = Ref{Dl_info}()
        ccall(:dladdr, Cint, (Ptr{Cvoid}, Ptr{Dl_info}), some_address_in_libpython,
            dlinfo1)
        ccall(:dladdr, Cint, (Ptr{Cvoid}, Ptr{Dl_info}), some_address_in_main_exe,
            dlinfo2)
        if dlinfo1[].dli_fbase == dlinfo2[].dli_fbase
            const libpython = nothing
        else
            const libpython = unsafe_string(dlinfo1[].dli_fname)
        end
    end
    # If we're not in charge, assume the user is installing necessary python
    # libraries rather than messing with their configuration
    const conda = false
    # Top-level code (`_current_python`) needs `pyprogramname`.  We
    # don't need its value but we need to assign something to it:
    const pyprogramname = ""
end

const pyversion = vparse(split(Py_GetVersion(libpy_handle))[1])

# PyUnicode_* may actually be a #define for another symbol, so
# we cache the correct dlsym
const PyUnicode_AsUTF8String =
    findsym(libpy_handle, :PyUnicode_AsUTF8String, :PyUnicodeUCS4_AsUTF8String, :PyUnicodeUCS2_AsUTF8String)
const PyUnicode_DecodeUTF8 =
    findsym(libpy_handle, :PyUnicode_DecodeUTF8, :PyUnicodeUCS4_DecodeUTF8, :PyUnicodeUCS2_DecodeUTF8)

# Python 2/3 compatibility: cache symbols for renamed functions
if hassym(libpy_handle, :PyString_FromStringAndSize)
    const PyString_FromStringAndSize = :PyString_FromStringAndSize
    const PyString_AsStringAndSize = :PyString_AsStringAndSize
    const PyString_Size = :PyString_Size
    const PyString_Type = :PyString_Type
else
    const PyString_FromStringAndSize = :PyBytes_FromStringAndSize
    const PyString_AsStringAndSize = :PyBytes_AsStringAndSize
    const PyString_Size = :PyBytes_Size
    const PyString_Type = :PyBytes_Type
end

# hashes changed from long to intptr_t in Python 3.2
const Py_hash_t = pyversion < v"3.2" ? Clong : Int

# whether to use unicode for strings by default, ala Python 3
const pyunicode_literals = pyversion >= v"3.0"

if libpython == nothing
    macro pysym(func)
        esc(func)
    end
    macro pyglobal(name)
        :(cglobal($(esc(name))))
    end
    macro pyglobalobj(name)
        :(cglobal($(esc(name)), PyObject_struct))
    end
    macro pyglobalobjptr(name)
        :(unsafe_load(cglobal($(esc(name)), Ptr{PyObject_struct})))
    end
else
    macro pysym(func)
        :(($(esc(func)), libpython))
    end
    macro pyglobal(name)
        :(cglobal(($(esc(name)), libpython)))
    end
    macro pyglobalobj(name)
        :(cglobal(($(esc(name)), libpython), PyObject_struct))
    end
    macro pyglobalobjptr(name)
        :(unsafe_load(cglobal(($(esc(name)), libpython), Ptr{PyObject_struct})))
    end
end


text:
using RobustLeastSquares
using BaseTestNext

Logging.configure(level=Logging.DEBUG) # Useful to check convergence

# Test solve()
@testset "RobustLeastSqures" begin

    @testset "Solvers" begin
        s1 = 20
        s2 = 10

        rng = srand(0)
        A = rand(rng,s1,s2)
        b = rand(rng,s1)
        weights = 1 + 0.2(rand(rng,s1)-0.5)

        sol_qr = RobustLeastSquares.solve(A,b,weights,:qr)
        sol_normal = RobustLeastSquares.solve(A,b,weights,:normal)
        sol_cg = RobustLeastSquares.solve(A,b,weights,:cg)

        @test norm(sol_qr - sol_normal) < 1e-10
        @test norm(sol_qr - sol_cg) < 1e-10
    end

    # Test the MEstimators
    @testset "MEstimators" begin
        width = 2.0
        r = [5.0]

        s1 = 20
        s2 = 10

        rng = srand(0)
        A = rand(rng,s1,s2)
        b = rand(rng,s1)
        weights = 1 + 0.2(rand(rng,s1)-0.5)

        @test RobustLeastSquares.estimator_rho(r,L2Estimator())        ≈ [12.5]
        @test RobustLeastSquares.estimator_psi(r,L2Estimator())        ≈ [5.0]
        @test RobustLeastSquares.estimator_weight(r,L2Estimator())     ≈ [1.0]
        @test RobustLeastSquares.estimator_sqrtweight(r,L2Estimator()) ≈ [1.0]

        @test RobustLeastSquares.estimator_rho(r,L1Estimator())        ≈ [5.0]
        @test RobustLeastSquares.estimator_psi(r,L1Estimator())        ≈ [1.0]
        @test RobustLeastSquares.estimator_weight(r,L1Estimator())     ≈ [0.20]
        @test RobustLeastSquares.estimator_sqrtweight(r,L1Estimator()) ≈ [0.4472135954999579]

        # problem here: width not implemented correctly for L1L2Estimator
        @test RobustLeastSquares.estimator_rho(r,L1L2Estimator(1.0))        ≈ [5.3484692283495345]
        @test RobustLeastSquares.estimator_psi(r,L1L2Estimator(1.0))        ≈ [1.3608276348795434]
        @test RobustLeastSquares.estimator_weight(r,L1L2Estimator(1.0))     ≈ [0.2721655269759087]
        @test RobustLeastSquares.estimator_sqrtweight(r,L1L2Estimator(1.0)) ≈ [0.5216948600244291]

        @test RobustLeastSquares.estimator_rho(r,HuberEstimator(width))        ≈ [8.0]
        @test RobustLeastSquares.estimator_psi(r,HuberEstimator(width))        ≈ [2.0]
        @test RobustLeastSquares.estimator_weight(r,HuberEstimator(width))     ≈ [0.4]
        @test RobustLeastSquares.estimator_sqrtweight(r,HuberEstimator(width)) ≈ [0.6324555320336759]

        @test RobustLeastSquares.estimator_rho(r,FairEstimator(width))        ≈ [4.988948126018528]
        @test RobustLeastSquares.estimator_psi(r,FairEstimator(width))        ≈ [1.4285714285714286]
        @test RobustLeastSquares.estimator_weight(r,FairEstimator(width))     ≈ [0.2857142857142857]
        @test RobustLeastSquares.estimator_sqrtweight(r,FairEstimator(width)) ≈ [0.5345224838248488]

        @test RobustLeastSquares.estimator_rho(r,CauchyEstimator(width))        ≈ [3.9620029377331667]
        @test RobustLeastSquares.estimator_psi(r,CauchyEstimator(width))        ≈ [0.6896551724137931]
        @test RobustLeastSquares.estimator_weight(r,CauchyEstimator(width))     ≈ [0.13793103448275862]
        @test RobustLeastSquares.estimator_sqrtweight(r,CauchyEstimator(width)) ≈ [0.3713906763541037]

        @test RobustLeastSquares.estimator_rho(r,GemanEstimator())        ≈ [0.4807692307692308]
        @test RobustLeastSquares.estimator_psi(r,GemanEstimator())        ≈ [0.0073964497041420114]
        @test RobustLeastSquares.estimator_weight(r,GemanEstimator())     ≈ [0.0014792899408284023]
        @test RobustLeastSquares.estimator_sqrtweight(r,GemanEstimator()) ≈ [0.038461538461538464]

        @test RobustLeastSquares.estimator_rho(r,WelschEstimator(width))        ≈ [1.9961390917275446]
        @test RobustLeastSquares.estimator_psi(r,WelschEstimator(width))        ≈ [0.009652270681138546]
        @test RobustLeastSquares.estimator_weight(r,WelschEstimator(width))     ≈ [0.0019304541362277093]
        @test RobustLeastSquares.estimator_sqrtweight(r,WelschEstimator(width)) ≈ [0.04393693362340742]

        @test RobustLeastSquares.estimator_rho(r,TukeyEstimator(width))        ≈ [0.6666666666666666]
        @test RobustLeastSquares.estimator_psi(r,TukeyEstimator(width))        ≈ [0.0]
        @test RobustLeastSquares.estimator_weight(r,TukeyEstimator(width))     ≈ [0.0]
        @test RobustLeastSquares.estimator_sqrtweight(r,TukeyEstimator(width)) ≈ [0.0]

        @test RobustLeastSquares.estimator_sqrtweight(b,MultiEstimator(CauchyEstimator(width)=>1:10))[1:10] ≈ RobustLeastSquares.estimator_sqrtweight(b[1:10],CauchyEstimator(width))
        @test RobustLeastSquares.estimator_sqrtweight(b,MultiEstimator(CauchyEstimator(width)=>1:10))[11:20] ≈ RobustLeastSquares.estimator_sqrtweight(b[11:20],L2Estimator(width))
    end

    @testset "Reweighted least squares" begin
        width = 2.0
        r = [5.0]

        s1 = 20
        s2 = 10

        rng = srand(0)
        A = rand(rng,s1,s2)
        b = rand(rng,s1)
        weights = 1 + 0.2(rand(rng,s1)-0.5)

        # Test the overall solver works...
        est = MultiEstimator(CauchyEstimator(width)=>1:10)
        sol2_qr, res, weights = reweighted_lsqr(A,b,est;method=:qr, n_iter=30, refit = false)
        sol2_normal, res, weights = reweighted_lsqr(A,b,est;method=:normal, n_iter=30, refit = false)
        sol2_cg, res, weights = reweighted_lsqr(A,b,est;method=:cg, n_iter=30, refit = false)

        @test norm(sol2_qr - sol2_normal) < 1e-7
        @test norm(sol2_qr - sol2_cg) < 1e-7

        sol3_qr, res, weights = reweighted_lsqr(A,b,est;method=:qr, n_iter=30, refit = true)
        sol3_normal, res, weights = reweighted_lsqr(A,b,est;method=:normal, n_iter=30, refit = true)
        sol3_cg, res, weights = reweighted_lsqr(A,b,est;method=:cg, n_iter=30, refit = true)

        @test norm(sol3_qr - sol3_normal) < 1e-7
        @test norm(sol3_qr - sol3_cg) < 1e-7
    end
end


text:
# plot gaussian probabilistic distribution
# using cumulative distribution function

module GaussProbDistCdf
    using Plots, Distributions
    pyplot()

    function main(is_test=false)
        norm_dist = Normal(209.7, 4.84) # mu, std_dev

        zs = Array(193:229)
        ys = [cdf(norm_dist, z) for z in zs]
        plot(zs, ys, color=:red, label="CDF")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/prob_model/gauss_prob_dist_cdf/gauss_prob_dist_cdf.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
struct Model
  alphaPrior::Vector{Float64} # concentration parameter for the symmetric Dirichlet prior on document topics
  betaPrior::Vector{Float64} # concentration parameter for the symmetric Dirichlet prior on words
  corpus::AbstractCorpus

  # initialize an untrained model
  Model(alphaPrior::Vector{Float64},
        betaPrior::Vector{Float64},
        corpus::AbstractCorpus) = begin
    K = length(alphaPrior)
    m = new(
      alphaPrior,
      betaPrior,
      corpus)
    return m
  end

  # initialize a trained model
  Model(trainedModel::Model,
        corpus::AbstractCorpus) = begin
    m = new(
      trainedModel.alphaPrior,
      trainedModel.betaPrior,
      corpus
    )
    return m
  end
end

struct State
  topics::Array{Float64,2}
  topicSums::Vector{Float64}
  docSums::Array{Float64,2}
  assignments::Array{Array{Int64,1},1}
  conditionals::Array{Array{Float64,2},1} # the p paramter for the word assignment (cat/multinom) variable
  frozen::Bool

  # randomly initialize the state
  State(model::Model,
        corpus::AbstractCorpus) = begin # length of the lexicon
    K = length(model.alphaPrior)
    s = new(
      zeros(Float64, K, length(corpus.lexicon)), # topics
      zeros(Float64, K), # topicSums
      zeros(Float64, K, length(corpus.docs)), #docSums
      fill(Array{Int64,1}(undef,0), length(corpus.docs)), # assignments
      fill(Array{Int64,2}(undef,0,K), length(corpus.docs)),
      false
    )
    initializeAssignments(model,s,corpus)
    return s
  end

  # initialize the state from a trained model
  State(topics::Array{Float64,2},
        topicSums::Vector{Float64},
        docSums::Array{Float64,2},
        assignments::Array{Array{Int64,1},1},
        conditionals::Array{Array{Float64,2},1},
        frozen::Bool) = begin # length of the lexicon
    s = new(
      topics,
      topicSums,
      docSums,
      assignmens,
      conditionals,
      frozen
    )
    return s
  end
end

function AllTopics(state::State)
  alltopics = []
  for i in 1:length(state.assignments)
    append!(alltopics,state.assignments[i])
  end
  return convert(Array{Int,1},alltopics)
end

function initializeAssignments(model::Model,state::State,corpus::AbstractCorpus)
  for dd in 1:length(corpus)
    @inbounds words = corpus.docs[dd].terms
    @inbounds state.assignments[dd] = zeros(length(words))
    @inbounds state.conditionals[dd] = zeros(length(words), length(model.alphaPrior))
    for ww in 1:length(words)
      @inbounds word = words[ww]
      @inbounds state.conditionals[dd][ww,:] = model.alphaPrior
      topic = sampleMultinomial(ww,dd,state)
      @inbounds state.assignments[dd][ww] = topic
      updateSufficientStatistics(word, topic, dd,
                                  model.corpus.weights[dd][ww],
                                  state)
    end
  end
  return
end

function sampleMultinomial(word_ind::Int64,
                           document::Int64,
                           state::State)
  cond = state.conditionals[document][word_ind,:]
  pSum = sum(cond)
  r = rand() * pSum
  K = length(cond)
  for k in 1:K
    if r < cond[k]
      return k
    else
      @inbounds r -= cond[k]
    end
  end
  return 0
end

function cond_word(word::Int,
                   word_ind::Int,
                   document::Int,
                   model::Model,
                   state::State)
  V = size(state.topics, 2)
  for ii in 1:length(model.alphaPrior)
    @inbounds state.conditionals[document][word_ind,ii] =
        (state.docSums[ii, document] + model.alphaPrior[ii]) *
        (state.topics[ii, word] + model.betaPrior[word]) /
        (state.topicSums[ii] + V * model.betaPrior[word])
  end
  return
end

function log_beta(x::Vector{Float64})
  # compute natural log of the multivariate beta function
  lb = sum(loggamma.(x))
  lb -= loggamma(sum(x))
end

function joint_log_p(model::Model,
                     state::State)
  #calculate the full joint log likelihood, this is usefull for testing
  log_pz = 0
  for k in 1:length(model.alphaPrior)
    @inbounds log_pz += (log_beta(state.topics[k,:] .+ model.betaPrior) -
                log_beta(model.betaPrior))
  end
  for d in 1:length(model.corpus)
    @inbounds log_pz += (log_beta(state.docSums[:,d] .+ model.alphaPrior) -
                log_beta(model.alphaPrior))
  end
  return log_pz
end

function sampleWord(word::Int,
                    word_ind::Int,
                    document::Int,
                    model::Model,
                    state::State)
  cond_word(word, word_ind, document, model, state)
  sampleMultinomial(word_ind, document, state)
end

function updateSufficientStatistics(word::Int64,
                                    topic::Int64,
                                    document::Int64,
                                    scale::Float64,
                                    state::State)
  fr = Float64(!state.frozen)
  @inbounds state.docSums[topic, document] += scale
  @inbounds state.topicSums[topic] += scale * fr
  @inbounds state.topics[topic, word] += scale * fr
  return
end

@doc raw"""
    getTermDist(state::State, model::Model)

Compute ``\phi_{k,v} = \frac{\Psi_{k,v} + \beta_t}{\left( \sum^V_{v'=1} \Psi_{k,v'} + \beta_{v'}\right)}``

Where ``\vec{ϕ_v}`` parameterizes the V-dimensional categorical distribution of a word.

Updates the `termDist` attribute of `state`
"""
function getTermDist(state::State, model::Model)
  Phi = Array{Float64,2}(undef,length(model.alphaPrior),length(model.betaPrior))
  for topic in 1:length(model.alphaPrior)
    Phi[topic,:] = (state.topics[topic,:] .+ model.betaPrior) ./ (state.topicSums[topic] + sum(model.betaPrior))
  end
  return Phi
end

@doc raw"""
    getTopicDist(state::State, model::Model)

Compute ``\theta_{k,m} = \frac{\Omega_{k,m} + \alpha_k}{\left( \sum^K_{k'=1} \Omega_{k,m'} + \alpha_{k'}\right)}``

Where ``\vec{\theta_m}`` parameterizes the K-dimensional categorical distribution of a document.

Updates the `topicDist` attribute of `state`
"""
function getTopicDist(state::State, model::Model)
  Theta = Array{Float64,2}(undef,length(model.alphaPrior),length(model.corpus))
  for doc in 1:length(model.corpus)
    Theta[:,doc] = (state.docSums[:,doc] .+ model.alphaPrior) ./ (sum(state.docSums[:,doc]) + sum(model.alphaPrior))
  end
  return Theta
end

function sampleDocument(document::Int,
                        model::Model,
                        state::State)
  words = model.corpus.docs[document].terms
  Nw = length(words)
  @inbounds weights = model.corpus.weights[document]
  K = length(model.alphaPrior)
  @inbounds assignments = state.assignments[document]
  for ii in 1:Nw
    word = words[ii]
    oldTopic = assignments[ii]
    updateSufficientStatistics(word, oldTopic, document, -weights[ii], state)
    newTopic = sampleWord(word, ii, document, model, state)
    @inbounds assignments[ii] = newTopic
    updateSufficientStatistics(word, newTopic, document, weights[ii], state)
  end
  return
end

function sampleCorpus(model::Model, state::State)
  for ii in 1:length(model.corpus)
    sampleDocument(ii, model, state)
  end
  return
end

# The functions below are designed for public consumption
function trainModel(model::Model,
                    state::State,
                    numIterations::Int64)
  for ii in 1:numIterations
    println(string("Iteration ", ii, "..."))
    sampleCorpus(model, state)
  end
  return
end

function topTopicWords(model::Model,
                       state::State,
                       numWords::Int64)
  [model.corpus.lexicon[reverse(sortperm(state.topics'[1:end, row]))[1:numWords]]
   for row in 1:size(state.topics,1)]
end


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
tid = length(ARGS)==1 ? parse(Int64, ARGS[1]) : 3

np = pyimport("numpy")

"""
nx = 3 blocks, node 0 is at the origin
2     4     6     8
|  _  |  _  |  _  |
|  X  |  X  |  X  |
|  _  |  _  |  _  |
1     3     5     7

element number: left, up, bottom, diag, anti-diag, right
2i    2i+2
|  -  |
|  X  |
|  -  |
2i-1  2i+1

local node id: left to right, bottom to top
"""

ndofs = 2
nx = 3
nnodes, neles = 2*(nx + 1), 5*nx + 1

# PARAMETER: distance between nodes 0-1 and 0-2
l = 1.0
nodes = zeros(Float64, nnodes, ndofs)
for i = 1:nx+1
    nodes[2 * i - 1, 1], nodes[2 * i - 1, 2] = (i - 1)*l, 0.0
    nodes[2 * i, 1], nodes[2 * i, 2] = (i - 1)*l, l
end

EBC, g = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)
#fix node 1 and node 2
EBC[[1,2], :] .= -1

function ggt(t)
    return zeros(sum(EBC.==-2)), zeros(sum(EBC.==-2))
end
gt = ggt

#pull in the x direction
FBC, fext = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)

# todo PARAMETER
FORCE_TYPE = "non-constant"

if FORCE_TYPE == "constant"
    #pull in the y direction
    FBC[2*(nx+1), 1] = -1
    fext[2*(nx+1), 1] = 1.0e2
    FBC[2*nx+1, 1] = -1
    fext[2*nx+1, 1] = 1.0e2
else
    FBC[2*(nx+1), 2] = -2
end

#force load function
function fft(t)
    # return 20.0 * sin(2*pi*t)
    return (tid*0.2+0.6)
end
ft = fft

elements = []
ngp = 2
for i = 1:nx
    if i == 1
        elnodes = [2*i-1, 2*i]; coords = nodes[elnodes,:];
        push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
    end
    elnodes = [2*i+1, 2*i+2]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))

    elnodes = [2*i-1, 2*i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))

    elnodes = [2*i, 2*i+2]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))

    elnodes = [2*i-1, 2*i+2]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))

    elnodes = [2*i, 2*i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
end

# T = 0.5
# NT = 100
T = 0.2
NT = 200

Δt = T/NT
stress_scale = 100.0


text:
using RvSpectML
 using Statistics
 using Dates

make_plots = false
include("expres_1_read.jl")
order_list_timeseries = RvSpectML.make_order_list_timeseries(expres_data)
order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
lambda_range_with_data = (min = maximum(d->minimum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data), max = minimum(d->maximum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data) )
RvSpectML.normalize_spectra!(order_list_timeseries,expres_data);

espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks","G2.espresso.mas")
 espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
 line_list_df = espresso_df |>
     @filter(lambda_range_with_data.min <= _.lambda ) |>
    @filter( _.lambda < lambda_range_with_data.max) |>
#    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
#    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
    DataFrame

# Setup to run CCF
mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=1.6)
 line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
 ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=-5e3)
 v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)

# Compute CCF's & measure RVs
tstart = now()    # Compute CCFs for each order
 @time order_ccfs = RvSpectML.CCF.calc_order_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan) # line_list, mask_shape=mask_shape, plan=ccf_plan)
 println("# Order CCFs runtime: ", now()-tstart)

save_to_file = false
if save_to_file
   using CSV
   inst = EXPRES.EXPRES2D()
   for (i, order) in enumerate(EXPRES.orders_to_use_default(inst))
      t = Tables.table(order_ccfs[:,i,:]',header=Symbol.(v_grid))
      CSV.write("101501_ccf_order=" * string(order) * ".csv",t)
   end
end

make_plots = true
 if make_plots
  using Plots
end

if make_plots
 t_idx = 1
 order_idx = 1:33
 plot(v_grid,order_ccfs[:,order_idx,t_idx]./maximum(order_ccfs[:,order_idx,t_idx],dims=1),label=:none)
 xlabel!("v (m/s)")
 ylabel!("CCF")
 end

ccf_contrast = vec(minimum(mean(order_ccfs,dims=3),dims=1).-maximum(mean(order_ccfs,dims=3),dims=1))

rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,order_ccfs[:,j,i],fit_type = "gaussian")  for i in 1:length(order_list_timeseries), j in 1:size(order_ccfs,2) ]
  for i in 1:size(rvs_ccf_gauss,2)
    rvs_ccf_gauss[:,i] .-= mean(rvs_ccf_gauss[:,i])
  end

if make_plots
     using Plots
    nbin = 4
    plt_t = (order_list_timeseries.times .- minimum(order_list_timeseries.times) )
    times_binned = RvSpectML.bin_times_nightly(plt_t)
    plt = plot()
    scatter!(plt,plt_t,rvs_ccf_gauss[:,1],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,5],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,10],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,15],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,20],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,25],label=:none)
    scatter!(plt,plt_t,rvs_ccf_gauss[:,30],label=:none)
    xlabel!("Time (days)")
    ylabel!("v (m/s)")
 end

if make_plots
    order_idx = 1:33;
    plt = plot()
    for i in order_idx
       rms = std(rvs_ccf_gauss[:,i])
       rvs_binned = RvSpectML.bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_gauss[:,i])
       rms_binned = std(rvs_binned)
       println("  order = ", i, " RMS(order) = ", rms, "  RMS_binned = ", rms_binned)
       scatter!(plt,times_binned, rvs_binned,label=:none)
       plot!(plt,times_binned, rvs_binned,label=:none)
    end
    xlabel!("Time (hours)")
    ylabel!("v (m/s)")
    display(plt)
 end


text:
include("../notebook/PhaseLock.jl") #Defines PhaseLock.* & BodePlots.* utilities.
import InspectDR
using NumericIO
using Interact
using Blink

#Convenience aliases:
SI(x) = formatted(x, :SI, ndigits=3)
_throt(x) = throttle(0.1, x) #Default throttling for this example.

#When MIME"image/svg+xml" is disabled, Jupyter eventually requests PNG inline graphics:
#(Not certain if PNG is faster than SVG)
InspectDR.defaults.rendersvg = false

#Create base Observable-s that control plot output:
fmin = Observable(Float64(1e5)) #Minimum plot frequency
fmax = Observable(Float64(1e9)) #Maximum plot frequency
tmax = Observable(Float64(10e-6)) #Maximum plot time
ωnperiods = Observable(Float64(10)) #Maximum plot time (relative to ωn)

KdB = Observable(Float64(150)) #Gain
fp = Observable(Float64(500e3)) #Pole frequency
fz = Observable(Float64(50e6)) #Zero frequency

dispinline = Observable(Bool(true)) #Display inline?

#Compute open-loop characteristics
#Calculations are relatively expensive... throttle all inputs:
G = map(_throt(KdB), _throt(fp), _throt(fz)) do _KdB, _fp, _fz
	PhaseLock.Xf1(K=10.0^(_KdB/20), ωp=2pi*_fp, ωz=2pi*_fz)
end

function frequency_slider(sig::Observable; label::String="")
	rng = 10 .^ range(log10(1), stop=log10(1e9), length=100)
	return slider(rng, value=sig, label=label)
end

#Construct widgets used to control plots:
txt_fmin = spinbox(value=fmin, label="min freq")
txt_fmax = spinbox(value=fmax, label="max freq")
txt_tmax = spinbox(value=tmax, label="max time")
txt_ωnperiods = spinbox(value=ωnperiods, label="ωn periods")

sld_KdB = slider(-10:1.0:400, value=KdB, label="OL gain (dB)")
sld_fp = frequency_slider(fp, label="pole (Hz)")
sld_fz = frequency_slider(fz, label="zero (Hz)")

chk_inline = checkbox(dispinline, label="Display inline")
btn_creategtkplot = button("Create GTK plot")

info_poles = map(G) do _G
	txt(x::String) = "$x" #"\\textrm{$x}" #Not currently using LaTeX
	return HTML(txt("G (open-loop): ") *
		"K=" * SI(20*log10(_G.K)) * txt(" dB") *
		", 𝑓p=" * txt(SI(_G.ωp/(2pi)) * "Hz") *
		", 𝑓z=" * txt(SI(_G.ωz/(2pi)) * "Hz") *
		txt(" (") * "𝑓0=" * txt(SI(PhaseLock.ω_0(_G)/(2pi))) * "Hz)"
	)
end

info_sys = map(G) do _G
	s = PhaseLock.Sys2(_G)
	α = SI(s.α); ζ = SI(s.ζ); ωn = SI(s.ωn)
    return HTML("H (closed-loop): α=$α, ζ=$ζ, <b>ωn=$(ωn)rad/s</b>")
end

#Build widget lists to control plot:
widgetlistA = [txt_fmin, txt_fmax, txt_tmax, txt_ωnperiods, chk_inline, btn_creategtkplot]
widgetlistB = [sld_KdB, sld_fp, sld_fz, info_poles, info_sys]

#Widget-controlled plot:
pobj = PhaseLock.newplot() #Create initial plot object
#Calculations are relatively expensive... throttle all inputs:
updater_plot = map(_throt(fmin), _throt(fmax), _throt(tmax), _throt(ωnperiods), G) do _fmin, _fmax, _tmax, _ωnperiods,_G
	s = PhaseLock.Sys2(_G)
	return PhaseLock.update(pobj, _fmin, _fmax, _tmax, _ωnperiods, _G)
end

#Conditionally display inline plot with throttling to reduce traffic:
inline_plot = map(updater_plot, dispinline) do _plot, _inline
	if _inline
		return _plot
	else
		return nothing
	end
end

#Provide means to spawn new InspectDR-GTK plot windows:
creator_gtkplot = map(observe(btn_creategtkplot)) do _s
	#Display plot in Gtk GUI:
	gtkgui = display(InspectDR.GtkDisplay(), pobj) #Use same plot object

	#Refresh GUI when plot changes:
	updater_gtkgui = map(updater_plot) do _up
		InspectDR.refresh(gtkgui)
	end
end

#Display main window:
wnd = Window()
body!(wnd, vbox()) #first call to body! does not always work.
body!(wnd, vbox(hbox(vbox(widgetlistA...), vbox(widgetlistB...)), vbox(inline_plot)))

:DEMO_READY


text:
@enum(LintCodes,
MissingRef,
IncorrectCallArgs,
IncorrectIterSpec,
NothingEquality,
NothingNotEq,
ConstIfCondition,
EqInIfConditional,
PointlessOR,
PointlessAND,
UnusedBinding,
InvalidTypeDeclaration,
UnusedTypeParameter,
IncludeLoop,
MissingFile,
InvalidModuleName,
TypePiracy)

const LintCodeDescriptions = Dict{LintCodes,String}(IncorrectCallArgs => "Possible method call error.",
    IncorrectIterSpec => "A loop iterator has been used that will likely error.",
    NothingEquality => "Compare against `nothing` using `===`",
    NothingNotEq => "Compare against `nothing` using `!==`",
    ConstIfCondition => "A boolean literal has been used as the conditional of an if statement - it will either always or never run.",
    EqInIfConditional => "Unbracketed assignment in if conditional statements is not allowed, did you mean to use ==?",
    PointlessOR => "The first argument of a `||` call is a boolean literal.",
    PointlessAND => "The first argument of a `&&` call is `false`.",
    UnusedBinding => "The variable name has been bound but not used.",
    InvalidTypeDeclaration => "A non-DataType has been used in a type declaration statement.",
    UnusedTypeParameter => "A DataType parameter has been specified but not used.",
    IncludeLoop => "Loop detected, this file has already been included.",
    MissingFile => "The included file can not be found.",
    InvalidModuleName => "Module name matches that of its parent.",
    TypePiracy => "An imported function has been extended without using module defined typed arguments.")

haserror(m::Meta) = m.error !== nothing
haserror(x::EXPR) = hasmeta(x) && haserror(x.meta)
errorof(x::EXPR) = hasmeta(x) ? x.meta.error : nothing
function seterror!(x::EXPR, e)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.error = e
end

function _typeof(x, state)
    if typof(x) in (CSTParser.Abstract, CSTParser.Primitive, CSTParser.Struct, CSTParser.Mutable)
        return CoreTypes.DataType
    elseif typof(x) in (CSTParser.ModuleH, CSTParser.BareModule)
        return CoreTypes.Module
    elseif CSTParser.defines_function(x)
        return CoreTypes.Function
    end
end

# Call
function struct_nargs(x::EXPR)
    # struct defs wrapped in macros are likely to have some arbirtary additional constructors, so lets allow anything
    parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.MacroCall && return 0, typemax(Int), String[], true
    minargs, maxargs, kws, kwsplat = 0, 0, String[], false
    args = typof(x) === CSTParser.Mutable ? x.args[4] : x.args[3]
    args.args === nothing && return 0, typemax(Int), kws, kwsplat
    inner_constructor = findfirst(a->CSTParser.defines_function(a), args.args)
    if inner_constructor !== nothing
        return func_nargs(args.args[inner_constructor])
    else
        minargs = maxargs = length(args.args)
    end
    return minargs, maxargs, kws, kwsplat
end

function func_nargs(x::EXPR)
    minargs, maxargs, kws, kwsplat = 0, 0, String[], false
    sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
    if sig.args isa Vector{EXPR}
        for i = 2:length(sig.args)
            arg = sig.args[i]
            if typof(arg) === CSTParser.MacroCall && arg.args isa Vector{EXPR} && length(arg.args) > 1 &&
                _expr_assert(arg.args[1], MacroName, 2) && valof(arg.args[1].args[2]) == "nospecialize"
                if length(arg.args) == 2
                    arg = arg.args[2]
                end
            end
            if typof(arg) === CSTParser.PUNCTUATION
                # skip
            elseif typof(arg) === CSTParser.Parameters
                for j = 1:length(arg)
                    arg1 = arg.args[j]
                    if typof(arg1) === CSTParser.Kw
                        push!(kws, CSTParser.str_value(CSTParser.get_arg_name(arg1.args[1])))
                    elseif typof(arg1) === CSTParser.BinaryOpCall && kindof(arg1.args[2]) === CSTParser.Tokens.DDDOT
                        kwsplat = true
                    end
                end
            elseif typof(arg) === CSTParser.Kw
                if typof(arg.args[1]) === UnaryOpCall && kindof(arg.args[1].args[2]) === CSTParser.Tokens.DDDOT
                    maxargs = typemax(Int)
                else
                    maxargs !== typemax(Int) && (maxargs += 1)
                end
            elseif (typof(arg) === UnaryOpCall && kindof(arg.args[2]) === CSTParser.Tokens.DDDOT) ||
                (_binary_assert(arg, CSTParser.Tokens.DECLARATION) &&
                ((typof(arg.args[3]) === CSTParser.IDENTIFIER && valof(arg.args[3]) == "Vararg") ||
                (typof(arg.args[3]) === CSTParser.Curly && typof(arg.args[3].args[1]) === CSTParser.IDENTIFIER && valof(arg.args[3].args[1]) == "Vararg")))
                maxargs = typemax(Int)
            else
                minargs += 1
                maxargs !== typemax(Int) && (maxargs += 1)
            end
        end
    end
    return minargs, maxargs, kws, kwsplat
end

function func_nargs(m::SymbolServer.MethodStore)
    minargs, maxargs, kws, kwsplat = 0, 0, String[], false
    for a in m.args
        if last(a) == ".KW"
            if endswith(first(a), "...")
                kwsplat = true
            else
                push!(kws, first(a))
            end
        elseif startswith(last(a), "Vararg") || endswith(first(a), "...") || endswith(last(a), "...")
            maxargs = typemax(Int)
        else
            minargs += 1
            maxargs !== typemax(Int) && (maxargs += 1)
        end
    end
    return minargs, maxargs, kws, kwsplat
end

function call_nargs(x::EXPR)
    minargs, maxargs, kws = 0, 0, String[]
    if x.args isa Vector{EXPR}
        for i = 2:length(x.args)
            arg = x.args[i]
            if typof(arg) === CSTParser.PUNCTUATION
                # skip
            elseif typof(x.args[i]) === CSTParser.Parameters
                for j = 1:length(x.args[i])
                    arg = x.args[i].args[j]
                    if typof(arg) === CSTParser.Kw
                        push!(kws, CSTParser.str_value(CSTParser.get_arg_name(arg.args[1])))
                    end
                end
            elseif typof(arg) === CSTParser.Kw
                push!(kws, CSTParser.str_value(CSTParser.get_arg_name(arg.args[1])))
            elseif typof(arg) === CSTParser.UnaryOpCall && kindof(arg.args[2]) === CSTParser.Tokens.DDDOT
                maxargs = typemax(Int)
            else
                minargs += 1
                maxargs !== typemax(Int) && (maxargs += 1)
            end
        end
    else
        @info string("call_nargs: ", Expr(x))
    end

    return minargs, maxargs, kws
end

# compare_f_call(m_counts, call_counts) = true # fallback method

function compare_f_call(m_counts::Tuple{Int,Int,Array{String},Bool}, call_counts::Tuple{Int,Int,Array{String}})
    if call_counts[2] == typemax(Int)
        call_counts[1] <= call_counts[2] < m_counts[1] && return false
    else
        !(m_counts[1] <= call_counts[1] <= call_counts[2] <= m_counts[2]) && return false
    end
    if !m_counts[4] # no splatted kw in method sig
        length(call_counts[3]) > length(m_counts[3]) && return false # call has more kws than method accepts
        !all(kw in m_counts[3] for kw in call_counts[3]) && return false # call supplies a kw that isn't defined in the method
    else # splatted kw in method so accept any kw in call
    end
    return true
end

function check_call(x, server)
    if typof(x) === Call
        parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.Do && return # TODO: add number of args specified in do block.
        x.args === nothing || isempty(x.args) && return
        if typof(first(x.args)) === IDENTIFIER && hasref(first(x.args))
            func_ref = refof(first(x.args))
        elseif _binary_assert(first(x), CSTParser.Tokens.DOT) && typof(first(x).args[3]) === CSTParser.Quotenode && first(x).args[3].args !== nothing && !isempty(first(x).args[3].args) && typof(first(x).args[3].args[1]) === IDENTIFIER && hasref(first(x).args[3].args[1])
            func_ref = refof(first(last(first(x.args).args)))
        else
            return
        end

        if func_ref isa SymbolServer.FunctionStore || func_ref isa SymbolServer.DataTypeStore
            call_counts = call_nargs(x)
            for m in func_ref.methods
                m_counts = func_nargs(m)
                if compare_f_call(m_counts, call_counts)
                    return
                end
            end
            seterror!(x, IncorrectCallArgs)
        elseif func_ref isa Binding && (func_ref.type === CoreTypes.Function || func_ref.type === CoreTypes.DataType)
            call_counts = call_nargs(x)
            b = func_ref
            while b.next isa Binding && b.next.type == CoreTypes.Function
                b = b.next
            end
            while true
                if !(b isa Binding) # Needs to be cleaned up
                    if b isa SymbolServer.FunctionStore || b isa SymbolServer.DataTypeStore
                        for m in b.methods
                            m_counts = func_nargs(m)
                            if compare_f_call(m_counts, call_counts)
                                return
                            end
                        end
                        break
                    else
                        return
                    end
                elseif b.type == CoreTypes.Function
                    m_counts = func_nargs(b.val)
                elseif b.type == CoreTypes.DataType && CSTParser.defines_struct(b.val)
                    m_counts = struct_nargs(b.val)
                elseif b.val isa SymbolServer.FunctionStore || b.val isa SymbolServer.DataTypeStore
                    for m in b.val.methods
                        m_counts = func_nargs(m)
                        if compare_f_call(m_counts, call_counts)
                            return
                        end
                    end
                    break
                else
                    break
                end
                if compare_f_call(m_counts, call_counts)
                    return
                elseif CSTParser.rem_where_decl(CSTParser.get_sig(b.val)) == x
                    return
                end
                b.prev === nothing && break
                b = b.prev
            end
            seterror!(x, IncorrectCallArgs)
        end
    end
end

function check_loop_iter(x::EXPR, server)
    if typof(x) === CSTParser.For
        if length(x.args) > 1 && typof(x.args[2]) === CSTParser.BinaryOpCall && refof(x.args[2]) === nothing
            rng = x.args[2].args[3]
            if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                seterror!(x.args[2], IncorrectIterSpec)
            elseif typof(rng) === CSTParser.Call && refof(rng.args[1]) === getsymbolserver(server)["Base"].vals["length"]
                seterror!(x.args[2], IncorrectIterSpec)
            end
        end
    elseif typof(x) === CSTParser.Generator
        for i = 3:length(x.args)
            if typof(x.args[i]) === CSTParser.BinaryOpCall && refof(x.args[i]) === nothing
                rng = x.args[i].args[3]
                if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                    seterror!(x.args[i], IncorrectIterSpec)
                elseif typof(rng) === CSTParser.Call && valof(rng.args[1]) == "length" && refof(rng.args[1]) === getsymbolserver(server)["Base"].vals["length"]
                    seterror!(x.args[i], IncorrectIterSpec)
                end
            end
        end
    end
end

function check_nothing_equality(x::EXPR, server)
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x.args[2]) === CSTParser.Tokens.EQEQ && valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbolserver(server)["Core"].vals["nothing"]
            seterror!(x.args[2], NothingEquality)
        elseif kindof(x.args[2]) === CSTParser.Tokens.NOT_EQ && valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbolserver(server)["Core"].vals["nothing"]
            seterror!(x.args[2], NothingNotEq)
        end
    end
end

function _get_top_binding(x::EXPR, name::String)
    if scopeof(x) isa Scope
        return _get_top_binding(scopeof(x), name)
    elseif parentof(x) isa EXPR
        return _get_top_binding(parentof(x), name)
    else
        return nothing
    end
end
# Given the name of a function binding lookup the last declared method (i.e. the one
# stored in the scope)
_get_last_method(x, b::Binding, server) = b

function _get_last_method(x::EXPR, b::Binding, server)
    if scopeof(x) isa Scope
        if haskey(scopeof(x).names, b.name)
            if scopeof(x).names[b.name] isa Binding && scopeof(x).names[b.name].type == b.type
                return scopeof(x).names[b.name]
            end
        end
    elseif parentof(x) isa EXPR
        return _get_last_method(parentof(x), b.name, server)
    end
    return b
end

function _get_top_binding(s::Scope, name::String)
    if haskey(s.names, name)
        return s.names[name]
    elseif parentof(s) isa Scope
        return _get_top_binding(parentof(s), name)
    else
        return nothing
    end
end

function _get_global_scope(s::Scope)
    if !s.ismodule && parentof(s) isa Scope && parentof(s) != s
        return _get_global_scope(parentof(s))
    else
        return s
    end
end

function check_if_conds(x::EXPR)
    if typof(x) === CSTParser.If && length(x.args) > 2
        cond = typof(first(x.args)) == CSTParser.KEYWORD ? x.args[2] : x.args[1]
        if typof(cond) === CSTParser.LITERAL && (kindof(cond) === CSTParser.Tokens.TRUE || kindof(cond) === CSTParser.Tokens.FALSE)
            seterror!(cond, ConstIfCondition)
        elseif _binary_assert(cond, CSTParser.Tokens.EQ)
            seterror!(cond, EqInIfConditional)
        end
    end
end

function check_lazy(x::EXPR)
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x.args[2]) === CSTParser.Tokens.LAZY_OR
            if (typof(x.args[1]) === CSTParser.LITERAL && (kindof(x.args[1]) === CSTParser.Tokens.TRUE || kindof(x.args[1]) === CSTParser.Tokens.FALSE))
                seterror!(x, PointlessOR)
            end
        elseif kindof(x.args[2]) === CSTParser.Tokens.LAZY_AND
            if (typof(x.args[1]) === CSTParser.LITERAL && kindof(x.args[1]) === CSTParser.Tokens.FALSE) || (typof(x.args[3]) === CSTParser.LITERAL && kindof(x.args[3]) === CSTParser.Tokens.FALSE)
                seterror!(x, PointlessAND)
            end
        end
    end
end

function check_is_callable(x::EXPR, server)
    if typof(x) === CSTParser.Call
        func = first(x.args)

        if hasref(func)
        end
    end
end

function check_datatype_decl(x::EXPR, server)
    # Only call in function signatures
    if typof(x) === CSTParser.BinaryOpCall && length(x.args) === 3 && kindof(x.args[2]) === CSTParser.Tokens.DECLARATION &&
        parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.Call
        if hasref(last(x.args))
            dt = refof(last(x.args))
            if dt isa SymbolServer.DataTypeStore || dt isa Binding && dt.val isa SymbolServer.DataTypeStore
            elseif dt isa Binding && dt.type !== nothing
                safety_trip = 0
                while dt.type == CoreTypes.Function
                    safety_trip += 1
                    dt.prev === nothing && break
                    dt = dt.prev
                    if safety_trip > 50
                        @warn string(Expr(dt.name))
                        return
                    end
                end
                if dt.type !== nothing && dt.type !== CoreTypes.DataType
                    seterror!(x, InvalidTypeDeclaration)
                end
            end
        elseif typof(last(x.args)) === CSTParser.LITERAL
            seterror!(x, InvalidTypeDeclaration)
        end
    end
end

function check_modulename(x::EXPR)
    if (typof(x) === CSTParser.ModuleH || typof(x) === CSTParser.BareModule) && # x is a module
        scopeof(x) isa Scope && parentof(scopeof(x)) isa Scope && # it has a scope and a parent scope
        (typof(parentof(scopeof(x)).expr) === CSTParser.ModuleH ||
        typof(parentof(scopeof(x)).expr) === CSTParser.BareModule) && # the parent scope is a module
        valof(CSTParser.get_name(x)) == valof(CSTParser.get_name(parentof(scopeof(x)).expr)) # their names match
        seterror!(CSTParser.get_name(x), InvalidModuleName)
    end
end

mutable struct LintOptions
    call::Bool
    iter::Bool
    nothingcomp::Bool
    constif::Bool
    lazy::Bool
    datadecl::Bool
    typeparam::Bool
    modname::Bool
    pirates::Bool
end
LintOptions() = LintOptions(true, true, true, true, true, false, true, true, true)

function check_all(x::EXPR, opts::LintOptions, server)
    # Do checks
    opts.call && check_call(x, server)
    opts.iter && check_loop_iter(x, server)
    opts.nothingcomp && check_nothing_equality(x, server)
    opts.constif && check_if_conds(x)
    opts.lazy && check_lazy(x)
    # check_is_callable(x, server)
    opts.datadecl && check_datatype_decl(x, server)
    opts.typeparam && check_typeparams(x)
    opts.modname && check_modulename(x)
    opts.pirates && check_for_pirates(x)
    if x.args !== nothing
        for i in 1:length(x.args)
            check_all(x.args[i], opts, server)
        end
    end
end

function collect_hints(x::EXPR, missing = true, isquoted = false, errs = Tuple{Int,EXPR}[], pos = 0)
    if quoted(x)
        isquoted = true
    elseif isquoted && unquoted(x)
        isquoted = false
    end
    if typof(x) === CSTParser.ErrorToken
        # collect parse errors
        push!(errs, (pos, x))
    elseif !isquoted
        if missing && CSTParser.isidentifier(x) && !hasref(x) &&
            !(valof(x) == "var" && parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.NONSTDIDENTIFIER) &&
            !((valof(x) == "stdcall" || valof(x) == "cdecl" || valof(x) == "fastcall" || valof(x) == "thiscall" || valof(x) == "llvmcall") && is_in_fexpr(x, x->typof(x) === CSTParser.Call && isidentifier(x.args[1]) && valof(x.args[1]) == "ccall"))
            push!(errs, (pos, x))
        elseif haserror(x) && errorof(x) isa StaticLint.LintCodes
            # collect lint hints
            push!(errs, (pos, x))
        end
    end

    if x.args !== nothing
        for i in 1:length(x.args)
            collect_hints(x.args[i], missing, isquoted, errs, pos)
            pos += x.args[i].fullspan
        end
    end
    errs
end

function check_typeparams(x::EXPR)
    if typof(x) === CSTParser.WhereOpCall
        for i = 3:length(x.args)
            if hasbinding(x.args[i])
                if !(bindingof(x.args[i]).refs isa Vector)
                    seterror!(x.args[i], UnusedTypeParameter)
                elseif length(bindingof(x.args[i]).refs) == 1
                    # there should (will?) always be at least one reference in the where declaration
                    seterror!(x.args[i], UnusedTypeParameter)
                end
            end
        end
    end
end

function check_for_pirates(x::EXPR)
    if CSTParser.defines_function(x) && hasbinding(x) && overwrites_imported_function(bindingof(x))
        sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
        if typof(sig) == CSTParser.Call
            for i = 2:length(sig.args)
                if hasbinding(sig.args[i]) && bindingof(sig.args[i]).type isa Binding
                    return
                elseif refers_to_nonimported_type(sig.args[i])
                    return
                end
            end
            seterror!(x, TypePiracy)
        end
    end
end

function refers_to_nonimported_type(arg::EXPR)
    if hasref(arg) && refof(arg) isa Binding
        return true
    elseif typof(arg) === CSTParser.UnaryOpCall && length(arg.args) == 2 && kindof(arg.args[1]) === CSTParser.Tokens.DECLARATION
        return refers_to_nonimported_type(arg.args[2])
    elseif _binary_assert(arg, CSTParser.Tokens.DECLARATION)
        return refers_to_nonimported_type(arg.args[3])
    elseif typof(arg) === CSTParser.Curly
        for i = 1:length(arg.args)
            if refers_to_nonimported_type(arg.args[i])
                return true
            end
        end
        return false
    end
    return false
end

overwrites_imported_function(b, visited_bindings = Binding[]) = false
function overwrites_imported_function(b::Binding, visited_bindings = Binding[])
    if b in visited_bindings
        return false
    else
        push!(visited_bindings, b)
    end

    if b.prev isa Binding
        if b.prev.val isa EXPR
            # overwrites a within source bindig so lets check that
            return overwrites_imported_function(b.prev, visited_bindings)
        elseif (b.prev.val isa SymbolServer.FunctionStore || b.prev.val isa SymbolServer.DataTypeStore) && parentof(b.prev.name) isa EXPR && typof(parentof(b.prev.name)) === CSTParser.Import
            # explicitly imported, e.g. import ModuleName: somefunction
            return true
        end
    elseif b.prev isa SymbolServer.FunctionStore &&
        parentof(b.name) isa EXPR && typof(parentof(b.name)) === CSTParser.Quotenode &&
        parentof(parentof(b.name)) isa EXPR && typof(parentof(parentof(b.name))) === CSTParser.BinaryOpCall && kindof(parentof(parentof(b.name))[2]) === CSTParser.Tokens.DOT
        fullname = parentof(parentof(b.name))
        # overwrites imported function by declaring full name, e.g. ModuleName.FunctionName
        if CSTParser.isidentifier(fullname.args[1]) && refof(fullname.args[1]) isa SymbolServer.ModuleStore
            return true
        elseif typof(fullname.args[1]) === CSTParser.BinaryOpCall && kindof(fullname.args[1].args[2]) == CSTParser.Tokens.DOT && typof(fullname.args[1].args[3]) === CSTParser.Quotenode && refof(fullname.args[1].args[3].args[1]) isa SymbolServer.ModuleStore
            return true
        end
    end
    return false
end