text:
# This defines a macro `mlj_model` which is a simpler version than the
# @sk_model macro defined to help import sklearn models.
# The difference is that the `mlj_model` macro only defines the constructor and the `clean!`
# and does not automatically define the `fit` and `predict` methods
#
# NOTE: it does NOT handle parametric types (yet).

"""
_process_model_def(ex)

Take an expression defining a model (`mutable struct Model ...`) and unpack key elements for
further processing:

- Model name (`modelname`)
- Names of parameters (`params`)
- Default values (`defaults`)
- Constraints (`constraints`)
"""
function _process_model_def(ex)
    defaults    = Dict{Symbol,Any}()
    constraints = Dict{Symbol,Any}()
    modelname   = ex.args[2] isa Symbol ? ex.args[2] : ex.args[2].args[1]
    params      = Symbol[]

    # inspect all lines which may define parameters, retrieve their names,
    # default values and constraints on values that can be given to them
    for i in 1:length(ex.args[3].args)
        # retrieve meaningful lines
        line = ex.args[3].args[i]
        line isa LineNumberNode && continue

        # line without information (e.g. just a name "a")
        if line isa Symbol
            param = line
            push!(params, param)
            defaults[param] = missing
        else
            # A meaningful line will look like
            #   line.args[1] = line.args[2]
            #
            # where line.args[1] will either be just `name`  or `name::Type`
            # and   line.args[2] will either be just `value` or `value::constraint`
            # ---------------------------------------------------------
            # 1. decompose `line.args[1]` appropriately (name and type)
            if line.args[1] isa Symbol # case :a
                param = line.args[1]
                type  = length(line.args) > 1 ? line.args[2] : :Any
            else                       # case :(a::Int)
                param, type = line.args[1].args[1:2] # (:a, Int)
            end
            push!(params, param)
            # ------------------------------------------------------------------
            # 2. decompose `line.args[2]` appropriately (values and constraints)
            if line.head == :(=) # assignment for default
                default = line.args[2]
                # if a constraint is given (value::constraint)
                if default isa Expr && length(default.args) > 1
                    constraints[param] = default.args[2]
                    # now discard the constraint to keep only the value
                    default = default.args[1]
                end
                defaults[param]    = default       # this will be a value not an expr
                ex.args[3].args[i] = line.args[1]  # name or name::Type (for the constructor)
            else
                # these are simple heuristics when no default value is given for the
                # field but an "obvious" one can be provided implicitly (ideally this should
                # not be used as it's not very clear that the intention matches the usage)
                eff_type = eval(type)
                if eff_type <: Number
                    defaults[param] = zero(eff_type)
                elseif eff_type <: AbstractString
                    defaults[param] = ""
                elseif eff_type == Any         # e.g. Any or no type given
                    defaults[param] = missing
                elseif eff_type >: Nothing     # e.g. Union{Nothing, ...}
                    defaults[param] = nothing
                elseif eff_type >: Missing     # e.g. Union{Missing, ...} (unlikely)
                    defaults[param] = missing
                else
                    @error "A default value for parameter '$param' (type '$type') must be given"
                end
            end
        end
    end
    return ex, modelname, params, defaults, constraints
end

"""
_unpack!(ex, rep)

Internal function to allow to read a constraint given after a default value for a parameter
and transform it in an executable condition (which is returned to be executed later).
For instance if we have

    alpha::Int = 0.5::(arg > 0.0)

Then it would transform the `(arg > 0.0)` in `(alpha > 0.0)` which is executable.
"""
function _unpack!(ex::Expr, rep)
    for i in eachindex(ex.args)
        if ex.args[i] ∈ (:_, :arg)
            ex.args[i] = rep
        end
        _unpack!(ex.args[i], rep)
    end
    return ex
end
_unpack!(ex, _) = ex # when it's been unpacked, it's not an expression anymore

"""
_model_constructor(modelname, params, defaults)

Build the expression of the keyword constructor associated with a model definition.
When the constructor is called, the `clean!` function is called as well to check that
parameter assignments are valid.
"""
function _model_constructor(modelname, params, defaults)
    Expr(:function, Expr(:call, modelname, Expr(:parameters, (Expr(:kw, p, defaults[p]) for p in params)...)),
        # body of the function
        Expr(:block,
             Expr(:(=), :model, Expr(:call, :new, params...)),
             :(message = clean!(model)),
			 :(isempty(message) || @warn message),
			 :(return model)
			 )
	 	)
end

"""
_model_cleaner(modelname, defaults, constraints)

Build the expression of the cleaner associated with the constraints specified in a model def.
"""
function _model_cleaner(modelname, defaults, constraints)
    Expr(:function, :(clean!(model::$modelname)),
        # body of the function
        Expr(:block,
             :(warning = ""),
             # condition and action for each constraint
             # each parameter is given as field::Type = default::constraint
             # here we recuperate the constraint and express it as an if statement
             # for instance if we had
             #     alpha::Real = 0.0::(arg > 0.0)
             # this would become
             #     if !(alpha > 0.0)
             (Expr(:if, Expr(:call, :!, _unpack!(constr, :(model.$param))),
                   # action of the constraint is violated:
                   # add a message and use default for the parameter
                   Expr(:block,
                        :(warning *= $("Constraint `$constr` failed; using default: $param=$(defaults[param]).")),
                        :(model.$param = $(defaults[param]))
                        )
                   ) for (param, constr) in constraints)...,
             # return full message
             :(return warning)
            )
        )
end

"""
mlj_model

Macro to help define MLJ models with constraints on the default parameters, this can be seen as
a tweaked version of the `@with_kw` macro from `Parameters`.
"""
macro mlj_model(ex)
    ex, modelname, params, defaults, constraints = _process_model_def(ex)
	# keyword constructor
    const_ex = _model_constructor(modelname, params, defaults)
	# associate the constructor with the definition of the struct
    push!(ex.args[3].args, const_ex)
	# cleaner
    clean_ex = _model_cleaner(modelname, defaults, constraints)
    esc(
        quote
            Base.@__doc__ $ex
            export $modelname
            $ex
            $clean_ex
        end
    )
end


text:
# ---
# title: 1325. Delete Leaves With a Given Value
# id: problem1325
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/delete-leaves-with-a-given-value/description/>
# hidden: true
# ---
#
# Given a binary tree `root` and an integer `target`, delete all the  **leaf
# nodes**  with value `target`.
#
# Note that once you delete a leaf node with value `target` **,  **if it's
# parent node becomes a leaf node and has the value `target`, it should also be
# deleted (you need to continue doing that until you can't).
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/09/sample_1_1684.png)**
#
#
#
#     Input: root = [1,2,3,2,null,2,4], target = 2
#     Output: [1,null,3,null,4]
#     Explanation: Leaf nodes in green with value (target = 2) are removed (Picture in left).
#     After removing, new nodes become leaf nodes with value (target = 2) (Picture in center).
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/09/sample_2_1684.png)**
#
#
#
#     Input: root = [1,3,3,3,2], target = 3
#     Output: [1,3,null,null,2]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/01/15/sample_3_1684.png)**
#
#
#
#     Input: root = [1,2,null,2,null,2], target = 2
#     Output: [1]
#     Explanation: Leaf nodes in green with value (target = 2) are removed at each step.
#
#
# **Example 4:**
#
#
#
#     Input: root = [1,1,1], target = 1
#     Output: []
#
#
# **Example 5:**
#
#
#
#     Input: root = [1,2,3], target = 1
#     Output: [1,2,3]
#
#
#
#
# **Constraints:**
#
#   * `1 <= target <= 1000`
#   * The given binary tree will have between `1` and `3000` nodes.
#   * Each node's value is between `[1, 1000]`.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
"""
Return snaps stored in `SeisForwExpt` after using `mod!`.

# Arguments

* `iss::Int64=1` : supersource index
"""
function Base.getindex(pa::Param, s::Symbol, iss::Int=1)
	@assert iss≤pa.c.acqgeom[1].nss
	@assert s in [:snaps]
	if(s==:snaps)
		nzd,nxd=length.(pa.c.model.mgrid)
		snaps_all=zeros(nzd,nxd,length(pa.c.itsnaps))
		ip=findall(in.(iss,pa.sschunks))[1]
		issp=findall(x->x==iss,pa.sschunks[ip])[]
		@sync begin
			p=procs(pa.p)[ip]
			@sync remotecall_wait(p) do
				pap=localpart(pa.p)
				for i in eachindex(snaps_all)
					snaps_all[i]=pap.ss[issp].snaps[i]
				end
			end
		end
		return snaps_all
	end
end


text:
# This file was generated, do not modify it. # hide
using MLJ, PrettyPrinting, DataFrames, Statistics

X, y = @load_reduced_ames
X = DataFrame(X)
@show size(X)
first(X, 3) |> pretty


text:
using LinearAlgebra
BLAS.set_num_threads(1)
using VecchiaFactorization
import VecchiaFactorization as VF
using BlockArrays
using Random
using Test
using LBblocks

@testset "vecchia_approx.jl and sparse_matrix_show.jl" begin

    n = 320

    Σ = @sblock let θtru = 1.5, n
        K = (x,y,θtru) -> exp(- θtru * abs(x - y) ^ 0.8 )
        x = range(0,4,n)
        K.(x, x', θtru)
    end

    block_sizes = [100, 150, 20, 50]
    perm = sortperm(@. sin((1:n)*2*π/n))

    R1, M1, P1 = VF.R_M_P(Σ, block_sizes, perm)
    R2, M2, P2 = VF.R_M_P(Σ[perm, perm], block_sizes)
    R3, M3, P3 = VF.R_M_P((i,j) -> Σ[i,j], block_sizes, perm)

    @test VF.sparse(R1) ≈ VF.sparse(R2) ≈ VF.sparse(R3)
    @test VF.sparse(M1) ≈ VF.sparse(M2) ≈ VF.sparse(M3)
    @test VF.sparse(P1) ≈ VF.sparse(P3)

    @test VF.Matrix(R1) ≈ VF.Matrix(R2) ≈ VF.Matrix(R3)
    @test VF.Matrix(M1) ≈ VF.Matrix(M2) ≈ VF.Matrix(M3)
    @test VF.Matrix(P1) ≈ VF.Matrix(P3)

    @test VF.sparse(R1') ≈ VF.sparse(R1)';
    @test VF.sparse(M1') ≈ VF.sparse(M1)';
    @test VF.sparse(P1') ≈ VF.sparse(P1)';

    VΣ = VF.vecchia(Σ, block_sizes, perm)
    VF.sparse(inv(VΣ))
    VF.sparse(inv(adjoint(VΣ)))
    VF.sparse(adjoint(inv(VΣ)))

    VF.Matrix(VΣ)
    VF.Matrix(inv(VΣ))
    VF.Matrix(adjoint(VΣ))
    VF.Matrix(inv(adjoint(VΣ)))
    VF.Matrix(adjoint(inv(VΣ)))

end

@testset "Ridiagonal and Qidiagonal" begin

    @sblock let

        bs = [2, 4, 3]
        @test VF.sizes_from_blocksides(Ridiagonal, bs) == [(4,2), (3,4)]
        @test VF.sizes_from_blocksides(Qidiagonal, bs) == [(2,4), (4,3)]
        @test VF.sizes_from_blocksides(Midiagonal, bs) == [(2,2), (4,4), (3,3)]

        Ri = rand(Ridiagonal{Float64}, bs)
        Qi = rand(Qidiagonal{Float64}, bs)
        Mi = rand(Midiagonal{Float64}, bs)

        @test Ri isa Ridiagonal
        @test Qi isa Qidiagonal
        @test Mi isa Midiagonal

        @test VF.diag_block_dlengths(Ri) == [2,4,3]
        @test VF.diag_block_dlengths(Qi) == [2,4,3]
        @test VF.diag_block_dlengths(Mi) == [2,4,3]

        @test size(Ri) == (9,9)
        @test size(Qi) == (9,9)
        @test size(Mi) == (9,9)

        @test (size(Ri,1), size(Ri,2)) == size(Ri)
        @test (size(Qi,1), size(Qi,2)) == size(Qi)
        @test (size(Mi,1), size(Mi,2)) == size(Mi)

    end

    @sblock let
        ## bs = vcat(fill(60,7), fill(80,10))
        ## bs = fill(80,20)
        bs = [2, 4, 3, 10]
        ## bs = fill(10,20)
        ## bs = fill(10,200)

        Ri = randn(Ridiagonal{Float64}, bs)
        Qi = randn(Qidiagonal{Float64}, bs)
        Mi = randn(Midiagonal{Float64}, bs)

        v = rand(Float64, sum(bs))

        @inferred Ri * v
        @inferred Ri' * v
        @inferred Ri \ v
        @test (Ri \ v) ≈ ldiv!(Ri, copy(v)) rtol = 1e-10
        @inferred Qi * v
        @inferred Qi' * v
        @inferred Qi \ v
        @test (Qi \ v) ≈ ldiv!(Qi, copy(v)) rtol = 1e-10

        ## -----
        x = 1:size(Ri,1)
        v = zeros(Float64, length(x));
        for i = 1:2  # it is strange how unstable this inversion is
            v[rand(x)] = 1
        end
        ## -----
        #
        ## τ = 10; v = sin.(τ .* 2 .* π .* x ./ x[end])
        ## -----
        ## v = rand(Float64, sum(bs))
        ## -----
        v1 = Ri \ Ri * v
        v2 = Qi \ Qi * v
        w1 = Ri * (Ri \ v)
        w2 = Qi * (Qi \ v)
        @test v ≈ v1 rtol=1e-5
        @test v ≈ v2 rtol=1e-5
        @test v ≈ w1 rtol=1e-5
        @test v ≈ w2 rtol=1e-5

    end

    #=
    using BenchmarkTools

    bs = fill(50,20)
    Ri = randn(Ridiagonal{Float64}, bs)
    Qi = randn(Qidiagonal{Float64}, bs)
    Mi = randn(Midiagonal{Float64}, bs)
    v  = randn(size(Ri,1))
    w  = randn(size(Ri,1))
    M = randn(size(Ri,1), size(Ri,1))
    @benchmark mul!(w, Mi, v, true, false)  # 9 μs
    @benchmark Mi * v  # 9 μs
    @benchmark Ri * v  # 9 μs
    @benchmark Qi * v  # 9 μs

    @benchmark Mi' * v # 9 μs
    @benchmark Ri' * v # 9 μs
    @benchmark Qi' * v # 9 μs

    @benchmark Mi \ v  # 716.999 μs
    @benchmark Ri \ v  # 9 μs
    @benchmark Qi \ v  # 9 μs

    @benchmark Mi' \ v # 734 μs
    @benchmark Ri' \ v # 9 μs
    @benchmark Qi' \ v # 9 μs

    @benchmark M * v # 112 μs
    =#

end

@testset "Vecchia, InvVecchia, VecchiaPivoted, InvVecchiaPivoted" begin

    K(x,y,θ) = exp(- θ * abs(x - y) ^ 0.8 )
    θtru    = 1.5
    nblocks, blocksz = 6, 150
    bsn    = fill(blocksz, nblocks)
    n      = sum(bsn)
    x      = vcat(range(0,4,length=n-100), 4*rand(100))
    prm    = sortperm(x)
    Σ      = PseudoBlockArray(K.(x, x', θtru), bsn, bsn)

    V = Vecchia(;
        diag_blocks=[Σ[Block(i,i)] for i = 1:nblocks],
        subdiag_blocks=[Σ[Block(i+1,i)] for i = 1:nblocks-1],
        )
    V′  = inv(V.R) * V.M * inv(V.R)'
    iV  = inv(V)
    iV′ = inv(V′)

    Vᴾ     = VecchiaPivoted(Vecchia(Σ[prm,prm], bsn), prm)
    Vᴾ′    = Piv(prm)' * inv(Vᴾ.R) * Vᴾ.M * inv(Vᴾ.R)' * Piv(prm)
    iVᴾ    = inv(Vᴾ)
    iVᴾ′   = inv(Vᴾ′)

    matV   = Matrix(V)
    matVᴾ  = Matrix(Vᴾ)
    matiV  = Matrix(iV)
    matiVᴾ = Matrix(iVᴾ)

    v    = randn(size(V,1))

    @test V  * v ≈ matV  * v rtol=1e-5
    @test V  * v ≈ V′    * v rtol=1e-5
    @test Vᴾ * v ≈ matVᴾ * v rtol=1e-5
    @test Vᴾ * v ≈ Vᴾ′   * v rtol=1e-5

    @test iV  * v ≈ matiV  * v rtol=1e-5
    @test iV  * v ≈ iV′    * v rtol=1e-5
    @test iVᴾ * v ≈ matiVᴾ * v rtol=1e-5
    @test iVᴾ * v ≈ iVᴾ′   * v rtol=1e-5

    #=
    using BenchmarkTools

    v  = randn(size(V,1))

    @benchmark V   * v   # 50  μs
    @benchmark V′  * v   # 100 μs
    @benchmark Vᴾ  * v   # 50  μs
    @benchmark Vᴾ′ * v   # 100 μs
    @benchmark iV   * v  # 50  μs
    @benchmark iV′  * v  # 50  μs
    @benchmark iVᴾ  * v  # 50  μs
    @benchmark iVᴾ′ * v  # 50  μs

    R = V.R
    @benchmark R * v          # 15 μs
    @benchmark $(R') * v      # 15 μs
    @benchmark $(inv(R)) * v  # 15 μs
    @benchmark $(inv(R)') * v # 15 μs

    Q = Qidiagonal(map(x->copy(x'), V.R.data))
    @benchmark Q * v          # 15 μs
    @benchmark $(Q') * v      # 15 μs
    @benchmark $(inv(Q)) * v  # 15 μs
    @benchmark $(inv(Q)') * v # 15 ms

    M = V.M
    @benchmark M * v          # 19 μs
    @benchmark $(M') * v      # 20 μs
    @benchmark $(inv(M)) * v  # 20 μs
    @benchmark $(inv(M)') * v # 20 μs
    @benchmark mul!(v, M, $(copy(v))) # 17 μs

    @benchmark Σ * v    # 83 μs
    =#

    #=
    using PyPlot
    fig, ax = subplots(2)
    ax[1].semilogy(eigen(Symmetric(matV)).values,label="matV eigen")
    ax[1].semilogy(eigen(Symmetric(Matrix(Σ))).values,"--",label="Σ eigen")
    ax[1].legend()
    ax[2].semilogy(eigen(Symmetric(matVᴾ)).values,label="matVᴾ eigen")
    ax[2].semilogy(eigen(Symmetric(Matrix(Σ))).values,"--",label="Σ eigen")
    ax[2].legend()
    =#

    #=
    using PyPlot
    fig, ax = subplots(2)
    ax[1].plot(Σv1,label="Σv1")
    ax[1].plot(Σv ,"--",label="Σv")
    ax[1].legend()
    ax[2].plot(Σv3,label="Σv3")
    ax[2].plot(Σv ,"--",label="Σv")
    ax[2].legend()
    =#

    #=
    using PyPlot

    @sblock let V = Vᴾ, Σ = Σ
        mV  = Matrix(V)
        imV = Matrix(inv(V))
        iΣ  = inv(Σ)

        fig,ax = subplots(2,3)
        mV     |> ax[1,1].imshow
        Σ      |> ax[1,2].imshow
        mV - Σ |> ax[1,3].imshow

        imV      .|> abs .|> log |> ax[2,1].imshow
        iΣ       .|> abs .|> log |> ax[2,2].imshow
        imV - iΣ .|> abs .|> log |> ax[2,3].imshow
    end

    @sblock let V
        matRᴴ = VF.Rᴴmat(V)
        matR  = VF.Rmat(V)
        matRᴴ .|> abs .|> log |> matshow; colorbar()
        matR  .|> abs .|> log |> matshow; colorbar()
    end

    =#

    #=
    using PyPlot
    fig, ax = subplots(2)
    ax[1].plot(w1,label="w1")
    ax[1].plot(w ,"--",label="w")
    ax[1].legend()
    ax[2].plot(w3,label="w3")
    ax[2].plot(w ,"--",label="w")
    ax[2].legend()
    =#

    ## invcholΣ = inv(cholesky(Hermitian(Σ,:L)).L)
    ## invcholV = VecchiaFactorization.inv_cholesky(V)
    ## invcholVᴾ = VecchiaFactorization.inv_cholesky(Vecchia(Vᴾ.R,Vᴾ.M, Vᴾ.bsds)) # [invperm(Vᴾ.piv), invperm(Vᴾ.piv)]

    #=
    using PyPlot
    v = randn(size(V,1))

    (invcholV \ v)[prm] |> plot
    (invcholVᴾ \ v[prm]) |> plot
    (invcholΣ \ v)[prm] |> plot
    =#

    #=
    using BenchmarkTools
    ## without threads ...
    @benchmark V * v   # 48 μs
    @benchmark iV * v  # 50 μs

    @benchmark $(Matrix(matV)) * v    # 93 μs
    @benchmark $(Matrix(imatV)) * v   # 90 μs

    @benchmark $matV * v    # 108 μs
    @benchmark $imatV * v   # 108 μs

    (@belapsed $(Matrix(matV)) * v) / (@belapsed V * v)     # 1.7
    Base.summarysize((Matrix(matV))) / Base.summarysize(V)  # 3.27
    =#

end

@testset "Vecchia approx inversion" begin

    K(x,y,θ) = exp(- θ * abs(x - y) ^ 0.8 )
    θtru    = 0.2
    nblocks, blocksz = 50, 50
    bsn    = fill(blocksz, nblocks)
    n      = sum(bsn)
    x      = sort(vcat(range(0,4,length=n-100), 4*rand(100)))

    Σ  = K.(x, x', θtru)
    iΣ = inv(Σ)
    V  = Vecchia(Σ,bsn)
    iV = inv(V)

    IVapx = hcat(map(x->iV*x, eachcol(Σ))...)

    #=
    using PyPlot
    fig, ax = subplots(2)
    ax[1].plot(real.(eigen(IVapx).values))
    ax[1].plot(zeros(size(IVapx,1)),":k")
    ax[2].plot(imag.(eigen(IVapx).values))
    ax[2].plot(zeros(size(IVapx,1)),":k")
    =#

end

## @testset "Inv" begin
##
##     @test VF.Inv(R') == VF.Inv(R)'
##
##
## end

@testset "util.jl" begin

    @test VF.block_split(20, 10) == [10,10]
    @test VF.block_split(20, 5) == [5,5,5,5]
    @test VF.block_split(20, 11) == [11,9]
    @test VF.block_split(11, 11) == [11]
    @test VF.block_split(11, 10) == [10,1]

end


text:
@testset "$(kind) operators" for (kind,Vfun) in [
    ("Linear", (R,fx) -> LinearOperator(R'*QuasiDiagonal(fx)*R, R)),
    ("Diagonal", (R,fx) -> DiagonalOperator(applied(*, R, R \ fx)))
]
    @testset "T = $(T)" for (T,o,γ) in [(Float64,x -> one(Float64),1.0),(ComplexF64,x -> one(ComplexF64),exp(im*π/4))]
        L = 10.0
        f = x -> γ * sin(2π*x/L)

        N = 100

        @testset "$(label)" for (label,R) in [
            ("Finite-differences", FiniteDifferences(N, L/(N+1))),
            ("Uniform, staggered finite-differences", StaggeredFiniteDifferences(N, L/(N-1/2))),
            ("Non-uniform, staggered finite-differences", StaggeredFiniteDifferences(0.1, 0.2, 0.01, L)),
            ("Implicit finite-differences", ImplicitFiniteDifferences(N, L/(N+1))),
            ("FE-DVR", FEDVR(range(0, stop=L, length=N), 4)),
            ("B-splines", BSpline(LinearKnotSet(5, 0, L, N)))
        ]
            x = axes(R,1)
            fx = f.(x)
            fc = R \ fx
            oc = R \ o.(x)

            V = Vfun(R, fx)
            m = size(R,2)
            @test size(V) == (m,m)
            @test size(V,1) == size(V,2) == m
            @test eltype(V) == T

            u = similar(oc)
            mul!(u, V, oc)
            @test u ≈ fc
        end
    end
end

@testset "Matricization of DiagonalOperators" begin
    L = 1.0

    N = 100

    @testset "$(label)" for (label,R) in [
        ("Finite-differences", FiniteDifferences(N, L/(N+1))),
        ("Uniform, staggered finite-differences", StaggeredFiniteDifferences(N, L/(N-1/2))),
        ("Non-uniform, staggered finite-differences", StaggeredFiniteDifferences(0.1, 0.2, 0.01, L)),
        ("Implicit finite-differences", ImplicitFiniteDifferences(N, L/(N+1))),
        ("FE-DVR", FEDVR(range(0, stop=L, length=N), 4)),
        ("B-splines", BSpline(LinearKnotSet(7, 0, L, N)))
    ]
        r = axes(R, 1)

        f = R \ exp.(-(r .- 0.5).^2/(2*0.04^2))
        g = R \ (exp.(-(r .- 0.5).^2/(2*0.04^2)) .* r)

        d = DiagonalOperator(applied(*, R, f))

        o = R \ (r -> 1).(r)
        ir = R \ identity.(r)

        m = Matrix(d)
        @test m isa (R isa BSpline ? BandedMatrix : Diagonal)
        f′ = LinearOperator(m, R)*o
        @test f ≈ f′

        n = Matrix(d, ir)
        @test n isa (R isa BSpline ? BandedMatrix : Diagonal)
        g′ = LinearOperator(n, R)*o
        @test g ≈ g′
    end
end

@testset "ShiftAndInvert" begin
    L = 10.0
    N = 100

    @testset "$(label)" for (label,R,tol) in [
        ("Finite-differences", FiniteDifferences(N, L/(N+1)), 0.004),
        ("Uniform, staggered finite-differences", StaggeredFiniteDifferences(N, L/(N-1/2)), 0.05),
        ("Non-uniform, staggered finite-differences", StaggeredFiniteDifferences(0.1, 0.2, 0.01, L), 0.05),
        ("Implicit finite-differences", ImplicitFiniteDifferences(N, L/(N+1)), 5e-6),
        ("FE-DVR", FEDVR(range(0, stop=L, length=N), 4)[:,2:end-1], 8e-11),
        ("B-splines", BSpline(LinearKnotSet(5, 0, L, N))[:,2:end-1], 5e-13)
    ]
        λ,ϕ,r̃,R,δλ = test_particle_in_a_box(R, L, 5)
        @test norm(δλ) < tol
    end
end


text:
using Base.Test
using Meganet
using LinearOperators

nImg = [8,10]
sK   = [3,3,4]

K = getConvCircKernel(Float64,nImg,sK)

@testset  "adjoint test" begin
theta = initTheta(K)
A     = getOp(K,theta);
v     = randn(nFeatIn(K))
w     = randn(nFeatIn(K))

t1    = dot(w,A*v)
t2    = dot(v,A'*w)
# println("adjointTest t1=$t1\t t2=$t2")
@test norm(t1-t2)/norm(t1) < 1e-10
end

@testset "derivative Test" begin
th  = initTheta(K);
dth = initTheta(K);
nex = 1;
Y  = randn(nFeatIn(K),nex)+nex;
Z  = randn(nFeatIn(K),nex)-nex;

t1 = vec(Z)'*vec(Jthetamv(K,dth,th,Y));
t2 = vec(dth)'*vec(JthetaTmv(K,Z,th,Y));
# println("derivativeTest t1=$t1\t t2=$t2")
@test norm(t1-t2)/norm(t2) < 1e-10
end


text:
"""
Something that is SSA-indexable.
"""
abstract type SSAIndexable end

struct SSADict{T}
    dict::Dictionary{ID,T}
    SSADict{T}() where {T} = new{T}(Dictionary{ID,T}())
    SSADict(pairs::AbstractVector{Pair{ID,T}}) where {T} = new{T}(Dictionary(pairs...))
    SSADict{T}(dict::Dictionary{ID,T}) where {T} = new{T}(dict)
    SSADict(dict::Dictionary{ID,T}) where {T} = SSADict{T}(dict)
end

SSADict(args::AbstractVector{<:SSAIndexable}) = SSADict(ID.(args) .=> args)
SSADict(args::Vararg) = SSADict(collect(args))
SSADict() = SSADict{Any}()

Base.convert(::Type{SSADict{T}}, dict::SSADict) where {T} = SSADict{T}(convert(Dictionary{ID,T}, dict.dict))

@forward SSADict.dict Base.getindex, Base.insert!, Dictionaries.set!, Base.get!, Base.get, Base.setindex!, Base.pop!, Base.first, Base.last, Base.broadcastable, Base.length, Base.iterate, Base.keys, Base.values, Base.haskey

Base.merge!(vec::SSADict, others::SSADict...) = merge!(vec.dict, getproperty.(others, :dict)...)


text:
##### calculate diffusivities of each individual
@kernel function calc_diffusion_kernel!(plank, rnd, κx, κy, κz, ΔT)
    i = @index(Global)
    @inbounds plank.x[i] = plank.x[i] + rnd.x[i] * √(1*κx*ΔT) * plank.ac[i]
    @inbounds plank.y[i] = plank.y[i] + rnd.y[i] * √(1*κy*ΔT) * plank.ac[i]
    # @inbounds plank.z[i] = plank.z[i] + rnd.z[i] * √(1*κz*ΔT) * plank.ac[i]
end
function calc_diffusion!(plank, rnd, κx, κy, κz, ΔT, arch::Architecture)
    kernel! = calc_diffusion_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, rnd, κx, κy, κz, ΔT)
    wait(device(arch), event)

    return nothing
end

function plankton_diffusion!(plank, rnd, κx, κy, κz, ΔT, g::AbstractGrid, arch::Architecture)
    ##### generate random numbers for grazing, mortality and division (0,1)
    randn!(rng_type(arch), rnd.x)
    randn!(rng_type(arch), rnd.y)
    # randn!(rng_type(arch), rnd.z)

    ##### shift the random numbers to (-1,1)
    @inbounds rnd.x .= rnd.x
    @inbounds rnd.y .= rnd.y
    # @inbounds rnd.z .= rnd.z .* 0.3

    ##### calculate diffusion
    calc_diffusion!(plank, rnd, κx, κy, κz, ΔT, arch)

    ##### keep individuals in the domain
    particle_boundaries!(plank, plank.ac, g, arch)
end

plankton_diffusion!(plank, rnd, κ, ΔT, g, arch) = plankton_diffusion!(plank, rnd, κ, κ, κ, ΔT, g, arch)


text:
using NNlib
import NNlib: softmax, ∇softmax, logsoftmax, ∇logsoftmax, conv, maxpool, meanpool

@grad softmax(xs) = softmax(xs), Δ -> (∇softmax(Δ, xs),)

@grad logsoftmax(xs) = logsoftmax(xs), Δ -> (∇logsoftmax(Δ, xs),)

@grad conv(x, w; kw...) =
  conv(x, w; kw...),
    Δ ->
      (NNlib.∇conv_data(Δ, x, w; kw...),
       NNlib.∇conv_filter(Δ, x, w; kw...))

@grad function maxpool(x, k; kw...)
  y = maxpool(x, k; kw...)
  y, Δ -> (NNlib.∇maxpool(Δ, y, x, k; kw...), nothing)
end

@grad function meanpool(x, k; kw...)
  y = meanpool(x, k; kw...)
  y, Δ -> (NNlib.∇meanpool(Δ, y, x, k; kw...), nothing)
end


text:
const PTRACE_TRACEME      = 0
const PTRACE_PEEKTEXT     = 1
const PTRACE_PEEKDATA     = 2
const PTRACE_PEEKUSER     = 3
const PTRACE_POKETEXT     = 4
const PTRACE_POKEDATA     = 5
const PTRACE_POKEUSER     = 6
const PTRACE_CONT         = 7
const PTRACE_KILL         = 8
const PTRACE_SINGLESTEP   = 9
const PTRACE_GETREGS      = 12
const PTRACE_SETREGS      = 13
const PTRACE_GETFPREGS    = 14
const PTRACE_SETFPREGS    = 15
const PTRACE_ATTACH       = 16
const PTRACE_DETACH       = 17
const PTRACE_GETFPXREGS   = 18
const PTRACE_SETFPXREGS   = 19
const PTRACE_SYSCALL      = 24
const PTRACE_SETOPTIONS   = 0x4200
const PTRACE_GETEVENTMSG  = 0x4201
const PTRACE_GETSIGINFO   = 0x4202
const PTRACE_SETSIGINFO   = 0x4203
const PTRACE_GETREGSET    = 0x4204
const PTRACE_SETREGSET    = 0x4205
const PTRACE_SEIZE        = 0x4206
const PTRACE_INTERRUPT    = 0x4207


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith
import REPL: LineEdit, REPLCompletions
import REPL

ans = nothing

"""
    finished(s)

Examine the buffer in the repl to see if the input is complete
"""
function finished(s)
    str = String(take!(copy(LineEdit.buffer(s))))
    if length(str) == 0
        return false
    elseif str[end] == ';' || str[end] == '$'
        return true
    else
        return false
    end
end

"""
    respond(repl, main)

Something dark and magic
"""
function respond(repl, main)
    (s, buf, ok) -> begin
        if !ok
            return REPL.transition(s, :abort)
        end
        input = String(take!(buf))
        if '%' in input
            tail = last(input)                              # chop the tail
            input = subst("$(ans)", "%", input[1:end-1])   # substitute
            input = string(input, tail)                     # add the tail
        end
        if !isempty(strip(input))
            try
                global ans = MExpr(input[1:end-1]) |> mcall
                REPL.reset(repl)
                if input[end] == ';'
                    REPL.print_response(repl, ans, nothing, true, Base.have_color)
                else
                    REPL.print_response(repl, nothing, nothing, true, Base.have_color)
                end
            catch err
                REPL.reset(repl)
                REPL.print_response(repl, err, catch_backtrace(), true, Base.have_color)
            end
        end
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main)
    end
end

"""
    MaximaCompletionProvider

Basic completion provider, just latex completions
"""
mutable struct MaximaCompletionProvider <: LineEdit.CompletionProvider
    r::REPL.LineEditREPL
end

function LineEdit.complete_line(c::MaximaCompletionProvider, s)
    buf = s.input_buffer
    partial = String(buf.data[1:buf.ptr-1])
    # complete latex
    full = LineEdit.input_string(s)
    ret, range, should_complete = REPLCompletions.bslash_completions(full, endof(partial))[2]

    if length(ret) > 0 && should_complete
        return ret, partial[range], true
    end

    return String[], 0:-1, false
end

function create_maxima_repl(repl, main)
    maxima_mode = LineEdit.Prompt(
        "maxima> ";
        prompt_prefix = Base.text_colors[:cyan],
        prompt_suffix = main.prompt_suffix,
        on_enter = finished,
        on_done = respond(repl, main),
        sticky = true)

    hp = main.hist
    hp.mode_mapping[:maxima] = maxima_mode
    maxima_mode.hist = hp
    maxima_mode.complete = MaximaCompletionProvider(repl)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, maxima_mode)

    mk = REPL.mode_keymap(main)

    b = Dict{Any,Any}[
        skeymap, mk, prefix_keymap, LineEdit.history_keymap,
        LineEdit.default_keymap, LineEdit.escape_defaults
    ]
    maxima_mode.keymap_dict = LineEdit.keymap(b)

    maxima_mode
end

function repl_init(repl)
    mirepl = isdefined(repl, :mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]
    maxima_mode = create_maxima_repl(mirepl, main_mode)
    push!(mirepl.interface.modes, maxima_mode)

    maxima_prompt_keymap = Dict{Any,Any}(
        ')' => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, maxima_mode) do
                    LineEdit.state(s, maxima_mode).input_buffer = buf
                end
            else
                if !isdefined(Main,:OhMyREPL)
                    LineEdit.edit_insert(s, ')')
                else
                    if Main.OhMyREPL.BracketInserter.AUTOMATIC_BRACKET_MATCH[] &&
                        !eof(LineEdit.buffer(s)) &&
                        Main.OhMyREPL.BracketInserter.peek(LineEdit.buffer(s)) == ')'
                        LineEdit.edit_move_right(LineEdit.buffer(s))
                    else
                        LineEdit.edit_insert(LineEdit.buffer(s), ')')
                    end
                    Main.OhMyREPL.Prompt.rewrite_with_ANSI(s)
                end
            end
        end
    )
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, maxima_prompt_keymap);
    nothing
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mc01wd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    DPMAX = 10
    P = Array{Float64,1}(undef, DPMAX+1)
    Q = Array{Float64,1}(undef, DPMAX+1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    DP = parse(BlasInt, vs[1])
    if ( DP<=-1 || DP>DPMAX )
        @error "Illegal DP=$DP"
    end

    vs = String[]
    _isz = DP+1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    P[1:_isz] .= parsex.(Float64, vs)

    vs = split(readline(f))
    U1 = parse(Float64, replace(vs[1],'D'=>'E'))
    U2 = parse(Float64, replace(vs[2],'D'=>'E'))
    close(f)
    # interp call 1

    INFO = SLICOT.mc01wd!(DP, P, U1, U2, Q)
    @test INFO == 0
    INFO == 0 || return

    println(io, "Q (power,coefft):")
    show(io, "text/plain", hcat(0:DP-2,Q[3:DP+1]))
    println(io)
    println(io, "R (power, coefft):")
    show(io, "text/plain", hcat(0:1,[Q[1]+Q[2]*U2,Q[2]]))
    println(io)
end # run_mc01wd()


text:
function main(input)

end

max_stars_repo_name: KrymoreAkbota/julia

max_issues_repo_name: KrymoreAkbota/julia

max_forks_repo_name: KrymoreAkbota/julia


text:
module Stages

import DiffEqBase: solve!, solve
using DiffEqBase: step!, init

using ProgressMeter: ProgressMeter, Progress

abstract type Stageable end
abstract type AbstractSource <: Stageable end
abstract type AbstractStage <: Stageable end

abstract type ExecPolicy end
struct NullPolicy <: ExecPolicy end
struct WithProgress{P} <: ExecPolicy
    progress_factory::P
    WithProgress(progress_factory::P) where P = new{P}(progress_factory)
end

get_exec_policy(progress::Real) =
    if progress > 0
        WithProgress((n, args...) -> Progress(n, progress, args...))
    else
        NullPolicy()
    end

pre_loop(::AbstractStage, ::NullPolicy) = nothing
post_loop(::AbstractStage, ::NullPolicy, ::Any) = nothing
post_step(::AbstractStage, ::NullPolicy, ::Any) = nothing
function remaining_steps end

"""
    is_finished(stage::Stageable) :: Bool

Ask if `stage`'s computation is done.
"""
function is_finished end

"""
    finish!(stage::Stageable) :: Stageable

Finish whatever the computation `stage` has to do.
"""
function finish!(stage::AbstractStage, policy::ExecPolicy = NullPolicy())
    policy_state = pre_loop(stage, policy)
    while ! is_finished(stage)
        step!(stage)
        policy_state = post_step(stage, policy, policy_state)
    end
    post_loop(stage, policy, policy_state)
    return stage
end

shortname(::T) where T = string(T.name.name)

pre_loop(stage::AbstractStage, policy::WithProgress) =
    policy.progress_factory(remaining_steps(stage), shortname(stage))

function post_step(::AbstractStage, ::WithProgress, prog::Progress)
    ProgressMeter.next!(prog)
    return prog
end

post_loop(::AbstractStage, ::WithProgress, prog::Progress) =
    ProgressMeter.finish!(prog)

function finish_if_not!(stage::Stageable, policy::ExecPolicy = NullPolicy())
    is_finished(stage) || finish!(stage, policy)
    return stage
end

finish!(::AbstractSource, ::ExecPolicy) = nothing
is_finished(::AbstractSource) = true

# Some "mix-in" methods for generating default `is_finished`:
stage_index(stage::AbstractStage) = stage.i
is_finished(stage::AbstractStage) = stage_index(stage) >= length(stage)
remaining_steps(stage::AbstractStage) = length(stage) - stage_index(stage)
# TODO: maybe this shouldn't be defined for generic stage.  Maybe
# define AbstractIndexedStage?

# Utility method. Define it and call it from step! for type-based
# recording injection.
"""
    record!(stage, Val{key})

Record some info (labeled by `key`) into `stage`.
"""
record!(::AbstractStage, ::Any) = nothing
# See: [[./clv/core_stages.jl::function record!]]

# Default iterator interface ("mix-in")
current_result(stage::AbstractStage) = nothing

Base.start(stage::AbstractStage) = nothing
Base.done(stage::AbstractStage, _state) = is_finished(stage)
@inline function Base.next(stage::AbstractStage, _state)
    step!(stage)
    return (current_result(stage), nothing)
end

struct StageIterator
    source::Stageable
    stage_types::Vector{Type{<: AbstractStage}}
    args::Tuple
end

struct StageState
    i::Int
    stage::Stageable
end

Base.start(iter::StageIterator) = StageState(1, iter.source)

function Base.next(iter::StageIterator, state::StageState)
    next_state = StageState(
        state.i + 1,
        iter.stage_types[state.i](
            finish_if_not!(state.stage),
            iter.args...))
    return (next_state.stage, next_state)
end

Base.done(iter::StageIterator, state::StageState) =
    length(iter.stage_types) < state.i

is_reachable(iter::StageIterator, stage_type::Type, i::Int = 1) =
    any(t <: stage_type for t in iter.stage_types[i:end])

is_reachable(iter::StageIterator, stage_type::Type, state::StageState) =
    is_reachable(iter, stage_type, state.i)

"""
Stateful solver based on [`StageIterator`](@ref).
"""
mutable struct StagedSolver{P, S}
    prob::P
    iter::StageIterator
    state::StageState
    sol::S

    function StagedSolver(prob::P, sol::S,
                          stage_types::AbstractVector, args::Tuple,
                          ) where {P, S}
        iter = StageIterator(
            prob,
            stage_types,
            args,  # additional arguments to each `stage_types`
        )
        state = start(iter)
        return new{P, S}(prob, iter, state, sol)
    end
end

function advance!(solver::StagedSolver, policy::ExecPolicy = NullPolicy())
    if done(solver.iter, solver.state)
        return nothing
    end
    finish_if_not!(solver.state.stage, policy)
    _, solver.state = next(solver.iter, solver.state)
    return solver.state.stage
end

function goto_last!(solver::StagedSolver, policy::ExecPolicy = NullPolicy())
    while advance!(solver, policy) !== nothing end
    return solver.state.stage
end

function solve!(solver::StagedSolver; progress=-1)
    if done(solver.iter, solver.state) && is_finished(solver.state.stage)
        error("No further computation is required.")
        # Should it be just a no-op?
    end
    policy = get_exec_policy(progress)
    goto_last!(solver, policy)
    finish_if_not!(solver.state.stage, policy)
    return solver
end

solve(prob::AbstractSource; progress=-1, kwargs...) =
    solve!(init(prob; kwargs...); progress=progress).sol
# This requires `init(prob)` to be defined elsewhere and to return a
# StagedSolver.  Not sure defining `solve` by default is a good idea.
# But repeating this definition all over the places is also not good.

"""
    goto!(solver::StagedSolver, stage_type::Type{T}) :: T

Advance the `solver` up to the stage of type `stage_type` and return
it.
"""
function goto!(solver::StagedSolver, stage_type::Type; progress=-1)
    @assert is_reachable(solver.iter, stage_type, solver.state)
    if solver.state.stage isa stage_type
        return solver.state.stage
    end
    policy = get_exec_policy(progress)
    while true
        stage = advance!(solver, policy)
        if stage isa stage_type
            return stage
        end
    end
    error("Unreachable!")
end

end


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

function GitRebase(repo::GitRepo, branch::GitAnnotated, upstream::GitAnnotated;
                   onto::Union{GitAnnotated, Nothing}=nothing,
                   opts::RebaseOptions = RebaseOptions())
    rebase_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_rebase_init, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
                   Ptr{Cvoid}, Ptr{RebaseOptions}),
                   rebase_ptr_ptr, repo.ptr, branch.ptr, upstream.ptr,
                   onto === nothing ? C_NULL : onto.ptr, Ref(opts))
    return GitRebase(repo, rebase_ptr_ptr[])
end

function Base.count(rb::GitRebase)
    return ccall((:git_rebase_operation_entrycount, :libgit2), Csize_t, (Ptr{Cvoid},), rb.ptr)
end

"""
    current(rb::GitRebase) -> Csize_t

Return the index of the current [`RebaseOperation`](@ref). If no operation has
yet been applied (because the [`GitRebase`](@ref) has been constructed but `next`
has not yet been called or iteration over `rb` has not yet begun), return
`GIT_REBASE_NO_OPERATION`, which is equal to `typemax(Csize_t)`.
"""
function current(rb::GitRebase)
    return ccall((:git_rebase_operation_current, :libgit2), Csize_t, (Ptr{Cvoid},), rb.ptr)
end

function Base.getindex(rb::GitRebase, i::Integer)
    if !(1 <= i <= count(rb))
        throw(BoundsError(rb, (i,)))
    end
    GC.@preserve rb begin
        rb_op_ptr = ccall((:git_rebase_operation_byindex, :libgit2),
                          Ptr{RebaseOperation},
                          (Ptr{Cvoid}, Csize_t), rb.ptr, i-1)
        rb_op = unsafe_load(rb_op_ptr)
    end
    return rb_op
end

function Base.next(rb::GitRebase)
    rb_op_ptr_ptr = Ref{Ptr{RebaseOperation}}(C_NULL)
    GC.@preserve rb begin
        try
            @check ccall((:git_rebase_next, :libgit2), Cint,
                          (Ptr{Ptr{RebaseOperation}}, Ptr{Cvoid}),
                           rb_op_ptr_ptr, rb.ptr)
        catch err
            err.code == Error.ITEROVER && return nothing
            rethrow(err)
        end
        rb_op_ptr = unsafe_load(rb_op_ptr_ptr[])
    end
    return rb_op_ptr
end

function Base.show(io::IO, rb::GitRebase)
    println(io, "GitRebase:")
    println(io, "Number: ", count(rb))
    println(io, "Currently performing operation: ", current(rb)+1)
end

"""
    LibGit2.commit(rb::GitRebase, sig::GitSignature)

Commit the current patch to the rebase `rb`, using `sig` as the committer. Is silent if
the commit has already been applied.
"""
function commit(rb::GitRebase, sig::GitSignature)
    oid_ptr = Ref(GitHash())
    try
        @check ccall((:git_rebase_commit, :libgit2), Error.Code,
                     (Ptr{GitHash}, Ptr{Cvoid}, Ptr{SignatureStruct}, Ptr{SignatureStruct}, Ptr{UInt8}, Ptr{UInt8}),
                      oid_ptr, rb.ptr, C_NULL, sig.ptr, C_NULL, C_NULL)
    catch err
        # TODO: return current HEAD instead
        err.code == Error.EAPPLIED && return nothing
        rethrow(err)
    end
    return oid_ptr[]
end

"""
    abort(rb::GitRebase) -> Csize_t

Cancel the in-progress rebase, undoing all changes made so far and returning
the parent repository of `rb` and its working directory to their state before
the rebase was initiated. Return `0` if the abort is successful,
`LibGit2.Error.ENOTFOUND` if no rebase is in progress (for example, if the
rebase had completed), and `-1` for other errors.
"""
function abort(rb::GitRebase)
    return ccall((:git_rebase_abort, :libgit2), Csize_t,
                      (Ptr{Cvoid},), rb.ptr)
end

"""
    finish(rb::GitRebase, sig::GitSignature) -> Csize_t

Complete the rebase described by `rb`. `sig` is a [`GitSignature`](@ref)
to specify the identity of the user finishing the rebase. Return `0` if the
rebase finishes successfully, `-1` if there is an error.
"""
function finish(rb::GitRebase, sig::GitSignature)
    return ccall((:git_rebase_finish, :libgit2), Csize_t,
                  (Ptr{Cvoid}, Ptr{SignatureStruct}),
                   rb.ptr, sig.ptr)
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
##################################################################################
#### DEFINTION OF NEW VARIABLES FOR STORAGE INVESTMENTS ACCODING TO FlexPlan MODEL
##################################################################################

function variable_absorbed_energy(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    e_abs = _PM.var(pm, nw)[:e_abs] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :storage_bounded_absorption)], base_name="$(nw)_e_abs",
    start = 0)

    if bounded
        for (s, storage) in _PM.ref(pm, nw, :storage_bounded_absorption)
            JuMP.set_lower_bound(e_abs[s],  0)
            JuMP.set_upper_bound(e_abs[s],  storage["max_energy_absorption"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :storage, :e_abs, _PM.ids(pm, nw, :storage_bounded_absorption), e_abs)
end

function variable_absorbed_energy_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    e_abs = _PM.var(pm, nw)[:e_abs_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :ne_storage_bounded_absorption)], base_name="$(nw)_e_abs_ne",
    start = 0)

    if bounded
        for (s, storage) in _PM.ref(pm, nw, :ne_storage_bounded_absorption)
            JuMP.set_lower_bound(e_abs[s],  0)
            JuMP.set_upper_bound(e_abs[s],  storage["max_energy_absorption"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :e_abs_ne, _PM.ids(pm, nw, :ne_storage_bounded_absorption), e_abs)
end

function variable_storage_power_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_storage_power_real_ne(pm; kwargs...)
    variable_storage_power_imaginary_ne(pm; kwargs...)
    variable_storage_power_control_imaginary_ne(pm; kwargs...)
    variable_storage_current_ne(pm; kwargs...)
    variable_storage_energy_ne(pm; kwargs...)
    variable_storage_charge_ne(pm; kwargs...)
    variable_storage_discharge_ne(pm; kwargs...)
    variable_storage_indicator(pm; kwargs..., relax=true)
    variable_storage_investment(pm; kwargs...)
end

function variable_storage_power_real_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    ps = _PM.var(pm, nw)[:ps_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_ps_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "ps_start")
    )

    if bounded
        inj_lb, inj_ub = _PM.ref_calc_storage_injection_bounds(_PM.ref(pm, nw, :ne_storage), _PM.ref(pm, nw, :bus))

        for i in _PM.ids(pm, nw, :ne_storage)
            if !isinf(inj_lb[i])
                JuMP.set_lower_bound(ps[i], inj_lb[i])
            end
            if !isinf(inj_ub[i])
                JuMP.set_upper_bound(ps[i], inj_ub[i])
            end
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :ps_ne, _PM.ids(pm, nw, :ne_storage), ps)
end

function variable_storage_power_imaginary_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    qs = _PM.var(pm, nw)[:qs_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_qs_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "qs_start")
    )

    if bounded
        inj_lb, inj_ub = _PM.ref_calc_storage_injection_bounds(_PM.ref(pm, nw, :ne_storage), _PM.ref(pm, nw, :bus))

        for (i, storage) in _PM.ref(pm, nw, :ne_storage)
            JuMP.set_lower_bound(qs[i], max(inj_lb[i], storage["qmin"]))
            JuMP.set_upper_bound(qs[i], min(inj_ub[i], storage["qmax"]))
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :qs_ne, _PM.ids(pm, nw, :ne_storage), qs)
end

"apo models ignore reactive power flows"
function variable_storage_power_imaginary_ne(pm::_PM.AbstractActivePowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    report && _IM.sol_component_fixed(pm, _PM.pm_it_sym, nw, :ne_storage, :qs_ne, _PM.ids(pm, nw, :ne_storage), NaN)
end

function variable_storage_power_control_imaginary_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    qsc = _PM.var(pm, nw)[:qsc_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_qsc_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "qsc_start")
    )

    if bounded
        inj_lb, inj_ub = _PM.ref_calc_storage_injection_bounds(_PM.ref(pm, nw, :ne_storage), _PM.ref(pm, nw, :bus))

        for (i,storage) in _PM.ref(pm, nw, :ne_storage)

            if !isinf(inj_lb[i]) || haskey(storage, "qmin")
                JuMP.set_lower_bound(qsc[i], max(inj_lb[i], get(storage, "qmin", -Inf)))
            end
            if !isinf(inj_ub[i]) || haskey(storage, "qmax")
                JuMP.set_upper_bound(qsc[i], min(inj_ub[i], get(storage, "qmax",  Inf)))
            end
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :qsc_ne, _PM.ids(pm, nw, :ne_storage), qsc)
end

"apo models ignore reactive power flows"
function variable_storage_power_control_imaginary_ne(pm::_PM.AbstractActivePowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    report && _IM.sol_component_fixed(pm, _PM.pm_it_sym, nw, :ne_storage, :qsc_ne, _PM.ids(pm, nw, :ne_storage), NaN)
end

"do nothing by default but some formulations require this"
function variable_storage_current_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
end

function variable_storage_energy_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    se = _PM.var(pm, nw)[:se_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_se_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "se_start", 1)
    )

    if bounded
        for (i, storage) in _PM.ref(pm, nw, :ne_storage)
            JuMP.set_lower_bound(se[i], 0)
            JuMP.set_upper_bound(se[i], storage["energy_rating"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :se_ne, _PM.ids(pm, nw, :ne_storage), se)
end

function variable_storage_charge_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    sc = _PM.var(pm, nw)[:sc_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_sc_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "sc_start", 1)
    )

    if bounded
        for (i, storage) in _PM.ref(pm, nw, :ne_storage)
            JuMP.set_lower_bound(sc[i], 0)
            JuMP.set_upper_bound(sc[i], storage["charge_rating"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :sc_ne, _PM.ids(pm, nw, :ne_storage), sc)
end

function variable_storage_discharge_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    sd = _PM.var(pm, nw)[:sd_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_sd_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_storage, i), "sd_start", 1)
    )

    if bounded
        for (i, storage) in _PM.ref(pm, nw, :ne_storage)
            JuMP.set_lower_bound(sd[i], 0)
            JuMP.set_upper_bound(sd[i], storage["discharge_rating"])
        end
    end

    report && _PM.sol_component_value(pm, nw, :ne_storage, :sd_ne, _PM.ids(pm, nw, :ne_storage), sd)
end

function variable_storage_indicator(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            z = _PM.var(pm, nw)[:z_strg_ne] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_z_strg_ne",
                binary = true,
                start = 0
            )
        else
            z = _PM.var(pm, nw)[:z_strg_ne] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_z_strg_ne",
                lower_bound = 0,
                upper_bound = 1,
                start = 0
            )
        end
    else
        z = _PM.var(pm, nw)[:z_strg_ne] = _PM.var(pm, first_n)[:z_strg_ne]
    end
    report && _PM.sol_component_value(pm, nw, :ne_storage, :isbuilt, _PM.ids(pm, nw, :ne_storage), z)
end

function variable_storage_investment(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            investment = _PM.var(pm, nw)[:z_strg_ne_investment] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_z_strg_ne_investment",
                binary = true,
                start = 0
            )
        else
            investment = _PM.var(pm, nw)[:z_strg_ne_investment] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :ne_storage)], base_name="$(nw)_z_strg_ne_investment",
                lower_bound = 0,
                upper_bound = 1,
                start = 0
            )
        end
    else
        investment = _PM.var(pm, nw)[:z_strg_ne_investment] = _PM.var(pm, first_n)[:z_strg_ne_investment]
    end
    report && _PM.sol_component_value(pm, nw, :ne_storage, :investment, _PM.ids(pm, nw, :ne_storage), investment)
end

# ####################################################
# Constraint Templates: They are used to do all data manipuations and return a function with the same name,
# this way the constraint itself only containts the mathematical formulation
# ###################################################

function constraint_storage_thermal_limit_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :ne_storage, i)
    constraint_storage_thermal_limit_ne(pm, nw, i, storage["thermal_rating"])
end

function constraint_storage_losses_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :ne_storage, i)

    constraint_storage_losses_ne(pm, nw, i, storage["storage_bus"], storage["r"], storage["x"], storage["p_loss"], storage["q_loss"])
end

function constraint_storage_bounds_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    constraint_storage_bounds_ne(pm, nw, i)
end

function constraint_storage_state(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :storage, i)

    if haskey(_PM.ref(pm, nw), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network data should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end
    constraint_storage_state_initial(pm, nw, i, storage["energy"], storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
end

function constraint_storage_state_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :ne_storage, i)

    if haskey(_PM.ref(pm, nw), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network data should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end
    constraint_storage_state_initial_ne(pm, nw, i, storage["energy"], storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
end

function constraint_storage_state(pm::_PM.AbstractPowerModel, i::Int, nw_1::Int, nw_2::Int)
    storage = _PM.ref(pm, nw_2, :storage, i)

    if haskey(_PM.ref(pm, nw_2), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw_2, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network $(nw_2) should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end

    if haskey(_PM.ref(pm, nw_1, :storage), i)
        constraint_storage_state(pm, nw_1, nw_2, i, storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
    else
        # if the storage device has status=0 in nw_1, then the stored energy variable will not exist. Initialize storage from data model instead.
        Memento.warn(_LOGGER, "storage component $(i) was not found in network $(nw_1) while building constraint_storage_state between networks $(nw_1) and $(nw_2). Using the energy value from the storage component in network $(nw_2) instead")
        constraint_storage_state_initial(pm, nw_2, i, storage["energy"], storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
    end
end

function constraint_storage_state_ne(pm::_PM.AbstractPowerModel, i::Int, nw_1::Int, nw_2::Int)
    storage = _PM.ref(pm, nw_2, :ne_storage, i)

    if haskey(_PM.ref(pm, nw_2), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw_2, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network $(nw_2) should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end

    if haskey(_PM.ref(pm, nw_1, :ne_storage), i)
        constraint_storage_state_ne(pm, nw_1, nw_2, i, storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
    else
        # if the storage device has status=0 in nw_1, then the stored energy variable will not exist. Initialize storage from data model instead.
        Memento.warn(_LOGGER, "storage component $(i) was not found in network $(nw_1) while building constraint_storage_state between networks $(nw_1) and $(nw_2). Using the energy value from the storage component in network $(nw_2) instead")
        constraint_storage_state_initial_ne(pm, nw_2, i, storage["energy"], storage["charge_efficiency"], storage["discharge_efficiency"], storage["stationary_energy_inflow"], storage["stationary_energy_outflow"], storage["self_discharge_rate"], time_elapsed)
    end
end

function constraint_storage_state_final(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :storage, i)
    constraint_storage_state_final(pm, nw, i, storage["energy"])
end

function constraint_storage_state_final_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :ne_storage, i)
    constraint_storage_state_final_ne(pm, nw, i, storage["energy"])
end

function constraint_storage_excl_slack(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    constraint_storage_excl_slack(pm, nw, i)
end

function constraint_storage_excl_slack_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    constraint_storage_excl_slack_ne(pm, nw, i)
end

function constraint_maximum_absorption(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :storage, i)

    if haskey(_PM.ref(pm, nw), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network data should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end
    constraint_maximum_absorption_initial(pm, nw, i, time_elapsed)
end

function constraint_maximum_absorption_ne(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    storage = _PM.ref(pm, nw, :ne_storage, i)

    if haskey(_PM.ref(pm, nw), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network data should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end
    constraint_maximum_absorption_initial_ne(pm, nw, i, time_elapsed)
end

function  constraint_maximum_absorption(pm::_PM.AbstractPowerModel, i::Int, nw_1::Int, nw_2::Int)
    storage = _PM.ref(pm, nw_2, :storage, i)

    if haskey(_PM.ref(pm, nw_2), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw_2, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network $(nw_2) should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end

    if haskey(_PM.ref(pm, nw_1, :storage), i)
        constraint_maximum_absorption(pm, nw_1, nw_2, i, time_elapsed)
    else
        # if the storage device has status=0 in nw_1, then the stored energy variable will not exist. Initialize storage from data model instead.
        Memento.warn(_LOGGER, "storage component $(i) was not found in network $(nw_1) while building constraint_storage_state between networks $(nw_1) and $(nw_2). Using the energy value from the storage component in network $(nw_2) instead")
        constraint_maximum_absorption_initial(pm, nw_2, i, time_elapsed)
    end
end

function  constraint_maximum_absorption_ne(pm::_PM.AbstractPowerModel, i::Int, nw_1::Int, nw_2::Int)
    storage = _PM.ref(pm, nw_2, :ne_storage, i)

    if haskey(_PM.ref(pm, nw_2), :time_elapsed)
        time_elapsed = _PM.ref(pm, nw_2, :time_elapsed)
    else
        Memento.warn(_LOGGER, "network $(nw_2) should specify time_elapsed, using 1.0 as a default")
        time_elapsed = 1.0
    end

    if haskey(_PM.ref(pm, nw_1, :ne_storage), i)
        constraint_maximum_absorption_ne(pm, nw_1, nw_2, i, time_elapsed)
    else
        # if the storage device has status=0 in nw_1, then the stored energy variable will not exist. Initialize storage from data model instead.
        Memento.warn(_LOGGER, "storage component $(i) was not found in network $(nw_1) while building constraint_storage_state between networks $(nw_1) and $(nw_2). Using the energy value from the storage component in network $(nw_2) instead")
        constraint_maximum_absorption_initial_ne(pm, nw_2, i, time_elapsed)
    end
end

function constraint_ne_storage_activation(pm::_PM.AbstractPowerModel, i::Int, prev_nws::Vector{Int}, nw::Int)
    investment_horizon = [nw]
    lifetime = _PM.ref(pm, nw, :ne_storage, i, "lifetime")
    for n in Iterators.reverse(prev_nws[max(end-lifetime+2,1):end])
        i in _PM.ids(pm, n, :ne_storage) ? push!(investment_horizon, n) : break
    end
    constraint_ne_storage_activation(pm, nw, i, investment_horizon)
end

####################################################
############### Constraints
###################################################

function _PM.constraint_storage_thermal_limit(pm::BFARadPowerModel, n::Int, i, rating)
    ps = _PM.var(pm, n, :ps, i)
    qs = _PM.var(pm, n, :qs, i)

    c_perp = cos(π/8) # ~0.92
    c_diag = sin(π/8) + cos(π/8) # == cos(π/8) * sqrt(2), ~1.31

    JuMP.@constraint(pm.model, -c_perp*rating <= ps      <= c_perp*rating)
    JuMP.@constraint(pm.model, -c_perp*rating <=      qs <= c_perp*rating)
    JuMP.@constraint(pm.model, -c_diag*rating <= ps + qs <= c_diag*rating)
    JuMP.@constraint(pm.model, -c_diag*rating <= ps - qs <= c_diag*rating)
end

function constraint_storage_thermal_limit_ne(pm::_PM.AbstractActivePowerModel, n::Int, i, rating)
    ps = _PM.var(pm, n, :ps_ne, i)

    JuMP.lower_bound(ps) < -rating && JuMP.set_lower_bound(ps, -rating)
    JuMP.upper_bound(ps) >  rating && JuMP.set_upper_bound(ps,  rating)
end

function constraint_storage_thermal_limit_ne(pm::BFARadPowerModel, n::Int, i, rating)
    ps = _PM.var(pm, n, :ps_ne, i)
    qs = _PM.var(pm, n, :qs_ne, i)

    c_perp = cos(π/8) # ~0.92
    c_diag = sin(π/8) + cos(π/8) # == cos(π/8) * sqrt(2), ~1.31

    JuMP.@constraint(pm.model, -c_perp*rating <= ps      <= c_perp*rating)
    JuMP.@constraint(pm.model, -c_perp*rating <=      qs <= c_perp*rating)
    JuMP.@constraint(pm.model, -c_diag*rating <= ps + qs <= c_diag*rating)
    JuMP.@constraint(pm.model, -c_diag*rating <= ps - qs <= c_diag*rating)
end

function constraint_storage_losses_ne(pm::_PM.AbstractAPLossLessModels, n::Int, i, bus, r, x, p_loss, q_loss)
    ps = _PM.var(pm, n, :ps_ne, i)
    sc = _PM.var(pm, n, :sc_ne, i)
    sd = _PM.var(pm, n, :sd_ne, i)

    JuMP.@constraint(pm.model, ps + (sd - sc) == p_loss)
end

"Neglects the active and reactive loss terms associated with the squared current magnitude."
function constraint_storage_losses_ne(pm::_PM.AbstractBFAModel, n::Int, i, bus, r, x, p_loss, q_loss)
    ps  = _PM.var(pm, n, :ps_ne, i)
    qs  = _PM.var(pm, n, :qs_ne, i)
    sc  = _PM.var(pm, n, :sc_ne, i)
    sd  = _PM.var(pm, n, :sd_ne, i)
    qsc = _PM.var(pm, n, :qsc_ne, i)

    JuMP.@constraint(pm.model, ps + (sd - sc) == p_loss)
    JuMP.@constraint(pm.model, qs == qsc + q_loss)
end

function constraint_storage_state_initial(pm::_PM.AbstractPowerModel, n::Int, i::Int, energy, charge_eff, discharge_eff, inflow, outflow, self_discharge_rate, time_elapsed)
    sc = _PM.var(pm, n, :sc, i)
    sd = _PM.var(pm, n, :sd, i)
    se = _PM.var(pm, n, :se, i)

    JuMP.@constraint(pm.model, se == ((1-self_discharge_rate)^time_elapsed)*energy + time_elapsed*(charge_eff*sc - sd/discharge_eff + inflow - outflow))
end

function constraint_storage_state_initial_ne(pm::_PM.AbstractPowerModel, n::Int, i::Int, energy, charge_eff, discharge_eff, inflow, outflow, self_discharge_rate, time_elapsed)
    sc = _PM.var(pm, n, :sc_ne, i)
    sd = _PM.var(pm, n, :sd_ne, i)
    se = _PM.var(pm, n, :se_ne, i)
    z = _PM.var(pm, n, :z_strg_ne, i)

    JuMP.@constraint(pm.model, se == ((1-self_discharge_rate)^time_elapsed)*energy*z + time_elapsed*(charge_eff*sc - sd/discharge_eff + inflow * z - outflow * z))
end

function constraint_storage_state(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int, charge_eff, discharge_eff, inflow, outflow, self_discharge_rate, time_elapsed)
    sc_2 = _PM.var(pm, n_2, :sc, i)
    sd_2 = _PM.var(pm, n_2, :sd, i)
    se_2 = _PM.var(pm, n_2, :se, i)
    se_1 = _PM.var(pm, n_1, :se, i)

    JuMP.@constraint(pm.model, se_2 == ((1-self_discharge_rate)^time_elapsed)*se_1 + time_elapsed*(charge_eff*sc_2 - sd_2/discharge_eff + inflow - outflow))
end

function constraint_storage_state_ne(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int, charge_eff, discharge_eff, inflow, outflow, self_discharge_rate, time_elapsed)
    sc_2 = _PM.var(pm, n_2, :sc_ne, i)
    sd_2 = _PM.var(pm, n_2, :sd_ne, i)
    se_2 = _PM.var(pm, n_2, :se_ne, i)
    se_1 = _PM.var(pm, n_1, :se_ne, i)
    z = _PM.var(pm, n_2, :z_strg_ne, i)

    JuMP.@constraint(pm.model, se_2 == ((1-self_discharge_rate)^time_elapsed)*se_1 + time_elapsed*(charge_eff*sc_2 - sd_2/discharge_eff + inflow * z - outflow * z))
end

function constraint_storage_state_final(pm::_PM.AbstractPowerModel, n::Int, i::Int, energy)
    se = _PM.var(pm, n, :se, i)

    JuMP.@constraint(pm.model, se >= energy)
end

function constraint_storage_state_final_ne(pm::_PM.AbstractPowerModel, n::Int, i::Int, energy)
    se = _PM.var(pm, n, :se_ne, i)
    z = _PM.var(pm, n, :z_strg_ne, i)

    JuMP.@constraint(pm.model, se >= energy * z)
end

function constraint_maximum_absorption_initial(pm::_PM.AbstractPowerModel, n::Int, i::Int, time_elapsed)
    sc = _PM.var(pm, n, :sc, i)
    e_abs = _PM.var(pm, n, :e_abs, i)

    JuMP.@constraint(pm.model, e_abs == time_elapsed * sc)
end

function constraint_maximum_absorption_initial_ne(pm::_PM.AbstractPowerModel, n::Int, i::Int, time_elapsed)
    sc = _PM.var(pm, n, :sc_ne, i)
    e_abs = _PM.var(pm, n, :e_abs_ne, i)

    JuMP.@constraint(pm.model, e_abs == time_elapsed * sc)
end

function constraint_maximum_absorption(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int, time_elapsed)
    sc_2 = _PM.var(pm, n_2, :sc, i)
    e_abs_2 = _PM.var(pm, n_2, :e_abs, i)
    e_abs_1 = _PM.var(pm, n_1, :e_abs, i)

    JuMP.@constraint(pm.model, e_abs_2 - e_abs_1 == time_elapsed * sc_2)
end

function constraint_maximum_absorption_ne(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int, time_elapsed)
    sc_2 = _PM.var(pm, n_2, :sc_ne, i)
    e_abs_2 = _PM.var(pm, n_2, :e_abs_ne, i)
    e_abs_1 = _PM.var(pm, n_1, :e_abs_ne, i)

    JuMP.@constraint(pm.model, e_abs_2 - e_abs_1 == time_elapsed * sc_2)
end

function constraint_storage_excl_slack(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    sc = _PM.var(pm, n, :sc, i)
    sd = _PM.var(pm, n, :sd, i)
    sc_max = JuMP.upper_bound(sc)
    sd_max = JuMP.upper_bound(sd)
    s_bound = max(sc_max, sd_max)

    JuMP.@constraint(pm.model, sc + sd <= s_bound)
end

function constraint_storage_excl_slack_ne(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    sc = _PM.var(pm, n, :sc_ne, i)
    sd = _PM.var(pm, n, :sd_ne, i)
    sc_max = JuMP.upper_bound(sc)
    sd_max = JuMP.upper_bound(sd)
    s_bound = max(sc_max, sd_max)

    JuMP.@constraint(pm.model, sc + sd <= s_bound)
end

function constraint_storage_bounds_ne(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    se = _PM.var(pm, n, :se_ne, i)
    sc = _PM.var(pm, n, :sc_ne, i)
    sd = _PM.var(pm, n, :sd_ne, i)
    ps = _PM.var(pm, n, :ps_ne, i)
    qs = _PM.var(pm, n, :qs_ne, i)
    z = _PM.var(pm, n, :z_strg_ne, i)

    se_min = JuMP.lower_bound(se)
    se_max = JuMP.upper_bound(se)
    sc_min = JuMP.lower_bound(sc)
    sc_max = JuMP.upper_bound(sc)
    sd_min = JuMP.lower_bound(sd)
    sd_max = JuMP.upper_bound(sd)
    ps_min = JuMP.lower_bound(ps)
    ps_max = JuMP.upper_bound(ps)
    qs_min = JuMP.lower_bound(qs)
    qs_max = JuMP.upper_bound(qs)

    JuMP.@constraint(pm.model, se  <= se_max * z)
    JuMP.@constraint(pm.model, se  >= se_min * z)
    JuMP.@constraint(pm.model, sc  <= sc_max * z)
    JuMP.@constraint(pm.model, sc  >= sc_min * z)
    JuMP.@constraint(pm.model, sd  <= sd_max * z)
    JuMP.@constraint(pm.model, sd  >= sd_min * z)
    JuMP.@constraint(pm.model, ps  <= ps_max * z)
    JuMP.@constraint(pm.model, ps  >= ps_min * z)
    JuMP.@constraint(pm.model, qs  <= qs_max * z)
    JuMP.@constraint(pm.model, qs  >= qs_min * z)
end

function constraint_storage_bounds_ne(pm::_PM.AbstractActivePowerModel, n::Int, i::Int)
    se = _PM.var(pm, n, :se_ne, i)
    sc = _PM.var(pm, n, :sc_ne, i)
    sd = _PM.var(pm, n, :sd_ne, i)
    ps = _PM.var(pm, n, :ps_ne, i)
    z = _PM.var(pm, n, :z_strg_ne, i)

    se_min = JuMP.lower_bound(se)
    se_max = JuMP.upper_bound(se)
    sc_min = JuMP.lower_bound(sc)
    sc_max = JuMP.upper_bound(sc)
    sd_min = JuMP.lower_bound(sd)
    sd_max = JuMP.upper_bound(sd)
    ps_min = JuMP.lower_bound(ps)
    ps_max = JuMP.upper_bound(ps)

    JuMP.@constraint(pm.model, se  <= se_max * z)
    JuMP.@constraint(pm.model, se  >= se_min * z)
    JuMP.@constraint(pm.model, sc  <= sc_max * z)
    JuMP.@constraint(pm.model, sc  >= sc_min * z)
    JuMP.@constraint(pm.model, sd  <= sd_max * z)
    JuMP.@constraint(pm.model, sd  >= sd_min * z)
    JuMP.@constraint(pm.model, ps  <= ps_max * z)
    JuMP.@constraint(pm.model, ps  >= ps_min * z)
end

function constraint_ne_storage_activation(pm::_PM.AbstractPowerModel, n::Int, i::Int, horizon::Vector{Int})
    indicator = _PM.var(pm, n, :z_strg_ne, i)
    investments = _PM.var.(Ref(pm), horizon, :z_strg_ne_investment, i)

    JuMP.@constraint(pm.model, indicator == sum(investments))
end

# Force investment decisions in two different nws to be the same.
function constraint_ne_storage_investment_same(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int)
    investment_1 = _PM.var(pm, n_1, :z_strg_ne_investment, i)
    investment_2 = _PM.var(pm, n_2, :z_strg_ne_investment, i)

    JuMP.@constraint(pm.model, investment_1 == investment_2)
end


text:
# Adaptation of BayesLogit C code by Nicholas Polson, James Scott, Jesse Windle, 2012-2019

function jj_m1(b::Float64, z::Float64)
    z = abs(z)
    m1 = 0.0
    if (z > 1e-12)
        m1 = b * tanh(z) / z
    else
        m1 = b * (1 - (1.0 / 3) * z^2 + (2.0 / 15) * z^4 - (17.0 / 315) * z^6)
    end
    return m1::Float64
end

function jj_m2(b::Float64, z::Float64)
    z = abs(z)
    m2 = 0.0
    if (z > 1e-12)
        m2 = (b + 1) * b * (tanh(z) / z)^2 + b * ((tanh(z) - z) / (z^3))
    else
        m2 =
            (b + 1) * b * (1 - (1.0 / 3) * z^2 + (2.0 / 15) * z^4 - (17.0 / 315) * z^6)^2 +
            b * ((-1.0 / 3) + (2.0 / 15) * z^2 - (17.0 / 315) * z^4)
    end
    return m2
end

function pg_m1(b::Float64, z::Float64)
    return jj_m1(b, 0.5 * z) * 0.25
end

function pg_m2(b::Float64, z::Float64)
    return jj_m2(b, 0.5 * z) * 0.0625
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@info "Running functionality tests..."
@testset "Stochastic Programs: Functionality" begin
    tol = 1e-2
    for (model, _scenarios, res, name) in problems
        sp = instantiate(model, _scenarios, optimizer = GLPK.Optimizer)
        @testset "SP Constructs: $name" begin
            optimize!(sp, cache = true)
            @test termination_status(sp) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp), res.x̄, rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp, i), res.ȳ[i], rtol = tol)
            end
            @test isapprox(objective_value(sp), res.VRP, rtol = tol)
            @test isapprox(EWS(sp), res.EWS, rtol = tol)
            @test isapprox(EVPI(sp), res.EVPI, rtol = tol)
            @test isapprox(VSS(sp), res.VSS, rtol = tol)
            @test isapprox(EV(sp), res.EV, rtol = tol)
            @test isapprox(EEV(sp), res.EEV, rtol = tol)
        end
        @testset "Inequalities: $name" begin
            @test EWS(sp) <= VRP(sp)
            @test VRP(sp) <= EEV(sp)
            @test VSS(sp) >= 0
            @test EVPI(sp) >= 0
            @test VSS(sp) <= EEV(sp) - EV(sp)
            @test EVPI(sp) <= EEV(sp) - EV(sp)
        end
        @testset "Copying: $name" begin
            sp_copy = copy(sp, optimizer = GLPK.Optimizer)
            add_scenarios!(sp_copy, scenarios(sp))
            @test num_scenarios(sp_copy) == num_scenarios(sp)
            generate!(sp_copy)
            @test num_subproblems(sp_copy) == num_subproblems(sp)
            optimize!(sp)
            optimize!(sp_copy)
            @test termination_status(sp_copy) == MOI.OPTIMAL
            @test isapprox(optimal_decision(sp_copy), optimal_decision(sp), rtol = tol)
            for i in 1:num_scenarios(sp)
                @test isapprox(optimal_recourse_decision(sp_copy, i), optimal_recourse_decision(sp, i), rtol = tol)
            end
            @test isapprox(objective_value(sp_copy), objective_value(sp), rtol = tol)
            @test isapprox(EWS(sp_copy), EWS(sp), rtol = tol)
            @test isapprox(EVPI(sp_copy), EVPI(sp), rtol = tol)
            @test isapprox(VSS(sp_copy), VSS(sp), rtol = tol)
            @test isapprox(EV(sp_copy), EV(sp), rtol = tol)
            @test isapprox(EEV(sp_copy), EEV(sp), rtol = tol)
        end
    end
    @testset "Sampling" begin
        sampled_sp = instantiate(simple, sampler, 100, optimizer = GLPK.Optimizer)
        generate!(sampled_sp)
        @test num_scenarios(sampled_sp) == 100
        @test isapprox(stage_probability(sampled_sp), 1.0)
        StochasticPrograms.sample!(sampled_sp, sampler, 100)
        generate!(sampled_sp)
        @test num_scenarios(sampled_sp) == 200
        @test isapprox(stage_probability(sampled_sp), 1.0)
    end
    @testset "Instant" begin
        optimize!(simple_sp)
        @test termination_status(simple_sp) == MOI.OPTIMAL
        @test isapprox(optimal_decision(simple_sp), simple_res.x̄, rtol = tol)
        for i in 1:num_scenarios(simple_sp)
            @test isapprox(optimal_recourse_decision(simple_sp, i), simple_res.ȳ[i], rtol = tol)
        end
        @test isapprox(objective_value(simple_sp), simple_res.VRP, rtol = tol)
        @test isapprox(EWS(simple_sp), simple_res.EWS, rtol = tol)
        @test isapprox(EVPI(simple_sp), simple_res.EVPI, rtol = tol)
        @test isapprox(VSS(simple_sp), simple_res.VSS, rtol = tol)
        @test isapprox(EV(simple_sp), simple_res.EV, rtol = tol)
        @test isapprox(EEV(simple_sp), simple_res.EEV, rtol = tol)
    end
    @testset "SMPS" begin
        simple_smps = read("io/smps/simple.smps", StochasticProgram, optimizer = GLPK.Optimizer)
        optimize!(simple_smps)
        @test termination_status(simple_smps) == MOI.OPTIMAL
        @test isapprox(optimal_decision(simple_smps), simple_res.x̄, rtol = tol)
        for i in 1:num_scenarios(simple_smps)
            @test isapprox(optimal_recourse_decision(simple_smps, i), simple_res.ȳ[i], rtol = tol)
        end
        @test isapprox(objective_value(simple_smps), simple_res.VRP, rtol = tol)
        @test isapprox(EWS(simple_smps), simple_res.EWS, rtol = tol)
        @test isapprox(EVPI(simple_smps), simple_res.EVPI, rtol = tol)
        @test isapprox(VSS(simple_smps), simple_res.VSS, rtol = tol)
        @test isapprox(EV(simple_smps), simple_res.EV, rtol = tol)
        @test isapprox(EEV(simple_smps), simple_res.EEV, rtol = tol)
    end
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
"""2M4
card1 = bb
card2 = bw
card3 = ww

ways to get b with card1: 2
ways to get b with card2: 1
ways to get b with card3: 0
total ways to get b: 3
otherside also black means this is card1
probability of card1 = (ways to get b with card1) / (total ways to get b) = 2/3
"""

# %%
"""2M5
card1 = bb
card2 = bw
card3 = ww
card4 = bb

ways to get b with card1: 2
ways to get b with card2: 1
ways to get b with card3: 0
ways to get b with card4: 2
total ways to get b: 5
otherside also black means this is card1 or card4
probability of card1 = (ways to get b with card1) / (total ways to get b) = 2/5
probability of card4 = (ways to get b with card4) / (total ways to get b) = 2/5
probability of either card1 or card4 = (probability of card1) + (probability of card4)
                                     = 2/5 + 2/5 = 4/5
"""

# %%
"""2M6
card1 = bb
card2 = bw
card3 = ww

ways to get b with card1: 2  |  ways to get card1: 1
ways to get b with card2: 1  |  ways to get card2: 2
ways to get b with card3: 0  |  ways to get card3: 3
total ways to get b: 4
otherside also black means this is card1
probability of card1 = (ways to get b with card1) / (total ways to get b) = 2/4
"""

# %%
"""2M7
card1 = bb
card2 = bw
card3 = ww

scenarios:

first draw: side of first draw, second draw: side of first draw -> ok/no

card1: b1, card2: b  -> no, second card not white
card1: b1, card2: w  -> ok
card1: b1, card3: w1 -> ok
card1: b1, card3: w2 -> ok

card1: b2, card2: b  -> no, second card not white
card1: b2, card2: w  -> ok
card1: b2, card3: w1 -> ok
card1: b2, card3: w2 -> ok

card2: b, card1: b1  -> no, second card not white
card2: b, card1: b2  -> no, second card not white
card2: b, card3: w1  -> ok
card2: b, card3: w2  -> ok

card2: w, card1: b1  -> no, first card not black
card2: w, card1: b2  -> no, first card not black
card2: w, card3: w1  -> no, first card not black
card2: w, card3: w2  -> no, first card not black

card3: w1, card1: b1 -> no, first card not black
card3: w1, card1: b2 -> no, first card not black
card3: w1, card2: b  -> no, first card not black
card3: w1, card2: w  -> no, first card not black

card3: w2, card1: b1 -> no, first card not black
card3: w2, card1: b2 -> no, first card not black
card3: w2, card2: b  -> no, first card not black
card3: w2, card2: w  -> no, first card not black

number of possible/ok scenarios: 8
possible scenarios with card1 as first card: 6
probability for card1 = 6/8 = 0.75
"""


text:
"""
    train!(gen_fn::GenerativeFunction, data_generator::Function,
           update::ParamUpdate,
           num_epoch, epoch_size, num_minibatch, minibatch_size; verbose::Bool=false)

Train the given generative function to maximize the expected conditional log probability (density) that `gen_fn` generates the assignment `constraints` given inputs, where the expectation is taken under the output distribution of `data_generator`.

The function `data_generator` is a function of no arguments that returns a tuple `(inputs, constraints)` where `inputs` is a `Tuple` of inputs (arguments) to `gen_fn`, and `constraints` is an `ChoiceMap`.
`conf` configures the optimization algorithm used.
`param_lists` is a map from generative function to lists of its parameters.
This is equivalent to minimizing the expected KL divergence from the conditional distribution `constraints | inputs` of the data generator to the distribution represented by the generative function, where the expectation is taken under the marginal distribution on `inputs` determined by the data generator.
"""
function train!(gen_fn::GenerativeFunction, data_generator::Function,
                update::ParamUpdate;
                num_epoch=1, epoch_size=1, num_minibatch=1, minibatch_size=1,
                evaluation_size=epoch_size, verbose=false)

    history = Vector{Float64}(undef, num_epoch)
    for epoch=1:num_epoch

        # generate data for epoch
        if verbose
            println("epoch $epoch: generating $epoch_size training examples...")
        end
        epoch_inputs = Vector{Tuple}(undef, epoch_size)
        epoch_choice_maps = Vector{ChoiceMap}(undef, epoch_size)
        for i=1:epoch_size
            (epoch_inputs[i], epoch_choice_maps[i]) = data_generator()
        end

        # train on epoch data
        if verbose
            println("epoch $epoch: training using $num_minibatch minibatches of size $minibatch_size...")
        end
        for minibatch=1:num_minibatch
            permuted = Random.randperm(epoch_size)
            minibatch_idx = permuted[1:minibatch_size]
            minibatch_inputs = epoch_inputs[minibatch_idx]
            minibatch_choice_maps = epoch_choice_maps[minibatch_idx]
            for (inputs, constraints) in zip(minibatch_inputs, minibatch_choice_maps)
                (trace, _) = generate(gen_fn, inputs, constraints)
                retval_grad = accepts_output_grad(gen_fn) ? zero(get_retval(trace)) : nothing
                accumulate_param_gradients!(trace, retval_grad)
            end
            apply!(update)
        end

        # evaluate score on held out data
        if verbose
            println("epoch $epoch: evaluating on $evaluation_size examples...")
        end
        avg_score = 0.
        for i=1:evaluation_size
            (inputs, constraints) = data_generator()
            (_, weight) = generate(gen_fn, inputs, constraints)
            avg_score += weight
        end
        avg_score /= evaluation_size

        history[epoch] = avg_score

        if verbose
            println("epoch $epoch: est. objective value: $avg_score")
        end
    end
    return history
end

export train!

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


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

mutable struct IoK8sApiCoreV1AzureFilePersistentVolumeSource <: SwaggerModel
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretName::Any # spec type: Union{ Nothing, String } # spec name: secretName
    secretNamespace::Any # spec type: Union{ Nothing, String } # spec name: secretNamespace
    shareName::Any # spec type: Union{ Nothing, String } # spec name: shareName

    function IoK8sApiCoreV1AzureFilePersistentVolumeSource(;readOnly=nothing, secretName=nothing, secretNamespace=nothing, shareName=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1AzureFilePersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1AzureFilePersistentVolumeSource, Symbol("secretName"), secretName)
        setfield!(o, Symbol("secretName"), secretName)
        validate_property(IoK8sApiCoreV1AzureFilePersistentVolumeSource, Symbol("secretNamespace"), secretNamespace)
        setfield!(o, Symbol("secretNamespace"), secretNamespace)
        validate_property(IoK8sApiCoreV1AzureFilePersistentVolumeSource, Symbol("shareName"), shareName)
        setfield!(o, Symbol("shareName"), shareName)
        o
    end
end # type IoK8sApiCoreV1AzureFilePersistentVolumeSource

const _property_map_IoK8sApiCoreV1AzureFilePersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretName")=>Symbol("secretName"), Symbol("secretNamespace")=>Symbol("secretNamespace"), Symbol("shareName")=>Symbol("shareName"))
const _property_types_IoK8sApiCoreV1AzureFilePersistentVolumeSource = Dict{Symbol,String}(Symbol("readOnly")=>"Bool", Symbol("secretName")=>"String", Symbol("secretNamespace")=>"String", Symbol("shareName")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1AzureFilePersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1AzureFilePersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1AzureFilePersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1AzureFilePersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1AzureFilePersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1AzureFilePersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1AzureFilePersistentVolumeSource)
    (getproperty(o, Symbol("secretName")) === nothing) && (return false)
    (getproperty(o, Symbol("shareName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1AzureFilePersistentVolumeSource }, name::Symbol, val)
end


text:
mutable struct SparseSolver{dim} <: KKTSolver{SqrScaling{dim}}
	k0::Vector{Float64}
	k1::Vector{Float64}
	k2::Vector{Float64}

	n0::Vector{Float64}
	n1::Vector{Float64}

	m0::Vector{Float64}
	AA::SparseMatrixCSC{Float64, Int}
	Gt::SparseMatrixCSC{Float64, Int}
	GiW::SuiteSparse.CHOLMOD.Sparse{Float64}
	AtS::SuiteSparse.CHOLMOD.Sparse{Float64}
	Gint::SparseMatrixCSC{Float64, Int}
	Gfact::SuiteSparse.CHOLMOD.Factor{Float64}
	Afact::SuiteSparse.CHOLMOD.Factor{Float64}

	ws1::LdivWorkspace{Float64}
	dv1::SuiteSparse.CHOLMOD.Dense{Float64}
	dv2::SuiteSparse.CHOLMOD.Dense{Float64}

	ws2::LdivWorkspace{Float64}
	dv3::SuiteSparse.CHOLMOD.Dense{Float64}
	dv4::SuiteSparse.CHOLMOD.Dense{Float64}
	function SparseSolver(pr::Problem{C,n,m,k,sing}) where {C, n, m, k, sing}
		k0, k1, k2 = zeros(k), zeros(k), zeros(k)
		n0, n1 = zeros(n), zeros(n)
		m0 = zeros(m)

		ws1 = LdivWorkspace(Float64)
		dv1, dv2 = SuiteSparse.CHOLMOD.Dense(zeros(pr.n)), SuiteSparse.CHOLMOD.Dense(zeros(pr.n))
		ws2 = LdivWorkspace(Float64)
		dv3, dv4 = SuiteSparse.CHOLMOD.Dense(zeros(pr.m)), SuiteSparse.CHOLMOD.Dense(zeros(pr.m))

		AA = pr.A'*pr.A
		if !sing
			Gt = sparse(pr.G')
			desired_stype = 0
		else
			Gt = sparse(pr.G'*pr.G + pr.A'*pr.A)
			desired_stype = -1
		end
		cm = SuiteSparse.CHOLMOD.defaults(SuiteSparse.CHOLMOD.common_struct[Threads.threadid()])
		GtS = SuiteSparse.CHOLMOD.Sparse(Gt, desired_stype)
		AtCSC = sparse(pr.A')
		AtS = SuiteSparse.CHOLMOD.Sparse(AtCSC)
		#note that we're setting up both Gfact and Afact to be unsymmetric, so CHOLMOD will factorize AA'
		Gfact = SuiteSparse.CHOLMOD.analyze(GtS, cm)
		cholesky!(Gfact, GtS) # get the sparsity pattern right
		examplemat = SuiteSparse.CHOLMOD.Sparse(sparse((Gfact.UP\AtCSC)'), 0)
		Afact = SuiteSparse.CHOLMOD.analyze(examplemat, cm) # factorize W^t W, where L W = A or W = L^-1 A; thus W^t = A^t L^-t
		Gint = spzeros(pr.k, pr.n)
		return new{n}(k0, k1, k2, n0, n1, m0,
			AA, Gt, GtS, AtS, Gint, Gfact, Afact,
			ws1, dv1, dv2,
			ws2, dv3, dv4)
	end
end

function setup_iter(ss::SparseSolver, pr::Problem{C,n,m,k,sing}, s::State, scaling::SqrScaling) where {C, n, m, k, sing}
	if !sing
		copyto!(ss.GiW, ss.Gt)
		lmul!(scaling.iW, ss.GiW)
		cholesky!(ss.Gfact, ss.GiW)
		# do the low rank updates
		modify_factors!(pr.cones, ss.Gfact, scaling, pr.G)
	else
		Gint = ss.Gint
		copyto!(Gint, pr.G)
		lmul!(scaling.iWiW, Gint)
		i1 = pr.G' * Gint :: SparseMatrixCSC{Float64, Int}
		i1 .+= ss.AA
		cholesky!(ss.Gfact, i1)
		# do the low rank updates
		modify_factors!(pr.cones, ss.Gfact, scaling, pr.G)
	end
	# convert the factorization to an LLt one (is LDLt if we had to modify the factor)
	SuiteSparse.CHOLMOD.change_factor!(ss.Gfact, true, false, false, false)
	# solve L C^t = A^t giving C^t = L^-1 A^t, then transpose to get C = A L^-T
	i2 = SuiteSparse.CHOLMOD.spsolve(SuiteSparse.CHOLMOD.CHOLMOD_P, ss.Gfact, ss.AtS)
	Ct = SuiteSparse.CHOLMOD.spsolve(SuiteSparse.CHOLMOD.CHOLMOD_L, ss.Gfact, i2)
	Ctt = SuiteSparse.CHOLMOD.transpose_(Ct, 1)
	cholesky!(ss.Afact, Ctt)
end

function solve_kkt(ss::SparseSolver, pr::Problem{C,n,m,k,sing}, s::State, scaling::SqrScaling,
				   dx::Vector{Float64}, dy::Vector{Float64}, dz::Vector{Float64}, ds::Vector{Float64},
				   cx::Vector{Float64}, cy::Vector{Float64}, cz::Vector{Float64}, cs::Vector{Float64}) where {C, n, m, k, sing}
	iprod!(pr.cones, ss.k0, scaling.l, ds)
	scale!(pr.cones, scaling, ss.k0, ss.k1)
	bx = dx
	by = dy
	ss.k2 .= dz .- ss.k1
	iscale!(pr.cones, scaling, ss.k2, ss.k1)
	iscale!(pr.cones, scaling, ss.k1, ss.k1)
	mul!(ss.n0, pr.G', ss.k1) #ss.n0 = G'*W^-1 * W^-T*(dz - ss.k1)
	ss.n0 .+= bx
	At = pr.A' :: Adjoint{Float64, SparseMatrixCSC{Float64, Int}}
	if sing
		ss.n0 .+= At*by
	end

	# ss.m0 = A L-TL-1 ss.n0
	copyto!(ss.dv1, ss.n0)
	div!(ss.Gfact, ss.dv2, ss.dv1, ss.ws1)
	copyto!(ss.n1, ss.dv2)
	mul!(ss.m0, pr.A, ss.n1)
	ss.m0 .-= by
	copyto!(ss.dv3, ss.m0)
	div!(ss.Afact, ss.dv4, ss.dv3, ss.ws2)
	copyto!(cy, ss.dv4)
	if sing
		ss.m0 .= by .- cy
	else
		ss.m0 .= .-cy
	end
	mul!(ss.n1, At, ss.m0)
	ss.n0 .+= ss.n1
	# cx = L-TL-1 ss.n0
	copyto!(ss.dv1, ss.n0)
	div!(ss.Gfact, ss.dv2, ss.dv1, ss.ws1)
	copyto!(cx, ss.dv2)
	mul!(ss.k1, pr.G, cx)
	ss.k1 .-= ss.k2
	iscale!(pr.cones, scaling, ss.k1, cz)
	iscale!(pr.cones, scaling, cz, cz)
	scale!(pr.cones, scaling, cz, ss.k1)
	ss.k0 .-= ss.k1
	scale!(pr.cones, scaling, ss.k0, cs)
end


text:
@testset "Selections" begin

    @testset "Rank" begin
        Random.seed!(2);
        fitness = [5, 1, 3, 4, 2.0]
        s = ranklinear(1.0)
        @test s(fitness, 2) == [2,3]

        s = ranklinear(2.0)
        @test s([1,2], 2) == [2,2]
    end

    @testset "Uniform" begin
        s = uniformranking(2)
        @test sort(unique(s([1.0,2.0,3.0], 10))) == [1,2]
        @test sort(unique(s([5,2,3], 5))) == [3]
        @test_throws AssertionError s([1.,2.], 2)
    end

    @testset "Roulette" begin
        @test roulette([0.0,0.0,3.0], 2) == [3, 3]
        @test roulette([0,2,0], 2) == [2,2]
        @test rouletteinv([0.0,0.0,3.0], 2) == [1, 1]
        @test rouletteinv([1,2,0], 2) == [3,3]
    end

    @testset "Truncation" begin
        @test truncation([0.0,2.0,3.0], 2) == [1,2]
        @test truncation([0,2,1], 1) == [1]
        @test_throws AssertionError truncation([1.,2.], 5)
    end

    @testset "Tournament" begin
        Random.seed!(2);
        @test_throws AssertionError tournament(0)
        t = tournament(2)
        @test t([0,2,0],2) == [2,2]
        @test t([0.0,0.0,1.0],2) == [3,3]
    end

    @testset "SUS" begin
        @test sort(unique(sus([1.0,0.0,2.0], 5))) == [1,3]
        @test sort(unique(sus([0,1,2], 5))) == [2,3]
    end

    @testset "DE" begin
        N = 5
        rng = collect(1:10)
        idxs = random(rng, N)
        @test length(idxs) == N
        @test all(i-> i ∈ rng, idxs)
        @test !all(length(unique(random(rng, N))) == N for i in 1:100)

        idxs = randomoffset(rng, N)
        @test length(idxs) == N
        @test all(i-> i ∈ rng, idxs)
        @test all( let idxs = randomoffset(rng, N); let l = (idxs[1] - 1 + N)%10; idxs[N] == (l == 0 ? 10 : l) end end for i in 1:100 )

        idxs = best(reverse(rng), N)
        @test length(idxs) == N
        @test all(i-> i ∈ rng, idxs)
        @test idxs == fill(10, N)

        idxs = permutation(rng, N)
        @test length(idxs) == N
        @test all(i-> i ∈ rng, idxs)
        @test all(length(unique(permutation(rng, N))) == N for i in 1:100)

        @testset "Mutually Exclusive Indices" for i in idxs
            targets = Evolutionary.randexcl(rng, [i], 4)
            @test length(unique(push!(targets, i))) == N
        end
    end

end


text:
@testset "WorldClim Weather" begin

    using RasterDataSources: rastername, zipurl, zipname, zippath

    raster_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Weather", "prec", "wc2.1_2.5m_prec_2001-01.tif")
    @test rasterpath(WorldClim{Weather}, :prec; date=Date(2001, 1)) == raster_file
    @test rastername(WorldClim{Weather}, :prec; date=Date(2001, 1)) == "wc2.1_2.5m_prec_2001-01.tif"

    zip_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "Weather", "zips", "wc2.1_2.5m_prec_2010-2018.zip")
    @test zippath(WorldClim{Weather}, :prec; decade=Date(2010)) == zip_file
    @test zipurl(WorldClim{Weather}, :prec; decade=Date(2010)) ==
        URI(scheme="https", host="biogeo.ucdavis.edu", path="/data/worldclim/v2.1/hist/wc2.1_2.5m_prec_2010-2018.zip")
    @test zipname(WorldClim{Weather}, :prec; decade=Date(2010)) ==
        "wc2.1_2.5m_prec_2010-2018.zip"

    # These files are 3GB each. Probably too big to test in CI.
    # Not sure what to do about that.
    # getraster(WorldClim{Weather}, :prec; dates=Date(2001):Month(1):Date(2001, 12))
end


text:
@everywhere function test_function1(n, p=nothing)
    ProgressMeter.@info "Starting the stuff..."
    ProgressMeter.@parallelprogress p for i in 1:n
        sleep(0.1)
    end
    ProgressMeter.@info "Done with the stuff..."
end

println("Testing pmap progress *without* pmap")
test_function1(15)

println("Testing pmap progress")
ProgressMeter.pmap(test_function1, [25, 5, 10, 15, 15, 20])

@everywhere function test_function2(n, p=nothing)
    ProgressMeter.@info "Starting the stuff..."
    ProgressMeter.@parallelprogress p "LABEL HERE!!!  " for i in 1:n
        sleep(0.1)
    end
    ProgressMeter.@info "Done with the stuff..."
end

println("Testing pmap progress *without* pmap and with label")
test_function2(15)

println("Testing pmap progress with label")
ProgressMeter.pmap(test_function2, [25, 5, 10, 15, 15, 20])


text:
struct TruncatedPlotContext{T₁ <: AbstractControl, T₂ <: PlotContext,   T₃ <: NestedIntervalContext} <: AbstractPlotContext
    control::T₁
    plot_context::T₂
    interval_context::T₃
end

function (context::TruncatedPlotContext{<: PlotlinesControl, <:  PlotContext,   <: NestedIntervalContext})()
        control = context.control
        plot_context = context.plot_context
        interval_context = context.interval_context
        #CImGui.Begin("Main",C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus)
        isenabled(control) ? control(plot_context, interval_context) : nothing
        #CImGui.End()
end


text:
using StatsBase

# SVRG solver for large scale DMD problem
mutable struct DMD_SVRG_Options{T<:AbstractFloat}
    tau::Integer	# sample size
    eta::T		# step size
    itm::Integer
    tol::T
    ptf::Integer
    # pre-allocate variables
    fars::Vector{T}
    Gars::Matrix{T}
    gars::Array{Vector{T},1}
    prox::Function
end

function prox_null(x)
    return
end

# constructor of the SVRG options
function DMD_SVRG_Options(tau, eta, params; itm=1000, tol=1e-5, ptf=100,
                          prox=prox_null::Function)
    k = params.k;
    T = typeof(real(params.X[1]));

    tol  = T(tol);

    fars = zeros(T, n);
    Gars = zeros(T, 2*k, n);
    gars = col_view(Gars);

    return DMD_SVRG_Options(tau, eta, itm, tol, ptf, fars,
                            Gars, gars, prox)
end

# SVRG solver
function solveDMD_withSVRG(params, opts)
	# load all the variables
    tau    = opts.tau;
    eta    = opts.eta;
    prox = opts.prox;
    n    = params.n;
    k    = params.k;
    T    = typeof(real(params.X[1]));

    ar   = params.ar;
    arold = copy(ar);
    prox(ar);
    dar  = zeros(T, 2*k);
    ind  = collect(1:n);

    fars = opts.fars;
    gars = opts.gars;

    tfar = T(0.0); tgar = zeros(T, 2*k);
    dfar = T(0.0); dgar = zeros(T, 2*k);
    rfar = T(0.0); rgar = zeros(T, 2*k);

    itm  = opts.itm;
    tol  = opts.tol;
    ptf  = opts.ptf;

    obj_his = zeros(T, itm);
    err_his = zeros(T, itm);

    # first full step
    # initialize the gradient (might take a while)
    for id = 1:n
	fars[id] = abFunc(params, id);
	abGrad(gars[id], params, id);
    end
    tfar = sum(fars);
    sum!(tgar, opts.Gars);
    # update alpha
    copyto!(dar, tgar); dar .*= eta/n;
    copyto!(arold,ar);
    ar .-= dar;

    prox(ar);
    BLAS.axpy!(-1.0,ar,arold)
    err  = sqrt(sum(abs2,arold));
    noi  = 0;

    ind2 = collect(1:tau);

    while err ≥ tol
	# random sample columns
	sample!(ind,ind2); fill!(dgar, T(0.0)); dfar = T(0.0);
	for i = 1:tau
	    id = ind2[i];
	    # calcualte the objecitve
	    rfar = abFunc(params, id);
	    dfar = dfar + rfar - fars[id];
	    fars[id] = rfar;
	    # calcualte the gradient
	    dgar .-= gars[id];
	    abGrad(gars[id], params, id)
	    dgar .+= gars[id];
	end
	# update alpha
	copyto!(dar, tgar); dar.*= eta/n;
	BLAS.axpy!(eta/tau, dgar, dar);
        copyto!(arold,ar);
	ar .-= dar;
        prox(ar);
	# update tfar and tgar
	tfar  += dfar;
	tgar .+= dgar;

	# update information
        BLAS.axpy!(-1.0,ar,arold)
	err = sqrt(sum(abs2, arold));
	noi = noi + 1;
	obj_his[noi] = tfar;
	err_his[noi] = err;

	# print information
	noi % ptf == 0 && @printf("iter %5d, obj %1.2e, err %1.2e\n",
			          noi, tfar, err);
	noi ≥ itm && break;
    end
    return obj_his[1:noi], err_his[1:noi]
end


text:
using AsyPlots, LaTeXStrings, SpecialFunctions

x = range(0, stop=1, length = 400)

p = plot(collect(x),x.^44 .* (1 .- x).^11 / beta(44,11),xticks=Ticks([1]),yticks=Ticks([1])) + Plot(Path([4/5 -0.04; 4/5 0.04])) + Plot(Label(L"\frac{\alpha}{\alpha+\beta}",(4/5,-0.2));width=128)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
# a free moving object will be gripped
function changeJointFromFreeMotionToFix!(obj)
    obj.joint = fixedJoint
#    obj.parent.hasChildJoint = true
end

# a fixed object will be released
function changeJointFromFixToFreeMotion!(obj)
    r = obj.r_rel
    q = Modia3D.from_R(obj.R_rel)
    v = Modia3D.resolve2(obj.parent.R_abs,
        (obj.v0 - obj.parent.v0) )
    w = Modia3D.resolve2(obj.parent.R_abs,
        (obj.w - obj.parent.w) )
    #=
    println("r ", r)
    println("q ", q)
    println("v ", v)
    println("w ", w)
    println("obj.v0 ", obj.v0)
    println("obj.parent.v0 ", obj.parent.v0)
    println("obj.w ", obj.w)
    println("obj.parent.w ", obj.parent.w)
    =#
    obj.joint = FreeMotion(obj;
                r_start = r,
                q_start = q,
                v_start = v,
                w_start = w)
#    obj.parent.hasChildJoint = true

# resolve2 R*v
# resolve1 R'*v

#    obj.joint = FreeMotion(obj;
#                    r_start = obj.r_abs,
#                    q_start = Modia3D.from_R(obj.R_abs) )
                   # v_start=v_start,
                   # w_start=w_start) #, v_start = v_start, q_start = obj.R_abs,

end