text:
using Random
export customop,
xavier_init,
load_op_and_grad,
load_op,
use_gpu,
test_jacobian,
install,
load_system_op,
install_adept,
register,
debug

"""
    xavier_init(size, dtype=Float64)

Returns a matrix of size `size` and its values are from Xavier initialization.
"""
function xavier_init(size, dtype=Float64)
    in_dim = size[1]
    xavier_stddev = 1. / sqrt(in_dim / 2.)
    return randn(dtype, size...)*xavier_stddev
end

############### custom operators ##################
function cmake(DIR::String="..")
    ENV_ = copy(ENV)
    if haskey(ENV_, "LD_LIBRARY_PATH")
        ENV_["LD_LIBRARY_PATH"] = ENV["LD_LIBRARY_PATH"]*":$LIBDIR"
    else
        ENV_["LD_LIBRARY_PATH"] = LIBDIR
    end
    if Sys.islinux()
        run(setenv(`$CMAKE -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX $DIR`, ENV_))
    else
        run(setenv(`$CMAKE $DIR`, ENV_))
    end
end

function make()
    ENV_ = copy(ENV)
    if haskey(ENV_, "LD_LIBRARY_PATH")
        ENV_["LD_LIBRARY_PATH"] = ENV["LD_LIBRARY_PATH"]*":$LIBDIR"
    else
        ENV_["LD_LIBRARY_PATH"] = LIBDIR
    end
    run(setenv(`$MAKE -j`, ENV_))
end

function clean()
    run(`$MAKE clean`)
end

load_op_dict = Dict{Tuple{String, String}, PyObject}()
load_op_grad_dict = Dict{Tuple{String, String}, PyObject}()

@doc """
    load_op(oplibpath::String, opname::String)

Loads the operator `opname` from library `oplibpath`.
"""
function load_op(oplibpath::String, opname::String)
    if splitext(oplibpath)[2]==""
        oplibpath = abspath(oplibpath * (Sys.islinux() ?
                        ".so" : Sys.isapple() ? ".dylib" : ".dll"))
    end
    oplibpath = abspath(oplibpath)
    if haskey(load_op_dict, (oplibpath,opname))
        return load_op_dict[(oplibpath,opname)]
    end

    if !isfile(oplibpath)
        error("File $oplibpath does not exist. Instruction:\nRunning `compile(oplibpath)` to compile the library first.")
    end
    fn_name = opname*randstring(8)
py"""
import tensorflow as tf
lib$$fn_name = tf.load_op_library($oplibpath)
"""
    lib = py"lib$$fn_name"
    s = getproperty(lib, opname)
    load_op_dict[(oplibpath,opname)] = s
    printstyled("Load library operator: $oplibpath ==> $opname\n", color=:green)
    return s
end

@doc """
    load_op_and_grad(oplibpath::String, opname::String; multiple::Bool=false)

Loads the operator `opname` from library `oplibpath`; gradients are also imported.
If `multiple` is true, the operator is assumed to have multiple outputs.
"""
function load_op_and_grad(oplibpath::String, opname::String; multiple::Bool=false)
    if splitext(oplibpath)[2]==""
        oplibpath = oplibpath * (Sys.islinux() ?
                        ".so" : Sys.isapple() ? ".dylib" : ".dll")
    end
    oplibpath = abspath(oplibpath)
    if haskey(load_op_grad_dict, (oplibpath,opname))
        return load_op_grad_dict[(oplibpath,opname)]
    end
    if !isfile(oplibpath)
        error("File $oplibpath does not exist. Instruction:\nRunning `compile(oplibpath)` to compile the library first.")
    end

    opname_grad = opname*"_grad"
    fn_name = opname*randstring(8)
if !multiple
py"""
import tensorflow as tf
lib$$fn_name = tf.load_op_library($oplibpath)
@tf.custom_gradient
def $$fn_name(*args):
    u = lib$$fn_name.$$opname(*args)
    def grad(dy):
        return lib$$fn_name.$$opname_grad(dy, u, *args)
    return u, grad
"""
else
py"""
import tensorflow as tf
lib$$fn_name = tf.load_op_library($oplibpath)
@tf.custom_gradient
def $$fn_name(*args):
    u = lib$$fn_name.$$opname(*args)
    def grad(*dy):
        dy = [y for y in dy if y is not None and y.dtype in [tf.float64, tf.float32]] # only float64 and float32 can backpropagate gradients
        return lib$$fn_name.$$opname_grad(*dy, *u, *args)
    return u, grad
"""
end
        s = py"$$fn_name"
        load_op_grad_dict[(oplibpath,opname)] = s
        printstyled("Load library operator (with gradient): $oplibpath ==> $opname\n", color=:green)
        return s
end

function refresh_cmake()
    NEWLIB = String[]
    for c in COLIB
        push!(NEWLIB, c.second[1])
    end
    NEWLIB = unique(NEWLIB)
    cmakecnt = read(joinpath(@__DIR__, "../deps/CustomOps/CMakeListsTemplate.txt"),String)
    cmakecnt = replace(cmakecnt, "set(LIBDIR_NAME"=>"set(LIBDIR_NAME "*join(NEWLIB, ' '))
    open(joinpath(@__DIR__, "../deps/CustomOps/CMakeLists.txt"), "w") do io
        write(io, cmakecnt)
    end
end

"""
    load_system_op(s::String, oplib::String, grad::Bool=true)

Loads custom operator from CustomOps directory (shipped with ADCME instead of TensorFlow)
For example
```
s = "SparseOperator"
oplib = "libSO"
grad = true
```
this will direct Julia to find library `CustomOps/SparseOperator/libSO.dylib` on MACOSX
"""
function load_system_op(s::String, oplib::String, opname::String, grad::Bool=true;
    return_str::Bool=false, multiple::Bool=false)
    dir = joinpath(joinpath("$(@__DIR__)", "../deps/CustomOps"), s)
    if !isdir(dir)
        error("Folder for the operator $s does not exist: $dir")
    end
    oplibpath = joinpath(joinpath(dir, "build"), oplib)
    # check if the library exists
    libfile = oplibpath * (Sys.islinux() ?
                        ".so" : Sys.isapple() ? ".dylib" : ".dll")
    # @show libfile
    if !isfile(libfile)
        @info "Lib $s exists in registery but was not initialized. Compiling..."
        ADCME.precompile()
    end
    if return_str
        return oplibpath
    end
    if grad
        load_op_and_grad(oplibpath, opname; multiple=multiple)
    else
        load_op(oplibpath, opname)
    end
end

load_system_op(s::String; kwargs...) = load_system_op(COLIB[s]...; kwargs...)

"""
    compile(s::String; force::Bool=false)

Compiles the library given by path `deps/s`. If `force` is false, `compile` first check whether
the binary product exists. If the binary product exists, return 2. Otherwise, `compile` tries to
compile the binary product, and returns 0 if successful; it return 1 otherwise.
"""
function compile(s::String; force::Bool=false, customdir::Bool = false)
    PWD = pwd()
    dir = s
    if !customdir
        dir = joinpath(joinpath("$(@__DIR__)", "../deps/CustomOps"), s)
    end
    if !isdir(dir)
        @warn("Folder for the operator $s does not exist: $dir")
        return 1
    end
    cd(dir)

    local surfix
    if Sys.isapple()
        surfix = ".dylib"
    elseif Sys.islinux()
        surfix = ".so"
    elseif Sys.iswindows()
        surfix = ".dll"
    end
    if !force && isdir("build") # check if product exists
        files = readdir("build")
        if any([endswith(x, surfix) for x in files])
            @warn("The binary product exists.")
            cd(PWD)
            return 2
        end
    end
    rm("build",force=true,recursive=true)
    mkdir("build")
    cd("build")
    try
        cmake()
        make()
        cd(PWD)
        return 0
    catch e
        error("Compilation error: $e")
        cd(PWD)
        return 1
    end
end

"""
    precompile(force::Bool=true)

Compiles all the operators in `formulas.txt`.
"""
function Base.:precompile(force::Bool=true)
    try
        run(`which julia`)
    catch
        printstyled(
"""Julia cannot be found using `which julia`. This will break custom operator.
To fix the error, add the julia binary path to your PATH environment variable:

export PATH=$(Sys.BINDIR):\$PATH

""",
        color=:red)
        error("Waiting for fixing the Julia binary path.")
    end

    if !(ADCME.LIBDIR in split(ENV["LD_LIBRARY_PATH"], ':'))
        @warn "$(ADCME.LIBDIR) is not in LD_LIBRARY_PATH; this may break the custom operator utilties.
You could add the path to LD_LIBRARY_PATH:

export LD_LIBRARY_PATH = $(ADCME.LIBDIR):\$LD_LIBRARY_PATH

"
    end
    PWD = pwd()
    if (!force) && isfile("$(@__DIR__)/../deps/CustomOps/CMakeLists.txt") &&
            isdir("$(@__DIR__)/../deps/CustomOps/build")
        @info "Reuse existing cmake files"
        PWD = pwd()
        cd("$(@__DIR__)/../deps/CustomOps/build")
        ADCME.cmake()
        ADCME.make()
        cd(PWD)
        return
    end
    refresh_cmake()
    cd("$(@__DIR__)/../deps/CustomOps")
    !isdir("build") && mkdir("build")
    cd("build")
    ADCME.cmake()
    ADCME.make()
    cd(PWD)
end

"""
    customop(simple::Bool=false)

Create a new custom operator. If `simple=true`, the custom operator only supports CPU and does not have gradients.

# Example

```julia-repl
julia> customop() # create an editable `customop.txt` file
[ Info: Edit custom_op.txt for custom operators
julia> customop() # after editing `customop.txt`, call it again to generate interface files.
```
"""
function customop(simple::Bool=false)
    # install_custom_op_dependency()
    py_dir = "$(@__DIR__)/../examples/custom_op/template"
    if !("custom_op.txt" in readdir("."))
        cp("$(py_dir)/custom_op.example", "custom_op.txt")
        @info "Edit custom_op.txt for custom operators"
        return
    else
        python = PyCall.python
        run(`$python $(py_dir)/customop.py custom_op.txt $py_dir $simple`)
    end
end

function use_gpu(i::Union{Nothing,Int64}=nothing)
    dl = pyimport("tensorflow.python.client.device_lib")
    if !isnothing(i) && i>=1
        i = join(collect(0:i-1),',')
        ENV["CUDA_VISIBLE_DEVICES"] = i
    elseif !isnothing(i) && i==0
        ENV["CUDA_VISIBLE_DEVICES"] = ""
    end
    local_device_protos = dl.list_local_devices()
    return [x.name for x in local_device_protos if x.device_type == "GPU"]
end

"""
    test_jacobian(f::Function, x0::Array{Float64}; scale::Float64 = 1.0)

Testing the gradients of a vector function `f`:
`y, J = f(x)` where `y` is a vector output and `J` is the Jacobian.
"""
function test_jacobian(f::Function, x0::Array{Float64}; scale::Float64 = 1.0)
    v0 = rand(Float64,size(x0))
    γs = scale ./10 .^(1:5)
    err2 = []
    err1 = []
    f0, J = f(x0)
    for i = 1:5
        f1, _ = f(x0+γs[i]*v0)
        push!(err1, norm(f1-f0))
        @show f1, f0, 2γs[i]*J*v0
        push!(err2, norm(f1-f0-γs[i]*J*v0))
        # push!(err2, norm((f1-f2)/(2γs[i])-J*v0))
        # #@show "test ", f1, f2, f1-f2
    end
    close("all")
    loglog(γs, err2, label="Automatic Differentiation")
    loglog(γs, err1, label="Finite Difference")
    loglog(γs, γs.^2 * 0.5*abs(err2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(γs, γs * 0.5*abs(err1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
    plt.gca().invert_xaxis()
    legend()
    println("Finite difference: $err1")
    println("Automatic differentiation: $err2")
    return err1, err2
end

"""
    install(s::String; force::Bool = false)

Install a custom operator via URL. `s` can be
- A URL. ADCME will download the directory through `git`
- A string. ADCME will search for the associated package on https://github.com/ADCMEMarket
"""
function install(s::String; force::Bool = false)
    global COLIB
    codir = "$(@__DIR__)/../deps/CustomOps"
    if !startswith(s, "https://github.com")
        s = "https://github.com/ADCMEMarket/"*s
    end
    _, name = splitdir(s)
    if name in readdir(codir)
        if force
            rm(joinpath(codir, name), recursive=true, force=true)
        else
            error("$name already in $codir, fix it with\n\n\tinstall(\"$s\", force=true)\n")
        end
    end
    try
        LibGit2.clone(s, joinpath(codir, name))
    catch
        LibGit2.clone("git://$(s[9:end]).git", joinpath(codir, name))
    end

    # If there is a build script `build.sh`, run the build script
    if isfile(joinpath(joinpath(codir, name), "build.sh"))
        PWD = pwd()
        cd(joinpath(codir, name))
        try
            run(`bash build.sh`)
        catch
        end
        cd(PWD)
    end

    formula = eval(Meta.parse(read(joinpath(joinpath(codir, name),"formula.txt"), String)))
    if isnothing(formula)
        error("Broken package: $s does not have formula.txt.")
    else
        @info "Add formula $formula"
        push!(COLIB, formula)
    end

    rm("$(@__DIR__)/../deps/CustomOps/formulas.jl", force=true)
    open("$(@__DIR__)/../deps/CustomOps/formulas.jl", "a") do io
        n = length(COLIB)
        for (k, c) in enumerate(COLIB)
            if k == n
                write(io, string(c)*"\n")
            else
                write(io, string(c)*",\n")
            end
        end
    end
end

"""
    install_adept(force::Bool=false)

Install adept-2 library: https://github.com/rjhogan/Adept-2
"""
function install_adept(force::Bool=false)
    PWD = pwd()
    cd(ADCME.LIBDIR)
    if !isdir("Adept-2")
        LibGit2.clone("https://github.com/ADCMEMarket/Adept-2", "Adept-2")
    end
    cd("Adept-2/adept")
    if !("openblas" in Conda._installed_packages())
        Conda.add("openblas", channel="anaconda")

        required_file = Sys.isapple() ? ".dylib" : (Sys.islinux() ? ".so" : ".dll")
        required_file = joinpath(ADCME.LIBDIR, "libopenblas")*required_file
        if !isfile(required_file)
            files = readdir(ADCME.LIBDIR)
            files = filter(x->!isnothing(x), match.(r"(libopenblas\S*.dylib)", files))[1]
            target = joinpath(ADCME.LIBDIR, files[1])
            symlink(target, required_file)
            @info "Symlink $(required_file) --> $(files[1])"
        end
    end
    try
        if force==true
            @info "remove .libs ..."
            rm(".libs", force=true, recursive=true)
        end
        if !isdir(".libs")
            @info """Copy "$(@__DIR__)/../deps/AdeptCMakeLists.txt" to "$(joinpath(pwd(), "CMakeLists.txt"))" ... """
            cp("$(@__DIR__)/../deps/AdeptCMakeLists.txt", "./CMakeLists.txt", force=true)
            @info """Remove $(joinpath(pwd(), "build")) ... """
            rm("build", force=true, recursive=true)
            @info "Make $(joinpath(pwd(), "build")) ... "
            mkdir("build")
            @info "Change directory into $(joinpath(pwd(), "build")) ... "
            cd("build")
            @info "Cmake ... "
            ADCME.cmake()
            @info "Make ... "
            ADCME.make()
        end
        printstyled("""
∘ Add the following lines to CMakeLists.txt

include_directories(\${LIBDIR}/Adept-2/include)
link_directories(\${LIBDIR}/Adept-2/adept/.libs)

∘ Add `adept` to `target_link_libraries`

∘ Add `$LIBDIR` to `LD_LIBRARY_PATH` environment variable
""", color=:green)
    catch
        printstyled("Compliation failed\n", color=:red)
    finally
        cd(PWD)
    end
end

# function register(forward::Function, backward::Function, fd_args::Int64, bd_args::Int64; multiple::Bool=false)
#     fn_name = "customgrad_"*randstring(8)

# fd = join(["o$k" for k = 1:fd_args], ",")
# bd = join(["o$k" for k = 1:bd_args], ",")

# py"""
# forward_$$fn_name = lambda $$fd: $forward($$fd)
# backward_$$fn_name = lambda $$bd: $backward($$bd)
# """
#     if !multiple
# py"""
# import tensorflow as tf
# @tf.custom_gradient
# def $$fn_name(*args):
#     u = forward_$$fn_name(*args)
#     def grad(dy):
#         return backward_$$fn_name(dy, u, *args)
#     return u, grad
# """
#     else
# py"""
# import tensorflow as tf
# @tf.custom_gradient
# def $$fn_name(*args):
#     u = forward_$$fn_name(*args)
#     def grad(*dy):
#         dy = [y for y in dy if y is not None and y.dtype in [tf.float64, tf.float32]] # only float64 and float32 can backpropagate gradients
#         return backward_$$fn_name(*dy, *u, *args)
#     return u, grad
# """
#     end
#     return py"$$fn_name"
# end

@doc raw"""
    register(forward::Function, backward::Function; multiple::Bool=false)

Register a function `forward` with back-propagated gradients rule `backward` to the backward.
∘ `forward`: it takes $n$ inputs and outputs $m$ tensors. When $m>1$, the keyword `multiple` must be true.
∘ `backward`: it takes $\tilde m$ top gradients from float/double output tensors of `forward`, $m$ outputs of the `forward`,
   and $n$ inputs of the `forward`. `backward` outputs $n$ gradients for each input of `forward`. When input $i$ of
   `forward` is not float/double, `backward` should return `nothing` for the corresponding gradients.

# Example
```julia
forward = x->log(1+exp(x))
backward = (dy, y, x)->dy*(1-1/(1+y))
f = register(forward, backward)
```
"""
function register(forward::Function, backward::Function; multiple::Bool=false)
    fn_name = "customgrad_"*randstring(8)
    if !multiple
py"""
import tensorflow as tf
@tf.custom_gradient
def $$fn_name(*args):
    u = $forward(*args)
    def grad(dy):
        return $backward(dy, u, *args)
    return u, grad
"""
    else
py"""
import tensorflow as tf
@tf.custom_gradient
def $$fn_name(*args):
    u = forward_$$fn_name(*args)
    def grad(*dy):
        dy = [y for y in dy if y is not None and y.dtype in [tf.float64, tf.float32]] # only float64 and float32 can backpropagate gradients
        return backward_$$fn_name(*dy, *u, *args)
    return u, grad
"""
    end
    return py"$$fn_name"
end

"""
    debug(sess::PyObject, o::PyObject)

In the case a session run yields `InvalidArgumentError()`, this function can help print the exact error.
"""
function debug(sess::PyObject, o::PyObject)
py"""
import tensorflow as tf
import traceback
try:
    $sess.run($o)
except Exception:
    print(traceback.format_exc())
"""
end


text:
# Potential contributions to MOI.

# # --- operate_coefficient with three arguments: coefficient,
# # first variable or nothing, second variable or nothing.

# function operate_index_coefficient(f, term::MOI.ScalarAffineTerm)
#     return MOI.ScalarAffineTerm(f(term.coefficient, term.variable, nothing), term.variable)
# end
# function operate_index_coefficient(f, term::MOI.ScalarQuadraticTerm)
#     return MOI.ScalarQuadraticTerm(
#         f(term.coefficient, term.variable_1, term.variable_2),
#         term.variable_1,
#         term.variable_2,
#     )
# end
# function operate_index_coefficient(f, term::MOI.VectorAffineTerm)
#     return MOI.VectorAffineTerm(
#         term.output_index,
#         operate_index_coefficient(f, term.scalar_term),
#     )
# end
# function operate_index_coefficient(f, term::MOI.VectorQuadraticTerm)
#     return MOI.VectorQuadraticTerm(
#         term.output_index,
#         operate_index_coefficient(f, term.scalar_term),
#     )
# end

# function operate_index_coefficients(f, func::MOI.ScalarAffineFunction)
#     return MOI.ScalarAffineFunction(
#         [operate_coefficient(f, term) for term in func.terms],
#         f(func.constant, nothing, nothing),
#     )
# end
# function operate_index_coefficients(f, func::MOI.ScalarQuadraticFunction)
#     return MOI.ScalarQuadraticFunction(
#         [operate_index_coefficient(f, term) for term in func.affine_terms],
#         [operate_index_coefficient(f, term) for term in func.quadratic_terms],
#         f(func.constant, nothing, nothing),
#     )
# end
# function operate_index_coefficients(f, func::MOI.VectorAffineFunction)
#     return MOI.VectorAffineFunction(
#         [operate_index_coefficient(f, term) for term in func.terms],
#         map(f, func.constants),
#     )
# end
# function operate_index_coefficients(f, func::MOI.VectorQuadraticFunction)
#     return MOI.VectorQuadraticFunction(
#         [operate_index_coefficient(f, term) for term in func.affine_terms],
#         [operate_index_coefficient(f, term) for term in func.quadratic_terms],
#         map(f, func.constants),
#     )
# end

# # --- operate_coefficient with two arguments: coefficient,
# # number of the output dimension. Only for vector functions.

# function operate_dimension_coefficient(f, term::MOI.VectorAffineTerm)
#     return MOI.VectorAffineTerm(
#         term.output_index,
#         operate_dimension_coefficient(
#             (v) -> f(v, term.output_index),
#             term.scalar_term,
#         ),
#     )
# end
# function operate_dimension_coefficient(f, term::MOI.VectorQuadraticTerm)
#     return MOI.VectorQuadraticTerm(
#         term.output_index,
#         operate_dimension_coefficient(
#             (v) -> f(v, term.output_index),
#             term.scalar_term,
#         ),
#     )
# end

# function operate_dimension_coefficients(f, func::MOI.VectorAffineFunction)
#     return MOI.VectorAffineFunction(
#         [operate_dimension_coefficient(f, term) for term in func.terms],
#         map(f, func.constants),
#     )
# end
# function operate_dimension_coefficients(f, func::MOI.VectorQuadraticFunction)
#     return MOI.VectorQuadraticFunction(
#         [operate_dimension_coefficient(f, term) for term in func.affine_terms],
#         [operate_dimension_coefficient(f, term) for term in func.quadratic_terms],
#         map(f, func.constants),
#     )
# end

# # --- from a vector function to a scalar function by summing all components.

# function Base.sum(func::MOI.VectorAffineFunction)
#     return MOIU.canonicalize(
#         MOI.ScalarAffineFunction(
#             [term.scalar_term for term in func.terms],
#             sum(func.constants)
#         )
#     )
# end
# function Base.sum(func::MOI.VectorQuadraticFunction)
#     return MOIU.canonicalize(
#         MOI.ScalarAffineFunction(
#             [term.scalar_term for term in func.affine_terms],
#             [term.scalar_term for term in func.quadratic_terms],
#             sum(func.constants)
#         )
#     )
# end

# --- casting from Bool to Int

# function MOI.ScalarAffineTerm{Int}(t::MOI.ScalarAffineTerm{Bool})
#     return MOI.ScalarAffineTerm{Int}(
#         Int(t.coefficient),
#         t.variable,
#     )
# end

# function MOI.VectorAffineTerm{Int}(t::MOI.VectorAffineTerm{Bool})
#     return MOI.VectorAffineTerm{Int}(
#         t.output_index,
#         MOI.ScalarAffineTerm{Int}(t.scalar_term),
#     )
# end

# function MOI.ScalarAffineFunction{Int}(f::MOI.ScalarAffineFunction{Bool})
#     return MOI.ScalarAffineFunction{Int}(
#         MOI.ScalarAffineTerm{Int}.(f.terms),
#         Int(f.constant),
#     )
# end

# function MOI.VectorAffineFunction{Int}(f::MOI.VectorAffineFunction{Bool})
#     return MOI.VectorAffineFunction{Int}(
#         MOI.VectorAffineTerm{Int}.(f.terms),
#         Int.(f.constants),
#     )
# end

ext: jl


text:
# AbstractGridComponent

export AbstractGridComponent
"""
    abstract type AbstractGridComponent{S,T,D} <: AbstractGridMesh{S,T,D} end
"""
abstract type AbstractGridComponent{S,T,D} <: AbstractGridMesh{S,T,D} end

export box
box(::AbstractGridComponent) = error("undefined")

export test_AbstractGridComponent
function test_AbstractGridComponent(gf::AbstractGridComponent)
    test_AbstractGridMesh(gf)

    @testset "$(typeof(gf)) <: AbstractGridComponent" begin
        D = ndims(gf)
        gfbox = box(gf)
        @test gfbox isa BBox{D}

        T = eltype(gf)
        for n in 1:10
            ind = map((a, b) -> rand(a:b), min(gfbox), max(gfbox))
            val = gf[CartesianIndex(Tuple(ind))]
            @test val isa T
        end
    end
end


text:
class Node {
  Shape elem;
  Node next;

  void setElem(Shape c) { elem = c; }

  void setNext(Node n) { next = n; }

  Shape getElem() { return elem; }

  Node getNext() { return next; }
}

class Stack {
  Node head;

  void push(Shape c) {
    Node newHead = new Node;
    newHead.setElem(c);
    newHead.setNext(head);
    head = newHead;
  }

  boolean isEmpty() {
    return head==(Node)null;
  }

  Shape top() {
    return head.getElem();
  }

  void pop() {
    head = head.getNext();
  }
}

class Shape {
  void tell () {
    printString("I'm a shape");
  }

  void tellAgain() {
     printString("I'm just a shape");
  }
}

class Rectangle extends Shape {
  void tellAgain() {
    printString("I'm really a rectangle");
  }
}

class Circle extends Shape {
  void tellAgain() {
    printString("I'm really a circle");
  }
}

class Square extends Rectangle {
  void tellAgain() {
    printString("I'm really a square");
  }
}

int main() {
  Stack stk = new Stack;
  Shape s = new Shape;
  stk.push(s);
  s = new Rectangle;
  stk.push(s);
  s = new Square;
  stk.push(s);
  s = new Circle;
  stk.push(s);
  while (!stk.isEmpty()) {
    s = stk.top();
    s.tell();
    s.tellAgain();
    stk.pop();
  }
  return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
"""
   compute the source-side wave field as (p[it+1/2]-p[it-1/2]-f[it]) / k
"""
function apply_image_condition!(g::Vector{Tv}, spt::Snapshot, wfd1::Wavefield,
         wfd2::Wavefield, params::TdParams) where {Tv <: AbstractFloat}

    # number of elements
    N = params.nz * params.nx

    # update gradient by current adjoint snapshot
    for i = 1 : N
        dpdt = 2.0 * (wfd2.p[i] - wfd1.p[i]) / params.vel[i]
        j    = params.spt2wfd[i]
        g[i] = g[i] + dpdt * spt.px[j]
    end

    return nothing
end

"""
   reconstruct source-side backward, which can be used for the adjoint of born
approximation
"""
function sourceside_reconstruct_backward(path_bnd::Tp, path_wfd::Tp,
         src::Ts, params::TdParams) where {Tp<:String, Ts<:Union{Source,Vector{Source}}}

    # make the source term as a vector has one element
    if Ts <: Source
       srcs    = Vector{Source}(undef,1)
       srcs[1] = src
    else
       srcs    = src
    end
    ns = length(srcs)

    # length of one-step pressure field
    N = params.nz * params.nx

    # initialize intermediate variables
    wfd1 = Wavefield(params)
    wfd2 = read_one_wavefield(path_wfd, 1)
    tmp_z1 = zeros(params.data_format, params.nz)
    tmp_z2 = zeros(params.data_format, params.nz)
    tmp_x1 = zeros(params.data_format, params.nx)
    tmp_x2 = zeros(params.data_format, params.nx)

    # memory for saving the sourceside wavefield at one time step
    dpdt   = zeros(params.data_format, N)

    # initialize boundary value
    fid_bnd = open(path_bnd, "r")
    bnd     = WavefieldBound(params)

    # pre-allocate memory for saving the reconstructed pressure field
    pre = zeros(params.data_format, N * (params.nt-1))
    idx_o = N * (params.nt-2) + 1

    # backward iterations
    for it = params.nt-1 : -1 : 1

        # remove dt * src[it+1]
        subtract_source!(wfd2, srcs, it+1)

        # read the boundary value
        read_one_boundary!(bnd, fid_bnd, it, params)

        # one step backward reconstruction
        one_step_backward!(wfd1, wfd2, bnd, params,
                           tmp_z1, tmp_z2, tmp_x1, tmp_x2)

        # compute the source-side wavefield
        for i = 1 : N
            dpdt[i] = 2.0 * (wfd2.p[i] - wfd1.p[i]) / params.vel[i]
        end

        # save it
        copyto!(pre, idx_o, dpdt, 1, N)
        idx_o = idx_o - N

        # prepare for next step
        copy_wavefield!(wfd2, wfd1)
    end

    close(fid_bnd)
    return reshape(pre, params.nz, params.nx, params.nt-1)
end

"""
   the adjoint operator of born approximation of simultaneouse source
"""
function born_approximation_adjoint(rec::Recordings, path_bnd::Tp, path_wfd::Tp,
                           src::Ts, params::TdParams) where {Tp<:String, Ts<:Union{Source,Vector{Source}}}

    # make the source term as a vector has one element
    if Ts <: Source
       srcs    = Vector{Source}(undef,1)
       srcs[1] = src
    else
       srcs    = src
    end
    ns = length(srcs)

    # model length
    N = params.nz * params.nx

    # allocate memory for computing adjoint wavefield
    spt1 = Snapshot(params)
    spt2 = Snapshot(params)
    tmp    = zeros(params.data_format, params.Nz * params.Nx)
    tmp_z1 = zeros(params.data_format, params.Nz)
    tmp_z2 = zeros(params.data_format, params.Nz)
    tmp_x1 = zeros(params.data_format, params.Nx)
    tmp_x2 = zeros(params.data_format, params.Nx)

    # allocate memory for reconstructing source-side wavefield backward
    wfd1 = Wavefield(params)
    wfd2 = read_one_wavefield(path_wfd, 1)        # read the last wavefield
    wfd_z1 = zeros(params.data_format, params.nz)
    wfd_z2 = zeros(params.data_format, params.nz)
    wfd_x1 = zeros(params.data_format, params.nx)
    wfd_x2 = zeros(params.data_format, params.nx)

    # initialize the boundary value as zero
    bnd     = WavefieldBound(params)
    fid_bnd = open(path_bnd, "r")

    # initialize gradient to be zeros
    m = zeros(params.data_format, N)

    # backward time stepping
    for it = params.nt : -1 : 2

        # compute adjoint wavefield
        one_step_adjoint!(spt1, spt2, params, tmp, tmp_z1, tmp_z2, tmp_x1, tmp_x2)
        inject_rec2spt!(spt1, rec, it)

        # reconstruct source-side wavefield
        subtract_source!(wfd2, srcs, it)
        read_one_boundary!(bnd, fid_bnd, it-1, params)
        one_step_backward!(wfd1, wfd2, bnd, params, wfd_z1, wfd_z2, wfd_x1, wfd_x2)

        # compute the gradient
        apply_image_condition!(m, spt1, wfd1, wfd2, params)

        # prepare for the next step backward reconstruction
        copy_snapshot!(spt2, spt1)
        copy_wavefield!(wfd2, wfd1)

    end

    # close the boundary value file
    close(fid_bnd)

    # return the gradient of velocity model
    return m
end

"""
   the adjoint operator of born approximation
"""
function born_approximation_adjoint(path_m::Ts, dir_rec::Ts, dir_sourceside::Ts, fidiff::TdParams;
                          normalization_flag=true, remove_flag=true, mute_index::Int64=0) where {Ts<:String, T<:Union{Source, Vector{Source}}}

    function wrap_born_adjoint(params::NamedTuple)

        rec = read_recordings(params.path_rec)
        src = read_source(params.path_src)

        m   = born_approximation_adjoint(rec, params.path_bnd, params.path_wfd, src, params.fidiff)
        hdr = RegularSampleHeader(m; title="temporary file")
        write_RSdata(params.path_tmp, hdr, m)

        return nothing
    end

    # set up directory
    dir_bnd = joinpath(dir_sourceside, "boundary")
    dir_wfd = joinpath(dir_sourceside, "wavefield")
    dir_src = joinpath(dir_sourceside, "source")
    isdir(dir_rec) || error("reflections doesn't exist")
    isdir(dir_bnd) || error("boundary file doesn't exist")
    isdir(dir_wfd) || error("wavefield file doesn't exist")
    isdir(dir_src) || error("source wavelet file doesn't exist")

    # the directory to save model
    dir_m = splitdir(path_m)[1]
    if !isdir(dir_m)
       mkdir(dir_m)
       if !isdir(dir_m) # check the directory is created
          error("can't create directory for the variable in model space")
       end
    end

    dir_tmp = joinpath(dir_m, "temporary")
    if !isdir(dir_tmp)
       mkdir(dir_tmp)
       if !isdir(dir_tmp) # check the directory is created
          error("can't create directory for one shot adjoint born approximation")
       end
    end

    # determine the number of shot
    file_rec = readdir(dir_rec)
    ns       = length(file_rec)
    if file_rec[1] == ".DS_Store"
       ns = ns - 1
    end

    file_src = readdir(dir_src)
    nb       = length(file_src)
    if file_src[1] == ".DS_Store"
       nb = nb - 1
    end
    ns == nb || error("number of shot doesn't match")

    argument_collection = Vector{NamedTuple}(undef, ns)
    for i = 1 : ns

        file_name = join(["temporary_" "$i" ".bin"])
        path_tmp  = joinpath(dir_tmp, file_name)

        file_name = join(["recordings_" "$i" ".bin"])
        path_rec  = joinpath(dir_rec, file_name)

        file_name = join(["boundary_" "$i" ".rsf"])
        path_bnd  = joinpath(dir_bnd, file_name)

        file_name = join(["wavefield_" "$i" ".rsf"])
        path_wfd  = joinpath(dir_wfd, file_name)

        file_name = join(["source_" "$i" ".bin"])
        path_src  = joinpath(dir_src, file_name)

        argument_collection[i] = (path_tmp=path_tmp, path_rec=path_rec, path_bnd=path_bnd,
                                  path_wfd=path_wfd, path_src=path_src, fidiff=fidiff)
    end

    # do simulation parallel
    if nprocs() == 1
       for i = 1 : ns
           wrap_born_adjoint(argument_collection[i])
       end
    else
       pmap(wrap_born_adjoint, argument_collection)
    end

    # read the adjoint result
    (hdr, p) = read_RSdata(argument_collection[1].path_tmp)
    for i = 2 : ns
        (hdr, tmp) = read_RSdata(argument_collection[i].path_tmp)
        p .= p .+ tmp
    end

    # remove the temporary folder
    remove_flag && rm(dir_tmp, force=true, recursive=true)

    # apply muting
    if mute_index > 0
       p  = reshape(p, fidiff.nz, fidiff.nx)
       p[1:mute_index,:] .= 0.0
       p  = vec(p)
    end

    # apply preconditioner to model parameter
    if normalization_flag
       path_normalization = joinpath(dir_sourceside, "normalization.rsf")
       (hdr, scale) = read_RSdata(path_normalization)
       p .*= vec(scale)
    end

    hdr = RegularSampleHeader(p, title="image")
    write_RSdata(path_m, hdr, p)

    return nothing
end


text:
export evaluate, evaluate!, weylnorm, jacobian, jacobian!, dt, dt!,
    nvariables, homogenize, dehomogenize, ishomogenized, ishomogenous

"""
    evaluate(H::AbstractPolynomialHomotopy, x, t)

Evaluate the homotopy `H` at `x` to time `t`, i.e. ``H(x,t)``.

    evaluate(H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Evaluate the homotopy `H` at `x` to time `t` using the precompuated values in `cfg`.
Note that this is significantly faster than `evaluate(H, x, t)`.
"""
function evaluate end

"""
    evaluate!(u::Vector, H::AbstractPolynomialHomotopy, x, t)

Evaluate the homotopy `H` at `x` to time `t`, i.e. ``H(x,t)``, and store the result in `u`.

    evaluate!(u::AbstractVector, H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Evaluate the homotopy `H` at `x` to time `t` using the precompuated values in `cfg` and store
the result in `u`.
"""
function evaluate! end

"""
    weylnorm(H::AbstractPolynomialHomotopy)

Creates a function with variable `t` that computes the Weyl norm (or Bombieri norm) of ``H(x,t)``.
See [here](https://en.wikipedia.org/wiki/Bombieri_norm) for details about the Weyl norm.
"""
function weylnorm end

"""
    jacobian(H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute the jacobian of `H` at `x` and `t` using the precomputed values in `cfg`.
The jacobian is constructed w.r.t. `x`, i.e. it doesn't contain the partial derivatives
w.r.t. `t`.
"""
function jacobian end

"""
    jacobian!(u, H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute the jacobian of `H` at `x` and `t` using the precomputed values in `cfg` and
store the result in `u`.

    jacobian!(r::JacobianDiffResult, H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute ``H(x, t)`` and the jacobian of `H` at `x` and `t` at once using the precomputated values in `cfg`
and store thre result in `r`. This is faster than computing both values separetely.

### Example
```julia
cfg = PolynomialHomotopyConfig(H)
r = JacobianDiffResult(cfg)
jacobian!(r, H, x, t, cfg)

value(r) == H(x, t)
jacobian(r) == jacobian(H, x, t, cfg)
```
"""
function jacobian! end

"""
    dt(H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute the derivative of `H` w.r.t. ``t`` at `x` and `t` using the precomputed values in `cfg`.
"""
function dt end

"""
    dt!(u, H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute the derivative of `H` w.r.t. ``t`` at `x` and `t` using the precomputed values in `cfg`
and store the result in `u`.

    dt!(r::DtDiffResult, H::AbstractPolynomialHomotopy, x, t, cfg::PolynomialHomotopyConfig)

Compute the derivative of `H` w.r.t. ``t`` at `x` and `t` using the precomputed values in `cfg`
and store the result in `r`. This is faster than computing both values separetely.

### Example
```julia
cfg = PolynomialHomotopyConfig(H)
r = DtDiffResult(cfg)
dt!(r, H, x, t, cfg)

value(r) == H(x, t)
dt(r) == dt(H, x, t, cfg)
```
"""
function dt! end

"""
    nvariables(H::AbstractPolynomialHomotopy)

The number of variables which `H` expects as input, i.e. to evaluate `H(x,t)` `x` has to be a
vector of length `nvariables(H)`.
"""
function nvariables end

"""
    homogenize(H::AbstractPolynomialHomotopy)

Homogenize the homotopy `H`. This adds an additional variable.
If `H` is already homogenized, this is the identity.
"""
function homogenize end

"""
    dehomogenize(H::AbstractPolynomialHomotopy)

Dehomogenize the homotopy `H`. This removes the first variable.
If `H` is not homogenized, this is the identity.
"""
function dehomogenize end

"""
    ishomogenized(H::AbstractPolynomialHomotopy)

Check whether the homotopy `H` was homogenized.
"""
function ishomogenized end

"""
    ishomogenous(H::AbstractPolynomialHomotopy)

Check whether the homotopy `H` is homogenous. This does not imply that `H` was homogenized.
"""
function ishomogenous end


text:
using WordCloud
mask = rendertext("World", 1000, border=10, color=0.9, backgroundcolor=0.98, type=:svg, font="Sans Bold")
words = repeat(["we", "are", "the", "world"], 150)
weights = repeat([1], length(words))
wc = wordcloud(
        words, weights,
        mask=mask,
        angles=0,
        colors=("#006BB0", "#EFA90D", "#1D1815", "#059341", "#DC2F1F"),
        density=0.55,
        ) |> generate!
println("results are saved to lettermask.svg")
paint(wc, "lettermask.svg" , background=false)
wc
#eval# runexample(:lettermask)
#md# ![](lettermask.svg)


text:
# # Example gallery
#
# Semi-curated collection of examples.
#
# ## Lines and markers
#
# ### A simple scatter plot

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

df = (x=rand(100), y=rand(100))
xy = data(df) * mapping(:x, :y)
draw(xy)

# ### A simple lines plot

x = range(-π, π, length=100)
y = sin.(x)
df = (; x, y)
xy = data(df) * mapping(:x, :y)
layer = visual(Lines)
draw(layer * xy)

# ### Lines and scatter combined plot

x = range(-π, π, length=100)
y = sin.(x)
df = (; x, y)
xy = data(df) * mapping(:x, :y)
layers = visual(Scatter) + visual(Lines)
draw(layers * xy)

#

x = range(-π, π, length=100)
y = sin.(x)
df1 = (; x, y)
df2 = (x=rand(10), y=rand(10))
layers = data(df1) * visual(Lines) + data(df2) * visual(Scatter)
draw(layers * mapping(:x, :y))

# ### Linear regression on a scatter plot

df = (x=rand(100), y=rand(100), z=rand(100))
xy = data(df) * mapping(:x, :y)
layers = linear() + visual(Scatter) * mapping(color=:z)
draw(layers * xy)

# ## Faceting
#
# ### Facet grid

df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
plt = data(df) * mapping(:x, :y, col=:i, row=:j)
draw(plt)

# ### Facet wrap

df = (x=rand(100), y=rand(100), l=rand(["a", "b", "c", "d", "e"], 100))
plt = data(df) * mapping(:x, :y, layout=:l)
draw(plt)

# ### Embedding facets
#
# All AlgebraOfGraphics plots can be inserted in any figure position, where the rest
# of the figure is managed by vanilla Makie.
# For example

resolution = (800, 400)
fig = Figure(; resolution)
ax = Axis(fig[1, 1], title="Some plot")

df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
plt = data(df) * mapping(:x, :y, col=:i, row=:j)

subfig = fig[1, 2:3]
ag = draw!(subfig, plt)
for ae in ag
    Axis(ae).xticklabelrotation[] = π/2
end
fig

# ### Adding traces to only some subplots

df1 = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
df2 = (x=[0, 1], y=[0.5, 0.5], i=fill("a", 2), j=fill("e", 2))
layers = data(df1) * visual(Scatter) + data(df2) * visual(Lines)
draw(layers * mapping(:x, :y, col=:i, row=:j))

# ## Statistical analyses
#
# ### Density plot

using AlgebraOfGraphics: density
df = (x=randn(1000), c=rand(["a", "b"], 1000))
plt = data(df) * mapping(:x, color=:c) * density(bandwidth=0.5)
draw(plt)

#

df = (x=randn(1000), c=rand(["a", "b"], 1000))
plt = data(df) * mapping(:x, color=:c) * density(bandwidth=0.5) * visual(orientation=:vertical)
"Not yet supported" # hide

# ## Discrete scales
#
# By default categorical ticks, as well as names from legend entries, are taken from the
# value of the variable converted to a string. Scales can be equipped with labels to
# overwrite that

df = (x=rand(["a", "b", "c"], 100), y=rand(100))
plt = data(df) * mapping(:x, :y) * visual(BoxPlot)
draw(plt)

#

df = (x=rand(["a", "b", "c"], 100), y=rand(100))
plt = data(df) *
    mapping(
        :x => renamer("a" => "label1", "b" => "label2", "c" => "label3"),
        :y
    ) * visual(BoxPlot)
draw(plt)

# The order can also be changed by tweaking the scale

plt = data(df) *
    mapping(
        :x => renamer("b" => "label b", "a" => "label a", "c" => "label c"),
        :y
    ) * visual(BoxPlot)
draw(plt)

# ## Continuous scales

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => log => "√x + 20x + 100 (log scale)",
    ) * visual(Lines)
draw(plt)

#

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => "√x + 20x + 100 (log scale)",
    ) * visual(Lines)
draw(plt, axis=(yscale=log,))

#

x = 0:100
y = @. 0.01 + x/1000
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => "y",
    ) * visual(Lines)
draw(plt, axis=(yscale=log,))

# ## Custom scales
#
# Sometimes, there is no default palettes for a specific attribute. In that
# case, the user can pass their own. TODO: allow legend to use custom attribute
# of plot, such as the arrowhead or the arrowcolor and pass correct legend symbol.

using Colors
x=repeat(1:20, inner=20)
y=repeat(1:20, outer=20)
u=cos.(x)
v=sin.(y)
c=rand(Bool, length(x))
d=rand(Bool, length(x))
df = (; x, y, u, v, c, d)
colors = [colorant"#E24A33", colorant"#348ABD"]
heads = ['▲', '●']
plt = data(df) *
    mapping(:x, :y, :u, :v) *
    mapping(arrowhead=:c => nonnumeric) *
    mapping(arrowcolor=:d => nonnumeric) *
    visual(Arrows, arrowsize=10, lengthscale=0.3)
draw(plt; palettes=(arrowcolor=colors, arrowhead=heads))

# ## Axis and figure keywords
#
# ### Axis tweaking
#
# To tweak one or more axes, simply use the `axis` keyword when plotting. For example

df = (x=rand(100), y=rand(100), z=rand(100))
layers = linear() + mapping(color=:z)
plt = data(df) * layers * mapping(:x, :y)
draw(plt, axis=(aspect=1,))

# ### Figure tweaking

df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
xyc = data(df) * mapping(:x, :y, layout=:c)
layers = linear() + mapping(color=:z)
plt = xyc * layers
draw(plt, axis=(aspect=1,), figure=(resolution=(800, 400),))

# ## Wide data

df = (a=randn(100), b=randn(100), c=randn(100))
labels = ["Trace 1", "Trace 2", "Trace 3"]
plt = data(df) *
    density() *
    mapping((:a, :b, :c) .=> "some label") *
    mapping(color=dims(1) => c -> labels[c])
draw(plt)

#

df = (a=rand(100), b=rand(100), c=rand(100), d=rand(100))
labels = ["Trace 1", "Trace 2", "Trace 3"]
layers = linear() + visual(Scatter)
plt = data(df) * layers * mapping(1, 2:4, color=dims(1) => c -> labels[c])
draw(plt)

# The wide format is combined with broadcast semantics.

df = (sepal_length=rand(100), sepal_width=rand(100), petal_length=rand(100), petal_width=rand(100))
xvars = ["sepal_length", "sepal_width"]
yvars = ["petal_length" "petal_width"]
layers = linear() + visual(Scatter)
plt = data(df) *
    layers *
    mapping(
        xvars .=> "sepal",
        yvars .=> "petal",
        row=dims(1) => c -> split(xvars[c], '_')[2],
        col=dims(2) => c -> split(yvars[c], '_')[2],
    )
draw(plt)

# ## Time series

using Dates

x = today() - Year(1) : Day(1) : today()
y = cumsum(randn(length(x)))
z = cumsum(randn(length(x)))
df = (; x, y, z)
labels = ["series 1", "series 2", "series 3", "series 4", "series 5"]
plt = data(df) *
    mapping(:x, [:y, :z], color=dims(1)=>(c -> labels[c])=>"series ") *
    visual(Lines)
draw(plt)

#

x = now() - Hour(6) : Minute(1) : now()
y = cumsum(randn(length(x)))
z = cumsum(randn(length(x)))
df = (; x, y, z)
plt = data(df) *
    mapping(:x, [:y, :z], color=dims(1)=>(c -> labels[c])=>"series ") *
    visual(Lines)
draw(plt)

# ## Geometries

using GeometryBasics

geometry = [Rect(Vec(i, j), Vec(1, 1)) for i in 0:7 for j in 0:7]
group = [isodd(i + j) ? "light square" : "dark square" for i in 0:7 for j in 0:7]
df = (; geometry, group)

plt = data(df) * visual(Poly) * mapping(:geometry, color = :group)
draw(plt; axis=(aspect=1,))

# ## New columns on the fly

df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
layers = linear() + mapping(color=:z)
plt = data(df) * layers * mapping(:x, (:x, :y, :z) => (+) => "x + y + z", layout=:c)
draw(plt)

# ## Legend merging

N = 40

x = [1:N; 1:N]
y = [cumsum(randn(N)); cumsum(randn(N))]
grp = [fill("a", N); fill("b", N)]

df = (; x, y, grp)

layers = visual(Lines) + visual(Scatter) * mapping(marker = :grp)
plt = data(df) * layers * mapping(:x, :y, color = :grp)

draw(plt)

#


text:
"""
    Compressor

Representation of compressor supprted by blosc
"""
struct Compressor
    name ::Symbol
    libname ::String
    code ::Int8
    version ::String
end

function fill_available_compressors!(dict::Dict{Symbol, Compressor})
    names = split(
        unsafe_string(
            Lib.blosc_list_compressors()
        ), ","
    )
    for n in names
        complib = Ref{Cstring}()
        version = Ref{Cstring}()
        code = Lib.blosc_get_complib_info(n, complib, version)
        code < 0 && continue
        dict[Symbol(n)] = Compressor(Symbol(n), unsafe_string(complib[]), code, unsafe_string(version[]))
    end
end

function compname_by_code(code)
    compname = Ref{Cstring}()
    r = Lib.blosc_compcode_to_compname(code, compname)
    (r == -1 || compname[] == C_NULL) && error("compressor code $compcode is not recognized, or there is not support for it in this build")
    return unsafe_string(compname[])
end

function _get_default_compressor_name()
    name = compname_by_code(0)
    return Symbol(name)
end

"""
    available_compressors()::Dict{Symbol, Compressor}

Get dictionary of available compressors
"""
available_compressors()::Dict{Symbol, Compressor} = _available_compressors

"""
    available_compressors_names()::Vector{Symbol}

Get vector of available compressors names
"""
available_compressors_names() = collect(keys(_available_compressors))

"""
    has_compressor(name::Symbol)

Check if compressor `name` available
"""
has_compressor(name::Symbol) = haskey(_available_compressors, name)

"""
    compressor(name::Symbol)

Get compressor by name
"""
function compressor_by_name(name::Symbol)
    !has_compressor(name) && throw(KeyError(name))
    _available_compressors[name]
end

"""
    default_compressor_name()

Name of default compressor
"""
default_compressor_name() = _default_compressor_name

"""
    default_compressor()

Default compressor
"""
default_compressor() = _available_compressors[_default_compressor_name]


text:
using JuliaFormatter

format(file; kwargs...) = JuliaFormatter.format(joinpath(@__DIR__, file); kwargs...)

format("src"; verbose = true)
format("test"; verbose = true)


text:
using Test
using LSH
using DelimitedFiles
using SparseArrays

function test_typed(T)
    @info("testing with type $T")

    # testinput.csv prepared from document-term matrix of few SPDX files:
    # BSD-1-Clause.txt, GPL-1.0-only.txt, LGPL-2.0-only.txt, MIT.txt, Vim.txt
    dtm = convert(SparseMatrixCSC{T,T}, sparse(convert(Array{T,2}, readdlm("testinput.csv", ','))))
    model = LSHModel(20, dtm)
    iob = IOBuffer()
    show(iob, model)
    @test String(take!(iob)) == "LSHModel{$T}(nhashes=20, ndocs=5, nterms=9554)"

    @info("testing signatures...")
    signatures = signature(model, dtm)
    @test signatures == model.dochashes

    nfailures = LSH.selftest(model)
    @test nfailures == 0

    allmatches = LSH.match(model, model.dochashes)
    for idx in 1:length(allmatches)
        @test first(allmatches[idx]).id == idx
    end
end

test_typed(Int64)
test_typed(Int32)


text:
type PartialQuotients{T <: Union{QuadraticIrrational}}
  number::T
  partialquotients::Vector{Int}
  accuracy::Int
  periodlength::Int
end

function PartialQuotients(q::QuadraticIrrational; accuracy = 100)
  continued_frac = PartialQuotients(q, Vector{Int}(), accuracy, 0)
  computepartialquotients!(continued_frac)
  return continued_frac
end

function show(io::IO, cf::PartialQuotients)
  println(io, "Partial quotients of $(cf.number)")
  println(io, " > $(cf.partialquotients)")
  print(io, " > period length is $(cf.periodlength)")
end

function computepartialquotients!(cf::PartialQuotients{QuadraticIrrational})
  x = cf.number
  accuracy = cf.accuracy
  while true
    period = findperiod(cf)
    if accuracy > 0 && period == 0
      partialquotient = Int(floor(eval(x)))
      if eval(x) > partialquotient
        push!(cf.partialquotients, partialquotient)
        next_x = 1 / (x - partialquotient)
      else
        push!(cf.partialquotients, partialquotient)
        cf.periodlength = 0
        return nothing
      end
    else
      cf.periodlength = period
      return nothing
    end
    accuracy -= 1
    x = next_x
  end
end


text:
@testset "Solvers" begin
    binary = ""
    if Sys.iswindows()
        binary = joinpath(pwd(), "solvers/bin/radius.exe")
    elseif Sys.isapple()
        binary = joinpath(pwd(), "solvers/bin/radius-mac")
    else
        binary = joinpath(pwd(), "solvers/bin/radius")
    end

    solver = Solver(binary, "", "in.txt")

    tmp = tempdir()

    open(joinpath(tmp, "in.txt"), "w") do input
        println(input, "0.5")
        println(input, "0.5")
    end

    run(solver, tmp)

    radius = map(x -> parse(Float64, x), readlines(joinpath(tmp, "out.txt")))

    @test radius[1] == sqrt(0.5^2 + 0.5^2)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
function find_dipoles(vi::VolumeImage; kwargs...)

    Logging.info("Finding dipoles for volume image")

    if size(vi.data, 4) > 1
        warn("Can not squeeze 4d volume image to 3d. Please reduce first.")
    end

    x = [xi / (1 * Meter) for xi in vi.x]
    y = [yi / (1 * Meter) for yi in vi.y]
    z = [zi / (1 * Meter) for zi in vi.z]

    # List comprehension returns type any which needs to be changed
    x = convert(Array{AbstractFloat}, x)
    y = convert(Array{AbstractFloat}, y)
    z = convert(Array{AbstractFloat}, z)

    unique_dipoles(find_dipoles(squeeze(vi.data, 4), x=x, y=y, z=z; kwargs...))
end

function new_dipole_method(vi::VolumeImage; min_size::Real = 1, kwargs...)

    old_dips = find_dipoles(vi; kwargs...)

    new_dips = Dipole[]

    for dip in old_dips

        threshold = 0.9 * dip.size

        tmp_vi = deepcopy(vi)

        tmp_vi.data[tmp_vi.data .< threshold] = 0

        val_x = abs(tmp_vi.x .- (dip.x) ) .> (0.015 * Meter)
        val_y = abs(tmp_vi.y .- (dip.y) ) .> (0.015 * Meter)
        val_z = abs(tmp_vi.z .- (dip.z) ) .> (0.015 * Meter)

        tmp_vi.data[val_x, :, :] = 0
        tmp_vi.data[:, val_y, :] = 0
        tmp_vi.data[:, :, val_z] = 0

        valid = tmp_vi.data .> threshold;

        x_loc = mean(vi.x[find(squeeze(sum(valid, [2, 3]), (2, 3)))]  / (1.0 * SIUnits.Meter))
        y_loc = mean(vi.y[find(squeeze(sum(valid, [1, 3]), (1, 3)))]  / (1.0 * SIUnits.Meter))
        z_loc = mean(vi.z[find(squeeze(sum(valid, [1, 2]), (1, 2)))]  / (1.0 * SIUnits.Meter))

        x, y, z, t = find_location(vi, x_loc, y_loc, z_loc)
        s = vi.data[x, y, z, t]

        push!(new_dips, Dipole("Talairach", x_loc, y_loc, z_loc, 0, 0, 0, 0, 0, s))
    end

    new_dips = new_dips[find([d.size > min_size for d in new_dips])]

    unique_dipoles(new_dips)
end

unique_dipoles(dips = Array{Dipoles}) = dips[find(![false; diff([d.size for d in dips]) .== 0])]

lowest_dipole(dips = Array{Dipoles}) = dips[find([d.z for d in dips] .== minimum([d.z for d in dips]))]


text:
mutable struct QM_FloatData_ref{T <: Real, S, Ssp} <: Abstract_QM_FloatData{T, S, Ssp}
  Q::Ssp
  A::Ssp # using AT is easier to form systems
  b::S
  c::S
  c0::T
  lvar::S
  uvar::S
  Δref::T
  x_approx::S
  x_approxQx_approx_2::T
  b_init::S
  bTy_approx::T
  c_init::S
  pri_obj_approx::T
  uplo::Symbol
end

function fd_refinement(
  fd::QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  Δxy::Vector{T},
  pt::Point{T},
  itd::IterData{T},
  ϵ::Tolerances{T},
  dda::DescentDirectionAllocs{T},
  pad::PreallocatedData{T},
  spd::StartingPointData{T},
  cnts::Counters,
  T0::DataType,
  refinement::Symbol;
  centering::Bool = false,
) where {T <: Real}

  # center Points before zoom
  if centering
    # if pad.fact_fail # re-factorize if it failed in a lower precision system
    #     out = factorize_K2!(pad.K, pad.K_fact, pad.D, pad.diag_Q, pad.diagind_K, pad.regu,
    #                         pt.s_l, pt.s_u, itd.x_m_lvar, itd.uvar_m_x, id.ilow, id.iupp,
    #                         id.ncon, id.nvar, cnts, itd.qp, T, T0)
    # end
    # dda.rxs_l .= -itd.μ
    # dda.rxs_u .= itd.μ
    # solve_augmented_system_cc!(pad.K_fact, itd.Δxy, itd.Δs_l, itd.Δs_u, itd.x_m_lvar, itd.uvar_m_x,
    #                         pad.rxs_l, pad.rxs_u, pt.s_l, pt.s_u, id.ilow, id.iupp)
    itd.Δxy .= 0
    itd.Δxy[id.ilow] .-= itd.μ ./ itd.x_m_lvar
    itd.Δxy[id.iupp] .+= itd.μ ./ itd.uvar_m_x
    if pad.fact_fail
      out = update_pad!(pad, dda, pt, itd, fd, id, res, cnts, T0)
      out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :IPF)
    else
      out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :cc)
    end
    itd.Δs_l .= @views .-(.-itd.μ .+ pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar
    itd.Δs_u .= @views (itd.μ .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x

    α_pri, α_dual =
      compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)
    update_data!(pt, α_pri, α_dual, itd, pad, res, fd, id)
  end

  c_ref = fd.c .- itd.ATy .+ itd.Qx
  if refinement == :zoom || refinement == :multizoom
    # zoom parameter
    Δref = one(T) / res.rbNorm
  elseif refinement == :ref || refinement == :multiref
    Δref = one(T)
    αref = T(1.0e12)
    δd = norm(c_ref, Inf)
    if id.nlow == 0 && id.nupp > 0
      δp = max(res.rbNorm, maximum(itd.uvar_m_x))
    elseif id.nlow > 0 && id.nupp == 0
      δp = max(res.rbNorm, maximum(itd.x_m_lvar))
    elseif id.nlow == 0 && id.nupp == 0
      δp = max(res.rbNorm)
    else
      δp = max(res.rbNorm, maximum(itd.x_m_lvar), maximum(itd.uvar_m_x))
    end
    Δref = max(αref / Δref, one(T) / δp, one(T) / δd)
  end

  c_ref .*= Δref
  fd_ref = QM_FloatData_ref(
    fd.Q,
    fd.A,
    .-res.rb .* Δref,
    c_ref,
    fd.c0,
    Δref .* (fd.lvar .- pt.x),
    Δref .* (fd.uvar .- pt.x),
    Δref,
    pt.x,
    copy(itd.xTQx_2),
    fd.b,
    dot(fd.b, pt.y),
    fd.c,
    copy(itd.pri_obj),
    fd.uplo,
  )

  # init zoom Point
  S = typeof(fd.c)
  pt_z = Point(
    fill!(S(undef, id.nvar), zero(T)),
    fill!(S(undef, id.ncon), zero(T)),
    max.(abs.(c_ref[id.ilow]), eps(T)),
    max.(abs.(c_ref[id.iupp]), eps(T)),
  )
  starting_points!(pt_z, fd_ref, id, itd, spd)

  # update residuals
  res.rb .= itd.Ax .- fd_ref.b
  res.rc .= itd.ATy .- itd.Qx .- fd_ref.c
  res.rc[id.ilow] .+= pt_z.s_l
  res.rc[id.iupp] .-= pt_z.s_u
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)

  return fd_ref, pt_z
end

function update_pt_ref!(
  Δref::T,
  pt::Point{T},
  pt_z::Point{T},
  res::AbstractResiduals{T},
  id::QM_IntData,
  fd::QM_FloatData{T},
  itd::IterData{T},
) where {T <: Real}

  # update Point
  pt.x .+= pt_z.x ./ Δref
  pt.y .+= pt_z.y ./ Δref
  pt.s_l .= pt_z.s_l ./ Δref
  pt.s_u .= pt_z.s_u ./ Δref

  # update IterData
  itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, fd.uplo), pt.x)
  itd.xTQx_2 = dot(pt.x, itd.Qx) / 2
  fd.uplo == :U ? mul!(itd.ATy, fd.A, pt.y) : mul!(itd.ATy, fd.A', pt.y)
  fd.uplo == :U ? mul!(itd.Ax, fd.A', pt.x) : mul!(itd.Ax, fd.A, pt.x)
  itd.cTx = dot(fd.c, pt.x)
  itd.pri_obj = itd.xTQx_2 + itd.cTx + fd.c0
  itd.dual_obj = @views dot(fd.b, pt.y) - itd.xTQx_2 + dot(pt.s_l, fd.lvar[id.ilow]) -
         dot(pt.s_u, fd.uvar[id.iupp]) + fd.c0
  itd.pdd = abs(itd.pri_obj - itd.dual_obj) / (one(T) + abs(itd.pri_obj))

  # update Residuals
  res.rb .= itd.Ax .- fd.b
  res.rc .= itd.ATy .- itd.Qx .- fd.c
  res.rc[id.ilow] .+= pt.s_l
  res.rc[id.iupp] .-= pt.s_u
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
end

function update_data!(
  pt::Point{T},
  α_pri::T,
  α_dual::T,
  itd::IterData{T},
  pad::PreallocatedData{T},
  res::AbstractResiduals{T},
  fd::QM_FloatData_ref{T},
  id::QM_IntData,
) where {T <: Real}

  # (x, y, s_l, s_u) += α * Δ
  update_pt!(
    pt.x,
    pt.y,
    pt.s_l,
    pt.s_u,
    α_pri,
    α_dual,
    itd.Δxy,
    itd.Δs_l,
    itd.Δs_u,
    id.ncon,
    id.nvar,
  )

  # update IterData
  itd.x_m_lvar .= @views pt.x[id.ilow] .- fd.lvar[id.ilow]
  itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt.x[id.iupp]
  boundary_safety!(itd.x_m_lvar, itd.uvar_m_x)

  itd.μ = compute_μ(itd.x_m_lvar, itd.uvar_m_x, pt.s_l, pt.s_u, id.nlow, id.nupp)
  itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, :U), pt.x)
  x_approxQx = dot(fd.x_approx, itd.Qx)
  itd.xTQx_2 = dot(pt.x, itd.Qx) / 2
  fd.uplo == :U ? mul!(itd.ATy, fd.A, pt.y) : mul!(itd.ATy, fd.A', pt.y)
  fd.uplo == :U ? mul!(itd.Ax, fd.A', pt.x) : mul!(itd.Ax, fd.A, pt.x)
  itd.cTx = dot(fd.c_init, pt.x)
  itd.pri_obj =
    fd.pri_obj_approx +
    x_approxQx / fd.Δref +
    itd.xTQx_2 / fd.Δref^2 +
    dot(fd.c_init, pt.x) / fd.Δref
  itd.dual_obj = @views fd.bTy_approx + dot(fd.b_init, pt.y) / fd.Δref - fd.x_approxQx_approx_2 -
         dot(fd.x_approx, itd.Qx) / fd.Δref - itd.xTQx_2 / fd.Δref^2 +
         dot(pt.s_l, fd.lvar[id.ilow]) / fd.Δref^2 +
         dot(pt.s_l, fd.x_approx[id.ilow]) / fd.Δref - dot(pt.s_u, fd.uvar[id.iupp]) / fd.Δref^2 -
         dot(pt.s_u, fd.x_approx[id.iupp]) / fd.Δref + fd.c0
  itd.pdd = abs(itd.pri_obj - itd.dual_obj) / (one(T) + abs(itd.pri_obj))

  #update Residuals
  res.rb .= itd.Ax .- fd.b
  res.rc .= itd.ATy .- itd.Qx .- fd.c
  res.rc[id.ilow] .+= pt.s_l
  res.rc[id.iupp] .-= pt.s_u
  res.rcNorm, res.rbNorm = norm(res.rc, Inf) / fd.Δref, norm(res.rb, Inf) / fd.Δref
end


text: hswish(x) = x * (relu6(x + 3) / 6)


text:
using ThArrays
using Random
using Test

@testset "Tensor Gradient" begin

    Random.seed!(0);

    a = rand(3, 2)
    b = rand(3, 2)

    ad_reseult = [4.647295015954825 3.354657692932529;
                  4.820713075852873 3.557760218662402;
                  3.3291315962673704 3.4069531160838453]

    @testset "Simple AD" begin
        ta = Tensor(a, requires_grad=true)
        tb = Tensor(b)
        tc = ta^2 + 3ta + sin(tb) - tb

        tg = Tensor(ones(3, 2))

        ThAD.backward(tc, tg)

        @test ThC.grad(ta) == Tensor(ad_reseult)
    end

    @testset "ThAD.gradient" begin
        f(x, y) = x^2 + 3x + sin(y) - y
        grads = ThAD.gradient(f, a, b; d=Tensor(ones(3,2)))
        @test grads[1] == Tensor(ad_reseult)
    end

    @testset "Reset gradient" begin
        t = Tensor(a, requires_grad=true)
        grads = ThAD.gradient(x -> sum(2x), t)
        @test grads[1] == Tensor(ones(3, 2)) * 2
        grads = ThAD.gradient(x -> sum(2x), t)
        @test grads[1] == Tensor(ones(3, 2)) * 4

        ThAD.reset_grad!(t)
        grads = ThAD.gradient(x -> sum(2x), t)
        @test grads[1] == Tensor(ones(3, 2)) * 2
    end

    @testset "ThAD.forward" begin
        f(x, y) = sum(2x + 2y)
        y, back = ThAD.forward(f, a, b)
        grads = back(1)
        @test grads[1] == grads[2] == Tensor(ones(3, 2)) * 2
    end

end


text:
# ## [Time Domain Simulation of a Typical Section](@id section-simulation)
#
# In this example, we simulate the response of two degree of freedom typical section model.
# We use the same parameters as in the [previous example](@ref section-stability)
#
# ![](../assets/section-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`section-simulation.ipynb`](@__NBVIEWER_ROOT_URL__/examples/section-simulation.ipynb).
#-

using Aeroelasticity, DifferentialEquations, LinearAlgebra

## define non-dimensional parameters
V = range(1e-6, 3.1, length=1000) # = U/(b*ωθ) (reduced velocity)
a = -1/5 # reference point normalized location
e = -1/10 # center of mass normalized location
μ = 20 # = m/(ρ*pi*b^2) (mass ratio)
r2 = 6/25 # = Iθ/(m*b^2) (radius of gyration about P)
σ = 2/5 # = ωh/ωθ (natural frequency ratio)
xθ = e - a # distance from center of mass to reference point
a0 = 2*pi # lift curve slope
α0 = 0 # zero lift angle
cd0 = 0 # drag coefficient
cm0 = 0 # moment coefficient

## choose dimensional parameters
b = 1 # semichord
ρ = 1 # air density
ωθ = 1 # pitch natural frequency
c = 343 # air speed of sound

## calculate dimensionalized parameters
U = V*b*ωθ # freestrean velocity
m = μ*ρ*pi*b^2 # mass
Sθ = m*xθ*b # mass imbalance
Iθ = r2*m*b^2 # inertia
ωh = σ*ωθ # plunge natural frequency
kh = m*ωh^2 # plunge spring constant
kθ = Iθ*ωθ^2 # pitch spring constant

## reduced velocity
V = 1.0 # = U/(b*ωθ)

## dimensionalized velocity
U = V*b*ωθ

## define coupled model
model = assemble_model(;
    aerodynamic_model = Peters{6}(),
    structural_model = Section())

## define parameter vector
p = assemble_parameters(model;
    aerodynamic_parameters = (; a = a, b = b, a0 = a0, alpha0 = α0, cd0 = cd0, cm0 = cm0),
    structural_parameters = (; kh = kh, ktheta = kθ, m = m, Stheta = Sθ, Itheta = Iθ),
    additional_parameters = (; U = U, rho = ρ, c = c))

## construct ODE function
f = ODEFunction(model)

## initial states
x0 = assemble_states(model;
    aerodynamic_states = (;lambda=zeros(6)),
    structural_states = (;h=0.5, theta=0, hdot=0, thetadot=0))

## simulate for 100 seconds
tspan = (0.0, 100.0)

## construct ODE problem
prob = ODEProblem(f, x0, tspan, p)

## solve ODE
sol = solve(prob)

#!jl nothing #hide

# We can then plot the solution using DifferentialEquations' built-in interface with the
# [Plots](https://github.com/JuliaPlots/Plots.jl) package.

using Plots
pyplot()

plot(sol,
    vars = [7,8,9,10],
    xlabel = "t",
    ylabel = permutedims([
        "\$h\$",
        "\$\\theta\$",
        "\$\\dot{h}\$",
        "\$\\dot{\\theta}\$",
        ]),
    label = "",
    layout = (4, 1),
    size = (600,1200)
    )

#md savefig("../assets/section-simulation-solution.svg") #hide
#md nothing #hide

#md # ![](../assets/section-simulation-solution.svg)
#-
# For aeroelastic models based on a typical section, we can also easily visualize the
# section's behavior.

## animation parameters
a = -1/5
b = 0.5

## create animation
#md anim = @animate for t in range(tspan[1], tspan[2], length=200)
#!md @gif for t in range(tspan[1], tspan[2], length=200)
    h, θ = sol(t, idxs=7:8)
    xplot, yplot = section_coordinates(h, θ; a, b)
    plot(xplot, yplot;
        framestyle = :origin,
        grid = :false,
        xlims = (-1.0, 1.0),
        ylims = (-0.75, 0.75),
        aspect_ratio = 1.0,
        label = "t = $(round(t, digits=1))")
end

#md gif(anim, "../assets/section-simulation.gif")
#md nothing #hide

#md # ![](../assets/section-simulation.gif)


text:
using Tethys
using Test
using Printf

@testset "Tethys.jl" begin
# Write your own tests here.

include("BareExpansion.jl")
end


text:
"""
This file is for creating a virtual environment reflecting the learned model.
For our work, we use obs -> obs as a prediction and use obs for reward computation etc.
So, statespace and obsspace are the same; however, we can easily make this to be the case for state -> state
prediction while obs is still used for other purposes.
"""

using LyceumBase
using Shapes
using Base: @propagate_inbounds

# we assume the data stored in state, obs, action are vectors of some float type
struct LearnedEnv{D,E,K,P,
    SP,OP,AP,RP} <: AbstractEnvironment
    state::D
    obs::D
    action::D
    time::D
    mjenv::E
    PredictMat::K
    phi::P
    features::D # scratch space
    obsact::D
    nextobs::D
    sp::SP
    op::OP
    ap::AP
    rp::RP
    function LearnedEnv{T}(PredictMat, phi, mjenv) where T<:AbstractFloat
        nobs = length(obsspace(mjenv))
        naction = length(actionspace(mjenv))
        sp = VectorShape(T, nobs)
        op = VectorShape(T, nobs)
        ap = VectorShape(T, naction)
        rp = ScalarShape(T)
        s = zeros(T, nobs)
        o = zeros(T, nobs)
        a = zeros(T, naction)
        env = new{typeof(s), typeof(mjenv), typeof(PredictMat), typeof(phi),
        typeof(sp), typeof(op), typeof(ap), typeof(rp)}(s,o,a,zeros(T,1),
                                                        mjenv,
                                                        copy(PredictMat), phi,
                                                        zeros(T, size(PredictMat, 2)),
                                                        zeros(T, nobs + naction),
                                                        zeros(T, nobs),
                                                        sp,op,ap,rp)
        reset!(env)
    end
end

function LyceumBase.tconstruct(::Type{LearnedEnv}, PredictMat, phi, mjenvs, n)
    @assert length(mjenvs) >= n
    Tuple(LearnedEnv{Float64}(PredictMat, phi, mjenv) for mjenv in mjenvs)
end
LearnedEnv(PredMat, phi, mjenv) = first(tconstruct(LearnedEnv, PredMat, phi, [mjenv], 1))

@inline LyceumBase.statespace(env::LearnedEnv) = env.sp
@inline LyceumBase.obsspace(env::LearnedEnv) = env.op
@inline LyceumBase.actionspace(env::LearnedEnv) = env.ap
@inline LyceumBase.rewardspace(env::LearnedEnv) = env.rp

@inline function LyceumBase.getstate!(state, env::LearnedEnv)
    state .= 0.0
    copyto!(state, env.state)
end
@inline function LyceumBase.setstate!(env::LearnedEnv, state)
    if length(env.state) >= length(state)
        env.state .= state
    else
        nstate = length(env.state) # drops warmstart
        for i=1:nstate
            env.state[i] = state[i]
        end
    end
end

@inline LyceumBase.getaction!(action, env::LearnedEnv) = action .= env.action
@inline LyceumBase.setaction!(env::LearnedEnv, action) = env.action .= action

@propagate_inbounds function LyceumBase.isdone(state, action, obs, env::LearnedEnv)
    isdone(state, action, obs, env.mjenv)
end
@propagate_inbounds function LyceumBase.getobs!(obs, env::LearnedEnv)
    copyto!(env.obs, env.state)
    obs .= env.obs
end
@propagate_inbounds function LyceumBase.getreward(s, a, o, env::LearnedEnv)
    getreward(s, a, o, env.mjenv) # uses original env's reward function
end
@propagate_inbounds function LyceumBase.geteval(s, a, o, env::LearnedEnv)
    geteval(s, a, o, env.mjenv)
end
@propagate_inbounds function LyceumBase.reset!(env::LearnedEnv)
    reset!(env.mjenv)
    getobs!(env.obs, env.mjenv)
    nstate = length(env.state) # drops warmstart
    copyto!(env.state, 1:nstate, env.obs, 1:nstate)
    env.time[1] = 0.0
    env
end
@propagate_inbounds function LyceumBase.randreset!(env::LearnedEnv)
    randreset!(env.mjenv)
    getobs!(env.obs, env.mjenv)
    nstate = length(env.state) # drops warmstart
    copyto!(env.state, 1:nstate, env.obs, 1:nstate)
    env.time[1] = 0.0
    env
end

@propagate_inbounds function LyceumBase.step!(env::LearnedEnv)
    # the code performs the following in a memory allocation free way
    # state += PredictMat * phi(vcat(obs, act))
    obsact, obs, act = env.obsact, env.obs, env.action
    copyto!(obsact, obs)                   # copy obs to start of vector
    r = (length(obs)+1):length(obsact)
    copyto!(obsact, r, act, 1:length(act)) # copy act to end of vector
    env.phi(env.features, obsact)
    mul!(env.nextobs, env.PredictMat, env.features)
    env.state .+= env.nextobs
    setstate!(env, env.state)
    env.time[1] += timestep(env)
    env
end

@inline LyceumBase.timestep(env::LearnedEnv) = timestep(env.mjenv)
@inline Base.time(env::LearnedEnv) = env.time[1]


text:
#    Problem :
#    *********
#	 The separable extension of Rosenbrock's function.
#
#    Origonal SIF Source:problem 21 in
#    J.J. More', B.S. Garbow and K.E. Hillstrom,
#    "Testing Unconstrained Optimization Software",
#    ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
#
#    SROSENBR.SIF classification SUR2-AN-V-0
#
#    Number of variables is variable
#
# Daniel Henderson, 08/2021

f = (x) -> begin
	fx = 0.0
	for i in firstindex(x):2:lastindex(x)
		t1     = 1 - x[i]
		t2     = 10(x[i+1] - x[i]^2)
		fx    += t1^2 + t2^2
	end
    return fx
end

g! = (g, x) -> begin
	for i in firstindex(x):2:lastindex(x)
		t1     = 1 - x[i]
		t2     = 10(x[i+1] - x[i]^2)
		g[i+1] = 20t2
		g[i]   = -2(x[i] * g[i+1] + t1)
	end
    return g
end

fg! = (g, x) -> begin
	fx = 0.0
	for i in firstindex(x):2:lastindex(x)
		t1     = 1 - x[i]
		t2     = 10(x[i+1] - x[i]^2)
		g[i+1] = 20t2
		g[i]   = -2(x[i] * g[i+1] + t1)
		fx    += t1^2 + t2^2
	end
    return fx, g
end

init = (n::Int=5000) -> begin
	mod(n, 2) > 0 && @warn "SROSENBR: number of variables must be even"
	q = max(1, div(n, 2))
	n = 2q

    x0 = [j % 2 == 1 ? 1.2 : 1.0 for j in 1:n]
    return n, x0
end

TestSet["SROSENBR"] = UncProgram("SROSENBR",  f, g!, fg!, init)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_Blocks.jl")
include("model_BlocksInner.jl")
include("model_DefsAppId.jl")
include("model_DefsBotId.jl")
include("model_DefsChannel.jl")
include("model_DefsChannelId.jl")
include("model_DefsChannelName.jl")
include("model_DefsCommentId.jl")
include("model_DefsDmId.jl")
include("model_DefsEnterpriseId.jl")
include("model_DefsEnterpriseName.jl")
include("model_DefsEnterpriseUserId.jl")
include("model_DefsFileId.jl")
include("model_DefsGroupId.jl")
include("model_DefsOkFalse.jl")
include("model_DefsOkTrue.jl")
include("model_DefsOptionalAppId.jl")
include("model_DefsPinnedInfo.jl")
include("model_DefsReminderId.jl")
include("model_DefsSubteamId.jl")
include("model_DefsTeam.jl")
include("model_DefsTopicPurposeCreator.jl")
include("model_DefsTs.jl")
include("model_DefsUserId.jl")
include("model_DefsWorkspaceId.jl")
include("model_ObjsBotProfile.jl")
include("model_ObjsBotProfileIcons.jl")
include("model_ObjsChannel.jl")
include("model_ObjsChannelPurpose.jl")
include("model_ObjsComment.jl")
include("model_ObjsComments.jl")
include("model_ObjsConversation.jl")
include("model_ObjsEnterpriseUser.jl")
include("model_ObjsExternalOrgMigrations.jl")
include("model_ObjsExternalOrgMigrationsCurrent.jl")
include("model_ObjsFile.jl")
include("model_ObjsFileShares.jl")
include("model_ObjsIcon.jl")
include("model_ObjsMessage.jl")
include("model_ObjsMessageAttachments.jl")
include("model_ObjsMessageIcons.jl")
include("model_ObjsPaging.jl")
include("model_ObjsPrimaryOwner.jl")
include("model_ObjsReaction.jl")
include("model_ObjsReminder.jl")
include("model_ObjsResources.jl")
include("model_ObjsResponseMetadata.jl")
include("model_ObjsScopes.jl")
include("model_ObjsSubteam.jl")
include("model_ObjsSubteamPrefs.jl")
include("model_ObjsTeam.jl")
include("model_ObjsTeamProfileField.jl")
include("model_ObjsTeamProfileFieldOption.jl")
include("model_ObjsTeamSsoProvider.jl")
include("model_ObjsUser.jl")
include("model_ObjsUserProfile.jl")
include("model_ObjsUserProfileShort.jl")


text:
import DataFrames: DataFrame

function rast2df(r::GeoArray)
    LON, LAT = st_coords(r)
    DataFrame(id = seq_along(LON), value = r.A[:], lon = LON[:], lat = LAT[:])
end

function rast2df(r::Array{GeoArray})
    df = [];
    for i in 1:length(r)
        d = rast2df(r[i])
        push!(df, d)
    end
    df
end

"""
    shrink_bbox(ga, mask::AbstractArray{Bool, 2} = nothing)

only true values in `mask` will be kept.
"""
function shrink_bbox(ga::GeoArray, mask::Union{Nothing, AbstractArray{Bool, 2}} = nothing)
    if mask === nothing; mask = ga.A[:, :, 1] .!= 0; end
    ind = findall(mask)
    # ind_vec = LinearIndices(mat)[ind]
    rows = map(x -> x[1], ind) # x, long
    cols = map(x -> x[2], ind) # y, lat

    I_x = seq(Range(rows)...)
    I_y = seq(Range(cols)...)
    ga[I_x, I_y]
end

"""
    st_as_sf(ga::GeoArray, mask::Union{Nothing, AbstractArray{Bool, 2}} = nothing)
    st_as_sf(file::AbstractString, mask = nothing)

If the input is a file path, `shrink_bbox` will be applied.
"""
function st_as_sf(ga::GeoArray, mask::Union{Nothing, AbstractArray{Bool, 2}} = nothing)
    if mask === nothing; mask = ga.A[:, :, 1] .!= 0; end
    ind = findall(mask)
    # ind_vec = LinearIndices(mat)[ind]
    ntime = size(ga, 3)
    res = zeros(eltype(ga.A), length(ind), ntime)
    @views for i = 1:ntime
        mat = ga.A[:, :, i]
        res[:, i] = mat[ind]
    end
    # keep enough information for the reverse operation
    res, rast(ga, vals = mask) # mat, mask
end

function st_as_sf(file::AbstractString)
    ga = rast(file)
    ga2 = shrink_bbox(ga)
    st_as_sf(ga2)
end

# mask is 3d boolean array
function st_as_sf(file::AbstractString, mask_shrink::AbstractArray{Bool}, mask::AbstractArray{Bool})
    ga = rast(file)
    ga2 = shrink_bbox(ga, mask_shrink)
    st_as_sf(ga2, mask)
end

function maskCoords(mask::GeoArray{Bool})
    LON, LAT = st_coords(mask)
    ind = findall(mask.A)
    DataFrame(id = seq_along(ind), lon = LON[ind], lat = LAT[ind])
end

export rast2df, shrink_bbox, st_as_sf, maskCoords


text:
function _nabla(ex::Expr; mod)
    def = splitdef(ex)
    symfuncs, exfuncs, exrules = Symbol[], Expr[], Expr[]

    def[:body] = MacroTools.postwalk(def[:body]) do x
        lhs, lhsind, rhs = Ref{Symbol}(), Ref{Vector{Any}}(), Ref{Expr}()
        which, opt = Ref{Symbol}(), Ref{Expr}()

        if @capture(x, @tensor _lhs_[_lhsind__] := _rhs_) ||
           @capture(x, @tensoropt _lhs_[_lhsind__] := _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[] = :assign
        elseif @capture(x, @tensor _lhs_[_lhsind__] += _rhs_) ||
               @capture(x, @tensoropt _lhs_[_lhsind__] += _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[] = :pluseq
        elseif @capture(x, @tensor _lhs_[_lhsind__] -= _rhs_) ||
               @capture(x, @tensoropt _lhs_[_lhsind__] -= _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[] = :subteq
        elseif @capture(x, @tensoropt _opt_ _lhs_[_lhsind__] := _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[], opt[] = :assign, _opt
        elseif @capture(x, @tensoropt _opt_ _lhs_[_lhsind__] += _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[], opt[] = :pluseq, _opt
        elseif @capture(x, @tensoropt _opt_ _lhs_[_lhsind__] -= _rhs_)
            lhs[], lhsind[], rhs[] = _lhs, _lhsind, _rhs
            which[], opt[] = :subteq, _opt
        elseif @capture(x, @tensor _[__] = _) ||
               @capture(x, @tensoropt _[__] = _) ||
               @capture(x, @tensoropt _ _[__] = _)
            error("use assignment with `:=` instead of inplace operation with `=`")
        else
            return x
        end

        lhs, lhsind, rhs, which = lhs[], lhsind[], rhs[], which[]

        rhsreplace, argsorig, argsdummy, isconjs, indsall = rhs_to_args(rhs)
        @gensym symfunc
        push!(symfuncs, symfunc)

        genargs = (symfunc, argsdummy, lhsind, rhsreplace, opt)

        exfunc = genfunc(genargs...)
        exfrule = genfrule(genargs...)
        exrrule = genrrule(genargs..., isconjs, indsall)

        @eval mod $(macroexpand(TensorOperations, exfunc))
        @eval mod $(macroexpand(TensorOperations, exfrule))
        @eval mod $(macroexpand(TensorOperations, exrrule))

        if which == :assign
            return :($lhs = $(Core.eval(mod, symfunc))($(argsorig...)))
        elseif which == :pluseq  # use x += y instead of x .+= y for Zygote
            return :($lhs += $(Core.eval(mod, symfunc))($(argsorig...)))
        elseif which == :subteq  # use x -= y instead of x .-= y for Zygote
            return :($lhs -= $(Core.eval(mod, symfunc))($(argsorig...)))
        end
    end

    return esc(MacroTools.combinedef(def)), symfuncs
end

macro ∇(ex)
    ex, _ = _nabla(ex; mod=@__MODULE__)
    return ex
end

macro fn∇(i, ex)
    _, fn = _nabla(ex; mod=@__MODULE__)
    return fn[i]
end


text:
import HTTP

function compat_helper()
    name = "CompatHelper"
    files_to_create = Dict{String, String}()
    files_to_create["CompatHelper.yml"] = String(HTTP.get("https://raw.githubusercontent.com/JuliaRegistries/CompatHelper.jl/master/.github/workflows/CompatHelper.yml").body)
    files_to_delete = Set{String}()
    return Workflow(name, files_to_create, files_to_delete)
end

function tag_bot()
    name = "TagBot"
    files_to_create = Dict{String, String}()
    files_to_create["TagBot.yml"] = String(HTTP.get("https://raw.githubusercontent.com/JuliaRegistries/TagBot/master/example.yml").body)
    files_to_delete = Set{String}()
    return Workflow(name, files_to_create, files_to_delete)
end


text:
mutable struct AccumulatedIntegral{Tt,Ft,Func,Wt}
    t::Tt
    ∫f::Ft
    f::Func
    c::Vector{Tt} # Quadrature roots
    w::Vector{Wt} # Quadrature weights

    function AccumulatedIntegral(f::Func, ::Type{Ft}, ::Type{Tt}=Ft;
                                 k = 3,
                                 t::Tt=zero(Tt), init=zero(Ft)) where {Ft,Func,Tt}
        c,w = gausslegendre(k)
        new{Tt,Ft,Func,eltype(w)}(t, init, f, c, w)
    end
end

function (F::AccumulatedIntegral)(a,b)
    μ = (a+b)/2
    δ = (b-a)/2
    v = zero(F.∫f)
    for (c,w) in zip(F.c,F.w)
        v += w*F.f(δ*c + μ)
    end
    δ*v
end

(F::AccumulatedIntegral)(t) = F.∫f + F(F.t, t)

function set!(F::AccumulatedIntegral, t)
    F.∫f = F(t)
    F.t = t
end

function field(F, ndt; kwargs...)
    t = timeaxis(F, ndt)
    tau = atomic_units.(t)

    Ev = atomic_units.(F.(t))
    A = AccumulatedIntegral(t -> -atomic_units(F(24.2u"as"*t)), Float64;
                            t=first(tau), kwargs...)
    nt = length(t)
    Av = zeros(Float64, nt)
    for i = 1:nt
        Av[i] = A(tau[i])
        set!(A, tau[i])
    end

    tau,Ev,Av
end


text:
module EntityComponentSystem

"""
    ECSComponent

A representation of an entity component system component.
"""
abstract type ECSComponent end

"""
    EntityKey

A generational index to an entity that will fail in usage if the entity has been destroyed.
"""
mutable struct EntityKey
    index::Int64
    generation::Int64
end

"""
    World

A collection of entities and their components.
"""
mutable struct World
    entity_keys::Array{Union{EntityKey,Nothing}}
    max_entity::Int64
    free_entities::Array{Int64,1}
    components::Dict{DataType,Any}
end

"""
    World()

Returns a new World instance with defaults.
"""
function World()
    World(Array{Union{Int64,Nothing}}(undef,100),0,[],Dict())
end

"""
    getentity(world::World,key::EntityKey)

Returns entity index if entity key is still valid.
"""
function getentity(world::World,key::EntityKey)
    entity = world.entity_keys[key.index]
    if entity == nothing || entity.generation != key.generation
        return nothing
    end
    key.index
end

"""
    register!(world::World,::Type{C}) where C <: ECSComponent

Registers storage space for a component type.
"""
function register!(world::World,::Type{C}) where C <: ECSComponent
    world.components[C] = Array{Union{C,Nothing}}(undef,1000)
    nothing
end

"""
    destroyentity!(world::World,entity::EntityKey)

Deallocates entity and memory if entity key is valid.
"""
function destroyentity!(world::World,entity::EntityKey)
    if getentity(world,entity) == nothing
        println("tried to destroy non existant entity")
        return
    end
    for x in keys(world.components)
        world.components[x][entity.index] = nothing
    end
    push!(world.free_entities,entity.index)
    # makes sure entity cant be retrieved again
    world.entity_keys[entity.index] = EntityKey(entity.index,world.entity_keys[entity.index].generation + 1)
    nothing
end

"""
    createentity!(world::World)

Allocates an entity and returns an entity key.
"""
function createentity!(world::World)
    i = if length(world.free_entities) > 0
        # If there's a free entity, use it
        pop!(world.free_entities)
    else
        # Otherwise add a new one
        world.max_entity += 1
    end

    entity_key = world.entity_keys[i]
    if entity_key == nothing
        entity_key = EntityKey(i,1)
        world.entity_keys[i] = entity_key
    else
        world.entity_keys[i] = EntityKey(i,world.entity_keys[i].generation + 1)
    end

    # clear out components for the new entity
    for x in keys(world.components)
        componentList = world.components[x]
        if i > length(componentList)
            resize!(componentList,2*length(componentList))
        end
        componentList[i] = nothing
    end
    entity_key
end

"""
    addcomponent!(world::World,entity::EntityKey,component::C) where C <: ECSComponent

Associate a component with an entity key.
"""
function addcomponent!(world::World,entity::EntityKey,component::C) where C <: ECSComponent
    if getentity(world,entity) == nothing
        println("tried to add component to non existant entity")
        return
    end
    world.components[C][entity.index] = component
    nothing
end

"""
    removecomponent!(world::World,entity::EntityKey,::Type{C}) where C <: ECSComponent

Dissociate a component of a specific type with an entity key.
"""
function removecomponent!(world::World,entity::EntityKey,::Type{C}) where C <: ECSComponent
    if getentity(world,entity) == nothing
        println("tried to remove component to non existant entity")
        return
    end
    world.components[C][entity.index] = nothing
    nothing
end

"""
    getcomponent(world::World,entity::EntityKey,::Type{C}) where C <: ECSComponent

Get a component of a specific type of an entity key.
"""
function getcomponent(world::World,entity::EntityKey,::Type{C}) where C <: ECSComponent
    if getentity(world,entity) == nothing
        println("tried to remove component to non existant entity")
        return
    end
    world.components[C][entity.index]
end

"""
    runsystem!(f,world::World,types::Array{DataType,1})

For each entity with components of passed in types, call provided function with that entity and components.
"""
function runsystem!(f,world::World,types::Array{DataType,1})
    num_types = length(types)
    for e = 1:world.max_entity
        components = []
        for t in types
            component_list = world.components[t]
            if e<=length(component_list) && component_list[e] != nothing
                # if we found a component, store it
                push!(components,component_list[e])
            else
                # if we find nothing we can stop right now
                break
            end
        end
        # if entity has all components call function
        if length(components) == num_types
            f(world.entity_keys[e],components)
        end
    end
    nothing
end

export ECSComponent,World,EntityKey,getentity,createentity!,register!,destroyentity!,
runsystem!,addcomponent!,removecomponent!,getcomponent
end


text:
@testset "Homology" begin
    cplx = SimplicialComplex(Simplex(1,2,3),
                        Simplex(2,4),
                        Simplex(3,4),
                        Simplex(5,4),
                        Simplex(6))

    h = homology(Int, cplx)
    @test eltype(h) == Tuple{Int64,Int64,Int64}
    @test grouptype(supertype(typeof(h))) == Nothing
    @test grouptype(typeof(h)) == Int
    @test length(h) == 2

    b, t, _ = group(h,0);
    @test b == 2
    @test t == 0

    itm, st = iterate(h)
    @test itm[1] == 0
    @test itm[2] == 2
    @test st[1] == 1

    itm, st = iterate(h, st)
    @test itm[1] == 1
    @test itm[2] == 1
    @test st[1] == 2

    @test iterate(h, st) === nothing

    @test betti(h) == (2,1)
    @test euler(betti(h)) == 1

    g = withgenerators(h)
    @test eltype(g) == Tuple{Int64,Int64,Int64,Dict{Chain,Int64}}
    @test length(g) == 2

    itm, st = iterate(g)
    @test itm[1] == 0
    @test itm[2] == 2
    @testset for (ch, d) in itm[4]
        @test d == 0
        @test all(k->k == 1, values(ch))
        @test all(k->k ∈ map(hash, map(Simplex, [5,6])), keys(ch))
    end

    itm, st = iterate(g, st)
    @test itm[1] == 1
    @test itm[2] == 1
    @test first(itm[4])[2] == 0
    gn =  Chain(1, Int) + (hash(Simplex(3,4)),1) + (hash(Simplex(2,4)),-1) + (hash(Simplex(2,3)),1)
    @testset for (a,b) in zip( simplify(first(itm[4])[1]), gn)
        @test a == b
    end

    @test iterate(g, st) === nothing

    @test length(generators(g)) == 2
end


text:
function UnixCompressFilter()
    return UnixCompressFilter(())
end

function filter(obj::UnixCompressFilter, arg0::NamedData)
    return jcall(obj, "filter", NamedData, (NamedData,), arg0)
end


text:
using HorizonSideRobots

function counter_peregorodok(r)
    goinNWcorner(r)
    signsign = true
    sum = 0
    while !isborder(r,Nord) || signsign
        sum += count(r)
        if !isborder(r,Nord)
            move!(r,Nord)
        else
            signsign = peregorodka(r,0)
        end
    end
    return sum
end

function goinNWcorner(r)
    ssign = true
    while !isborder(r,Sud) || ssign
        if !isborder(r,Sud)
            move!(r,Sud)
        end
        if isborder(r,Sud)
            ssign = peregorodka(r,2)
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    if !isborder(r,Nord)
        move!(r,Nord)
    else
        peregorodka(r,0)
    end
end

function peregorodka(r,side)
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide((side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide(side))
    end
    for i = 1:x
        move!(r,HorizonSide((side+1) % 4))
    end
    return sign
end

function count(r)
    c1 = 0
    c2 = 0
    StopCount = true
    while  !isborder(r,Ost) || StopCount
        if isborder(r,Sud)
            c1 += 1
        else
            if c1 != 0
                c2 += 1
            end
            c1 = 0
        end
        if !isborder(r,Ost)
            move!(r,Ost)
        else
            StopCount = false
            if c1 != 0
                c2 += 1
            end
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    return(c2)
end

counter_peregorodok(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs