text:
export audit, collect_files

using ObjectFile

# AUDITOR TODO LIST:
#
# * Auto-determine minimum glibc version (to sate our own curiosity)
# * Detect instruction sets that are non-portable

"""
    audit(prefix::Prefix; platform::Symbol = platform_key();
                          verbose::Bool = false,
                          autofix::Bool = false)

Audits a prefix to attempt to find deployability issues with the binary objects
that have been installed within.  This auditing will check for relocatability
issues such as dependencies on libraries outside of the current `prefix`,
usage of advanced instruction sets such as AVX2 that may not be usable on many
platforms, linkage against newer glibc symbols, etc...

This method is still a work in progress, only some of the above list is
actually implemented, be sure to actually inspect `Auditor.jl` to see what is
and is not currently in the realm of fantasy.
"""
function audit(prefix::Prefix; platform::Symbol = platform_key(),
                               verbose::Bool = false,
                               autofix::Bool = false)
    if verbose
        info("Beginning audit of $(prefix.path)")
    end

    # If this is false then it's bedtime for bonzo boy
    all_ok = true

    # Inspect binary files, looking for improper linkage
    predicate = f -> (filemode(f) & 0o111) != 0 || valid_dl_path(f, platform)
    bin_files = collect_files(prefix, predicate)
    for f in bin_files
        # Peel this binary file open like a delicious tangerine
        oh = try
            readmeta(f)
        catch
            # If this isn't an actual binary file, skip it
            if verbose
                info("Skipping binary analysis of $(relpath(f, prefix.path))")
            end
            continue
        end
        rp = RPath(oh)

        if verbose
            msg = strip("""
            Checking $(relpath(f, prefix.path)) with RPath list $(rpaths(rp))
            """)
            info(msg)
        end

        # Look at every non-default dynamic link
        libs = filter_default_linkages(find_libraries(oh), oh)
        for libname in keys(libs)
            if !isfile(libs[libname])
                # If we couldn't resolve this library, let's try autofixing,
                # if we're allowed to by the user
                if autofix
                    # First, is this a library that we already know about?
                    known_bins = lowercase.(basename.(bin_files))
                    kidx = findfirst(known_bins .== lowercase(basename(libname)))
                    if kidx > 0
                        # If it is, point to that file instead!
                        new_link = update_linkage(prefix, platform, path(oh), libs[libname], bin_files[kidx]; verbose=verbose)

                        if verbose
                            msg = replace("""
                            Linked library $(libname) has been auto-mapped to
                            $(new_link)
                            """, '\n', ' ')
                            info(strip(msg))
                        end
                    else
                        msg = replace("""
                        Linked library $(libname) could not be resolved and
                        could not be auto-mapped
                        """, '\n', ' ')
                        warn(strip(msg))
                        all_ok = false
                    end
                else
                    msg = replace("""
                    Linked library $(libname) could not be resolved within
                    the given prefix
                    """, '\n', ' ')
                    warn(strip(msg))
                    all_ok = false
                end
            elseif !startswith(libs[libname], prefix.path)
                msg = replace("""
                Linked library $(libname) (resolved path $(libs[libname]))
                is not within the given prefix
                """, '\n', ' ')
                warn(strip(msg))
                all_ok = false
            end
        end
    end

    # Inspect all shared library files for our platform (but only if we're
    # running native, don't try to load library files from other platforms)
    if platform == platform_key()
        # Find all dynamic libraries
        predicate = f -> valid_dl_path(f, platform)
        shlib_files = collect_files(prefix, predicate)

        for f in shlib_files
            if verbose
                info("Checking shared library $(relpath(f, prefix.path))")
            end
            hdl = Libdl.dlopen_e(f)
            if hdl == C_NULL
                # TODO: Use the relevant ObjFileBase packages to inspect why this
                # file is being nasty to us.

                warn("$(relpath(f, prefix.path)) cannot be dlopen()'ed")
                all_ok = false
            else
                Libdl.dlclose(hdl)
            end
        end
    end

    # Search _every_ file in the prefix path to find hardcoded paths
    predicate = f -> !startswith(f, joinpath(prefix, "logs")) &&
                     !startswith(f, joinpath(prefix, "manifests"))
    all_files = collect_files(prefix, predicate)

    # Finally, check for absolute paths in any files.  This is not a "fatal"
    # offense, as many files have absolute paths.  We want to know about it
    # though, so we'll still warn the user.
    for f in all_files
        file_contents = readstring(f)
        if contains(file_contents, prefix.path)
            warn("$(relpath(f, prefix.path)) contains an absolute path")
        end
    end

    return all_ok
end

"""
    collect_files(prefix::Prefix, predicate::Function = f -> true)

Find all files that satisfy `predicate()` when the full path to that file is
passed in, returning the list of file paths.
"""
function collect_files(prefix::Prefix, predicate::Function = f -> true)
    collected = String[]
    real_paths = String[]
    for (root, dirs, files) in walkdir(prefix.path)
        for f in files
            f_path = joinpath(root, f)

            # Calculate the realpath, but keep the nicely formatted path too
            f_real_path = realpath(f_path)

            # Only add this file into our list if it is not already contained.
            # This removes duplicate symlinks
            if !(f_real_path in real_paths) && predicate(f_path)
                push!(collected, f_path)
                push!(real_paths, realpath(f_path))
            end
        end
    end
    return collected
end

"""
    filter_default_linkages(libs::Dict, oh::ObjectHandle)

Given libraries obtained through `ObjectFile.find_libraries()`, filter out
libraries that are "default" libraries and should be available on any system.
"""
function filter_default_linkages(libs::Dict, oh::ObjectHandle)
    return Dict(k => libs[k] for k in keys(libs) if !should_ignore_lib(k, oh))
end

function should_ignore_lib(lib, ::ELFHandle)
    default_libs = [
        "libc.so.6",
        "libgcc_s.1.so",
    ]
    return lowercase(basename(lib)) in default_libs
end

function should_ignore_lib(lib, ::MachOHandle)
    default_libs = [
        "libsystem.b.dylib",
        "libgcc_s.1.dylib",
    ]
    return lowercase(basename(lib)) in default_libs
end

function should_ignore_lib(lib, ::COFFHandle)
    default_libs = [
        "msvcrt.dll",
        "kernel32.dll",
        "user32.dll",
        "libgcc_s_sjlj-1.dll",
    ]
    return lowercase(basename(lib)) in default_libs
end

"""
    update_linkage(prefix::Prefix, platform::Symbol, path::AbstractString,
                   old_libpath, new_libpath; verbose::Bool = false)

Given a binary object located at `path` within `prefix`, update its dynamic
linkage to point to `new_libpath` instead of `old_libpath`.  This is done using
a tool within the cross-compilation environment such as `install_name_tool` on
MacOS or `patchelf` on Linux.  Windows platforms are completely skipped, as
they do not encode paths or RPaths within their executables.
"""
function update_linkage(prefix::Prefix, platform::Symbol, path::AbstractString,
                        old_libpath, new_libpath; verbose::Bool = false)
    # Windows doesn't do updating of linkage
    if is_windows(platform)
        return
    end

    dr = DockerRunner(;prefix=prefix, platform=platform)

    add_rpath = x -> ``
    relink = (x, y) -> ``
    origin = ""
    if is_apple(platform)
        origin = "@loader_path"
        add_rpath = rp -> `install_name_tool -add_rpath $(rp) $(path)`
        relink = (op, np) -> `install_name_tool -change $(op) $(np) $(path)`
    elseif is_linux(platform)
        origin = "\$ORIGIN"
        full_rpath = join(':', rpaths(RPath(readmeta(path))))
        add_rpath = rp -> `patchelf --set-rpath $(full_rpath):$(rp) $(path)`
        relink = (op, np) -> `patchelf --replace-needed $(op) $(np) $(path)`
    end

    if !(dirname(new_libpath) in canonical_rpaths(RPath(readmeta(path))))
        libname = basename(old_libpath)
        logpath = joinpath(logdir(prefix), "update_rpath_$(libname).log")
        cmd = add_rpath(relpath(dirname(new_libpath), dirname(path)))
        run(dr, cmd, logpath; verbose=verbose)
    end

    # Create a new linkage that looks like $ORIGIN/../lib, or similar
    libname = basename(old_libpath)
    logpath = joinpath(logdir(prefix), "update_linkage_$(libname).log")
    origin_relpath = joinpath(origin, relpath(new_libpath, dirname(path)))
    cmd = relink(old_libpath, origin_relpath)
    run(dr, cmd, logpath; verbose=verbose)

    return origin_relpath
end


text:
"""
    Thermometer{T <: Real} <: AbstractEncoder

Thermometer encoder, used to represent real values as binary patterns.
"""
struct Thermometer{T <: Real} <: AbstractEncoder
    min::Union{T,AbstractVector{T}}
    max::Union{T,AbstractVector{T}}
    resolution::Int # TODO: This could be unsigned

    """
        Thermometer{T}(min::T, max::T, resolution::Int) where {T <: Real}

    Constructs a thermometer encoder. `min` and `max` denotes the admissible
    numerical limits of encoded values and `resolution` specifies the number of
    bits dedicated to the encoding of a scalar. Values to be encoded must be of
    type T.
    """
    function Thermometer{T}(min::T, max::T, resolution::Int) where {T <: Real}
        max ≤ min && throw(ArgumentError(
            "the maximum of the thermometer may not be lesser or equal to its minimum"
        ))

        resolution < 1 && throw(DomainError(
            resolution,
            "the thermometer's resolution must be greater then zero"
        ))

        new(min, max, resolution)
    end

    """
        Thermometer{T}(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real}

    Constructs a thermometer encoder. `min` and `max` denotes the admissible
    numerical limits of encoded values and `resolution` specifies the number of
    bits dedicated to the encoding of a scalar. Values to be encoded are
    expected to have the same shape as `min` and `max`. Values to be encoded
    must be of type T.
    """
    function Thermometer{T}(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real}
        length(min) != length(max) && throw(DimensionMismatch(
            "min and max vectors must have the same length"
        ))

        !all(min .≤ max) && throw(ArgumentError(
            "the maximum of the thermometer may not be lesser or equal to its minimum"
        ))

        resolution < 1 && throw(DomainError(
            resolution,
            "the thermometer's resolution must be greater then zero"
        ))

        new(min, max, resolution)
    end
end

"""
    Thermometer(min::T, max::T, resolution::Int) where {T <: Real}

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. Values to be encoded must be of type
T.
"""
Thermometer(min::T, max::T, resolution::Int) where {T <: Real} = Thermometer{T}(min, max, resolution)

"""
    Thermometer(min::Real, max::Real, resolution::Int)

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. Values to be encoded must be same as
the promotion of `min` and `max` types.
"""
Thermometer(min::Real, max::Real, resolution::Int) = Thermometer(promote(min, max)..., resolution)

"""
    Thermometer(::Type{T}, min::Real, max::Real, resolution) where {T <: Real}

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of
bits dedicated to the encoding of a scalar. First argument determines the type
of values to be encoded.
"""
Thermometer(::Type{T}, min::Real, max::Real, resolution) where {T <: Real} = Thermometer(convert.(T, (min, max))..., resolution)

"""
    Thermometer(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real}

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of bits
dedicated to the encoding of a scalar. Values to be encoded are expected to have
the same shape as `min` and `max`. Values to be encoded must be of type T.
"""
Thermometer(min::AbstractVector{T}, max::AbstractVector{T}, resolution::Int) where {T <: Real} = Thermometer{T}(min, max, resolution)

"""
    Thermometer(min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution::Int)

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of bits
dedicated to the encoding of a scalar. Values to be encoded are expected to have
the same shape as `min` and `max`. Values to be encoded must be same as the
promotion of `min` and `max` types.
"""
Thermometer(min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution::Int) = Thermometer(promote(min, max)..., resolution)

"""
    Thermometer(::Type{T}, min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution)

Constructs a thermometer encoder. `min` and `max` denotes the admissible
numerical limits of encoded values and `resolution` specifies the number of bits
dedicated to the encoding of a scalar. Values to be encoded are expected to have
the same shape as `min` and `max`. First argument determines the type of values
to be encoded.
"""
Thermometer(::Type{T}, min::AbstractVector{<:Real}, max::AbstractVector{<:Real}, resolution) where {T <: Real} = Thermometer(convert.(Vector{T}, (min, max))..., resolution)

function encode!(encoder::Thermometer{T}, min::T, max::T, x::T, pattern::AbstractVector{Bool}) where {T <: Real}
    for i in eachindex(pattern)
        pattern[i] = x > min + (i - 1) * (max - min) / encoder.resolution
    end

    pattern
end

function encode!(encoder::Thermometer{T}, x::T, pattern::AbstractVector{Bool}) where {T <: Real}
    validate_input_output(encoder, x, pattern)

    encode!(encoder, encoder.min[1], encoder.max[1], x, pattern)
end

ext: jl

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
export
	GoalConstraint,
	BoundConstraint,
	CircleConstraint,
	SphereConstraint,
	NormConstraint,
	LinearConstraint,
	VariableBoundConstraint,
	QuatNormConstraint,
	QuatSlackConstraint

############################################################################################
#                              GOAL CONSTRAINTS 										   #
############################################################################################

"""
$(TYPEDEF)
Constraint of the form ``x_g = a``, where ``x_g`` can be only part of the state
vector.

# Constructors:
```julia
GoalConstraint(xf::AbstractVector)
GoalConstraint(xf::AbstractVector, inds)
```
where `xf` is an n-dimensional goal state. If `inds` is provided,
only `xf[inds]` will be used.
"""
struct GoalConstraint{T,P,N,L} <: AbstractConstraint{Equality,State,P}
	xf::SVector{P,T}
	Ix::SMatrix{P,N,T,L}
	inds::SVector{P,Int}
end

function GoalConstraint(xf::AbstractVector, inds=1:length(xf))
	n = length(xf)
	p = length(inds)
	xf = SVector{n}(xf)
	inds = SVector{p}(inds)
	GoalConstraint(xf, inds)
end

function GoalConstraint(xf::SVector{n}, inds::SVector{p,Int}) where {n,p}
	Ix = SMatrix{n,n}(Matrix(1.0I,n,n))
	Ix = Ix[inds,:]
	GoalConstraint(SVector{p}(xf[inds]), Ix, inds)
end

state_dim(::GoalConstraint{T,P,N}) where {T,P,N} = N
evaluate(con::GoalConstraint, x::SVector) = x[con.inds] - con.xf
jacobian(con::GoalConstraint, z::KnotPoint) = con.Ix

############################################################################################
#                              LINEAR CONSTRAINTS 										   #
############################################################################################
"""
$(TYPEDEF)
Linear constraint of the form ``Ay - b \\{\\leq,=\\} 0`` where ``y`` may be either the
state or controls (but not a combination of both).

# Constructor: ```julia
LinearConstraint{S,W}(n,m,A,b)
```
where `W <: Union{State,Control}`.
"""
struct LinearConstraint{S,W<:Union{State,Control},P,N,L,T} <: AbstractConstraint{S,W,P}
	n::Int
	m::Int
	A::SMatrix{P,N,T,L}
	b::SVector{P,T}
end

function LinearConstraint{S,W}(n,m, A::AbstractMatrix, b::AbstractVector) where {S,W}
	@assert size(A,1) == length(b)
	p,q = size(A)
	A = SMatrix{p,q}(A)
	b = SVector{p}(b)
	LinearConstraint{S,W}(n,m, A, b)
end

function LinearConstraint{S,W}(n::Int,m::Int, A::SMatrix{P,N,T}, b::SVector{P,T}) where {S,W,P,N,T}
	con = LinearConstraint{S,W,P,N,P*N,T}(n,m,A,b)
	@assert check_dims(con,n,m) "Dimensions of LinearConstraint are inconsistent"
	return con
end

state_dim(::LinearConstraint{S,State,P,N}) where {S,P,N} = N
control_dim(::LinearConstraint{S,Control,P,N}) where {S,P,N} = N
evaluate(con::LinearConstraint,x::SVector) = con.A*x - con.b
jacobian(con::LinearConstraint,x::SVector) = con.A

############################################################################################
#                              CIRCLE/SPHERE CONSTRAINTS 								   #
############################################################################################
"""
$(TYPEDEF)
Constraint of the form
`` (x - x_c)^2 + (y - y_c)^2 \\leq r^2 ``
where ``x``, ``y`` are given by `x[xi]`,`x[yi]`, ``(x_c,y_c)`` is the center
of the circle, and ``r`` is the radius.

# Constructor:
```julia
CircleConstraint(n, xc::SVector{P}, yc::SVector{P}, radius::SVector{P}, xi=1, yi=2)
```
"""
struct CircleConstraint{T,P} <: AbstractConstraint{Inequality,State,P}
	n::Int
	x::SVector{P,T}
	y::SVector{P,T}
	radius::SVector{P,T}
	xi::Int  # index of x-state
	yi::Int  # index of y-state
	CircleConstraint(n::Int, xc::SVector{P,T}, yc::SVector{P,T}, radius::SVector{P,T},
			xi=1, yi=2) where {T,P} =
		 new{T,P}(n,xc,yc,radius,xi,yi)
end
state_dim(con::CircleConstraint) = con.n

function evaluate(con::CircleConstraint{T,P}, x::SVector) where {T,P}
	xc = con.x; xi = con.xi
	yc = con.y; yi = con.yi
	r = con.radius
	-(x[xi] .- xc).^2 - (x[yi] .- yc).^2 + r.^2
end

"""
$(TYPEDEF)
Constraint of the form
`` (x - x_c)^2 + (y - y_c)^2 + (z - z_c)^2 \\leq r^2 ``
where ``x``, ``y``, ``z`` are given by `x[xi]`,`x[yi]`,`x[zi]`, ``(x_c,y_c,z_c)`` is the center
of the sphere, and ``r`` is the radius.

# Constructor:
```
SphereConstraint(n, xc::SVector{P}, yc::SVector{P}, zc::SVector{P},
	radius::SVector{P}, xi=1, yi=2, zi=3)
```
"""
struct SphereConstraint{T,P} <: AbstractConstraint{Inequality,State,P}
	n::Int
	x::SVector{P,T}
	y::SVector{P,T}
	z::SVector{P,T}
	xi::Int
	yi::Int
	zi::Int
	radius::SVector{P,T}
	SphereConstraint(n::Int, xc::SVector{P,T}, yc::SVector{P,T}, zc::SVector{P,T},
			radius::SVector{P,T}, xi=1, yi=2, zi=3) where {T,P} =
			new{T,P}(n,xc,yc,zc,xi,yi,zi,radius)
end

state_dim(con::SphereConstraint) = con.n

function evaluate(con::SphereConstraint{T,P}, x::SVector) where {T,P}
	xc = con.x; xi = con.xi
	yc = con.y; yi = con.yi
	zc = con.z; zi = con.zi
	r = con.radius

	-((x[xi] .- xc).^2 + (x[yi] .- yc).^2 + (x[zi] .- zc).^2 - r.^2)
end

############################################################################################
#  								SELF-COLLISION CONSTRAINT 								   #
############################################################################################

struct CollisionConstraint{D} <: AbstractConstraint{Inequality,State,1}
	n::Int
    x1::SVector{D,Int}
    x2::SVector{D,Int}
    radius::Float64
end

state_dim(con::CollisionConstraint) = con.n

function evaluate(con::CollisionConstraint, x::SVector)
    x1 = x[con.x1]
    x2 = x[con.x2]
    d = x1 - x2
    @SVector [con.radius^2 - d'd]
end

############################################################################################
#								NORM CONSTRAINT											   #
############################################################################################

"""
$(TYPEDEF)
Constraint of the form
``\\|y\\|^2 \\{\\leq,=\\} a``
where ``y`` is either a state or a control vector (but not both)

# Constructors:
```
NormConstraint{S,State}(n,a)
NormConstraint{S,Control}(m,a)
```
where `a` is the constant on the right-hand side of the equation.

# Examples:
```julia
NormConstraint{Equality,Control}(2,4.0)
```
creates a constraint equivalent to
``\\|u\\|^2 = 4.0`` for a problem with 2 controls.

```julia
NormConstraint{Inequality,State}(3, 2.3)
```
creates a constraint equivalent to
``\\|x\\|^2 \\leq 2.3`` for a problem with 3 states.
"""
struct NormConstraint{S,W<:Union{State,Control},T} <: AbstractConstraint{S,W,1}
	dim::Int
	val::T
	function NormConstraint{S,W,T}(dim::Int, val::T) where {S,W<:Union{State,Control},T}
		@assert val ≥ 0 "Value must be greater than or equal to zero"
		new{S,W,T}(dim, val)
	end
end
NormConstraint{S,W}(n::Int, val::T) where {S,W,T} = NormConstraint{S,W,T}(n, val)

state_dim(con::NormConstraint{S,State}) where S = con.dim
control_dim(con::NormConstraint{S,Control}) where S = con.dim

function evaluate(con::NormConstraint, x::SVector)
	return @SVector [norm(x)^2 - con.val]
end

struct QuatNormConstraint <: AbstractConstraint{Equality,State,1}
	n::Int
	qinds::SVector{4,Int}
end

QuatNormConstraint(n::Int=13, qinds=(@SVector [4,5,6,7])) = QuatNormConstraint(n, qinds)

state_dim(con::QuatNormConstraint) = con.n

function evaluate(con::QuatNormConstraint, x::SVector)
	q = x[con.qinds]
	return @SVector [norm(q) - 1.0]
end

struct QuatSlackConstraint <: AbstractConstraint{Equality,Stage,1}
	qinds::SVector{4,Int}
end
(::Type{<:QuatSlackConstraint})(qinds=(@SVector [4,5,6,7])) = QuatSlackConstraint(qinds)

state_dim(::QuatSlackConstraint) = 13
control_dim(::QuatSlackConstraint) = 5  # special cased for quadrotor

function evaluate(con::QuatSlackConstraint, x::SVector, u::SVector)
	s = u[end]
	q = x[con.qinds]
	return @SVector [norm(q)*s - 1.0]
end

function jacobian(con::QuatSlackConstraint, x::SVector, u::SVector)
	s = u[end]
	q = x[con.qinds]
	nq = norm(q)
	M = s/nq
	return @SMatrix [0 0 0 q[1]*M  q[2]*M  q[3]*M  q[4]*M  0 0 0  0 0 0  0 0 0 0 nq]
end

############################################################################################
# 								COPY CONSTRAINT 										   #
############################################################################################

# struct CopyConstraint{K,W,S,P,N,M} <: AbstractConstraint{W,S,P}
# 	con::AbstractConstraint{W,S,P}
#     xinds::Vector{SVector{N,Int}}
#     uinds::Vector{SVector{M,Int}}
# end
#
# function evaluate(con::CopyConstraint{K}, z::KnotPoint)
# 	c = evaluate(con,)
# 	for 2 = 1:K
# 	end
# end

############################################################################################
# 								BOUND CONSTRAINTS 										   #
############################################################################################
"""$(TYPEDEF) Linear bound constraint on states and controls
# Constructors
```julia
BoundConstraint(n, m; x_min, x_max, u_min, u_max)
```
Any of the bounds can be ±∞. The bound can also be specifed as a single scalar, which applies the bound to all state/controls.
"""
struct BoundConstraint{T,P,NM,PNM} <: AbstractConstraint{Inequality,Stage,P}
	n::Int
	m::Int
	z_max::SVector{NM,T}
	z_min::SVector{NM,T}
	b::SVector{P,T}
	B::SMatrix{P,NM,T,PNM}
	inds::SVector{P,Int}
end

function BoundConstraint(n, m; x_max=Inf*(@SVector ones(n)), x_min=-Inf*(@SVector ones(n)),
		u_max=Inf*(@SVector ones(m)), u_min=-Inf*(@SVector ones(m)))
	# Check and convert bounds
	x_max, x_min = checkBounds(Val(n), x_max, x_min)
	u_max, u_min = checkBounds(Val(m), u_max, u_min)

	# Concatenate bounds
	z_max = [x_max; u_max]
	z_min = [x_min; u_min]
	b = [-z_max; z_min]
	bN = [x_max; u_max*Inf; x_min; -u_min*Inf]

	active = isfinite.(b)
	p = sum(active)
	inds = SVector{p}(findall(active))

	B = SMatrix{2(n+m), n+m}([1.0I(n+m); -1.0I(n+m)])

	BoundConstraint(n, m, z_max, z_min, b[inds], B[inds,:], inds)
end

function con_label(con::BoundConstraint, ind::Int)
	i = con.inds[ind]
	n,m = state_dim(con), control_dim(con)
	if 1 <= i <= n
		return "x max $i"
	elseif n < i <= n + m
		j = i - n
		return "u max $j"
	elseif n + m < i <= 2n+m
		j = i - (n+m)
		return "x min $j"
	elseif 2n+m < i <= 2n+2m
		j = i - (2n+m)
		return "u min $j"
	else
		throw(BoundsError())
	end
end

function checkBounds(::Val{N}, u::AbstractVector, l::AbstractVector) where N
	if all(u .>= l)
		return SVector{N}(u), SVector{N}(l)
	else
		throw(ArgumentError("Upper bounds must be greater than or equal to lower bounds"))
	end
end

checkBounds(sze::Val{N}, u::Real, l::Real) where N =
	checkBounds(sze, (@SVector fill(u,N)), (@SVector fill(l,N)))
checkBounds(sze::Val{N}, u::AbstractVector, l::Real) where N =
	checkBounds(sze, u, (@SVector fill(l,N)))
checkBounds(sze::Val{N}, u::Real, l::AbstractVector) where N =
	checkBounds(sze, (@SVector fill(u,N)), l)

state_dim(con::BoundConstraint) = con.n
control_dim(con::BoundConstraint) = con.m
is_bound(::BoundConstraint) = true
lower_bound(bnd::BoundConstraint) = bnd.z_min
upper_bound(bnd::BoundConstraint) = bnd.z_max

function evaluate(bnd::BoundConstraint{T,P,NM}, x, u) where {T,P,NM}
	bnd.B*SVector{NM}([x; u]) + bnd.b
end

function jacobian!(∇c, bnd::BoundConstraint, z::AbstractKnotPoint)
	∇c .= bnd.B
end

############################################################################################
#  							VARIABLE BOUND CONSTRAINT 									   #
############################################################################################

struct VariableBoundConstraint{T,P,NM,PNM} <: AbstractConstraint{Inequality,Stage,P}
	n::Int
	m::Int
	z_max::Vector{SVector{NM,T}}
	z_min::Vector{SVector{NM,T}}
	b::Vector{SVector{P,T}}
	B::SMatrix{P,NM,T,PNM}
	function VariableBoundConstraint(n::Int,m::Int,
			z_max::Vector{<:SVector{NM,T}}, z_min::Vector{<:SVector{NM,T}},
			b::Vector{<:SVector{P}}, B::SMatrix{P,NM,T,PNM}) where {T,P,PN,NM,PNM}
		new{T,P,NM,PNM}(n,m,z_max,z_min,b,B)
	end
end

state_dim(con::VariableBoundConstraint) = con.n
control_dim(con::VariableBoundConstraint) = con.m
is_bound(::VariableBoundConstraint) = true

function evaluate!(vals::Vector{<:AbstractVector},
		con::VariableBoundConstraint, Z::Traj, inds=1:length(Z)-1)
	for (i,k) in enumerate(inds)
		vals[i] = con.B*Z[k].z + con.b[k]
	end
end

function jacobian(con::VariableBoundConstraint, z::KnotPoint)
	return con.B
end

function VariableBoundConstraint(n, m, N;
		x_max=[Inf*(@SVector ones(n)) for k = 1:N], x_min=[-Inf*(@SVector ones(n)) for k = 1:N],
		u_max=[Inf*(@SVector ones(m)) for k = 1:N], u_min=[-Inf*(@SVector ones(m)) for k = 1:N])
	@assert length(x_max) == N
	@assert length(u_max) == N
	@assert length(x_min) == N
	@assert length(u_min) == N

	# Check and convert bounds
	for k = 1:N
		x_max[k], x_min[k] = checkBounds(Val(n), x_max[k], x_min[k])
		u_max[k], u_min[k] = checkBounds(Val(m), u_max[k], u_min[k])
	end

	# Concatenate bounds
	z_max = [SVector{n+m}([x_max[k]; u_max[k]]) for k = 1:N]
	z_min = [SVector{n+m}([x_min[k]; u_min[k]]) for k = 1:N]
	b = [[-z_max[k]; z_min[k]] for k = 1:N]

	active = map(x->isfinite.(x), b)
	equal_active = all(1:N-2) do k
		active[k] == active[k+1]
	end
	if !equal_active
		throw(ArgumentError("All bounds must have the same active constraints"))
	end
	active = active[1]
	p = sum(active)

	inds = SVector{p}(findall(active))

	b = [bi[inds] for bi in b]
	B = SMatrix{2(n+m), n+m}([1.0I(n+m); -1.0I(n+m)])

	VariableBoundConstraint(n, m, z_max, z_min, b, B[inds,:])
end

############################################################################################
#  								INFEASIBLE CONSTRAINT 									   #
############################################################################################
""" $(TYPEDEF) Constraints additional ``infeasible'' controls to be zero.
Constructors: ```julia
InfeasibleConstraint(model::InfeasibleModel)
InfeasibleConstraint(n,m)
```
"""
struct InfeasibleConstraint{N,M} <: AbstractConstraint{Equality, Control, N} end

InfeasibleConstraint(model::InfeasibleModel{N,M}) where {N,M} = InfeasibleConstraint{N,M}()
InfeasibleConstraint(n::Int, m::Int) = InfeasibleConstraint{n,m}()
control_dim(::InfeasibleConstraint{N,M}) where {N,M} = N+M

@generated function evaluate(con::InfeasibleConstraint{N,M}, u::SVector) where {N,M}
    _u = SVector{M}(1:M)
    _ui = SVector{N}((1:N) .+ M)
	quote
        ui = u[$_ui] # infeasible controls
	end
end

@generated function jacobian(con::InfeasibleConstraint{N,M}, u::SVector) where {N,M}
	Iu = [(@SMatrix zeros(N,M)) Diagonal(@SVector ones(N))]
	return :($Iu)
end

############################################################################################
#  								INDEXED CONSTRAINT 	 									   #
############################################################################################
""" $(TYPEDEF) Compute a constraint on an arbitrary portion of either the state or control,
or both. Useful for dynamics augmentation. e.g. you are controlling two models, and have
individual constraints on each. You can define constraints as if they applied to the individual
model, and then wrap it in an `IndexedConstraint` to apply it to the appropriate portion of
the concatenated state. Assumes the indexed state portion is contiguous.

Type params:
* S - Inequality or Equality
* W - ConstraintType
* P - Constraint length
* N,M - original state and control dimensions
* NM - N+M
* Bx - location of the first element in the state index
* Bu - location of the first element in the control index
* C - type of original constraint

Constructors:
```julia
IndexedConstraint(n, m, con)
IndexedConstraint(n, m, con, ix::SVector, iu::SVector)
```
where the arguments `n` and `m` are the state and control dimensions of the new dynamics.
`ix` and `iu` are the indices into the state and control vectors. If left out, they are
assumed to start at the beginning of the vector.

NOTE: Only part of this functionality has been tested. Use with caution!
"""
struct IndexedConstraint{S,W,P,N,M,w,C} <: AbstractConstraint{S,W,P}
	n::Int  # new dimension
	m::Int  # new dimension
	con::C
	ix::SVector{N,Int}
	iu::SVector{M,Int}
	∇c::SizedMatrix{P,w,Float64,2}
	A::SubArray{Float64,2,SizedMatrix{P,w,Float64,2},Tuple{UnitRange{Int},UnitRange{Int}},false}
	B::SubArray{Float64,2,SizedMatrix{P,w,Float64,2},Tuple{UnitRange{Int},UnitRange{Int}},false}
end

state_dim(con::IndexedConstraint{<:Any,<:Union{Stage,State}}) = con.n
control_dim(con::IndexedConstraint{<:Any,<:Union{Stage,Control}}) = con.m
Base.length(::IndexedConstraint{S,W,P}) where {S,W,P} = P

function IndexedConstraint(n,m,con::AbstractConstraint{S,W,P},
		ix::SVector{N}, iu::SVector{M}) where {S,W,P,N,M}
	x = @SVector rand(N)
	u = @SVector rand(M)
	w = width(con)
	∇c = SizedMatrix{P,w}(zeros(P,w))
	if W == Stage
		A = view(∇c, 1:P, 1:N)
		B = view(∇c, 1:P, N .+ (1:M))
	else
		A = view(∇c, 1:0, 1:0)
		B = view(∇c, 1:0, 1:0)
	end
	IndexedConstraint{S,W,P,N,M,w,typeof(con)}(n,m,con,ix,iu,∇c,A,B)
end

function IndexedConstraint(n,m,con::AbstractConstraint{S,W}) where {S,W}
	if W <: Union{State,CoupledState}
		m0 = m
	else
		m0 = control_dim(con)
	end
	if W<: Union{Control,CoupledControl}
		n0 = n
	else
		n0 = state_dim(con)
	end
	ix = SVector{n0}(1:n0)
	iu = SVector{m0}(1:m0)
	IndexedConstraint(n,m,con, ix, iu)
end

# TODO: define higher-level evaluate! function instead
@generated function evaluate(con::IndexedConstraint{<:Any,<:Stage,<:Any,N,M}, z::KnotPoint) where {N,M}
	ix = SVector{N}(1:N)
	iu = N .+ SVector{M}(1:M)
	return quote
		x0 = state(z)[con.ix]
		u0 = control(z)[con.iu]
		z_ = StaticKnotPoint([x0; u0], $ix, $iu, z.dt, z.t)
		evaluate(con.con, z_)
	end
end

# TODO: define higher-leel jacobian! function instead
@generated function jacobian!(∇c, con::IndexedConstraint{<:Any,Stage,P,N0,M0},
		z::KnotPoint{<:Any,N}) where {P,N0,M0,N}
	iP = 1:P
	ix = SVector{N0}(1:N0)
	iu = SVector{M0}(N0 .+ (1:M0))
	if eltype(∇c) <: SizedMatrix
		assignment = quote
			uview(∇c.data,$iP,iA) .= con.A
			uview(∇c.data,$iP,iB) .= con.B
		end
	else
		assignment = quote
			uview(∇c,$iP,iA) .= con.A
			uview(∇c,$iP,iB) .= con.B
		end
	end
	quote
		x0 = state(z)[con.ix]
		u0 = control(z)[con.iu]
		z_ = StaticKnotPoint([x0;u0], $ix, $iu, z.dt, z.t)
		jacobian!(con.∇c, con.con, z_)
		iA = con.ix
		iB = N .+ con.iu
		$assignment
	end
end

@generated function jacobian!(∇c, con::IndexedConstraint{<:Any,State,P,N0,M0},
		z::KnotPoint{<:Any,N}) where {P,N0,M0,N}
	iP = 1:P
	ix = SVector{N0}(1:N0)
	iu = SVector{M0}(N0 .+ (1:M0))
	if eltype(∇c) <: SizedArray
		assignment = :(uview(∇c.data,$iP,iA) .= con.∇c)
	else
		assignment = :(uview(∇c,$iP,iA) .= con.∇c)
	end
	quote
		x0 = state(z)[con.ix]
		u0 = control(z)[con.iu]
		z_ = StaticKnotPoint([x0;u0], $ix, $iu, z.dt, z.t)
		jacobian!(con.∇c, con.con, z_)
		iA = con.ix
		$assignment
	end
end

@generated function jacobian!(∇c, con::IndexedConstraint{<:Any,Control,P,N0,M0},
		z::KnotPoint{<:Any,N}) where {P,N0,M0,N}
	iP = 1:P
	ix = SVector{N0}(1:N0)
	iu = SVector{M0}(N0 .+ (1:M0))
	if eltype(∇c) <: SizedArray
		assignment = :(uview(∇c.data,$iP,iB) .= con.∇c)
	else
		assignment = :(uview(∇c,$iP,iB) .= con.∇c)
	end
	quote
		x0 = state(z)[con.ix]
		u0 = control(z)[con.iu]
		z_ = StaticKnotPoint([x0;u0], $ix, $iu, z.dt, z.t)
		jacobian!(con.∇c, con.con, z_)
		iB = con.iu
		$assignment
	end
end


text:
# by Lazaro Alonso # HIDE
using Pkg # HIDE
Pkg.activate(@__DIR__) # HIDE #
Pkg.add(["CairoMakie","PalmerPenguins", "DataFrames", "AlgebraOfGraphics"]) # HIDE
using PalmerPenguins, DataFrames, AlgebraOfGraphics, CairoMakie
let
    penguins = dropmissing(DataFrame(PalmerPenguins.load()))
    ychi = filter(:species => x -> x == "Chinstrap", penguins)[!,:bill_depth_mm]
    yade = filter(:species => x -> x == "Adelie", penguins)[!,:bill_depth_mm]
    ygen = filter(:species => x -> x == "Gentoo", penguins)[!,:bill_depth_mm]

    penguin_bill = data(penguins) * mapping(:bill_length_mm  => "bill length mm",
        :bill_depth_mm  => "bill depth mm")

    layers = AlgebraOfGraphics.density() * visual(Contour) + linear() +
                visual(alpha = 0.5, markersize = 15)
    aesPoints = penguin_bill * layers * mapping(color = :species, marker = :species,)
    aesHist = data(penguins) * mapping(:bill_length_mm  => "bill length mm",
        color = :species, stack = :species) * AlgebraOfGraphics.histogram(;bins = 20)

    estilo = (color=["#FC7808", "#8C00EC", "#107A78"],
                marker=[:circle, :utriangle, :rect])
    # the actual plot with AoG and Makie
    fig = Figure(; resolution = (700, 700))
    axP = draw!(fig[2,1], aesPoints, palettes = estilo)
    axH = draw!(fig[1,1], aesHist, palettes = estilo)
    ax1 = Axis(axH[1,1])
    ax2 = Axis(axP[1,1])

    ax3 = Axis(fig[2,2],xgridstyle=:dash,ygridstyle=:dash,xtickalign=1, ytickalign=1)
    chin = density!(ax3, ychi, direction = :y, color =("#8C00EC", 0.25),
            strokewidth = 1, strokecolor = "#8C00EC")
    ade = density!(ax3, yade, direction = :y, color =("#FC7808", 0.25),
            strokewidth = 1, strokecolor = "#FC7808")
    gen = density!(ax3, ygen, direction = :y, color =("#107A78", 0.25),
            strokewidth = 1, strokecolor = "#107A78")
    leg = Legend(fig, [chin,ade,gen], ["Chinstrap", "Adelie", "Gentoo"],)
    fig[1,2] = leg
    limits!(ax2, 30, 60, 10, 25)
    xlims!(ax1, 30,60)
    ylims!(ax3, 10, 25)
    hidexdecorations!(ax1, ticks = false, grid = false)
    hideydecorations!(ax3, ticks = false, grid = false)
    hidespines!(ax3, :b,:r,:t)
    hidespines!(ax1, :l,:r,:t)
    colsize!(fig.layout, 1, Relative(3/4))
    rowsize!(fig.layout, 1, Relative(1/4))
    fig
    save(joinpath(@__DIR__, "output", "penguinsAoGsides.png"), fig) # HIDE
end
Pkg.status() # HIDE
Pkg.activate() # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Pkg
using JuMP;
using Plots;
using Distributions
using Random

#a
function simul_a()
    a_1 = 0
    p_1 = 1
    mean = [0.,0.]
    C = [0.2 0; 0 0.3]
    L = 1000
    mu = zeros(L)
    y = zeros(L)

    d = Normal(a_1, p_1)
    mu[1] = rand(d, 1)[1]

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + x[2,1]
    for i in 1:L-1
        mu[i+1] = mu[i] + x[1,i+1]
        y[i+1] = mu[i] + x[2,i+1]
    end

    return y
end
function plot_a()
    p = plot(simul_a())
    plot!(simul_a())
    plot!(simul_a())
    plot!(simul_a())
    plot!(simul_a())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\a.png")
end

plot_a()

#b
function simul_b()
    phi = [0.1, -0.1, 0.1]
    theta = [0.1, -0.2, 0.3, -0.4]
    c = 8
    L = 1000
    y = zeros(L)

    d = Normal(0, 0.1)
    x = rand(d, L*2)

    y[1] = 8
    y[2] = 8.1
    y[3] = 8
    for i in 4:L
        y[i] = c + phi[1]*y[i-1]+ phi[2]*y[i-2] + phi[3]*y[i-3]
        + theta[1]*x[i+3] + theta[2]*x[i+2] + theta[2]*x[i+1] + theta[2]*x[i] + x[i+4]
    end

    return y
end
function plot_b()
    p = plot(simul_b())
    plot!(simul_b())
    plot!(simul_b())
    plot!(simul_b())
    plot!(simul_b())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\b.png")
end

plot_b()

#c
function simul_c()
    phi = [0.1, -0.1, 0.1]
    theta = [0.1, -0.2, 0.3, -0.4]
    c = 8
    L = 1000
    y = zeros(L)

    d = Normal(0, 0.1)
    x = rand(d, L*2)

    y[1] = 8
    y[2] = 8.1
    y[3] = 8
    for i in 4:L
        y[i] = c + (1+phi[1])*y[i-1]+ (phi[2]-phi[2])*y[i-2] - phi[3]*y[i-3]
        + theta[1]*x[i-1] + x[i]
    end

    return y
end
function plot_c()
    p = plot(simul_c())
    plot!(simul_c())
    plot!(simul_c())
    plot!(simul_c())
    plot!(simul_c())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\c.png")
end

plot_c()

#d
function simul_d()
    mean = [0.,0.]
    C = [0.2 0;0 0.1]
    L = 1000
    phi = zeros(L)
    y = zeros(L)
    beta = 0.05
    alpha = 0.2
    phi[1] = 1
    y[1] = 1

    d = MvNormal(mean, C)
    x = rand(d, L)
    for i in 1:L-1
        phi[i+1] = beta + alpha*phi[i] + beta + x[1,i+1]
        y[i+1] = phi[i]*y[i] + x[2,i]
    end

    return y
end
function plot_d()
    p = plot(simul_d())
    plot!(simul_d())
    plot!(simul_d())
    plot!(simul_d())
    plot!(simul_d())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\d.png")
end

plot_d()

#e
function simul_e()
    mean = [0.,0., 0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.25]
    L = 1000
    mu = zeros(L)
    psi = zeros(L)
    y = zeros(L)
    beta = 0.05
    phi = 0.2
    theta = -0.5
    mu[1] = 1
    psi[1] = 0.9

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + psi[1] + x[3,1]
    for i in 1:L-1
        mu[i+1] = mu[i] + beta + x[1,i+1]
        psi[i+1] = phi*psi[i] + theta*x[2,i] + x[2,i+1]
        y[i+1] = mu[i] + psi[i] + x[3,i]
    end

    return y
end
function plot_e()
    p = plot(simul_e())
    plot!(simul_e())
    plot!(simul_e())
    plot!(simul_e())
    plot!(simul_e())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\e.png")
end

plot_e()

#f
function simul_f()
    mean = [0.,0., 0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.25]
    L = 1000
    mu = zeros(L)
    psi = zeros(L)
    y = zeros(L)
    mu[1] = 1
    psi[1] = 0.9

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + psi[1] + x[3,1]
    for i in 1:L-1
        mu[i+1] = mu[i] + x[1,i+1]
        psi[i+1] = psi[i] + x[2,i+1]
        y[i+1] = mu[i] + psi[i] + x[3,i]
    end

    return y
end
function plot_f()
    p = plot(simul_f())
    plot!(simul_f())
    plot!(simul_f())
    plot!(simul_f())
    plot!(simul_f())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\f.png")
end

plot_f()

#g
function simul_g()
    mean = [0.,0., 0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.25]
    L = 1000
    mu = zeros(L)
    beta = zeros(L)
    y = zeros(L)
    mu[1] = 1
    beta[1] = 0.9

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + beta[1] + x[3,1]
    for i in 1:L-1
        beta[i+1] = beta[i] + x[2,i+1]
        mu[i+1] = mu[i] + beta[i] + x[1,i+1]
        y[i+1] = mu[i] + x[3,i]
    end

    return y
end
function plot_g()
    p = plot(simul_g())
    plot!(simul_g())
    plot!(simul_g())
    plot!(simul_g())
    plot!(simul_g())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\g.png")
end

plot_g()

#h
function simul_h()
    mean = [0.,0., 0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.25]
    L = 1000
    mu = zeros(L)
    beta = zeros(L)
    y = zeros(L)
    mu[1] = 1
    beta[1] = 0.9

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + beta[1] + x[3,1]
    for i in 1:L-1
        beta[i+1] = beta[i] + x[2,i+1]
        mu[i+1] = mu[i] + x[1,i+1]
        y[i+1] = mu[i] + x[3,i]
    end

    return y
end
function plot_h()
    p = plot(simul_h())
    plot!(simul_h())
    plot!(simul_h())
    plot!(simul_h())
    plot!(simul_h())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\h.png")
end

plot_h()

#i
function simul_i()
    mean = [0.,0., 0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.25]
    L = 1000
    mu = zeros(L)
    beta = zeros(L)
    y = zeros(L)
    mu[1] = 1
    beta[1] = 0.9

    phi = 0.7
    gamma = -0.9

    d = MvNormal(mean, C)
    x = rand(d, L)
    y[1] = mu[1] + beta[1] + x[3,1]
    for i in 1:L-1
        beta[i+1] = gamma*beta[i] + x[2,i+1]
        mu[i+1] = phi*mu[i] + x[1,i+1]
        y[i+1] = mu[i] + x[3,i]
    end

    return y
end
function plot_i()
    p = plot(simul_i())
    plot!(simul_i())
    plot!(simul_i())
    plot!(simul_i())
    plot!(simul_i())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\i.png")
end

plot_i()

#j
function simul_j()
    mean = [0.,0.]
    C = [0.2 0;0 0.3]
    L = 1000
    mu = zeros(L)
    epsilon = zeros(L)
    y = zeros(L)
    mu[1] = 1

    phi = 0.7

    d = MvNormal(mean, C)
    x = rand(d, L*2)
    # y[1] = mu[1] + beta[1] + x[3,1]
    for i in 1:L-1
        epsilon[i+1] = phi*x[2,i] + x[2,i+1]
        mu[i+1] = mu[i] + x[1,i+1]
        y[i+1] = mu[i] + epsilon[i]
    end

    return y
end
function plot_j()
    p = plot(simul_j())
    plot!(simul_j())
    plot!(simul_j())
    plot!(simul_j())
    plot!(simul_j())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\j.png")
end

plot_j()

#k
function simul_k()
    mean = [0.,0.]
    C = [0.2 0;0 0.3]
    L = 1000
    mu = zeros(L)
    epsilon = zeros(L)
    y = zeros(L)
    mu[1] = 1

    phi = 0.7

    d = MvNormal(mean, C)
    x = rand(d, L*2)
    # y[1] = mu[1] + beta[1] + x[3,1]
    for i in 1:L-1
        epsilon[i+1] = phi*epsilon[i] + x[2,i+1]
        mu[i+1] = mu[i] + x[1,i+1]
        y[i+1] = mu[i] + epsilon[i]
    end

    return y
end
function plot_k()
    p = plot(simul_k())
    plot!(simul_k())
    plot!(simul_k())
    plot!(simul_k())
    plot!(simul_k())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\k.png")
end

plot_k()

#l
function simul_l()
    mean = [0.,0.,0.]
    C = [0.2 0 0;0 0.3 0;0 0 0.1]
    L = 1000
    mu = zeros(L)
    beta = zeros(L+2)
    y = zeros(L)
    mu[1] = 1
    beta[1] = 0.5
    beta[2] = 0.7
    beta[3] = 0.7

    phi = 0.7

    d = MvNormal(mean, C)
    x = rand(d, L*2)
    for i in 1:L-1
        beta[i+3] = 3phi*beta[i+2] - 3phi^2*beta[i+1] + phi^3*beta[i] + x[1,i]
        mu[i+1] = mu[i] + beta[i+3] + x[2,i]
        y[i+1] = mu[i+1] + x[3,i]
    end

    return y
end
function plot_l()
    p = plot(simul_l())
    plot!(simul_l())
    plot!(simul_l())
    plot!(simul_l())
    plot!(simul_l())
    savefig(p,"Modelos Estruturais\\Lista 1\\imagens\\l.png")
end

plot_l()

max_stars_repo_path: Modelos Estruturais/Lista 1/simul.jl

max_issues_repo_path: Modelos Estruturais/Lista 1/simul.jl

max_forks_repo_path: Modelos Estruturais/Lista 1/simul.jl


text:
module EasyGit

using Dates

export gitstaginginfo, gitisclean, githead, gitshorthead,
    githeadtime, gitroot

function gitverify(repfolder::String=".")
    try
        read(`command -v git`)
    catch
        error("git is not installed.")
    end
    try
        read(`git -C $repfolder rev-parse --is-inside-work-tree`)
    catch
        error("Not a git repository.")
    end
    true
end

"""
    gitstaginginfo(repfolder::String=".")
Get modified and untracked files of a git repository as String vectors.
"""
function gitstaginginfo(repfolder::String=".")
    gitverify(repfolder)
    status = split(readchomp(`git -C $repfolder status --porcelain`), "\n")
    modified = String[]
    untracked = String[]
    for entry in status
        if length(entry) < 2
            continue
        elseif entry[1:2] == " M"
            push!(modified, entry[4:end])
        elseif entry[1:2] == "??"
            push!(untracked, entry[4:end])
        end
    end
    modified, untracked
end

"""
    gitisclean(repfolder::String=".")
Is true if there are no modified and no untracked files in the repository that are not
yet committed (or ignored).
"""
function gitisclean(repfolder::String=".")
    modified, untracked = gitstaginginfo(repfolder)
    isempty(modified) && isempty(untracked)
end

"""
    githead(repfolder::String=".")
Get hash of the repository head.
"""
function githead(repfolder::String=".")
    gitverify(repfolder)
    readchomp(`git -C $repfolder rev-parse HEAD`)
end

"""
    githead(repfolder::String=".")
Get short hash of the repository head.
"""
function gitshorthead(repfolder::String=".")
    hid = githead(repfolder)
    hid[1:7]
end

"""
    githeadtime(repfolder::String=".")
Get date and time of the last commit.
"""
function githeadtime(repfolder::String=".")
    gitverify(repfolder)
    htime_unix = readchomp(`git show -s --format=%ct HEAD`)
    unix2datetime(parse(Int64, htime_unix))
end

"""
    gitroot(repfolder::String=".")
Get root folder of the repository.
"""
function gitroot(repfolder::String=".")
    gitverify(repfolder)
    readchomp(`git rev-parse --show-toplevel`)
end

end


text:
"""
    BurgersPeriodicSemidiscretization(D, Di, split_form)

A semidiscretization of Burgers' equation
    \$\\partial_t u(t,x) + \\partial_x \\frac{u(t,x)^2}{2} = 0\$
with periodic boundary conditions.

`D` is a first-derivative SBP operator, `Di` an associated dissipation operator
or `nothing`, and `split_form::Union{Val(true), Val(false)}` determines whether
the canonical split form or the conservative form is used.
"""
struct BurgersPeriodicSemidiscretization{T,Derivative<:AbstractDerivativeOperator{T},
                                            Dissipation,
                                            SplitForm<:Union{Val{false}, Val{true}}} <: AbstractSemidiscretization
    derivative::Derivative
    dissipation::Dissipation
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm

    function BurgersPeriodicSemidiscretization(derivative::Derivative, dissipation::Dissipation, split_form::SplitForm=Val{false}()) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}}
        if dissipation !== nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm}(derivative, dissipation, tmp1, tmp2, split_form)
    end
end

function Base.show(io::IO, semi::BurgersPeriodicSemidiscretization)
    if get(io, :compact, false)
        print(io, "Semidiscretization of Burgers' equation (periodic)")
    else
        print(io, "Semidiscretization of Burgers' equation\n")
        print(io, "  \$ \\partial_t u(t,x) + \\partial_x \\frac{u(t,x)^2}{2} = 0 \$ \n")
        print(io, "with periodic boundaries using")
        if semi.split_form == Val{true}()
            print(io, " a split form and: \n")
        else
            print(io, " no split form and: \n")
        end
        print(io, semi.derivative)
        print(io, semi.dissipation)
    end
end

function (disc::BurgersPeriodicSemidiscretization)(du, u, p, t)
    @unpack tmp1, tmp2, derivative, dissipation, split_form = disc
    @boundscheck begin
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        m1_3 = -one(eltype(u)) / 3

        ## u * D * u
        mul!(tmp1, derivative, u)
        @. du = m1_3 * u * tmp1
        ## D * u^2
        @. tmp2 = u^2
        mul!(tmp1, derivative, tmp2)
        @. du += m1_3 * tmp1
    else
        m1_2 = -one(eltype(u)) / 2
        @. tmp2 = u^2
        mul!(tmp1, derivative, tmp2)
        @. du = m1_2 * tmp1
    end

    # dissipation
    if dissipation !== nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    nothing
end

"""
    BurgersNonperiodicSemidiscretization(D, Di, split_form, left_bc, right_bc)

A semidiscretization of Burgers' equation
    \$\\partial_t u(t,x) + \\partial_x \\frac{u(t,x)^2}{2} = 0\$
with boundary conditions `left_bc(t)`, `right_bc(t)`.

`D` is a first-derivative SBP operator, `Di` an associated dissipation operator
or `nothing`, and `split_form::Union{Val(true), Val(false)}` determines whether
the canonical split form or the conservative form is used.
"""
struct BurgersNonperiodicSemidiscretization{T,Derivative<:AbstractDerivativeOperator{T},
                                            Dissipation,
                                            SplitForm<:Union{Val{false}, Val{true}},
                                            LeftBC, RightBC} <: AbstractSemidiscretization
    derivative::Derivative
    dissipation::Dissipation
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm
    left_bc::LeftBC
    right_bc::RightBC

    function BurgersNonperiodicSemidiscretization(derivative::Derivative, dissipation::Dissipation, split_form::SplitForm, left_bc::LeftBC, right_bc::RightBC) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}, LeftBC, RightBC}
        if dissipation !== nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm,LeftBC,RightBC}(derivative, dissipation, tmp1, tmp2, split_form, left_bc, right_bc)
    end
end

function Base.show(io::IO, semi::BurgersNonperiodicSemidiscretization)
    if get(io, :compact, false)
        print(io, "Semidiscretization of Burgers' equation (non-periodic)")
    else
        print(io, "Semidiscretization of Burgers' equation\n")
        print(io, "  \$ \\partial_t u(t,x) + \\partial_x \\frac{u(t,x)^2}{2} = 0 \$ \n")
        print(io, "with nonperiodic boundaries using")
        if semi.split_form == Val{true}()
            print(io, " a split form and: \n")
        else
            print(io, " no split form and: \n")
        end
        print(io, semi.derivative)
        print(io, semi.dissipation)
    end
end

function godunov_flux_burgers(uₗ::T, uᵣ::T) where {T<:Real}
    if uₗ < uᵣ
        if uₗ < 0 && 0 < uᵣ
            zero(T)
        else
            min(uₗ^2/2, uᵣ^2/2)
        end
    else
        max(uₗ^2/2, uᵣ^2/2)
    end
end

function (disc::BurgersNonperiodicSemidiscretization)(du, u, p, t)
    @unpack tmp1, tmp2, derivative, dissipation, split_form, left_bc, right_bc = disc
    @boundscheck begin
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        m1_3 = -one(eltype(u)) / 3

        ## u * D * u
        mul!(tmp1, derivative, u)
        @. du = m1_3 * u * tmp1
        ## D * u^2
        @. tmp2 = u^2
        mul!(tmp1, derivative, tmp2)
        @. du += m1_3 * tmp1
    else
        m1_2 = -one(eltype(u)) / 2
        @. tmp2 = u^2
        mul!(tmp1, derivative, tmp2)
        @. du = m1_2 * tmp1
    end

    # dissipation
    if dissipation !== nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    # boundary conditions via Godunov's flux
    @inbounds fnum_left = godunov_flux_burgers(left_bc(t), u[1])
    @inbounds du[1] += (fnum_left - u[1]^2/2) / left_boundary_weight(derivative)
    @inbounds fnum_right = godunov_flux_burgers(u[end], right_bc(t))
    @inbounds du[end] -= (fnum_right - u[end]^2/2) / right_boundary_weight(derivative)

    nothing
end


text:
import Pkg
import UUIDs

function git_clone(tmp_dir, previous_directory, url, name)
    cd(tmp_dir)
    run(`git clone $(url) $(name)`)
    cd(previous_directory)
end

function download_or_clone(tmp_dir, previous_director, reg_url, registry_path, url, name)
    try
        Pkg.Types.download_verify_unpack(reg_url, nothing, registry_path, ignore_existence = true)
    catch err
        # in case of fail, fallback to clone
        git_clone(tmp_dir, previous_director, url, name)
    end
end

function get_latest_version_from_registries!(dep_to_latest_version::Dict{Package, Union{VersionNumber, Nothing}},
                                             registry_list::Vector{Pkg.Types.RegistrySpec};
                                             use_pkg_server::Bool=false)
    original_directory = pwd()
    num_registries = length(registry_list)
    registry_temp_dirs = Vector{String}(undef, num_registries)
    registry_urls = nothing
    for (i, registry) in enumerate(registry_list)
        tmp_dir = mktempdir()
        atexit(() -> rm(tmp_dir; force = true, recursive = true))
        registry_temp_dirs[i] = tmp_dir
        name = registry.name
        url = registry.url
        uuid = registry.uuid
        previous_directory = pwd()
        if use_pkg_server
            registry_path = joinpath(tmp_dir, name)
            @static if VERSION >= v"1.5.0"
                # copied from https://github.com/JuliaLang/Pkg.jl/blob/release-1.5/src/Types.jl#L981
                reg_url, registry_urls = Pkg.Types.pkg_server_registry_url(uuid, registry_urls)
            elseif v"1.4.0" <= VERSION < v"1.5.0"
                # copied from https://github.com/JuliaLang/Pkg.jl/blob/v1.4.2/src/Types.jl#L921
                reg_url = Pkg.Types.pkg_server_registry_url(uuid)
            else
                reg_url = nothing
            end

            if reg_url !== nothing
                download_or_clone(tmp_dir, previous_directory, reg_url, registry_path, url, name)
            elseif url !== nothing # clone from url
                git_clone(tmp_dir, previous_directory, url, name)
            else
                error("Could not download registry, nor pkg server not git download works.")
            end
        else
            git_clone(tmp_dir, previous_directory, url, name)
        end
    end
    for (registry_temp_dir, registry) in zip(registry_temp_dirs, registry_list)
        previous_directory = pwd()
        name = registry.name
        registry_path = joinpath(registry_temp_dir, name)
        cd(registry_path)
        registry_parsed = Pkg.TOML.parsefile(joinpath(registry_path, "Registry.toml"))
        packages = registry_parsed["packages"]
        for p in packages
            name = p[2]["name"]
            uuid = UUIDs.UUID(p[1])
            package = Package(name, uuid)
            path = p[2]["path"]
            if package in keys(dep_to_latest_version)
                versions = VersionNumber.(collect(keys(Pkg.TOML.parsefile(joinpath(registry_path, path, "Versions.toml")))))
                old_value = dep_to_latest_version[package]
                if isnothing(old_value)
                    dep_to_latest_version[package] = maximum(versions)
                else
                    dep_to_latest_version[package] = max(old_value, maximum(versions))
                end
            end
        end
        cd(previous_directory)
    end
    cd(original_directory)
    for tmp_dir in registry_temp_dirs
        rm(tmp_dir; force = true, recursive = true)
    end
    return dep_to_latest_version
end


text:
export NonInteractiveEnv

abstract type NonInteractiveEnv <: AbstractEnv end
(env::NonInteractiveEnv)() = env(nothing)
RLBase.action_space(::NonInteractiveEnv) = (nothing,)

include("pendulum.jl")


text:
@testset "Contains" begin
    @test_throws TypeError(:typeassert, "", Bool, 1) Contains("a", case_sensitive=1)
end

@testset "Fragment" begin
    @test_throws TypeError(:typeassert, "", Bool, 1) Fragment("a", case_sensitive=1)
end

@testset "SearchF" begin
    @test_throws TypeError(:typeassert, "", ExtractionFunction, 1) SearchF("dimensions", Contains("a"), extractionFn=1)
end

@testset "Like" begin
    @test_throws TypeError(:typeassert, "", AbstractString, 1) Like("dimension", "pattern", escape=1)
    @test_throws TypeError(:typeassert, "", ExtractionFunction, 1) Like("dimension", "pattern", extractionFn=1)
end

@testset "Bound" begin
    @test_throws TypeError(:typeassert, "", String, 1) Bound("dimension", lower=1)
    @test_throws TypeError(:typeassert, "", String, 1) Bound("dimension", upper=1)
    @test_throws ErrorException("At least one of lower and upper must be specified") Bound("diimension")
    @test_throws TypeError(:typeassert, "", Bool, 1) Bound("dimension", lower="a", lowerStrict=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) Bound("dimension", lower="a", upperStrict=1)
    @test_throws TypeError(:typeassert, "", String, 1) Bound("dimension", lower="a", ordering=1)
    @test_throws TypeError(:typeassert, "", ExtractionFunction, 1) Bound("dimension", lower="a", extractionFn=1)
end

@testset "IntervalF" begin
    @test_throws TypeError(:typeassert, "", ExtractionFunction, 1) IntervalF("dimension", [Interval("a", "b")], extractionFn=1)
end


text:
using Documenter, Compose
import Cairo

struct SVGJSWritable{T}
    x :: T
end
Base.show(io::IO, m::MIME"text/html", x::SVGJSWritable) = show(io, m, x.x)

makedocs(
    modules = [Compose],
    clean = false,
    sitename = "Compose.jl",
    pages = Any[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Gallery" => Any[
            "Forms" => "gallery/forms.md",
            "Properties"=> "gallery/properties.md",
            "Transformations"=> "gallery/transforms.md",
            ],
        "Library" => "library.md"
    ]
)

deploydocs(
    repo   = "github.com/GiovineItalia/Compose.jl.git",
)


text:
include("/Users/ksb/Catlab.jl/src/atp/ATP.jl")
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Theories
using Catlab.Programs
using Catlab.Present

"""
Co-diagonals only are created when they are used in an output, so sometimes we
must create 'fake' outputs and then trim them away (replacing with a junction).
We keep the first `O` outputs and `I` inputs, defaulting to keeping all inputs
and no outputs.
"""
function trim(p::WiringDiagram{BiproductCategory, Any, Any, Any}, O::Int=0,
              I::Int=-1)
  pd = p.diagram
  I = I < 0 ? nparts(pd, :OuterInPort) : I # default to keeping all inputs
  N_O = nparts(pd, :OuterOutPort)
  rem_ows, rem_pws = Int[], Int[]
  for oop in O+1:N_O
    typ = pd[:outer_out_port_type][oop]
    b = add_part!(pd, :Box, value=typ, box_type=Junction{nothing, Symbol})
    for pw in incident(pd, oop, :pass_tgt)
      ip = add_part!(pd, :InPort, in_port_box=b, in_port_type=typ)
      add_part!(pd, :InWire, in_src=pd[:pass_src][pw], in_tgt=ip,
                in_wire_value=nothing)
      push!(rem_pws, pw)
    end
    for ow in incident(pd, oop, :out_tgt)
      ip = add_part!(pd, :InPort, in_port_box=b, in_port_type=typ)
      add_part!(pd, :Wire, src=pd[:out_src][ow], tgt=ip, wire_value=nothing)
      push!(rem_ows, ow)
    end
  end
  rem_parts!(pd, :OutWire, sort(rem_ows))
  rem_parts!(pd, :PassWire, sort(rem_pws))
  rem_parts!(pd, :OuterOutPort, O+1:N_O)
  # Repeat, swapping in for out
  N_I = nparts(pd, :OuterInPort)
  rem_iws, rem_pws = Int[], Int[]
  for oip in I+1:N_I
    typ = pd[:outer_in_port_type][oip]
    b = add_part!(pd, :Box, value=typ, box_type=Junction{nothing, Symbol})
    for pw in incident(pd, oip, :pass_src)
      op = add_part!(pd, :OutPort, out_port_box=b, out_port_type=typ)
      add_part!(pd, :OutWire, out_tgt=pd[:pass_src][pw], out_src=op,
                out_wire_value=nothing)
      push!(rem_pws, pw)
    end
    for iw in incident(pd, oip, :in_src)
      op = add_part!(pd, :OutPort, out_port_box=b, out_port_type=typ)
      add_part!(pd, :Wire, tgt=pd[:in_tgt][iw], src=op, wire_value=nothing)
      push!(rem_iws, iw)
    end
  end
  rem_parts!(pd, :InWire, sort(rem_iws))
  rem_parts!(pd, :PassWire, sort(rem_pws))
  rem_parts!(pd, :OuterInPort, I+1:N_I)
  return p
end

# Single sorted inputs/outputs
Zero, One, Two, Three, Four, Five, Six = [repeat([:X], i) for i in 0:6]

# Boxes
mmul = Box(:mul, Two,   One);
e  = Box(:e,   Zero,  One);
inv  = Box(:inv, One,   One);
gens = Box(:s_,   Zero,  One);
genr = Box(:r_,   Zero,  One);
R  = Box(:R,   One,   One);
f  = Box(:f,   One,   One);
g  = Box(:g,   One,   One);
s  = Box(:s,   [:A],  [:O]);
t  = Box(:t,   [:A],  [:O]);
Id   = Box(:Id,  [:O],  [:A]);
cmp  = Box(:cmp, [:A,:A], [:A]);

o_o  = Box(:o_o, [:O,:O], [:O]); # ⊗ₒ
o_a  = Box(:o_a, [:A,:A], [:A]); # ⊗ₐ

ounit= Box(:⊤, Symbol[], [:O])
σ  = Box(:σ,  [:O, :O], [:A])
del  = Box(:δ,  [:O],   [:A])
eps  = Box(:ϵ,  [:O],   [:A])
mu   = Box(:μ,  [:O],   [:A])
ev   = Box(:ev,   [:O, :O], [:A])
Exp   = Box(:exp, [:O, :O], [:O])

lam   = Box(:lam, [:O, :O, :O, :A], [:A])

ϵ(x::Symbol=:X)   = Junction(x, 1, 0)
η(x::Symbol=:X)   = Junction(x, 0, 1)
δ(x::Symbol=:X)   = Junction(x, 1, 2)
μ(x::Symbol=:X)   = Junction(x, 2, 1)
dot(x::Symbol=:X) = Junction(x, 1, 1)
cap(x::Symbol=:X) = Junction(x, 2, 0)
δ_2(x::Symbol=:X)   = Junction(x, 2, 2)

# Generator sets
Σ0     = Set{Box{Symbol}}()
Σ_monoid   = Set([mmul, e]);
Σ_group  = Set([inv]);
Σ_dihedral = Set([gens, genr]);
Σ_rel    = Set([R, f, g]);
Σ_reflG  = Set([s, t, Id]);
Σ_cat    = Set([cmp]);
Σ_mc     = Set([o_o, o_a, ounit])
Σ_smc    = Set([σ])
Σ_crc    = Set([del, eps])
Σ_dcr    = Set([mu])
Σ_ccc    = Set([ev, Exp, lam])

@present C(FreeBiproductCategory) begin
  X::Ob
  A::Ob
  O::Ob
  mul::Hom(otimes(X, X), X)
  e::Hom(munit(), X)
  inv::Hom(X,X)
  s_::Hom(munit(), X)
  r_::Hom(munit(), X)
  R::Hom(X, X)
  f::Hom(X,X)
  g::Hom(X,X)
  s::Hom(A,O)
  t::Hom(A,O)
  Id::Hom(O,A)
  cmp::Hom(otimes(A, A),A)
  o_o::Hom(otimes(O, O), O)
  o_a::Hom(otimes(A, A), A)

  ⊤::Hom(munit(), O)
  σ::Hom(otimes(O,O), A)
  δ::Hom(O,A)
  ϵ::Hom(O,A)
  μ::Hom(O,A)
  ev::Hom(otimes(O,O), A)
  exp::Hom(otimes(O,O), O)
  lam::Hom(otimes(O, O, O, A), A)
end

X, O, A = (generator(C, name) for name in [:X,:O,:A])
mul_,e_,lam_ = (generator(C, name) for name in [:mul, :e, :lam])

passx = @program C (x::X) -> x
passo = @program C (x::O) -> x
passa = @program C (x::A) -> x

epass =  trim(@program(C, (x::X,y::X), (x, [y,e()])), 1)

e11 = @program C (x::X) (e(),x)
ll = @program(C, (x::X), let y=e(); (y, mul(y, x)) end)
leftid = Eq(:leftid, ll, e11, false);

rr = @program(C, (x::X), let y=e(); (y, mul(x, y)) end)
rightid = Eq(:rightid, rr, e11, false);

ma1 = @program C (x::X,y::X,z::X) -> mul(x, mul(y, z))
ma2 = @program C (x::X,y::X,z::X) -> mul(mul(x, y), z)
mul_assoc = Eq(:mul_assoc, ma1, ma2);

ma1 = @program C (x::A,y::A,z::A) -> cmp(x, cmp(y, z))
ma2 = @program C (x::A,y::A,z::A) -> cmp(cmp(x, y), z)
cmp_assoc = Eq(:cmp_assoc, ma1, ma2);

ma1 = @program C (x::O,y::O,z::O) -> o_o(x, o_o(y, z))
ma2 = @program C (x::O,y::O,z::O) -> o_o(o_o(x, y), z)
o_o_assoc = Eq(:o_o_assoc, ma1, ma2);

ma1 = @program C (x::A,y::A,z::A) -> o_a(x, o_a(y, z))
ma2 = @program C (x::A,y::A,z::A) -> o_a(o_a(x, y), z)
o_a_assoc = Eq(:o_a_assoc, ma1, ma2);

e2x = trim(@program(C, (x::X), let y=e(); (y, y, x) end), 2)
rightinv = Eq(:rightinv, @program(C, (x::X), (e(),mul(x, inv(x)))), e2x, false);
leftinv  = Eq(:leftinv,  @program(C, (x::X), (e(),mul(inv(x),x))), e2x, false);
if 1+1==2
posl = trim(@program(C, (x::X,y::X), [mul(x,y), e()]))
posr = trim(@program(C, (x::X,y::X), [x, e(), y]))
pos = Eq(:pos,  posl, posr, false)

cancell = trim(@program(C, (x::X,y::X), (e(), [mul(x,y), y])), 1)
cancelr = @program(C, (x::X,y::X), [x, e()])
cancel = Eq(:cancel, cancell, cancelr, false)

sym = Eq(:sym, @program(C, (x::X,y::X) -> mul(x,y)),
               @program(C, (x::X,y::X) -> mul(y,x)))

uniq_inv_1 = trim(@program(C,(a::X,b::X,c::X),[mul(a, b),mul(b, c),e()]))
uniq_inv_2 = trim(@program(C,(a::X,b::X,c::X),[a,c]))
uniq_inv   = Eq(:uniq_inv,   uniq_inv_1,   uniq_inv_2);

div_1 =  @program C (a::X,b::X,c::X) [mul(a,b),c]
div_2 =  @program C (a::X,b::X,c::X) [mul(inv(a),c),b]
gdiv     = Eq(:gdiv,     div_1,    div_2);

leftcancel_1 = @program C (a::X,b::X,c::X) [mul(a, b), mul(a, c)]
leftcancel_2 = @program C (a::X,b::X,c::X) [b, c]
leftcancel   = Eq(:leftcancel, leftcancel_1, leftcancel_2);

e_wd = @program C () e()

ss = @program C () let x=s_(); mul(x,x) end
s_order_2 = Eq(:s_ord_2, ss, e_wd);

rrr = @program C () let x=r_(); mul(x,mul(x,x)) end
r_order_3 = Eq(:r_ord_3,  rrr, e_wd);

srs = @program C () let x=r_(), y=s_(); mul(mul(y, x), inv(y)) end
rinv = @program C () inv(r_())
sr2 =  @program C () let x=mul(s_(), r_()); mul(x,x) end

R_copy = @program C (i::X) let x=R(i); (x,x) end
fig26 =  @program C (i::X) let x=g(f(i)); (x,x) end

refl_s = @program C (a::O) s(Id(a))
refls = Eq(:refls, refl_s, passo);

refl_t = @program C (a::O) t(Id(a))
reflt = Eq(:reflt, refl_t, passo);

cmp_1 = trim(@program(C,(ff::A,gg::A),(s(ff),[t(ff),s(gg)],t(gg))))

cmp_2 = trim(@program(C,(ff::A,gg::A),begin
  cm = cmp(ff,gg)
  ([s(cm),s(ff)],[t(gg),t(cm)])
  end ))
cmp_intro = Eq(:cmp, cmp_1, cmp_2, false);

id_s = @program C (a::A) cmp(Id(s(a)),a)
l_unit = Eq(:l_unit, id_s, passa);

id_t = @program C (a::A) cmp(a,Id(t(a)))
r_unit = Eq(:r_unit, id_t, passa);

o_func_sl = @program C (a::A, b::A) s(o_a(a,b))
o_func_sr = @program C (a::A, b::A) o_o(t(a),t(b))
o_func_s = Eq(:o_func_s,o_func_sl,o_func_sr)

o_func_tl = @program C (a::A, b::A) t(o_a(a,b))
o_func_tr = @program C (a::A, b::A) o_o(t(a),t(b))
o_func_t = Eq(:o_func_t,o_func_tl,o_func_tr)

o_f_c_1 =  @program C (a::A,b::A,c::A,d::A) cmp(o_a(a,b),o_a(c,d))
o_f_c_2 =  @program C (a::A,b::A,c::A,d::A) o_a(cmp(a,c),cmp(b,d))
o_func_cmp = Eq(:o_func_cmp, o_f_c_1,    o_f_c_2);

lunit_o_a_ = @program C (a::A) o_a(Id(⊤()), a)
lunit_o_a  = Eq(:lunit_o_a,  lunit_o_a_,   passa);

runit_o_a_ = @program C (a::A) o_a(a, Id(⊤()))
runit_o_a  = Eq(:runit_o_a,  runit_o_a_,   passa);

lunit_o_o_ = @program C (a::O) o_o(⊤(), a)
lunit_o_o  = Eq(:lunit_o_o,  lunit_o_o_,   passo);

runit_o_o_ = @program C (a::O) o_o(a, ⊤())
runit_o_o  = Eq(:runit_o_o,  runit_o_o_,   passo);

bb_1 = @program C (a::O,b::O) cmp(σ(a,b),σ(b,a))
bb_2 = @program C (a::O,b::O) Id(o_o(a,b))
braidbraid = Eq(:braidbraid, bb_1,     bb_2);

braid_tl = @program C (a::O,b::O) t(σ(a,b))
braid_tr = @program C (a::O,b::O) o_o(b,a)
braid_t = Eq(:braid_t, braid_tl, braid_tr)

braid_sl = @program C (a::O,b::O) s(σ(a,b))
braid_sr = @program C (a::O,b::O) o_o(a,b)
braid_s = Eq(:braid_s, braid_sl, braid_sr)

braid_ots_t = @program C (a::A,b::A) cmp(o_a(a,b), σ(t(a),t(b)))
braid_ots_s = @program C (a::A,b::A) cmp(σ(s(a),s(b)),o_a(b,a))
braid_ots = Eq(:braid_ots, braid_ots_t, braid_ots_s);

eps_s_l = @program C (a::O) s(ϵ(a))
eps_s   = Eq(:eps_s, eps_s_l, passo)

del_s_l =@program C (a::O) s(δ(a))
del_s   = Eq(:del_s, del_s_l, passo)

e_t = @program C (a::O) t(ϵ(a))
eps_t_2 = @program C (_::O) ⊤()
eps_t   = Eq(:eps_o,   e_t, eps_t_2);

d_t = @program C (a::O) t(δ(a))
del_t_2 = @program C (a::O) o_o(a,a)
del_t   = Eq(:del_t,   d_t, del_t_2);

eps_coh_1 = @program C (a::O,b::O) ϵ(o_o(a,b))
eps_coh_2 = @program C (a::O,b::O) o_a(ϵ(a),ϵ(b))
eps_coh = Eq(:eps_coh, eps_coh_1, eps_coh_2);

del_coh_1 = @program C (a::O,b::O) δ(o_o(a,b))
del_coh_2 = @program C (a::O,b::O) [o_a(δ(a),δ(b)),o_a(o_a(Id(a),σ(a,b)),Id(b))]
del_coh = Eq(:del_coh, del_coh_1, del_coh_2);

del_nat_1 = @program C (a::O) cmp(δ(a),σ(a,a))
del_nat_2 =@program C (a::O) δ(a)
del_nat = Eq(:del_nat, del_nat_1, del_nat_2);

cc1_tt = @program C (a::O) let x=δ(a); cmp(x,o_a(x,Id(a))) end
cc1_ft = @program C (a::O) let x=δ(a); cmp(x,o_a(Id(a),x)) end
cc1   = Eq(:cc1, cc1_tt, cc1_ft);

cc1_tf = @program C (a::O) cmp(δ(a),o_a(ϵ(a),Id(a)))
cc1_ff = @program C (a::O) cmp(δ(a),o_a(Id(a),ϵ(a)))
cc2   = Eq(:cc2, cc1_tf, cc1_ff);

cc3_1 = @program C (a::A) cmp(δ(s(a)), o_a(a,a))
cc3_2 = @program C (a::A) cmp(a, δ(t(a)))
cc3   = Eq(:cc3, cc3_1, cc3_2);

delmu_t = @program C (a::O) cmp(δ(a), μ(a))
idbox = @program C (a::O) Id(a)
bone  = Eq(:bone, delmu_t, idbox)

proj1   = Eq(:mu_s, @program(C,(a::O),s(μ(a))), passo)
proj2   = Eq(:mu_t, @program(C,(a::O),t(μ(a))), passo)

frobr = @program C (a::O) cmp(μ(a),δ(a))
frobll = @program C (a::O) let x=Id(a); cmp(o_a(x,δ(a)),o_a(μ(a),x)) end
frob1   = Eq(:frob_l, frobll, frobr)

frobrl = @program C (a::O) let x=Id(a); cmp(o_a(δ(a),x),o_a(x,μ(a))) end
frob2   = Eq(:frob_r, frobrl, frobr)

epsnat1 = @program C (a::A) cmp(a,ϵ(t(a)))
epsnat2 = @program C (a::A) ϵ(s(a))
eps_nat = Eq(:epsnat, epsnat1, epsnat2)

ev_s = @program C (a::O,b::O) s(ev(a,b))
evs2 = @program C (a::O,b::O) o_o(exp(a,b),a)
evs   = Eq(:evs, ev_s, evs2)

ev_t = @program C (a::O,b::O) t(ev(a,b))
evt2 = @program C (_::O,b::O) (b,)
evt   = Eq(:evt, ev_t, evt2)

lam_intro1 = trim(@program C (x::O,y::O,z::O,q::O) ([o_o(x,y),s(q)],[z,t(q)]))
lam_intro2 = trim(@program C (x::O,y::O,z::O,q::O) lam(x,y,z,q))
λ_intro = Eq(:λ_intro, lam_intro1,  lam_intro2, false)

lam_s_ =  @program C (a::O,b::O,c::O,d::A) s(lam(a,b,c,d))
lam_s2 = @program C (x::O,_::O,_::O,_::A) (x,)
lam_s   = Eq(:lam_s,   lam_s_,  lam_s2)

lam_t_ =  @program C (a::O,b::O,c::O,d::A) t(lam(a,b,c,d))
lam_t2 = @program C (_::O,x::O,y::O,_::A) exp(x,y)
lam_t   = Eq(:lam_t,   lam_t_, lam_t2)

lam_eqf1 = @program C (a::O,b::O,c::O,d::A) cmp(o_a(lam(a,b,c,d),Id(b)),ev(b,c))
lam_eqf2 = trim(@program(C, (a::O,b::O,c::O,d::A), (d, [o_o(a,b),s(d)],[c,t(d)])), 1)
lam_eqf = Eq(:lam_eqf, lam_eqf1,    lam_eqf2)

lam_eqg1 = trim(@program(C, (a::O,b::O,c::O,d::A), (d, [a,s(d)], [exp(b,c),t(d)])), 1)
lam_eqg2 =  @program C (a::O,b::O,c::O,d::A) lam(a,b,c,cmp(o_a(d, Id(b)), ev(b,c)))
lam_eqg = Eq(:lam_eqg, lam_eqg1,    lam_eqg2)

# Equation sets
I_monoid = Set([mul_assoc, leftid, rightid]);

I_group  = Set([leftinv, rightinv]);
I_d3h  = Set([s_order_2, r_order_3]);
I_reflG  = Set([refls, reflt])
I_cat  = Set([cmp_assoc, l_unit, r_unit])
I_mc   = Set([o_o_assoc, o_a_assoc, o_func_s, o_func_t, o_func_cmp,
        lunit_o_a, runit_o_a, lunit_o_o, runit_o_o])
I_smc  = Set([braidbraid, braid_t, braid_s, braid_ots])
I_crc  = Set([eps_s, del_s, eps_t, del_t, eps_coh, del_coh, del_nat,
        cc1, cc2, cc3])
I_dcr  = Set([bone, proj1, proj2, frob1, frob2])
I_cc   = Set([eps_nat])
I_ccc  = Set([evs, evt, λ_intro, lam_s, lam_t, lam_eqf, lam_eqg])
I_pcs_mon = Set([mul_assoc, leftid, pos,cancel,sym])

# Theories
T_monoid = EqTheory(Σ_monoid, I_monoid);

T_pcs_mon = EqTheory(Σ_monoid, I_pcs_mon);
T_group  = union(T_monoid, Σ_group,  I_group);
T_d3h  = union(T_group,  Σ_dihedral, I_d3h);

T_reflG  = EqTheory(Σ_reflG, I_reflG);
T_cat  = union(T_reflG, Σ_cat, I_cat);
T_mc   = union(T_cat,   Σ_mc,  I_mc);
T_smc  = union(T_mc,  Σ_smc, I_smc);
T_crc  = union(T_smc,   Σ_crc, I_crc);
T_dcr  = union(T_crc,   Σ_dcr, I_dcr);
T_cc   = union(T_dcr,   Σ0,  I_cc)
T_ccc  = union(T_cc,  Σ_ccc, I_ccc)
end


text:
"""
    writeMps!

Write the deterministic equivalent formulation in MPS format.
"""
function writeMps!(m::SJ.StructuredModel, filename::AbstractString; options...)
    # free any existing model pointer
    freeModel(dspenv)

    # set options
    setoptions!(options)

    # classify the second stage constraints into linear ones and quadratic ones
    classifyConstrs(m)

    # load problem
    load_problem!(m)

    # write problem into MPS file
    writeMps!(dspenv, filename)

    return
end


text:
VERSION <= v"0.5.0" && __precompile__()

module DataTrees

import Base: start, next, done, get, getindex, print
# import isnull?
export TreeNode, addchild!, addsubtree!, buildtree, print, getroot, getdepth, getindex, getsiblings, isleaf, isroot, traverse

include("TreeNode.jl")
include("indexing.jl")
include("buildtree.jl")
include("printing.jl")
include("traversing.jl")

# interface for AbstractTrees
children(tn::TreeNode) = tn.children
start(tn::TreeNode) = start(tn.children)
next(tn::TreeNode, state) = next(tn.children, state)
done(tn::TreeNode, state) = done(tn.children, state)

# TODO: all relevant methods also for Nullable{TreeNode}?

end # module


text:
using ArgParse, MAT

include("BBTNet/models/retinaface.jl")
include("./BBTNet/utils/ImageReader.jl")
include("configs.jl")

function parse_cmd()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--load_dir"
            help = "If there is a pretrained model, then the path of it."
            arg_type = String
            required = true
        "--mode", "-m"
            help = "Training mode: 0 for only baseline until FPN, 1 for full model, 2 for no cascaded structure."
            arg_type = Int
            default = 1
            required = false
        "--laterals", "-l"
            help = "How many lateral connections will be processed, either 3 or 5. 5 is needed for full model."
            arg_type = Int
            default = 5
            required = false
        "--eval_dir"
            help = "WIDER FACE Easy/Medium/Hard .mat files directory path.."
            arg_type = String
            required = true
        "--save_dir"
            help = "Directory to save all prediction files."
            arg_type = String
            default = "./data/wider_txt/"
        "--nms"
            help = "NMS threshold, for eliminating redundant boxes for computational effiency."
            arg_type = Float64
            default = 0.4
        "--conf"
            help = "Confidence threshold, for eliminating really low probability boxes for computational effiency."
            arg_type = Float64
            default = 0.02
        "--topk"
            help = "Keep top k many elements after confidence filtering."
            arg_type = Int
            default = 5000
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_cmd()
    nms_threshold = parsed_args["nms"]; conf_threshold = parsed_args["conf"];
    topk = parsed_args["topk"]
    mode = parsed_args["mode"]; scale_cnt = parsed_args["laterals"];
    num_anchors = scale_cnt == 3 ? 2 : 3
    anchor_info = scale_cnt == 3 ? lat_3_anchors : lat_5_anchors

    load_path = parsed_args["load_dir"]; eval_dir = parsed_args["eval_dir"];
    save_dir = parsed_args["save_dir"]; mkdir(save_dir);
    model = RetinaFace(
        mode=mode, num_anchors=num_anchors, anchor_info=anchor_info, load_path=load_path, dtype=atype
    )
    print("[INFO] Model is loaded!\n")

    val_data = matread(eval_dir * "wider_face_val.mat")
    events = val_data["event_list"]; esize = size(events, 1); files = val_data["file_list"];

    iter_no = 0;
    for e in 1:esize
        event = events[e, 1]
        for f in 1:length(files[e])
            file = files[e][f]
            path = event * "/" * file * ".jpg"
            img_path = wf_path  * "val/images/" * path;
            # reading the image
            img = convert(Array{Int64}, rawview(channelview(load(img_path))))
            c, h, w = size(img)
            len = max(w, h) # makes image square without changing its size, for evaluation
            len = 64 * ceil(Int, len / 64)
            img = convert(Array{Float32}, img) ./ 255
            img = squaritize_img(img, len)
            img .*= 255
            img .-= avg_img
            reversed_img = reshape(reverse(img, dims=1), (size(img)..., 1))
            imgs_arr = convert(atype, permutedims(reversed_img, (3,2,1,4)))
            c, b, l = predict_image(model, imgs_arr, mode=mode, verbose=false, nms_thold=nms_threshold, conf_thold=conf_threshold, topk=topk)

            invalids = 0
            for person in 1:size(b, 2)
                x1, y1, x2, y2 = floor.(Int, b[:,person])
                if x1 >= w || y1 >= h
                    invalids += 1
                end
            end

            logfile = save_dir * event * "/" * file * ".txt"
            mkpath(save_dir * event);
            open(logfile, "w") do io
                write(io, file * "\n" * string(size(b, 2) - invalids) * "\n")
                for person in 1:size(b, 2)
                    x1, y1, x2, y2 = floor.(Int, b[:,person]); wid = x2 - x1; hei = y2 - y1;
                    if x1 >= w || y1 >= h
                        continue
                    end
                    conf = c[2,person]
                    to_print = string(x1) * " " * string(y1) * " " * string(wid) * " " * string(hei) * " " * string(conf) * "\n"
                    write(io, to_print)
                end
            end
            iter_no += 1; println("Iteration: ", iter_no, " finished! --> Found:", size(b, 2)-invalids);
        end
    end
end

main()

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
struct DSASignature{T} <: AbstractSignature where T<:Integer
    hash::T
    r::T
    s::T
    pubkey::T
end

==(x::DSASignature,y::DSASignature) = x.hash==y.hash && x.r==y.r && x.s==y.s && x.pubkey==y.pubkey

import Base.mod
mod(G::CyclicGroup) = mod(value(G),order(G))
mod(n::Integer,G::CyclicGroup) = mod(n,order(G))

# https://github.com/zydeon/DSA/blob/master/DSA.py
function modinv_euclidean(z::Integer,a::Integer)
    if 0<z<a
	i = a
	j = z
	y1 = BigInt(1)
	y2 = BigInt(0)
	while j > 0
	    q = div(i,j)
	    r = i-j*q
	    y = y2 - y1*q
	    i, j = j, r
	    y2, y1 = y1, y
	    if i == 1
		return mod(y2,a)
            end
        end
    else
        error("Inverse Error")
    end
end

modinv_fermat(k::Integer,a::Integer) = powermod(k,a-2,a)

modinv(k::Integer,G::CyclicGroup) = modinv_euclidean(k,order(G))

function DSASignature(hash,signer::AbstractSigner)

    h = BigInt(hash)
    x = BigInt(signer.privkey)
    G = signer.G

    k = rand(1:order(G)) ### chooses a number from 0 to q
    r = mod(G^k)
    kinv = modinv(k,G)
    s = mod(kinv*(h + x*r),G)

    # This always works for modinv_fermat
    # s < order(G)/2 || s = order(G) - s

    if s==0
        return DSASignature(hash,signer,G)
    else
        return DSASignature(hash,r,s,signer.pubkey)
    end
end

function verify(sr::DSASignature,G::CyclicGroup)
    r,s = sr.r,sr.s
    h = BigInt(sr.hash)
    Y = typeof(G)(sr.pubkey,G)

    0<r<order(G) || return false
    0<s<order(G) || return false

    sinv = modinv(s,G)
    w = mod(sinv,G)
    u1 = mod(h*w,G)
    u2 = mod(r*w,G)
    v = mod(G^u1 * Y^u2)

    if v==r
        return true
    else
        return false
    end
end

import Base.Dict
function Dict(signature::DSASignature)
    config = Dict()
    config["hash"] = string(signature.hash,base=16)
    config["pub"] = string(signature.pubkey,base=16)
    config["r"] = string(signature.r,base=16)
    config["s"] = string(signature.s,base=16)
    return config
end

function DSASignature{BigInt}(dict::Dict)
    hash = parse(BigInt,dict["hash"],base=16)
    r = parse(BigInt,dict["r"],base=16)
    s = parse(BigInt,dict["s"],base=16)
    pub = parse(BigInt,dict["pub"],base=16)
    DSASignature(hash,r,s,pub)
end


text:
__precompile__()
module LanczosInterpolation_mod

#using Yeppp

export get_interp_lag,lanczos_interp_weight!, lanczos_interp_trace!

const NLAG=5

function sinc(x::Float64)::Float64
if abs(x) > 0
   return sin(pi*x)/(pi*x)
else
   return 1.0
end
end # sinc

function dsinc(x::Float64)::Float64
if abs(x) > 0
   return (cos(pi*x)-sinc(pi*x))/x
else
   return 0.0
end
end # dsinc

function lanczos(x::Float64)::Float64
   a::Float64=Float64(NLAG)
   if x < a
      return sinc(x)*sinc(x/a)
   else
      return 0.0
   end
end # lanczos

function get_interp_lag()
lag::Int64 = NLAG
return lag
end

function lanczos_interp_weight!( t::Float64,
                                 n::Int64,
                                 t0::Float64,
                                 dt::Float64,
                                 it::Int64,
                                 wlanczos::Array{Float64,1} )
        if t0 <= t &&  t <= t0+Float64(n-1)*dt
           x::Float64 = (t-t0)/dt
           it = floor(Int64,x)
           dx::Float64 = x-Float64(it)
           it += 1
           for lag=-NLAG:NLAG
                   wlanczos[lag+NLAG+1]=lanczos(dx+Float64(lag))
           end
        else
           it = 0
           wlanczos[1:2*NLAG+1] .= Float64(0.0)
        end
nothing
end # lanczos_interp_weight!

function lanczos_interp_trace!( t::Float64,
                                n::Int64,
                                t0::Float64,
                                dt::Float64,
                                it::Int64,
                                trace::Array{Float64,1} )
        interp::Float64=0.0
        if t0 <= t && t <= t0+Float64(n-1)*dt
           x::Float64 = (t-t0)/dt
           it = floor(Int64,x)
           dx::Float64 = x-Float64(it)
           it += 1
           for lag=-NLAG:NLAG
               is::Int64=Int64(min(max(it+lag,1),n))
               interp += trace[is]*lanczos(dx+Float64(lag))
           end
           return interp
        end
        return Float64(0.0)
end # lanczos_interp_trace!

end


text:
# Solve NY Times spelling bee puzzles

function load_all_words(file_name::String = "/usr/share/dict/words")::Set{String}
    F = open(file_name)
    wlist = readlines(F)

    S = Set{String}()

    for w in wlist
        if length(w) < 4 || any(isuppercase(c) for c in w)
            continue
        end
        push!(S, uppercase(w))
    end

    return S
end

function bee_filter(word::String, middle::Char, surround::String)::Bool
    if middle ∉ word
        return false
    end
    valids = surround * middle
    return all(c ∈ valids for c ∈ word)
end

"""
    bee_solver('a', "bcdefg")
Find all words with 4 or more letters that must use `a`
and any of the letters in "bcdefg".
"""
function bee_solver(middle::Char, surround::String)::Set{String}
    wlist = load_all_words()
    middle = uppercase(middle)
    surround = uppercase(surround)

    Set(w for w ∈ wlist if bee_filter(w, middle, surround))
end

function bee_solver(middle::String, surround::String)::Set{String}
    if length(middle) != 1
        error("First argument should be a single character")
    end
    bee_solver(middle[1], surround)
end

function bingo_filter(word::String, letters::String)::Bool
    all(c ∈ letters for c in word)
end

"""
    bingo(letters)
Find all words that include all the characters in `letters`.
"""
function bingo(letters::String)
    wlist = load_all_words()
    letters = uppercase(letters)
    S = Set(
        word for
        word in wlist if bingo_filter(word, letters) && bingo_filter(letters, word)
    )
    return S
end


text:
module MyTools
using OffsetArrays
using FFTW
include("HHG.jl")
include("fdcoeff.jl")
include("dictandstruct.jl")
include("plotrecipe.jl")
end


text:
const ParameterName=Union{Symbol,Expr}

"""
    Grid(goal=nothing, resolution=10, rng=Random.GLOBAL_RNG, shuffle=true)

Instantiate a Cartesian grid-based hyperparameter tuning strategy with
a specified number of grid points as `goal`, or using a specified
default `resolution` in each numeric dimension.

### Supported ranges:

A single one-dimensional range or vector of one-dimensioinal ranges
can be specified. Specifically, in `Grid` search, the `range` field
of a `TunedModel` instance can be:

- A single one-dimensional range - ie, `ParamRange` object - `r`, or
  pair of the form `(r, res)` where `res` specifies a resolution to
  override the default `resolution`.

- Any vector of objects of the above form

Two elements of a `range` vector may share the same `field` attribute,
with the effect that their grids are combined, as in Example 3 below.

`ParamRange` objects are constructed using the `range` method.

Example 1:

    range(model, :hyper1, lower=1, origin=2, unit=1)

Example 2:

    [(range(model, :hyper1, lower=1, upper=10), 15),
      range(model, :hyper2, lower=2, upper=4),
      range(model, :hyper3, values=[:ball, :tree])]

Example 3:

    # a range generating the grid `[1, 2, 10, 20, 30]` for `:hyper1`:
    [range(model, :hyper1, values=[1, 2]),
     (range(model, :hyper1, lower= 10, upper=30), 3)]

Note: All the `field` values of the `ParamRange` objects (`:hyper1`,
`:hyper2`, `:hyper3` in the preceding example) must refer to field
names a of single model (the `model` specified during `TunedModel`
construction).

### Algorithm

This is a standard grid search with the following specifics: In all
cases all `values` of each specified `NominalRange` are exhausted. If
`goal` is specified, then all resolutions are ignored, and a global
resolution is applied to the `NumericRange` objects that maximizes the
number of grid points, subject to the restriction that this not exceed
`goal`. (This assumes no field appears twice in the `range` vector.)
Otherwise the default `resolution` and any parameter-specific
resolutions apply.

In all cases the models generated are shuffled using `rng`, unless
`shuffle=false`.

See also [`TunedModel`](@ref), [`range`](@ref).

"""
mutable struct Grid <: TuningStrategy
    goal::Union{Nothing,Int}
    resolution::Int
    shuffle::Bool
    rng::Random.AbstractRNG
end

# Constructor with keywords
Grid(; goal=nothing, resolution=10, shuffle=true,
     rng=Random.GLOBAL_RNG) =
    Grid(goal, resolution, MLJBase.shuffle_and_rng(shuffle, rng)...)

isnumeric(::Any) = false
isnumeric(::NumericRange) = true

# To replace resolutions for numeric ranges with goal-adjusted ones if
# a goal is specified:
adjusted_resolutions(::Nothing,  ranges, resolutions) = resolutions
function adjusted_resolutions(goal, ranges, resolutions)
    # get the product Π of the lengths of the NominalRanges:
    len(::NumericRange) = 1
    len(r::NominalRange) = length(r.values)
    Π = prod(len.(ranges))

    n_numeric = sum(isnumeric.(ranges))

    # compute new uniform resolution:
    goal = goal/Π
    res = round(Int, goal^(1/n_numeric))
    return  map(eachindex(resolutions)) do j
        isnumeric(ranges[j]) ? res : resolutions[j]
    end
end

# For deciding scale for duplicated fields:
_merge(s1, s2) = (s1 == :none ? s2 : s1)

function fields_iterators_and_scales(ranges, resolutions)

    # following could have non-unique entries:
    fields = map(r -> r.field, ranges)

    iterator_given_field = Dict{Union{Symbol,Expr},Vector}()
    scale_given_field = Dict{Union{Symbol,Expr},Any}()
    for i in eachindex(ranges)
        fld = fields[i]
        r = ranges[i]
        if haskey(iterator_given_field, fld)
            iterator_given_field[fld] =
                vcat(iterator_given_field[fld], iterator(r, resolutions[i]))
            scale_given_field[fld] =
                _merge(scale_given_field[fld], scale(r))
        else
            iterator_given_field[fld] = iterator(r, resolutions[i])
            scale_given_field[fld] = scale(r)
        end
    end
    fields = unique(fields)
    iterators = map(fld->iterator_given_field[fld], fields)
    scales = map(fld->scale_given_field[fld], fields)

    return fields, iterators, scales

end

function setup(tuning::Grid, model, user_range, verbosity)
    ranges, resolutions =
        process_grid_range(user_range, tuning.resolution, verbosity)

    resolutions = adjusted_resolutions(tuning.goal, ranges, resolutions)

    fields, iterators, parameter_scales =
        fields_iterators_and_scales(ranges, resolutions)

    if tuning.shuffle
        models = grid(tuning.rng, model, fields, iterators)
    else
        models = grid(model, fields, iterators)
    end

    state = (models=models,
             fields=fields,
             parameter_scales=parameter_scales)

    return state

end

MLJTuning.models(tuning::Grid,
                 model,
                 history,
                 state,
                 n_remaining,
                 verbosity) =
                     state.models[_length(history) + 1:end], state

tuning_report(tuning::Grid, history, state) =
    (plotting = plotting_report(state.fields, state.parameter_scales, history),)

function default_n(tuning::Grid, user_range)

    ranges, resolutions =
        process_grid_range(user_range, tuning.resolution, -1)

    resolutions = adjusted_resolutions(tuning.goal, ranges, resolutions)

    fields, iterators, parameter_scales =
        fields_iterators_and_scales(ranges, resolutions)

    return prod(length.(iterators))

end


text:
using JLD
using HTTP
using ..Config: Configuration
using ..Serde: serialize_payload, deserialize_payload, unpack
using ..Fed: STATS, update_stats!, initialize_stats

"""
    CentralNode(host, port, weights, strategy, evaluate, config)

Orchestrate the entire federated learning process. It handles the clients
registration and requests them to train on their local data.
"""
struct CentralNode
    # networking
    host::String
    port::Int
    client_manager::ClientManager

    # machine learning
    weights::Vector{Float32}
    strategy::Function

    # proxy
    evaluate::Function

    config::Configuration

    CentralNode(
        host::String,
        port::Int,
        weights::Vector{Float32},
        strategy::Function,
        evaluate::Function,
        config::Configuration
    ) = new(host, port, ClientManager(), weights, strategy, evaluate, config)
end

function get_store_url(url)
    protocol, url, port = split(url, ":")
    return protocol * ":" * url * ":" * string(parse(Int, port) + 1010)
end

"""
    fit(central_node::CentralNode)

Fit the global model by orchestrating the federated training.
"""
function fit(central_node::CentralNode)
    global_weights = central_node.weights

    initialize_stats(central_node.config, length(global_weights))

    jldopen("data/request.jld", "w") do req_file
        jldopen("data/response.jld", "w") do res_file
            for round_num in 1:central_node.config.base.num_comm_rounds
                @info "Communication round $round_num"

                # uplink communication
                # 1. chose the clients subset for the round.
                clients = sample_clients(central_node.client_manager, central_node.config.base.fraction_clients)
                # 2. serialize the global model.
                payload = serialize_payload(central_node.config.payload_serde, global_weights)
                # 3. ask clients to train on the global model.
                clients_payload = fit_clients(central_node.config.base.fit_node, clients, payload)

                write(req_file, "$round_num", global_weights)
                client_w = map(p -> p.data, unpack.(clients_payload))
                write(res_file, "$round_num", client_w)

                # downlink communication
                # 4. deserialize the results.
                round_weights = [
                    deserialize_payload(central_node.config.payload_serde, payload, get_store_url(clients[i]))
                    for (i, payload) in enumerate(clients_payload)
                ]

                # update, convergeance tracking, ...
                # 5. update the global model.
                global_weights = central_node.strategy(round_weights)
                # 6. evaluate global model.
                loss, acc = central_node.evaluate(global_weights)
                @info "loss: $loss, acc: $acc"

                # 7. record statistics.
                update_stats!(STATS, round_num, loss, acc)

                save("stats.jld", "stats", STATS)
            end
        end
    end
end

"""
    fit_clients(endpoint::String, clients::Vector{String}, payload::Vector{UInt8})

Ask asynchronously each client from `clients` to train on their local data by
using the model weights contained in the `payload`. Return a `Vector` where each
element is the serialized weights of one client.
"""
function fit_clients(
    endpoint::String,
    clients::Vector{String},
    payload::Vector{UInt8}
)::Vector{Vector{UInt8}}
    # asynchronously ask the clients subset to train, and wait on the results.
    tasks = [@async fit_client(client * endpoint, payload) for client in clients]
    wait.(tasks)

    # extract the results.
    round_weights = map(task -> task.result, tasks)

    return round_weights
end

"""
    fit_client(url::String, payload::Vector{UInt8})

Ask one `client` to train on its local data, by using the model weights
contained in the `payload`. Return the serialized updated `weights`.
"""
function fit_client(
    url::String,
    payload::Vector{UInt8}
)::Vector{UInt8}
    response = HTTP.request("POST", url, [], payload)
    return response.body
end