text:
module OshiroDoor

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Oshiro Door" => Ahorn.EntityPlacement(
        Maple.OshiroDoor
    )
)

clutterDoorColor = (74, 71, 135, 153) ./ 255.0

function Ahorn.selection(entity::Maple.OshiroDoor)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 16, y - 16, 32, 32)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.OshiroDoor, room::Maple.Room)
    Ahorn.drawSprite(ctx, "characters/oshiro/oshiro24", 0, 0, alpha=0.8)
    Ahorn.drawRectangle(ctx, -16, -16, 32, 32, clutterDoorColor, (1.0, 1.0, 1.0, 8.0))
end

end


text:
# this is much faster if the sizeof(T) <=1
function sumby_smallT{T, S<:Number}(by::AbstractVector{T},  val::AbstractVector{S}; j = 1)::Dict{T,S}
  by_sim = similar(by)
  val1 = similar(val)
  hi = length(by)

  if hi == 0
      return Dict{T,S}();
  elseif hi == 1
      return Dict{T,S}(by[1] => val[1])
  end
  o = Forward
  lo = 1

  # Make sure we're sorting a bits type
  #TT = Base.Order.ordtype(o, by)
  # only use this in case of small T & number of uniques in by is smaller
  if sizeof(T) > 1
    error("only use this sort if the sizeof(eltype(by)) <= 1")
  elseif !isbits(T)
    error("Radix sort only sorts bits types (got $T)")
  end

  # Init
  bin1 = zeros(UInt32, 2^RADIX_SIZE) # just one bin for the initial run

  # Histogram for each element, radix
  for i = lo:hi
    v = uint_mapping(o, by[i])
    # perform counting then sort them
    idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
    @inbounds bin1[idx] += 1
  end
  # swaps  = 0
  #len = hi-lo+1
  v = uint_mapping(o, by[hi])
  idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

  # are all values the same at this radix?
  iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)

  if bin1[idx] == hi
    #println(by[1])
    #return Dict{T,S}(by[1] => sum(val))
    return sumby(by, val)
  end

  cbin = cumsum(bin1)
  non_empty_cnt = Int.(cbin[vcat(cbin[1] != 0, diff(cbin) .!= 0)])

  ci = cbin[idx]
  by_sim[ci] = by[hi]
  val1[ci] = val[hi]
  cbin[idx] -= 1

  # Finish the loop...
  @inbounds for i in hi-1:-1:lo
      v = uint_mapping(o, by[i])
      idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
      ci = cbin[idx]
      by_sim[ci] = by[i]
      val1[ci] = val[i]
      cbin[idx] -= 1
  end

  if iters == j
    #println(countmap(by_sim))
    return sumby_contiguous(by_sim, val1)
  else
    indices = vcat([1:non_empty_cnt[1]],
      [(a+1):b for (a,b) in zip(non_empty_cnt[1:end-1], non_empty_cnt[2:end])]
      )
    return mapreduce(merge, indices) do ii
      sumby_smallT(view(by_sim,ii), view(val1,ii); j = j + 1)
    end
  end
end

ext: jl


text:
# This is a test of the MRIReco package. In the future I plan to change the reconstruction step
# to use this. I will need to change the raw data output to a ISMRMRDFile.

using MRIReco

f = ISMRMRDFile("/home/ccp/Downloads/knee_3dFSE_slice170.h5")
acqData = AcquisitionData(f);

recoParams = Dict{Symbol, Any}()
recoParams[:reco] = "direct"
recoParams[:reconSize] = (320,320) # this size is also contained in acqData.encodingSize
img = reconstruction(acqData, recoParams)
# #samplingDensity.jl, v \ u uses too much memory

# ##
img = sqrt.(sum(img.^2,dims=5))
using Plots
heatmap(reverse( abs.(img[:,:,1,1,1]),dims=1), c=:viridis)

# multiCoil
# sampling
redFac= 4.0
acqDataSub = sample_kspace(acqData,redFac,"poisson",calsize=30,profiles=false);
# show sampling pattern
msk = zeros(acqDataSub.encodingSize[1],acqDataSub.encodingSize[2])
msk[acqDataSub.subsampleIndices[1]] .= 1
# #sampling
# heatmap(reverse( msk,dims=1), c=:grays)

#ESPIRIT
# coil sensitivities
smaps = espirit(acqData,(6,6),30,eigThresh_1=0.035,eigThresh_2=0.98)
# SENSE reconstruction
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (320,320) #previously it says :shape
params[:senseMaps] = smaps
params[:solver] = "cgnr"
params[:regularization] = "L2"
params[:λ] = 1e-4
params[:iterations] = 1
params[:normalizeReg] = true

img_reco = reconstruction(acqDataSub, params)

## show image
heatmap(reverse( abs.(img_reco[:,:,1,1,1]),dims=1), c=:viridis)

##
using MRIsim

imcoils = []
Nc = numChannels(data)
for coil = 1:Nc
    raw = reshape(data.kdata[1][:,coil], sz[1:2]...)
    img = MRIsim.ifftc(raw)
    append!(imcoils, [img])
end
imgss = sqrt.( sum([abs.(imcoils[i]).^2 for i=1:Nc]) ) #sum of squares

##

l = PlotlyJS.Layout(;height=500, width=500, colorbar=false)
p = PlotlyJS.heatmap(z = abs.(imgss), colorscale="Greys")
PlotlyJS.plot(p,l)


text:
module FiscomTools
include("test.jl")
include("incremental.jl")
include("bisection.jl")
include("ndiff.jl")
include("newton.jl")
include("gauss.jl")
include("intsimp.jl")
include("lucrout.jl")
include("trapz.jl")
include("diff_operators.jl")

export incremental, bisection_it, bisection_err, newtonR # For root approximation
export ndiff # Numerical differentiation
export redGauss, lgsolve, lucrout # Matrices and systems of equations
export trapz,intsimp # Integration techniques
export d1hat, d2hat # Differential operators of the first and second numerical derivative
end


text:
function reserve_slacks(optimization_container::OptimizationContainer, reserve_name::String)
    time_steps = model_time_steps(optimization_container)
    var_name = make_variable_name(reserve_name, SLACK_UP)
    variable = add_var_container!(optimization_container, var_name, time_steps)

    for jx in time_steps
        variable[jx] = JuMP.@variable(
            optimization_container.JuMPmodel,
            base_name = "$(var_name)_{$(jx)}",
            lower_bound = 0.0
        )
        JuMP.add_to_expression!(
            optimization_container.cost_function,
            variable[jx] * SERVICES_SLACK_COST,
        )
    end
    return variable
end


text:
"""
    struct CSF

Represents a single _configuration state function_ (CSF) adapted to angular momentum
symmetries. Depending on the type parameters, it can represent both non-relativistic CSFs in
[``LS``-coupling](@ref) and relativistic CSFs in [``jj``-coupling](@ref).

A CSF is defined by the following information:

* The configuration, i.e. an ordered list of subshells together with their occupations.
* A list of intermediate coupling terms (including the seniority quantum number to label
  states in degenerate subspaces) for each subshell in the configuration.
* A list of coupling terms for the coupling tree. The coupling is assumed to be done by,
  first, coupling two orbitals together, then coupling that to the next orbital, and so on.

An instance of a [`CSF`](@ref) object does not specify the ``J_z`` or ``L_z``/``S_z``
quantum number(s).

# Constructors

    CSF(configuration::Configuration, subshell_terms::Vector, terms::Vector)

Constructs an instance of a [`CSF`](@ref) from the information provided. The arguments have
different requirements depending on whether it is `configuration` is based on relativistic
or non-relativistic orbitals.

* If the configuration is based on [`Orbital`](@ref)s, `subshell_terms` must be a list of
  [`IntermediateTerm`](@ref)s and `terms` a list of [`Term`](@ref)s.
* If it is a configuration of [`RelativisticOrbital`](@ref)s, both `subshell_terms` and
  `terms` should both be a list of half-integer values.
"""
struct CSF{O<:AbstractOrbital, T<:Union{Term,HalfInteger}, S}
    config::Configuration{<:O}
    subshell_terms::Vector{IntermediateTerm{T,S}}
    terms::Vector{T}

    function CSF(config::Configuration{O},
                 subshell_terms::Vector{IntermediateTerm{T,S}},
                 terms::Vector{T}) where {O<:Union{<:Orbital,<:RelativisticOrbital},
                                          T<:Union{Term,HalfInt}, S}
        length(subshell_terms) == length(peel(config)) ||
            throw(ArgumentError("Need to provide $(length(peel(config))) subshell terms for $(config)"))
        length(terms) == length(peel(config)) ||
            throw(ArgumentError("Need to provide $(length(peel(config))) terms for $(config)"))

        for (i,(orb,occ,_)) in enumerate(peel(config))
            st = subshell_terms[i]
            assert_unique_classification(orb, occ, st) ||
                throw(ArgumentError("$(st) is not a unique classification of $(orb)$(to_superscript(occ))"))
        end

        new{O,T,S}(config, subshell_terms, terms)
    end

    CSF(config, subshell_terms::Vector{<:IntermediateTerm}, terms::Vector{<:Real}) =
        CSF(config, subshell_terms, convert.(HalfInt, terms))
end

"""
    const NonRelativisticCSF = CSF{<:Orbital,Term}
"""
const NonRelativisticCSF = CSF{<:Orbital,Term}

"""
    const RelativisticCSF = CSF{<:RelativisticOrbital,HalfInt}
"""
const RelativisticCSF = CSF{<:RelativisticOrbital,HalfInt}

Base.:(==)(a::CSF{O,T}, b::CSF{O,T}) where {O,T} =
    (a.config == b.config) && (a.subshell_terms == b.subshell_terms) && (a.terms == b.terms)

# We possibly want to sort on configuration as well
Base.isless(a::CSF, b::CSF) = last(a.terms) < last(b.terms)

num_electrons(csf::CSF) = num_electrons(csf.config)

"""
    orbitals(csf::CSF{O}) -> Vector

Access the underlying list of orbitals

```jldoctest
julia> csf = first(csfs(rc"1s2 2p-2"))
1s²   2p-²
     0     0
      0     0+

julia> orbitals(csf)
2-element Vector{RelativisticOrbital{Int64}}:
 1s
 2p-
```
"""
orbitals(csf::CSF) = orbitals(csf.config)

"""
    csfs(::Configuration) -> Vector{CSF}
    csfs(::Vector{Configuration}) -> Vector{CSF}

Generate all [`CSF`](@ref)s corresponding to a particular configuration or a set of
configurations.
"""
function csfs end

function csfs(config::Configuration)
    map(allchoices(intermediate_terms(peel(config)))) do subshell_terms
        map(intermediate_couplings(subshell_terms)) do coupled_terms
            CSF(config, subshell_terms, coupled_terms[2:end])
        end
    end |> c -> vcat(c...) |> sort
end

csfs(configs::Vector{Configuration{O}}) where O = vcat(map(csfs, configs)...)

Base.length(csf::CSF) = length(peel(csf.config))
Base.getindex(csf::CSF, i::Integer) =
    (peel(csf.config)[i],csf.subshell_terms[i],csf.terms[i])

Base.iterate(csf::CSF, (el, i)=(length(csf)>0 ? csf[1] : nothing,1)) =
    i > length(csf) ? nothing : (el, (csf[i==length(csf) ? i : i+1],i+1))

function Base.show(io::IO, csf::CSF)
    c = core(csf.config)
    p = peel(csf.config)
    nc = length(c)
    if nc > 0
        show(io, c)
        write(io, " ")
    end

    for (i,((orb,occ,state),st,t)) in enumerate(csf)
        show(io, orb)
        occ > 1 && write(io, to_superscript(occ))
        write(io, "($(st)|$(t))")
        i != lastindex(p) && write(io, " ")
    end
    print(io, iseven(parity(csf.config)) ? "+" : "-")
end

function Base.show(io::IO, ::MIME"text/plain", csf::RelativisticCSF)
    c = core(csf.config)
    nc = length(c)
    cw = length(string(c))

    w = 0
    p = peel(csf.config)
    for (orb,occ,state) in p
        w = max(w, length(string(orb))+length(digits(occ)))
    end

    w += 1

    cfmt = FormatExpr("{1:$(cw)s} ")
    ofmt = FormatExpr("{1:<$(w+1)s}")
    ifmt = FormatExpr("{1:>$(w+1)d}")
    rfmt = FormatExpr("{1:>$(w-1)d}/2")

    nc > 0 && printfmt(io, cfmt, c)
    for (orb,occ,state) in p
        printfmt(io, ofmt, join([string(orb), occ > 1 ? to_superscript(occ) : ""], ""))
    end
    println(io)

    nc > 0 && printfmt(io, cfmt, "")
    for it in csf.subshell_terms
        j = it.term
        if denominator(j) == 1
            printfmt(io, ifmt, numerator(j))
        else
            printfmt(io, rfmt, numerator(j))
        end
    end
    println(io)

    nc > 0 && printfmt(io, cfmt, "")
    print(io, " ")
    for j in csf.terms
        if denominator(j) == 1
            printfmt(io, ifmt, numerator(j))
        else
            printfmt(io, rfmt, numerator(j))
        end
    end
    print(io, iseven(parity(csf.config)) ? "+" : "-")
end

export CSF, NonRelativisticCSF, RelativisticCSF, csfs


text:
"""
```julia
produce_orbitdiagram(
	ds::ContinuousDynamicalSystem, plane, i::Int, p_index, pvalues; kwargs...
)
```
Produce an orbit diagram (sometimes wrongly called bifurcation diagram)
for the `i` variable(s) of the given continuous
system by computing Poincaré surfaces of section using `plane`
for the given parameter values (see [`poincaresos`](@ref)).

`i` can be `Int` or `AbstractVector{Int}`.
If `i` is `Int`, returns a vector of vectors. Else
it returns a vector of vectors of vectors.
Each entry are the points at each parameter value.

## Keyword Arguments
* `printparams::Bool = false` : Whether to print the parameter used
  during computation in order to keep track of running time.
* `direction, warning, Ttr, rootkw, diffeq` :
  Propagated into [`poincaresos`](@ref).
* `u0 = nothing` : Specify an initial state. If `nothing`, the previous state after each
  parameter is used to seed the new initial condition at the new parameter
  (with the very first state being the system's state). This makes convergence to the
  attractor faster, necessitating smaller `Ttr`. Otherwise `u0` can be a standard state,
  or a vector of states, so that a specific state is used for each parameter.

## Description
For each parameter, a PSOS reduces the system from a flow to a map. This then allows
the formal computation of an "orbit diagram" for the `i` variable
of the system, just like it is done in [`orbitdiagram`](@ref).

The parameter change is done as `p[p_index] = value` taking values from `pvalues`
and thus you must use a parameter container that supports this
(either `Array`, `LMArray`, dictionary or other).

See also [`poincaresos`](@ref), [`orbitdiagram`](@ref).
"""
function produce_orbitdiagram(
        ds::CDS{IIP, S, D}, plane, idxs, p_index, pvalues;
        tfinal::Real = 100.0, direction = -1, printparams = false, warning = true,
        Ttr = 0.0, u0 = nothing, rootkw = (xrtol = 1e-6, atol = 1e-6),
        diffeq = NamedTuple(), kwargs...
    ) where {IIP, S, D}

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    i = typeof(idxs) <: Int ? idxs : SVector{length(idxs), Int}(idxs...)

    _check_plane(plane, D)
    typeof(u0) <: Vector{<:AbstractVector} && @assert length(u0)==length(p)
    integ = integrator(ds; diffeq)
    planecrossing = PlaneCrossing(plane, direction > 0)
    p0 = ds.p[p_index]
    output = Vector{typeof(ds.u0[i])}[]
	plane_distance = (t) -> planecrossing(integ(t))

    for (n, p) in enumerate(pvalues)
        integ.p[p_index] = p
        printparams && println("parameter = $p")
        if typeof(u0) <: Vector{<:AbstractVector}
            st = u0[n]
        elseif isnothing(u0)
            st = integ.u
        else
            st = u0
        end
        reinit!(integ, st)
		Ttr != 0 && step!(integ, Ttr)
		data = _poincaresos(integ, plane_distance, planecrossing, tfinal+Ttr, i, rootkw)
		push!(output, data)
        warning && length(output[end]) == 0 && @warn "For parameter $p $PSOS_ERROR"
    end
    # Reset the parameter of the system:
    ds.p[p_index] = p0
    return output
end


text:
"""
    isscopedname(ex, modpath, name::Symbol)

Test whether expression `ex` is a module-scoped `name`. Both of these return `true`:

```
isscopedname(:(Base.OneTo), :Base, :OneTo)
isscopedname(:(Base.Checked.checked_add), (:Base, :Checked), :checked_add)
```
"""
function isscopedname(ex, modpath, name::Symbol)
  isexpr(ex, :(.), 2) &&
    (a = ex.args[2]; isa(a, QuoteNode) && a.value === name) &&
    hasscope(ex.args[1], modpath)
end
hasscope(modex, mod::Symbol) = modex === mod
hasscope(modex, mod::Tuple{Symbol}) = hasscope(modex, mod[1])
hasscope(modex, modpath::Tuple{Vararg{Symbol}}) =
  isscopedname(modex, Base.front(modpath), modpath[end])

"""
    isglobalref(g, mod, name)

Return true if `g` is equal to `GlobalRef(mod, name)`.
"""
isglobalref(g, mod, name) = isa(g, GlobalRef) && g.mod === mod && g.name === name


text:
@enum(
    LintCodes,

    MissingRef,
    IncorrectCallArgs,
    IncorrectIterSpec,
    NothingEquality,
    NothingNotEq,
    ConstIfCondition,
    EqInIfConditional,
    PointlessOR,
    PointlessAND,
    UnusedBinding,
    InvalidTypeDeclaration,
    UnusedTypeParameter,
    IncludeLoop,
    MissingFile,
    InvalidModuleName,
    TypePiracy,
    UnusedFunctionArgument,
    CannotDeclareConst,
    InvalidRedefofConst,
    NotEqDef,
    KwDefaultMismatch,
    InappropriateUseOfLiteral,
    ShouldBeInALoop,
    TypeDeclOnGlobalVariable,
    UnsupportedConstLocalVariable,
    UnassignedKeywordArgument,
    CannotDefineFuncAlreadyHasValue,
    DuplicateFuncArgName,
    IncludePathContainsNULL,
    IndexFromLength
)

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 a boolean literal.",
    UnusedBinding => "Variable has been assigned 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.",
    UnusedFunctionArgument => "An argument is included in a function signature but not used within its body.",
    CannotDeclareConst => "Cannot declare constant; it already has a value.",
    InvalidRedefofConst => "Invalid redefinition of constant.",
    NotEqDef => "`!=` is defined as `const != = !(==)` and should not be overloaded. Overload `==` instead.",
    KwDefaultMismatch => "The default value provided does not match the specified argument type.",
    InappropriateUseOfLiteral => "You really shouldn't be using a literal value here.",
    ShouldBeInALoop => "`break` or `continue` used outside loop.",
    TypeDeclOnGlobalVariable => "Type declarations on global variables are not yet supported.",
    UnsupportedConstLocalVariable => "Unsupported `const` declaration on local variable.",
    UnassignedKeywordArgument => "Keyword argument not assigned.",
    CannotDefineFuncAlreadyHasValue => "Cannot define function ; it already has a value.",
    DuplicateFuncArgName => "Function argument name not unique.",
    IncludePathContainsNULL => "Cannot include file, path contains NULL characters.",
    IndexFromLength => "Indexing with indices obtained from `length`, `size` etc is discouraged. Use `eachindex` or `axes` instead."
)

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

const default_options = (true, true, true, true, true, true, true, true, true, true)

struct LintOptions
    call::Bool
    iter::Bool
    nothingcomp::Bool
    constif::Bool
    lazy::Bool
    datadecl::Bool
    typeparam::Bool
    modname::Bool
    pirates::Bool
    useoffuncargs::Bool
end
LintOptions() = LintOptions(default_options...)
LintOptions(::Colon) = LintOptions(fill(true, length(default_options))...)

LintOptions(options::Vararg{Union{Bool,Nothing},length(default_options)}) =
    LintOptions(something.(options, default_options)...)

function check_all(x::EXPR, opts::LintOptions, env::ExternalEnv)
    # Do checks
    opts.call && check_call(x, env)
    opts.iter && check_loop_iter(x, env)
    opts.nothingcomp && check_nothing_equality(x, env)
    opts.constif && check_if_conds(x)
    opts.lazy && check_lazy(x)
    opts.datadecl && check_datatype_decl(x, env)
    opts.typeparam && check_typeparams(x)
    opts.modname && check_modulename(x)
    opts.pirates && check_for_pirates(x)
    opts.useoffuncargs && check_farg_unused(x)
    check_kw_default(x, env)
    check_use_of_literal(x)
    check_break_continue(x)
    check_const(x)

    if x.args !== nothing
        for i in 1:length(x.args)
            check_all(x.args[i], opts, env)
        end
    end
end

function _typeof(x, state)
    if x isa EXPR
        if headof(x) in (:abstract, :primitive, :struct)
            return CoreTypes.DataType
        elseif CSTParser.defines_module(x)
            return CoreTypes.Module
        elseif CSTParser.defines_function(x)
            return CoreTypes.Function
        end
    elseif x isa SymbolServer.DataTypeStore
        return CoreTypes.DataType
    elseif x isa SymbolServer.FunctionStore
        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 && CSTParser.ismacrocall(parentof(x)) && return 0, typemax(Int), Symbol[], true
    minargs, maxargs, kws, kwsplat = 0, 0, Symbol[], false
    args = x.args[3]
    length(args.args) == 0 && 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, Symbol[], false
    sig = CSTParser.rem_wheres_decls(CSTParser.get_sig(x))

    if sig.args !== nothing
        for i = 2:length(sig.args)
            arg = unwrap_nospecialize(sig.args[i])
            if isparameters(arg)
                for j = 1:length(arg.args)
                    arg1 = arg.args[j]
                    if iskwarg(arg1)
                        push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg1.args[1]))))
                    elseif isidentifier(arg1) || isdeclaration(arg1)
                        push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg1))))
                    elseif issplat(arg1)
                        kwsplat = true
                    end
                end
            elseif iskwarg(arg)
                if issplat(arg.args[1])
                    maxargs = typemax(Int)
                else
                    maxargs !== typemax(Int) && (maxargs += 1)
                end
            elseif issplat(arg) ||
                (isdeclaration(arg) &&
                ((isidentifier(arg.args[2]) && valofid(arg.args[2]) == "Vararg") ||
                (iscurly(arg.args[2]) && isidentifier(arg.args[2].args[1]) && valofid(arg.args[2].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, Symbol[], false

    for arg in m.sig
        if CoreTypes.isva(last(arg))
            maxargs = typemax(Int)
        else
            minargs += 1
            maxargs !== typemax(Int) && (maxargs += 1)
        end
    end
    for kw in m.kws
        if endswith(String(kw), "...")
            kwsplat = true
        else
            push!(kws, kw)
        end
    end
    return minargs, maxargs, kws, kwsplat
end

function call_nargs(x::EXPR)
    minargs, maxargs, kws = 0, 0, Symbol[]
    if length(x.args) > 0
        for i = 2:length(x.args)
            arg = x.args[i]
            if isparameters(arg)
                for j = 1:length(arg.args)
                    arg1 = arg.args[j]
                    if iskwarg(arg1)
                        push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg1.args[1]))))
                    end
                end
            elseif iskwarg(arg)
                push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg.args[1]))))
            elseif issplat(arg)
                maxargs = typemax(Int)
            else
                minargs += 1
                maxargs !== typemax(Int) && (maxargs += 1)
            end
        end
    else
        @info string("call_nargs: ", to_codeobject(x))
    end

    return minargs, maxargs, kws
end

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

function compare_f_call(
        (ref_minargs, ref_maxargs, ref_kws, kwsplat),
        (act_minargs, act_maxargs, act_kws),
    )
    # check matching on positional arguments
    if act_maxargs == typemax(Int)
        act_minargs <= act_maxargs < ref_minargs && return false
    else
        !(ref_minargs <= act_minargs <= act_maxargs <= ref_maxargs) && return false
    end

    # check matching on keyword arguments
    kwsplat && return true # splatted kw in method so accept any kw in call

    # no splatted kw in method sig
    length(act_kws) > length(ref_kws) && return false # call has more kws than method accepts
    !all(kw in ref_kws for kw in act_kws) && return false # call supplies a kw that isn't defined in the method

    return true
end

function is_something_with_methods(x::Binding)
    (CoreTypes.isfunction(x.type) && x.val isa EXPR) ||
    (CoreTypes.isdatatype(x.type) && x.val isa EXPR && CSTParser.defines_struct(x.val)) ||
    (x.val isa SymbolServer.FunctionStore || x.val isa SymbolServer.DataTypeStore)
end
is_something_with_methods(x::T) where T <: Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore} = true
is_something_with_methods(x) = false

function check_call(x, env::ExternalEnv)
    if iscall(x)
        parentof(x) isa EXPR && headof(parentof(x)) === :do && return # TODO: add number of args specified in do block.
        length(x.args) == 0 && return
        # find the function we're dealing with
        func_ref = refof_call_func(x)
        func_ref === nothing && return

        if is_something_with_methods(func_ref) && !(func_ref isa Binding && func_ref.val isa EXPR && func_ref.val.head === :macro)
            # intentionally empty
            if func_ref isa Binding && func_ref.val isa EXPR && isassignment(func_ref.val) && isidentifier(func_ref.val.args[1]) && isidentifier(func_ref.val.args[2])
                # if func_ref is a shadow binding (for these purposes, an assignment that just changes the name of a mehtod), redirect to the rhs of the assignment.
                func_ref = refof(func_ref.val.args[2])
            end
        else
            return
        end
        call_counts = call_nargs(x)
        tls = retrieve_toplevel_scope(x)
        tls === nothing && return @warn "Couldn't get top-level scope." # General check, this means something has gone wrong.
        func_ref === nothing && return
        !sig_match_any(func_ref, x, call_counts, tls, env) && seterror!(x, IncorrectCallArgs)
    end
end

function sig_match_any(func_ref::Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore}, x, call_counts, tls::Scope, env::ExternalEnv)
    iterate_over_ss_methods(func_ref, tls, env, m -> compare_f_call(func_nargs(m), call_counts))
end

function sig_match_any(func_ref::Binding, x, call_counts, tls::Scope, env::ExternalEnv)
    if func_ref.val isa SymbolServer.FunctionStore || func_ref.val isa SymbolServer.DataTypeStore
        match = sig_match_any(func_ref.val, x, call_counts, tls, env)
        match && return true
    end

    has_at_least_one_method = func_ref.val isa EXPR && defines_function(func_ref.val)
    # handle case where func_ref is typed as Function and yet has no methods

    for r in func_ref.refs
        method = get_method(r)
        method === nothing && continue
        has_at_least_one_method = true
        sig_match_any(method, x, call_counts, tls, env) && return true
    end
    return !has_at_least_one_method
end

function sig_match_any(func::EXPR, x, call_counts, tls::Scope, env::ExternalEnv)
    if CSTParser.defines_function(func)
        m_counts = func_nargs(func)
    elseif CSTParser.defines_struct(func)
        m_counts = struct_nargs(func)
    else
        return true # We shouldn't get here
    end
    if compare_f_call(m_counts, call_counts) || (CSTParser.rem_where_decl(CSTParser.get_sig(func)) == x)
        return true
    end
    return false
end

function get_method(name::EXPR)
    f = maybe_get_parent_fexpr(name, x -> CSTParser.defines_function(x) || CSTParser.defines_struct(x))
    if f !== nothing && CSTParser.get_name(f) == name
        return f
    end
end
function get_method(x::Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore})
    x
end
get_method(x) = nothing

isdocumented(x::EXPR) = parentof(x) isa EXPR && CSTParser.ismacrocall(parentof(x)) && headof(parentof(x).args[1]) === :globalrefdoc

function check_loop_iter(x::EXPR, env::ExternalEnv)
    if headof(x) === :for
        if length(x.args) > 1
            body = x.args[2]
            if headof(x.args[1]) === :block && x.args[1].args !== nothing
                for arg in x.args[1].args
                    check_incorrect_iter_spec(arg, body, env)
                end
            else
                check_incorrect_iter_spec(x.args[1], body, env)
            end
        end
    elseif headof(x) === :generator
        body = x.args[1]
        for i = 2:length(x.args)
            check_incorrect_iter_spec(x.args[i], body, env)
        end
    end
end

function check_incorrect_iter_spec(x, body, env)
    if x.args !== nothing && CSTParser.is_range(x)
        rng = rhs_of_iterator(x)

        if headof(rng) === :FLOAT || headof(rng) === :INTEGER || (iscall(rng) && refof(rng.args[1]) === getsymbols(env)[:Base][:length])
            seterror!(x, IncorrectIterSpec)
        elseif iscall(rng) && valof(rng.args[1]) == ":" &&
            length(rng.args) === 3 &&
            headof(rng.args[2]) === :INTEGER &&
            iscall(rng.args[3]) &&
            length(rng.args[3].args) > 1 && (
                refof(rng.args[3].args[1]) === getsymbols(env)[:Base][:length] ||
                refof(rng.args[3].args[1]) === getsymbols(env)[:Base][:size]
            )
            if length(x.args) >= 1
                lhs = x.args[1]
                arr = rng.args[3].args[2]
                b = refof(arr)

                # 1:length(arr) indexing is ok for Vector and Array specifically
                if b isa Binding && (CoreTypes.isarray(b.type) || CoreTypes.isvector(b.type))
                    return
                end
                if !all_underscore(valof(lhs))
                    if check_is_used_in_getindex(body, lhs, arr)
                        seterror!(x, IndexFromLength)
                    end
                end
            end
        end
    end
end

function check_is_used_in_getindex(expr, lhs, arr)
    if headof(expr) === :ref && expr.args !== nothing && length(expr.args) > 1
        this_arr = expr.args[1]
        if hasref(this_arr) && hasref(arr) && refof(this_arr) == refof(arr)
            for index_arg in expr.args[2:end]
                if hasref(index_arg) && hasref(lhs) && refof(index_arg) == refof(lhs)
                    seterror!(expr, IndexFromLength)
                    return true
                end
            end
        end
    end
    if expr.args !== nothing
        for arg in expr.args
            check_is_used_in_getindex(arg, lhs, arr) && return true
        end
    end
    return false
end

function check_nothing_equality(x::EXPR, env::ExternalEnv)
    if isbinarycall(x) && length(x.args) == 3
        if valof(x.args[1]) == "==" && (
                (valof(x.args[2]) == "nothing" && refof(x.args[2]) === getsymbols(env)[:Core][:nothing]) ||
                (valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbols(env)[:Core][:nothing])
            )
            seterror!(x.args[1], NothingEquality)
        elseif valof(x.args[1]) == "!=" && (
                (valof(x.args[2]) == "nothing" && refof(x.args[2]) === getsymbols(env)[:Core][:nothing]) ||
                (valof(x.args[3]) == "nothing" && refof(x.args[3]) === getsymbols(env)[:Core][:nothing])
            )
            seterror!(x.args[1], 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

function _get_top_binding(s::Scope, name::String)
    if scopehasbinding(s, 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 !CSTParser.defines_module(s.expr) && 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 headof(x) === :if
        cond = x.args[1]
        if headof(cond) === :TRUE || headof(cond) === :FALSE
            seterror!(cond, ConstIfCondition)
        elseif isassignment(cond)
            seterror!(cond, EqInIfConditional)
        end
    end
end

function check_lazy(x::EXPR)
    if isbinarysyntax(x)
        if valof(headof(x)) == "||"
            if headof(x.args[1]) === :TRUE || headof(x.args[1]) === :FALSE
                seterror!(x, PointlessOR)
            end
        elseif valof(headof(x)) == "&&"
            if headof(x.args[1]) === :TRUE || headof(x.args[1]) === :FALSE || headof(x.args[2]) === :TRUE || headof(x.args[2]) === :FALSE
                seterror!(x, PointlessAND)
            end
        end
    end
end

is_never_datatype(b, env::ExternalEnv) = false
is_never_datatype(b::SymbolServer.DataTypeStore, env::ExternalEnv) = false
function is_never_datatype(b::SymbolServer.FunctionStore, env::ExternalEnv)
    !(SymbolServer._lookup(b.extends, getsymbols(env)) isa SymbolServer.DataTypeStore)
end
function is_never_datatype(b::Binding, env::ExternalEnv)
    if b.val isa Binding
        return is_never_datatype(b.val, env)
    elseif b.val isa SymbolServer.FunctionStore
        return is_never_datatype(b.val, env)
    elseif CoreTypes.isdatatype(b.type)
        return false
    elseif b.type !== nothing
        return true
    end
    return false
end

function check_datatype_decl(x::EXPR, env::ExternalEnv)
    # Only call in function signatures?
    if isdeclaration(x) && parentof(x) isa EXPR && iscall(parentof(x))
        if (dt = refof_maybe_getfield(last(x.args))) !== nothing
            if is_never_datatype(dt, env)
                seterror!(x, InvalidTypeDeclaration)
            end
        elseif CSTParser.isliteral(last(x.args))
            seterror!(x, InvalidTypeDeclaration)
        end
    end
end

function check_modulename(x::EXPR)
    if CSTParser.defines_module(x) && # x is a module
        scopeof(x) isa Scope && parentof(scopeof(x)) isa Scope && # it has a scope and a parent scope
        CSTParser.defines_module(parentof(scopeof(x)).expr) && # 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

# Check whether function arguments are unused
function check_farg_unused(x::EXPR)
    if CSTParser.defines_function(x)
        sig = CSTParser.rem_wheres_decls(CSTParser.get_sig(x))
        if (headof(x) === :function && length(x.args) == 2 && x.args[2] isa EXPR && length(x.args[2].args) == 1 && CSTParser.isliteral(x.args[2].args[1])) ||
            (length(x.args) > 1 && headof(x.args[2]) === :block && length(x.args[2].args) == 1 && CSTParser.isliteral(x.args[2].args[1]))
            return # Allow functions that return constants
        end
        if iscall(sig)
            arg_names = Set{String}()
            for i = 2:length(sig.args)
                arg = sig.args[i]
                if arg.head === :parameters
                    for arg2 in arg.args
                        !check_farg_unused_(arg2, arg_names) && return
                    end
                else
                    !check_farg_unused_(arg, arg_names) && return
                end
            end
        end
    end
end

function check_farg_unused_(arg, arg_names)
    if hasbinding(arg)
    elseif iskwarg(arg) && hasbinding(arg.args[1])
        arg = arg.args[1]
    elseif is_nospecialize_call(arg) && hasbinding(unwrap_nospecialize(arg))
        arg = unwrap_nospecialize(arg)
    else
        return false
    end
    b = bindingof(arg)

    # We don't care about these
    valof(b.name) isa String && all_underscore(valof(b.name)) && return false

    if b === nothing ||
        # no refs:
       isempty(b.refs) ||
        # only self ref:
       (length(b.refs) == 1 && first(b.refs) == b.name) ||
        # first usage has binding:
        (length(b.refs) > 1 && b.refs[2] isa EXPR && hasbinding(b.refs[2]))
        seterror!(arg, UnusedFunctionArgument)
    end

    if valof(b.name) === nothing
    elseif valof(b.name) in arg_names
        seterror!(arg, DuplicateFuncArgName)
    else
        push!(arg_names, valof(b.name))
    end
    true
end

function unwrap_nospecialize(x)
    is_nospecialize_call(x) || return x
    x.args[3]
end

function is_nospecialize_call(x)
    CSTParser.ismacrocall(x) &&
    CSTParser.ismacroname(x.args[1]) &&
    is_nospecialize(x.args[1])
end

"""
collect_hints(x::EXPR, env, missingrefs = :all, isquoted = false, errs = Tuple{Int,EXPR}[], pos = 0)

Collect hints and errors from an expression. `missingrefs` = (:none, :id, :all) determines whether unresolved
identifiers are marked, the :all option will mark identifiers used in getfield calls."
"""
function collect_hints(x::EXPR, env, missingrefs=:all, isquoted=false, errs=Tuple{Int,EXPR}[], pos=0)
    if quoted(x)
        isquoted = true
    elseif isquoted && unquoted(x)
        isquoted = false
    end
    if headof(x) === :errortoken
        # collect parse errors
        push!(errs, (pos, x))
    elseif !isquoted
        if missingrefs != :none && isidentifier(x) && !hasref(x) &&
            !(valof(x) == "var" && parentof(x) isa EXPR && isnonstdid(parentof(x))) &&
            !((valof(x) == "stdcall" || valof(x) == "cdecl" || valof(x) == "fastcall" || valof(x) == "thiscall" || valof(x) == "llvmcall") && is_in_fexpr(x, x -> iscall(x) && 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
    elseif isquoted && missingrefs == :all && should_mark_missing_getfield_ref(x, env)
        push!(errs, (pos, x))
    end

    for i in 1:length(x)
        collect_hints(x[i], env, missingrefs, isquoted, errs, pos)
        pos += x[i].fullspan
    end

    errs
end

function refof_maybe_getfield(x::EXPR)
    if isidentifier(x)
        return refof(x)
    elseif is_getfield_w_quotenode(x)
        return refof(x.args[2].args[1])
    end
end

function should_mark_missing_getfield_ref(x, env)
    if isidentifier(x) && !hasref(x) && # x has no ref
    parentof(x) isa EXPR && headof(parentof(x)) === :quotenode && parentof(parentof(x)) isa EXPR && is_getfield(parentof(parentof(x)))  # x is the rhs of a getproperty
        lhsref = refof_maybe_getfield(parentof(parentof(x)).args[1])
        hasref(x) && return false # We've resolved
        if lhsref isa SymbolServer.ModuleStore || (lhsref isa Binding && lhsref.val isa SymbolServer.ModuleStore)
            # a module, we should know this.
            return true
        elseif lhsref isa Binding
            # by-use type inference runs after we've resolved references so we may not have known lhsref's type first time round, lets try and find `x` again
            resolve_getfield(x, lhsref, ResolveOnly(retrieve_scope(x), env, nothing)) # FIXME: Setting `server` to nothing might be sketchy?
            hasref(x) && return false # We've resolved
            if lhsref.val isa Binding
                lhsref = lhsref.val
            end
            lhsref = get_root_method(lhsref, nothing)
            if lhsref isa EXPR
                # Not clear what is happening here.
                return false
            elseif lhsref.type isa SymbolServer.DataTypeStore && !(isempty(lhsref.type.fieldnames) || isunionfaketype(lhsref.type.name) || has_getproperty_method(lhsref.type, env))
                return true
            elseif lhsref.type isa Binding && lhsref.type.val isa EXPR && CSTParser.defines_struct(lhsref.type.val) && !has_getproperty_method(lhsref.type)
                # We may have infered the lhs type after the semantic pass that was resolving references. Copied from `resolve_getfield(x::EXPR, parent_type::EXPR, state::State)::Bool`.
                if scopehasbinding(scopeof(lhsref.type.val), valof(x))
                    setref!(x, scopeof(lhsref.type.val).names[valof(x)])
                    return false
                end
                return true
            end
        end
    end
    return false
end

unwrap_fakeunionall(x) = x isa SymbolServer.FakeUnionAll ? unwrap_fakeunionall(x.body) : x
function has_getproperty_method(b::SymbolServer.DataTypeStore, env)
    getprop_vr = SymbolServer.VarRef(SymbolServer.VarRef(nothing, :Base), :getproperty)
    if haskey(getsymbolextendeds(env), getprop_vr)
        for ext in getsymbolextendeds(env)[getprop_vr]
            for m in SymbolServer._lookup(ext, getsymbols(env))[:getproperty].methods
                t = unwrap_fakeunionall(m.sig[1][2])
                !(t isa SymbolServer.FakeUnion) && t.name == b.name.name && return true
            end
        end
    else
        for m in getsymbols(env)[:Base][:getproperty].methods
            t = unwrap_fakeunionall(m.sig[1][2])
            !(t isa SymbolServer.FakeUnion) && t.name == b.name.name && return true
        end
    end
    return false
end

function has_getproperty_method(b::Binding)
    if b.val isa Binding || b.val isa SymbolServer.DataTypeStore
        return has_getproperty_method(b.val)
    elseif b isa Binding && CoreTypes.isdatatype(b.type)
        for ref in b.refs
            if ref isa EXPR && is_type_of_call_to_getproperty(ref)
                return true
            end
        end
    end
    return false
end

function is_type_of_call_to_getproperty(x::EXPR)
    function is_call_to_getproperty(x::EXPR)
        if iscall(x)
            func_name = x.args[1]
            return (isidentifier(func_name) && valof(func_name) == "getproperty") || # getproperty()
            (is_getfield_w_quotenode(func_name) && isidentifier(func_name.args[2].args[1]) && valof(func_name.args[2].args[1]) == "getproperty") # Base.getproperty()
        end
        return false
    end

    return parentof(x) isa EXPR && parentof(parentof(x)) isa EXPR &&
        ((isdeclaration(parentof(x)) && x === parentof(x).args[2] && is_call_to_getproperty(parentof(parentof(x)))) ||
        (iscurly(parentof(x)) && x === parentof(x).args[1] && isdeclaration(parentof(parentof(x))) &&  parentof(parentof(parentof(x))) isa EXPR && is_call_to_getproperty(parentof(parentof(parentof(x))))))
end

isunionfaketype(t::SymbolServer.FakeTypeName) = t.name.name === :Union && t.name.parent isa SymbolServer.VarRef && t.name.parent.name === :Core

function check_typeparams(x::EXPR)
    if iswhere(x)
        for i in 2:length(x.args)
            a = x.args[i]
            if hasbinding(a) && (bindingof(a).refs === nothing || length(bindingof(a).refs) < 2)
                seterror!(a, UnusedTypeParameter)
            end
        end
    end
end

function check_for_pirates(x::EXPR)
    if CSTParser.defines_function(x)
        sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
        fname = CSTParser.get_name(sig)
        if fname_is_noteq(fname)
            seterror!(x, NotEqDef)
        elseif iscall(sig) && hasbinding(x) && overwrites_imported_function(refof(fname))
            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 fname_is_noteq(x)
    if x isa EXPR
        if isoperator(x) && valof(x) == "!="
            return true
        elseif is_getfield_w_quotenode(x)
            return fname_is_noteq(x.args[2].args[1])
        end
    end
    return false
end

function refers_to_nonimported_type(arg::EXPR)
    arg = CSTParser.rem_wheres(arg)
    if hasref(arg) && refof(arg) isa Binding
        return true
    elseif isunarysyntax(arg) && (valof(headof(arg)) == "::" || valof(headof(arg)) == "<:")
        return refers_to_nonimported_type(arg.args[1])
    elseif isdeclaration(arg)
        return refers_to_nonimported_type(arg.args[2])
    elseif iscurly(arg)
        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) = false
function overwrites_imported_function(b::Binding)
    if ((b.val isa SymbolServer.FunctionStore || b.val isa SymbolServer.DataTypeStore) &&
        (is_in_fexpr(b.name, x -> headof(x) === :import)) || (b.refs isa Vector && length(b.refs) > 0 && (first(b.refs) isa SymbolServer.FunctionStore || first(b.refs) isa SymbolServer.DataTypeStore)))
        return true
    end
    return false
end

# Now called from add_binding
# Should return true/false indicating whether the binding should actually be added?
function check_const_decl(name::String, b::Binding, scope)
    # assumes `scopehasbinding(scope, name)`
    b.val isa Binding && return check_const_decl(name, b.val, scope)
    if b.val isa EXPR && (CSTParser.defines_datatype(b.val) || is_const(bind))
        seterror!(b.val, CannotDeclareConst)
    else
        prev = scope.names[name]
        if (CoreTypes.isdatatype(prev.type) && !is_mask_binding_of_datatype(prev)) || is_const(prev)
            if b.val isa EXPR && prev.val isa EXPR && !in_same_if_branch(b.val, prev.val)
                return
            end
            if b.val isa EXPR
                seterror!(b.val, InvalidRedefofConst)
            else
                # TODO check what's going on here
                seterror!(b.name, InvalidRedefofConst)
            end
        end
    end
end

function is_mask_binding_of_datatype(b::Binding)
    b.val isa EXPR && CSTParser.isassignment(b.val) && (rhsref = refof(b.val.args[2])) !== nothing && (rhsref isa SymbolServer.DataTypeStore || (rhsref.val isa EXPR && rhsref.val isa SymbolServer.DataTypeStore) || (rhsref.val isa EXPR && CSTParser.defines_datatype(rhsref.val)))
end

# check whether a and b are in all the same :if blocks and in the same branches
in_same_if_branch(a::EXPR, b::EXPR) = in_same_if_branch(find_if_parents(a), find_if_parents(b))
in_same_if_branch(a::Dict, b::EXPR) = in_same_if_branch(a, find_if_parents(b))
function in_same_if_branch(a::Dict, b::Dict)
    return length(a) == length(b) && all(k in keys(b) for k in keys(a)) && all(a[k] == b[k] for k in keys(a))
end

# find any parent nodes that are :if blocks and a pseudo-index of which branch
# x is in
function find_if_parents(x::EXPR, current=Int[], list=Dict{EXPR,Vector{Int}}())
    if x.head in (:block, :elseif) && parentof(x) isa EXPR && headof(parentof(x)) in (:if, :elseif)
        i = 1
        while i <= length(parentof(x).args)
            if parentof(x).args[i] == x
                pushfirst!(current, i)
                break
            end
            i += 1
        end
        if headof(parentof(x)) == :if
            list[parentof(x)] = current
            current = []
        end
    end
    return parentof(x) isa EXPR ? find_if_parents(parentof(x), current, list) : list
end

is_const(x) = false
is_const(b::Binding) = is_const(b.val)
is_const(x::EXPR) = is_const_expr(parentof(x))

is_const_expr(x) = false
is_const_expr(x::EXPR) = headof(x) === :const

"""
    check_kw_default(x::EXPR, server)

Check that the default value matches the type for keyword arguments. Following types are
checked: `String, Symbol, Int, Char, Bool, Float32, Float64, UInt8, UInt16, UInt32,
UInt64, UInt128`.
"""
function check_kw_default(x::EXPR, env::ExternalEnv)
    if headof(x) == :kw && isdeclaration(x.args[1]) && CSTParser.isliteral(x.args[2]) && hasref(x.args[1].args[2])
        decl_T = refof(x.args[1].args[2])
        rhs = x.args[2]
        rhsval = valof(rhs)
        if decl_T == getsymbols(env)[:Core][:String] && !CSTParser.isstringliteral(rhs)
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Symbol] && headof(rhs) !== :IDENTIFIER
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Int] && headof(rhs) !== :INTEGER
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][Sys.WORD_SIZE == 64 ? :Int64 : :Int32] && headof(rhs) !== :INTEGER
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Bool] && !(headof(rhs) === :TRUE || headof(rhs) === :FALSE)
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Char] && headof(rhs) !== :CHAR
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Float64] && headof(rhs) !== :FLOAT
            seterror!(rhs, KwDefaultMismatch)
        elseif decl_T == getsymbols(env)[:Core][:Float32] && !(headof(rhs) === :FLOAT && occursin("f", rhsval))
            seterror!(rhs, KwDefaultMismatch)
        else
            for T in (UInt8, UInt16, UInt32, UInt64, UInt128)
                if decl_T == getsymbols(env)[:Core][Symbol(T)]
                    # count the digits without prefix (=0x, 0o, 0b) and make sure it fits
                    # between upper and lower literal boundaries for `T` where the boundaries
                    # depend on the type of literal (binary, octal, hex)
                    n = count(x -> x != '_', rhsval) - 2
                    ub = sizeof(T)
                    lb = ub ÷ 2
                    if headof(rhs) == :BININT
                        8lb < n <= 8ub || seterror!(rhs, KwDefaultMismatch)
                    elseif headof(rhs) == :OCTINT
                        3lb < n <= 3ub || seterror!(rhs, KwDefaultMismatch)
                    elseif headof(rhs) == :HEXINT
                        2lb < n <= 2ub || seterror!(rhs, KwDefaultMismatch)
                    else
                        seterror!(rhs, KwDefaultMismatch)
                    end
                end
            end
            # signed integers of non native size can't be declared as literal
            for T in (Int8, Int16, Sys.WORD_SIZE == 64 ? Int32 : Int64, Int128)
                if decl_T == getsymbols(env)[:Core][Symbol(T)]
                    seterror!(rhs, KwDefaultMismatch)
                end
            end

        end
    end
end

function check_use_of_literal(x::EXPR)
    if CSTParser.defines_module(x) && length(x.args) > 1 && isbadliteral(x.args[2])
        seterror!(x.args[2], InappropriateUseOfLiteral)
    elseif (CSTParser.defines_abstract(x) || CSTParser.defines_primitive(x)) && isbadliteral(x.args[1])
        seterror!(x.args[1], InappropriateUseOfLiteral)
    elseif CSTParser.defines_struct(x) && isbadliteral(x.args[2])
        seterror!(x.args[2], InappropriateUseOfLiteral)
    elseif (isassignment(x) || iskwarg(x)) && isbadliteral(x.args[1])
        seterror!(x.args[1], InappropriateUseOfLiteral)
    elseif isdeclaration(x) && isbadliteral(x.args[2])
        seterror!(x.args[2], InappropriateUseOfLiteral)
    elseif isbinarycall(x, "isa") && isbadliteral(x.args[3])
        seterror!(x.args[3], InappropriateUseOfLiteral)
    end
end

isbadliteral(x::EXPR) = CSTParser.isliteral(x) && (CSTParser.isstringliteral(x) || headof(x) === :INTEGER || headof(x) === :FLOAT || headof(x) === :CHAR || headof(x) === :TRUE || headof(x) === :FALSE)

function check_break_continue(x::EXPR)
    if iskeyword(x) && (headof(x) === :CONTINUE || headof(x) === :BREAK) && !is_in_fexpr(x, x -> headof(x) in (:for, :while))
        seterror!(x, ShouldBeInALoop)
    end
end

function check_const(x::EXPR)
    if headof(x) === :const
        if VERSION < v"1.8.0-DEV.1500" && CSTParser.isassignment(x.args[1]) && CSTParser.isdeclaration(x.args[1].args[1])
            seterror!(x, TypeDeclOnGlobalVariable)
        elseif headof(x.args[1]) === :local
            seterror!(x, UnsupportedConstLocalVariable)
        end
    end
end

function check_unused_binding(b::Binding, scope::Scope)
    if headof(scope.expr) !== :struct && headof(scope.expr) !== :tuple && !all_underscore(valof(b.name))
        refs = loose_refs(b)
        if (isempty(refs) || length(refs) == 1 && refs[1] == b.name) &&
                !is_sig_arg(b.name) && !is_overwritten_in_loop(b.name) &&
                !is_overwritten_subsequently(b, scope) && !is_kw_of_macrocall(b)
            seterror!(b.name, UnusedBinding)
        end
    end
end

all_underscore(s) = false
all_underscore(s::String) = all(==(0x5f), codeunits(s))

function is_sig_arg(x)
    is_in_fexpr(x, CSTParser.iscall)
end

function is_kw_of_macrocall(b::Binding)
    b.val isa EXPR && isassignment(b.val) && parentof(b.val) isa EXPR && CSTParser.ismacrocall(parentof(b.val))
end

function is_overwritten_in_loop(x)
    # Cuts out false positives for check_unused_binding - the linear nature of our
    # semantic passes mean a variable declared at the end of a loop's block but used at
    # the start won't appear to be referenced.

    # Cheap version:
    # is_in_fexpr(x, x -> x.head === :while || x.head === :for)

    # We really want to check whether the enclosing scope(s) of the loop has a binding
    # with matching name.
    # Is this too expensive?
    loop = maybe_get_parent_fexpr(x, x -> x.head === :while || x.head === :for)
    if loop !== nothing
        s = scopeof(loop)
        if s isa Scope && parentof(s) isa Scope
            s2 = check_parent_scopes_for(s, valof(x))
            if s2 isa Scope
                prev_binding = parentof(s2).names[valof(x)]
                if prev_binding isa Binding
                    return true
                    # s = ComesBefore(prev_binding.name, s2.expr, 0)
                    # traverse(parentof(s2).expr, s)
                    # return s.result == 1
                    # for r in prev_binding.refs
                    #     if r isa EXPR && is_in_fexpr(r, x -> x === loop)
                    #         return true
                    #     end
                    # end
                else
                    return false
                end
            end
        else
            return false
        end
    else
        false
    end
    false
end

"""
    ComesBefore

Check whether x1 comes before x2
"""
mutable struct ComesBefore
    x1::EXPR
    x2::EXPR
    result::Int
end

function (state::ComesBefore)(x::EXPR)
    state.result > 0 && return
    if x == state.x1
        state.result = 1
        return
    elseif x == state.x2
        state.result = 2
        return
    end
    if !hasscope(x)
        traverse(x, state)
        state.result > 0 && return
    end
end

"""
    check_parent_scopes_for(s::Scope, name)

Checks whether the parent scope of `s` has the name `name`.
"""
function check_parent_scopes_for(s::Scope, name)
    # This returns `s` rather than the parent so that s.expr can be used in the linear
    # search (e.g. `bound_before`)
    if s.expr.head !== :module && parentof(s) isa Scope && haskey(parentof(s).names, name)
        s
    elseif s.parent isa Scope
        check_parent_scopes_for(parentof(s), name)
    end
end

function is_overwritten_subsequently(b::Binding, scope::Scope)
    valof(b.name) === nothing && return false
    s = BoundAfter(b.name, valof(b.name), 0)
    traverse(scope.expr, s)
    return s.result == 2
end

"""
    ComesBefore

Check whether x1 comes before x2
"""
mutable struct BoundAfter
    x1::EXPR
    name::String
    result::Int
end

function (state::BoundAfter)(x::EXPR)
    state.result > 1 && return
    if x == state.x1
        state.result = 1
        return
    end
    if scopeof(x) isa Scope && haskey(scopeof(x).names, state.name)
        state.result = 2
        return
    end
    traverse(x, state)
end


text:
module JSom

using DataFrames
using DataStructures
using Distances
using Hexagons
using StatsBase

import Base.size

export
    GridSOM,
    HexSOM,

    # SOM objects methods
    get_unit_weight,
    set_unit_weight,
    size,
    update,
    get_BMU,
    activate,
    sequential_random_epoch,
    sequential_epoch,
    quantize,
    quantization_error,
    bmu_map,
    activation_response,
    umatrix,
    reset_state,

    # Neighborhood functions
    Gaussian_Neighborhood,
    Ricker_Neighborhood,
    Triangular_Neighborhood

abstract SOM
abstract NeighborhoodFunction
type Gaussian_Neighborhood      <: NeighborhoodFunction end
type Ricker_Neighborhood        <: NeighborhoodFunction end
type Triangular_Neighborhood    <: NeighborhoodFunction end

function SOM__init__(this::SOM, x::Int, y::Int, input_len::Int;
                     σ::Float64=1.0, η::Float64=0.5, seed=0)
    this.t = 0
    this.epoch = 0
    this.λ = 0
    this.σ₀ = σ
    this.η₀ = η
    this.weights = Array{Float64}((x,y,input_len))
    this.activation_map = Array{Float64}((x,y))
    this.τ = τ_inverse
    this.ħ = Gaussian_Neighborhood()
    this.dist = euclidean
    this.seed = seed != 0 ? seed : round(Int, rand()*1e7)
    this.rng = MersenneTwister(seed)
    init_weights(this)
end

type GridSOM <: SOM
    weights::Array{Float64,3}
    activation_map::Array{Float64,2}
    η₀::Float64
    σ₀::Float64
    λ::Float64
    t::Int
    epoch::Int
    τ::Function
    ħ::NeighborhoodFunction
    dist::Function
    seed::Int
    rng::MersenneTwister

    function GridSOM(x::Int, y::Int, input_len::Int;
                 σ::Float64=1.0, η::Float64=0.5, seed=0)
        this = new()
        SOM__init__(this, x, y, input_len, σ=σ, η=η, seed=seed)
        return this
    end
end

type HexSOM <: SOM
    weights::Array{Float64,3}
    activation_map::Array{Float64,2}
    η₀::Float64
    σ₀::Float64
    λ::Float64
    t::Int
    epoch::Int
    τ::Function
    ħ::NeighborhoodFunction
    dist::Function
    seed::Int
    rng::MersenneTwister

    function HexSOM(x::Int, y::Int, input_len::Int;
                    σ::Float64=1.0, η::Float64=0.5, seed=0)
        this = new()
        SOM__init__(this, x, y, input_len, σ=σ, η=η, seed=seed)
        return this
    end
end

function size(som::SOM)
    return size(som.activation_map)
end

# --------------------------
# Neighboring units

function neighbor_units(som::GridSOM, x::Int, y::Int)
    r, c = size(som)
    neighbors = [
        (x-1,y-1), (x-1,y), (x-1,y+1),
         (x,y-1),            (x,y+1),
        (x+1,y-1), (x+1,y), (x+1,y+1)
    ]
    return neighbor_units(neighbors, r, c)
end

function neighbor_units(som::HexSOM, x::Int, y::Int)
    r, c = size(som)
    neighbors = [(n.x, n.y) for n in Hexagons.neighbors(HexagonAxial(x, y))]
    return neighbor_units(neighbors, r, c)
end

function neighbor_units(neighbors::Vector{Tuple{Int,Int}}, r, c)
    f(u) = 0 < u[1] ≤ r && 0 < u[2] ≤ c
    filter(f, neighbors)
end

neighbor_units(som::SOM, u::Tuple{Int,Int}) = neighbor_units(som, u...)

function hexagonal(u::Tuple{Int,Int}, v::Tuple{Int,Int})
    return float(Hexagons.distance(HexagonAxial(u...), HexagonAxial(v...)))
end

# --------------------------
# Decay functions

function τ_inverse(x::Float64, t::Int, λ::Float64)
    return x / (1 + t/λ)
end

function τ_exponential(x::Float64, t::Int, λ::Float64)
    return x * exp(-t/λ)
end

# --------------------------
# Neighborhood functions

function ħ(som::SOM, u::Tuple{Int,Int}, bmu::Tuple{Int,Int}, σ::Float64)
    return _ħ(som, som.ħ, u, bmu, σ)
end

# Gaussian
_ħ(som::SOM, ::Gaussian_Neighborhood, u, v, σ) = ħ_gaussian(som, u, v, σ)
ħ_gaussian(::GridSOM, u, v, σ) = ħ_gaussian(σ, euclidean(collect(u), collect(v)))
ħ_gaussian(::HexSOM, u, v, σ) = ħ_gaussian(σ, hexagonal(u, v))
ħ_gaussian(σ::Float64, d::Float64) = exp(-d^2 / (2 * σ^2))

# Mexican hat
_ħ(som::SOM, ::Ricker_Neighborhood, u, v, σ) = ħ_ricker(som, u, v, σ)
ħ_ricker(::GridSOM, u, v, σ) = ħ_ricker(σ, euclidean(collect(u), collect(v)))
ħ_ricker(::HexSOM, u, v, σ) = ħ_ricker(σ, hexagonal(u, v))
ħ_ricker(σ::Float64, d::Float64) = (1 - d^2 / σ^2) * ħ_gaussian(σ, d)

# Triangular
_ħ(som::SOM, ::Triangular_Neighborhood, u, v, σ) = ħ_triangular(som, u, v, σ)
ħ_triangular(::GridSOM, u, v, σ) = ħ_triangular(σ, euclidean(collect(u), collect(v)))
ħ_triangular(::HexSOM, u, v, σ) = ħ_triangular(σ, hexagonal(u, v))
ħ_triangular(σ::Float64, d::Float64) = abs(d) ≤ σ ? 1 - abs(d) / σ : 0.0

# --------------------------
# Getter / Setters

function get_unit_weight(som::SOM, i::Int, j::Int)
    return som.weights[i,j,:]
end

function get_unit_weight(som::SOM, c::Tuple{Int, Int})
    return get_unit_weight(som, c[1], c[2])
end

function set_unit_weight(som::SOM, i::Int, j::Int, w::Array{Float64})
    som.weights[i,j,:] = w
end

function set_unit_weight(som::SOM, i::Int, j::Int, w::Float64)
    som.weights[i,j,:] = [w]
end

function set_unit_weight(som::SOM, c::Tuple{Int, Int}, w)
    set_unit_weight(som, c[1], c[2], w)
end

# --------------------------
# Main API

function init_weights(som::SOM)
    l = size(som.weights, 3)
    for k in eachindex(som.activation_map)
        i, j = ind2sub(som.activation_map, k)
        w = rand(som.rng, l) * 2 - 1
        w /= norm(w)
        set_unit_weight(som, i, j, vec(w))
    end
end

function init_λ(som::SOM, num_iter::Int)
    som.λ = (som.t + num_iter) / 2
end

function update(som::SOM, input::Array, bmu::Tuple)
    som.t += 1
    input = vec(input)
    ηₜ = som.τ(som.η₀, som.t, som.λ)
    σₜ = som.τ(som.σ₀, som.t, som.λ)
    for k in eachindex(som.activation_map)
        u = ind2sub(som.activation_map, k)
        h = ħ(som, u, bmu, σₜ)
        weight = vec(get_unit_weight(som, u))
        weight = weight + h * ηₜ * (input - weight)
        set_unit_weight(som, u, weight)
    end
end

function get_BMU(som::SOM, input::Array)
    activate(som, input)
    return ind2sub(size(som), indmin(som.activation_map))
end

function activate(som::SOM, input::Array)
    input = reshape(input, (1,1,size(input, 2)))
    for k in eachindex(som.activation_map)
        i, j = ind2sub(som.activation_map, k)
        weight = get_unit_weight(som, i, j)
        som.activation_map[k] = som.dist(vec(input), vec(weight))
    end
end

function quantize(som::SOM, data::Array)
    q = similar(data)
    for i=1:size(data, 1)
        input = data[i,:]
        bmu = get_BMU(som, input)
        q[i,:] = get_unit_weight(som, bmu)
    end
    return q
end

function activation_response(som::SOM, data::Array)
    a = zeros(size(som))
    for i=1:size(data, 1)
        input = data[i,:]
        bmu = get_BMU(som, input)
        a[bmu...] += 1
    end
    return a
end

function umatrix(som::SOM)
    u_matrix = similar(som.activation_map)
    l = size(som)
    indices = [(x,y) for x in 1:l[1], y in 1:l[2]]
    for u in indices
        u_matrix[u...] = 0.0
        neighbors = neighbor_units(som, u)
        w = get_unit_weight(som, u)
        for n in neighbors
            v = get_unit_weight(som, n)
            u_matrix[u...] += som.dist(collect(v), collect(w))
        end
        u_matrix[u...] /= length(neighbors)
    end
    return u_matrix
end

function bmu_map(som::SOM, data::Array)
    bmus = DefaultDict(Tuple{Int,Int}, Vector{Vector}, Vector{Vector})
    for i=1:size(data, 1)
        input = data[i,:]
        bmu = get_BMU(som, input)
        push!(bmus[bmu...], vec(input))
    end
    return bmus
end

function sequential_random_epoch(som::SOM, data::Array, num_iter::Int)
    som.epoch += 1
    init_λ(som, num_iter)
    for t = 0:num_iter
        i = rand(som.rng, 1:size(data, 1))
        input = data[i,:]
        update(som, input, get_BMU(som, input))
    end
end

function sequential_epoch(som::SOM, data::Array, epochs=1)
    for i=1:epochs
        som.epoch += 1
        num_iter = size(data, 1)
        init_λ(som, num_iter)
        for t in shuffle(som.rng, collect(1:num_iter))
            input = data[t,:]
            update(som, input, get_BMU(som, input))
        end
    end
end

function quantization_error(som::SOM, data::Array)
    err = 0
    for i=1:size(data, 1)
        input = data[i,:]
        weight = get_unit_weight(som, get_BMU(som, input))
        err += som.dist(vec(weight), vec(input))
    end
    return err / size(data, 1)
end

function reset_state(som::SOM)
    som.t = 0
    som.epoch = 0
    som.rng = MersenneTwister(som.seed)
    init_weights(som)
end

end

max_stars_repo_name: agrimaldi/JSom

max_issues_repo_name: grmwld/JSom

max_forks_repo_name: agrimaldi/JSom


text:
"""
    CN0_monte_carlo(CN0, dopplers, doppler_rates, t_length, f_s, channel="I",
                    signal_type=define_l1ca_code_type(t_length); bins=1000,
                    iterations=100, prn=1, nADC=4, phase_noise_scaler=1/10,
                    include_phase_noise=true, include_thermal_noise=true,
                    include_databits_I=true, include_databits_Q=true,
                    f_if=0., include_adc=true, Tsys=535., include_carrier=true,
                    σω=1000., fd_range=5000., RLM=10, replica_t_length=1e-3,
                    cov_mult=1, q_a=1, q_mult=1, dynamickf=true, dll_b=5,
                    state_num=3, fine_acq_method=:fft, M=10)

Performs a benchmark test of signal simulation and processing by varying all
parameters, except the carrier-to-noise ratio, C/N₀, which is set by the user.

Required Arguments:

- `CN0`: carrier-to-noise ratio (C/N₀) of signal in dB
- `dopplers`: vector of Dopplers where duplicates are allowed
- `doppler_rates`: vector of Doppler rates where duplicates are allowed
- `t_length`: length of the signal to simulate and process in seconds
- `f_s`: sampling rate of signal in Hz

Optional Arguments:

- `channel`: channel to process signal on (must be either `"I"` or `"Q"`)
             `(default = "I")`
- `signal_type`: `SignalType` struct describing signal type
                 `(default = define_l1ca_code_type(t_length))`
- `bins`: number of bins that the Doppler vectors will be split into
          `(default = 1000)`
- `iterations`: number of times the simulation will be re-run and re-processed
                `(default = 100)`
- `prn`: PRN number to simulate `(default = 1)`
- `nADC`: bit depth of signal
- `phase_noise_scaler`: standard deviation of phase noise in rads
                       `(default = 1/10)`
- `include_phase_noise`: flag to include phase noise `(default = true)`
    * will be regenerated on each iteration
- `include_thermal_noise`: flag to include thermal noise `(default = true)`
    * will be regenerated on each iteration
- `include_databits_I`: flag to include databits on the I channel if they exist
                        `(default = true)`
- `include_databits_Q`: flag to include databits on the Q channel if they exist
                        `(default = true)`
- `f_if`: signal IF frequency in Hz `(default = 0Hz)`
- `include_adc`: flag to quantize signal `(default = true)`
- `Tsys`: system noise temperature in kelvin `(default = 535K)`
- `include_carrier`: flag to include carrier in signal `(default = true)`
- `σω`: initial uncertainty in Doppler rate `(default = 1000Hz/s)
- `fd_range`: range of Doppler above and below Doppler search center frequency
              in Hz `(default = 5000Hz)`
- `RLM`: a multiple of the replica signal length and used for FFT based fine
         acquisition `(default = 10)`
- `replica_t_length`: length of the replica signal in seconds `(default = 1e-3s)`
    * determines the length of the replica used in course acquisition and
      tracking
- `cov_mult`: scalar to inflate the initial covariance matrix `P₀`
              `(default = 1)`
- `q_a`: line of site platform dynamics in m²/s⁶ `(default = 1m²/s⁶)`
- `q_mult`: scalar to inflate the process noise matrix `Q` `(default = 1)`
- `dynamickf`: flag to specify if steady state KF gain is used or if KF can
               change from time step to time step `(default = true)`
- `dll_b`: bandwidth of the DLL filter in Hz `(default = 5Hz)`
- `state_num`: number of states to track `(default = 3)`
	* if set to `3`, carrier phase, Doppler, and Doppler rate are tracked
	* if set to `2`, only carrier phase and Doppler are tracked
- `fine_acq_method`:
- `M`: a multiple of the replica signal length and used for carrier based fine
       acquisition `(default = 10)`
- `h_parms`: 5 element vector containing oscillator h parameeters

Returns:

-
"""
function CN0_monte_carlo(CN0, dopplers, doppler_rates, t_length, f_s, channel="I",
                         signal_type=define_l1ca_code_type(t_length); bins=1000,
                         iterations=100, prn=1, nADC=4, phase_noise_scaler=1/10,
                         include_phase_noise=true, include_thermal_noise=true,
                         include_databits_I=true, include_databits_Q=true,
                         f_if=0., include_adc=true, Tsys=535., include_carrier=true,
                         σω=1000., fd_range=5000., RLM=10, replica_t_length=1e-3,
                         cov_mult=1, q_a=1, q_mult=1, dynamickf=true, dll_b=5,
                         state_num=3, fine_acq_method=:fft, M=10,
                         h_parms=h_parms_tcxo[1], σ_phi=π/2, fd_center=fd_center)
    if (channel != "I") && (channel != "Q")
        error("Invalid channel. Channel to be processed must be either I or Q.")
    end
    doppler_hist = fit(Histogram, dopplers, nbins=bins)
    doppler_weights = Weights(doppler_hist.weights)
    doppler_rate_hist = fit(Histogram, doppler_rates, nbins=bins)
    doppler_rate_weights = Weights(doppler_rate_hist.weights)
    Δfd = 1/replica_t_length
    if signal_type.include_I
        include_databits_I = signal_type.I_codes.databits
    else
        include_databits_I = false
    end
    if signal_type.include_Q
        include_databits_Q = signal_type.Q_codes.databits
    else
        include_databits_Q = false
    end
    include_databits = include_databits_I || include_databits_Q
    signal = definesignal(signal_type, f_s, t_length; prn=prn, CN0=CN0,
                          nADC=nADC, f_if=f_if, include_adc=include_adc,
                          Tsys=Tsys, include_carrier=include_carrier,
                          phase_noise_scaler=phase_noise_scaler,
                          include_phase_noise=include_phase_noise,
                          include_thermal_noise=include_thermal_noise,
                          include_databits_I=include_databits_I,
                          include_databits_Q=include_databits_Q,
                          skip_noise_generation=true, receiver_h_parms=h_parms)
    ϕ = rand(0:0.0001:2π, iterations)
    code_start_idx = rand(1:signal.sample_num, iterations)
    T_num = floor(Int, signal.t_length/replica_t_length)
    t = calctvector(T_num, Δfd)
    p = Progress(iterations, 1, "Processing...")
    for i in 1:iterations
        f_d = sample(doppler_hist.edges[1], doppler_weights)
        if ismissing(fd_center)
            fd_center = round(f_d/Δfd)*Δfd
        end
        fd_rate = sample(doppler_rate_hist.edges[1], doppler_rate_weights)
        definesignal!(signal; f_d=f_d, fd_rate=fd_rate, phi=ϕ[i],
                      code_start_idx=code_start_idx[i], new_databits=true,
                      new_thermal_noise=true, new_phase_noise=true,)
        generatesignal!(signal)
        acqresults, trackresults = process(signal, signal_type, prn, channel;
                                           show_plot=false, σω=σω, q_a=q_a,
                                           fd_range=fd_range, RLM=RLM,
                                           q_mult=q_mult, cov_mult=cov_mult,
                                           dynamickf=dynamickf, dll_b=dll_b,
                                           state_num=state_num, fd_rate=fd_rate,
                                           fine_acq_method=fine_acq_method, M=M,
                                           replica_t_length=replica_t_length,
                                           fd_center=fd_center,
                                           return_corrresult=false,
                                           fine_acq=fine_acq, σ_phi=σ_phi)
        # Truth data
        if channel == "I"
            init_code_phase = calcinitcodephase(signal_type.I_codes.code_lengths[1],
                                                signal.f_code_d_I[1],
                                                signal.f_code_dd_I[1],
                                                f_s, code_start_idx[i])
            truth_code = calccodeidx.(init_code_phase, signal.f_code_d_I[1],
                                      signal.f_code_dd_I[1], t,
                                      signal_type.I_codes.code_lengths[1])
            if signal.signal_type.I_codes.databits
                truth_databits = signal.signal_type.I_codes.codes[end][prn]
            else
                truth_databits = missing
            end
        elseif channel == "Q"
            init_code_phase = calcinitcodephase(signal_type.Q_codes.code_lengths[1],
                                                signal.f_code_d_Q[1],
                                                signal.f_code_dd_Q[1],
                                                f_s, code_start_idx[i])
            truth_code = calccodeidx.(init_code_phase, signal.f_code_d_Q[1],
                                      signal.f_code_dd_Q[1], t,
                                      signal_type.Q_codes.code_lengths[1])
            if signal.signal_type.Q_codes.databits
                truth_databits = signal.signal_type.Q_codes.codes[end][prn]
            else
                truth_databits = missing
            end
        else
            error("Invalid channel specified. Channel must be either I or Q.")
        end
        truth_doppler = f_d .+ fd_rate.*t
        next!(p)
    end
end


text:
using ArgParse
using EventStreamMatrix
using EventStreamGLM
using LinearAlgebra
using YAML
using GLM
using BSplines
import Dates: now, Millisecond

import Random: seed!
import IterativeSolvers: Identity
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--infolder"
            help="Folder to read in simulated dataframes from"
            default="/tmp/"
        "--outfolder"
            help="Folder to write out fitted models"
            default="/tmp/"
        "--seed"
            help="RNG Seed"
            default=564
            arg_type=Int
        "--name"
            help="Name to read in and write"
            default="sim"
        "--include_lag"
            help="Whether or not to use response variable lagged by one bin as additional input"
            action=:store_true
        "--dense"
            help="Whether to use a dense design matrix representation"
            action=:store_true
        "--nonconj"
            help="Whether to default to a non-CG fit"
            action=:store_true
        "--fineness"
            help="Level of discretization to use"
            default=0.1
            arg_type=Float64
        "--breakafter"
            help="Early stopping on file number"
            default=1000
            arg_type=Int
    end
    return parse_args(s)
end

function main()
    args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in args
        println("  $arg  =>  $val")
    end
    seed!(args["seed"])
    BLAS.set_num_threads(1)
    q=1
    for file in readdir(abspath(args["infolder"]))
        if q > args["breakafter"]
            break
        end
        if startswith(file, args["name"])
            println("Reading $file")
            dat = YAML.load_file(joinpath(abspath(args["infolder"]), file))
            fittime_start = now()
            outpoints = dat["output"]
            δ = args["fineness"]
            nsources=convert(Int, dat["nsources"][1])
            maxtime = convert(Float64, dat["maxtime"][1])
            order = convert(Int, dat["order"][1])
            eventstream = Tuple{Float64, Int}[]
            labels = ["source_$i" for i in 1:nsources]
            for (i, l) in enumerate(labels)
                eventstream = vcat(eventstream, [(t, i) for t in dat[l]]) :: Vector{Tuple{Float64, Int}}
            end
            if args["include_lag"]
                push!(labels, "response_lag")
                # Lag by
                lag_resp = outpoints .+ δ
                filter!(t -> t<=maxtime, lag_resp)
                eventstream = vcat(eventstream, [(t, nsources + 1) for t in lag_resp])
            end
            sort!(eventstream)

            # Construct  data
            basis = BSplineBasis(dat["order"], dat["breakpoints"])
            E = FirstOrderBSplineEventStreamMatrix(eventstream, labels, δ, maxtime, order, dat["breakpoints"], true)
            y_es = IdentityEventStreamVector(outpoints, δ, maxtime)
            y = Vector(y_es)
            y = y .>= 1 ## For binary response

            # Construct GLM
            cntrl = CGControl(0.001, 0.001, 40, false, true, Identity())
            if args["dense"] && !args["nonconj"]
                pp=DensePredConjGrad{Float64}(Matrix(E), zeros(Float64, size(E)[2]) , cntrl)
            elseif !args["dense"] && !args["nonconj"]
                pp = EventStreamPredConjGrad{Float64, String}(E, zeros(Float64, size(E)[2]), cntrl)
            else
                pp = GLM.DensePredChol(Matrix(E), false)
            end
            rr = GLM.GlmResp(y, Binomial(), LogitLink(), repeat([0.0], length(y)), repeat([1.0], length(y)))
            mod = GeneralizedLinearModel(rr, pp, false)

            # Fitting!
            println("Fitting $file")
            fit!(mod)
            fittime_end = now()
            fittime = Millisecond(fittime_end - fittime_start).value

            # Save output
            outdata = Dict{String, Vector{Float64}}("fineness" => [δ])
            outdata["coefs"] = coef(mod)
            outdata["fittime"] = [fittime]
            YAML.write_file(joinpath(abspath(args["outfolder"], file)), outdata)
            q+=1
        end
    end
end

main()


text:
# helpful bits for working with n \geq 2

"""
   uvec(x)

Helper to find a unit vector.
"""
function uvec(x)
    nm = norm(x)
    nm == 0 && return x
    return x/nm
end

"""
    `unzip(vs)`
    `unzip(v1, v2, ...)`
    `unzip(r::Function, a, b)`

Take a vector of points described by vectors (as returned by, say
`r(t)=[sin(t),cos(t)], r.([1,2,3])`, and return a tuple of collected x
values, y values, and optionally z values.

If the argument is specified as a comma separated collection of vectors, then these are combined and passed along.

If the argument is a function and two end point, then the function is
evaluated at 100 points between `a` and `b`.

This is useful for plotting when the data is more conveniently
represented in terms of vectors, but the plotting interface requires the x and y values collected.

Examples:
```
using Plots
r(t) = [sin(t), cos(t)]
rp(t) = [cos(t), -sin(t)]
plot(unzip(r, 0, 2pi)...)  # calls plot(xs, ys)

t0, t1 = pi/6, pi/4

p, v = r(t0), rp(t0)
plot!(unzip(p, p+v)...)  # connect p to p+v with line

p, v = r(t1), rp(t1)
quiver!(unzip([p])..., quiver=unzip([v]))
```

Based on `unzip` from the `Plots` package.
"""
unzip(vs::Vector) = Tuple([[vs[i][j] for i in eachindex(vs)] for j in eachindex(vs[1])])
function unzip(ws::Array; recursive=false)
    if recursive
        unzip([unzip(ws[:,j]) for j in 1:size(ws)[end]])
    else
        Tuple(eltype(first(ws))[xyz[j] for xyz in ws] for j in eachindex(first(ws)))
    end
end
#unzip(vs) = (A=hcat(vs...); Tuple([A[i,:] for i in eachindex(vs[1])]))
unzip(v,vs...) = unzip([v, vs...])
unzip(r::Function, a, b, n=100) = unzip(r.(range(a, stop=b, length=n)))

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

## The gradient in SymPy.
## We avoid name collision with  ForwardDiff.gradient we use `grad`
import ForwardDiff: gradient
gradient(ex::SymPy.Sym, vars::AbstractArray=free_symbols(ex)) = diff.(ex, [vars...])
gradient(f::Function) = (x, xs...) -> ForwardDiff.gradient(f, vcat(x, xs...))

"""
    curl(F)

Find curl of a 2 or 3-D vector field.
"""
function curl(J::Matrix)
    if size(J) == (2,2)
        Mx, Nx, My, Ny = J
        return Nx - My # a scalar
    elseif size(J) == (3,3)
        Mx, Nx, Px, My, Ny, Py, Mz, Nz, Pz = J
        return [Py-Nz, Mz-Px, Nx-My] # ∇×VF
    else
        throw(ArgumentError("Wrong size jacobian matrix for a curl"))
    end
end
curl(F::Vector{Sym}, vars=free_symbols(F)) = curl(F.jacobian(vars))
curl(F::Tuple) = curl(F[1], F[2])
curl(F::Function, pt) = curl(ForwardDiff.jacobian(F, float.(pt)))
curl(F::Function) = (pt, pts...) -> curl(F, vcat(pt, pts...))

"""
    divergence(F)

Find divergence of a 3-D vector vield.
"""
divergence(F::Vector{Sym}, vars=free_symbols(F)) = sum(diff.(F, vars))
divergence(F::Tuple) = divergence(F[1], F[2])
divergence(F::Function, pt) = sum(diag(ForwardDiff.jacobian(F, float.(pt))))
divergence(F::Function) = (pt, pts...) -> divergence(F, vcat(pt, pts...))

## Is this a bad idea?
## syntax is a bit heavy with parentheses

struct DelOperator end
const ∇ = DelOperator()

## For symbolic objects, these use `free_symbols(ex)` to find variables
## this *may* not be right if sorting isn't as desired, or variables don't
## appear in the formula
(::DelOperator)(f) = return gradient(f)
(::DelOperator)(f::Tuple) = gradient(f[1],f[2])
LinearAlgebra.dot(::DelOperator, F) = divergence(F)
LinearAlgebra.dot(::DelOperator, F, vars) = divergence(F, vars)

LinearAlgebra.cross(::DelOperator, F) = curl(F)
LinearAlgebra.cross(::DelOperator, F, vars) = curl(F, vars)

##


text:
function rotate(rot, chr::Char)
    if !isletter(chr)
        return chr
    end

    position = Int(lowercase(chr)) - 97
    newposition =  (position + rot) % 26
    rotated = Char(97 + newposition)
    chr == lowercase(chr) ? rotated : uppercase(rotated)
end

function rotate(rot, phrase::String)
    map(ch -> rotate(rot, ch), phrase)
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
using Documenter, KernelMachines, Literate

const example_dir = joinpath(@__DIR__, "..", "examples")
const generated_dir = joinpath(@__DIR__, "src", "generated")

for file in readdir(example_dir)
    Literate.markdown(joinpath(example_dir, file), generated_dir)
end

ENV["GKSwstype"] = "100"

makedocs(
    sitename="Kernel Machines",
    repo = "https://gitlab.com/VeosDigital/KernelMachines.jl/blob/{commit}{path}#{line}",
    pages = [
        "index.md",
        "machines.md",
        "kernels.md",
        "Examples" => [
            "generated/multiscale.md",
            "generated/polynomial.md",
        ]
    ])


text:
using YAXArrayBase, DimensionalData, AxisArrays, AxisIndices, Test

struct M
end
Base.ndims(::M) = 2
YAXArrayBase.getdata(::M)  = reshape(1:12,3,4)
YAXArrayBase.dimname(::M,i) = i==1 ? :x : :y
YAXArrayBase.dimvals(::M,i) = i==1 ? (0.5:1.0:2.5) : (1.5:0.5:3.0)
YAXArrayBase.getattributes(::M) = Dict{String,Any}("a1"=>5, "a2"=>"att")

@testset "AxisIndices" begin
    using AxisIndices: AxisIndices
    d = yaxconvert(AxisIndices.AxisArray,M())
    @test d isa AxisIndices.AxisArray
    @test getdata(d) == reshape(1:12,3,4)
    @test YAXArrayBase.dimnames(d) == (:Dim_1, :Dim_2)
    @test dimvals(d,1) == 0.5:1.0:2.5
    @test dimvals(d,2) == 1.5:0.5:3.0
end

@testset "AxisArrays" begin
    using AxisArrays: AxisArrays
    d = yaxconvert(AxisArrays.AxisArray,M())
    @test d isa AxisArrays.AxisArray
    @test getdata(d) == reshape(1:12,3,4)
    @test YAXArrayBase.dimnames(d) == (:x, :y)
    @test dimvals(d,1) == 0.5:1.0:2.5
    @test dimvals(d,2) == 1.5:0.5:3.0
end

@testset "NamedTuples" begin
    d = yaxconvert(NamedTuple,M())
    @test d isa NamedTuple
    @test getdata(d) == reshape(1:12,3,4)
    @test YAXArrayBase.dimnames(d) == (:x, :y)
    @test dimvals(d,1) == 0.5:1.0:2.5
    @test dimvals(d,2) == 1.5:0.5:3.0
end

@testset "NamedDims" begin
    using NamedDims: NamedDimsArray
    d = yaxconvert(NamedDimsArray,M())
    @test d isa NamedDimsArray
    @test getdata(d) == reshape(1:12,3,4)
    @test YAXArrayBase.dimnames(d) == (:x, :y)
end

@testset "DimensionalData" begin
    using DimensionalData
    d = yaxconvert(DimensionalArray,M())
    @test d isa DimensionalArray
    @test getdata(d) == reshape(1:12,3,4)
    @test YAXArrayBase.dimnames(d) == (:x, :y)
    @test dimvals(d,1) == 0.5:1.0:2.5
    @test dimvals(d,2) == 1.5:0.5:3.0
    @test getattributes(d) == Dict{String,Any}("a1"=>5, "a2"=>"att")
end

@testset "ArchGDAL" begin
  import Downloads
  p = Downloads.download("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")
  using ArchGDAL
  AG=ArchGDAL
  r = AG.readraster(p)
  @test YAXArrayBase.dimnames(r) == (:Y, :X, :Band)
  @test YAXArrayBase.dimname(r,1) == :Y
  @test YAXArrayBase.dimname(r,2) == :X
  @test YAXArrayBase.dimname(r,3) == :Band
  @test YAXArrayBase.dimvals(r,1) == -28493.166784412522:60.02213698319374:2298.189487965865
  @test YAXArrayBase.dimvals(r,2) == 4.2558845438021915e6:-60.02213698319374:4.22503316539283e6
  @test YAXArrayBase.dimvals(r,3) == ["Gray"]
  @test_throws Exception YAXArrayBase.dimname(r,4)
  @test_throws Exception YAXArrayBase.dimvals(r,4)
  @test YAXArrayBase.iscontdim(r,1) == true
  @test YAXArrayBase.iscontdim(r,2) == true
  @test YAXArrayBase.iscontdim(r,3) == true
  @test YAXArrayBase.getattributes(r)["projection_PROJ4"] == "+proj=cea +lat_ts=33.75 +lon_0=-117.333333333333 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs"
  b = AG.getband(r,1)
  @test YAXArrayBase.dimnames(b) == (:Y, :X)
  @test YAXArrayBase.dimname(b,1) == :Y
  @test YAXArrayBase.dimname(b,2) == :X
  @test YAXArrayBase.dimvals(b,1) == -28493.166784412522:60.02213698319374:2298.189487965865
  @test YAXArrayBase.dimvals(b,2) == 4.2558845438021915e6:-60.02213698319374:4.22503316539283e6
  @test_throws Exception YAXArrayBase.dimname(b,3)
  @test_throws Exception YAXArrayBase.dimvals(b,3)
  @test YAXArrayBase.iscontdim(b,1) == true
  @test YAXArrayBase.iscontdim(b,2) == true
  @test YAXArrayBase.getattributes(b)["projection_PROJ4"] == "+proj=cea +lat_ts=33.75 +lon_0=-117.333333333333 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs"
end


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

using Interpolations

# -----------------------------------------------------------------------------
#           AMPLITUDE POINT SPREAD FUNCTION BASE TYPE & METHODS
# -----------------------------------------------------------------------------

# Base Type
#
# All PointSpreadFunctionModel subtypes will inherit from this base type.
abstract type PointSpreadFunctionModel end

# In a lateral shift invariant models, the image at each z plane can
# be computed using a convolution with a PSF for that z-plane.
#
# Examples: DebyeModel(), VectorialDebyeModel(), and DefocusModel().
#
abstract type LateralShiftInvariantModel <: PointSpreadFunctionModel end

# In a periodic shift invariant models, the point spread function
# displays some periodic lateral shift invariance, but is not completey
# shift invariant.  Projection operators for such models can still be
# computed fairly efficiently using a sequence of convolution operations.
#
# Subtypes of this type must define 'shift_invariance_subpixels' and
# 'shift_invariance_dimensions' to indicate the periodicity of that
# particular model.
#
# Examples: LightFieldModel()
#
abstract type PeriodicShiftInvariantModel <: PointSpreadFunctionModel end

# Shift variant models cannot be projected using convolution, since
# they allow for PSFs that change with every 3D point location.
#
# Examples: CodedApertureLightFieldModel()
#
abstract type ShiftVariantModel <: PointSpreadFunctionModel end

# --------------
# Shared Methods
#
function num_observations(model::PointSpreadFunctionModel)
    return Int64(model.num_observations)
end

function atf(model::PointSpreadFunctionModel, p, theta = Float64[]; na_limit = true, apply_tilt = true)
    error("The atf() method is not implemented for this model type.")
end

function otf_2d(model::PointSpreadFunctionModel, p, theta = Float64[])
    error("The atf() method is not implemented for this model type.")
end

function apsf(model::PointSpreadFunctionModel, p, theta = Float64[])
    error("The apsf() method is not implemented for this model type.")
end

function psf(model::PointSpreadFunctionModel, p, theta = Float64[])
    wf = WaveOptics.apsf(model, p, theta)
    return abs2.(wf)
end

"""
Estimate the diameter (in meters) of an image large enough to contain the point
spread function for a wide field microscope without any phase mask.  Most
(but not all) point spread functions with back aperture masks will fit in
this diameter.

Here, p = (x,y,z) is a tuple containing the point generating this PSF.
"""
function estimate_psf_size(model::PointSpreadFunctionModel, p)

        # Check length and unpack p
        assert(length(p) == 3)
        x,y,z = p

        objective_theta = asin(model.recipe.objective_na / model.recipe.medium_index)
        aperture_diameter = abs( 2 * z * tan(objective_theta) )
        rayleigh_limit = model.recipe.wavelength / (2 * model.recipe.objective_na)

        # Choose a PSF size that comfortably fits the full PSF, with special
        # handling for when z = 0. (The min size is 8 * rayleigh_limit,
        # otherwise it is 2 * aperture_diameter.)
        final_diameter_m = max(8*rayleigh_limit, 2*aperture_diameter)
        return final_diameter_m
end

function optimal_sampling(model::PointSpreadFunctionModel, sim_size_m; oversampling = nothing)
        """
        Given a desired simulation size, this function returns a sampling
        rate (in pixels) that will yield an accurate wave optics simulation.

        Note that this sampling rate is only an estimate based on the abbe
        limit. It may be too high in some cases, and too low in others. Use it
        as a starting point.
        """

        if oversampling == nothing
            oversampling = 1
        end

        # We use the abbe limit as an estimate of the optimal sampling rate. We
        # recommend sampling at 2x the Abbe limit.
        abbe_limit = model.recipe.wavelength / (2. * model.recipe.objective_na);
    return sim_size_m / (2*abbe_limit) * oversampling
end

# --------------------------------------------------------------------------
#                COMMON PHASE MASKS & APSF UTILITIES
# --------------------------------------------------------------------------

"""
This function takes an APSF lookup up table computed using apsf_scalar_debye
or apsf_defocus_atf() and uses the translational invariance property of the
APSF to interpolate and place it at any location x, y at the depth z for
which the apsf_template was computed.

Note, the translation could have been accomplished with a tilt in the back
aperture plane via translation_mask(), but this would have potentially led
to sampling problems if a large tilt (due to a large x or y translation)
exceeded the nyquist sampling rate of the back aperture. For general
translations which may be large, this method is more robust.
"""
function translate_apsf(sim_size_m, sim_size_px, apsf::AbstractArray{T, 2}, dx, dy) where T<:Number
    sample_period = sim_size_m / sim_size_px

    # Use 2D interpolation to shift the translationally invariant microscope PSF
    apsf_real_interp = interpolate(real(apsf), BSpline(Interpolations.Linear()))
    apsf_imag_interp = interpolate(imag(apsf), BSpline(Interpolations.Linear()))

    result = zeros(eltype(apsf), size(apsf))
    for c in 1:size(apsf,2)
        for r in 1:size(apsf,1)
            result[r,c] = apsf_real_interp(r + dy/sample_period, c - dx/sample_period) + im * apsf_imag_interp(r + dy/sample_period, c - dx/sample_period)
        end
    end
    return result
end

"""
Creates a phase mask in the back aperture plane corresponding to a
translation in the image plane. This translation mask is based on the
Fourier shift theorem. For small shifts, it should match with the
translate_apsf() function below.

fx, fy must be in 'normalized' cartesian coordinates (see fourier_plane_meshgrid()).
"""
function translation_mask(dx, dy, fx::Array{T, 2}, fy::Array{T, 2},
                          objective_na::Float64,
                          wavelength::Float64) where T<:AbstractFloat
    result = Array{ComplexF32}(size(fx))
    omega = -2. * im * pi * objective_na/wavelength
    for i in eachindex(result)
        result[i] = exp(omega * (fx[i] * dx + fy[i] * (-dy)))
    end
    return result
end

"""
Implements the defocus phase function described in

Hanser, B. M., Gustafsson, M., & Agard, D. A. (2004). Phase-retrieved pupil
functions in wide-field fluorescence microscopy. Journal of Microscopy. Vol.
216. pp. 32-48. 2004.

This phase function models roughly what happens if the microscope objective
were to be refocused to a different depth 'dz' in the object.

r, theta must be in 'normalized' polar coordinates (see fourier_plane_meshgrid()).
"""
function defocus_mask(dz::S, r::AbstractArray{T,2}, theta::AbstractArray{T,2},
                      objective_na::Float64, wavelength::Float64, medium_index::Float64) where T<:AbstractFloat where S<:Number
    M,N = size(r)
    k = medium_index/wavelength
    omega = T(k^2) .- T(objective_na/wavelength)^2 * r.^2

    # Ensure we don't try to take the sqrt() of a negative number. This should
    # only happen outside the boundary of the aperture, so this check amounts to
    # zeroing out the phase of mask pixels outside the numerical aperture.
    binary_mask = omega .>= 0

    # Since CuArrays.jl doesn't have a function to compute exp(ComplexFloat),
    # we manually compute it.
    phi = T(-2π * dz) * sqrt.(max.(omega, 0))
    phase =  (cos.(phi) + 1im .* sin.(phi)) .* binary_mask
    return phase
end

# --------------------------------------------------------------------------
#                       APODIZATION FUNCTIONS
# --------------------------------------------------------------------------

"""
Apply a the apodization function for the Abbe sine condition. This
implements eqn. 6.3.7 from Min Gu's "Advanced Optical Imaging Theory"

r, theta must be in "normalized" polar coordinates (see fourier_plane_meshgrid()).
"""
function abbe_sine_apodization(rho::Array{T}, theta::Array{T}, objective_na::AbstractFloat, medium_index::AbstractFloat) where T<:AbstractFloat

    result = zeros(ComplexF32, size(rho))
    omega = objective_na/medium_index
    for i in eachindex(rho)
        if rho[i] <= 1.0
            alpha = asin(omega * rho[i])    # Gross eq. 21-3 and 21-4
            apodization = sqrt(cos(alpha))                # Gu eq. 6.3.10
            result[i] = apodization + 0. * im
        else
            result[i] = 0.
        end
    end
    return result
end


text:
using ProgressBars
using SparseArrays

struct Graph
    transition_matrix::SparseMatrixCSC{Int64, Int64}
    nodes::Vector{Int64}
    edges::Vector{Tuple{Int64, Int64}}
    degree_out::Vector{Int64}
    degree_in::Vector{Int64}
    n::Int32
    function Graph(T, nodes, edges, degree_out, degree_in, n)
        new(T, nodes, edges, degree_out, degree_in, n)
    end
    function Graph(filename)
        # Read edges and nodes from filename
        io = open(filename)
        lines = readlines(io)

        edges = Vector{Tuple{Int64, Int64}}()

        I = Vector{Int64}()
        J = Vector{Int64}()

        println("Loading graph in memory")
        maxI = 0
        maxJ = 0
        for (index, line) in ProgressBar(enumerate(lines))
            if !startswith(line, "#") && line != ""
                r = r"(\w+)[\t\s](\w+)"
                i1, i2 = match(r, line).captures
                i1 = parse(Int64, i1)
                i2 = parse(Int64, i2)
                push!(I, i1)
                push!(J, i2)

                push!(edges, (i1, i2))

                if i1 > maxI
                    maxI = i1
                end

                if i2 > maxJ
                    maxJ = i2
                end
            end
        end

        n = max(maxI, maxJ)

        nodes = Vector{Int64}(1:n)

        degree_out = zeros(n)
        degree_in = zeros(n)
        for i in 1:length(I)
            u = I[i]
            v = J[i]
            degree_out[u] += 1
            degree_in[v] += 1
        end

        V = ones(Int64, length(I))
        transition_matrix = sparse(I, J, V, n, n)
        new(transition_matrix, nodes, edges, degree_out, degree_in, n)
    end

end

function remove_node(graph::Graph, node::Int64)

end

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes