text:
using Test
using ConvDiffMIPDECO
using jInv.ForwardShare
using jInv.LinearSolvers
using KrylovMethods
using jInv.Mesh
using MUMPSjInv
using Statistics
using LinearAlgebra

domain = [0 1. 0 3. 0 2.]

N = [8;16;32;]
v = X -> [ones(size(X,1)) zeros(size(X,1)) zeros(size(X,1))]
bc = (:dir,:neu,:neu,:neu,:neu,:neu)

u  = X -> cos.(pi*X[:,1]).*cos.(pi*X[:,2]).*cos.(pi*X[:,3])
ux = X -> -pi*sin.(pi*X[:,1]).*cos.(pi*X[:,2]).*cos.(pi*X[:,3])
uy = X -> -pi*cos.(pi*X[:,1]).*sin.(pi*X[:,2]).*cos.(pi*X[:,3])
uz = X -> -pi*cos.(pi*X[:,1]).*cos.(pi*X[:,2]).*sin.(pi*X[:,3])
sig = .01;

f  = X->  3*sig*pi^2*u(X)+ux(X)

err = zeros(length(N))

times = zeros(length(N),3)
pFor=1;
for j=1:2
	for k=1:length(N)
		Ainvs = (getJuliaSolver(), getMUMPSsolver())

		nk  = [N[k];N[k];N[k]]
		Mk  = getRegularMesh(domain,nk)
		xck = getCellCenteredGrid(Mk)
		xnk = getNodalGrid(Mk)
		pFork = getConvDiffFEMParam(Mk,sig=sig,gd=u,bc=bc,Ainv=Ainvs[j])

		dobsk,pFork = getData(f(xck), pFork)
		pFork.Fields .-= mean(pFork.Fields)
		utruek = u(xnk)
		utruek .-= mean(utruek)

		 times[k,j]= @elapsed begin
			 err[k] = norm(pFork.Fields - utruek,Inf)/norm(utruek,Inf)
		 end
		println(err[k])
	end
	@test all(abs.(diff(log2.(err))).>1.4)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
"""
    couple_models(aero::QuasiSteady, stru::LiftingLineSection)

Create an aerostructural model using a quasi-steady aerodynamics model and a
lifting line section model.  The existence of this coupling allows
[`QuasiSteady`](@ref) to be used with [`LiftingLine`](@ref). This model
introduces the freestream air density ``\\rho`` as an additional parameter.
"""
couple_models(aero::QuasiSteady, stru::LiftingLineSection) = (aero, stru)

# --- Traits --- #

# steady
number_of_additional_parameters(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = 1
coupling_inplaceness(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = Nonlinear()
copuling_time_gradient_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection}) = Zeros()

# quasisteady without acceleration terms
number_of_additional_parameters(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = 1
coupling_inplaceness(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = Zeros()
coupling_state_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = Nonlinear()
copuling_time_gradient_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection}) = Zeros()

# quasisteady with acceleration terms
number_of_additional_parameters(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = 1
coupling_inplaceness(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = OutOfPlace()
coupling_rate_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = Linear()
coupling_state_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = Nonlinear()
copuling_time_gradient_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection}) = Zeros()

# --- Methods --- #

# steady
function get_coupling_inputs(aero::QuasiSteady{0}, stru::LiftingLineSection, dx, x, p, t)
    # extract rate variables
    dvx, dvy, dvz, dωx, dωy, dωz = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    # calculate loads
    N, A, M = quasisteady0_loads(a, b, ρ, a0, α0, cd0, cm0, u, v)
    # loads per unit span
    f = SVector(A, 0, N)
    m = SVector(0, M, 0)
    # return inputs
    return SVector(f..., m...)
end

# quasisteady without accelaration terms
function get_coupling_inputs(aero::QuasiSteady{1}, stru::LiftingLineSection, dx, x, p, t)
    # extract rate variables
    vx, vy, vz, ωx, ωy, ωz = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    # calculate aerodynamic loads
    N, A, M = quasisteady1_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    # loads per unit span
    f = SVector(A, 0, N)
    m = SVector(0, M, 0)
    # return inputs
    return SVector(f..., m...)
end

# quasisteady with acceleration terms
function get_coupling_inputs(aero::QuasiSteady{2}, stru::LiftingLineSection, dx, x, p, t)
    # extract rate variables
    dvx, dvy, dvz, dωx, dωy, dωz = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    udot, vdot, ωdot = liftingline_accelerations(dvx, dvz, dωy)
    # calculate aerodynamic loads
    N, A, M = quasisteady2_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # loads per unit span
    f = SVector(A, 0, N)
    m = SVector(0, M, 0)
    # return inputs
    return SVector(f..., m...)
end

# --- Performance Overloads --- #

# TODO: State rate, state, and parameter jacobians

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::QuasiSteady, stru::LiftingLineSection; rho)

    padd[1] = rho

    return padd
end

function separate_additional_parameters(aero::QuasiSteady, stru::LiftingLineSection, padd)

    return (rho = padd[1],)
end

ext: jl


text:
module CUDAint

using ..CUDA
using CUDA: CUDNN

import ..Flux: Flux
import Zygote
using Zygote: @adjoint

include("cudnn.jl")

end


text:
using Tokenize
import Tokenize.Tokens: kind, exactkind, startpos

# Defines a convenient type alias
const Tokenized = Array{Tokenize.Tokens.Token,1}
#+Tokenizer,Internal
# Defines a convenient method to tokenize a =String=
tokenized(s::String) = collect(tokenize(s))

# +Tokenizer,Internal
is_structure(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.STRUCT)
is_function(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.FUNCTION)
is_comment(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.COMMENT)
is_whitespace(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.WHITESPACE)
# +Tokenizer,Internal
# Checks for "\n"
is_whitespace_n(tok::Tokenized,idx::Int)::Bool = is_whitespace(tok,idx) && (untokenize(tok[idx])=="\n")
# +Tokenizer,Internal
# Checks for ENDMARKER (end of code)
is_endmarker(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.ENDMARKER)
is_identifier(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.IDENTIFIER)
is_declaration(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.DECLARATION)
is_dot(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.DOT)
is_comma(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.COMMA)
is_where(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.WHERE)
is_struct(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.STRUCT)
is_mutable(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.MUTABLE)
is_const(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.CONST)
is_global(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.GLOBAL)
is_local(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.LOCAL)
is_macro(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.MACRO)
is_issubtype(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.ISSUBTYPE)
is_export(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.EXPORT)
is_abstract(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.ABSTRACT)
is_type(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.TYPE)
is_end(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.END)
is_if(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.IF)
is_for(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.FOR)
is_while(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.WHILE)
is_begin(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.BEGIN)
is_try(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (exactkind(tok[idx])==Tokenize.Tokens.TRY)

# +Tokenizer,Internal
#
# Check if tok[i] points on
#
# #+BEGIN_SRC julia :eval never :exports code
# @enum ...
# #+END_SRC
#
function is_enum(tok::Tokenized,idx::Int)::Bool
    return (idx+1<=length(tok)) &&
        (exactkind(tok[idx])==Tokenize.Tokens.AT_SIGN) &&
        (untokenize(tok[idx+1])=="enum")
end

# +Tokenizer,Internal
is_opening_parenthesis(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.LPAREN)
is_closing_parenthesis(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.RPAREN)
is_opening_brace(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.LBRACE)
is_closing_brace(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.RBRACE)
is_opening_square(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.LSQUARE)
is_closing_square(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok)) && (kind(tok[idx])==Tokenize.Tokens.RSQUARE)

function is_opening_block(tok::Tokenized,idx::Int)::Bool
    return (idx<=length(tok))&&(is_structure(tok,idx)||
                                is_function(tok,idx)||
                                is_if(tok,idx)||
                                is_for(tok,idx)||
                                is_begin(tok,idx)||
                                is_try(tok,idx))
end

is_closing_block(tok::Tokenized,idx::Int)::Bool = (idx<=length(tok))&&(is_end(tok,idx))

#+Tokenizer,Internal,Obsolete
# Skip comment
# -> replaced by skip uninformative
function skip_comment(tok::Tokenized,idx::Int)
    while (idx<=length(tok))&&(is_comment(tok,idx)||is_whitespace(tok,idx))
        idx=idx+1
    end
    return idx
end

#+Tokenizer,Internal
# Skip whitespace
function skip_uninformative(tok::Tokenized,idx::Int)
    while (idx<=length(tok))&&(is_comment(tok,idx)||is_whitespace(tok,idx))
        idx=idx+1
    end
    return idx
end

#+Tokenizer,Internal
# Skip comment + whitespace
function skip_whitespace(tok::Tokenized,idx::Int)
    while (idx<=length(tok))&&is_whitespace(tok,idx)
        idx=idx+1
    end
    return idx
end

#<Tokenizer,Internal
# Check "strict" whitespace: "  \t    "
#
function is_whitespace_strict(s::String)::Bool
    n=length(s)
    for i = 1:n
        if (s[i]!=' ')&&(s[i]!='\t')
            return false
        end
    end
    return true
end
function is_whitespace_strict(tok::Tokenized,idx::Int)::Bool
    return is_whitespace(tok,idx)&&is_whitespace_strict(untokenize(tok[idx]))
end
#>

#+Tokenizer,Internal
# Skip whitespace
function skip_whitespace_strict(tok::Tokenized,idx::Int)
    while (idx<=length(tok))&&(is_whitespace_strict(tok,idx))
        idx=idx+1
    end
    return idx
end

#+Internal
#
function find_closing_X_helper(tok::Tokenized,idx::Int,
                               is_opening_X::Function,
                               is_closing_X::Function,
                               count_opening::Int = 0)::Int

    # @assert is_opening_X(tok,idx)

    # count_opening = 0
    while idx<=length(tok)
        if is_opening_X(tok,idx)
            count_opening =  count_opening+1
	    idx=idx+1
            continue
        end

        if is_closing_X(tok,idx)
            count_opening = count_opening-1
            if count_opening == 0
                return idx
            end
        end
        idx=idx+1
    end

    @assert false "no closing X"
end

#+
# assume an open parenthesis, find the closed one
function find_closing_parenthesis(tok::Tokenized,idx::Int)::Int
    @assert is_opening_parenthesis(tok,idx)
    return find_closing_X_helper(tok,idx,
                                 is_opening_parenthesis,
                                 is_closing_parenthesis)
end

#+
function find_closing_brace(tok::Tokenized,idx::Int)::Int
    @assert is_opening_brace(tok,idx)
    return find_closing_X_helper(tok,idx,
                                 is_opening_brace,
                                 is_closing_brace)
end

#+
function find_closing_square(tok::Tokenized,idx::Int)::Int
    @assert is_opening_square(tok,idx)
    return find_closing_X_helper(tok,idx,
                                 is_opening_square,
                                 is_closing_square)
end

#+Internal
#
# Find closing block (-> END)
#
# This function does not impose that the idx position is an opening
# block. If this is the case use the count_opening::Int = 0 default
# value.
#
# If you want to reach end of block starting from an already opened
# one (for instance from the middle of a function body, use the
# count_opening::Int = 1 value.
#
function find_closing_block(tok::Tokenized,idx::Int,
                            count_opening::Int = 0)::Int
    return find_closing_X_helper(tok,idx,
                                 is_opening_block,
                                 is_closing_block,
                                 count_opening)
end

#+Tokenizer,Internal
# Skips uniformative & identifier A or A.B.C or A.B{...}
# Returns (name,idx) if prod
# Does not move is identifier not found
#
function skip_identifier(tok::Tokenized,idx::Int;
                         identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)
    if !is_identifier(tok,idx)
        return idx_save
    end

    n=length(tok)
    identifier[] = untokenize(tok[idx])
    while idx<=n
        idx=idx+1
        idx_save = idx
        idx=skip_whitespace_strict(tok,idx)
        if is_dot(tok,idx)
            idx=idx+1
            idx=skip_whitespace_strict(tok,idx)
            if is_identifier(tok,idx)
                identifier[]=identifier[]*"."*untokenize(tok[idx])
                continue
            else
                identifier[]=""
                return idx_save # failure
            end
        end

        if is_opening_brace(tok,idx)
            idx=find_closing_brace(tok,idx)
            return idx+1
        end

        return idx_save
    end

    identifier[]=""
    return idx_save
end

# function skip_identifier(tok::Tokenized,idx::Int)::Int
#     return last(get_and_skip_identifier(tok,idx))
# end

#+Tokenizer,Internal
# Skips uninformative & a comma separated sequence of identifiers A,A.B.C,A.B{...}
# Does not move is comma separeted identifier not found
function skip_comma_separated_identifiers(tok::Tokenized,idx::Int)::Int
    idx_save = idx
    idx = skip_uninformative(tok,idx)
    if !is_identifier(tok,idx)
        return idx_save
    end

    n=length(tok)
    while idx<=n
        idx=skip_identifier(tok,idx)
        idx_save=idx
        idx=skip_uninformative(tok,idx)
        if !is_comma(tok,idx)
            return idx_save
        end
        idx=idx+1
        idx=skip_uninformative(tok,idx)
        if !is_identifier(tok,idx)
            return idx_save
        end
    end
    @assert false
end

#+Tokenizer,Internal
# Skips uninformative & a where block (where {...} or where A,B
# Does not move is comma separeted identifier not found
function skip_where_block(tok::Tokenized,idx::Int)::Int
    idx_save = idx
    idx = skip_uninformative(tok,idx)
    if !is_where(tok,idx)
        return idx_save
    end

    idx=idx+1
    idx = skip_uninformative(tok,idx)

    if is_opening_brace(tok,idx)
        return find_closing_brace(tok,idx)+1
    end

    if is_identifier(tok,idx)
        return skip_comma_separated_identifiers(tok,idx)
    end

    return idx_save
end

#+Tokenizer,Internal
# Skips uninformative & a declaration block ::identifier{...}
# Does not move in case of identification failure
function skip_declaration_block(tok::Tokenized,idx::Int)::Int
    idx_save = idx
    idx = skip_uninformative(tok,idx)
    if !is_declaration(tok,idx)
        return idx_save
    end

    idx=idx+1
    idx=skip_uninformative(tok,idx)

    if is_identifier(tok,idx)
        return skip_identifier(tok,idx)
    end

    return idx_save
end
#+Tokenizer,Internal
# Skips uninformative & a function call block identifier{...}(....)
# Does not move in case of identification failure
function skip_function_call_block(tok::Tokenized,idx::Int;
                         identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)
    if !is_identifier(tok,idx)
        return idx_save
    end

    idx=skip_identifier(tok,idx,identifier=identifier)
    idx=skip_uninformative(tok,idx)

    if is_opening_parenthesis(tok,idx)
        idx=find_closing_parenthesis(tok,idx)
        idx=idx+1
        return idx
    end

    return idx_save
end

#+Tokenizer,Internal
# Skips uniformative & issubtype indentifier
# Does not move is identifier not found
#
function skip_issubtype_block(tok::Tokenized,idx::Int)::Int
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    if !is_issubtype(tok,idx)
        return idx_save
    end
    idx=idx+1
    idx_check_success = idx
    idx = skip_identifier(tok,idx)

    if idx != idx_check_success
        return idx
    end

    return idx_save
end

#+Tokenizer,Internal
# Skips uniformative & structure block
# Does not move is identifier not found
#
function skip_struct_block(tok::Tokenized,idx::Int;
                          identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    check_is_structure = is_structure(tok,idx)
    check_is_mutable = is_mutable(tok,idx)

    if !(check_is_structure||check_is_mutable)
        return idx_save
    end

    if check_is_mutable
        idx=idx+1
        idx = skip_uninformative(tok,idx)

        if !is_structure(tok,idx)
            return idx_save
        end
    end

    @assert is_structure(tok,idx)

    idx=idx+1
    idx = skip_identifier(tok,idx,identifier=identifier)
    idx = skip_issubtype_block(tok,idx)

    return idx
end

#+Tokenizer,Internal
# Skips uniformative & abstract type block
# Does not move is identifier not found
#
function skip_abstract_block(tok::Tokenized,idx::Int;
                             identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    if is_abstract(tok,idx)

        idx = idx+1
        idx = skip_uninformative(tok,idx)

        if is_type(tok,idx)
            idx = idx + 1
            idx_check_success = idx
            idx = skip_identifier(tok,idx,identifier=identifier)

            if idx_check_success != idx
                # no need to check success as subtype is not mandatory
                idx = idx + 1
                idx = skip_issubtype_block(tok,idx)

                # success
                return idx
            end
        end
    end

    return idx_save
end

#+Tokenizer,Internal
# Moves idx until it reaches a comment #, end of line \n or of file ENDMARKER
#
function skip_line(tok::Tokenized,idx::Int)::Int
    while !(is_comment(tok,idx)||is_endmarker(tok,idx)||is_whitespace_n(tok,idx))
        idx=idx+1
    end
    return idx
end

#+Tokenizer,Internal
# Skips uniformative & @enum ... line
# Does not move is identifier not found
#
function skip_enum_block(tok::Tokenized,idx::Int;
                         identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    if is_enum(tok,idx)

        idx = idx+2 # caveat skip @ AND enum
        idx_check_success = idx
        idx = skip_identifier(tok,idx,identifier=identifier)

        if idx_check_success != idx
            # move until end of line (ignoring comment if any)
            idx = skip_line(tok,idx)
            return idx
        end
    end

    return idx_save
end

#+Tokenizer,Internal
# Skips uniformative & const A
# Does not move is identifier not found
#
function skip_variable_block(tok::Tokenized,idx::Int;
                         identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    if is_const(tok,idx)
        idx =  skip_uninformative(tok,idx+1)
    end

    if is_global(tok,idx)||is_local(tok,idx)
        idx =  skip_uninformative(tok,idx+1)
    end

    idx_check_success = idx
    idx = skip_identifier(tok,idx,identifier=identifier)

    if idx_check_success != idx
        return idx
    end

    return idx_save
end

#+Tokenizer,Internal
# Skips uniformative & @enum ... line
# Does not move is identifier not found
#
function skip_macro_block(tok::Tokenized,idx::Int;
                         identifier::Ref{String}=Ref{String}(""))::Int
    identifier[]=""
    idx_save = idx
    idx = skip_uninformative(tok,idx)

    if is_macro(tok,idx)
        idx=idx+1
        idx_check_success = idx
        idx = skip_identifier(tok,idx,identifier=identifier)

        if idx_check_success != idx
            idx=skip_uninformative(tok,idx)

            if is_opening_parenthesis(tok,idx)
                idx=find_closing_parenthesis(tok,idx)
                idx=idx+1
                return idx
            end
        end
    end
    return idx_save
end


text:
@safetestset "CanonicalSymplecticMatrix" begin
    using PoincareInvariants: CanonicalSymplecticMatrix
    using LinearAlgebra: dot

    @test CanonicalSymplecticMatrix{Int}(2) == [0 -1; 1 0]
    @test CanonicalSymplecticMatrix{Int}(4) == [0 0 -1 0; 0 0 0 -1; 1 0 0 0; 0 1 0 0]
    @test CanonicalSymplecticMatrix{Int}(6) == [0  0  0 -1  0  0;
                                                0  0  0  0 -1  0;
                                                0  0  0  0  0 -1;
                                                1  0  0  0  0  0;
                                                0  1  0  0  0  0;
                                                0  0  1  0  0  0]

    function test_canonicalmatrix(mat, n, ::Type{T}) where T
        mid = n ÷ 2

        for i in 1:n, j in 1:n
            if i == j - mid
                mat[i, j] === T(-1) || return false
            elseif i - mid == j
                mat[i, j] === T(1) || return false
            else
                mat[i, j] === T(0) || return false
            end
        end

        return true
    end

    @testset "CanonicalSymplecticMatrix{$T}($n)" for T in [Int, Float64], n in [10, 100, 1000]
        @test CanonicalSymplecticMatrix(n) === CanonicalSymplecticMatrix{Int}(n)

        C = CanonicalSymplecticMatrix{T}(n)
        @test size(C) == (n, n)
        @test eltype(C) == T

        mid = n ÷ 2

        @test test_canonicalmatrix(C, n, T)

        @test C * collect(1:n) == [(-mid-1:-1:-n)..., 1:mid...]

        @test dot(ones(n), C, ones(n)) == 0

        @test_throws ArgumentError CanonicalSymplecticMatrix{T}(n + 1)
    end
end


text:
# data/z-test.jl

using Test

@test extrema(ir_load_brainweb_t1_256()) == (0.0f0, 241.0f0)


text:
using Printf

abstract type NonlinearSolver{T} end

solve!(s::NonlinearSolver) = error("solve! not implemented for $(typeof(s))")

function solve!(s::NonlinearSolver{T}, x₀::Vector{T}) where {T}
    setInitialConditions!(s, x₀)
    solve!(s)
end

# default parameters for nonlinear solvers
const DEFAULT_nwarn=100

struct NonlinearSolverParameters{T}
    nmin::Int   # minimum number of iterations
    nmax::Int   # maximum number of iterations

    atol::T     # absolute tolerance
    rtol::T     # relative tolerance
    stol::T     # successive tolerance

    atol²::T
    rtol²::T
    stol²::T

    atol_break::T
    rtol_break::T
    stol_break::T

    function NonlinearSolverParameters{T}(nmin, nmax, atol, rtol, stol, atol_break, rtol_break, stol_break) where {T}
        @assert nmin ≥ 0
        @assert nmax > 0
        @assert atol > 0
        @assert rtol > 0
        @assert stol > 0
        @assert atol_break > 0
        @assert rtol_break > 0
        @assert stol_break > 0

        new(nmin, nmax, atol, rtol, stol, atol^2, rtol^2, stol^2, atol_break, rtol_break, stol_break)
    end
end

function NonlinearSolverParameters(T)
    NonlinearSolverParameters{T}(get_config(:nls_nmin),
                                 get_config(:nls_nmax),
                                 get_config(:nls_atol),
                                 get_config(:nls_rtol),
                                 get_config(:nls_stol),
                                 get_config(:nls_atol_break),
                                 get_config(:nls_rtol_break),
                                 get_config(:nls_stol_break))
end

mutable struct NonlinearSolverStatus{T}
    i::Int         # iteration number
    rₐ::T          # residual (absolute)
    rᵣ::T          # residual (relative)
    rₛ::T          # residual (successive)
    r₀::Vector{T}  # initial residual (absolute)

    x₀::Vector{T}  # initial solution
    xₚ::Vector{T}  # previous solution

    y₀::Vector{T}  # initial function
    yₚ::Vector{T}  # previous function

    NonlinearSolverStatus{T}(n) where {T} = new(0, 0, 0, 0, zeros(T,n),
                                                zeros(T,n), zeros(T,n),
                                                zeros(T,n), zeros(T,n),)
end

Base.show(io::IO, status::NonlinearSolverStatus) = print(io,
                        (@sprintf "    n=%4i" status.i),  ",   ", (@sprintf "rₐ=%14.8e" status.rₐ), ",   ",
                        (@sprintf "rᵣ=%14.8e" status.rᵣ), ",   ", (@sprintf "rₛ=%14.8e" status.rₛ))

function print_solver_status(status::NonlinearSolverStatus, params::NonlinearSolverParameters, n::Int)
    if (get_config(:verbosity) == 1 && !(check_solver_converged(status, params) && status.i ≤ params.nmax)) ||
        get_config(:verbosity) > 1
        println((@sprintf "  i=%07i" n), ",", status)
    end
end

function print_solver_status(status::NonlinearSolverStatus, params::NonlinearSolverParameters)
    if (get_config(:verbosity) == 1 && !(check_solver_converged(status, params) && status.i ≤ params.nmax)) ||
        get_config(:verbosity) > 1
        println(status)
    end
end

function check_solver_converged(status::NonlinearSolverStatus, params::NonlinearSolverParameters)
    return status.rₐ ≤ params.atol  ||
           status.rᵣ ≤ params.rtol  ||
           status.rₛ ≤ params.stol
        #    all(x -> x==0, status.Δx) ||
        #    maximum(status.Δx) / maximum(abs.(status.xₚ)) ≤ params.stol
end

function check_solver_status(status::NonlinearSolverStatus, params::NonlinearSolverParameters, n::Int)
    if any(x -> isnan(x), status.xₚ)
        error("Detected NaN in it=", n)
    end

    if status.rₐ > params.atol_break
        error("Absolute error of nonlinear solver ($(status.rₐ)) larger than allowed ($(params.atol_break)) in it=", n)
    end

    if status.rᵣ > params.rtol_break
        error("Relative error of nonlinear solver ($(status.rᵣ)) larger than allowed ($(params.rtol_break)) in it=", n)
    end

    if status.rₛ > params.stol_break
        error("Succesive error of nonlinear solver ($(status.rₛ)) larger than allowed ($(params.stol_break)) in it=", n)
    end
end

function getLinearSolver(T, n)
    linear_solver = get_config(:ls_solver)

    if linear_solver == nothing || linear_solver == :lapack
    # if linear_solver == :lapack
        linear_solver = LUSolverLAPACK{T}(BlasInt(n))
    elseif linear_solver == :julia
        linear_solver = LUSolver{T}(n)
    # elseif linear_solver == nothing
    #     linear_solver = IterativeSolver{T}(n)
    else
        @assert typeof(linear_solver) <: LinearSolver{T}
        @assert n == linear_solver.n
    end
    return linear_solver
end

function residual_initial!(status::NonlinearSolverStatus{T}, x::Vector{T}, y::Vector{T}) where {T}
    @assert length(x) == length(y) == length(status.r₀)

    status.r₀ .= y.^2
    status.rₐ  = maximum(status.r₀)
    status.rᵣ  = 1
    status.x₀ .= x
    status.xₚ .= x
    status.y₀ .= y
    status.yₚ .= y
end

function residual!(status::NonlinearSolverStatus{T}, δx::Vector{T}, x::Vector{T}, y::Vector{T}) where {T}
    residual_absolute!(status, y)
    residual_relative!(status, y)
    residual_successive!(status, δx, x)
    status.xₚ .= x
    status.yₚ .= y
end

function residual_absolute!(status::NonlinearSolverStatus{T}, y::Vector{T}) where {T}
    @assert length(y) == length(status.y₀)
    local rₐ²::T = 0
    @inbounds for yᵢ in y
        rₐ² = max(rₐ², yᵢ^2)
    end
    status.rₐ = sqrt(rₐ²)
end

function residual_relative!(status::NonlinearSolverStatus{T}, y::Vector{T}) where {T}
    @assert length(y) == length(status.y₀)
    status.rᵣ = norm(y) / norm(status.y₀)
end

function residual_successive!(status::NonlinearSolverStatus{T}, δx::Vector{T}, x::Vector{T}) where {T}
    @assert length(δx) == length(x)
    status.rₛ = norm(δx) / norm(x)
end


text:
"""
Convert an MPCCModel to an NLPModels as follows.
Definit le type NLMPCC :
min 	f(x)
s.t. 	l <= x <= u
	lcon(tb) <= cnl(x) <= ucon

with

cnl(x) := c(x),G(x),H(x),G(x).*H(x)
"""
mutable struct NLMPCC <: AbstractNLPModel

 meta     :: NLPModelMeta
 counters :: Counters

 x0       :: Vector

 mod      :: AbstractMPCCModel

 function NLMPCC(mod     :: AbstractMPCCModel;
                  x      :: Vector = mod.meta.x0,
				  name::String = "Generic",
				  lin::AbstractVector{<: Integer}=Int[])

  nvar, ncc = mod.meta.nvar, mod.meta.ncc
  new_lcon = vcat(mod.meta.lcon,
                  mod.meta.lccG,
                  mod.meta.lccH,
                  -Inf*ones(ncc))
  new_ucon = vcat(mod.meta.ucon,
                  Inf*ones(2*ncc),
                  zeros(ncc))
  ncon = maximum([length(new_lcon); length(new_ucon)])
  y0 = vcat(mod.meta.y0, zeros(3*ncc))

  nnzh = nvar * (nvar + 1) / 2
  nnzj = nvar * ncon
  nln = setdiff(1:ncon, lin)

  meta = NLPModelMeta(nvar, x0 = x, lvar = mod.meta.lvar, uvar = mod.meta.uvar,
                                    ncon = mod.meta.ncon+3*ncc, y0=y0,
                                    lcon = new_lcon, ucon = new_ucon,
								    nnzj=nnzj, nnzh=nnzh, lin=lin, nln=nln,
								    minimize=true, islp=false, name=name)

  return new(meta, Counters(), x, mod)
 end
end

function obj(nlp :: NLMPCC, x ::  AbstractVector)
 increment!(nlp, :neval_obj)
 return obj(nlp.mod, x)
end

function grad!(nlp :: NLMPCC, x :: Vector, gx :: AbstractVector)
 increment!(nlp, :neval_grad)
 return grad!(nlp.mod, x, gx)
end

function cons!(nlp :: NLMPCC, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_cons)
  Gx, Hx = consG(nlp.mod, x), consH(nlp.mod, x)
  c[1:nlp.meta.ncon] .= vcat(cons_nl(nlp.mod, x), Gx, Hx, Gx .* Hx)
  return c
end

function jac(nlp :: NLMPCC, x :: AbstractVector)
 increment!(nlp, :neval_jac)
 JGx, JHx = jacG(nlp.mod,x),   jacH(nlp.mod,x)
 Gx,  Hx  = consG(nlp.mod, x), consH(nlp.mod, x)
 A = vcat(jac_nl(nlp.mod,x), JGx, JHx, diagm(0 => Hx) * JGx + diagm(0 => Gx) * JHx)
 return A
end

function jac_coord!(nlp :: NLMPCC, x :: AbstractVector, vals ::AbstractVector)
 Jx = jac(nlp, x) #findnz(jac(nlp.mod, x))
 m, n = size(Jx)
 vals[1 : n*m] .= Jx[:]
 return vals
end

function jac_structure!(nlp :: NLMPCC, rows :: AbstractVector{<:Integer}, cols :: AbstractVector{<:Integer})
 m, n = nlp.meta.ncon, nlp.meta.nvar
 I = ((i,j) for i = 1:m, j = 1:n)
 rows[1 : n*m] .= getindex.(I, 1)[:]
 cols[1 : n*m] .= getindex.(I, 2)[:]
 return rows, cols
end

function jprod!(nlp :: NLMPCC, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  increment!(nlp, :neval_jprod)
  Jv .= jac(nlp, x) * v
  return Jv
end

function jtprod!(nlp :: NLMPCC, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  increment!(nlp, :neval_jtprod)
  Jtv .= jac(nlp,x)' * v
  return Jtv
end

function hess(nlp :: NLMPCC, x :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp.mod, x, obj_weight = obj_weight)
  return tril(Hx)
end

function hess(nlp :: NLMPCC, x :: AbstractVector, y :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  ncc, nlc, nvar = nlp.mod.meta.ncc, nlp.mod.meta.ncon, nlp.mod.meta.nvar
  JGx, JHx = jacG(nlp.mod,x),   jacH(nlp.mod,x)
  Gx,  Hx  = consG(nlp.mod, x), consH(nlp.mod, x)
  lphi = y[nlc+2*ncc+1:nlc+3*ncc]
  ny = vcat(y[1:nlc], y[nlc+1:nlc+ncc] - lphi .* Hx, y[nlc+ncc+1:nlc+2*ncc] - lphi .* Gx)
  Hx = hess(nlp.mod, x, ny, obj_weight = obj_weight)
  for i=1:ncc
   Hx += lphi[i] * (JGx[i,:] * JHx[i,:]' + JHx[i,:] * JGx[i,:]')
  end
  return tril(Hx)
end

function hess_structure!(nlp :: NLMPCC, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows .= getindex.(I, 1)
  cols .= getindex.(I, 2)
  return rows, cols
end

function hess_coord!(nlp :: NLMPCC, x :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hess_coord!(nlp :: NLMPCC, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, y, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hprod!(nlp :: NLMPCC, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, obj_weight=obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end

function hprod!(nlp :: NLMPCC, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, y, obj_weight=obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end

"""
Return the violation of the constraints
lb <= x <= ub,
lc <= c(x) <= uc,
lccG <= G(x),
lccH <= H(x),
G(x).*H(x) <= 0.
"""
function viol(nlp :: NLMPCC, x :: AbstractVector)

 mod, n, ncc = nlp.mod, nlp.mod.meta.nvar, nlp.mod.meta.ncc
 x = length(x) == n ? x : x[1:n]

 feas_x = vcat(max.(mod.meta.lvar-x, 0), max.(x-mod.meta.uvar, 0))

 if mod.meta.ncon !=0

  c = cons_nl(mod, x)
  feas_c = vcat(max.(mod.meta.lcon-c, 0), max.(c-mod.meta.ucon, 0))

 else

  feas_c = Float64[]

 end

 if ncc != 0

  G=consG(mod, x)
  H=consH(mod, x)

  feas_cp = vcat(max.(mod.meta.lccG-G, 0), max.(mod.meta.lccH-H, 0))
  feas_cc = max.(G.*H, 0)
 else
  feas_cp = Float64[]
  feas_cc = Float64[]
 end
 return vcat(feas_x, feas_c, feas_cp, feas_cc)
end


text:
### libs
using PyCall
using HDF5
using PyPlot, Colors
using Statistics
###
# path
file = "path\\to\\file.h5"
filesavepath = "path\\to\\save\\file\\"
# loading raw data
avgSpectrum = HDF5.h5read(file, "SPECdata/AverageSpec")
timebins = Array{Float64}(undef,length(avgSpectrum),1)
timebins = range(1,stop=length(avgSpectrum),length=length(avgSpectrum))
rawMassCalParams = HDF5.h5read(file, "CALdata/Spectrum")[1:2,2]
rawMassAxis = Array{Float64}(undef,length(avgSpectrum),1)
rawMassAxis = ((timebins .- rawMassCalParams[2])./rawMassCalParams[1]).^2
### plotting raw data
fig = figure(figsize=(12,8));
ax1 = fig.add_subplot(2,1,1);
semilogy(rawMassAxis,avgSpectrum)
xlabel("m/z [Th]")
ylabel("intensity [a.u.]")
title("To get the timebin interval of the desired peaks for mass axis calibration:\n click left and right of each peak in the lower panel.")
ax2 = fig.add_subplot(2,1,2);
semilogy(timebins,avgSpectrum)
xlabel("time bin [#]")
ylabel("intensity [a.u.]")
#savefig("$filesavepath plot_file.png", papertype="A4",orientation="landscape")
### processing data
using LsqFit
@. model(x, p) = p[1]*exp((x-p[2])^2/p[3]) # Gaussian distribution
#@. model(x, p) = p[1]/(p[1]^2 + (x-p[2])^2 ) # Lorentzian distribution, p[1]>0, -inf<p[2]<inf

# select TIMEBIN intervals for mass axis calibration
# -> click left and right of each peak to get the interval
nSteps = length(massCalCompounds) # number of peaks for mass axis calibration
gdata = Array{Float32}(undef,(2*nSteps),1)
for i =1:2:2*nSteps
    println("please zoom to the correct position of the masscal compound. When happy with the zoom stage, click a random key")
    readline()
    println("now click in the plot to set the peak boundaries")
    gin = ginput()
    gdata[i] = gin[1][1]
    ax2.plot([gdata[i], gdata[i]], [minimum(avgSpectrum), maximum(avgSpectrum)],"r-.")
    gin = ginput()
    gdata[i+1] = gin[1][1]
    ax2.plot([gdata[i+1], gdata[i+1]], [minimum(avgSpectrum), maximum(avgSpectrum)],"r-.")
end
# find interval indices
maxPeakVal = Any[]
IdxT = Array{Float64}(undef,length(avgSpectrum),1)
for l=1:1
    for i = 1:2:2*nSteps
        specInt = Any[]
        idxSpecInt = Any[]
        for m = 1:length(timebins)
            if gdata[i] < timebins[m] && timebins[m] < gdata[i+1]
                push!(specInt,avgSpectrum[m])
                push!(idxSpecInt,timebins[m])
            end
        end
        p0 = [maximum(specInt), idxSpecInt[findall(x->x==maximum(specInt),specInt)[1]], -maximum(specInt)*1000] # p0 for Guassian distribution
        #p0 = [maximum(specInt)*7, mean(idxSpecInt)] # p0 for Lorentzian distribution
        #lb = [maximum(specInt)-0.1*maximum(specInt), mean(idxSpecInt)-0.01*mean(idxSpecInt),-maximum(specInt)*10000]
        #ub = [maximum(specInt)+0.1*maximum(specInt), mean(idxSpecInt)+0.01*mean(idxSpecInt),Inf]
        fit = LsqFit.curve_fit(model, idxSpecInt, specInt, p0)#, lower=lb, upper=ub )
        push!(maxPeakVal, fit.param[2])
        figure(figsize=(10,7));
        plot(idxSpecInt,specInt)
        plot(idxSpecInt,model(idxSpecInt,fit.param) )
    end
    #return specInt, idxSpecInt, fit, maxPeakVal
end
# 3-params
@. model_tb2m(x, p) = ((x-p[2])/p[1])^(1/p[3])
p0_tb2m = [15000.0, -60000.0, 0.5] # [minimum(maxPeakVal)/1.4,maximum(maxPeakVal)]
fit_tb2m = LsqFit.curve_fit(model_tb2m, maxPeakVal, massCalCompounds, p0_tb2m ) # , p0_bounds, lower=lb, upper=ub)
tb2m = ((timebins.-fit_tb2m.param[2]) ./fit_tb2m.param[1]).^(1/fit_tb2m.param[3])
### plot result
figure(figsize=(12,8));
semilogy(rawMassAxis,avgSpectrum, label="old mass calibration")
semilogy(tb2m,avgSpectrum, label="new mass calibration")
xlabel("m/z [Th]")
ylabel("intensity [a.u.]")
legend()

### save data to APi HDF5 file
processedfilepath = file
fid = HDF5.h5open(processedfilepath, "cw")
HDF5.h5write(processedfilepath, "fitParams", fit_tb2m.param)
HDF5.close(fid)
### test if writing was succesful
fitiPara = HDF5.h5read(processedfilepath, "/fitParams")


text:
module TestBase

using Test
using ContextManagers.++

function test_iobuffer()
    @with(io = IOBuffer()) do
        @test io isa IOBuffer
    end
    @test true
end

function test_lock()
    l = ReentrantLock()
    @with(l) do
        @test l isa ReentrantLock
    end
    @sync begin
        @async lock(l) do
        end
    end
    @test true
end

function test_channel()
    ch = Channel(1)
    @with(ch) do
        put!(ch, 1)
    end
    @test !isopen(ch)
end

function test_opentemp_success()
    tmppath = nothing
    @with((path, io) = opentemp()) do
        tmppath = path
        @test io isa IO
        @test isfile(tmppath)
    end
    @test !isfile(tmppath)
end

function test_opentemp_error()
    tmppath = Ref{Any}(nothing)
    @with(IgnoreError(), (path, io) = opentemp()) do
        tmppath[] = path
        @test io isa IO
        @test isfile(tmppath[])
        error("error")
    end
    @test !isfile(tmppath[])
end

function test_opentempdir_success()
    tmppath = nothing
    @with(path = opentempdir()) do
        tmppath = path
        @test isdir(tmppath)
    end
    @test !isdir(tmppath)
end

function test_opentempdir_error()
    tmppath = Ref{Any}(nothing)
    @with(IgnoreError(), path = opentempdir()) do
        tmppath[] = path
        @test isdir(tmppath[])
        error("error")
    end
    @test !isdir(tmppath[])
end

end  # module


text:
# ---
# title: Symmetric monogenic filters
# id: demo_phasesymmono
# cover: assets/phasesymmono.gif
# author: Peter Kovesi
# date: 2018-10-26
# ---

# Phase symmetry responds well to line like features and circular objects.  The number of
# filter scales will affect the scale of features that are marked. Phase symmetry marks
# features independently of contrast (a bright circle is not more symmetric than a grey
# circle) and is a dimensionless quantity between 0 and 1.  However this may not be what one
# desires in which case the symmetry energy may be of greater interest.

using TestImages
using Images
using ImagePhaseCongruency

img = Gray.(testimage("blobs"))
## Detect regions of bright symmetry (polarity = 1)
phase_bright, = phasesymmono(img; nscale=5, polarity=1)

## Detect regions of dark symmetry (polarity = -1)
phase_dark, = phasesymmono(img; nscale=5, polarity=-1)

mosaic(img, phase_bright, phase_dark; nrow=1)

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "phasesymmono.gif"), Images.gif([phase_bright, phase_dark]); fps=1) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
# import libraries
using GeometricIntegrators
using ProgressMeter

# import ParticleMethods package
using ParticleMethods

# import HDF5 package for storing solution
using HDF5

# parameters
Δt = 0.1    # time step size
nt = 200    # number of time steps
np = 10000  # number of particles
nx = 16     # number of grid points

# output file
h5file = "vlasov_poisson.hdf5"

# random initial conditions
x₀ = randn(np)
v₀ = randn(np)

# shift x₀ to the interval [0,1]
xmax = ceil(maximum(abs.(x₀)))
x₀ .+= xmax
x₀ ./= 2*xmax

# concatenate x₀ and v₀
z₀ = collect(hcat(x₀, v₀)')

# vector field
function lorentz_force!(t, z, ż, p::PoissonSolver)
    Particles.solve!(p, z[1,:])
    for i in axes(ż, 2)
        ż[1,i] = z[2,i]
        ż[2,i] = eval_field(p, z[1,i])
    end
end

# splitting fields
function v_advection!(t, z, ż)
    for i in axes(ż, 2)
        ż[1,i] = z[2,i]
        ż[2,i] = 0
    end
end

function s_advection!(t, z, s, h)
    for i in axes(s, 2)
        s[1,i] = z[1,i] + h * z[2,i]
        s[2,i] = z[2,i]
    end
end

function v_lorentz_force!(t, z, ż, p::PoissonSolver)
    Particles.solve!(p, z[1,:])
    for i in axes(ż, 2)
        ż[1,i] = 0
        ż[2,i] = eval_field(p, z[1,i])
    end
end

function s_lorentz_force!(t, z, s, h, p::PoissonSolver)
    Particles.solve!(p, z[1,:])
    for i in axes(s, 2)
        s[1,i] = z[1,i]
        s[2,i] = z[2,i] + h * eval_field(p, z[1,i])
    end
end

# solution storage
function copy_to_hdf5(h5z, z, n)
    h5z[:,:,n+1] = z
end

# create Poisson solver
p = PoissonSolverFFT{eltype(z₀)}(nx)
Particles.solve!(p, x₀)

# create an ODE instance
# ode = ODE((t, z, ż) -> lorentz_force!(t, z, ż, p), z₀)

# create a splitting ODE instance
sode = SODE((v_advection!, (t, z, ż) -> v_lorentz_force!(t, z, ż, p)),
            (s_advection!, (t, z, s, h) -> s_lorentz_force!(t, z, s, h, p)),
            z₀)

# create HDF5 file and copy initial conditions
h5  = h5open(h5file, "w")
h5z = create_dataset(h5, "z", eltype(z₀), ((2, np, nt+1), (2, np, -1)), chunk=(2,np,1))
copy_to_hdf5(h5z, z₀, 0)

# create integrator
# int = IntegratorExplicitEuler(ode, Δt)
int = Integrator(sode, TableauStrang(), Δt)

# create atomic solution
# asol = AtomicSolution(ode)
asol = AtomicSolution(sode)

# copy initial conditons to atomic solution
# set_initial_conditions!(asol, ode)
set_initial_conditions!(asol, sode)

# initilize integrator
initialize!(int, asol)

# loop over time steps showing progress bar
try
    @showprogress 5 for n in 1:nt
        integrate_step!(int, asol)
        copy_to_hdf5(h5z, asol.q, n)
    end
finally
    # close HDF5 file
    close(h5)
end

# load Plots package
using Plots

# read array from HDF5 file
z = h5read(h5file, "z")

# compute plot ranges
vmax = ceil(maximum(abs.(v₀)))
xlim = (0, 1)
vlim = (-vmax, +vmax)

# plot initial condition
scatter(mod.(z[1,:,begin], 1), z[2,:,begin],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
)

# save figure to file
savefig("vlasov_poisson_z₀.png")

# create animation
anim = @animate for n in axes(z,3)
    scatter(mod.(z[1,:,n], 1), z[2,:,n],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
    )
end

# save animation to file
gif(anim, "vlasov_poisson_anim.gif", fps=10)


text:
using Documenter, JuliaInterpreter, Test, CodeTracking

DocMeta.setdocmeta!(JuliaInterpreter, :DocTestSetup, :(
    begin
        using JuliaInterpreter
        JuliaInterpreter.clear_caches()
        JuliaInterpreter.remove()
    end); recursive=true)

makedocs(
    modules = [JuliaInterpreter],
    clean = false,
    strict = true,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "JuliaInterpreter.jl",
    authors = "Keno Fischer, Tim Holy, and others",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "ast.md",
        "internals.md",
        "dev_reference.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaDebug/JuliaInterpreter.jl.git",
)


text:
module COVID19
# Usage:
#import COVID19
#(covid19g, covid19us) = COVID19.load_data();
#
# This code loads COVID19 data from https://github.com/CSSEGISandData/COVID-19
#
# COVID19:
# https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_covid_19_time_series
#
# FIPS:
# https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/UID_ISO_FIPS_LookUp_Table.csv

using CSV, DataFrames, Dates, Downloads

export load_data_from_web, load_data

mutable struct COVID19data
    country::AbstractString
    province_state::AbstractString
    county::AbstractString
    key::AbstractString
    fips::Int
    lat::AbstractFloat
    lon::AbstractFloat
    population::Int
    time::Array{Date}
    confirmed::Array{Int}
    death::Array{Int}
    #recovered::Array{Int}
end

mutable struct COVID19nyt
    time::Array{Date}
    state::Array{AbstractString}
    county::Array{AbstractString}
    fips::Array{Int}
    cases::Array{Int}
    death::Array{Int}
end

# define data path for non-windows and windows systems
if Sys.iswindows() == false
    datadir0 = ENV["HOME"]
    datadir1 = datadir0 * "/Box/Data/COVID19/"
else
    datadir0 = ENV["HOMEDRIVE"] * ENV["HOMEPATH"]
    datadir1 = datadir0 * "\\Box\\Data\\COVID19\\"
end

function load_data_from_nyt(datadir = datadir1)
    # data source: https://github.com/nytimes/covid-19-data
    covid19nyt_url_county = "https://github.com/nytimes/covid-19-data/raw/master/us-counties.csv";
    covid19nyt_url_state = "https://github.com/nytimes/covid-19-data/raw/master/us-states.csv";

    covid19nyt_path_county = datadir * "covid19nyt_county.csv";
    covid19nyt_path_state = datadir * "covid19nyt_state.csv";

    Downloads.download(covid19nyt_url_county, covid19nyt_path_county);
    Downloads.download(covid19nyt_url_state, covid19nyt_path_state);

    covid19nyt_county = DataFrame(CSV.File(covid19nyt_path_county));
    covid19nyt_state = DataFrame(CSV.File(covid19nyt_path_state));
    return covid19nyt_state, covid19nyt_county;
end

function load_data_from_web(datadir = datadir1)
    covid19url_confirmed_US = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
    covid19url_confirmed_global = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
    covid19url_death_US = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
    covid19url_death_global = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"
    covid19url_recovered_global = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv"
    covid19url_fips = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/UID_ISO_FIPS_LookUp_Table.csv"

    covid19path_confirmed_US = datadir * "covid19_confirmed_US.csv"
    covid19path_confirmed_global = datadir * "covid19_confirmed_global.csv"
    covid19path_death_US = datadir * "covid19_death_US.csv"
    covid19path_death_global = datadir * "covid19_death_global.csv"
    covid19path_recovered_global = datadir * "covid19_recovered_global.csv"
    covid19path_fips = datadir * "UID_ISO_FIPS_LookUp_Table.csv"

    Downloads.download(covid19url_confirmed_US, covid19path_confirmed_US)
    Downloads.download(covid19url_confirmed_global, covid19path_confirmed_global)
    Downloads.download(covid19url_death_US, covid19path_death_US)
    Downloads.download(covid19url_death_global, covid19path_death_global)
    Downloads.download(covid19url_recovered_global, covid19path_recovered_global)
    Downloads.download(covid19url_fips, covid19path_fips)

#=     io = open(covid19path_confirmed, "w")
    r = HTTP.request("GET", covid19url_confirmed, response_stream=io)
    close(io)
    io = open(covid19path_death, "w")
    r = HTTP.request("GET", covid19url_death, response_stream=io)
    close(io)
    io = open(covid19path_recovered, "w")
    r = HTTP.request("GET", covid19url_recovered, response_stream=io)
    close(io) =#

    confirmed_US = DataFrame(CSV.File(covid19path_confirmed_US));
    death_US = DataFrame(CSV.File(covid19path_death_US));
    confirmed_global = DataFrame(CSV.File(covid19path_confirmed_global));
    death_global = DataFrame(CSV.File(covid19path_death_global));
    recovered_global = DataFrame(CSV.File(covid19path_recovered_global));
    fipstable = DataFrame(CSV.File(covid19path_fips));

    return confirmed_US, death_US, confirmed_global, death_global, recovered_global, fipstable
end

#function load_covid19(confirmed_US::Any,death_US::Any,confirmed_global::Any,death_global::Any,fipstable::Any)
function load_data(datadir = datadir1)
    #(confirmed_US, death_US, confirmed_global, death_global, recovered_global, ipstable) = load_data_from_web(datadir);
    #(covid19nyt_state, covid19nyt_county) = load_data_from_nyt(datadir);

    covid19path_confirmed_US = datadir * "covid19_confirmed_US.csv"
    covid19path_confirmed_global = datadir * "covid19_confirmed_global.csv"
    covid19path_death_US = datadir * "covid19_death_US.csv"
    covid19path_death_global = datadir * "covid19_death_global.csv"
    covid19path_recovered_global = datadir * "covid19_recovered_global.csv"
    covid19path_fips = datadir * "UID_ISO_FIPS_LookUp_Table.csv"

    covid19nyt_path_county = datadir * "covid19nyt_county.csv";
    covid19nyt_path_state = datadir * "covid19nyt_state.csv";

    confirmed_US = DataFrame(CSV.File(covid19path_confirmed_US));
    death_US = DataFrame(CSV.File(covid19path_death_US));
    confirmed_global = DataFrame(CSV.File(covid19path_confirmed_global));
    death_global = DataFrame(CSV.File(covid19path_death_global));
    recovered_global = DataFrame(CSV.File(covid19path_recovered_global));
    fipstable = DataFrame(CSV.File(covid19path_fips));

    covid19nyt_county = DataFrame(CSV.File(covid19nyt_path_county));
    covid19nyt_state = DataFrame(CSV.File(covid19nyt_path_state));

    # defining the time
    t0 = Date("2020-01-22")
    ndays = Dates.today() - t0;
    t = collect(t0 : Day(1) : Dates.today());

    # specifying the first column in the data where the time series start
    c0_global = 5;
    c0_US = 12;

    # loading global COVID19 time series data
    latg = collect(Missings.replace(confirmed_global[:,3],NaN));
    long = collect(Missings.replace(confirmed_global[:,4],NaN));
    country = confirmed_global[:,2];
    province_state = collect(Missings.replace(confirmed_global[:,1],""));

    covid19global = COVID19data[];
    for i = 1:size(latg)[1]
        covid19global = push!(covid19global, COVID19data(country[i], province_state[i], "", province_state[i] * " ," * country[i], -9999, latg[i], long[i], -9999, t, collect(confirmed_global[i,5:end]), collect(death_global[i,5:end])));
    end

    # loading the US COVID19 time series
    latus = collect(Missings.replace(confirmed_US[:,9],NaN));
    lonus = collect(Missings.replace(confirmed_US[:,10],NaN));
    countryUS = confirmed_US[:,8];
    province_stateUS = collect(Missings.replace(confirmed_US[:,7],""));
    fips = collect(Missings.replace(confirmed_US[:,5],-99999));
    county = collect(Missings.replace(confirmed_US[:,6],""));
    key = confirmed_US[:,11];
    population = collect(Missings.replace(death_US[:,12],0));

    covid19us = COVID19data[];
    for i = 1:size(latus)[1]
        #display(i)
        covid19us = push!(covid19us, COVID19data(countryUS[i], province_stateUS[i], county[i], key[i], fips[i], latus[i], lonus[i], population[i], t, collect(confirmed_US[i,12:end]), collect(death_US[i,13:end])));
    end

    return covid19global, covid19us, covid19nyt_state, covid19nyt_county
end

end


text:
__precompile__()
module Map_removeEQ

include("get_kurtosis.jl")
include("remove_eq.jl")

using .Get_kurtosis, .Remove_eq
using SeisIO, Dates, JLD2, Printf, FileIO
using Suppressor

export map_removeEQ

"""
loaddata(finame, path)
"""
function loaddata(finame, path)
    try
        @suppress_err return FileIO.load(finame, path)
    catch
        return false;
    end
end

"""
    ParallelEQremoval(dlid, InputDict::Dict)
    remove earthquake and save it into jld2 file.

"""
function map_removeEQ(dlid, InputDict::Dict)

    #store data

    finame                 =InputDict["finame"]
    IsKurtosisRemoval      =InputDict["IsKurtosisRemoval"]
    max_edgetaper_duration =InputDict["max_edgetaper_duration"]
    kurtosis_tw_sparse     =InputDict["kurtosis_tw_sparse"]
    kurtosis_timewindow    =InputDict["kurtosis_timewindow"]
    kurtosis_threshold     =InputDict["kurtosis_threshold"]
    IsSTALTARemoval        =InputDict["IsSTALTARemoval"]
    stalta_longtimewindow  =InputDict["stalta_longtimewindow"]
    stalta_threshold       =InputDict["stalta_threshold"]
    stalta_absoluteclip    =InputDict["stalta_absoluteclip"]
    max_wintaper_duration  =InputDict["max_wintaper_duration"]
    removal_shorttimewindow=InputDict["removal_shorttimewindow"]
    overlap                =InputDict["overlap"]
    Iswhiten               =InputDict["Iswhiten"]
    freqmin_whiten         =InputDict["freqmin_whiten"]
    freqmax_whiten         =InputDict["freqmax_whiten"]
    plot_kurtosis_α        =InputDict["plot_kurtosis_α"]
    plot_boxheight         =InputDict["plot_boxheight"]
    plot_span              =InputDict["plot_span"]
    plot_fmt               =InputDict["plot_fmt"]
    fodir                  =InputDict["fodir"]
    foname                 =InputDict["foname"]
    fopath                 =InputDict["fopath"]
    IsSaveFig              =InputDict["IsSaveFig"]

    DLtimestamplist        =InputDict["DLtimestamplist"]
    stationlist            =InputDict["stationlist"]
    NumofTimestamp         =InputDict["NumofTimestamp"]

    tstamp = DLtimestamplist[dlid]

    if mod(dlid, round(0.1*NumofTimestamp)+1) == 0
        println(@sprintf("start process %s", tstamp))
    end

    bt_getkurtosis = 0.0
    bt_removeeq = 0.0

    for st = stationlist
        #S = t[joinpath(tstamp, st)]
        st1 = replace(st, "-"=>"")

        # load raw data
        Stemp = loaddata(finame, joinpath(tstamp, st1))

        # skip if it does not exist
        if Stemp == false
            continue;
        end

        #convert if it's SeisChannel
        if Stemp isa SeisIO.SeisChannel
            Stemp = SeisData(Stemp)
        end

        SremEQ = SeisData()

        #loop by SeisChannel and save into temporal wile
        for Seisid = 1:Stemp.n

            Schan = Stemp[Seisid]

            if !haskey(Schan.misc, "dlerror")
                if iszero(Schan.x)
                    Schan.misc["dlerror"] = true
                else
                    Schan.misc["dlerror"] = false
                end
            end

            #---Transient error on combination between download margin and resampling---#
            # S.t may contain index 0 in t[:, 1]. in that case , this causes the error on later part, so discard that data.
            if any(Schan.t[:,1] .== 0)
                println("data gap index found.")
                println(Schan)
                Schan.misc["dlerror"] = true
            end
            #---------------------------------------------------------------------------#

            if Schan.misc["dlerror"] == 0

                dt = 1/Schan.fs
                tvec = collect(0:length(Schan.x)-1) * dt ./ 60 ./ 60

                #tapering to avoid instrumental edge artifacts
                #SeisIO.taper!(Schan,  t_max = max_edgetaper_duration, α=0.05) this is already done by download margin

                S1 = deepcopy(Schan)

                #set long window length to user input since last window of previous channel will have been adjusted
                S1.misc["eqtimewindow"] = fill(true, length(S1.x))

                if IsKurtosisRemoval
                    # compute kurtosis and detect
                    bt_1 = @elapsed S1 = Get_kurtosis.get_kurtosis(S1, float(kurtosis_timewindow), float(kurtosis_tw_sparse))

                    bt_2 = @elapsed S1 = Remove_eq.detect_eq_kurtosis(S1, tw=float(removal_shorttimewindow), kurtosis_threshold=float(kurtosis_threshold), overlap=float(overlap))

                    btsta_1 = 0

                    if IsSTALTARemoval
                        # detect earthquake and tremors by STA/LTA
                        btsta_1 = @elapsed S1 = Remove_eq.detect_eq_stalta(S1, float(stalta_longtimewindow), float(removal_shorttimewindow),
                                            float(stalta_threshold), float(overlap), float(stalta_absoluteclip),  InputDict, tstamp, st1)
                    end

                    if Iswhiten
                        Remove_eq.s_whiten!(S1, freqmin_whiten, freqmax_whiten)# apply spectral whitening on this channel
                    end

                    bt_3 = @elapsed S1 = Remove_eq.remove_eq(S1, Schan, float(plot_kurtosis_α), float(max_wintaper_duration),
                                    plot_boxheight, trunc(Int, plot_span), plot_fmt, fodir, tstamp, tvec, IsSaveFig)

                    bt_getkurtosis += bt_1
                    bt_removeeq += bt_2 + bt_3 + btsta_1

                    #println([bt_2, bt_3, btsta_1])
                    #if mod(dlid, round(0.1*NumofTimestamp)+1) == 0
                    #    println([bt_1, bt_2, btsta_1, bt_3])
                    #end

                else
                    #only STA/LTA
                    if IsSTALTARemoval

                        data.misc["kurtosis"] = zeros(length(S1.x))

                        bt_2 = @elapsed S1 = Remove_eq.detect_eq_stalta(S1, float(stalta_longtimewindow), float(removal_shorttimewindow),
                                            float(stalta_threshold), float(overlap), InputDict, tstamp, st1)

                        if Iswhiten
                            Remove_eq.s_whiten!(S1, freqmin_whiten, freqmax_whiten) # apply spectral whitening on this channel
                        end

                        bt_3 = @elapsed S1 = Remove_eq.remove_eq(S1, S, float(plot_kurtosis_α), float(max_wintaper_duration),
                                        plot_boxheight, trunc(Int, plot_span), plot_fmt, fodir, tstamp, tvec, IsSaveFig)

                        bt_getkurtosis += 0.0
                        bt_removeeq += bt_2 + bt_3

                    else
                        #no removal process.
                        @warn "Both 'IsKurtosisRemoval' and 'IsKurtosisRemoval' are false. No removal process is executed. Abort."
                        exit(0)
                    end

                end

                if InputDict["IsOutputRemovalFrac"]
                    #output removal fraction on this channel
                    eqidlist = S1.misc["eqtimewindow"][:]
                    numofremoval = sum(x->x==false, eqidlist, dims=1)
                    fractionofremoval = numofremoval[1] / length(eqidlist)

                    y, jd = parse.(Int64, split(InputDict["DLtimestamplist"][dlid], ".")[1:2])
                    tstamp_fname = replace(tstamp, ":" => "-")
                    fname_out = join([tstamp_fname, st1, "removalfraction","dat"], '.')
                    fopath = InputDict["removal_fractionpath"]*"/"*fname_out
                    open(fopath, "w") do io
                       write(io, @sprintf("%f\n", fractionofremoval))
                    end
                end

                if InputDict["dumptraces"]
                    #dump raw trace
                    tstamp_fname = replace(tstamp, ":" => "-")
                    fname_out = join([tstamp_fname, st1, "rawdata","dat"], '.')
                    fopath = InputDict["dumppath"]*"/"*fname_out
                    open(fopath, "w") do io
                        for i in 1:length(Schan.x)
                            write(io, @sprintf("%12.8f\n", Schan.x[i]))
                        end
                    end

                    #dump kurtosis trace
                    mkpath(InputDict["dumppath"])
                    tstamp_fname = replace(tstamp, ":" => "-")
                    fname_out = join([tstamp_fname, st1, "kurtosis","dat"], '.')
                    fopath = InputDict["dumppath"]*"/"*fname_out
                    open(fopath, "w") do io
                        for i in 1:length(S1.misc["kurtosis"])
                            write(io, @sprintf("%12.8f\n", S1.misc["kurtosis"][i]))
                        end
                    end

                    #dump removed trace
                    tstamp_fname = replace(tstamp, ":" => "-")
                    fname_out = join([tstamp_fname, st1, "remdata","dat"], '.')
                    fopath = InputDict["dumppath"]*"/"*fname_out
                    open(fopath, "w") do io
                        for i in 1:length(S1.x)
                            write(io, @sprintf("%12.8f\n", S1.x[i]))
                        end
                    end
                end

                #it's not allowed to save this into binary;
                delete!(S1.misc, "eqtimewindow")
                delete!(S1.misc, "kurtosis")

            else
                #download error found: save as it is.
                S1 = Schan
            end

            SremEQ += S1

        end

        # if some of SeisChannels in Stemp have a data, save temp file
        y, jd = parse.(Int64, split(InputDict["DLtimestamplist"][dlid], ".")[1:2])
        tstamp_fname = replace(tstamp, ":" => "-")
        fname_out = join([tstamp_fname, st1, "FDSNWS","dat"], '.')
        # save as intermediate binary file
        t_write = @elapsed SeisIO.wseis(InputDict["tmppath"]*"/"*fname_out, SremEQ)

    end

    return nothing
end

end

max_stars_repo_name: kura-okubo/SeisRemoveEQ

max_issues_repo_name: kura-okubo/SeisRemoveEQ

max_forks_repo_name: kura-okubo/SeisRemoveEQ


text:
@doc raw"""
Penning trap with asymmetric magnetic field in (x,y,z) coordinates.
Based on Yanyan Shi, Yajuan Sun, Yulei Wang, Jian Liu, Study of adaptive symplectic methods for
      simulating charged particle dynamics, Journal of Computational Dynamics 6, 429-448, 2019.

The covariant components of the vector potential are given by
```math
A (x,y,z) = B_0 / 2 \, ( -y , x-z/6, y / 6)^T + B_1 / 2 \, (z^2 - y^2, z^2 - x^2, y^2 - x^2)^T,
```
resulting in the magnetic field with covariant components
```math
B (x,y,z) = B_0 \, ( 1/3, 0, 1)^T + B_1 \, ( y-z, x+z, y-x )^T ,
```
and the electrostatic potential given by
```math
\varphi (x,y,z) = E_0 \, ( x^2 / 2 + y^2 / 2 - z^2) ,
```
resulting in the electric field with covariant components
```math
E (x,y,z) = E_0 \, ( x, y, - 2 z)^T .
```

Parameters:
* `B₀`: B-field strength
* `Bₚ`: B-field perturbation strength
* `E₀`: E-field strength
"""
module PenningTrapAsymmetric

    using RecipesBase

    import ..ElectromagneticFields
    import ..ElectromagneticFields: CartesianEquilibrium, code
    import ..AnalyticCartesianField: X, Y, Z

    export  PenningTrapAsymmetricEquilibrium

    const DEFAULT_B₀ = 100.0
    const DEFAULT_Bₚ = 50.0
    const DEFAULT_E₀ = 10.0

    struct PenningTrapAsymmetricEquilibrium{T <: Number} <: CartesianEquilibrium
        name::String
        B₀::T
        Bₚ::T
        E₀::T

        function PenningTrapAsymmetricEquilibrium{T}(B₀::T, Bₚ::T, E₀::T) where T <: Number
            new("PenningTrapAsymmetricEquilibrium", B₀, Bₚ, E₀)
        end
    end

    PenningTrapAsymmetricEquilibrium(B₀::T=DEFAULT_B₀, Bₚ::T=DEFAULT_Bₚ, E₀::T=DEFAULT_E₀) where T <: Number = PenningTrapAsymmetricEquilibrium{T}(B₀, Bₚ, E₀)

    function init(B₀=DEFAULT_B₀, Bₚ=DEFAULT_Bₚ, E₀=DEFAULT_E₀)
        PenningTrapAsymmetricEquilibrium(B₀, E₀, Bₚ)
    end

    macro code(B₀=DEFAULT_B₀, Bₚ=DEFAULT_Bₚ, E₀=DEFAULT_E₀)
        code(init(B₀, Bₚ, E₀); escape=true)
    end

    function Base.show(io::IO, equ::PenningTrapAsymmetricEquilibrium)
        print(io, "Penning trap with asymmetric magnetic field in (x,y,z) coordinates with\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  Bₚ = ", equ.Bₚ, "\n")
        print(io, "  E₀ = ", equ.E₀)
    end

    ElectromagneticFields.A₁(x::AbstractVector, equ::PenningTrapAsymmetricEquilibrium) = - equ.B₀ / 2 *  Y(x,equ)                  + equ.Bₚ / 2* (Z(x,equ)^2 - Y(x,equ)^2)
    ElectromagneticFields.A₂(x::AbstractVector, equ::PenningTrapAsymmetricEquilibrium) = + equ.B₀ / 2 * (X(x,equ) - Z(x,equ) / 6 ) + equ.Bₚ / 2* (Z(x,equ)^2 - X(x,equ)^2)
    ElectromagneticFields.A₃(x::AbstractVector, equ::PenningTrapAsymmetricEquilibrium) = + equ.B₀ / 2 *             Y(x,equ) / 6   + equ.Bₚ / 2* (Y(x,equ)^2 - X(x,equ)^2)
    ElectromagneticFields.φ(x::AbstractVector, equ::PenningTrapAsymmetricEquilibrium) = - equ.E₀ * (X(x,equ)^2 / 2 + Y(x,equ)^2 / 2 - Z(x,equ)^2)

    ElectromagneticFields.get_functions(::PenningTrapAsymmetricEquilibrium) = (X=X, Y=Y, Z=Z)

end


text:
type QPTrace
    x
    xc
    value
    search
    clamped
    nfactor
end

"""
 Minimize `0.5*x'*H*x + x'*g`  s.t. lower<=x<=upper

  inputs:
     `H`            - positive definite matrix   (n * n)
     `g`            - bias vector                (n)
     `lower`        - lower bounds               (n)
     `upper`        - upper bounds               (n)

   optional inputs:
     `x0`           - initial state              (n)
     `options`      - see below                  (7)

  outputs:
     `x`            - solution                   (n)
     `result`       - result type (roughly, higher is better, see below)
     `Hfree`        - subspace cholesky factor   (n_free * n_free)
     `free`         - set of free dimensions     (n)
"""
function boxQP(H,g,lower,upper,x0;
               maxIter        = 100,
               minGrad        = 1e-8,
               minRelImprove  = 1e-8,
               stepDec        = 0.6,
               minStep        = 1e-22,
               Armijo         = 0.1,
               print          = 0)
    #     maxIter        = 100        # maximum number of iterations
    #     minGrad        = 1e-8       # minimum norm of non-fixed gradient
    #     minRelImprove  = 1e-8       # minimum relative improvement
    #     stepDec        = 0.6        # factor for decreasing stepsize
    #     minStep        = 1e-22      # minimal stepsize for linesearch
    #     Armijo         = 0.1    	# Armijo parameter (fraction of linear improvement required)
    #     print          = 0 			# verbosity

    n        = size(H,1)
    clamped  = falses(n)
    free     = trues(n)
    oldvalue = 0.
    result   = 0
    gnorm    = 0.
    nfactor  = 0
    trace    = Array(QPTrace,maxIter)
    Hfree    = zeros(n,n)
    clamp(x) = max(lower, min(upper, x))

    debug("# initial state")
    x = clamp(x0[:])
    LU = [lower upper]
    LU[!isfinite(LU)] = NaN

    debug("# initial objective value")
    value    = (x'*g + 0.5*x'*H*x )[1]

    if print > 0
        @printf("==========\nStarting box-QP, dimension %-3d, initial value: %-12.3f\n",n, value)
    end

    debug("# main loop")
    iter = 1
    while iter <= maxIter

        if result != 0
            break
        end

        debug("# check relative improvement")
        if iter>1 && (oldvalue - value) < minRelImprove*abs(oldvalue)
            result = 4
            break
        end
        oldvalue = value

        debug("# get gradient")
        grad     = g + H*x

        debug("# find clamped dimensions")
        old_clamped                        = clamped
        clamped                            = falses(n)
        #         clamped[(x[:,1] .== lower)&(grad[:,1].>0)]   = true
        #         clamped[(x[:,1] .== upper)&(grad[:,1].<0)]   = true
        for i = 1:n
            clamped[i]   = ((x[i,1] == lower[i])&&(grad[i,1]>0)) || ((x[i,1] == upper[i])&&(grad[i,1]<0))
        end
        free                               = !clamped

        debug("# check for all clamped")
        if all(clamped)
            result = 6
            break
        end

        debug("# factorize if clamped has changed")
        if iter == 1
            factorize    = true
        else
            factorize    = any(old_clamped != clamped)
        end

        if factorize
            Hfree  = chol(H[free,free])  # was   (Hfree, indef)  = chol(H[free,free])
            #             if indef
            #                 result = -1
            #                 break
            #             end
            nfactor += 1
        end

        debug("# check gradient norm")
        gnorm  = norm(grad[free])
        if gnorm < minGrad
            result = 5
            break
        end

        debug("# get search direction")
        grad_clamped   = g  + H*(x.*clamped)
        search         = zeros(n,1)
        search[free]   = -Hfree\(Hfree'\grad_clamped[free]) - x[free]

        debug("# check for descent direction")
        sdotg          = sum(search.*grad)
        if sdotg >= 0 # (should not happen)
            break
        end

        debug("# armijo linesearch")
        step  = 1
        nstep = 0
        xc    = clamp(x+step*search)
        vc    = (xc'*g + 0.5*xc'*H*xc)[1]
        while (vc - oldvalue)/(step*sdotg) < Armijo
            step  = step*stepDec
            nstep += 1
            xc    = clamp(x+step*search)
            vc    = (xc'*g + 0.5*xc'*H*xc)[1]
            if step<minStep
                result = 2
                break
            end
        end

        if print > 1
            @printf("iter %-3d  value % -9.5g |g| %-9.3g  reduction %-9.3g  linesearch %g^%-2d  n_clamped %d\n",
                    iter, vc, gnorm, oldvalue-vc, stepDec, nstep, sum(clamped))
        end

        trace[iter] = QPTrace(x,xc,value,search,clamped,nfactor )

        debug("# accept candidate")
        x     = xc
        value = vc
        iter += 1

    end

    if iter == maxIter
        result = 1
    end

    results = ["Hessian is not positive definite",          # result = -1
               "No descent direction found",                # result = 0    SHOULD NOT OCCUR
               "Maximum main iterations exceeded",          # result = 1
               "Maximum line-search iterations exceeded",   # result = 2
               "No bounds, returning Newton point",         # result = 3
               "Improvement smaller than tolerance",        # result = 4
               "Gradient norm smaller than tolerance",      # result = 5
               "All dimensions are clamped"]                # result = 6

    if print > 0
        @printf("RESULT: %s.\niterations %d  gradient %-12.6g final value %-12.6g  factorizations %d\n",
                results[result+2], iter, gnorm, value, nfactor)
    end

    return x,result,Hfree,free,trace
end

function demoQP()
    n 		= 500
    g 		= randn(n,1)
    H 		= randn(n,n)
    H 		= H*H'
    lower 	= -ones(n,1)
    upper 	=  ones(n,1)
    @time boxQP(H, g, lower, upper, randn(n,1),print=1)

end


text:
#########################################################################
#
# GET Functions, No Required Arguments, Returns JSON Julia data structure
#
#########################################################################

funcname = (:get_help_configuration,
            :get_help_languages,
            :get_help_privacy,
            :get_help_tos,
            :get_application_rate_limit_status,
            :get_profile_banner,
            :get_blocks_ids,
            :get_account_settings,
            :get_oembed,
            :get_trends_available,
            :get_user_suggestions,
            :get_saved_searches_list,
            :get_lists_subscribers_show,
            :get_friendships_no_retweets,
            :get_friends_ids,
            :get_followers_ids,
            :get_friendships_lookup,
            :get_friendships_incoming,
            :get_friendships_outgoing,
            :get_friendships_show)

endpoint = ("help/configuration.json",
            "help/languages.json",
            "help/privacy.json",
            "help/tos.json",
            "application/rate_limit_status.json",
            "users/profile_banner.json",
            "blocks/ids.json",
            "account/settings.json",
            "statuses/oembed.json",
            "trends/available.json",
            "users/suggestions.json",
            "saved_searches/list.json",
            "lists/subscribers/show.json",
            "friendships/no_retweets/ids.json",
            "friends/ids.json",
            "followers/ids.json",
            "friendships/lookup.json",
            "friendships/incoming.json",
            "friendships/outgoing.json",
            "friendships/show.json")

for (func, endp) in zip(funcname, endpoint)
@eval begin function ($func)(; options=Dict{AbstractString, AbstractString}())

        r = get_oauth($"https://api.twitter.com/1.1/$endp", options)

        return r.status == 200 ? Requests.json(r) : r

        end
    end
end


text:
# Element Click
"""
    click!(element::Element)::Nothing

The Element Click command scrolls into view the element if it is not already pointer-interactable, and clicks its in-view center point.
If the element’s center point is obscured by another element, an element click intercepted error is returned. If the element is outside the viewport, an element not interactable error is returned.
"""
function click!(element::Element)::Nothing
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.post(
        "$addr/session/$id/element/$element_id/click",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
end


text:
module RaviartThomasRefFEsTest

using Test
using Gridap.Polynomials
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Fields: MockField
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData
using Gridap.Arrays
using Gridap.Visualization

using Gridap
using Gridap.FESpaces

p = QUAD
D = num_dims(QUAD)
et = Float64
order = 0

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test Conformity(reffe) == DivConformity()
p = QUAD
D = num_dims(QUAD)
et = Float64
order = 1

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 12
@test num_dofs(reffe) == 12
@test get_order(get_prebasis(reffe)) == 1

prebasis = get_prebasis(reffe)
dof_basis = get_dof_basis(reffe)

v = VectorValue(3.0,0.0)
field = GenericField(x->v*x[1])

cache = return_cache(dof_basis,field)
r = evaluate!(cache, dof_basis, field)
test_dof_array(dof_basis,field,r)

cache = return_cache(dof_basis,prebasis)
r = evaluate!(cache, dof_basis, prebasis)
test_dof_array(dof_basis,prebasis,r)

p = TET
D = num_dims(TET)
et = Float64
order = 0

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 4
@test num_dofs(reffe) == 4
@test get_order(get_prebasis(reffe)) == 0
@test Conformity(reffe) == DivConformity()

p = TET
D = num_dims(p)
et = Float64
order = 2

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 36
@test num_dofs(reffe) == 36
@test get_order(get_prebasis(reffe)) == 2
@test Conformity(reffe) == DivConformity()

prebasis = get_prebasis(reffe)
dof_basis = get_dof_basis(reffe)

v = VectorValue(0.0,3.0,0.0)
field = GenericField(x->v)

cache = return_cache(dof_basis,field)
r = evaluate!(cache, dof_basis, field)
test_dof_array(dof_basis,field,r)

cache = return_cache(dof_basis,prebasis)
r = evaluate!(cache, dof_basis, prebasis)
test_dof_array(dof_basis,prebasis,r)

# Factory function
reffe = ReferenceFE(QUAD,raviart_thomas,0)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test Conformity(reffe) == DivConformity()

reffe = ReferenceFE(QUAD,raviart_thomas,Float64,0)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test Conformity(reffe) == DivConformity()

@test Conformity(reffe,:L2) == L2Conformity()
@test Conformity(reffe,:Hdiv) == DivConformity()
@test Conformity(reffe,:HDiv) == DivConformity()

@test RaviartThomas() == raviart_thomas

end # module


text:
module SummarizedExperiments

using Reexport
@reexport using DataStructures: OrderedDict
@reexport using DataFrames: DataFrame, insertcols!

include("class.jl")
export SummarizedExperiment

include("subset.jl")
include("assign.jl")

include("getters.jl")
export size, rowdata, coldata, assay, assays, metadata

include("setters.jl")
export setrowdata!, setcoldata!, setassay!, setassays!, setmetadata!

include("combine.jl")
include("miscellaneous.jl")

include("example.jl")
export exampleobject

end


text:
using Recommenders
using Test, Documenter

const tests = [
    "dataset/movielens.jl",
    "dataset/data_split.jl",
    "dataset/data_utils.jl",
    "core/item_knn.jl",
    "core/loss_function.jl",
    "core/randomwalk_bipartite.jl",
    "metric.jl",
    "model/most_popular.jl",
    "model/item_knn.jl",
    "model/implicit_mf.jl",
    "model/bpr.jl",
    "model/slim.jl",
    "model/randomwalk.jl",
]

for t in tests
    @testset "Test $t" begin
        include(t)
    end
end


text:
useproc("CPU")

StartCode()
x=ADvariable() # global memory
y=x[:,1]+x[:,2]
loss=sum(y)
net=EndCode() # defines the graph

# parameters:
net.value[x]=randn(3,2)
net=compile(net,debug=true) # compile and preallocate memory
@gpu CUDArt.init([0])
@gpu net=convert(net,"GPU")

gradcheck(net;showgrad=true) # use a small number of datapoints and small network to check the gradient, otherwise

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Test
using SymDiff

@test subs(:a, (:a => 1)) == 1
@test subs(:a, (:b => 1)) == :a


text:
using Literate
using Printf

repo_src = joinpath(@__DIR__,"..","src")
notebooks_dir = joinpath(@__DIR__,"..","notebooks")

files = [
  "Example of a fixed point algorithm"=>"fixed-point.jl",
  "The overfitting problem (handle unrelated stopping criteria)"=>"overfitting.jl",
  "Checkpointing"=>"checkpointing.jl",
  "Mixed-algorithms, a ListofStates tutorial"=>"gradient-lbfgs.jl",
  "Benchmark solvers"=>"benchmark.jl",
]

Sys.rm(notebooks_dir;recursive=true,force=true)
for (i,(title,filename)) in enumerate(files)
  notebook_prefix = string("t",@sprintf "%03d_" i)
  notebook = string(notebook_prefix,splitext(filename)[1])
  notebook_title = string("# # Tutorial ", i, ": ", title)
  function preprocess_notebook(content)
    return string(notebook_title, "\n\n", content)
  end
  Literate.notebook(joinpath(repo_src,filename), notebooks_dir; name=notebook, preprocess=preprocess_notebook, documenter=false, execute=false)
end

deps_jl = joinpath(@__DIR__, "deps.jl")

if isfile(deps_jl)
  rm(deps_jl)
end

open(deps_jl,"w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const files = $files))
end


text:
using Documenter, ForTheBadge

makedocs(;
    modules=[ForTheBadge],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/asinghvi17/ForTheBadge.jl/blob/{commit}{path}#L{line}",
    sitename="ForTheBadge.jl",
    authors="Anshul Singhvi",
    assets=String[],
)

deploydocs(;
    repo="github.com/asinghvi17/ForTheBadge.jl",
)


text:
export matching_poly

"""
`matching_poly(G)` returns the matching polynomial of the
`SimpleGraph` `G`.
"""
function matching_poly(G::SimpleGraph, cache_flag::Bool = true)
    if cache_flag && cache_check(G, :matching_poly)
        return cache_recall(G, :matching_poly)
    end
    if NE(G) == 0
        x = getx()
        p = x^NV(G)
        # p = Polynomial([0,1])^NV(G)
        if cache_flag
            SimpleGraphs.cache_save(G, :matching_poly, p)
        end
        return p
    end

    if is_connected(G)
        e = first(G.E)  # get any edge
        v, w = e

        G1 = deepcopy(G)
        SimpleGraphs.delete!(G1, v, w)
        p1 = matching_poly(G1, false)

        SimpleGraphs.delete!(G1, v)
        SimpleGraphs.delete!(G1, w)
        p2 = matching_poly(G1, false)

        p = p1 - p2
        if cache_flag
            cache_save(G, :matching_poly, p)
        end
        return p
    end

    comps = parts(components(G))
    p = SimplePolynomial(1)
    for S in comps
        H = induce(G, S)
        pH = matching_poly(H, false)
        p *= pH
    end
    if cache_flag
        cache_save(G, :matching_poly, p)
    end
    return p
end


text:
function get_component(forecast::Forecast)
    return IS.get_component(forecast)
end

function get_forecast_component_name(forecast::Forecast)
    return IS.get_forecast_component_name(forecast)
end
function get_forecast_value(forecast::Forecast, ix)
    return IS.get_forecast_value(forecast, ix)
end
function get_horizon(forecast::Forecast)
    return IS.get_horizon(forecast)
end
function get_timeseries(forecast::Forecast)
    return IS.get_timeseries(forecast)
end

function get_data(forecast::Forecast)
    return IS.get_data(forecast)
end


text:
#! /usr/bin/env julia

using Base64: base64encode, base64decode
using BenchmarkTools
using Test
using CUDA

function getFrequencies()::Tuple{UInt32, UInt32, Vector{Float32}}
    letterScores  = [
    ['A',50.0778966131907],
    ['B',30.2092691622103],
    ['C',40.8846702317291],
    ['D',23.1073083778966],
    ['E',24.6778966131907],
    ['F',17.9591800356506],
    ['G',16.6153297682709],
    ['H',22.0377896613191],
    ['I',39.8060606060606],
    ['J',14.0295900178253],
    ['K',8.3030303030303],
    ['L',19.0702317290553],
    ['M',46.2520499108734],
    ['N',36.6147950089127],
    ['O',18.8413547237077],
    ['P',25.7110516934046],
    ['Q',2.07825311942959],
    ['R',26.1048128342246],
    ['S',54.3620320855615],
    ['T',58.0146167557932],
    ['U',10.2474153297683],
    ['V',5.53529411764706],
    ['W',19.1078431372549],
    ['X',1.35080213903743],
    ['Y',16.8087344028521],
    ['Z',1.0],
    ['a',938.285026737968],
    ['b',154.395008912656],
    ['c',349.449554367201],
    ['d',422.427807486631],
    ['e',1380.00748663102],
    ['f',231.180926916221],
    ['g',215.106417112299],
    ['h',526.890909090909],
    ['I',807.011051693405],
    ['j',11.7390374331551],
    ['k',82.1368983957219],
    ['l',455.107308377897],
    ['m',261.564349376114],
    ['n',808.47504456328],
    ['o',843.006417112299],
    ['p',223.810873440285],
    ['q',9.6650623885918],
    ['r',737.602317290553],
    ['s',746.204991087344],
    ['t',981.763279857397],
    ['u',287.579857397504],
    ['v',116.465240641711],
    ['w',181.04385026738],
    ['x',22.0279857397504],
    ['y',189.311942959002],
    ['z',11.8401069518717],
    [' ',2279.00118835413]]
    minChar = typemax(UInt32)
    maxChar = typemin(UInt32)

    for letterScore in letterScores
        letter, score = letterScore
        charVal = UInt32(letter)
        if charVal < minChar
            minChar = charVal
        end
        if charVal > maxChar
            maxChar = charVal
        end
    end

    scores = Float32[0.0 for i in minChar:maxChar]
    for letterScore in letterScores
        letter, score = letterScore
        scores[Int(letter)-minChar+1] = score
    end

    return minChar, maxChar, scores
end

function getKeyOracle():Char
    return '🔥'
end

global const MAX_VALUE_CHECKED = 0x80000000

function getKey(secretText::Vector{UInt32})::Char
    lowestChar, highestChar, frequencies = getFrequencies()
    maxScore = Float32(0)
    maxKey = Char(0)
    for key in typemin(UInt32):MAX_VALUE_CHECKED
        score = 0
        for c in secretText[1:64]
            cGuess = c ⊻ key
            if (lowestChar <= cGuess) && (highestChar >= cGuess)
                score += frequencies[(cGuess - lowestChar) + 1]
            end
        end
        if score > maxScore
            maxScore = score
            maxKey = key
        end
    end
    return Char(maxKey)
end

function getKeyThreaded(secretText::Vector{UInt32})::Char
    lowestChar, highestChar, frequencies = getFrequencies()
    maxScore = Float32(0)
    maxKey = Char(0)
    lk = ReentrantLock()
    threadLoops = Threads.nthreads()*128
    Threads.@threads for thrdIdx in 1:threadLoops
        threadMaxScore = maxScore
        threadMaxKey = maxKey
        for key in typemin(UInt32):threadLoops:MAX_VALUE_CHECKED
            key = key + thrdIdx - 1
            score = UInt32(0)
            for c in secretText
                cGuess = c ⊻ key
                if (lowestChar <= cGuess) && (highestChar >= cGuess)
                    @inbounds score += frequencies[(cGuess - lowestChar) + 1]
                end
            end
            if score > threadMaxScore
                threadMaxScore = score
                threadMaxKey = key
            end
        end
        if threadMaxScore > maxScore
            lock(lk) do
                if threadMaxScore > maxScore
                    maxScore = threadMaxScore
                    maxKey = threadMaxKey
                end
            end
        end
    end
    return Char(maxKey)
end

function getKeyGPUKernel!(secretText::CuDeviceVector{UInt32, 1},
                            keyScores::CuDeviceVector{Float32, 1},
                            frequencies::CuDeviceVector{Float32,1},
                            lowestChar::UInt32,
                            highestChar::UInt32)

    keyIndex = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    score = Float32(0)
    for c in secretText
        key = keyIndex - 1
        cGuess = c ⊻ key
        if (lowestChar <= cGuess) && (highestChar >= cGuess)
            @inbounds score += frequencies[(cGuess - lowestChar) + 1]
        end
    end

    @inbounds keyScores[keyIndex] = score
    return nothing
end

function getKeyGPU(secretText::Vector{UInt32})::Char
    lowestChar, highestChar, frequencies = getFrequencies()
    keyScores = CUDA.fill(0.0f0, MAX_VALUE_CHECKED+1)
    @cuda  threads=1024 blocks=cld(length(keyScores), 1024) getKeyGPUKernel!(
                            cu(secretText),
                            keyScores,
                            cu(frequencies),
                            lowestChar,
                            highestChar)
    maxScore, maxCharIndex = findmax(keyScores)
    return Char(maxCharIndex - 1)
end

function main()
    open("secret.txt") do f
        secretCodedText = ""
        for line in eachline(f)
            secretCodedText *= line
        end
        secretText = base64decode(secretCodedText)
        secretText = Vector(reinterpret(UInt32, secretText))
        secretTextSub = secretText[1:64]
        key = getKeyOracle()
        print("Normal Loop          ")
        @test key == @btime getKey($secretTextSub)
        print("Threaded (", Threads.nthreads(), " threads)")
        @test key == @btime getKeyThreaded($secretTextSub)
        print("GPU                  ")
        @test key == @btime getKeyGPU($secretTextSub)
        text = [Char(UInt32(c) ⊻ UInt32(key)) for c in secretText]
        println("Bounds:")
        println("\tMin: 0x00000000")
        println("\tMax: 0x80000000")
        println("Key: ", key)
        println("Text:")
        println(join(text[1:57]))
    end
end

main()

max_stars_repo_name: billybrickner/juliaDemo

max_issues_repo_name: billybrickner/juliaDemo

max_forks_repo_name: billybrickner/juliaDemo