text:
using IntervalArithmetic: AbstractInterval, mince
import Base: +, -, *, ∈, ⊆, rand, min, max

export Interval,
       dim, σ, center,
       vertices_list,
       isflat,
       linear_map
       constraints_list

"""
    Interval{N, IN<:AbstractInterval{N}} <: AbstractHyperrectangle{N}

Type representing an interval on the real line.
Mathematically, it is of the form

```math
[a, b] := \\{ a ≤ x ≤ b \\} ⊆ \\mathbb{R}.
```

### Fields

- `dat` -- data container for the given interval

### Notes

This type relies on the
[IntervalArithmetic.jl](https://juliaintervals.github.io/IntervalArithmetic.jl/stable/)
library for representation of intervals and arithmetic operations.

### Examples

Unidimensional intervals are symbolic representations of a real closed interval.

We can create intervals in different ways, the simpler way is to pass a pair
of numbers:

```jldoctest interval_constructor
julia> x = Interval(0.0, 1.0)
Interval{Float64, IntervalArithmetic.Interval{Float64}}([0, 1])
```
or a 2-vector:

```jldoctest interval_constructor
julia> x = Interval([0.0, 1.0])
Interval{Float64, IntervalArithmetic.Interval{Float64}}([0, 1])
```

Note that if the package `IntervalArithmetic` is loaded in the current scope,
you have to prepend `LazySets.` to the `Interval` type since there is a name
conflict otherwise.

```jldoctest interval_constructor
julia> using IntervalArithmetic
WARNING: using IntervalArithmetic.Interval in module Main conflicts with an existing identifier.

julia> x = Interval(IntervalArithmetic.Interval(0.0, 1.0))
Interval{Float64, IntervalArithmetic.Interval{Float64}}([0, 1])

julia> dim(x)
1

julia> center(x)
1-element Vector{Float64}:
 0.5
```

This type is such that the usual pairwise arithmetic operators `+`, `-`, `*` trigger
the corresponding interval arithmetic backend method, and return a new
`Interval` object. For the symbolic Minkowksi sum, use `MinkowskiSum` or `⊕`.

Interval of other numeric types can be created as well, eg. a rational interval:

```jldoctest interval_constructor
julia> Interval(0//1, 2//1)
Interval{Rational{Int64}, IntervalArithmetic.Interval{Rational{Int64}}}([0//1, 2//1])
```
"""
struct Interval{N, IN<:AbstractInterval{N}} <: AbstractHyperrectangle{N}
    dat::IN

    function Interval(dat::IN) where {N, IN<:AbstractInterval{N}}
        @assert isfinite(dat.lo) && isfinite(dat.hi) "intervals must be bounded"

        return new{N, IN}(dat)
    end
end

# constructor from two numbers with type promotion
function Interval(lo::N1, hi::N2) where {N1, N2}
    N = promote_type(N1, N2)
    Interval(IntervalArithmetic.Interval(N(lo), N(hi)))
end

# constructor from a vector
function Interval(x::AbstractVector)
    @assert length(x) == 2 "vector for Interval constructor has to be 2D"
    Interval(x[1], x[2])
end

# constructor from a single number
function Interval(x::Real)
    Interval(IntervalArithmetic.Interval(x))
end

isoperationtype(::Type{<:Interval}) = false
isconvextype(::Type{<:Interval}) = true

"""
    dim(x::Interval)

Return the ambient dimension of an interval.

### Input

- `x` -- interval

### Output

The integer 1.
"""
dim(x::Interval) = 1

"""
    σ(d::AbstractVector, x::Interval)

Return the support vector of an interval in a given direction.

### Input

- `d` -- direction
- `x` -- interval

### Output

Support vector in the given direction.
"""
function σ(d::AbstractVector, x::Interval)
    @assert length(d) == dim(x)
    N = promote_type(eltype(d), eltype(x))
    return d[1] > zero(N) ? high(x) : low(x)
end

"""
    ρ(d::AbstractVector, x::Interval)

Evaluate the support function of an interval in a given direction.

### Input

- `d` -- direction
- `x` -- interval

### Output

Evaluation of the support function in the given direction.
"""
function ρ(d::AbstractVector, x::Interval)
    @assert length(d) == dim(x)
    N = promote_type(eltype(d), eltype(x))
    return d[1] * (d[1] > zero(N) ? max(x) : min(x))
end

"""
    center(x::Interval)

Return the interval's center.

### Input

- `x` -- interval

### Output

The center, or midpoint, of ``x``.
"""
function center(x::Interval)
    return [IntervalArithmetic.mid(x.dat)]
end

"""
    -(x::Interval, y::Interval)

Return the difference of the intervals.

### Input

- `x` -- interval
- `y` -- interval

### Output

The difference of the intervals as a new `Interval` set.
"""
function -(x::Interval, y::Interval)
    return Interval(x.dat - y.dat)
end

"""
```
    *(x::Interval, y::Interval)
```

Return the product of the intervals.

### Input

- `x` -- interval
- `y` -- interval

### Output

The product of the intervals as a new `Interval` set.
"""
function *(x::Interval, y::Interval)
    return Interval(x.dat * y.dat)
end

"""
    ∈(v::AbstractVector, x::Interval))

Return whether a vector is contained in the interval.

### Input

- `v` -- one-dimensional vector
- `x` -- interval

### Output

`true` iff `x` contains `v`'s first component.
"""
function ∈(v::AbstractVector, x::Interval)
    return v[1] ∈ x.dat
end

"""
    ∈(v::Number, x::Interval)

Return whether a number is contained in the interval.

### Input

- `v` -- scalar
- `x` -- interval

### Output

`true` iff `x` contains `v`.
"""
function ∈(v::Number, x::Interval)
    return v ∈ x.dat
end

"""
    min(x::Interval)

Return the lower component of an interval.

### Input

- `x` -- interval

### Output

The lower (`lo`) component of the interval.
"""
function min(x::Interval)
    return x.dat.lo
end

"""
    max(x::Interval)

Return the higher or upper component of an interval.

### Input

- `x` -- interval

### Output

The higher (`hi`) component of the interval.
"""
function max(x::Interval)
    return x.dat.hi
end

"""
    low(x::Interval{N}) where {N}

Return the lower coordinate of an interval set.

### Input

- `x` -- interval

### Output

A vector with the lower coordinate of the interval.
"""
function low(x::Interval{N}) where {N}
    return N[x.dat.lo]
end

"""
    high(x::Interval{N}) where {N}

Return the higher coordinate of an interval set.

### Input

- `x` -- interval

### Output

A vector with the higher coordinate of the interval.
"""
function high(x::Interval{N}) where {N}
    return N[x.dat.hi]
end

"""
    an_element(x::Interval)

Return some element of an interval.

### Input

- `x` -- interval

### Output

The left border (`min(x)`) of the interval.
"""
function an_element(x::Interval)
    return [min(x)]
end

"""
    rand(::Type{Interval}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random interval.

### Input

- `Interval` -- type for dispatch
- `N`        -- (optional, default: `Float64`) numeric type
- `dim`      -- (optional, default: 1) dimension
- `rng`      -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`     -- (optional, default: `nothing`) seed for reseeding

### Output

A random interval.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
"""
function rand(::Type{Interval};
              N::Type{<:Real}=Float64,
              dim::Int=1,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    @assert dim == 1 "cannot create a random Interval of dimension $dim"
    rng = reseed(rng, seed)
    x = randn(rng, N)
    y = randn(rng, N)
    return x < y ? Interval(x, y) : Interval(y, x)
end

"""
    vertices_list(x::Interval)

Return the list of vertices of this interval.

### Input

- `x` -- interval

### Output

The list of vertices of the interval represented as two one-dimensional vectors,
or just one if the interval is degenerate (the endpoints match within the tolerance).
"""
function vertices_list(x::Interval)
    a = min(x)
    b = max(x)
    return _isapprox(a, b) ? [[a]] : [[a], [b]]
end

"""
    constraints_list(x::Interval{N}) where {N}

Return the list of constraints of the given interval.

### Input

- `x` -- interval

### Output

The list of constraints of the interval represented as two one-dimensional
half-spaces.
"""
function constraints_list(x::Interval{N}) where {N}
    constraints = Vector{LinearConstraint{N, SingleEntryVector{N}}}(undef, 2)
    e₁ = SingleEntryVector(1, 1, one(N))
    constraints[1] = HalfSpace(e₁, max(x))
    constraints[2] = HalfSpace(-e₁, -min(x))
    return constraints
end
"""
    translate(x::Interval, v::AbstractVector)

Translate (i.e., shift) an interval by a given vector.

### Input

- `x` -- interval
- `v` -- translation vector

### Output

A translated interval.

### Algorithm

We add the vector to the left and right of the interval.
"""
function translate(x::Interval, v::AbstractVector)
    @assert length(v) == dim(x) "cannot translate a $(dim(x))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return Interval(x.dat + v[1])
end

# --- AbstractHyperrectangle interface functions ---

"""
    radius_hyperrectangle(x::Interval{N}, i::Int) where {N}

Return the box radius of an interval in a given dimension.

### Input

- `x` -- interval
- `i` -- dimension index (must be `1`)

### Output

The box radius in the given dimension.
"""
function radius_hyperrectangle(x::Interval{N}, i::Int) where {N}
    @assert i == 1 "an interval is one-dimensional"
    return (max(x) - min(x)) / N(2)
end

"""
    radius_hyperrectangle(x::Interval)

Return the box radius of an interval in every dimension.

### Input

- `x` -- interval

### Output

The box radius of the interval (a one-dimensional vector).
"""
function radius_hyperrectangle(x::Interval)
    return [radius_hyperrectangle(x, 1)]
end

"""
    isflat(I::Interval)

Determine whether an interval is flat, i.e. whether its extreme values coincide.

### Input

- `I` -- interval

### Output

A boolean which is `true` if the interval is flat and `false` otherwise.

### Notes

For robustness with respect to floating-point inputs, this function relies on
the result of `isapproxzero` when applied to the diameter of the interval.
Hence, this function depends on the absolute zero tolerance `ABSZTOL`.
"""
function isflat(I::Interval)
    return isapproxzero(IntervalArithmetic.diam(I.dat))
end

"""
    plot_recipe(I::Interval{N}, [ε]=zero(N)) where {N}

Convert an interval to a pair `(x, y)` of points for plotting.

### Input

- `I` -- interval
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

A pair `(x, y)` of two points that can be plotted.

### Notes

We consider the interval as a line segment with y coordinate equal to zero.
"""
function plot_recipe(I::Interval{N}, ε=zero(N)) where {N}
    return [min(I), max(I)], zeros(N, 2)
end

"""
    linear_map(M::AbstractMatrix, x::Interval)

Concrete linear map of an interval.

### Input

- `M` -- matrix
- `x` -- interval

### Output

Either an interval or a zonotope, depending on the leading dimension (i.e. the
number of rows) of `M`:

- If `size(M, 1) == 1`, the output is an interval obtained by scaling `x` by the
  matrix `M`.
- If `size(M, 1) > 1`, the output is a zonotope whose center is `M * center(x)`
  and it has the single generator, `M * g`, where `g = (high(x)-low(x))/2`.
"""
function linear_map(M::AbstractMatrix, x::Interval)
    @assert size(M, 2) == 1 "a linear map of size $(size(M)) " *
                            "cannot be applied to an interval"
    nout = size(M, 1)
    if nout == 1
        return _linear_map_interval(M, x)
    else
        return _linear_map_zonotope(M, x)
    end
end

function _linear_map_interval(M::AbstractMatrix, x::Interval)
    α = M[1, 1]
    return Interval(α * x.dat)
end

function _linear_map_zonotope(M::AbstractMatrix, x::Interval)
    nout = size(M, 1)
    cx = IntervalArithmetic.mid(x.dat)
    gx = cx - min(x)
    N = promote_type(eltype(M), eltype(x))
    c = Vector{N}(undef, nout)
    gen = Matrix{N}(undef, nout, 1)
    @inbounds for i in 1:nout
        c[i] = M[i, 1] * cx
        gen[i] = M[i, 1] * gx
    end
    return Zonotope(c, gen)
end

"""
    scale(α::Real, x::Interval)

Concrete scaling of an interval.

### Input

- `α` -- scalar
- `x` -- interval

### Output

The interval obtained by applying the numerical scale to the given interval.
"""
function scale(α::Real, x::Interval)
    return Interval(α * x.dat)
end

"""
    center(H::Interval, i::Int)

Return the center along a given dimension of a interval.

### Input

- `x` -- interval
- `i` -- dimension of interest

### Output

The center along a given dimension of the interval.
"""
@inline function center(x::Interval, i::Int)
    @boundscheck i == 1 || throw(ArgumentError("an interval has dimension one, but the index is $i"))
    return IntervalArithmetic.mid(x.dat)
end

"""
    rectify(x::Interval{N}) where {N}

Concrete rectification of an interval.

### Input

- `x` -- interval

### Output

The `Interval` that corresponds to the rectification of `x`.

### Notes

Note that the result is an `Interval` even if the set becomes a singleton (which
is the case if the original interval was nonpositive).
"""
function rectify(x::Interval{N}) where {N}
    if x.dat.lo >= zero(N)
        # interval is already nonnegative
        return x
    else
        # lower end is negative
        return Interval(zero(N), max(x.dat.hi, zero(N)))
    end
end

"""
    diameter(x::Interval, [p]::Real=Inf)

Compute the diameter of an interval, defined as ``\\Vert b - a\\Vert`` in the ``p`-norm, where ``a`` (resp. ``b``) are the minimum (resp. maximum) of the given interval.

### Input

- `x` -- interval
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the diameter.

### Notes

In one dimension all norms are the same.
"""
function diameter(x::Interval, p::Real=Inf)
    return IntervalArithmetic.diam(x.dat)
end

"""
    split(x::Interval, k)

Partition an interval into `k` uniform sub-intervals.

### Input

- `x` -- interval
- `k` -- number of sub-intervals, possibly wrapped in a vector of length 1

### Output

A list of `k` `Interval`s.
"""
function split(x::Interval, k::AbstractVector{Int})
    @assert length(k) == 1 "an interval can only be split along one dimension"
    return split(x, k[1])
end

function split(x::Interval, k::Int)
    @assert k > 0 "can only split into a positive number of intervals"
    return [Interval(x2) for x2 in mince(x.dat, k)]
end

function vertices_list(x::IntervalArithmetic.Interval{N}) where {N}
    a = IntervalArithmetic.inf(x)
    b = IntervalArithmetic.sup(x)
    ST = IntervalArithmetic.SVector{1, N}
    if _isapprox(a, b)
        vlist = [ST(a)]
    else
        vlist = [ST(a), ST(b)]
    end
    return vlist
end

function vertices_list(H::IntervalArithmetic.IntervalBox)
    return vertices_list(convert(Hyperrectangle, H))
end

function chebyshev_center(x::Interval{N}; compute_radius::Bool=false) where {N}
    if compute_radius
        return center(x), zero(N)
    end
    return center(x)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

datadir = joinpath(Pkg.dir("FEMBeam"), "test", first(splitext(basename(@__FILE__))))

"""
    read_mtx(fn)

Read matrix file written using ABAQUS:

*MATRIX GENERATE, STIFFNESS, MASS, LOAD
*MATRIX OUTPUT, STIFFNESS, MASS, LOAD

"""
function read_mtx(fn; dim=error("dofs / node missing: pass argument dim::Int64"))
    I = Int64[]
    J = Int64[]
    V = Float64[]
    open(fn) do fid
        for ln in eachline(fid)
            i,idof,j,jdof,value = map(parse, split(ln, ','))
            if i < 1 || j < 1
                continue
            end
            push!(I, (i-1)*dim+idof)
            push!(J, (j-1)*dim+jdof)
            push!(V, value)
        end
    end
    return sparse(I,J,V)
end

function read_mtx!(A, fn)
    open(fn) do fid
        for ln in eachline(fid)
            i,idof,j,jdof,value = map(parse, split(ln, ','))
            A[i,idof,j,jdof] = value
        end
    end
    return A
end

K = read_mtx(joinpath(datadir, "eb3arxs4_STIF1.mtx"); dim=6)
# K1 = read_mtx(joinpath(datadir, "stif.mtx"); dim=2)
# K2 = zeros(2, 2, 2, 2)
# read_mtx!(K2, joinpath(datadir, "stif.mtx"))

println(K)

ext: jl


text:
cd(dirname(@__FILE__))

using ProgressMeter, JLD, Statistics, MAT, Plots

fileIn = matopen("i140703-001_lfp-spikes.mat")
Data = read(fileIn)
close(fileIn)

events = Data["block"]["segments"][1]["events"][1]

TS_ON = events["times"][findall(x -> x == "65296", events["labels"])] ./ 30000
WS_ON = events["times"][findall(x -> x == "65360", events["labels"])][1:2:end] ./ 30000
CUE_ON = events["times"][findall(x -> x == "65365" || x == "65370", events["labels"])] ./ 30000
CUE_OFF = events["times"][findall(x -> x == "65360", events["labels"])][2:2:end] ./ 30000
GO_ON = events["times"][findall(x -> x == "65369" || x == "65366", events["labels"])] ./ 30000
SR = events["times"][findall(x -> x == "65385" || x == "65382", events["labels"])] ./ 30000
RW_ON = events["times"][findall(x -> x == "65509" || x == "65514", events["labels"])] ./ 30000
GO_RW_OFF = events["times"][findall(x -> x == "65376", events["labels"])] ./ 30000
STOP = events["times"][findall(x -> x == "65312", events["labels"])] ./ 30000
spiketrains = [Data["block"]["segments"][1, 1]["spiketrains"][i]["times"][:] for i ∈ 1:271] / 30000

####################################

Prog = Progress(271 * length(TS_ON))
cutspiketrains = []
for i ∈ 1:271
    for j ∈ 1:length(TS_ON)
        indices = findall(x -> x > [[0.0]; TS_ON][j] && x < [[0.0]; TS_ON][j+1], spiketrains[i])
        push!(cutspiketrains, [(i, TS_ON[j] - spiketrains[i][n]) for n ∈ indices]...)
        next!(Prog)
    end
end

save("cutspiketrains.jld", "cutspiketrains", cutspiketrains)

cutspiketrains

####################################

cut_TS_ON = []

for j ∈ 1:271
    push!(cut_TS_ON, (j, 0))
end

cut_TS_ON

####################################

cut_WS_ON = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], WS_ON)
    for j ∈ 1:271
        push!(cut_WS_ON, [(j, TS_ON[i] - WS_ON[n]) for n ∈ indices]...)
    end
end

cut_WS_ON

####################################

cut_CUE_ON = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], CUE_ON)
    for j ∈ 1:271
        push!(cut_CUE_ON, [(j, TS_ON[i] - CUE_ON[n]) for n ∈ indices]...)
    end
end

cut_CUE_ON

####################################

cut_CUE_OFF = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], CUE_OFF)
    for j ∈ 1:271
        push!(cut_CUE_OFF, [(j, TS_ON[i] - CUE_OFF[n]) for n ∈ indices]...)
    end
end

cut_CUE_OFF

####################################

cut_GO_ON = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], GO_ON)
    for j ∈ 1:271
        push!(cut_GO_ON, [(j, TS_ON[i] - GO_ON[n]) for n ∈ indices]...)
    end
end

cut_GO_ON

####################################

cut_CUE_OFF = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], CUE_OFF)
    for j ∈ 1:271
        push!(cut_CUE_OFF, [(j, TS_ON[i] - CUE_OFF[n]) for n ∈ indices]...)
    end
end

cut_CUE_OFF

####################################

cut_SR = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], SR)
    for j ∈ 1:271
        push!(cut_SR, [(j, TS_ON[i] - SR[n]) for n ∈ indices]...)
    end
end

cut_SR

####################################

cut_RW_ON = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], RW_ON)
    for j ∈ 1:271
        push!(cut_RW_ON, [(j, TS_ON[i] - RW_ON[n]) for n ∈ indices]...)
    end
end

cut_RW_ON

####################################

cut_GO_RW_OFF = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], GO_RW_OFF)
    for j ∈ 1:271
        push!(cut_GO_RW_OFF, [(j, TS_ON[i] - GO_RW_OFF[n]) for n ∈ indices]...)
    end
end

cut_GO_RW_OFF

####################################

cut_STOP = []

for i ∈ 1:length(TS_ON)
    indices = findall(x -> x > [[0.0]; TS_ON][i] && x < [[0.0]; TS_ON][i+1], STOP)
    for j ∈ 1:271
        push!(cut_STOP, [(j, TS_ON[i] - STOP[n]) for n ∈ indices]...)
    end
end

cut_STOP

####################################
[[1]; [2, 3]]

scatter(getindex.(cutspiketrains, 2), getindex.(cutspiketrains, 1), marker = :vline, c = :black, alpha = 0.01)

scatter!(getindex.(cut_TS_ON, 2), getindex.(cut_TS_ON, 1), marker = :vline, c = :red, ms = 5)

scatter!(getindex.(cut_WS_ON, 2), getindex.(cut_WS_ON, 1), marker = :vline, c = :green)

scatter!(getindex.(cut_CUE_OFF, 2), getindex.(cut_CUE_OFF, 1), marker = :vline, c = :blue)

scatter!(getindex.(cut_GO_ON, 2), getindex.(cut_GO_ON, 1), marker = :vline, c = :purple)

scatter!(getindex.(cut_SR, 2), getindex.(cut_SR, 1), marker = :vline, c = :brown)

scatter!(getindex.(cut_RW_ON, 2), getindex.(cut_RW_ON, 1), marker = :vline, c = :orange)

scatter!(getindex.(cut_GO_RW_OFF, 2), getindex.(cut_GO_RW_OFF, 1), marker = :vline, c = :yellow)

scatter!(getindex.(cut_STOP, 2), getindex.(cut_STOP, 1), marker = :vline, c = :black)

max_stars_repo_name: ali-mahani/neuro-project

max_issues_repo_name: ali-mahani/neuro-project

max_forks_repo_name: ali-mahani/neuro-project


text:
@testset "new macro notation" begin
    A = [i^2 for i in 1:12]

    # index ranges
    @test reshape(A,3,4) == @cast C[i,j] := A[(i,j)]  i in 1:3
    @test_logs min_level=Logging.Warn @reduce D[i] := sum(j:4) A[i⊗j]  # no warning for this

    # underscores
    @test A' == @cast _[_,i] := A[i]
    @test (@cast _[_,i] := A[i]) isa Matrix  # transmute is reshape here

    # ... use in recursion, reads from the _ just written to:
    A = rand(4);
    B = randn(4,4);
    R1 = @reduce sum(i) A[i] * log( @reduce _[i] := sum(j) A[j] * exp(B[i,j]) ) # new
    R2 = @reduce sum(i) A[i] * log( @reduce [i] := sum(j) A[j] * exp(B[i,j]) ) # old, with warning
    @reduce inner[i] := sum(j) A[j] * exp(B[i,j])
    S = @reduce sum(i) A[i] * log(inner[i])
    @test S == R1 == R2

    # lazier
    Cs = [fill(i^2,7) for i in 1:12]
    using LazyStack
    @test (@cast _[i,o] := Cs[o][i]) isa LazyStack.Stacked
    @test (@cast _[i,o] |= Cs[o][i]) isa Matrix  # collects afterwards
    @test (@cast _[i,o] := Cs[o][i] lazy=false) isa Matrix  # uses stack_iter

    using LinearAlgebra, TransmuteDims
    @test (@cast _[o,i] := Cs[o][i]) isa Transpose{Int, <:LazyStack.Stacked}
    @test (@cast _[o,_,i] := Cs[o][i]) isa TransmutedDimsArray

    # broadcasting
    using Strided, LazyArrays
    @test A .+ 1 == @cast @strided B[i] := A[i] + 1
    # @test (@cast @lazy B[i] := A[i] + 1) isa BroadcastVector  # fails on 1.3?
end

@testset "ternary operator" begin
    CNT = Ref(0)
    add(x) = begin CNT[] += 1; x+1 end
    xs = collect(1:12)

    @cast ys[i] := iseven(xs[i]) ? add(xs[i]) : 0
    @test CNT[] == 6
    @test ys[1:6] == [0, 3, 0, 5, 0, 7]

    @cast zs[i,j] := iseven(xs[i⊗j]+1) ? add(ys[j⊗i] / 2) : ((1:3)[i] + 100)
    @test CNT[] == 12
    @test size(zs) == (3, 4)

    rs = randn(10,10)
    @reduce ps[i,_] := sum(j) if rs[i,j] > 0
        sqrt(rs[i,j])
    else
        0 * cbrt(rs[i,j])
    end
    @test ps ≈ sum(x -> x>0 ? sqrt(x) : 0.0, rs, dims=2)
end

@testset "naked indices" begin
    @cast A[i,j] := i+10j  (i in 1:3, j in 1:4)  # explicit size
    @test A == (1:3) .+ (10:10:40)'

    j = 99
    @cast A[i,j] = div(i, j)  # inferred from LHS
    @test A[3,1] == 3

    @reduce B[i] := sum(j) A[i,j] + i + $j  # interpolate j
    @test B[3] > 400

    @test [1, 4, 9] .- im == @cast _[i'] := i'^2 + im'  (i' in 1:3)  # primes

    @cast C[i,j,k] := 0 * A[i,(j,k)] + j  (k in 1:2)  # used to infer sz_j = (:)
    @test all(==(2), C[:,2,:])

    @reduce D[k] := sum(i) B[i]/k (k in 1:4)  # no indexing by k on RHS
    @test D ≈ vec(sum(B ./ (1:4)', dims=1))

    @reduce E := sum(i,k) i/k (i in 1:2, k in 1:4)  # no indexing on RHS of reduction
    @test E ≈ sum((1:2) ./ (1:4)')
end

@testset "offset handling" begin
    using OffsetArrays
    @cast A[i,j] := i+10j  (i in 0:1, j in 7:15)
    @test axes(A) == (0:1, 7:15)
    @test A[1,7] == 71

    # reshape
    @cast B[(i,k),j] := A[i,(j,k)]  k in 1:3
    @test axes(B) === (Base.OneTo(6), Base.OneTo(3))

    # reduction
    @reduce C[_,j] := sum(i) A[i,j]
    @test axes(C) == (1:1, 7:15)
    @test extrema(C) == (141, 301)  # was reading out of bounds, TransmuteDims bug

    # slicing
    @test axes(@cast _[j] := A[:,j]) == (7:15,)
    @test axes(@cast _[j][i] := A[i,j]) == (7:15,)
    @test axes(first(@cast _[j] := A[:,j])) == (0:1,)
    @test axes(first(@cast _[j][i] := A[i,j])) == (0:1,)

    using StaticArrays
    @test_throws Exception @cast _[j] := A{:,j}  # similar error to reinterpret(reshape, SVector{2,Int}, A)
    @cast D[i,j] := i+10j  (i in 1:3, j in 7:15)
    @test axes(@cast _[j] := D{:,j}) == (7:15,)
    @test first(@cast _[j] := D{:,j}) === SVector(71, 72, 73)
    @test first(@cast _[j] := D{:3,j}) === SVector(71, 72, 73)
    @test_throws Exception @cast _[j] := D{:2,j}  # wrong size
end

@testset "tuples" begin
    # stack acts on vectors of tuples
    x = rand(3, 5)
    @cast vi[j,k] := findmax(x[:, k])[j]
    @test vi[1,1] == maximum(x[:,1])
    @test vi[2,2] == argmax(x[:,2])

    # and on tuples of vectors!
    vs = (ones(4), rand(4), zeros(4))
    @cast xs[i,j] := vs[j][i]
    @test xs[4,3] == 0
    @test (@cast _[j,i] := vs[j][i] lazy=false) isa Matrix

    # construction from vectors
    y = rand(Int8, 5,4,3)
    @cast z[j,i] := Tuple(y[i,j,:])
    @test z[1,2] == Tuple(y[2,1,:])

    # construction with if else
    @cast tri[i,j,c] := tuple(c==1 ? i : 0, c==2 ? j : 0)  (i in 1:3, j in 1:5, c in 1:2)
    @test all(x -> x[2]==0, tri[:,:,1])
    @test (x->x[2]).(tri[:,:,2]) == [1,1,1] .* (1:5)'
    @cast tri[i,j,c] = (c==1 ? i^2 : 0, ifelse(c==2, j^2, 0))  # no "tuple", sizes inferred, ifelse
    @test tri[3,2,1] == (3^2, 0)
end

@testset "splats" begin
    y = rand(Int8, 5,4,3)

    @cast z′[j,i] := tuple(y[i,j,:]...)
    @test z′[3,4] == Tuple(y[4,3,:])

    # with other arguments, they have to come first at the moment:
    @cast z2[j,i] := tuple(i, j, y[i,j,:]...)
    @test z2[4,5] == (5, 4, y[5,4,:]...)

    struct Quad x; y; z; t; end
    @cast z3[i,j] := Quad(y[i,:,j]...)
    @test z3[2,3] == Quad(y[2,:,3]...)

    # stupid operation:
    @test_broken y == @cast _[i,j,k] := (y[i,:,k]...,)[j]  # (,) does not work with ...
    @test y == @cast _[i,j,k] := tuple(y[i,:,k]...)[j]
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function declared_floor(x::Float64)
    return ccall("llvm.floor.f64", llvmcall, Float64, (Float64,), x)
end
@test declared_floor(4.2) == 4.0
ir = sprint(io->code_llvm(io, declared_floor, Tuple{Float64}))
@test contains(ir, "call double @llvm.floor.f64") # should be inlined

function doubly_declared_floor(x::Float64)
    a = ccall("llvm.floor.f64", llvmcall, Float64, (Float64,), x)
    b = ccall("llvm.floor.f64", llvmcall, Float64, (Float64,), x + 1) - 1
    return a + b
end
@test doubly_declared_floor(4.2) == 8.0

function doubly_declared2_trunc(x::Float64)
    a = ccall("llvm.trunc.f64", llvmcall, Float64, (Float64,), x)
    b = ccall("llvm.trunc.f64", llvmcall, Float64, (Float64,), x + 1) - 1
    return a + b
end
@test doubly_declared2_trunc(4.2) == 8.0

# Test for single line
function declared_ceil(x::Float64)
    return ccall("llvm.ceil.f64", llvmcall, Float64, (Float64,), x)
end
@test declared_ceil(4.2) == 5.0

# Test for multiple lines
function ceilfloor(x::Float64)
    a = ccall("llvm.ceil.f64", llvmcall, Float64, (Float64,), x)
    b = ccall("llvm.floor.f64", llvmcall, Float64, (Float64,), a)
    return b
end
@test ceilfloor(7.4) == 8.0


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: July 2020

### Model
model, model0, dm = setup_model(tti, viscoacoustic, nlayer)
q, srcGeometry, recGeometry, f0 = setup_geom(model)
dt = srcGeometry.dt[1]

@testset "Gradient options test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
        ##################################ISIC########################################################
        println("Testing isic")
        @timeit TIMEROUTPUT "ISIC" begin
                opt = Options(sum_padding=true, free_surface=fs, isic=true, f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)
                @test norm(J*(0f0.*dm)) == 0

                y0 = F*q
                y_hat = J*dm
                x_hat1 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat1)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test isapprox(c, d, rtol=5f-2)
                @test !isnan(norm(y0))
                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat1))
        end

        ##################################checkpointing###############################################
        println("Testing checkpointing")
        @timeit TIMEROUTPUT "Checkpointing" begin
                opt = Options(sum_padding=true, free_surface=fs, optimal_checkpointing=true, f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)

                y_hat = J*dm
                x_hat2 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat2)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test isapprox(c, d, rtol=1f-2)

                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat2))
        end

        ##################################DFT#########################################################
        println("Testing DFT")
        @timeit TIMEROUTPUT "DFT" begin
                opt = Options(sum_padding=true, free_surface=fs, frequencies=[2.5, 4.5], f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)
                @test norm(J*(0f0.*dm)) == 0

                y_hat = J*dm
                x_hat3 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat3)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat3))
        end

        ################################## DFT time subsampled#########################################
        println("Testing subsampled in time DFT")
        @timeit TIMEROUTPUT "Subsampled DFT" begin
                opt = Options(sum_padding=true, free_surface=fs, frequencies=[2.5, 4.5], dft_subsampling_factor=4, f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)
                @test norm(J*(0f0.*dm)) == 0

                y_hat = J*dm
                x_hat3 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat3)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat3))
        end

        ##################################subsampling#################################################
        println("Testing subsampling")
        @timeit TIMEROUTPUT "Subsampling" begin
                opt = Options(sum_padding=true, free_surface=fs, subsampling_factor=4, f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)

                y_hat = J*dm
                x_hat4 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat4)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test isapprox(c, d, rtol=1f-2)
                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat4))
        end
        ##################################ISIC + DFT #########################################################
        println("Testing isic+dft")
        @timeit TIMEROUTPUT "ISIC+DFT" begin
                opt = Options(sum_padding=true, free_surface=fs, isic=true, frequencies=[2.5, 4.5], f0=f0)
                F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

                # Linearized modeling
                J = judiJacobian(F, q)
                @test norm(J*(0f0.*dm)) == 0

                y_hat = J*dm
                x_hat5 = adjoint(J)*y0

                c = dot(y0, y_hat)
                d = dot(dm, x_hat5)
                @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, c/d - 1)
                @test !isnan(norm(y_hat))
                @test !isnan(norm(x_hat5))
        end
end


text:
###############################################################################
### ShapefileProvider
#
# ShapefileProvider's must have the following fields:
#
# basename::AbstractString -- basename of the shapefile; also used as cache
#                             directory
# url::AbstractString -- canonical url of the shapefile
#
# polygons:Array -- will store the polygons extracted from the .shp
#
# records::Dict -- will store the records extracted from the .dbf (if available)

function polygons(provider::ShapefileProvider)
  return provider.polygons
end

function records(provider::ShapefileProvider)
  return provider.records
end

function names(provider::ShapefileProvider)
  return collect(keys((records(provider))))
end

# I don't really like this Pythonic naming scheme, but since I specifically
# don't want to export the Provider pseudo-constructors (eg. Provider.STATE())
# and Julia doesn't provide a straightforward way to completely hide functions,
# this seems like a decent compromise.

function _localize(provider::ShapefileProvider)
  _iscached(provider) && return
  cachedir = joinpath(Pkg.dir("ChoroplethMaps"), "cache")
  zipname = joinpath(cachedir, string(provider.basename, ".zip"))
  mkpath(dirname(zipname))

  # Maybe use Requests.jl here instead?
  download(provider.url, zipname)

  _unzip(zipname, joinpath(cachedir, provider.basename))
end

function _iscached(provider::ShapefileProvider)
  if isdir(_getcachedir(provider))
    return isfile(_getshpfilepath(provider))
  end
  return false
end

function _getcachedir(provider::ShapefileProvider)
  return joinpath(Pkg.dir("ChoroplethMaps"), "cache", provider.basename)
end

function _getshpfilepath(provider::ShapefileProvider)
  return joinpath(_getcachedir(provider), "$(provider.basename).shp")
end

function _getdbffilepath(provider::ShapefileProvider)
  return joinpath(_getcachedir(provider), "$(provider.basename).dbf")
end

function _unzip(cfile, dir)
  reader = ZipFile.Reader(cfile)
  for file in reader.files
    ucname = joinpath(dir, file.name)
    mkpath(dirname(ucname))
    ucfile = open(ucname, "w")
    write(ucfile, readall(file))
    close(ucfile)
  end
  close(reader)
end

function _parseshapes(provider::ShapefileProvider)
  shpf = open(_getshpfilepath(provider)) do fd
    read(fd, Shapefile.Handle)
  end
  # ChoroplethMaps module expects polygon data to be a
  # Vector of Vector of Vector of Float64
  polygons = Vector[]
  for shape in shpf.shapes
    polygon = Vector[]
    for pt in shape.points
      push!(polygon, [pt.x, pt.y])
    end
    push!(polygons, polygon)
  end
  return polygons
end

function _parserecords(provider::ShapefileProvider)
  dbffile = _getdbffilepath(provider)
  if !isfile(dbffile)
    warn("No .dbf associated with shapefile $(provider.basename)!")
    return Dict()
  end
  return readdbf(dbffile)
end

###############################################################################


text:
using NiLang.AD
export treeverse, treeverse_solve

struct TreeverseAction
    action::Symbol
    τ::Int
    δ::Int
    step::Int
    depth::Int
end

struct TreeverseLog
    actions::Vector{TreeverseAction}
    depth::Base.RefValue{Int}
    peak_mem::Base.RefValue{Int}  # should be `n*(k-1)+2`
end
TreeverseLog() = TreeverseLog(TreeverseAction[], Ref(0), Ref(0))
Base.push!(tlog::TreeverseLog, args...) = push!(tlog.actions, TreeverseAction(args..., tlog.depth[]))

Base.show(io::IO, ::MIME"text/plain", logger::TreeverseLog) = Base.show(io, logger)
function Base.show(io::IO, logger::TreeverseLog)
    print(io, """Treeverse log
| peak memory usage = $(logger.peak_mem[])
| number of function calls = $(count(x->x.action==:call, logger.actions))
| number of gradient calls = $(count(x->x.action==:grad, logger.actions))
| number of stack push/pop = $(count(x->x.action==:store, logger.actions))/$(count(x->x.action==:fetch, logger.actions))""")
end

function binomial_fit(N::Int, δ::Int)
    τ = 1
    while N > binomial(τ+δ, τ)
        τ += 1
    end
    return τ
end

function mid(δ, τ, σ, ϕ)
    κ = ceil(Int, (δ*σ + τ*ϕ)/(τ+δ))
    if κ >= ϕ && δ > 0
        κ = max(σ+1, ϕ-1)
    end
    return κ
end

"""
    treeverse(f, gf, s; δ, N, τ=binomial_fit(N,δ), f_inplace=true, logger = TreeverseLog())

Treeverse algorithm for back-propagating a program memory efficiently.

Positional arguments
* `f`, the step function that ``s_{i+1} = f(s_i)``,
* `gf`, the single step gradient function that ``g_i = gf(s_i, g_{i+1})``
   !!! When ``g_{i+1}`` is `nothing`, it should return the gradient of the loss function,
* `s`, the initial state ``s_0``,

Keyword arguments
* `δ`, the number of checkpoints,
* `N`, the number of time steps,
* `τ`, the number of sweeps, it is chosen as the smallest integer that `binomial(τ+δ, τ) >= N` by default,
* `f_inplace = false`, whether `f` is inplace,
* `logger = TreeverseLog()`, the logger.

Ref: https://www.tandfonline.com/doi/abs/10.1080/10556789208805505
"""
function treeverse(f, gf, s::T; δ, N, τ=binomial_fit(N,δ), f_inplace=false, logger = TreeverseLog()) where T
    state = Dict(0=>s)
    if N > binomial(τ+δ, τ)
        error("please input a larger `τ` and `δ` so that `binomial(τ+δ, τ) >= N`!")
    end
    g = treeverse!(f, gf, state, nothing, δ, τ, 0, 0, N, logger, f_inplace)
    return g
end

function treeverse!(f, gf, state::Dict{Int,T}, g, δ, τ, β, σ, ϕ, logger, f_inplace) where T
    logger.depth[] += 1
    # cache sσ
    if σ > β
        δ -= 1
        s = state[β]
        for j=β:σ-1
            s = NiLang.getf(f, j)(s)
            push!(logger, :call, τ, δ, j)
        end
        store_state!(state,σ, f_inplace ? copy(s) : s)
        s = nothing
        push!(logger, :store, τ, δ, σ)
        logger.peak_mem[] = max(logger.peak_mem[], length(state))
    elseif σ < β
        error("treeverse fails! σ < β")
    end

    κ = mid(δ, τ, σ, ϕ)
    while τ>0 && κ < ϕ
        g = treeverse!(f, gf, state, g, δ, τ, σ, κ, ϕ, logger, f_inplace)
        τ -= 1
        ϕ = κ
        κ = mid(δ, τ, σ, ϕ)
    end

    g = NiLang.getf(gf, σ)(state[σ], g)
    push!(logger, :grad, τ, δ, σ)
    if σ>β
        # remove state[σ]
        delete_state!(state, σ)
        push!(logger, :fetch, τ, δ, σ)
    end
    return g
end

@inline function store_state!(state::Dict, i::Int, x)
    state[i] = x
end

@inline function delete_state!(state::Dict, i::Int)
    pop!(state, i)
end

function treeverse_step(s, param, src, srcv, c)
    unext, φ, ψ = zero(s.u), copy(s.φ), copy(s.ψ)
    ReversibleSeismic.one_step!(param, unext, s.u, s.upre, φ, ψ, param.Σx, param.Σy, c)
    s2 = SeismicState(s.u, unext, φ, ψ, Ref(s.step[]+1))
    s2.u[src...] += srcv[s2.step[]]*param.DELTAT^2
    return s2
end

function treeverse_grad(x, g, param, src, srcv, gsrcv, c, gc)
    y = treeverse_step(x, param, src, srcv, c)
    gt = SeismicState([GVar(getfield(y, field), getfield(g, field)) for field in fieldnames(SeismicState)[1:end-1]]..., Ref(y.step[]))
    _, gs, _, _, gv, gc2 = (~bennett_step!)(gt, GVar(x), param, src, GVar(srcv, gsrcv), GVar(c, gc))
    (grad(gs), grad(gv), grad(gc2))
end

"""
    treeverse_solve(s0; param, src, srcv, c, δ=20, logger=TreeverseLog())

* `s0` is the initial state,
"""
function treeverse_solve(s0, gnf; param, src, srcv, c, δ=20, logger=TreeverseLog())
    f = x->treeverse_step(x, param, src, srcv, c)
    res = []
    function gf(x, g)
        if g === nothing
            y = f(x)
            push!(res, y)
            g = gnf(y)
        end
        treeverse_grad(x, g[1], param, src, srcv, g[2], c, g[3])
    end
    g = treeverse(f, gf,
        copy(s0); δ=δ, N=param.NSTEP-1, f_inplace=false, logger=logger)
    return res[], g
end


text:
include("task_func.jl")

function hel(r, side, lockSide)
    count = 0
    flag = 0
    while(!isborder(r, lockSide))
        listOfShags = []
        while(!isborder(r, side) || isRect(r, side))
            if(isRect(r, side))
                append!(listOfShags, aroundBorder!(r, side, save = true))
            else
                push!(listOfShags, goToSideWithShag!(r, side, 1, save = true))
            end
            if(flag == 0 && isborder(r, lockSide))
                flag = 1
                count += 1
            elseif(flag == 1 && !isborder(r, lockSide))
                flag = 0
            end
        end
        goToStart!(r, listOfShags)
        move!(r, lockSide)
    end
    return count
end

function draw(r)
    while !isborder(r, Nord) || !isborder(r, West)
        goToSide!(r, Nord)
        goToSide!(r, West)
    end
    answer = 0
    answer += hel(r, Ost, Sud)
    answer += hel(r, Nord, Ost)
    print(answer)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
module GingerGUI

include("main.jl")

end # module

max_stars_repo_name: pacific-zitter/GingerGUI

max_issues_repo_name: pacific-zitter/GingerGUI

max_forks_repo_name: pacific-zitter/GingerGUI


text:
"""
    data = get_permissions(db::Database)

Fetch all current permissions. Note: this is Cloudant-specific.

[API reference](https://docs.cloudant.com/authorization.html#viewing-permissions)
"""
function get_permissions(db::Database)
  relax(get, endpoint(db.client.url, "_api/v2/db/$(db.name)/_security"); cookies=db.client.cookies)
end

"""
    result = set_permissions(db::Database, current::Dict=Dict{AbstractString, Any}(); key="", roles=[])

Modify permissions. Note: this is Cloudant-specific.

[API reference](https://docs.cloudant.com/authorization.html#modifying-permissions)
"""
function set_permissions(db::Database, current::Dict=Dict{AbstractString, Any}(); key="", roles=[])
  if key != "" && length(roles) == 0
    error("No roles provided for key")
  end
  if !haskey(current, "cloudant") && key == ""
    error("Expected a non-empty permissions dict or a non-empty key")
  end

  if !haskey(current, "cloudant")
    current["cloudant"] = Dict(key => roles)
  else
    current["cloudant"][key] = roles
  end

  relax(put, endpoint(db.client.url, "_api/v2/db/$(db.name)/_security"); json=current, cookies=db.client.cookies)
end

"""
    data = make_api_key(client::Client)

Generate a new API key. Note: this is Cloudant-specific.

Note also that API keys take a long time to propagate around a cluster. It's unsafe to rely
on a newly created key to be immediately available. The reason for this is that Cloudant keeps
its auth-related documents centrally, and replicate out to all clusters.

[API reference](https://docs.cloudant.com/authorization.html#creating-api-keys)
"""
function make_api_key(client::Client)
  relax(post, endpoint(client.url, "_api/v2/api_keys"); cookies=client.cookies)
end

"""
    result = delete_api_key(db::Database, key::AbstractString)

Remove an existing API key. Note: this is Cloudant-specific. This is implemented via `set_permissions()`.

[API reference](https://docs.cloudant.com/authorization.html#deleting-api-keys)
"""
function delete_api_key(db::Database, key::AbstractString)
  current = get_permissions(db)
  if haskey(current, "cloudant") && haskey(current["cloudant"], key)
    delete!(current["cloudant"], key)
    set_permissions(db, current)
    return true
  end
  false
end


text:
using HorizonSideRobots

function mark_field_with_horizontal_obstacles(r)
    back_track = go_to_bottom_left_corner_return_path(r)

    while isborder(r, Nord) == false
        mark_n_back(r, Ost)
        putmarker!(r)
        move!(r, Nord)
    end

    mark_n_back(r, Ost)
    putmarker!(r)
    while isborder(r, Sud) == false
        move!(r, Sud)
    end

    back_by_path_o_sides(r, back_track)
end

function go_to_bottom_left_corner_return_path(r::Robot)::Array
    path = []
    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        if isborder(r, Sud) == false
            move!(r, Sud)
            push!(path, Nord)
        end
        if isborder(r, West) == false
            move!(r, West)
            push!(path, Ost)
        end
    end
    return path
end

function back_by_path_o_sides(r, path::Array)
    n = length(path)
    while n > 0
        move!(r, path[n])
        n = n - 1
    end
end

function mark_n_back(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
        putmarker!(r)
    end
    while ismarker(r)==true
        move!(r, HorizonSide(mod(Int(side) + 2, 4))) # В противоположную сторону
    end
end


text:
function constraint_model_voltage_active(pm::_PM.AbstractDCPModel, n::Int, c::Int)
end

""
function variable_voltage_active(pm::_PM.AbstractDCPModel; kwargs...)
    _PM.variable_voltage_angle(pm; kwargs...)
end

"no vm values to turn off"
function constraint_bus_active(pm::_PM.AbstractDCPModel, i::Int; nw::Int=nw_id_default, kwargs...)
end


text:
# maximum(p::_Poly.Polynomial) = Base.maximum(p)
# minimum(p::_Poly.Polynomial) = Base.minimum(p)
pdf(p::_Poly.Polynomial) = error("`pdf` of a polynomial is not definied because it is expected that the polynomial used is an approximation of the `logpdf` instead!")
logpdf(p::_Poly.Polynomial, x::Real) = p(x)
gradlogpdf(p::_Poly.Polynomial, x::Real) = _Poly.derivative(p, 1)(x)
heslogpdf(p::_Poly.Polynomial, x::Real) = _Poly.derivative(p, 2)(x)


text:
# ---
# title: 956. Tallest Billboard
# id: problem956
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/tallest-billboard/description/>
# hidden: true
# ---
#
# You are installing a billboard and want it to have the largest height.  The
# billboard will have two steel supports, one on each side.  Each steel support
# must be an equal height.
#
# You have a collection of `rods` which can be welded together.  For example, if
# you have rods of lengths 1, 2, and 3, you can weld them together to make a
# support of length 6.
#
# Return the largest possible height of your billboard installation.  If you
# cannot support the billboard, return 0.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,6]
#     Output: 6
#     Explanation: We have two disjoint subsets {1,2,3} and {6}, which have the same sum = 6.
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,4,5,6]
#     Output: 10
#     Explanation: We have two disjoint subsets {2,3,5} and {4,6}, which have the same sum = 10.
#
#
# **Example 3:**
#
#
#
#     Input: [1,2]
#     Output: 0
#     Explanation: The billboard cannot be supported, so we return 0.
#
#
#
#
# **Note:**
#
#   1. `0 <= rods.length <= 20`
#   2. `1 <= rods[i] <= 1000`
#   3. `The sum of rods is at most 5000.`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using GadJet, Test, DelimitedFiles, Unitful, UnitfulAstro

@testset "GadJet" begin

    @testset "Read Snapshot" begin

        snap_file = joinpath(dirname(@__FILE__), "snap_050")

        @test_nowarn read_snap(snap_file)
        @test_nowarn read_header(snap_file)
        @test_nowarn read_info(snap_file)

        d = read_snap(snap_file, "POS", 0)

        ideal_file = joinpath(dirname(@__FILE__), "pos.dat")
        d_ideal = Float32.(readdlm(ideal_file))

        @test d == d_ideal
    end

    @testset "Write Snapshot" begin

        # read in reference file
        ref_file = joinpath(dirname(@__FILE__), "snap_050")
        head = head_to_obj(ref_file)
        x = read_snap(ref_file, "POS", 0)

        # specify output file for testing
        output_file = joinpath(dirname(@__FILE__), "write_test.dat")
        f = open(output_file, "w")

        @test_nowarn write_header(f, head)
        @test_nowarn write_block(f, x, "POS")
        @test_nowarn write_block(f, x, "", snap_format=1)

        # @test_warn "Please specify blockname!" write_block(f, x, "")

        close(f)
    end

    @testset "Unit Conversion" begin

        GU = GadgetPhysicalUnits()

        @test GU.t_s ≈ 3.085678e16u"s"
        @test GU.E_cgs ≈ 1.989e53u"erg"

        d = strip_unit(1.0u"g")
        @test d == 1.0

        GU = GadgetPhysical()

        @test GU.t_s ≈ 3.085678e16
        @test GU.E_cgs ≈ 1.989e53

    end

    # @testset "Riemann Parameters" begin
    #
    #     @test_warn "Both Ul and Pl are zero!" RiemannParameters()
    # end

    @testset "Riemann Sod-Shock" begin

        par = RiemannParameters(Ul=100.0, Mach=10.0, t=1.5)

        sol = solve([86.0], par)

        @test sol.rho3 ≈ 0.37667409437005994
        @test sol.rho4 ≈ 0.4854368932038835
        @test sol.P34  ≈ 13.097397114653086

    end

    @testset "Riemann CR-Shock" begin

        par = RiemannParameters(Pl=63.400, Pr=0.1, Mach=10.0, dsa_model=4, t=1.5)

        sol = solve([86.0], par)

        @test sol.rho3    ≈ 0.3699882303652922
        @test sol.rho4    ≈ 0.5925766595991485
        @test sol.P34_tot ≈ 12.089335761741005
        @test sol.P4_cr   ≈ 3.583255961157783

    end

    @testset "DSA models" begin

        #@test_warn "Invalid DSA model selection!" RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=5, t = 1.5)

        # KR07
        par = RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=0, t = 1.5)
        @test par.acc_function(5.0) ≈ 0.25185919999999995
        # KR13
        par = RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=1, t = 1.5)
        @test par.acc_function(5.0) ≈ 0.09999999999999998
        # Ryu+19
        par = RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=2, t = 1.5)
        @test par.acc_function(5.0) ≈ 0.017286554080677037
        # CS14
        par = RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=3, t = 1.5)
        @test par.acc_function(5.0) ≈ 0.04999999999999999
        # Pfrommer+16
        par = RiemannParameters( Ul = 100.0, Ur = 0.1, dsa_model=4, t = 1.5)
        @test par.acc_function(5.0) == 0.5
    end

    @testset "Synchrotron Kernel" begin

        @test synchrotron_kernel(  1.0) ≈ 0.6514228153553652
        @test synchrotron_kernel(  3.5) ≈ 0.08268719536694746
        @test synchrotron_kernel(  4.5) ≈ 0.03357441971502366
        @test synchrotron_kernel( 10.0) ≈ 0.00019223826430086885
        @test synchrotron_kernel( 50.0) ≈ 1.734785203976593e-21
        @test synchrotron_kernel(100.0) ≈ 4.697593665922202e-43

    end

    @testset "SPH Kernels" begin

        k = Cubic()
        # < 0.5
        d = kernel_value_2D(k, 0.4, 0.5)
        @test d ≈ 0.26992678348385446
        d = kernel_value_3D(k, 0.4, 0.5)
        @test d ≈ 0.13496339174192723
        # < 1.0
        d = kernel_value_2D(k, 0.5, 0.5)
        @test d ≈ 0.15915494309189535
        d = kernel_value_3D(k, 0.5, 0.5)
        @test d ≈ 0.07957747154594767
        # > 1.0
        d = kernel_value_2D(k, 1.5, 0.5)
        @test d == 0.0
        d = kernel_value_3D(k, 1.5, 0.5)
        @test d == 0.0

        k = Quintic()
        # < 1/3
        d = kernel_value_2D(k, 0.3, 0.5)
        @test d ≈ 0.32700352517410614
        d = kernel_value_3D(k, 0.3, 0.5)
        @test d ≈ 0.2791208661307549
        # 2/3
        d = kernel_value_2D(k, 0.5, 0.5)
        @test d ≈ 0.07767855829318415
        d = kernel_value_3D(k, 0.5, 0.5)
        @test d ≈ 0.06630419797168217
        # < 1.0
        d = kernel_value_2D(k, 0.8, 0.5)
        @test d ≈ 0.0008155658657050455
        d = kernel_value_3D(k, 0.8, 0.5)
        @test d ≈ 0.0006961437210839494
        # > 1.0
        d = kernel_value_2D(k, 1.5, 0.5)
        @test d == 0.0
        d = kernel_value_3D(k, 1.5, 0.5)
        @test d == 0.0

        k = WendlandC4()
        # < 1.0
        d = kernel_value_2D(k, 0.5, 0.5)
        @test d ≈ 0.30960610023345264
        d = kernel_value_3D(k, 0.5, 0.5)
        @test d ≈ 0.26606774238812336
        # > 1.0
        d = kernel_value_2D(k, 1.5, 0.5)
        @test d == 0.0
        d = kernel_value_3D(k, 1.5, 0.5)
        @test d == 0.0

        k = WendlandC6()
        # < 1.0
        d = kernel_value_2D(k, 0.5, 0.5)
        @test d ≈ 0.052822211162893276
        d = kernel_value_3D(k, 0.5, 0.5)
        @test d ≈ 0.05055250677698771
        # > 1.0
        d = kernel_value_2D(k, 1.5, 0.5)
        @test d == 0.0
        d = kernel_value_3D(k, 1.5, 0.5)
        @test d == 0.0

    end

    @testset "SPH mappingParameters" begin

        # @test_warn "Giving a center position requires extent in x, y and z direction." mappingParameters()
        #
        # @test_warn "Please specify pixelSideLenght or number of pixels!" mappingParameters(center=[0.0, 0.0, 0.0],
        #                                                                                xlim = [-1.0, 1.0],
        #                                                                                ylim = [-1.0, 1.0],
        #                                                                                zlim = [-1.0, 1.0])

        @test_nowarn mappingParameters(center=[0.0, 0.0, 0.0],
                                       x_lim = [-1.0, 1.0],
                                       y_lim = [-1.0, 1.0],
                                       z_lim = [-1.0, 1.0],
                                       Npixels=100)

        @test_nowarn mappingParameters(center=[0.0, 0.0, 0.0],
                                       x_lim = [-1.0, 1.0],
                                       y_lim = [-1.0, 1.0],
                                       z_lim = [-1.0, 1.0],
                                       pixelSideLength=0.2)
    end

    @testset "glimpse" begin

        snap_file = joinpath(dirname(@__FILE__), "snap_050")
        d = glimpse(snap_file, "RHO", [3.0, 3.0, 3.0], 6.0, 6.0, 6.0,
                    parallel=false, verbose=false)

        ideal_file = joinpath(dirname(@__FILE__), "image.dat")
        d_ideal = readdlm(ideal_file)

        @test d[  1,  1] ≈ d_ideal[1, 1]
        @test d[ 30, 32] ≈ d_ideal[30, 32]
        @test d[117, 92] ≈ d_ideal[117, 92]

    end

    @testset "Smac utility" begin

        @test_nowarn write_smac1_par("", 0, "", "", "", "",
                                     0, 0, 0, 4, 3,
                                     20.0, 10.0, 1, 1,
                                     24, 1.0, 1.e6, 10,
                                     1, 0.0, 0.0, 0.0)

        # @test_warn "Read error: Incorrect image format!" read_smac1_binary_image(joinpath(dirname(@__FILE__), "image.dat"))

        filename = joinpath(dirname(@__FILE__), "GadJet_test.140.a.z.pix")
        info = read_smac1_binary_info(filename)

        @test info.snap == 140

        image = read_smac1_binary_image(filename)
        @test length(image[:,1]) == 128
        @test image[1,1] ≈ 0.000120693

        @test_nowarn write_smac2_par(1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                                     1.0, 1.0, 1024,
                                     "", "", "")
    end

    @testset "COM" begin

        m = ones(5)

        x = [  3.2  4.1 -2.4;
               2.9 -2.3  1.7;
              -2.3  4.2  3.2;
              -1.1 -2.4 -1.1;
               2.8 -3.2 -3.2 ]

        result = calculate_center_of_mass(x, m)

        @test result[1] ≈  1.1
        @test result[2] ≈  0.07999999999999999
        @test result[3] ≈ -0.36000000000000004

    end

    @testset "Spectral CRs types" begin

        d = CRShockData(1)
        @test d.dt[1] ==  0.0

        d = CRMomentumDistributionConfig()
        @test d.mp == 1.0

        d = CRMomentumDistribution(1)
        @test d.CRp_dis[1] == 0.0

    end

end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
"""
    PolygonMesh(nodes[, edges])

Polygon given by `nodes` and `edges`. Several formats are supported:
 * nodes as nx2 matrix
 * nodes as n-vector of 2-vectors or 2-tuples of x-y-coordinates
 * edges as nx2 matrix of indices
 * edges as n-vector of indices`
 * edges empty

Incomplete edges are filles up by best guess. The columns must be permutations of `1:n`.
"""
struct PolygonMesh{A,U,E<:Union{Nothing,AbstractArray{<:Integer}}}
    nodes::U
    edges::E
    function PolygonMesh(n::U, ed::E) where {U,E}
        A = ed isa AbstractMatrix ? max(size(ed, 2) - 2, 0) : 0
        # the number of stored area indices per edge
        checkarguments(n, ed)
        new{A,U,E}(n, ed)
    end
    PolygonMesh(n::U) where U = PolygonMesh(n, nothing)
end

"""
    PointsInbound(p)

Vector of points. Several formats are supported for points:
 * as 2-vector of real numbers (set of single point)
 * as nx2 matrix of real numbers - each row represents `(x, y)`
 * as n-vector of 2-vectors or 2-tuples of x-y-coordinates
"""
struct PointsInbound{V}
    points::V
    PointsInbound(p::V) where V = new{V}(p)
end

"""
    see `InOnOut` and `InOnBit`
"""
abstract type AbstractOutputFormat end
"""
    InOnOut{IN,ON,OUT}

Expected format of output is one of three constant values.
`InOnOut{1,0,-1}` returns `1`, `0`, `-1` if point is in, on, or out respectively.
"""
struct InOnOut{IN,ON,OUT} <: AbstractOutputFormat end

"""
    InOnBit

Expected output is a pair of booleans, the first one, `inout`, indicating the point
is inside, the second one, `bnds`, indicating the point is close to the border.
Vectors are stored in `BitArray`s of 3 dimensions.

Depending on the used algorithm, the value of `inout` may be invalid, if `bnds` is true.
"""
struct InOnBit <: AbstractOutputFormat end

# interface for above structures standard cases

# access functions for PolygonMesh

function nodecount(poly::PolygonMesh)
    size(poly.nodes, 1)
end
function Base.length(poly::PolygonMesh)
    nodecount(poly)
end
function edgecount(poly::PolygonMesh{<:Any,<:Any,<:Union{Nothing,AbstractArray{<:Any,0}}})
    length(poly)
end
function edgecount(poly::PolygonMesh{<:Any,<:Any,<:AbstractArray})
    size(poly.edges, 1)
end

function edgeindex(poly::PolygonMesh{0,<:Any,<:Union{Nothing,AbstractArray{<:Any,0}}}, i::Integer, n::Integer)
    n == 1 ? i : i % nodecount(poly) + 1
end
function edgeindex(poly::PolygonMesh{<:Any,<:Any,<:AbstractVector}, i::Integer, n::Integer)
    n == 1 ? i : poly.edges[i]
end
function edgeindex(poly::PolygonMesh{<:Any,<:Any,<:AbstractMatrix}, i::Integer, n::Integer)
    poly.edges[i,n]
end

function areaindex(poly::PolygonMesh{A}, i::Integer, a::Integer) where A
    A >= 1 ? poly.edges[i,a+2] : 1
end
function areacount(poly::PolygonMesh{A}) where A
    if A >= 1
        maximum(view(poly.edges, :, 3:A+2))
    else
        1
    end
end

function all_areas(poly::PolygonMesh{A}) where A
    if A >= 1
        edges = poly.edges
        m = size(edges, 2)
        m <= 2 ? [1] : [a for a in unique(vec(view(edges,:, 3:m))) if a != 0]
    else
        [1]
    end
end

function edges_dict_for_area(poly::PolygonMesh{A}, a::Integer) where A
    edges = poly.edges
    n, m = size(edges)
    if A >= 1
        has(edges, i, a) = any(view(edges, i, 3:m) .== a)
        Dict((edges[k,1], edges[k,2]) for k in 1:n if has(edges, k, a))
    else
        Dict((edges[k,1], edges[k,2]) for k in 1:n)
    end
end

function all_edges(poly::PolygonMesh)
    ((edgeindex(poly, i, 1), edgeindex(poly, i, 2)) for i in 1:edgecount(poly))
end

function all_boxlengths(poly::PolygonMesh, xy::Integer)
    (abs(vertex(poly,v,xy) - vertex(poly,w,xy)) for (v, w) in all_edges(poly))
end

function vertex(poly::PolygonMesh{<:Any,<:AbstractMatrix{<:Real}}, v::Integer, xy::Integer)
    poly.nodes[v,xy]
end
function vertex(poly::PolygonMesh{<:Any,<:AbstractVector}, v::Integer, xy::Integer)
    poly.nodes[v][xy]
end

# access functions for PointsInbounds

function Base.length(p::PointsInbound{<:AbstractVector{<:Real}})
    1
end
function Base.length(p::PointsInbound{<:AbstractMatrix{<:Real}})
    size(p.points, 1)
end
function Base.length(p::PointsInbound{<:AbstractVector})
    length(p.points)
end

function vertex(p::PointsInbound{<:AbstractVector{<:Real}}, v::Integer, xy::Integer)
    p.points[xy]
end
function vertex(p::PointsInbound{<:AbstractMatrix{<:Real}}, v::Integer, xy::Integer)
    p.points[v,xy]
end
function vertex(p::PointsInbound{<:AbstractVector}, v::Integer, xy::Integer)
    p.points[v][xy]
end

# standard functions for PointsInbound

function Base.minimum(p::Union{PointsInbound,PolygonMesh})
    n = length(p)
    minimum([[vertex(p, k, 1) for k = 1:n] [vertex(p, k, 2) for k = 1:n]], dims = 1)
end
function Base.maximum(p::Union{PointsInbound,PolygonMesh})
    n = length(p)
    maximum([[vertex(p, k, 1) for k = 1:n] [vertex(p, k, 2) for k = 1:n]], dims = 1)
end
function Base.sortperm(p::PointsInbound, ix::Integer)
    n = length(p)
    sortperm([vertex(p, k, ix) for k in 1:n])
end
function Base.sortperm(p::PointsInbound{<:AbstractMatrix}, ix::Integer)
    sortperm(view(p.points, :, ix))
end

# conversion between different output formats

function convertout(to::Type{T}, from::Type{S}, args...) where {S<:AbstractOutputFormat,T<:AbstractOutputFormat}
    T == S && return length(args) == 1 ? args[1] : args
    convertto(to, convertfrom(from, args...)...)
end

function convertfrom(::Type{<:InOnBit}, stat::T, bnds::T) where T<:Bool
    stat, bnds
end
function convertfrom(::Type{<:InOnBit}, stat::T) where T<:AbstractArray{Bool}
    stat
end
function convertfrom(::Type{InOnOut{IN,ON,OUT}}, stat::Integer) where {IN,ON,OUT}
    if stat == IN
        true, false
    elseif stat == OUT
        false, false
    elseif stat == ON
        true, true
    else
        false, true
    end
end
function convertfrom(::Type{OT}, stat::Union{AbstractVector{T},AbstractMatrix{T}}) where {T<:Integer,OT<:InOnOut}
    n, m = size(stat)
    inout = BitArray(undef, n, 2, m)
    for j = 1:m
        for i in 1:n
            a, b = convertfrom(OT, stat[i,j])
            inout[i,1,j] = a
            inout[i,2,j] = a
        end
    end
    inout, bnds
end

function convertto(::Type{<:InOnBit}, stat::T, bnds::T) where T<:Bool
    stat, bnds
end
function convertto(::Type{<:InOnBit}, stat::T, bnds::T) where T<:BitArray
    stat
end
function convertto(::Type{InOnOut{IN,ON,OUT}}, stat::Bool, bnds::Bool) where {IN,ON,OUT}
    bnds ? ON : stat ? IN : OUT
end
function convertto(::Type{OT}, stat::T) where {OT<:InOnOut,T<:AbstractArray{Bool}}
    convertto.(OT, view(stat, :,1,:), view(stat, :, 2, :))
end

# check functions for PolygonMesh constructor

function checknodes(nodes::AbstractMatrix)
    size(nodes, 2) >= 2 || throw(ArgumentError("nodes require at least 2 coordinates"))
    nothing
end
function checknodes(nodes::AbstractVector)
    nothing
end
function checknodes(nodes)
    nothing
end

function checkarguments(nodes, edges::Union{Nothing,AbstractArray{<:Integer,0}})
    checknodes(nodes)
    nothing
end
function checkarguments(nodes, edges::AbstractVector{<:Integer})
    checknodes(nodes)
    n = size(nodes, 1)
    mi, ma = extrema(edges)
    0 < mi <= ma <= n || throw(ArgumentError("invalid edge indices"))
    nothing
end
function checkarguments(nodes, edges::AbstractMatrix{<:Integer})
    checknodes(nodes)
    n, m = size(edges)
    m >= 2 || throw(ArgumentError("edges matrix requires at least 2 columns"))
    m <= 5 || throw(ArgumentError("edges matrix accepts at most 3 area columns"))
    mi, ma = extrema(view(edges, :, 1:2))
    0 < mi <= ma <= n || throw(ArgumentError("invalid edge indices"))
    if m > 2
        mi, ma = extrema(view(edges, :, 3:m))
        0 <= mi <= 4096 || throw(ArgumentError("invalid area indices ($mi, $ma)"))
    end
    if m <= 2 || m == 3 && mi == ma
        check(i) = isperm(view(edges, :, i)) ||
        throw(ArgumentError("edges column $i is no permutation"))
        check(1)
        check(2)
    end
    nothing
end


text:
module LineNumbers

export SourceFile, compute_line, LineBreaking

# Offsets are 0 based
struct SourceFile
    data::Vector{UInt8}
    offsets::Vector{UInt64}
end

function SourceFile(data::AbstractString)
    offsets = UInt64[0]
    buf = IOBuffer(data)
    local line
    while !eof(buf)
        line = readuntil(buf,'\n')
        !eof(buf) && push!(offsets, position(buf))
    end
    if !isempty(offsets) && !isempty(line) && line[end] == '\n'
        push!(offsets, position(buf))
    end
    SourceFile(copy(codeunits(data)), offsets)
end

function compute_line(file::SourceFile, offset::Integer)
    ind = searchsortedfirst(file.offsets, offset)
    ind <= length(file.offsets) && file.offsets[ind] == offset ? ind : ind - 1
end

end # module


text:
using OpenSCADPolyhedron
using OpenSCADCube

function generateBaseFace (d::Float64)
  return [
    # four base points
    Point3{Float64}([ 1.0,  1.0, -1.0]),
    Point3{Float64}([-1.0,  1.0, -1.0]),
    Point3{Float64}([-1.0, -1.0, -1.0]),
    Point3{Float64}([ 1.0, -1.0, -1.0]),
    # base points taken in by d
    Point3{Float64}([ 1.0 - d,  1.0 - d, -1.0]),
    Point3{Float64}([-1.0 + d,  1.0 - d, -1.0]),
    Point3{Float64}([-1.0 + d, -1.0 + d, -1.0]),
    Point3{Float64}([ 1.0 - d, -1.0 + d, -1.0]),
    # now taken up by d as well
    Point3{Float64}([ 1.0 - d,  1.0 - d, -1.0 + d]),
    Point3{Float64}([-1.0 + d,  1.0 - d, -1.0 + d]),
    Point3{Float64}([-1.0 + d, -1.0 + d, -1.0 + d]),
    Point3{Float64}([ 1.0 - d, -1.0 + d, -1.0 + d]),
  ]
end

cubeFrameBase = Polyhedron{Float64}(generateBaseFace(0.2))
# bottom faces
append!(cubeFrameBase, Int[1, 2, 6, 5])
append!(cubeFrameBase, Int[2, 3, 7, 6])
append!(cubeFrameBase, Int[3, 4, 8, 7])
append!(cubeFrameBase, Int[4, 1, 5, 8])
# inside faces
append!(cubeFrameBase, Int[10, 9, 5, 6])
append!(cubeFrameBase, Int[11, 10, 6, 7])
append!(cubeFrameBase, Int[12, 11, 7, 8])
append!(cubeFrameBase, Int[9, 12, 8, 5])

# These are the points that need symmetry relations.
# Note that they must be specified in the order that lines them up with cubeFace
cubeFrameLayers = Set{Array{Int,1}}()
push!(cubeFrameLayers, Int[1, 2, 3, 4])
push!(cubeFrameLayers, Int[9, 10, 11, 12])

# Here are our symmetries, generated from the basic cube symmetries
cubeFrameSymmetries = expandSymmetries(cubeSymmetries, cubeFace, cubeFrameLayers)

cube = generatePolyhedron(cubeFrameBase, cubeTransforms, cubeFrameSymmetries)
println(stringify(cube))

max_stars_repo_name: erikstrand/OpenSCADPolyhedron

max_issues_repo_name: erikstrand/OpenSCADPolyhedron

max_forks_repo_name: erikstrand/OpenSCADPolyhedron


text:
"""
    getNominalReturn
calculates the return of a "profit and loss"-series as absolute value
per definition: the net profit or loss of an investment defined by the change
in value over a stated period including all costs and incomes
Formula is:     return = (endValue - initialValue)/initialValue
"""
function getNominalReturn(arr::Array{Float64,1})::Float64
    return Base.sum(arr)
end

"""
    getReturn_rate
calculate the return of a tradeSeries as a rate
formula: returnRate = p(t)/p(0)-1 = (p(t) - p(0))/p(0)
"""
function getReturn_rate(profitAndLossArray::Array{Float64,1}, initialValue::Float64)::Float64
    # if the initial value is 0 or less abort the programm and inform user
    if initialValue <= 0
        Error("Initial Value is 0 or less: See function getReturnRate");
        exit(0);
    end
    endValue = getNominalReturn(profitAndLossArray)+initialValue;
    return ((endValue - initialValue) / initialValue)
end

"""
    getReturn_percentage
calculate the return of a tradeSeries as percentage
"""
function getReturn_percentage(profitAndLossArray::Array{Float64,1}, initialValue::Float64)::Float64
    return getReturn_rate(profitAndLossArray,initialValue)*100
end

"""
    convertReturnToAnnulizedReturn
convert the return Rate to annulized Return with periods
"""
function convertReturnToAnnulizedReturn(returnRate::Float64, amountOfPeriods::Float64)::Float64
    return round(((1+returnRate)^(1/amountOfPeriods))-1,digits=3)
end

"""
    getAnnulizedReturnRateBetweenDates
calculates the annulized Return between twoDates for yearly scale
"""
function getAnnulizedReturnRateBetweenDatesYearly(profitAndLossArray::Array{Float64,1}, capitalAtBeginn::Float64, beginDate::Date, endDate::Date)::Float64
    averageAmountOfDaysPerYear::Float64 = 365.25                                #Average amount of days per year taking in account leap years
    returnRate::Float64 = getReturnRate(profitAndLossArray,capitalAtBeginn)
    amountOfDays = Dates.days((endDate + Dates.Day(1)) - beginDate)             # +1 to set the end-day to a closed-interval: 01.01.2018 to 31.12.2018 includes the 31.12.2018
    return ((1+returnRate)^(365.25/amountOfDays))-1                             # 365.25 amount of days per year taking in account leap years
end

"""
    getYearlyAverageProfit_betweenDates
calculates the average yearly return for a given array of profits and losses
between two given Dates
"""

function getYearlyAverageProfit_betweenDates(profitAndLossArray::Array{Float64,1},beginDate::Date, endDate::Date)::Float64
    nominalReturn = getNominalReturn(profitAndLossArray);
    amountOfDaysBetween = getDaysBetweenDates_includeLast(beginDate, endDate);
    amountOfYearsBetween = getDaysBetweenDates_includeLast(beginDate,endDate)/365.25
    return nominalReturn / amountOfYearsBetween
end

"""
    getYearlyAverageProfit_percentage
calculates the average yearly return for a given array of profits and losses
and returns the percentage values
"""

function getYearlyAverageProfit_percentage(profitAndLossArray::Array{Float64,1},beginDate::Date, endDate::Date, initialCapital::Float64)::Float64
    yearlyAverageProfit::Float64 = getYearlyAverageProfit(profitAndLossArray, beginDate,endDate)
    return (yearlyAverageProfit/initialCapital)*100
end

max_stars_repo_name: RainerNeu/tradeStatistics

max_issues_repo_name: RainerNeu/tradeStatistics

max_forks_repo_name: RainerNeu/tradeStatistics


text:
module ReverseDiffFlux

using Flux, ReverseDiff
import Flux.NNlib: softmax, ∇softmax, logsoftmax, ∇logsoftmax, conv, ∇conv_data, depthwiseconv, maxpool, meanpool
import Flux.NNlib: DenseConvDims, DepthwiseConvDims, PoolDims
using MacroTools: @forward

include("nnlib.jl")
include("flux.jl")

function ReverseDiff.gradient(f, ps::Flux.Params)
    input = tuple(ps...)
    @assert !isempty(input)
    tape = input[1].tape
    @assert all(x -> x.tape === tape, input)
    cfg = ReverseDiff._GradientConfig(input, tape)
    output = ReverseDiff.gradient((a...) -> f(), input, cfg)
    return Flux.Zygote.Grads(Base.IdDict(zip(input, output)), ps)
end

function overload_gradient()
    @eval Flux.gradient(f, args...) = ReverseDiff.gradient(f, args...)
end

end


text:
@doc raw"""
    GeneralizedStiefel{n,k,𝔽,B} <: AbstractDecoratorManifold{𝔽}

The Generalized Stiefel manifold consists of all $n\times k$, $n\geq k$ orthonormal
matrices w.r.t. an arbitrary scalar product with symmetric positive definite matrix
$B\in R^{n × n}$, i.e.

````math
\operatorname{St}(n,k,B) = \bigl\{ p \in \mathbb F^{n × k}\ \big|\ p^{\mathrm{H}} B p = I_k \bigr\},
````

where $𝔽 ∈ \{ℝ, ℂ\}$,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and
$I_k \in \mathbb R^{k × k}$ denotes the $k × k$ identity matrix.

In the case $B=I_k$ one gets the usual [`Stiefel`](@ref) manifold.

The tangent space at a point $p\in\mathcal M=\operatorname{St}(n,k,B)$ is given by

````math
T_p\mathcal M = \{ X \in 𝔽^{n × k} : p^{\mathrm{H}}BX + X^{\mathrm{H}}Bp=0_n\},
````
where $0_k$ is the $k × k$ zero matrix.

This manifold is modeled as an embedded manifold to the [`Euclidean`](@ref), i.e.
several functions like the [`zero_vector`](@ref) are inherited from the embedding.

The manifold is named after
[Eduard L. Stiefel](https://en.wikipedia.org/wiki/Eduard_Stiefel) (1909–1978).

# Constructor
    GeneralizedStiefel(n, k, B=I_n, F=ℝ)

Generate the (real-valued) Generalized Stiefel manifold of $n\times k$ dimensional
orthonormal matrices with scalar product `B`.
"""
struct GeneralizedStiefel{n,k,𝔽,TB<:AbstractMatrix} <: AbstractDecoratorManifold{𝔽}
    B::TB
end

function GeneralizedStiefel(
    n::Int,
    k::Int,
    B::AbstractMatrix=Matrix{Float64}(I, n, n),
    𝔽::AbstractNumbers=ℝ,
)
    return GeneralizedStiefel{n,k,𝔽,typeof(B)}(B)
end

active_traits(f, ::GeneralizedStiefel, args...) = merge_traits(IsEmbeddedManifold())

@doc raw"""
    check_point(M::GeneralizedStiefel, p; kwargs...)

Check whether `p` is a valid point on the [`GeneralizedStiefel`](@ref) `M`=$\operatorname{St}(n,k,B)$,
i.e. that it has the right [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system) type and $x^{\mathrm{H}}Bx$
is (approximately) the identity, where $\cdot^{\mathrm{H}}$ is the complex conjugate
transpose. The settings for approximately can be set with `kwargs...`.
"""
function check_point(M::GeneralizedStiefel{n,k,𝔽}, p; kwargs...) where {n,k,𝔽}
    c = p' * M.B * p
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(p) does not lie on $(M), because x'Bx is not the unit matrix.",
        )
    end
    return nothing
end

# overwrite passing to embedding
function check_size(M::GeneralizedStiefel{n,k,𝔽}, p) where {n,k,𝔽}
    return check_size(get_embedding(M), p) #avoid embed, since it uses copyto!
end
function check_size(M::GeneralizedStiefel{n,k,𝔽}, p, X) where {n,k,𝔽}
    return check_size(get_embedding(M), p, X) #avoid embed, since it uses copyto!
end

@doc raw"""
    check_vector(M::GeneralizedStiefel, p, X; kwargs...)

Check whether `X` is a valid tangent vector at `p` on the [`GeneralizedStiefel`](@ref)
`M`=$\operatorname{St}(n,k,B)$, i.e. the [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system) fits,
`p` is a valid point on `M` and
it (approximately) holds that $p^{\mathrm{H}}BX + \overline{X^{\mathrm{H}}Bp} = 0$, where
`kwargs...` is passed to the `isapprox`.
"""
function check_vector(M::GeneralizedStiefel{n,k,𝔽}, p, X; kwargs...) where {n,k,B,𝔽}
    if !isapprox(p' * M.B * X, -conj(X' * M.B * p); kwargs...)
        return DomainError(
            norm(p' * M.B * X + conj(X' * M.B * p)),
            "The matrix $(X) does not lie in the tangent space of $(p) on $(M), since x'Bv + v'Bx is not the zero matrix.",
        )
    end
    return nothing
end

get_embedding(::GeneralizedStiefel{N,K,𝔽}) where {N,K,𝔽} = Euclidean(N, K; field=𝔽)

@doc raw"""
    inner(M::GeneralizedStiefel, p, X, Y)

Compute the inner product for two tangent vectors `X`, `Y` from the tangent space of `p`
on the [`GeneralizedStiefel`](@ref) manifold `M`. The formula reads

````math
(X, Y)_p = \operatorname{trace}(v^{\mathrm{H}}Bw),
````
i.e. the metric induced by the scalar product `B` from the embedding, restricted to the
tangent space.
"""
inner(M::GeneralizedStiefel, p, X, Y) = dot(X, M.B * Y)

@doc raw"""
    manifold_dimension(M::GeneralizedStiefel)

Return the dimension of the [`GeneralizedStiefel`](@ref) manifold `M`=$\operatorname{St}(n,k,B,𝔽)$.
The dimension is given by

````math
\begin{aligned}
\dim \mathrm{St}(n, k, B, ℝ) &= nk - \frac{1}{2}k(k+1) \\
\dim \mathrm{St}(n, k, B, ℂ) &= 2nk - k^2\\
\dim \mathrm{St}(n, k, B, ℍ) &= 4nk - k(2k-1)
\end{aligned}
````
"""
function manifold_dimension(::GeneralizedStiefel{n,k,ℝ}) where {n,k}
    return n * k - div(k * (k + 1), 2)
end
manifold_dimension(::GeneralizedStiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k
manifold_dimension(::GeneralizedStiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1)

@doc raw"""
    project(M::GeneralizedStiefel,p)

Project `p` from the embedding onto the [`GeneralizedStiefel`](@ref) `M`, i.e. compute `q`
as the polar decomposition of $p$ such that $q^{\mathrm{H}}Bq$ is the identity,
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
project(::GeneralizedStiefel, ::Any)

function project!(M::GeneralizedStiefel, q, p)
    s = svd(p)
    e = eigen(s.U' * M.B * s.U)
    qsinv = e.vectors * Diagonal(1 ./ sqrt.(e.values))
    q .= s.U * qsinv * e.vectors' * s.V'
    return q
end

@doc raw"""
    project(M:GeneralizedStiefel, p, X)

Project `X` onto the tangent space of `p` to the [`GeneralizedStiefel`](@ref) manifold `M`.
The formula reads

````math
\operatorname{proj}_{\operatorname{St}(n,k)}(p,X) = X - p\operatorname{Sym}(p^{\mathrm{H}}BX),
````

where $\operatorname{Sym}(y)$ is the symmetrization of $y$, e.g. by
$\operatorname{Sym}(y) = \frac{y^{\mathrm{H}}+y}{2}$.
"""
project(::GeneralizedStiefel, ::Any, ::Any)

function project!(M::GeneralizedStiefel, Y, p, X)
    A = p' * M.B' * X
    copyto!(Y, X - p * Hermitian((A + A') / 2))
    return Y
end
@doc raw"""
    retract(M::GeneralizedStiefel, p, X)
    retract(M::GeneralizedStiefel, p, X, ::PolarRetraction)
    retract(M::GeneralizedStiefel, p, X, ::ProjectionRetraction)

Compute the SVD-based retraction [`PolarRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.PolarRetraction) on the
[`GeneralizedStiefel`](@ref) manifold `M`, which in this case is the same as
the projection based retraction employing the exponential map in the embedding
and projecting the result back to the manifold.

The default retraction for this manifold is the [`ProjectionRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.ProjectionRetraction).
"""
retract(::GeneralizedStiefel, ::Any...)

default_retraction_method(::GeneralizedStiefel) = ProjectionRetraction()

function retract_polar!(M::GeneralizedStiefel, q, p, X)
    project!(M, q, p + X)
    return q
end
function retract_project!(M::GeneralizedStiefel, q, p, X)
    project!(M, q, p + X)
    return q
end

function Base.show(io::IO, M::GeneralizedStiefel{n,k,𝔽}) where {n,k,𝔽}
    return print(io, "GeneralizedStiefel($(n), $(k), $(M.B), $(𝔽))")
end


text:
## Script to determine values for filtering the datasets

using Revise
using SwitchingAnalysis
using HypothesisTests
using BrowseTables
using Plots.PlotMeasures

## load
if ispath(linux_gdrive)
    ongoing_dir = linux_gdrive
else
    ongoing_dir = mac_gdrive
end

files_loc = joinpath(ongoing_dir,files_dir)
figs_loc = joinpath(ongoing_dir,figs_dir)
fullS =  CSV.read(joinpath(files_loc,"streaks.csv"), DataFrame; types = columns_types)
fullP =  CSV.read(joinpath(files_loc,"pokes.csv"), DataFrame)
gr(size=(600,600), tick_orientation = :out, grid = false, linecolor = :black,
markerstrokecolor = :black)
## removing irrelevant events and computes relevant variables for pokes dataframe

# turning Protocol from Float to string to allow for categorical analysis
fullP[!,:Protocol] = [ismissing(x) ? missing : string(x) for x in fullP[:,:Protocol]]
pokes = filter(r-> !ismissing(r.Protocol) &&
    r.ExpDay > 5 && #remove early training days
    r.Protocol in ["0.5","0.75","1.0"], #remove protocol "0.0": an exeption for wrong trials
    fullP)
#removing events when the animal travel to a patch but doesn't poke
pokes =  dropmissing(pokes, :Poke)
disallowmissing!(pokes, :Reward)
# computing the probability of next poke to be a reward and instanteneous reward rate
pokes[!,:Next_Prew] = [Prew(parse(Float64,prot) + 1,Int64(pok)) for (prot,pok) in zip(pokes.Protocol,pokes.Poke_within_Trial)]
gd = groupby(pokes, [:Day,:MouseID,:Trial])
transform!(gd,:Reward => Pnext => :Pnextrew)
# calculate time interval from previous poke out to current poke out
poke_time = [ismissing(i) ? d : i+d for (i,d) in zip( pokes.Pre_Interpoke, pokes.PokeDuration)]
pokes[!,:InstRewRate] = pokes.Pnextrew ./ poke_time
#calculate elapsed time and cumulative reward from the first poke in
gd = groupby(pokes, [:Day,:MouseID])
transform!(gd,[:PokeIn,:PokeOut,] => (i,o) -> o .- i[1])
rename!(pokes, :PokeIn_PokeOut_function => :ElapsedTime)
transform!(gd,:Reward => cumsum => :Cumulative_Reward)
# calculate average reward rate for each poke cumulative reward / elapsed time
pokes[!,:AverageRewRate] = 1 ./ (pokes.ElapsedTime ./ pokes.Cumulative_Reward)
# calculate time elapsed from the last pokeout of the previous trial
pokes[!,:TimeFromLeaving] = zeros(nrow(pokes))
gd1 = groupby(pokes,[:MouseID,:Day])
combine(gd1) do dd
    dd[:,:TimeFromLeaving] = TimeFromLeaving(dd)
end
nbins = 200
lr = LinRange(0,60,nbins)
pokes[!,:TimeFromLeaving] = [isnothing(findfirst( x .< lr)) ? round(nbins * step(lr), digits = 1) : round(findfirst( x .< lr) * step(lr), digits = 1) for x in pokes.TimeFromLeaving]
pokes[!,:CumRewTrial] = Vector{Float64}(undef,nrow(pokes))
gd2 = groupby(pokes,[:Trial,:MouseID,:Day])
combine(gd2) do dd
    dd[:,:CumRewTrial] = cumsum(dd.Reward)
end
gd3 = groupby(pokes,[:Protocol])
combine(gd3) do dd
    dd[:,:CumRewTrial] = dd[:,:CumRewTrial]./maximum(dd.CumRewTrial)
end
checkstim = combine(groupby(pokes,[:MouseID,:ExpDay]), :Stim => (s -> length(union(s)) > 1) => :StimDay)
checkstim[!,:StimDay] = [m in ["pc10", "pc1", "pc2", "pc3", "pc4", "pc5", "pc6", "pc7", "pc8", "pc9"] ? s : false for (m,s) in zip(checkstim.MouseID, checkstim.StimDay)]
expanded = leftjoin(pokes, checkstim, on = [:MouseID,:ExpDay]).StimDay
pokes.Stim_Day = expanded

# open_html_table(pokes[1:100,[:Side,:Poke,:Trial,:PokeIn,:PokeOut,:TimeFromLeaving]])
## Process trials dataframe
streaks = combine(groupby(pokes,[:MouseID,:Day,:Phase,:Group,:Treatment, :Injection])) do dd
        process_streaks(dd)
    end
##Process trials dataframe
count_bouts!(pokes)
bouts = combine(groupby(pokes,[:MouseID,:Day,:Phase,:Group,:Treatment, :Injection])) do dd
    process_bouts(dd)
end
##
#= Filter dataset cutting the 5th percentile tail (either bilaterally or to the right extreme)
if filtering results skewed a filter for values with a probability lower than 0.99 is applied =#

# trim_conf_ints!(pokes,:Pre_Interpoke)
# trim_conf_ints!(pokes,:PokeDuration)
#
# trim_conf_ints!(streaks,:Trial)
#= the distribution of number of pokes before leaving results bimodal;
    with early leaving seaprated from the rest. Analysis on trials
    are therefore performed with at least 2 or more pokes before leaving=#
# filter!(r-> 1< r.Num_pokes < 31, streaks)
# trim_conf_ints!(streaks, :Num_pokes; percent = 99)
# filter!(r-> r.Trial <61,streaks)
##
# @df streaks density(:Num_pokes)
# @df trim_conf_ints(streaks, :Num_pokes) density(:Num_pokes)


text:
module ADBCVUMPS

using Zygote
using OMEinsum
using BCVUMPS
using BCVUMPS: HamiltonianModel

export num_grad
export TFIsing, Heisenberg
export hamiltonian, HamiltonianModel
export init_ipeps, energy, optimiseipeps

include("hamiltonianmodels.jl")
include("autodiff.jl")
include("bcipeps.jl")
include("variationalipeps.jl")

end


text:
# each of these includes the bits it needs so don't do all includes here
include("harmonic_test.jl")
include("bernoulli_test.jl")
include("stieltjes_test.jl")
include("gamma_test.jl")
include("polylog_test1.jl")
include("polylog_test2.jl")


text:
module ConicProgram

using LinearAlgebra, SparseArrays

using MathOptInterface
const MOI = MathOptInterface

import BlockDiagonals
import IterativeSolvers

import DiffOpt

Base.@kwdef struct Cache
    M::SparseMatrixCSC{Float64, Int}
    vp::Vector{Float64}
    Dπv::BlockDiagonals.BlockDiagonal{Float64, Matrix{Float64}}
    A::SparseMatrixCSC{Float64, Int}
    b::Vector{Float64}
    c::Vector{Float64}
end

Base.@kwdef struct ForwCache
    du::Vector{Float64}
    dv::Vector{Float64}
    dw::Vector{Float64}
end
Base.@kwdef struct ReverseCache
    g::Vector{Float64}
    πz::Vector{Float64}
end

# Geometric conic standard form
const Form{T} = MOI.Utilities.GenericModel{
    T,
    MOI.Utilities.ObjectiveContainer{T},
    MOI.Utilities.FreeVariables,
    MOI.Utilities.MatrixOfConstraints{
        T,
        MOI.Utilities.MutableSparseMatrixCSC{
            T,
            Int,
            # We use `OneBasedIndexing` as it is the same indexing as used
            # by `SparseMatrixCSC` so we can do an allocation-free conversion to
            # `SparseMatrixCSC`.
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{T},
        DiffOpt.ProductOfSets{T},
    },
}

"""
    Diffopt.ConicProgram.Model <: DiffOpt.AbstractModel

Model to differentiate conic programs.

The forward differentiation computes the product of the derivative (Jacobian) at the
conic program parameters `A`, `b`, `c` to the perturbations `dA`, `db`, `dc`.

The reverse differentiation computes the product of the transpose of the derivative (Jacobian) at the
conic program parameters `A`, `b`, `c` to the perturbations `dx`, `dy`, `ds`.

For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043
"""
mutable struct Model <: DiffOpt.AbstractModel
    # storage for problem data in matrix form
    model::Form{Float64}
    # includes maps from matrix indices to problem data held in `optimizer`
    # also includes KKT matrices
    # also includes the solution
    gradient_cache::Union{Nothing,Cache}

    # caches for sensitivity output
    # result from solving KKT/residualmap linear systems
    # this allows keeping the same `gradient_cache`
    # if only sensitivy input changes
    forw_grad_cache::Union{Nothing,ForwCache}
    back_grad_cache::Union{Nothing,ReverseCache}

    # sensitivity input cache using MOI like sparse format
    input_cache::DiffOpt.InputCache

    x::Vector{Float64} # Primal
    s::Vector{Float64} # Slack
    y::Vector{Float64} # Dual
end
function Model()
    return Model(Form{Float64}(), nothing, nothing, nothing, DiffOpt.InputCache(), Float64[], Float64[], Float64[])
end

function MOI.is_empty(model::Model)
    return MOI.is_empty(model.model)
end

function MOI.empty!(model::Model)
    MOI.empty!(model.model)
    model.gradient_cache = nothing
    model.forw_grad_cache = nothing
    model.back_grad_cache = nothing
    empty!(model.input_cache)
    empty!(model.x)
    empty!(model.s)
    empty!(model.y)
    return
end

function MOI.supports_constraint(model::Model, F::Type{MOI.VectorAffineFunction{Float64}}, ::Type{S}) where {S<:MOI.AbstractVectorSet}
    if DiffOpt.add_set_types(model.model.constraints.sets, S)
        push!(model.model.constraints.caches, Tuple{F,S}[])
        push!(model.model.constraints.are_indices_mapped, BitSet())
    end
    return MOI.supports_constraint(model.model, F, S)
end

function MOI.set(model::Model, ::MOI.ConstraintPrimalStart, ci::MOI.ConstraintIndex, value)
    MOI.throw_if_not_valid(model, ci)
    DiffOpt._enlarge_set(model.s, MOI.Utilities.rows(model.model.constraints, ci), value)
end

function MOI.set(model::Model, ::MOI.ConstraintDualStart, ci::MOI.ConstraintIndex, value)
    MOI.throw_if_not_valid(model, ci)
    DiffOpt._enlarge_set(model.y, MOI.Utilities.rows(model.model.constraints, ci), value)
end

function _gradient_cache(model::Model)
    if model.gradient_cache !== nothing
        return model.gradient_cache
    end

    # For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043

    A = -convert(SparseMatrixCSC{Float64, Int}, model.model.constraints.coefficients)
    b = model.model.constraints.constants

    if MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
        c = spzeros(size(A, 2))
    else
        obj = MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        c = DiffOpt.sparse_array_representation(obj, size(A, 2)).terms
        if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
            c = -c
        end
    end

    # programs in tests were cross-checked against `diffcp`, which follows SCS format
    # hence, some arrays saved during `MOI.optimize!` are not same across all optimizers
    # specifically there's an extra preprocessing step for `PositiveSemidefiniteConeTriangle` constraint for SCS/Mosek

    # pre-compute quantities for the derivative
    m = A.m
    n = A.n
    N = m + n + 1
    # NOTE: w = 1.0 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (model.x, model.y - model.s, 1.0)

    # find gradient of projections on dual of the cones
    Dπv = DiffOpt.Dπ(v, model.model, model.model.constraints.sets)

    # Q = [
    #      0   A'   c;
    #     -A   0    b;
    #     -c' -b'   0;
    # ]
    # M = (Q- I) * B + I
    # with B =
    # [
    #  I    .   .    # Πx = x because x is a solution and hence satistfies the constraints
    #  .  Dπv   .
    #  .    .   1    # w >= 0, but in the solution x = 1
    # ]
    # see: https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
    # for the definition of Π and why we get I and 1 for x and w respectectively
    # K is defined in (5), Π in sect 2, and projection sin sect 3

    M = [
        spzeros(n,n)     (A' * Dπv)    c
        -A               -Dπv + I      b
        -c'              -b' * Dπv     0.0
    ]
    # find projections on dual of the cones
    vp = DiffOpt.π(v, model.model, model.model.constraints.sets)

    model.gradient_cache = Cache(
        M = M,
        vp = vp,
        Dπv = Dπv,
        A = A,
        b = b,
        c = c,
    )

    return model.gradient_cache
end

function DiffOpt.forward_differentiate!(model::Model)
    gradient_cache = _gradient_cache(model)
    M = gradient_cache.M
    vp = gradient_cache.vp
    Dπv = gradient_cache.Dπv
    x = model.x
    y = model.y
    s = model.s
    A = gradient_cache.A
    b = gradient_cache.b
    c = gradient_cache.c

    objective_function = DiffOpt._convert(MOI.ScalarAffineFunction{Float64}, model.input_cache.objective)
    sparse_array_obj = DiffOpt.sparse_array_representation(objective_function, length(c))
    dc = sparse_array_obj.terms

    db = zeros(length(b))
    DiffOpt._fill(S -> false, gradient_cache, model.input_cache, model.model.constraints.sets, db)
    (lines, cols) = size(A)
    nz = nnz(A)
    dAi = zeros(Int, 0)
    dAj = zeros(Int, 0)
    dAv = zeros(Float64, 0)
    sizehint!(dAi, nz)
    sizehint!(dAj, nz)
    sizehint!(dAv, nz)
    DiffOpt._fill(S -> false, gradient_cache, model.input_cache, model.model.constraints.sets, dAi, dAj, dAv)
    dA = sparse(dAi, dAj, dAv, lines, cols)

    m = size(A, 1)
    n = size(A, 2)
    N = m + n + 1
    # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # g = dQ * Π(z/|w|) = dQ * [u, vp, 1.0]
    RHS = [dA' * vp + dc; -dA * u + db; -dc ⋅ u - db ⋅ vp]

    dz = if norm(RHS) <= 1e-400 # TODO: parametrize or remove
        RHS .= 0 # because M is square
    else
        IterativeSolvers.lsqr(M, RHS)
    end

    du, dv, dw = dz[1:n], dz[n+1:n+m], dz[n+m+1]
    model.forw_grad_cache = ForwCache(du, dv, [dw])
    return nothing
    # dx = du - x * dw
    # dy = Dπv * dv - y * dw
    # ds = Dπv * dv - dv - s * dw
    # return -dx, -dy, -ds
end

function DiffOpt.reverse_differentiate!(model::Model)
    gradient_cache = _gradient_cache(model)
    M = gradient_cache.M
    vp = gradient_cache.vp
    Dπv = gradient_cache.Dπv
    x = model.x
    y = model.y
    s = model.s
    A = gradient_cache.A
    b = gradient_cache.b
    c = gradient_cache.c

    dx = zeros(length(c))
    for (vi, value) in model.input_cache.dx
        dx[vi.value] = value
    end
    dy = zeros(length(b))
    ds = zeros(length(b))

    m = size(A, 1)
    n = size(A, 2)
    N = m + n + 1
    # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # dz = D \phi (z)^T (dx,dy,dz)
    dz = [
        dx
        Dπv' * (dy + ds) - ds
        - x' * dx - y' * dy - s' * ds
    ]

    g = if norm(dz) <= 1e-4 # TODO: parametrize or remove
        dz .= 0 # because M is square
    else
        IterativeSolvers.lsqr(M, dz)
    end

    πz = [
        u
        vp
        1.0
    ]

    # TODO: very important
    # contrast with:
    # http://reports-archive.adm.cs.cmu.edu/anon/2019/CMU-CS-19-109.pdf
    # pg 97, cap 7.4.2

    model.back_grad_cache = ReverseCache(g, πz)
    return nothing
    # dQ = - g * πz'
    # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
    # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
    # dc = - dQ[1:n, end] + dQ[end, 1:n]'
    # return dA, db, dc
end

function MOI.get(model::Model, ::DiffOpt.ReverseObjectiveFunction)
    g = model.back_grad_cache.g
    πz = model.back_grad_cache.πz
    dc = DiffOpt.lazy_combination(-, πz, g, length(g))
    return DiffOpt.VectorScalarAffineFunction(dc, 0.0)
end

function MOI.get(model::Model, ::DiffOpt.ForwardVariablePrimal, vi::MOI.VariableIndex)
    i = vi.value
    du = model.forw_grad_cache.du
    dw = model.forw_grad_cache.dw
    return -(du[i] - model.x[i] * dw[])
end
function DiffOpt._get_db(model::Model, ci::MOI.ConstraintIndex{F,S}
) where {F<:MOI.AbstractVectorFunction,S}
    i = MOI.Utilities.rows(model.model.constraints, ci) # vector
    # i = ci.value
    n = length(model.x) # columns in A
    # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
    g = model.back_grad_cache.g
    πz = model.back_grad_cache.πz
    return DiffOpt.lazy_combination(-, πz, g, length(g), n .+ i)
end
function DiffOpt._get_dA(model::Model, ci::MOI.ConstraintIndex{<:MOI.AbstractVectorFunction})
    i = MOI.Utilities.rows(model.model.constraints, ci) # vector
    # i = ci.value
    n = length(model.x) # columns in A
    m = length(model.y) # lines in A
    # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
    g = model.back_grad_cache.g
    πz = model.back_grad_cache.πz
    return DiffOpt.lazy_combination(-, g, πz, i, n .+ (1:n))
end

end


text:
function calculUtil(var::Int64, x::Vector{Int64}, prob::Problem)
	return prob.C[var]^3 / (prob.occurConst[var])
end

function findVariable(x::Vector{Int64}, prob::Problem)::Int64
	varMax::Int64 = -1
	utilMax::Float64 = -Inf

	for var = 1:prob.nbVar
		if x[var] == -1
			util = calculUtil(var, x, prob)
			if util > utilMax
				varMax = var
				utilMax = util
			end
		end
	end
	return varMax
end

function glouton(prob::Problem)::Solution

	estLibre::BitArray = trues(prob.nbVar)

	x::Vector{Int64} = Vector{Int64}(undef, prob.nbVar)
	for ind = 1:prob.nbVar
		x[ind] = -1
	end

	z::Int64 = 0

	nbVarLibre::Int64 = prob.nbVar
	while nbVarLibre > 0
		varAdd = findVariable(x, prob)

		x[varAdd] = 1
		estLibre[varAdd] = false
		nbVarLibre -= 1

		z += prob.C[varAdd]

		for varBloque = 1:prob.nbVar
			if prob.listOccur[varAdd, varBloque] == 1 && x[varBloque] == -1
				x[varBloque] = 0
				estLibre[varBloque] = false
				nbVarLibre -= 1
			end
		end
	end
	return Solution(x, z, estLibre)

end


text:
# Samplers for Categorical distribution

abstract AbstractCategoricalSampler

include("drawtable.jl")
include("alias.jl")
include("huffman.jl")