text:
using Documenter, MIPVerify

makedocs(
    modules = [MIPVerify],
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "MIPVerify",
    authors = "Vincent Tjeng and contributors.",
    pages = [
        "Home" => "index.md",
        "Tutorials" => "tutorials.md",
        "Finding Adversarial Examples" => [
            "finding_adversarial_examples/single_image.md",
            "finding_adversarial_examples/batch_processing.md",
        ],
        "Importing" =>
            ["utils/import_example_nets.md", "utils/import_weights.md", "utils/import_datasets.md"],
        "Working with Neural Net Parameters" => [
            "net_components/overview.md",
            "net_components/layers.md",
            "net_components/nets.md",
            "net_components/core_ops.md",
        ],
    ],
    checkdocs = :exports,
)

deploydocs(repo = "github.com/vtjeng/MIPVerify.jl.git")


text:
# packing and unpacking of Graph related types

using IncrementalInference
using DistributedFactorGraphs
using Manifolds

using Test

##

@testset "Serialization of SamplableBelief types" begin
##

td = Uniform()

ptd = convert(PackedSamplableBelief, td)
utd = convert(SamplableBelief, td)

@test td.a - utd.a |> abs < 1e-10
@test td.b - utd.b |> abs < 1e-10

##
end

##

dfg = initfg()

@testset "hard-coded test of PackedPrior to Prior" begin

pt = PackedPrior("KDE:100:[1.5015]:[-98.8276 -101.803 -98.1296 -96.1897 -99.3076 -99.1881 -101.721 -101.332 -100.431 -97.7293 -96.7652 -99.3806 -95.5593 -104.237 -94.9318 -101.691 -102.255 -98.9559 -99.3386 -99.2361 -102.483 -102.896 -97.0244 -98.9643 -99.4457 -101.298 -103.263 -2.75251 5.14065 0.327863 3.60042 -0.604114 -0.0564047 -0.804898 3.05469 1.4974 1.34657 2.22745 4.78117 1.89485 -0.48091 6.63068 0.63013 -3.11422 1.73705 5.22904 -1.73223 2.47453 1.10584 -0.0179944 3.65585 4.50016 -1.95042 98.2664 98.9983 103.748 100.789 98.4127 101.397 104.364 102.125 96.3685 103.59 99.0581 100.913 101.461 105.211 103.513 99.3325 101.201 98.05 103.508 99.9785 104.624 100.202 100.258 101.579 96.6931 95.4181 299.02 296.804 301.322 298.127 299.578 298.36 296.339 300.156 299.641 297.731 299.822 296.941 295.857 299.482 302.531 301.875 302.192 301.999 300.634 294.084 300.44]")

tt = convert(BallTreeDensity, pt.Z)
@test isa(tt, BallTreeDensity)

F = Prior

upt = convert(F, pt)

@test isa(upt, Prior)

  # TODO add more tests
end

##

fg = initfg()
N=100

doors = zeros(1,4)
doors[1,1:4] = [-100.0;0.0;100.0;300.0]
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
bws2 = Array{Float64,2}(undef, length(bws),1)
bws2[:,1] = bws[:]
doors2 = getPoints(pd);
v1 = addVariable!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[:x1], Prior(kde!(doors2,bws2[:])))

v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
lc = LinearRelative(  Normal(50.0, 2.0) )
f2 = addFactor!(fg, [:x1; :x2], lc)

##

@testset "Testing conversion to packed function node data structure and back" begin

topack = getSolverData(f1)
dd = convert(PackedFunctionNodeData{PackedPrior},topack)
upd = convert(FunctionNodeData{CommonConvWrapper{Prior}}, dd)

@test compare(topack, upd)

topack = getSolverData(f2)
dd =  convert(IncrementalInference.PackedFunctionNodeData{PackedLinearRelative},topack)
upd = convert(IncrementalInference.FunctionNodeData{CommonConvWrapper{LinearRelative}}, dd)

@test compare(topack, upd)

end

@testset "Testing conversion to packed variable node data structure and back" begin

dat = getSolverData(getVariable(fg,:x1))

# dat.BayesNetVertID

pd = packVariableNodeData(dfg, dat)
unpckd = unpackVariableNodeData(dfg, pd)

@test compareFields(dat, unpckd, skip=[:variableType])
@test compareFields(getVariableType(dat), getVariableType(unpckd))
@test isa(getVariableType(dat), ContinuousScalar)
@test isa(getVariableType(unpckd), ContinuousScalar)

end

@testset "test serialization of ManifoldKernelDensity" begin
##

# create a basic manifoldkerneldensity
mkd = manikde!(TranslationGroup(2), [randn(2) for _ in 1:100])

# convert up and down
st = convert(PackedSamplableBelief, mkd)
upk = convert(SamplableBelief, st)

# and check the basics
@test isapprox( getPoints(mkd)[1], getPoints(upk)[1])
@test isapprox( getPoints(mkd)[end], getPoints(upk)[end])

@test mkd.manifold == upk.manifold
@test mkd._partial == upk._partial
@test mkd.infoPerCoord == upk.infoPerCoord

##
end

#

ext: jl


text:
module TokenType
@enum Token begin
    # Reserved word
    INCLUDE = 0                         # include
    DEFINE = 1                          # define
    MAIN = 2                            # main
    CHAR = 3                            # char
    INT = 4                             # int
    FLOAT = 5                           # float
    DOUBLE = 6                          # double
    IF = 7                              # if
    ELSE = 8                            # else
    ELSEIF = 9                          # elseif
    FOR = 10                            # for
    WHILE = 11                          # while
    DO = 12                             # do
    RETURN = 13                         # return
    SWITCH = 14                         # switch
    CASE = 15                           # case
    PRINTF = 16                         # printf
    SCANF = 17                          # scanf
    # Library name
    LIBRARY_NAME = 30                   # <libname.h>, "libname.h"
    # Comment
    SINGLE_COMMENT = 50                 # //
    MULTIPLE_COMMENT = 51               # /**/
    # Identifier
    IDENTIFIER = 60                     # a b c
    # Pointer
    POINTER_IDENTIFIER = 61             # *a *b *c
    # Address
    ADDRESS_IDENTIFIER = 62             # declared &a &b &c
    # Unidentifier
    UNIDENTIFIER = 69                   # undefine
    # Sequence token
    CHAR_SEQ = 70                       # "a b c" -> "a", "b", "c"
    STRING_SEQ = 71                     # "abc efg \n" -> DQUOTE, CHAR_SEQ("abc", "efg "), NEWLINE_ESC, DQUOTE
                                        # 'b' -> QUOTE, CHAR_SEQ("b"), QUOTE
                                        # “123” -> LDQUOTE, CHAR_SEQ("123"), RDQUOTE
    # Constant
    DEC_CON = 80                        # 1 2 3 -1 (-2) -3
    OCT_CON = 81                        # 0101 0234
    HEX_CON = 82                        # 0x12
    FLOAT_CON = 83                      # 1.2f 1.f .5f 1.2
    # Operator
    PLUS = 100                          # +
    MINUS = 101                         # -
    MULTIPLY = 102                      # *
    DIVIDE = 103                        # /
    MODULO = 104                        # %
    ASSIGN = 105                        # =
    NOT = 106                           # !
    BIT_AND = 107                       # &
    BIT_OR = 108                        # |
    BIT_XOR = 109                       # ^
    BIT_NOT = 110                       # ~
    DPLUS = 111                         # ++
    DMINUS = 112                        # --
    ACC_PLUS = 113                      # +=
    ACC_MINUS = 114                     # -=
    ACC_MULTIPLY = 115                  # *=
    ACC_DIVIDE = 116                    # /=
    ACC_MODULO = 117                    # %=
    ACC_BIT_AND = 118                   # &=
    ACC_BIT_OR = 119                    # |=
    ACC_BIT_XOR = 120                   # ^=
    LEFT_SHIFT = 121                    # <<
    RIGHT_SHIFT = 122                   # >>
    ACC_LEFT_SHIFT = 123                # <<=
    ACC_RIGHT_SHIFT = 124               # >>=
    AND = 125                           # &&
    OR = 126                            # ||
    # Comparator
    EQUAL = 150                         # ==
    IDENTIFIER_EQUAL = 151              # ===
    GREATER = 152                       # >
    LESS = 153                          # <
    GREATER_EQUAL = 154                 # >=
    LESS_EQUAL = 155                    # <=
    NOT_EQUAL = 156                     # !=
    # Bracket
    LEFT_PARENTHESE = 170               # (
    RIGHT_PARENTHESE = 171              # )
    LEFT_BRACE = 172                    # {
    RIGHT_BRACE = 173                   # }
    LEFT_BRACKET = 174                  # [
    RIGHT_BRACKET = 175                 # ]
    # Format specifier
    DECIMAL_FORMAT = 180                # %x.xd
    LONG_FORMAT = 181                   # %x.xl
    FLOAT_FORMAT = 182                  # %x.xf
    DOUBLE_FORMAT = 183                 # %x.xlf
    CHAR_FORMAT = 184                   # %c
    STRING_FORMAT = 185                 # %x.xs
    ALERT_ESC = 200                     # \a
    BACKSPACE_ESC = 201                 # \b
    ESCAPE_ESC = 202                    # \e
    FORMFEED_ESC = 203                  # \f
    NEWLINE_ESC = 204                   # \n
    RETURN_ESC = 205                    # \r
    HTAB_ESC = 206                      # \t
    VTAB_ESC = 207                      # \v
    BACKSLASH_ESC = 208                 # \\
    QUOTE_ESC = 209                     # \'
    DQUOTE_ESC = 210                    # \"
    QUESTION_ESC = 211                  # \?
    OCT_ESC = 212                       # \123 \2 \32
    HEX_ESC = 213                       # \u0014
    # Punctuation
    COMMA = 250                         # ,
    SEMICOLON = 251                     # ;
    COLON = 252                         # :
    SHARP = 253                         # #
    QUOTE = 254                         # '
    DQUOTE = 255                        # "
    # New line of stream
    LN_TOKEN = -1
    # End of stream
    EOS_TOKEN = -2
end

function toPunctuationType(symbol::Char)::Union{Token,Nothing}
    if symbol == ','
        return COMMA
    elseif symbol == ';'
        return SEMICOLON
    elseif symbol == ':'
        return COLON
    elseif symbol == '#'
        return SHARP
    elseif symbol == '\''
        return QUOTE
    elseif symbol == '"'
        return DQUOTE
    end
    return nothing
end

function toPreversedType(str::String)::Union{Token,Nothing}
    str = lowercase(str)
    if str == "include"
        return INCLUDE
    elseif str == "define"
        return DEFINE
    elseif str == "main"
        return MAIN
    elseif str == "char"
        return CHAR
    elseif str == "int"
        return INT
    elseif str == "float"
        return FLOAT
    elseif str == "double"
        return DOUBLE
    elseif str == "if"
        return IF
    elseif str == "elseif"
        return ELSEIF
    elseif str == "else"
        return ELSE
    elseif str == "for"
        return FOR
    elseif str == "while"
        return WHILE
    elseif str == "do"
        return DO
    elseif str == "return"
        return RETURN
    elseif str == "switch"
        return SWITCH
    elseif str == "case"
        return CASE
    elseif str == "printf"
        return PRINTF
    elseif str == "scanf"
        return SCANF
    end
    return nothing
end
end


text:
# @generated default_offsets(::Val{N}) where {N} = :(Base.Cartesian.@ntuple $N n -> One())
function default_stridedpointer_quote(::Type{T}, N, Ot) where {T}
    C = 1
    B = 0
    R = Expr(:tuple)
    o = Expr(:tuple)
    for n ∈ 1:N
        push!(R.args, n)
        push!(o.args, Expr(:call, Ot))
    end
    quote
        $(Expr(:meta,:inline))
        StridedPointer{$T,$N,$C,$B,$R}(ptr, x, $o)
    end
end

@generated function default_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :One)
end
@generated function default_zerobased_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :Zero)
end

function dense_quote(N::Int, b::Bool)
    d = Expr(:tuple)
    for n in 1:N
        push!(d.args, b)
    end
    Expr(:call, Expr(:curly, :DenseDims, d))
end
@generated all_dense(::Val{N}) where {N} = dense_quote(N, true)
@generated none_dense(::Val{N}) where {N} = dense_quote(N, false)

@inline function ptrarray0(ptr::Ptr{T}, s::Tuple{Vararg{Integer,N}}, x::Tuple{Vararg{Integer,N}}, ::DenseDims{D}) where {T,N,D}
    PtrArray(default_zerobased_stridedpointer(ptr, x), s, DenseDims{D}())
end
@inline function PtrArray(ptr::Ptr{T}, s::Tuple{Vararg{Integer,N}}, x::Tuple{Vararg{Integer,N}}, ::DenseDims{D}) where {T,N,D}
    PtrArray(default_stridedpointer(ptr, x), s, DenseDims{D}())
end
# @generated function PtrArray(ptr::Ptr{T}, s::Tuple{Vararg{Integer,N}}, x::Tuple{Vararg{Integer,N}}) where {T,N}
#     q = Expr(:block, Expr(:meta,:inline))

# end
@inline zeroindex(A::PtrArray{S,D}) where {S,D} = PtrArray(zstridedpointer(A), size(A), DenseDims{D}())

function ptrarray_densestride_quote(::Type{T}, N, stridedpointer_offsets) where {T}
    last_sx = :s_0
    q = Expr(:block, Expr(:meta,:inline), Expr(:(=), last_sx, static_expr(sizeof(T))))
    t = Expr(:tuple); d = Expr(:tuple);
    n = 0
    while true
        n += 1
        push!(t.args, last_sx)
        push!(d.args, true)
        n == N && break
        new_sx = Symbol(:s_,n)
        push!(q.args, Expr(:(=), new_sx, Expr(:call, :vmul, last_sx, Expr(:ref, :s, n))))
        last_sx = new_sx
    end
    push!(q.args, :(PtrArray($stridedpointer_offsets(ptr, $t), s, DenseDims{$d}())))
    q
end

@generated function PtrArray(ptr::Ptr{T}, s::Tuple{Vararg{Integer,N}}) where {T,N}
    ptrarray_densestride_quote(T, N, :default_stridedpointer)
end
@generated function ptrarray0(ptr::Ptr{T}, s::Tuple{Vararg{Integer,N}}) where {T,N}
    ptrarray_densestride_quote(T, N, :default_zerobased_stridedpointer)
end
# @inline function PtrArray(ptr::StridedPointer

static_expr(N::Int) = Expr(:call, Expr(:curly, :StaticInt, N))
@generated function calc_strides_len(::Type{T}, s::Tuple{Vararg{StaticInt,N}}) where {T, N}
    L = sizeof(T)
    t = Expr(:tuple)
    for n ∈ 1:N
        push!(t.args, static_expr(L))
        L *= s.parameters[n].parameters[1]
    end
    Expr(:tuple, t, static_expr(L))
end
@generated function calc_strides_len(::Type{T}, s::Tuple{Vararg{Any,N}}) where {T, N}
    last_sx = :s_0
    q = Expr(:block, Expr(:meta,:inline), Expr(:(=), last_sx, static_expr(sizeof(T))))
    t = Expr(:tuple)
    for n ∈ 1:N
        push!(t.args, last_sx)
        new_sx = Symbol(:s_,n)
        push!(q.args, Expr(:(=), new_sx, Expr(:call, :vmul, last_sx, Expr(:ref, :s, n))))
        last_sx = new_sx
    end
    push!(q.args, Expr(:tuple, t, last_sx))
    q
end

# @inlines are because we want to make sure the compiler has the chance to avoid the allocation
# @inline function FixedSizeArray{T}(::UndefInitializer, s::Tuple{Vararg{StaticInt,N}}) where {N,T}
#     x, L = calc_strides_len(s)
#     b = MemoryBuffer{T}(undef, L)
#     ptr = VectorizationBase.align(pointer(r))
#     FixedSizeArray(ptr, s, x, b)
# end
# @inline function FixedSizeArray(ptr::Ptr{T}, s::S, x::X, b::NTuple)
#     FixedSizeArray(PtrArray(ptr, s, x), b)
# end

@inline function StrideArray{T}(::UndefInitializer, s::Tuple{Vararg{Integer,N}}) where {N,T}
    x, L = calc_strides_len(T,s)
    b = undef_memory_buffer(T, L ÷ static_sizeof(T))
    # For now, just trust Julia's alignment heuristics are doing the right thing
    # to save us from having to over-allocate
    # ptr = VectorizationBase.align(pointer(b))
    ptr = pointer(b)
    StrideArray(ptr, s, x, b, all_dense(Val{N}()))
end
@inline function StrideArray(ptr::Ptr{T}, s::S, x::X, b, ::DenseDims{D}) where {S,X,T,D}
    StrideArray(PtrArray(ptr, s, x, DenseDims{D}()), b)
end

@inline PtrArray(A::StrideArray) = A.ptr
@inline PtrArray(A::AbstractArray) = PtrArray(stridedpointer(A), size(A), dense_dims(A))

@generated rank_to_sortperm_val(::Val{R}) where {R} = :(Val{$(ArrayInterface.rank_to_sortperm(R))}())
@inline function similar_layout(A::AbstractStrideArray{S,D,T,N,C,B,R}) where {S,D,T,N,C,B,R}
    permutedims(similar(permutedims(A, rank_to_sortperm_val(Val{R}()))), Val{R}())
end
@inline function similar_layout(A::AbstractArray)
    b = preserve_buffer(A)
    GC.@preserve b begin
        similar_layout(PtrArray(A))
    end
end
@inline function Base.similar(A::AbstractStrideArray{S,D,T}) where {S,D,T}
    StrideArray{T}(undef, size(A))
end

# function calc_padding(nrow::Int, T)
#     W = VectorizationBase.pick_vector_width(T)
#     W > nrow ? VectorizationBase.nextpow2(nrow) : VectorizationBase.align(nrow, T)
# end


text:
using Query
using DataFrames
using NamedTuples

# We want a DataFrame without Nullable's, so
# we need to construct a DataFrame with columns of Array
columns=[]
push!(columns, ["John", "Sally", "Kirk"])
push!(columns, [23., 42., 59.])
push!(columns, [3,5,2])
df = DataFrame(columns, [:name, :age, :children])

# It happens that the whole thing also works fairly well
# with standard julia generator syntax. Probably would make
# sense to think hard how that could be integrated better...
# In some way the generator syntax right now is a little bit
# like a "query syntax light".

x = collect(i.name for i in Query.query(df) if i.age>30.)
println(x)

x = collect(@NT(Name=>i.name) for i in Query.query(df) if i.age>30.)
println(x)

x = collect(Query.query(x), DataFrame)
println(x)


text:
"""
    KMarkovUpdater

Updater that stores the k most recent observations as the belief.

Example:

```julia
up = KMarkovUpdater(5)
s0 = initialstate(pomdp, rng)
initial_observation = generate_o(pomdp, s0, rng)
initial_obs_vec = fill(initial_observation, 5)
hr = HistoryRecorder(rng=rng, max_steps=100)
hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)
```
"""
struct KMarkovUpdater <: Updater
    k::Int
end

function initialize_belief(bu::KMarkovUpdater, obs_vec::AbstractVector)
    if length(obs_vec) != bu.k
        error("KMarkovUpdater: The length of the initial observation vector
               does not match the number of observation to stack\n"*throw_example(bu))
    end
    return obs_vec
end

function update(bu::KMarkovUpdater, old_b::AbstractVector{O}, action, obs) where {O}
    obs_stacked = Vector{O}(undef, bu.k)
    if !isa(obs, O)
        error("KMarkovUpdater: Observation did not match previous observation type.\n"*throw_example(bu))
    end
    for i=1:bu.k-1
        obs_stacked[i] = old_b[i+1]
    end
    obs_stacked[bu.k] = obs
    return obs_stacked
end

function initialize_belief(bu::KMarkovUpdater, obs_vec)
    error("KMarkovUpdater: To initialize the belief, pass in a vector of observation.\n"*throw_example(bu))
end

function throw_example(bu::KMarkovUpdater)
    example = """
    Did you forget to pass the initial observation to the simulator?
    Example:
    ```julia
    up = KMarkovUpdater(5)
    s0 = initialstate(pomdp, rng)
    initial_observation = generate_o(pomdp, s0, rng)
    initial_obs_vec = fill(initial_observation, 5)
    hr = HistoryRecorder(rng=rng, max_steps=100)
    hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)
    ```
    """
    return example
end


text:
using LockPatternComplexity
using Documenter

DocMeta.setdocmeta!(LockPatternComplexity, :DocTestSetup, :(using LockPatternComplexity); recursive=true)

makedocs(;
    modules=[LockPatternComplexity],
    authors="Jaan Tollander de Balsch",
    repo="https://github.com/jaantollander/LockPatternComplexity.jl/blob/{commit}{path}#{line}",
    sitename="LockPatternComplexity.jl",
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://jaantollander.github.io/LockPatternComplexity.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/jaantollander/LockPatternComplexity.jl",
)


text:
__precompile__()
module RosDataProcess

using LinearAlgebra, PyCall, PyPlot, Colors, Statistics, JLD2

const matplotlib2tikz = PyNULL()
const yaml = PyNULL()
const rosbag = PyNULL()

function __init__()
  try
    copy!(matplotlib2tikz, pyimport("tikzplotlib"))
  catch e
    println("Could not import tikzplotlib, trying matploblib2tikz instead")
    copy!(matplotlib2tikz, pyimport("matplotlib2tikz"))
  end

  copy!(yaml, pyimport_conda("yaml", "pyyaml"))

  try
    copy!(rosbag, pyimport("rosbag"))
  catch e
    println("Could not import rosbag.")
    println("Bags will not load, but other functions will work.")
  end

end

export TimeSeries, time_type, data_type, time_eltype, get_time, get_data,
mutate_time, mutate_data, map_time, map_data, indices_match

include("time_series.jl")

export normalize_start, intersect_intervals, intersect_regular, interpolate,
get_at_time, get_at_nearest_time, intersect_interpolate, intersect_series,
select_times, select_indices, standard_error, differentiate

include("data_process.jl")

export plot, plot_standard_error, plot_mean, plot_mean, plot_trials,
generate_colors, first_horizontal_intersection

include("plotting.jl")

export AnnotatedBag, load_yaml, load_yaml_data, load_bag, get_name,
get_bag_names, filter_extension, load_directory, get_param, filter_bags,
get_range, sort_bags, to_sec, get_topic_names, first_bag, read_topic,
read_series, to_file_name, save_latex

include("bag_process.jl")

end # module

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
"""
Submodule providing helper functions for generating plots.
"""
module Plotter
using ..CarrierCapture: Potential, conf_coord

export plot_pots, plot_pot!, plot_ccs, plot_cc!

using Plots, LaTeXStrings

"""
    plot_pot!(pot; lplt_wf = false, plt = nothing, color = Nothing, label = "", scale_factor = 2e-2)

Plots `Potential` with its' wave functions if `lplt_wf` is `true`.
"""
function plot_pot!(pot::Potential; lplt_wf = false, plt = nothing, color = Nothing, label = "", scale_factor = 2e-2)
    if plt == nothing
        plt = plot()
    end
    label = if label == "" pot.name else label end
    color = if color == nothing "black" else color end

    # plot wave functions
    if lplt_wf
        ϵ = pot.ϵ; χ = pot.χ
        for i = 1:length(ϵ)
            plot!(plt, pot.Q, χ[i, :] * scale_factor .+ ϵ[i], lw = 0,
                  fillrange = [χ[i, :] * 0 .+ ϵ[i], χ[i, :] * scale_factor .+ ϵ[i]],
                  color = color, alpha = 0.5, label = "")
        end
    end

    # plot function
    plot!(plt, pot.Q, pot.E, lw = 4, color = color, label = label)

    # plot data
    if size(pot.QE_data)[1] > 1
        scatter!(plt, pot.QE_data.Q, pot.QE_data.E, color = color, label = "")
    end
    xaxis!(L"\ Q (amu^{1\/2} Å) \ (^{}$$"); yaxis!(L"\ Energy  (eV) \ (^{}$$")

    return plt
end

"""
    plot_cc!(cc; plt = nothing, color = nothing, label = "")

Plots capture coefficient `cc.capt_coeff` as a function of `1000/T`.
"""
function plot_cc!(cc::conf_coord; plt = nothing, color = nothing, label = "")
    if plt == nothing
        plt = plot()
    end
    label = if label == "" cc.name else label end
    # color = if color==Nothing "black" else color end
    plot!(plt, 1000 ./ cc.temperature, cc.capt_coeff, lw = 4, label = label)
    xaxis!(L"\ 1000\/T (K^{-1}) \ (^{}$$"); yaxis!(L"C (cm^{3}\/s) \ (^{}$$", :log10)
    return plt
end

end


text:
"""
    loadlocalizations(source, inputformat)

Load localizations with the given format from the path or stream provided.
"""
function loadlocalizations(source::Union{String, IO}, inputformat::InputFormat)
    if inputformat == nikonelementstext
        loadfrom_nikonelementstext(source)
    elseif inputformat == nikonelementsbinary
        loadfrom_nikonelementsbinary(source)
    end
end

function loadfrom_nikonelementstext(source::Union{String, IO})
    file = CSV.File(source, delim = '\t', normalizenames = true, types = Dict(1 => String))

    localizations = Vector{Localization}()
    i = 1
    for row ∈ file
        loc =  Localization(i, row.Channel_Name, row.Xwc, row.Ywc, row.Zwc,
            row.Lateral_Localization_Accuracy, row.Frame, row.Length)
        push!(localizations, loc)
        i += 1
    end
    localizations
end

#=
function loadfrom_nikonelementsbinary(source::Union{String, IO})
    file = CSV.File(source, delim = '\t', normalizenames = true, types = Dict(1 => String))

    open(source) do f
        skip(f, 69) # 69-byte header
        nlocalizations = read(f, Int32)
        A = Matrix{Float32}(undef, 18, nlocalizations)
        A = read!(f, A) # table of 18 columns containing floating point values

        [Localization(i, "1", col[3], col[4], col[18],) (i, col) ∈ enumerate(eachrow(A))]
    end
    localizations
end

From MATLAB
mol.x = A(1:18:end);
mol.y = A(2:18:end);
mol.xc = A(3:18:end);
mol.yc = A(4:18:end);
mol.h = A(5:18:end);
mol.area = A(6:18:end);
mol.width = A(7:18:end);
mol.phi = A(8:18:end);
mol.Ax = A(9:18:end);
mol.bg = A(10:18:end);
mol.I = A(11:18:end);
mol.cat = floor(A(12:18:end)/1.401e-45); ==> reinterpret(Int32, A[12,:])
mol.valid = floor(A(13:18:end)/1.401e-45);
mol.frame = floor(A(14:18:end)/1.401e-45);
mol.length = floor(A(15:18:end)/1.401e-45);
mol.link = floor(A(16:18:end)/1.401e-45);
mol.z = A(17:18:num*18);
mol.zc = A(18:18:num*18);

# where is xwc, ywc, zwc? Channel name? Lateral Localization Accuracy?
# I believe these may all be calculated values that Elements populates when saving a text file
# LLA is calculated by a formula which I could reproduce here, but the warp corrections
# would be stored in the microscope. So... maybe I don't want to read Nikon binary files?

=#


text:
abstract type AbstractGrammar{R} end

chart_key(g::AbstractGrammar{R}) where {R} = chart_key(R)
rule_type(g::AbstractGrammar{R}) where {R} = R

function productions end
function terminal_productions end
function start_symbol end

abstract type AbstractRule{R} end
chart_key(::Type{<:AbstractRule{R}}) where {R} = R

lhs(p::Pair) = first(p)
rhs(p::Pair) = last(p)

struct SimpleRule <: AbstractRule{Symbol}
    rule::Pair{Symbol, Vector{Symbol}}
end

rhs(r::SimpleRule) = rhs(r.rule)
lhs(r::SimpleRule) = lhs(r.rule)
score(r::SimpleRule) = 1.0

Base.convert(::Type{SimpleRule}, p::Pair{Symbol, Vector{Symbol}}) = SimpleRule(p)

struct SimpleGrammar <: AbstractGrammar{SimpleRule}
    productions::Vector{SimpleRule}
    categories::Dict{String, Vector{Symbol}}
    start::Symbol
end

productions(g::SimpleGrammar) = g.productions
start_symbol(g::SimpleGrammar) = g.start

function terminal_productions(g::SimpleGrammar, tokens::AbstractVector{<:AbstractString})
    R = rule_type(g)
    result = Arc{R}[]
    for (i, token) in enumerate(tokens)
        for category in get(g.categories, token, Symbol[])
            push!(result, Arc{R}(i - 1, i, category => Symbol[], Arc{R}[], 1.0))
        end
    end
    result
end

struct SimpleWeightedRule <: AbstractRule{Symbol}
    rule::Pair{Symbol, Vector{Symbol}}
    score::Float64
end

rhs(r::SimpleWeightedRule) = rhs(r.rule)
lhs(r::SimpleWeightedRule) = lhs(r.rule)
score(r::SimpleWeightedRule) = r.score

Base.convert(::Type{SimpleWeightedRule}, t::Tuple{Pair{Symbol, Vector{Symbol}}, Real}) = SimpleWeightedRule(t...)

struct SimpleWeightedGrammar <: AbstractGrammar{SimpleWeightedRule}
    productions::Vector{SimpleWeightedRule}
    categories::Dict{String, Vector{Pair{Symbol, Float64}}}
    start::Symbol
end

productions(g::SimpleWeightedGrammar) = g.productions
start_symbol(g::SimpleWeightedGrammar) = g.start

function terminal_productions(g::SimpleWeightedGrammar, tokens::AbstractVector{<:AbstractString})
    R = rule_type(g)
    result = Arc{R}[]
    for (i, token) in enumerate(tokens)
        for (category, weight) in get(g.categories, token, Symbol[])
            push!(result, Arc{R}(i - 1, i, SimpleWeightedRule(category => Symbol[], weight), Arc{R}[], weight))
        end
    end
    result
end


text:
#!/bin/bash
#=
project=$(dirname $0)/..
sysimage=$project/sysimage.so

if [ ! -e "$sysimage" ]; then
	# echo not found
	julia --project="$project" "$project/src/precompile.jl"
fi

exec julia --project=$project --sysimage="$sysimage" --startup-file=no \
    --color=yes -e 'include(popfirst!(ARGS))' "${BASH_SOURCE[0]}" "$@"
=#

using Logging, Dates
include("./preprocess/preprocess.jl")
include("./training/training.jl")
include("./LazyDataLoader.jl")
include("./cliparser.jl")

function main(options)

	log_io = open("train-$(Dates.format(now(), "yyyy-mm-ddTHH-MM-SS")).log", "w+")
	global_logger(SimpleLogger(log_io))

	feature_shape = Dict(
		"chroma" => (12, 582),
		"melspectrogram" => (128, 582),
		"mfcc" => (26, 582),
	)[options["selectedfeature"]]

	if options["shufflesamples"]
		shuffle_rows(options["tagsfile"])
	end

	if options["preprocess"]
		@info "Preprocessing audio..."
		preprocess_audio(options["selectedfeature"], options["tagsfile"],
			options["pathcolumn"], options["samplespath"],
			options["audiofeaturesdir"])
	end

	if options["trainnetwork"]
		(train_set, validation_set, test_set), tagnames =
			preprocess_tags(options["tagsfile"], options["pathcolumn"],
				options["excludedcolumns"], partitions=(13, 1, 2))

		for set in (train_set, validation_set, test_set)
			# Append features path and remove file extension
			for (i,) in enumerate(set[1])
				set[1][i] = joinpath(options["audiofeaturesdir"], set[1][i])[1:end-4]
			end
		end

		# Create data loader to load data from disk to gpu on demand
		train_set = LazyDataLoader(create_feature_loader(feature_shape), train_set, batchsize=options["batchsize"])

		@info "Loading validation set..."
		validation_set = (load_features(validation_set[1], feature_shape), copy(validation_set[2]))

		train_cnn(feature_shape, train_set, validation_set,
			options["epochs"], options["modeloutputdir"], options["dumpoutputdir"],
			log_io, options["eta"], options["model"])
	end

	close(log_io)

end

function load_features(filepaths, feature_shape)
	features = Array{Float32}(undef, feature_shape..., length(filepaths))
	for (i, filepath) in enumerate(filepaths)
		features[:, :, i] = reshape(reinterpret(Float32, read(filepath)), feature_shape...)
	end

	return features
end

function create_feature_loader(feature_shape)
	return ((filepaths, tags),) -> gpu((load_features(filepaths, feature_shape), copy(tags)))
end

main(options)

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
col1 = open("col1.txt")
col2 = open("col2.txt")
merge = open("merge.txt", "w")
for line in readlines(col1)
  c1 = chomp(line)
  c2 = chomp(readline(col2))
  println(merge, join([c1,c2], "\t"))
end


text:
"""
    autocor(x,lag)

Returns autocorrelation function for a distance of lag, i.e. between values of
```x[t]``` and ```x[t+lag]```. The autocorrelation function is normalized such
that ```autocor(x,0) = 1```. This is equvivalent to the quantity
``\\Gamma_X(t)`` of equation (4.61) in [Gattringer/Lang](@cite).
"""
function autocor(x, lag)
    # (wasteful in terms of allocations but clear)
    z = x .- mean(x)
    a = sum(z[1+lag:end].*z[1:end-lag])/sum(z.*z)
    return a
end
"""
    autocor(x)

Returns autocorrelation function for a set of lags from ``0`` up to the closest
integer to ``10 \\text{log}_{10}(l)``. Here ``l`` is the number of elements in
```x```.
"""
function autocor(x)
    lx   = length(x)
    lags = collect(0:min(lx-1, round(Int,10*log10(lx))))
    a = zeros(eltype(x),length(lags))
    for i in 1:length(lags)
        a[i] = autocor(x, lags[i])
    end
    return a
end
"""
    autotimeexp(x)

Returns ``\\max(1,\\tau)`` where ``\\tau`` is the exponential autocorrelation
time of a series of measurements ``x``. This is obtained by fitting the
autocorrelation function to an exponential function of the form
``A \\exp(\\frac{t}{\\tau})``.
"""
function autotimeexp(O)
    a = autocor(O)
    @. modelτ(x,p) = abs.(p[2])*exp(-x/p[1])
    x = collect(1:length(a))
    c = curve_fit(modelτ, x, a, ones(2))
    τ = c.param[1]
    return max(one(τ),τ)
end


text:
export lbi!, set_live_nodes!

let state::Float64 = 0.
	global lbi_newstate(x) = (state=x)
	global lbi_getstate() = (state)
	global lbi_resetstate() = lbi_newstate(0)
end

"""
	lbi!(r::TreeNode{LBIData}, τ)
	lbi!(t::Tree{LBIData}, τ)
"""
function lbi!(r::TreeNode{LBIData}, τ::Float64 ; normalize=true)
	# set_to_zero!(r::TreeNode{LBIData})
	get_message_up!(r::TreeNode{LBIData}, τ)
	send_message_down!(r::TreeNode{LBIData}, τ)
	compute_lbi!(r::TreeNode{LBIData})
	if normalize
		normalize_lbi!(r)
	end
	return nothing
end
lbi!(r::TreeNode{LBIData}, τ; normalize=true) = lbi!(r, convert(Float64, τ), normalize=true)
function lbi!(t::Tree{LBIData}, τ; normalize=true)
	if length(findall(x->x.data.alive, t.lleaves)) < 1
		# @warn "Cannot compute LBI for tree with dead leaves."
		return true
	end
	lbi!(t.root, τ, normalize=normalize)
	return false
end
lbi!(t::Tree{EvoData}, τ; normalize=true) = error("Cannot compute LBI for type Tree{EvoData}. Read the tree with keyword `LBIData`.")

"""
	set_to_zero!(n::TreeNode{LBIData})
"""
function set_to_zero!(n::TreeNode{LBIData})
	n.data.lbi = 0.
	for c in n.child
		set_to_zero!(c)
	end
end

"""
	compute_lbi!(n::TreeNode{LBIData})
"""
function compute_lbi!(n::TreeNode{LBIData})
	if n.data.alive
		n.data.lbi = n.data.message_down
		for c in n.child
			compute_lbi!(c)
			n.data.lbi += c.data.message_up
		end
	else
		n.data.lbi = 0.
		for c in n.child
			set_to_zero!(c)
		end
	end
end

"""
	get_message_up!(n::TreeNode{LBIData}, τ)

Get message going up from node `n`. Ask for up messages `m_c` for all children `c` of `n`. Return `exp(-t/τ) * sum(m_c) + τ * (1 - exp(-t/τ)` where `t=n.data.tau`. Field `message_up` in `n`'s data is modified.
"""
function get_message_up!(n::TreeNode{LBIData}, τ::Float64)
	n.data.message_up = 0.
	n.data.message_down = 0.
	for c in n.child
		n.data.message_up += get_message_up!(c, τ)
	end
	n.data.message_up *= exp(-n.data.tau/τ)
	if n.data.alive
		(n.data.message_up += τ*(1 - exp(-n.data.tau/τ)))
	end
	return n.data.message_up
end

"""
	send_message_down!(n::TreeNode{LBIData}, τ)

Send message going down from node `n`. Field `c.message_down` is modified for all children `c` of `n`. No return value.

## Note
It's a bit easier to think of this function as operating on `c.anc` where `c` is a child of `n`. Maybe I should code it like this.
"""
function send_message_down!(n::TreeNode{LBIData}, τ::Float64)
	for c1 in n.child
		c1.data.message_down = n.data.message_down
		for c2 in n.child
			if c1 != c2
				c1.data.message_down += c2.data.message_up
			end
		end
		c1.data.message_down *= exp(-c1.data.tau/τ)
		if c1.data.alive
			c1.data.message_down += τ*(1 - exp(-c1.data.tau/τ))
		end
		send_message_down!(c1,τ)
	end
	return nothing
end

"""
	normalize_lbi!(r::TreeNode{LBIData})
	normalize_lbi!(t::Tree)
"""
function normalize_lbi!(r::TreeNode{LBIData})
	max_lbi = get_max_lbi(r::TreeNode{LBIData})
	normalize_lbi!(r::TreeNode{LBIData}, max_lbi)
end
normalize_lbi!(t::Tree) = normalize_lbi!(t.root)
function normalize_lbi!(n::TreeNode{LBIData}, max_lbi::Float64)
	if max_lbi == 0.
		@error "Maximum LBI is 0, cannot normalize."
	end
	for c in n.child
		normalize_lbi!(c, max_lbi)
	end
	n.data.lbi /= max_lbi
end

"""
	get_max_lbi(r::TreeNode{LBIData})

Return maximal LBI value of the subtree below `r`.
"""
function get_max_lbi(r::TreeNode{LBIData})
	lbi_resetstate()
	_get_max_lbi(r)
	return lbi_getstate()
end
function _get_max_lbi(n::TreeNode{LBIData})
	for c in n.child
		_get_max_lbi(c)
	end
	if n.data.lbi > lbi_getstate()
		lbi_newstate(n.data.lbi)
	end
end

"""
	set_live_nodes!(n::TreeNode{LBIData}; datemin=missing, datemax=missing)
	set_live_nodes!(t::Tree; datemin=missing, datemax=missing)

Recursively sets the `.data.alive` field of tree nodes.
- If `n` is a leaf node, `n.data.date` must strictly be between `datemin` and `datemax` (if not missing).
- If `n` is a non-terminal node, at least one of its children must be alive.
"""
function set_live_nodes!(n::TreeNode{LBIData}; datemin=missing, datemax=missing, set_leaves=true)
	for c in n.child
		set_live_nodes!(c, datemin=datemin, datemax=datemax, set_leaves=set_leaves)
	end
	if set_leaves && n.isleaf
		n.data.alive = (ismissing(datemin) || n.data.date > datemin) && (ismissing(datemax) || n.data.date < datemax)
	elseif !n.isleaf
		n.data.alive = mapreduce(x->x.data.alive, |, n.child, init=false)
	end
	return nothing
end
set_live_nodes!(t::Tree; datemin=missing, datemax=missing, set_leaves=true) = set_live_nodes!(t.root, datemin=datemin, datemax=datemax, set_leaves=set_leaves)

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
module FoldsCUDATests

import CUDA
using GPUArrays
using Test
using TestFunctionRunner

include("utils.jl")

module TestGeneric
using ..Utils: include_tests
include_tests(@__MODULE__, joinpath(@__DIR__, "generic"))
end

module TestGPU
using ..Utils: include_tests, should_test_gpu
include_tests(@__MODULE__)
should_test_module() = should_test_gpu()
end

""" A list of tests to be run from UnionArrays.jl """
const UNIONARRAYS_TESTS = [TestGPU.TestTypeChangingAccumulators]

function runtests_unionarrays(; kwargs...)
    TestFunctionRunner.run(UNIONARRAYS_TESTS; kwargs...)
end

function before_test_module()
    if !CUDA.functional()
        @info "CUDA.jl not functional"
        return
    end

    GPUArrays.allowscalar(false)

    if lowercase(get(ENV, "CI", "false")) == "true"
        CUDA.versioninfo()
        println()
    end
end

end # module


text:
function make_constraints(::Type{Val{nothing}}, args, kwargs, verbose)
    verbose && println("No constraints given.")
    return x -> [0.0]
end

function make_constraints(eqs::Vector, args, kwargs, verbose)
    verbose && println("Converting ieqcons to a single constraint function.")
    return x -> [f(x, args...; kwargs...) for f in eqs]
end

function make_constraints(eqs::Function, args, kwargs, verbose)
    verbose && println("Single constraint function given in f_ieqcons.")
    return x -> eqs(x, args...; kwargs...)
end

make_constraints(eqs, args, kwargs, verbose) = make_constraints(Val{eqs}, args, kwargs, verbose)

# atualiza melhores informações das partículas
function update_position!(particle::Particles)
    i_update = (particle.fx .< particle.fp) .& particle.fs
    particle.p[i_update] = copy(particle.x[i_update])
    particle.fp[i_update] = particle.fx[i_update]
end

# funçõa auxiliar para ajuda na busca local
verify(x, size) = (abs(x) % size == 0) ? size : mod(x,size)

# função que atualiza o valor dos melhores vizinhos com base nos dados
# dos n vizinhos mais próximos de cada partícula
function findBestLocal!(particle::Particles, n::Integer)
    sortperm!(particle.order, particle.fx, rev=true)
    for i in particle.order
        temp = verify.((i-n):(i+n), particle.swarmsize)
        val, ind = findmin(particle.fp[particle.order[temp]])
        particle.fg[i] = val
        particle.g[i] = particle.p[particle.order[temp]][ind]
    end
end

# função que atualiza os valores dos melhores vizinhos com base na melhor partícula
function findBestGlobal!(particle::Particles)
    val, ind = findmin(particle.fp)
    particle.fg = fill(val, particle.swarmsize)
    particle.g = fill(particle.p[ind], particle.swarmsize)
end


text:
#!/usr/bin/env julia

# Calculate Mandelbrot set using OpenCL

if VERSION < v"0.4-"
  typealias UInt16 Uint16
end

import OpenCL
const cl = OpenCL

import GR

const mandel_kernel = "
#pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable
__kernel void mandelbrot(__global double2 *q, __global ushort *output,
                         double const min_x, double const max_x,
                         double const min_y, double const max_y,
                         ushort const width, ushort const height,
                         ushort const iters)
{
    int ci = 0, inc = 1;
    int gid = get_global_id(0);
    double nreal, real = 0;
    double imag = 0;

    q[gid].x = min_x + (gid % width) * (max_x - min_x) / width;
    q[gid].y = min_y + (gid / width) * (max_y - min_y) / height;

    output[gid] = iters;

    for (int curiter = 0; curiter < iters; curiter++) {
        nreal = real * real - imag * imag + q[gid].x;
        imag = 2 * real * imag + q[gid].y;
        real = nreal;

        if (real * real + imag * imag >= 4) {
            output[gid] = ci;
            return;
        }
        ci += inc;
        if (ci == 0 || ci == 255)
            inc = -inc;
    }
}"

for device in reverse(cl.available_devices(cl.platforms()[1]))
    global ctx, queue, prg

    ctx = cl.Context(device)
    queue = cl.CmdQueue(ctx)
    try
        prg = cl.Program(ctx, source = mandel_kernel) |> cl.build!
        println(device[:name])
        break
    end
end

function calc_fractal(q, min_x, max_x, min_y, max_y, width, height, iters)
    global ctx, queue, prg

    output = Array(UInt16, size(q))

    q_opencl = cl.Buffer(Complex128, ctx, (:r, :copy), hostbuf=q)
    output_opencl = cl.Buffer(UInt16, ctx, :w, length(output))

    k = cl.Kernel(prg, "mandelbrot")
    cl.call(queue, k, length(q), nothing, q_opencl, output_opencl,
            min_x, max_x, min_y, max_y, UInt16(width), UInt16(height), UInt16(iters))
    cl.copy!(queue, output, output_opencl)

    return output
end

function create_fractal(min_x, max_x, min_y, max_y, width, height, iters)
    q = zeros(Complex128, (width, height))

    output = calc_fractal(q, min_x, max_x, min_y, max_y, width, height, iters)

    return output
end

x = -0.9223327810370947027656057193752719757635
y = 0.3102598350874576432708737495917724836010

f = 0.5
for i in 1:200
    tic()
    image = create_fractal(x-f, x+f, y-f, y+f, 500, 500, 400)
    dt = toq()

    @printf("Mandelbrot created in %f s\n", dt)

    GR.clearws()
    GR.setviewport(0, 1, 0, 1)
    GR.setcolormap(13)
    GR.cellarray(0, 1, 0, 1, 500, 500, image + 1000)
    GR.updatews()

    f *= 0.9
end


text:
add_verbose_scope(:StabSub)

###############################################################################
#
#  Tools for ZpnGModules
#
###############################################################################

#
#  Lifts a matrix from F_p to Z/p^nZ
#

function lift(M::fq_nmod_mat, R::Nemo.NmodRing)

  x=factor(fmpz(R.n))
  @assert length(x.fac)==1
  @assert order(parent(M[1,1]))==first(keys(x.fac))
  N=zero_matrix(R,rows(M),cols(M))
  for i=1:rows(M)
    for j=1:cols(M)
      N[i,j]=FlintZZ(coeff(M[i,j],0))
    end
  end
  return N

end

#
#  Action of a matrix on an element of the group
#

function *(x::GrpAbFinGenElem, M::nmod_mat)

  G=parent(x)
  @assert ngens(G)==rows(M)
  R=parent(M[1,1])
  coeff=zero_matrix(R,1,cols(M))
  for i=1:cols(M)
    coeff[1,i]=x.coeff[1,i]
  end
  y=coeff*M
  return G(lift(y))

end

#
#  Smith Normal Form for a ZpnGModule
#

function Nemo.snf(M::ZpnGModule)

  A=M.V
  G=M.G
  if issnf(A)
    return M, GrpAbFinGenMap(A,A, identity_matrix(FlintZZ, ngens(A)))
  end
  S,mS=snf(A)
  W=[mS(s) for s in gens(S)]
  R=M.R
  H=Array{nmod_mat,1}(length(G))
  for i=1:length(G)
    N=zero_matrix(R, ngens(S),ngens(S))
    for j=1:length(W)
      y=mS\(W[j]*G[i])
      for k=1:ngens(S)
        N[j,k]=y.coeff[1,k]
      end
    end
    H[i]=N
  end
  return ZpnGModule(S,H), mS

end

function is_stable(act::Array{T, 1}, mS::GrpAbFinGenMap) where T <: Map{GrpAbFinGen, GrpAbFinGen}

  S=mS.header.domain
  for s in gens(S)
    x=mS(s)
    for g in act
      if !haspreimage(mS,g(x))[1]
        return false
      end
    end
  end
  return true

end

function issubmodule(M::ZpnGModule, S::nmod_mat)

  @assert issnf(M.V)
  s, ms=subm_to_subg(M,S)
  for x in gens(s)
    el=ms(x)
    for g in M.G
      if !haspreimage(ms,el*g)[1]
        return false
      end
    end
  end
  return true

end

#
#  Given a group $G$ and a group of automorphisms of G, this function returns the corresponding ZpnGModule
#

function action(V::GrpAbFinGen, act::Array{T,1}) where T<: Map{GrpAbFinGen, GrpAbFinGen}

  expon=Int(exponent(V))
  @assert length(factor(order(V)).fac)==1
  RR=ResidueRing(FlintZZ, expon, cached=false)
  act_mat=Array{nmod_mat,1}(length(act))
  for z=1:length(act)
    A=zero_matrix(RR,ngens(V), ngens(V))
    for i=1:ngens(V)
      y=act[z](V[i])
      for j=1:ngens(V)
        A[i,j]=y.coeff[1,j]
      end
    end
    act_mat[z]=A
  end
  return ZpnGModule(V,act_mat)

end

#################################################################################
#
#  Duality
#
#################################################################################

#=
  To get the transpose map of a homomorphism, and so the action on the dual module,
  you need to transpose the matrix and multiply or the elements for a power of p, assuming that the group is in snf.
  Precisely, let p^e be the exponent of the group and let p^t_i be the powers of p appearing in the snf. Then
  ( p^(e-t_1)     )                  ( p^(e-t_1)     )
  (          ...  )  x ^t A  = A  x  (          ...  )
  (              1)                  (              1)
=#

function dual_module(M::ZpnGModule)

  @assert issnf(M.V)
  G1=deepcopy(M.G)
  for i=1:length(G1)
    for j=1:ngens(M.V)-1
      for k=j+1:ngens(M.V)
        x=divexact(M.V.snf[k], M.V.snf[j])
        G1[i][j,k]=divexact(G1[i][j,k].data, x)
        G1[i][k,j]*=x
      end
    end
    transpose!(G1[i])
  end
  return ZpnGModule(M.V, G1)

end

function _dualize(M::nmod_mat, V::GrpAbFinGen, v::Array{fmpz,1})
  #
  #  First, compute the kernel of the corresponding homomorphisms
  #
  K=DiagonalGroup([V.snf[end] for j=1:rows(M)])
  A=lift(transpose(M))
  for j=1:rows(A)
    for k=1:cols(A)
      A[j,k]*=v[j]
    end
  end
  mH=Hecke.GrpAbFinGenMap(V,K,A)
  newel=kernel_as_submodule(mH)
  return MatrixSpace(parent(M[1,1]),rows(newel), cols(newel), false)(newel)

end

function _dualize_1(M::nmod_mat, snf_struct::Array{fmpz,1})

  A=nullspace(transpose(M))
  B=vcat(transpose(A),zero_matrix(M[1,1].parent, cols(A),cols(A)))
  for j=1:cols(A)
    B[rows(A)+j,j]=snf_struct[j]
  end
  S=nullspace(B)
  C=vcat(transpose(A),zero_matrix(M[1,1].parent, cols(A),cols(A)))
  return S*C

end

function kernel_as_submodule(h::GrpAbFinGenMap)
  G = domain(h)
  H = codomain(h)
  hn, t = hnf_with_transform(vcat(h.map, rels(H)))
  for i=1:rows(hn)
    if iszero_row(hn, i)
      return sub(t, i:rows(t), 1:ngens(G))
    end
  end
  error("JH")
end

########################################################################################################
#
#  Quotients and subgroups of ZpnGModules
#
########################################################################################################

function quo(M::ZpnGModule, S:: nmod_mat)

  Q,mQ=quo(M.V,lift(S))
  return ZpnGModule(Q,M.G), mQ

end

function sub(M::ZpnGModule, S::nmod_mat)

  sg,msg=sub(M.V,lift(S))
  G=Array{nmod_mat,1}(length(M.G))
  for k=1:length(M.G)
    A=zero_matrix(M.R, ngens(sg), ngens(sg))
    for i=1:ngens(sg)
      x=msg(sg[i])*M.G[k]
      x=haspreimage(msg, x)[2].coeff
      for j=1:ngens(sg)
        A[i,j]=x[1,j]
      end
    end
    G[k]=A
  end
  return ZpnGModule(sg,G), msg

end

function sub(M::ZpnGModule, n::Int)

  sg,msg=sub(M.V,n)
  G=Array{nmod_mat,1}(length(M.G))
  for k=1:length(M.G)
    A=zero_matrix(M.R, ngens(sg), ngens(sg))
    for i=1:ngens(sg)
      x=msg(sg[i])*M.G[k]
      x=haspreimage(msg, x)[2].coeff
      for j=1:ngens(sg)
        A[i,j]=x[1,j]
      end
    end
    G[k]=A
  end
  return ZpnGModule(sg,G), msg

end

function _exponent_p_sub(M::ZpnGModule)

  @assert issnf(M.V)
  G=M.G
  V=M.V
  p=M.p
  F, a = Nemo.FiniteField(p, 1, "a", cached=false)
  v=fmpz[divexact(V.snf[i],p) for i=1:ngens(V)]
  G1=Array{fq_nmod_mat,1}(length(G))
  MS=MatrixSpace(F,ngens(V),ngens(V), false)
  for s=1:length(G1)
    G1[s]=MS(0)
    for i=1:ngens(V)
      for j=1:ngens(V)
        if G[s][i,j]!=0 && v[i] <= v[j]
          G1[s][i,j]=divexact((G[s][i,j].data)*v[i],v[j])
        end
      end
    end
  end
  return FqGModule(G1)

end

function subm_to_subg(M::ZpnGModule, S::nmod_mat; op=sub)

  return op(M.V,lift(S))

end

##########################################################################
#
#  Minimal Submodules function
#
##########################################################################

function minimal_submodules(M::ZpnGModule, ord::Int=-1)

  R=M.R
  S,mS=snf(M)
  N=_exponent_p_sub(S)
  if ord==-1
    list_sub=minimal_submodules(N)
  else
    list_sub=minimal_submodules(N, ngens(S.V)-ord)
  end
  list=Array{nmod_mat,1}(length(list_sub))
  v=[M.p^(valuation(S.V.snf[i], M.p)-1) for i=1:ngens(S.V)]
  W=MatrixSpace(R,1, ngens(M.V), false)
  for z=1:length(list)
    list[z]= vcat([W((mS(S.V([FlintZZ(coeff(list_sub[z][k,i],0))*v[i] for i=1:ngens(S.V)]))).coeff) for k=1:rows(list_sub[z])])
  end
  return list

end

###########################################################################################
#
#  Maximal Submodules function
#
###########################################################################################

function maximal_submodules(M::ZpnGModule, ind::Int=-1)

  R=M.R
  S,mS=snf(M)
  N=dual_module(S)
  if ind==-1
    minlist=minimal_submodules(N)
  else
    minlist=minimal_submodules(N,ind)
  end
  list=Array{nmod_mat,1}(length(minlist))
  W=MatrixSpace(R,1,ngens(S.V), false)
  v=[divexact(fmpz(R.n),S.V.snf[j]) for j=1:ngens(S.V) ]
  for x in minlist
    K=DiagonalGroup([fmpz(R.n) for j=1:rows(x)])
    A=lift(transpose(x))
    for j=1:rows(A)
      for k=1:cols(A)
        A[j,k]*=v[j]
      end
    end
    mH=Hecke.GrpAbFinGenMap(S.V,K,A)
    sg,msg=kernel(mH)
    push!(list, vcat([ (mS(msg(y))).coeff for y in gens(sg)]))
  end
  return list

end

##########################################################################################
#
#  Composition factors and series
#
##########################################################################################

#
#  Given a list of square matrices G, it returns a list of matrices given by the minors
#  (n-s) x (n-s) of the matrices G[i] mod p
#

function _change_ring(G::Array{nmod_mat,1}, F::Nemo.FqNmodFiniteField, s::Int)

  G1=Array{fq_nmod_mat,1}(length(G))
  n=rows(G[1])
  for i=1:length(G)
    M=zero_matrix(F,n-s+1,n-s+1)
    for j=s:n
      for k=s:n
        M[j-s+1,k-s+1]=(G[i][j,k]).data
      end
    end
    G1[i]=M
  end
  return G1

end

#
#  Cut the module in submodules with exponent p, returning the quotients p^i M /p^(i+1) M
#

function _mult_by_p(M::ZpnGModule)

  G=M.G
  V=M.V
  p=M.p
  @assert issnf(V)
  #
  #  First, the quotient M/pM
  #
  F,a=Nemo.FiniteField(p,1,"a", cached=false)
  n=ngens(V)
  Gq=_change_ring(G,F,1)
  spaces=FqGModule[FqGModule(Gq)]
  #
  #  Now, the others
  #
  s=valuation(fmpz(M.R.n),p)
  j=1
  for i=2:s
    while !divisible(V.snf[j],p^i)
      j+=1
    end
    GNew=_change_ring(G,F,j)
    push!(spaces, FqGModule(GNew))
  end
  return spaces

end

function composition_factors(M::ZpnGModule)

  N,_=snf(M)
  list_sub=_mult_by_p(N)
  list=[]
  for x in list_sub
    append!(list, composition_factors(x))
  end
  for i=1:length(list)
    j=i+1
    while j<=length(list)
      if Hecke.isisomorphic(list[i][1], list[j][1])
        list[i][2]+=list[j][2]
        deleteat!(list,j)
      else
        j+=1
      end
    end
  end
  return list

end

#######################################################################################
#
#  Submodules function
#
#######################################################################################

doc"""
    submodules(M::ZpnGModule; typequo, typesub, order)
> Given a ZpnGModule M, the function returns all the submodules of M.
>
"""

function submodules(M::ZpnGModule; typequo=Int[-1], typesub=Int[-1], ord=-1)

  if typequo!=[-1]
    return submodules_with_quo_struct(M,typequo)
  elseif typesub!=[-1]
    return submodules_with_struct(M,typesub)
  elseif ord!=-1
    return submodules_order(M,ord)
  else
    return submodules_all(M)
  end

end

function submodules_all(M::ZpnGModule)

  R=M.R
  S,mS=snf(M)
  minlist=minimal_submodules(S)
  list=nmod_mat[MatrixSpace(R,length(S.V.snf),length(S.V.snf), false)(1), MatrixSpace(R,1,length(S.V.snf), false)(0)]

  #
  #  Find minimal submodules, factor out and recursion on the quotients
  #

  for x in minlist
    N,_=quo(S,x)
    newlist=submodules_all(N)
    for y in newlist
      push!(list,vcat(y,x))
    end
  end

  append!(list,minlist)
  #
  #  Eliminate redundance via Howell form
  #
  w=fmpz[divexact(fmpz(R.n), S.V.snf[j]) for j=1:ngens(S.V)]
  list=_no_redundancy(list,w)

  #
  #  Writing the submodules in terms of the given generators and returning an iterator
  #
  W=MatrixSpace(R,rows(mS.map), cols(mS.map), false)
  MatSnf=W(mS.map)
  return (x*MatSnf for x in list)

end

function submodules_with_struct_cyclic(M::ZpnGModule, ord::Int)

  R=M.R
  #
  #  We search for an element in p^(ord-1)*G
  #
  s,ms=sub(M, M.p^(ord-1))
  S,mS=snf(s)
  N=_exponent_p_sub(S)
  submod=fq_nmod_mat[]
  if N.dim==1
    push!(submod, identity_matrix(N.K, 1))
  else
    @vtime :StabSub 1 submod=minimal_submodules(N,1,composition_factors(N))
  end
  list1=Array{nmod_mat,1}(length(submod))
  v=fmpz[(M.p)^(valuation(S.V.snf[i], M.p)-1) for i=1:ngens(S.V)]
  for i=1:length(submod)
    @views list1[i]=lift(submod[i],R)
    @assert rows(list1[i])==1
    for k=1:cols(list1[i])
      @inbounds list1[i][1,k]*=v[k]
    end
  end
  W=MatrixSpace(R,rows(mS.map), cols(ms.map), false)
  MatSnf=W(mS.map*ms.map)
  for j=1:length(list1)
    @inbounds @views list1[j]=list1[j]*MatSnf
  end
  if ord==1
    return list1
  end
  list=nmod_mat[]
  for x in list1
    L,_=quo(M,x)
    newlist=submodules_with_struct_cyclic(L,ord-1)
    i=1
    while i<=length(newlist)
      @views if iszero(M.V(lift(M.p^(ord-1)*newlist[i])))
        deleteat!(newlist,i)
      else
        i+=1
      end
    end
    append!(list, newlist)
  end
  return list

end

function submodule_with_struct_exp_p(M::ZpnGModule, l::Int)

    R=M.R
    S,mS=snf(M)
    N=Hecke._exponent_p_sub(S)
    lf=composition_factors(N)
    list=nmod_mat[]
    v=fmpz[divexact(S.V.snf[i], M.p) for i=1:ngens(S.V)]
    submod=submodules(N,ngens(S.V)-l,comp_factors=lf)
    list1=Array{nmod_mat,1}(length(submod))
    for i=1:length(submod)
      list1[i]=lift(submod[i],R)
      for s=1:rows(list1[i])
        for t=1:cols(list1[i])
          @inbounds list1[i][s,t]*=v[t]
        end
      end
    end
    W=MatrixSpace(R,rows(mS.map), cols(mS.map), false)
    MatSnf=W(mS.map)
    return ( x *MatSnf for x in list1)

end

function submodules_with_struct(M::ZpnGModule, typesub::Array{Int,1})

  # If the group is cyclic, it is easier
  if length(typesub)==1
    return submodules_with_struct_cyclic(M,typesub[1])
  end
  sort!(typesub)
  # If the subgroups we are searching for have exponent p, it is easier
  if typesub[end]==1
    return submodule_with_struct_exp_p(M,length(typesub))
  end
  R=M.R

  a=1
  while a<length(typesub) && typesub[end-a]==typesub[end]
    a+=1
  end

  S1,mS1=snf(M)
  s,ms=sub(S1, M.p^(typesub[end]-1))
  S,mS=snf(s)
  N=_exponent_p_sub(S)
  @vtime :StabSub 1 submod=submodules(N,(N.dim)-a)
  #
  #  Write the modules as elements of S
  #
  list1=Array{nmod_mat,1}(length(submod))
  v=[divexact(S.V.snf[i], M.p) for i=1:ngens(S.V)]
  for i=1:length(submod)
    @inbounds list1[i]=lift(submod[i],R)
    for j=1:rows(list1[i])
      for k=1:cols(list1[i])
        @inbounds list1[i][j,k]*=v[k]
      end
    end
  end

  auxmat=mS.imap*ms.map
  W=MatrixSpace(R,rows(auxmat), cols(auxmat), false)
  auxmat=W(auxmat)
  for j=1:length(list1)
    @inbounds list1[j]=list1[j]*auxmat
  end
  #
  #  I create the group to check if the candidates are isomorphic to it
  #

  l=[M.p^(typesub[j]) for j=1:length(typesub)]
  G=DiagonalGroup(l)
  G.issnf=true
  G.snf=l

  #
  #  I create the type of the group I am searching for in the quotient
  #
  new_subtype=deepcopy(typesub)
  for j=length(typesub)-a+1:length(typesub)
    @inbounds new_subtype[j]-=1
  end

  #
  #  Recursion on the quotient
  #

  list=nmod_mat[]
  for x in list1
    L, _=quo(S1,x)
    newlist=collect(submodules_with_struct(L,new_subtype))
    for j=1:length(newlist)
      @inbounds newlist[j]=vcat(newlist[j],x)
    end
    i=1
    while i<=length(newlist)
      t,mt=subm_to_subg(S1,newlist[i])
      if !isisomorphic(t,G)
        deleteat!(newlist,i)
      else
        i+=1
      end
    end
    if !isempty(newlist)
      append!(list, newlist)
    end
  end
  #
  #  Check for redundancy
  #
  w=fmpz[divexact(fmpz(R.n), S1.V.snf[j]) for j=1:ngens(S1.V)]
  list=_no_redundancy(list,w)

  #
  #  Write the submodules in terms of the set of given generators
  #  and return an iterator over the list
  #

  W=MatrixSpace(R,ngens(S1.V), ngens(M.V), false)
  MatSnf=W(mS1.map)

  return (el*MatSnf for el in list)

end

function _no_redundancy(list::Array{nmod_mat,1},w::Array{fmpz,1})

  #
  #  Howell form of every candidate, embedding them in a free module
  #
  for i=1:length(list)
    if cols(list[i])>=rows(list[i])
      @inbounds list[i]=vcat(list[i],zero_matrix(parent(list[i][1,1]),cols(list[i])-rows(list[i]), cols(list[i])))
    end
    for j=1:cols(list[i])
      for k=1:rows(list[i])
        @inbounds list[i][k,j]*=w[j]
      end
    end
    howell_form!(list[i])
    list[i]=sub(list[i],1:cols(list[i]),1:cols(list[i]))
  end
  #
  #  Now, check if they are equal
  #
  i=1
  while i<length(list)
    k=i+1
    while k<=length(list)
      if list[i]==list[k]
        deleteat!(list,k)
      else
        k+=1
      end
    end
    i+=1
  end

  #
  #  Write them again not embedded
  #
  for i=1:length(list)
    for j=1:cols(list[i])
      for k=1:j
        list[i][k,j]=divexact(list[i][k,j].data,w[j])
      end
    end
  end
  return list

end

function submodules_order(M::ZpnGModule, ord::Int)

  R=M.R
  S,mS=snf(M)
  @assert exponent(S.V)==fmpz(R.n)
  N=Hecke._exponent_p_sub(S)
  lf=composition_factors(N)
  list=nmod_mat[]
  v=fmpz[divexact(S.V.snf[i], M.p) for i=1:ngens(S.V)]
  for i=1:ord-1
    minlist=minimal_submodules(N,i,lf)
    for x in minlist
      A=lift(x,R)
      for s=1:rows(A)
        for t=1:cols(A)
          A[s,t]*=v[t]
        end
      end
      L, _=quo(S,A)
      newlist=Hecke.submodules_order(L,ord-i)
      for z=1:length(newlist)
        push!(list,vcat(newlist[z],A))
      end
    end
  end

  #
  #  Check for redundancy
  #
  w=fmpz[divexact(fmpz(R.n), S.V.snf[j]) for j=1:ngens(S.V)]
  list=_no_redundancy(list,w)

  #
  #  Write the submodules in terms of the set of given generators
  #

  W=MatrixSpace(R,rows(mS.map), cols(mS.map), false)
  MatSnf=W(mS.map)
  for j=1:length(list)
    list[j]=list[j]*MatSnf #vcat([W(( mS( S.V([list[j][k,i].data for i=1:ngens(S.V)]))).coeff)  for k=1:rows(list[j])])
  end

  #
  #  Add the minimal submodules
  #

  minlist=minimal_submodules(N,ord, lf)
  for x in minlist
    push!(list, vcat([W((mS( S.V([FlintZZ(coeff(x[k,i],0))*((M.p)^(v[i]-1)) for i=1:ngens(S.V)]))).coeff) for k=1:rows(x) ]))
  end
  return (x for x in list)

end

function submodules_with_quo_struct(M::ZpnGModule, typequo::Array{Int,1})

  R=M.R
  S,mS=snf(M)
  sort!(typequo)
  l=[(M.p)^typequo[i] for i=1:length(typequo)]
  wish=DiagonalGroup(l)
  wish.issnf=true
  wish.snf=l
  if isisomorphic(wish,S.V)
    return nmod_mat[zero_matrix(R, 1, ngens(M.V))]
  end
  if length(l)>length(S.V.snf)
    return nmod_mat[]
  end
  for i=1:length(typequo)
    if !divisible(S.V.snf[ngens(S.V)+1-i],fmpz((M.p)^typequo[length(typequo)+1-i]))
      return nmod_mat[]
    end
  end

  #
  #  Dual Module and candidate submodules
  #
  M_dual=dual_module(S)
  candidates=submodules_with_struct(M_dual,typequo)

  #
  #  Dualize the modules
  #
  v=[divexact(S.V.snf[end],S.V.snf[j]) for j=1:ngens(S.V) ]
  list=(_dualize(x, S.V, v) for x in candidates)
  #
  #  Write the submodules in terms of the given generators
  #
  W=MatrixSpace(R,rows(mS.map), cols(mS.map), false)
  MatSnf=W(mS.map)
  return (x*MatSnf for x in list)

end


text:
# sequences.jl

"""
A Sequence is a bi-infinite discrete set of values, indexed by all integers.

Each Sequence has an element type, given by `eltype`.
"""
abstract Sequence{T}

Base.eltype{T}(::Type{Sequence{T}}) = T
Base.eltype{S <: Sequence}(::Type{S}) = eltype(supertype(S))

# Return a range of the sequence as a vector
Base.getindex(s::Sequence, r::Range) = [s[k] for k in range]

iscompact(s::Sequence) = false

# We adopt the convention that the transpose of a sequence corresponds to time-reversal
Base.transpose(s::Sequence) = reverse(s)

# Hence, ctranspose is time-reversal plus taking conjugates
Base.ctranspose(s::Sequence) = conj(reverse(s))


text:
using FaultTolerantControl
using UnPack, ComponentArrays
using Transducers
using Plots

function make_env()
    multicopter = IslamQuadcopterEnv()
    Λ_func_compat = (x, p, t) -> ones(4) |> Diagonal |> Matrix  # effectiveness matrix
    @unpack m, g, B = multicopter
    allocator = PseudoInverseAllocator(B)
    x0_multicopter = State(multicopter)()
    pos0 = copy(x0_multicopter.p)
    vel0 = copy(x0_multicopter.v)
    helixCG = HelixCommandGenerator(pos0)
    cg = Command(helixCG)
    controller = BacksteppingPositionControllerEnv(m; pos_cmd_func=cg)
    x0_controller = State(controller)(pos0, m, g)
    x0 = ComponentArray(multicopter=x0_multicopter, controller=x0_controller)
    multicopter, controller, allocator, x0, cg, Λ_func_compat
end

function main()
    multicopter, controller, allocator, x0, cg, Λ_func_compat = make_env()
    prob, df = sim(
                    x0,
                    apply_inputs(Dynamics!(multicopter, controller, allocator); Λ=Λ_func_compat);
                    tf=40.0,
                   )
    t0, tf = prob.tspan
    Δt = 0.01  # data sampling period; not simulation time step
    ts = t0:Δt:tf
    poss = ts |> Map(t -> sol(t).multicopter.p) |> collect
    poss_ref = ts |> Map(t -> cg(t)) |> collect
    ## plot
    # 3d traj
    p_traj = plot3d(; title="position", legend=:outertopright)
    plot!(p_traj, hcat(poss...)'[:, 1], hcat(poss...)'[:, 2], hcat(poss...)'[:, 3]; label="position", color="red")
    plot!(p_traj, hcat(poss_ref...)'[:, 1], hcat(poss_ref...)'[:, 2], hcat(poss_ref...)'[:, 3]; label="position (ref)", color="black")
    # save
    dir_log = "figures"
    mkpath(dir_log)
    savefig(p_traj, joinpath(dir_log, "traj_multicopter_backstepping.png"))
end


text:
# Sampling from manifolds. Currently used for testing only.
abstract type Manifold{Int, T<:AbstractFloat} end

# ............................................................................ #
struct Torus{T} <: Manifold{3, T}
    R::T
    r::T
    center::Tuple{T, T, T}
end

function Torus(R, r)
    args = convert.(AbstractFloat, (R, r))
    z = zero(eltype(args))
    Torus(args..., (z, z, z))
end

(t::Torus)(θ, φ) =
    ((t.R + t.r*cos(θ)) * cos(φ), (t.R + t.r*cos(θ)) * sin(φ), t.r * sin(θ)) .+
    t.center

# ............................................................................ #
struct AsymTorus{T} <: Manifold{3, T}
    R::T
    r_large::T
    r_small::T
    center::Tuple{T, T, T}
end

function AsymTorus(R, r1, r2)
    args = convert.(AbstractFloat, (R, r1, r2))
    z = zero(eltype(args))
    AsymTorus(args..., (z, z, z))
end

function (at::AsymTorus)(θ, φ)
    R = at.R
    r_large = at.r_large
    r_small = at.r_small

    r = (r_large - r_small) * (cos(φ) + 1)/2 + r_small
    ((R + r*cos(θ)) * cos(φ), (R + r*cos(θ)) * sin(φ), r*sin(θ)) .+ at.center
end

# ............................................................................ #
struct Sphere{T} <: Manifold{3, T}
    r::T
    center::Tuple{T, T, T}
end

function Sphere(r)
    arg = convert(AbstractFloat, r)
    z = zero(typeof(arg))
    Sphere(arg, (z, z, z))
end

(s::Sphere)(θ, φ) =
    (s.r*cos(θ) * sin(φ), s.r*sin(θ) * sin(φ), s.r*cos(φ)) .+ s.center

# ............................................................................ #
struct Circle{T} <: Manifold{2, T}
    r::T
    center::Tuple{T, T}
end

function Circle(r)
    arg = convert(AbstractFloat, r)
    z = zero(typeof(arg))
    Circle(arg, (z, z))
end

(c::Circle)(θ) = (c.r * cos(θ), c.r * sin(θ)) .+ c.center

# ............................................................................ #
# Uniformly sampling n points from a manifold returning D by n matrix.
# TODO: Circular noise? Overenginnered as it is?
function Base.rand(m::Manifold{D, T}, n=1; noise=0) where {D, T}
    res = zeros(T, D, n)
    for i in 1:n
        pars = 2π .* rand(D-1)
        res[:, i] .= m(pars...) .+ (rand(D) - 0.5) * noise
    end
    res
end


text:
# minimal use of the package, user manages all the plotting
# PyPlot version

module M
using Pseudospectra
using PyPlot

A = randn(100,100)
# we want a true Schur decomposition, so force complexity
Tschur,U,eigA  = schur(A .+ 0im)
ax = [-12,12,-12,12]
npts = 100
opts = Dict{Symbol,Any}(:real_matrix => true)
Z,x,y,levels,err,Tproj,eigAproj,algo = psa_compute(Tschur,npts,
                                                   ax,eigA,opts)
plot(real(eigA),imag(eigA),"k.")
if isempty(levels)
    contour(x,y,log10.(Z),extend="both",linewidths=2)
else
    contour(x,y,log10.(Z),levels=levels,extend="both",linewidths=2)
end
colorbar()
end


text:
import J4Org: extract_tag, find_tag, tags, link, line, tokenized

@testset "extract_tag" begin

    r=extract_tag(tokenized("one_word"),1)
    @test r==nothing

    r=extract_tag(tokenized("#one_word"),1)
    @test r==nothing

    r=extract_tag(tokenized("#+one_word"),1)
    @test r!=nothing
    @test length(tags(r))==1
    @test line(r) == 1

    r=extract_tag(tokenized("#+one_word, two_Words L:my_link   "),1)
    @test r!=nothing
    @test length(tags(r))==2
    @test tags(r)[1] == "one_word"
    @test tags(r)[2] == "two_Words"
    @test link(r) == "my_link"
    @test line(r) == 1

    # also support "# +"
    r=extract_tag(tokenized("# +one_word, two_Words L:my_link   "),1)
    @test r!=nothing
    @test length(tags(r))==2
    @test tags(r)[1] == "one_word"
    @test tags(r)[2] == "two_Words"
    @test link(r) == "my_link"
    @test line(r) == 1

end;

@testset "extract_tag" begin

    filename = "$(dirname(@__FILE__))/code_examples/basic.jl"
    code=readcode(filename)
    t=tokenized(code)

    @test extract_tag(t,1)==nothing
    @test extract_tag(t,2)==nothing
    @test extract_tag(t,3)!=nothing
end;

@testset "find_tag" begin

    filename = "$(dirname(@__FILE__))/code_examples/basic.jl"

    code=readcode(filename)
    t=tokenized(code)

    r = find_tag(t,1)
    @test line(r) == 3
    @test tags(r) == ["MyTag1";"MyTag2"]

    r = find_tag(t,skip(r))
    @test r==nothing
end;


text:
t1 = 1
t2 = 2
sum = 2

while t1 + t2 < 4_000_000

    if (t1 + t2) % 2 == 0
        global sum += (t1+t2)
    else
    end
    global t1, t2 = t2, t1 + t2
end

println(sum)


text:
#test-warlpiri-addition.jl

###############################################################################
##COMPLETE ORDER TABLE

#column key

#none_ pfew_ pone_ psome ptwo_ pmany pall_ junk_ nnone n_few n_one nsome n_two nmany n_all
wcmp = Bool[
#none_ > line
 false false false false false false false false false true  true  true  true  true  true;
#pfew_ > line
 true  false false false false false false false true  true  true  true  true  true  true;
#pone_ > line
 true  true  false false false false false false true  true  true  true  true  true  true;
#psome > line
 true  true  true  false false false false false true  true  true  true  true  true  true;
#ptwo_ > line
 true  true  true  true  false false false false true  true  true  true  true  true  true;
#pmany > line
 true  true  true  true  true  false false false true  true  true  true  true  true  true;
#pall_ > line
 true  true  true  true  true  true  false false true  true  true  true  true  true  true;
#junk_ > line
 false false false false false false false false false false false false false false false;
#nnone > line
 false false false false false false false false false true  true  true  true  true  true;
#n_few > line
 false false false false false false false false false false true  true  true  true  true;
#n_one > line
 false false false false false false false false false false false true  true  true  true;
#nsome > line
 false false false false false false false false false false false false true  true  true;
#n_two > line
 false false false false false false false false false false false false false true  true;
#nmany > line
 false false false false false false false false false false false false false false true;
#n_all > line
 false false false false false false false false false false false false false false false;
]

testop(>, wcmp)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
@testset "lightweight gen function" begin

##########
# update #
##########

@testset "force update" begin

    @gen function bar()
        @addr(normal(0, 1), :a)
    end

    @gen function baz()
        @addr(normal(0, 1), :b)
    end

    @gen function foo()
        if @addr(bernoulli(0.4), :branch)
            @addr(normal(0, 1), :x)
            @addr(bar(), :u)
        else
            @addr(normal(0, 1), :y)
            @addr(baz(), :v)
        end
    end

    # get a trace which follows the first branch
    constraints = DynamicAssignment()
    constraints[:branch] = true
    (trace,) = generate(foo, (), constraints)
    x = get_assignment(trace)[:x]
    a = get_assignment(trace)[:u => :a]

    # force to follow the second branch
    y = 1.123
    b = -2.1
    constraints = DynamicAssignment()
    constraints[:branch] = false
    constraints[:y] = y
    constraints[:v => :b] = b
    (new_trace, weight, discard, retchange) = update(
        foo, (), nothing, trace, constraints)

    # test discard
    @test get_leaf_node(discard, :branch) == true
    @test get_leaf_node(discard, :x) == x
    @test get_leaf_node(discard, :u => :a) == a
    @test length(collect(get_leaf_nodes(discard))) == 2
    @test length(collect(get_internal_nodes(discard))) == 1

    # test new trace
    new_assignment = get_assignment(new_trace)
    @test get_leaf_node(new_assignment, :branch) == false
    @test get_leaf_node(new_assignment, :y) == y
    @test get_leaf_node(new_assignment, :v => :b) == b
    @test length(collect(get_leaf_nodes(new_assignment))) == 2
    @test length(collect(get_internal_nodes(new_assignment))) == 1

    # test score and weight
    prev_score = (
        logpdf(bernoulli, true, 0.4) +
        logpdf(normal, x, 0, 1) +
        logpdf(normal, a, 0, 1))
    expected_new_score = (
        logpdf(bernoulli, false, 0.4) +
        logpdf(normal, y, 0, 1) +
        logpdf(normal, b, 0, 1))
    expected_weight = expected_new_score - prev_score
    @test isapprox(expected_new_score, get_call_record(new_trace).score)
    @test isapprox(expected_weight, weight)

    # test retchange (should be nothing by default)
    @test retchange === nothing
end

##############
# fix_update #
##############

@testset "fix update" begin

    @gen function bar()
        @addr(normal(0, 1), :a)
    end

    @gen function baz()
        @addr(normal(0, 1), :b)
    end

    @gen function foo()
        if @addr(bernoulli(0.4), :branch)
            @addr(normal(0, 1), :x)
            @addr(bar(), :u)
        else
            @addr(normal(0, 1), :y)
            @addr(baz(), :v)
        end
        @addr(normal(0, 1), :z)
    end

    # get a trace which follows the first branch
    constraints = DynamicAssignment()
    constraints[:branch] = true
    (trace,) = generate(foo, (), constraints)
    x = get_assignment(trace)[:x]
    a = get_assignment(trace)[:u => :a]
    z = get_assignment(trace)[:z]

    # force to follow the second branch, and change z
    y = 1.123
    b = -2.1
    z_new = 0.4
    constraints = DynamicAssignment()
    constraints[:branch] = false
    constraints[:z] = z_new
    (new_trace, weight, discard, retchange) = fix_update(
        foo, (), nothing, trace, constraints)

    # test discard
    @test get_leaf_node(discard, :branch) == true
    @test get_leaf_node(discard, :z) == z
    @test length(collect(get_leaf_nodes(discard))) == 2
    @test length(collect(get_internal_nodes(discard))) == 0

    # test new trace
    new_assignment = get_assignment(new_trace)
    @test get_leaf_node(new_assignment, :branch) == false
    @test get_leaf_node(new_assignment, :z) == z_new
    y = get_leaf_node(new_assignment, :y)
    b = get_leaf_node(new_assignment, :v => :b)
    @test length(collect(get_leaf_nodes(new_assignment))) == 3
    @test length(collect(get_internal_nodes(new_assignment))) == 1

    # test score and weight
    expected_new_score = (
        logpdf(bernoulli, false, 0.4) +
        logpdf(normal, y, 0, 1) +
        logpdf(normal, b, 0, 1) +
        logpdf(normal, z_new, 0, 1))
    expected_weight = (
        logpdf(bernoulli, false, 0.4)
        - logpdf(bernoulli, true, 0.4)
        + logpdf(normal, z_new, 0, 1)
        - logpdf(normal, z, 0, 1))
    @test isapprox(expected_new_score, get_call_record(new_trace).score)
    @test isapprox(expected_weight, weight)

    # test retchange (should be nothing by default)
    @test retchange === nothing
end

###############
# free_update #
###############

@testset "regenerate" begin

    @gen function bar(mu)
        @addr(normal(mu, 1), :a)
    end

    @gen function baz(mu)
        @addr(normal(mu, 1), :b)
    end

    @gen function foo(mu)
        if @addr(bernoulli(0.4), :branch)
            @addr(normal(mu, 1), :x)
            @addr(bar(mu), :u)
        else
            @addr(normal(mu, 1), :y)
            @addr(baz(mu), :v)
        end
    end

    # get a trace which follows the first branch
    mu = 0.123
    constraints = DynamicAssignment()
    constraints[:branch] = true
    (trace,) = generate(foo, (mu,), constraints)
    x = get_assignment(trace)[:x]
    a = get_assignment(trace)[:u => :a]

    # resimulate branch
    selection = DynamicAddressSet()
    push_leaf_node!(selection, :branch)

    # try 10 times, so we are likely to get both a stay and a switch
    for i=1:10
        prev_assignment = get_assignment(trace)

        # change the argument so that the weights can be nonzer
        prev_mu = mu
        mu = rand()
        (trace, weight, retchange) = regenerate(foo, (mu,), nothing, trace, selection)
        assignment = get_assignment(trace)

        # test score
        if assignment[:branch]
            expected_score = logpdf(normal, assignment[:x], mu, 1) + logpdf(normal, assignment[:u => :a], mu, 1) + logpdf(bernoulli, true, 0.4)
        else
            expected_score = logpdf(normal, assignment[:y], mu, 1) + logpdf(normal, assignment[:v => :b], mu, 1) + logpdf(bernoulli, false, 0.4)
        end
        @test isapprox(expected_score, get_call_record(trace).score)

        # test values
        if assignment[:branch]
            @test has_leaf_node(assignment, :x)
            @test has_internal_node(assignment, :u)
        else
            @test has_leaf_node(assignment, :y)
            @test has_internal_node(assignment, :v)
        end
        @test length(collect(get_leaf_nodes(assignment))) == 2
        @test length(collect(get_internal_nodes(assignment))) == 1

        # test weight
        if assignment[:branch] == prev_assignment[:branch]
            if assignment[:branch]
                expected_weight = logpdf(normal, assignment[:x], mu, 1) + logpdf(normal, assignment[:u => :a], mu, 1)
                expected_weight -= logpdf(normal, assignment[:x], prev_mu, 1) + logpdf(normal, assignment[:u => :a], prev_mu, 1)
            else
                expected_weight = logpdf(normal, assignment[:y], mu, 1) + logpdf(normal, assignment[:v => :b], mu, 1)
                expected_weight -= logpdf(normal, assignment[:y], prev_mu, 1) + logpdf(normal, assignment[:v => :b], prev_mu, 1)
            end
        else
            expected_weight = 0.
        end
        @test isapprox(expected_weight, weight)

        # test retchange (should be nothing by default)
        @test retchange === nothing
    end

end

end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Logging
using JuMP, Gurobi
using DecisionProgramming
using CSV, DataFrames, PrettyTables

# Setting subproblem specific parameters
const chosen_risk_level = "12%"

# Reading tests' technical performance data (dummy data in this case)
data = CSV.read("CHD_preventative_care_data.csv", DataFrame)

# Bayes posterior risk probabilities calculation function
# prior = prior risk level for which the posterior risk distribution is calculated for,
# t = test done
# returns a 100x1 vector with the probabilities of getting CHD given the prior risk level and test result
# for no test done (i.e. t = 3) returns a zero vector
function update_risk_distribution(prior::Int64, t::Int64)
    if t == 1 # the test is TRS
        # P(TRS = result | sick) = P(TRS_if_sick = result) * P(sick) = P(TRS_if_sick = result) * P(prior_risk)
        numerators = data.TRS_if_sick .* data.risk_levels[prior]

        # P(TRS = result) = P(TRS_if_sick = result) * P(sick) + P(TRS_if_healthy = result) * P(healthy)
        denominators = data.TRS_if_sick .* data.risk_levels[prior]  + data.TRS_if_healthy .* (1 - data.risk_levels[prior])

        posterior_risks = numerators./denominators

        # if the denominator is zero, post_risk is NaN, changing those to 0
        for i = 1:101
            if isnan(posterior_risks[i])
                posterior_risks[i] = 0
            end
        end

        return posterior_risks

    elseif t == 2 #the test is GRS
        numerators = (data.GRS_if_sick .* data.risk_levels[prior])
        denominators = data.GRS_if_sick .* data.risk_levels[prior]  + data.GRS_if_healthy .* (1 .- data.risk_levels[prior])

        posterior_risks =  numerators./denominators

        # if the denominator is zero, post_risk is NaN, changing those to 0
        for i = 1:101
            if isnan(posterior_risks[i])
                posterior_risks[i] = 0
            end
        end

        return posterior_risks

    else # no test performed
        risks_unchanged = zeros(100,1)

        return risks_unchanged

    end
end

# State probabilites calculation function
# risk_p = the resulting array from update_risk_distribution
# t = test done
# h = CHD or no CHD
# returns the probability distribution in 101x1 vector for the states of the R node given the prior risk level (must be same as to function update_risk_distribution), test t and health h
function state_probabilities(risk_p::Array{Float64}, t::Int64, h::Int64, prior::Int64)

    #if no test is performed, then the probabilities of moving to states (other than the prior risk level) are 0 and to the prior risk element is 1
    if t == 3
        state_probabilites = zeros(101)
        state_probabilites[prior] = 1.0
        return state_probabilites
    end

    # return vector
    state_probabilites = zeros(101)

    # copying the probabilities of the scores for ease of readability
    if h == 1 && t == 1    # CHD and TRS
        p_scores = data.TRS_if_sick
    elseif t ==1    # no CHD and TRS
        p_scores = data.TRS_if_healthy
    elseif h == 1 && t == 2 # CHD and GRS
        p_scores = data.GRS_if_sick
    else # no CHD and GRS
        p_scores = data.GRS_if_healthy
    end

    for i = 1:101 #iterating through all risk levels 0%, 1%, ..., 99%, 100% in data.risk_levels
        for j = 1:101 #iterates through all risk estimates in risk_p
            #finding all risk estimates risk_p[j] within risk level i
            # risk_level[i] <= risk_p < risk_level[i]
            if i < 101 && data.risk_levels[i] <= risk_p[j] && risk_p[j] < data.risk_levels[i+1]
                state_probabilites[i] += p_scores[j]
            elseif i == 101 && data.risk_levels[i] <= risk_p[j] #special case: the highest risk level[101] = 100%
                state_probabilites[i] += p_scores[j]
            end
        end
    end

    return state_probabilites
end

@info("Creating the influence diagram.")
diagram = InfluenceDiagram()

const H_states = ["CHD", "no CHD"]
const T_states = ["TRS", "GRS", "no test"]
const TD_states = ["treatment", "no treatment"]
const R_states = [string(x) * "%" for x in [0:1:100;]]

add_node!(diagram, ChanceNode("R0", [], R_states))
add_node!(diagram, ChanceNode("R1", ["R0", "H", "T1"], R_states))
add_node!(diagram, ChanceNode("R2", ["R1", "H", "T2"], R_states))
add_node!(diagram, ChanceNode("H", ["R0"], H_states))

add_node!(diagram, DecisionNode("T1", ["R0"], T_states))
add_node!(diagram, DecisionNode("T2", ["R1"], T_states))
add_node!(diagram, DecisionNode("TD", ["R2"], TD_states))

add_node!(diagram, ValueNode("TC", ["T1", "T2"]))
add_node!(diagram, ValueNode("HB", ["H", "TD"]))

generate_arcs!(diagram)

X_R0 = ProbabilityMatrix(diagram, "R0")
X_R0[chosen_risk_level] = 1
add_probabilities!(diagram, "R0", X_R0)

X_H = ProbabilityMatrix(diagram, "H")
X_H[:, "CHD"] = data.risk_levels
X_H[:, "no CHD"] = 1 .- data.risk_levels
add_probabilities!(diagram, "H", X_H)

X_R = ProbabilityMatrix(diagram, "R1")
for s_R0 = 1:101, s_H = 1:2, s_T1 = 1:3
    X_R[s_R0, s_H,  s_T1, :] =  state_probabilities(update_risk_distribution(s_R0, s_T1), s_T1, s_H, s_R0)
end
add_probabilities!(diagram, "R1", X_R)
add_probabilities!(diagram, "R2", X_R)

cost_TRS = -0.0034645
cost_GRS = -0.004
forbidden = 0     #the cost of forbidden test combinations is negligible
Y_TC = UtilityMatrix(diagram, "TC")
Y_TC["TRS", "TRS"] = forbidden
Y_TC["TRS", "GRS"] = cost_TRS + cost_GRS
Y_TC["TRS", "no test"] = cost_TRS
Y_TC["GRS", "TRS"] = cost_TRS + cost_GRS
Y_TC["GRS", "GRS"] = forbidden
Y_TC["GRS", "no test"] = cost_GRS
Y_TC["no test", "TRS"] = cost_TRS
Y_TC["no test", "GRS"] = cost_GRS
Y_TC["no test", "no test"] = 0
add_utilities!(diagram, "TC", Y_TC)

Y_HB = UtilityMatrix(diagram, "HB")
Y_HB["CHD", "treatment"] = 6.89713671259061
Y_HB["CHD", "no treatment"] = 6.65436854256236
Y_HB["no CHD", "treatment"] = 7.64528451705134
Y_HB["no CHD", "no treatment"] = 7.70088349200034
add_utilities!(diagram, "HB", Y_HB)

generate_diagram!(diagram)

@info("Creating the decision model.")
model = Model()
z = DecisionVariables(model, diagram)

# Defining forbidden paths to include all those where a test is repeated twice
forbidden_tests = ForbiddenPath(diagram, ["T1","T2"], [("TRS", "TRS"),("GRS", "GRS"),("no test", "TRS"), ("no test", "GRS")])
fixed_R0 = FixedPath(diagram, Dict("R0" => chosen_risk_level))
scale_factor = 10000.0
x_s = PathCompatibilityVariables(model, diagram, z; fixed = fixed_R0, forbidden_paths = [forbidden_tests], probability_cut=false)

EV = expected_value(model, diagram, x_s, probability_scale_factor = scale_factor)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "MIPFocus" => 3,
    "MIPGap" => 1e-6,
)
set_optimizer(model, optimizer)

optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)
S_probabilities = StateProbabilities(diagram, Z)
U_distribution = UtilityDistribution(diagram, Z)

@info("Printing decision strategy using tailor made function:")
print_decision_strategy(diagram, Z, S_probabilities)

@info("Printing state probabilities:")
# Here we can see that the probability of having a CHD event is exactly that of the chosen risk level
print_state_probabilities(diagram, S_probabilities, ["R0", "R1", "R2"])

@info("Printing utility distribution.")
print_utility_distribution(U_distribution)

print_statistics(U_distribution)


text:
"""
    calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm, t)

Calculate element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
"""
@inline @muladd function calculate_residuals(ũ::Number, u₀::Number, u₁::Number,
                                    α, ρ, internalnorm,t)
    ũ / (α + max(internalnorm(u₀,t), internalnorm(u₁,t)) * ρ)
end

@inline function calculate_residuals(ũ::Array{T}, u₀::Array{T}, u₁::Array{T}, α::T2,
                                             ρ::Real, internalnorm,t) where
                                             {T<:Number,T2<:Number}
    out = similar(ũ)
    calculate_residuals!(out, ũ, u₀, u₁, α, ρ, internalnorm,t)
    out
end

@inline function calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm,t)
  @.. calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm,t)
end

"""
    calculate_residuals(u₀, u₁, α, ρ, internalnorm, t)

Calculate element-wise residuals
```math
\\frac{u₁ - u₀}{α+\\max{|u₀|,|u₁|}*ρ}
```
"""

@inline @muladd function calculate_residuals(u₀::Number, u₁::Number,
                                    α, ρ, internalnorm,t)
    (u₁ - u₀) / (α + max(internalnorm(u₀,t), internalnorm(u₁,t)) * ρ)
end

@inline function calculate_residuals(u₀::Array{T}, u₁::Array{T}, α::T2,
                                             ρ::Real, internalnorm,t) where
                                             {T<:Number,T2<:Number}
    out = similar(u₀)
    calculate_residuals!(out, u₀, u₁, α, ρ, internalnorm,t)
    out
end

@inline function calculate_residuals(u₀, u₁, α, ρ, internalnorm,t)
  @.. calculate_residuals(u₀, u₁, α, ρ, internalnorm,t)
end

"""
    calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)

Return element-wise residuals
```math
\\frac{δ E₁ + E₂}{α+\\max{scalarnorm(u₀),scalarnorm(u₁)}*ρ}.
```
"""
@inline @muladd function calculate_residuals(E₁::Number, E₂::Number, u₀::Number, u₁::Number,
                                             α::Real, ρ::Real, δ::Number, scalarnorm, t)
    (δ * E₁ + E₂) / (α + max(scalarnorm(u₀,t), scalarnorm(u₁,t)) * ρ)
end

@inline function calculate_residuals(E₁::Array{<:Number}, E₂::Array{<:Number},
                                     u₀::Array{<:Number}, u₁::Array{<:Number}, α::Real,
                                     ρ::Real, δ::Number, scalarnorm, t)
    out = similar(u₀)
    calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)
    out
end

@inline function calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)
  @.. calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)
end

# Inplace Versions

"""
    DiffEqBase.calculate_residuals!(out, ũ, u₀, u₁, α, ρ)

Save element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
in `out`.
"""
@inline function calculate_residuals!(out, ũ, u₀, u₁, α, ρ, internalnorm,t)
  @.. out = calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm,t)
  nothing
end

"""
    calculate_residuals!(out, u₀, u₁, α, ρ)

Save element-wise residuals
```math
\\frac{u₁ - u₀}{α+\\max{|u₀|,|u₁|}*ρ}
```
in `out`.
"""
@inline function calculate_residuals!(out, u₀, u₁, α, ρ, internalnorm,t)
  @.. out = calculate_residuals(u₀, u₁, α, ρ, internalnorm,t)
end

"""
    calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm)

Calculate element-wise residuals
```math
\\frac{δ E₁ + E₂}{α+\\max{scalarnorm(u₀),scalarnorm(u₁)}*ρ}.
```
"""
@inline function calculate_residuals!(out, E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)
  @.. out = calculate_residuals(E₁, E₂, u₀, u₁, α, ρ, δ, scalarnorm, t)
  out
end


text:
using MAT
using PyPlot
using ADCME
using DelimitedFiles

data = matread("data.mat")["data"]

kid = 3
tid = 1

if length(ARGS)==2
    global kid = parse(Int64, ARGS[1])
    global tid = parse(Int64, ARGS[2])
end

@info kid, tid

d = data[kid]
n = div(size(d, 1), 5)

H = Variable(1.0)^2
ε = d[1:n, 1]
σ = d[1:n, 2]

function prediction(ε, σ0)
    ε, σ0 = convert_to_tensor([ε, σ0], [Float64, Float64])
    N = length(ε)
    condition = (i, σ)-> (i<=N)
    function body(i, taσ)
        σl = read(taσ, i-1)
        εc = ε[i]
        εl = ε[i-1]
        σc = H * εc + squeeze(ae(reshape([εl;σl], (1,2)), [20,20,20,1]))
        taσ = write(taσ, i, σc)
        i+1, taσ
    end
    i = constant(2, dtype=Int32)
    taσ = TensorArray(N)
    taσ = write(taσ, 1, σ0)
    _, out = while_loop(condition, body, [i, taσ])
    set_shape(stack(out), (N,))
end

function visualize(iter)
    full = run(sess, σpred_full)

    close("all")
    plot(d[:,1], d[:,2])
    plot(d[1:end,1], full, "--")
    xlabel("Strain")
    ylabel("Stress (MPa)")
    savefig("figures/strainstress$kid$(tid)_$iter.png")

    close("all")
    ref = d[1:end,2]
    plot(1:n-1,ref[1:n-1], "C1", label="Reference")
    plot(1:n-1,full[1:n-1], "C2--", label="Training Prediction")
    plot(n:size(d,1),ref[n:end], "C1")
    plot(n:size(d,1),full[n:end], "C3--", label="Testing Prediction")
    xlabel("Index")
    ylabel("Stress (MPa)")
    legend()

    savefig("figures/stress$kid$(tid)_$iter.png")
    writedlm("figures/res$kid$(tid)_$iter.txt", full)
end

σpred_train = prediction(d[1:n], σ[1])
loss = sum((σpred_train - σ)^2)

σpred_full = prediction(d[1:end,1], σ[1])

sess = Session(); init(sess)

for i = 1:20
    visualize(i-1)
    BFGS!(sess, loss, 100)
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using ManualMemory: MemoryBuffer, load, store!, LazyPreserve, preserve
using Test

@testset "ManualMemory.jl" begin

  @test_throws AssertionError MemoryBuffer{4,String}(undef)
  m = MemoryBuffer{4,Float64}(undef);
  store!(pointer(m), 1.23)
  @test load(pointer(m)) == 1.23
  str = "Hello world!"
  GC.@preserve str begin
    store!(Base.unsafe_convert(Ptr{String}, m), str)
    @test load(Base.unsafe_convert(Ptr{String}, m)) === str
  end

  x = [0 0; 0 0];
  preserve(store!, LazyPreserve(x), 1)
  @test x[1] === 1
end

using ThreadingUtilities
include(joinpath(pkgdir(ThreadingUtilities), "test", "runtests.jl"))


text:
#weave(joinpath(dirname(pathof(MetidaNCA)), "..", "docs", "validation_report.jmd");
weave(joinpath(dirname(@__FILE__), "validation_report.jmd");
#doctype = "pandoc2pdf",
doctype = "pandoc2pdf",
#out_path = joinpath(dirname(pathof(MetidaNCA)), "..", "docs", "src"),
out_path = joinpath(dirname(@__FILE__), "src"),
pandoc_options=["--toc", "-V colorlinks=true" , "-V linkcolor=blue", "-V urlcolor=red",
"-V toccolor=gray", "--number-sections"])

rm(joinpath(dirname(@__FILE__), "src", "validation_report.aux"); force=true)
rm(joinpath(dirname(@__FILE__), "src", "validation_report.log"); force=true)
rm(joinpath(dirname(@__FILE__), "src", "validation_report.out"); force=true)

#=
using Documenter, MetidaNCA, Weave, PrettyTables, CSV, DataFrames
weave(joinpath(dirname(@__FILE__), "validation_report.jmd");
#doctype = "pandoc2pdf",
doctype = "pandoc2pdf",
out_path = ppath,
pandoc_options=["--toc", "-V colorlinks=true" , "-V linkcolor=blue", "-V urlcolor=red", "-V toccolor=gray"])

mainfont: romanuni.ttf
sansfont: NotoSans-Regular.ttf
monofont: NotoSansMono-Regular.ttf
mathfont: texgyredejavu-math.otf

- name: Install TeXlive and Pandoc
        run: sudo apt-get install pandoc texlive texlive-publishers texlive-science latexmk texlive-xetex texlive-latex-base texlive-fonts-recommended
=#


text:
"""
    ILP(optimizer, max_iter)

ILP iteratively solves a linearized primal optimization to compute maximum allowable disturbance.
It iteratively adds the linear constraint to the problem.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: PolytopeComplement

# Return
`AdversarialResult`

# Method
Iteratively solve a linear encoding of the problem.
It only considers the linear piece of the network that has the same activation pattern as the reference input.
Default `optimizer` is `GLPKSolverMIP()`.
We provide both iterative method and non-iterative method to solve the LP problem.
Default `iterative` is `true`.

# Property
Sound but not complete.

# Reference
[O. Bastani, Y. Ioannou, L. Lampropoulos, D. Vytiniotis, A. Nori, and A. Criminisi,
"Measuring Neural Net Robustness with Constraints,"
in *Advances in Neural Information Processing Systems*, 2016.](https://arxiv.org/abs/1605.07262)
"""
@with_kw struct ILP <: Solver
    optimizer = GLPK.Optimizer
    iterative::Bool = true
end

function solve(solver::ILP, problem::Problem)
    nnet = problem.network
    x = problem.input.center
    model = Model(solver)
    model[:δ] = δ = get_activation(nnet, x)
    z = init_vars(model, nnet, :z, with_input=true)
    add_complementary_set_constraint!(model, problem.output, last(z))
    o = max_disturbance!(model, first(z) - problem.input.center)

    if !solver.iterative
        encode_network!(model, nnet, StandardLP())
        optimize!(model)
        termination_status(model) != OPTIMAL && return AdversarialResult(:unknown)
        x = value(first(z))
        return interpret_result(solver, x, problem.input)
    end

    encode_network!(model, nnet, LinearRelaxedLP())
    while true
        optimize!(model)
        termination_status(model) != OPTIMAL && return AdversarialResult(:unknown)
        x = value(first(z))
        matched, index = match_activation(nnet, x, δ)
        if matched
            return interpret_result(solver, x, problem.input)
        end
        add_constraint!(model, nnet, z, δ, index)
    end
end

function interpret_result(solver::ILP, x, input)
    radius = abs.(x .- center(input))

    if all(radius .>= radius_hyperrectangle(input))
        return AdversarialResult(:holds, minimum(radius))
    else
        return AdversarialResult(:violated, minimum(radius))
    end
end

function add_constraint!(model::Model,
                         nnet::Network,
                         z::Vector{Vector{VariableRef}},
                         δ::Vector{Vector{Bool}},
                         (i, j)::Tuple{Int64, Int64})
    layer = nnet.layers[i]
    val = layer.weights[j, :]' * z[i] + layer.bias[j]
    if δ[i][j]
        @constraint(model, val >= 0.0)
    else
        @constraint(model, val <= 0.0)
    end
end

function match_activation(nnet::Network, x::Vector{Float64}, δ::Vector{Vector{Bool}})
    curr_value = x
    for (i, layer) in enumerate(nnet.layers)
        curr_value = layer.weights * curr_value + layer.bias
        matched, j = match_activation(layer, δ[i], curr_value)

        !matched && return false, (i,j)

        curr_value = layer.activation(curr_value)
    end
    return (true, nothing)
end

match_activation(L::Layer{Id}, args...) = (true, nothing)
function match_activation(L::Layer{ReLU}, δᵢ, val, tol = 1e-4)
    for (j, val) in enumerate(val)
        act = δᵢ[j]
        if act
            val < -tol && return false, j
        else
            val > tol  && return false, j
        end

    end
    return (true, nothing)
end


text:
LinearAlgebra.lq(a::IndicesArray, args...; kws...) = lq!(copy(a), args...; kws...)
function LinearAlgebra.lq!(a::IndicesArray, args...; kwargs...)
    inner = lq!(parent(nda), args...; kwargs...)
    return LQ(IndicesArray(getfield(inner, :factors), axes(a)), getfield(inner, :τ))
end
function Base.parent(fact::LQ{T,<:IndicesArray}) where {T}
    return LQ(parent(getfield(fact, :factors)), getfield(fact, :τ))
end

function Base.getproperty(fact::LQ{T,<:IndicesArray}, d::Symbol) where {T}
    if d === :L
        return IndicesArray{(n1, :_)}(getproperty(parent(fact), d), axes(fact, 1), nothing)
    elseif d === :Q
        return IndicesArray{(:_, n2)}(getproperty(parent(fact), d), nothing, axes(fact, 2))
    else
        return getproperty(parent(fact), d)
    end
end


text:
using BenchmarkTools
using Clustering, Distances, Random

Random.seed!(345678)

const SUITE = BenchmarkGroup()

SUITE["hclust"] = BenchmarkGroup()

function random_distance_matrix(n::Integer, m::Integer=10, dist::PreMetric=Euclidean())
    pts = rand(m, n)
    return Distances.pairwise(dist, pts, dims=2)
end

function hclust_benchmark(n::Integer, m::Integer=10, dist::PreMetric=Euclidean())
    res = BenchmarkGroup()
    for linkage in ("single", "average", "complete")
        res[linkage] = @benchmarkable hclust(D, linkage=Symbol($linkage)) setup=(D=random_distance_matrix($n, $m, $dist))
    end
    return res
end

for n in (10, 100, 1000, 10000)
    SUITE["hclust"][n] = hclust_benchmark(n)
end

SUITE["cutree"] = BenchmarkGroup()

for (n, k) in ((10, 3), (100, 10), (1000, 100), (10000, 1000))
    SUITE["cutree"][(n,k)] = @benchmarkable cutree(hclu, k=$k) setup=(D=random_distance_matrix($n, 5); hclu=hclust(D, linkage=:single))
end


text:
#module Bruker

export readBrukerFID

function readBrukerFID(f::IOStream; format=Float64)
  data = Array{format,1}([])
  while !eof(f)
    append!(data,(read(f,format)));
  end
  return data[1:2:end]+im*data[2:2:end]
end

"""
`Bruker.readBrukerFID(s::String; format=Float64)` reads a Bruker fid or ser file and returns
the data contained in it as a complex array. For files generated by TopSpin 4.0 and
higher, the default format (`Float64`) will work. For TopSpin 2.0, use
`format=Int32`.
"""
function readBrukerFID(s::String;opts...)
    f=open(s,"r")
    fid=readBrukerFID(f;opts...)
    close(f)
    return fid
end

#end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR