text:
"""
    create!(newIdx::Int, size::Int, rng=GLOBAL_RNG)

Propose to generate new index (uniformly) randomly in [1, size]

# Arguments
- `newIdx`:  index ∈ [1, size]
- `size` : up limit of the index
- `rng=GLOBAL_RNG` : random number generator
"""
@inline function create!(d::Discrete, idx::Int, config)
    (idx >= length(d.data) - 1) && error("$idx overflow!")
    d[idx] = rand(config.rng, d.lower:d.upper)
    return Float64(d.upper - d.lower + 1) # lower:upper has upper-lower+1 elements!
end

@inline createRollback!(d::Discrete, idx::Int, config) = nothing

"""
    remove!(newIdx::Int, size::Int, rng=GLOBAL_RNG)

Propose to remove the old index in [1, size]

# Arguments
- `oldIdx`:  index ∈ [1, size]
- `size` : up limit of the index
- `rng=GLOBAL_RNG` : random number generator
"""
@inline function remove!(d::Discrete, idx::Int, config)
    (idx >= length(d.data) - 1) && error("$idx overflow!")
    return 1.0 / Float64(d.upper - d.lower + 1)
end

@inline removeRollback!(d::Discrete, idx::Int, config) = nothing

"""
    shift!(d::Discrete, idx::Int, config)

Propose to shift the old index in [1, size] to a new index

# Arguments
- `oldIdx`:  old index ∈ [1, size]
- `newIdx`:  new index ∈ [1, size], will be modified!
- `size` : up limit of the index
- `rng=GLOBAL_RNG` : random number generator
"""
@inline function shift!(d::Discrete, idx::Int, config)
    (idx >= length(d.data) - 1) && error("$idx overflow!")
    d[end] = d[idx] # save the current variable
    d[idx] = rand(config.rng, d.lower:d.upper)
    return 1.0
end

@inline function shiftRollback!(d::Discrete, idx::Int, config)
    (idx >= length(d.data) - 1) && error("$idx overflow!")
    d[idx] = d[end]
end

"""
    swap!(d::Discrete, idx1::Int, idx2::Int, config)

 Swap the variables idx1 and idx2

"""
@inline function swap!(d::Discrete, idx1::Int, idx2::Int, config)
    ((idx1 >= length(d.data) - 1) || (idx2 >= length(d.data) - 1)) && error("$idx overflow!")
    d[idx1], d[idx2] = d[idx2], d[idx1]
    return 1.0
end

@inline function swapRollback!(d::Discrete, idx1::Int, idx2::Int, config)
    ((idx1 >= length(d.data) - 1) || (idx2 >= length(d.data) - 1)) && error("$idx overflow!")
    d[idx1], d[idx2] = d[idx2], d[idx1]
end

"""
    create!(K::FermiK{D}, idx::Int, rng=GLOBAL_RNG)

Propose to generate new Fermi K in [Kf-δK, Kf+δK)

# Arguments
- `newK`:  vector of dimension of d=2 or 3
"""
function create!(K::FermiK{D}, idx::Int, config) where {D}
    @assert idx > K.offset
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    rng = config.rng
    ############ Simple Way ########################
    # for i in 1:DIM
    #     newK[i] = Kf * (rand(rng) - 0.5) * 2.0
    # end
    # return (2.0 * Kf)^DIM
    ################################################

    Kamp = K.kF + (rand(rng) - 0.5) * 2.0 * K.δk
    (Kamp <= 0.0) && return 0.0
    # Kf-dK<Kamp<Kf+dK
    ϕ = 2π * rand(rng)
    if D == 3 # dimension 3
        θ = π * rand(rng)
        # newK .= Kamp .* Mom(cos(ϕ) * sin(θ), sin(ϕ) * sin(θ), cos(θ))
        # K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
        K.data[1, idx] = Kamp * cos(ϕ) * sin(θ)
        K.data[2, idx] = Kamp * sin(ϕ) * sin(θ)
        K.data[3, idx] = Kamp * cos(θ)
        return 2 * K.δk * 2π * π * (sin(θ) * Kamp^2)
        # prop density of KAmp in [Kf-dK, Kf+dK), prop density of Phi
        # prop density of Theta, Jacobian
    else  # DIM==2
        # K[idx] = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
        K.data[1, idx] = Kamp * cos(ϕ)
        K.data[2, idx] = Kamp * sin(ϕ)
        return 2 * K.δk * 2π * Kamp
        # prop density of KAmp in [Kf-dK, Kf+dK), prop density of Phi, Jacobian
    end
end
createRollback!(K::FermiK{D}, idx::Int, config) where {D} = nothing

"""
    removeFermiK!(oldK, Kf=1.0, δK=0.5, rng=GLOBAL_RNG)

Propose to remove an existing Fermi K in [Kf-δK, Kf+δK)

# Arguments
- `oldK`:  vector of dimension of d=2 or 3
"""
function remove!(K::FermiK{D}, idx::Int, config) where {D}
    @assert idx > K.offset
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    ############## Simple Way #########################
    # for i in 1:DIM
    #     if abs(oldK[i]) > Kf
    #         return 0.0
    #     end
    # end
    # return 1.0 / (2.0 * Kf)^DIM
    ####################################################

    oldK = K[idx]
    Kamp = sqrt(dot(oldK, oldK))
    if !(K.kF - K.δk < Kamp < K.kF + K.δk)
        return 0.0
    end
    # (Kamp < Kf - dK || Kamp > Kf + dK) && return 0.0
    if D == 3 # dimension 3
        sinθ = sqrt(oldK[1]^2 + oldK[2]^2) / Kamp
        sinθ < 1.0e-15 && return 0.0
        return 1.0 / (2 * K.δk * 2π * π * sinθ * Kamp^2)
    else  # DIM==2
        return 1.0 / (2 * K.δk * 2π * Kamp)
    end
end

removeRollback!(K::FermiK{D}, idx::Int, config) where {D} = nothing

"""
    shiftK!(oldK, newK, step, rng=GLOBAL_RNG)

Propose to shift oldK to newK. Work for generic momentum vector
"""
function shift!(K::FermiK{D}, idx::Int, config) where {D}
    @assert idx > K.offset
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    K[end] = K[idx]  # save current K

    rng = config.rng
    x = rand(rng)
    if x < 1.0 / 3
        λ = 1.5
        ratio = 1.0 / λ + rand(rng) * (λ - 1.0 / λ)
        K[idx] = K[idx] * ratio
        return (D == 2) ? 1.0 : ratio
    elseif x < 2.0 / 3
        ϕ = rand(rng) * 2π
        if (D == 3)
            # sample uniformly on sphere, check http://corysimon.github.io/articles/uniformdistn-on-sphere/
            θ = acos(1 - 2 * rand(rng))
            Kamp = sqrt(K[idx][1]^2 + K[idx][2]^2 + K[idx][3]^2)
            # K[idx] = @SVector [Kamp * cos(ϕ) * sin(θ), Kamp * sin(ϕ) * sin(θ), Kamp * cos(θ)]
            K.data[1, idx] = Kamp * cos(ϕ) * sin(θ)
            K.data[2, idx] = Kamp * sin(ϕ) * sin(θ)
            K.data[3, idx] = Kamp * cos(θ)
            return 1.0
        else # D=2
            Kamp = sqrt(K[idx][1]^2 + K[idx][2]^2)
            # K = @SVector [Kamp * cos(ϕ), Kamp * sin(ϕ)]
            K.data[1, idx] = Kamp * cos(ϕ)
            K.data[2, idx] = Kamp * sin(ϕ)
            return 1.0
        end
    else
        Kc, dk = K[idx], K.δk
        if (D == 3)
            # K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk, Kc[3] + (rand(rng) - 0.5) * dk]
            K.data[1, idx] = Kc[1] + (rand(rng) - 0.5) * dk
            K.data[2, idx] = Kc[2] + (rand(rng) - 0.5) * dk
            K.data[3, idx] = Kc[3] + (rand(rng) - 0.5) * dk
        else # D=2
            # K[idx] = @SVector [Kc[1] + (rand(rng) - 0.5) * dk, Kc[2] + (rand(rng) - 0.5) * dk]
            K.data[1, idx] = Kc[1] + (rand(rng) - 0.5) * dk
            K.data[2, idx] = Kc[2] + (rand(rng) - 0.5) * dk
        end
        # K[idx] += (rand(rng, D) .- 0.5) .* K.δk
        return 1.0
    end
end

function shiftRollback!(K::FermiK{D}, idx::Int, config) where {D}
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    K[idx] = K[end]
end

@inline function swap!(K::FermiK{D}, idx1::Int, idx2::Int, config) where {D}
    ((idx1 >= length(K.data) - 1) || (idx2 >= length(K.data) - 1)) && error("$idx1 or $idx2 overflow!")
    if D == 2
        K.data[1, idx1], K.data[1, idx2] = K.data[1, idx2], K.data[1, idx1]
        K.data[2, idx1], K.data[2, idx2] = K.data[2, idx2], K.data[2, idx1]
    elseif D == 3
        K.data[1, idx1], K.data[1, idx2] = K.data[1, idx2], K.data[1, idx1]
        K.data[2, idx1], K.data[2, idx2] = K.data[2, idx2], K.data[2, idx1]
        K.data[3, idx1], K.data[3, idx2] = K.data[3, idx2], K.data[3, idx1]
    else
        error("not implemented!")
    end
    return 1.0
end

@inline function swapRollback!(K::FermiK{D}, idx1::Int, idx2::Int, config) where {D}
    ((idx1 >= length(K.data) - 1) || (idx2 >= length(K.data) - 1)) && error("$idx1 or $idx2 overflow!")
    if D == 2
        K.data[1, idx1], K.data[1, idx2] = K.data[1, idx2], K.data[1, idx1]
        K.data[2, idx1], K.data[2, idx2] = K.data[2, idx2], K.data[2, idx1]
    elseif D == 3
        K.data[1, idx1], K.data[1, idx2] = K.data[1, idx2], K.data[1, idx1]
        K.data[2, idx1], K.data[2, idx2] = K.data[2, idx2], K.data[2, idx1]
        K.data[3, idx1], K.data[3, idx2] = K.data[3, idx2], K.data[3, idx1]
    else
        error("not implemented!")
    end
end

"""
    create!(T::Tau, idx::Int, rng=GLOBAL_RNG)

Propose to generate new tau (uniformly) randomly in [0, β), return proposal probability

# Arguments
- `T`:  Tau variable
- `idx`: T.data[idx] will be updated
"""
@inline function create!(T::Tau, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[idx] = rand(config.rng) * T.β
    return T.β
end
@inline createRollback!(T::Tau, idx::Int, config) = nothing

"""
    remove(T::Tau, idx::Int, rng=GLOBAL_RNG)

Propose to remove old tau in [0, β), return proposal probability

# Arguments
- `T`:  Tau variable
- `idx`: T.data[idx] will be updated
"""
@inline function remove!(T::Tau, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    return 1.0 / T.β
end
@inline removeRollback!(T::Tau, idx::Int, config) = nothing

"""
    shift!(T::Tau, idx::Int, rng=GLOBAL_RNG)

Propose to shift an existing tau to a new tau, both in [0, β), return proposal probability

# Arguments
- `T`:  Tau variable
- `idx`: T.data[idx] will be updated
"""
@inline function shift!(T::Tau, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[end] = T[idx]
    rng = config.rng
    x = rand(rng)
    if x < 1.0 / 3
        T[idx] = T[idx] + 2 * T.λ * (rand(rng) - 0.5)
    elseif x < 2.0 / 3
        T[idx] = T.β - T[idx]
    else
        T[idx] = rand(rng) * T.β
    end

    if T[idx] < 0.0
        T[idx] += T.β
    elseif T[idx] > T.β
        T[idx] -= T.β
    end

    return 1.0
end

@inline function shiftRollback!(T::Tau, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[idx] = T[end]
end

"""
    swap!(T::Tau, idx1::Int, idx2::Int, config)

 Swap the variables idx1 and idx2

"""
@inline function swap!(T::Tau, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
    return 1.0
end

@inline function swapRollback!(T::Tau, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
end

"""
    create!(T::TauPair, idx::Int, rng=GLOBAL_RNG)

Propose to generate a new pair of tau (uniformly) randomly in [0, β), return proposal probability

# Arguments
- `T`:  TauPair variable
- `idx`: T.data[idx] will be updated
"""
@inline function create!(T::TauPair, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    rng = config.rng
    T[idx][1] = rand(rng) * T.β
    T[idx][2] = rand(rng) * T.β
    return T.β * T.β
end

@inline createRollback!(T::TauPair, idx::Int, config) = nothing

"""
    remove(T::TauPair, idx::Int, rng=GLOBAL_RNG)

Propose to remove an existing pair of tau in [0, β), return proposal probability

# Arguments
- `T`:  Tau variable
- `idx`: T.data[idx] will be updated
"""
@inline function remove!(T::TauPair, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    return 1.0 / T.β / T.β
end
@inline removeRollback!(T::TauPair, idx::Int, config) = nothing

"""
    shift!(T::TauPair, idx::Int, rng=GLOBAL_RNG)

Propose to shift an existing tau pair to a new tau pair, both in [0, β), return proposal probability

# Arguments
- `T`:  Tau variable
- `idx`: T.t[idx] will be updated
"""
@inline function shift!(T::TauPair, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[end] .= T[idx]
    rng = config.rng
    x = rand(rng)
    if x < 1.0 / 3
        T[idx][1] += 2 * T.λ * (rand(rng) - 0.5)
        T[idx][2] += 2 * T.λ * (rand(rng) - 0.5)
    elseif x < 2.0 / 3
        T[idx][1] = T.β - T[idx][1]
        T[idx][2] = T.β - T[idx][2]
    else
        T[idx][1] = rand(rng) * T.β
        T[idx][2] = rand(rng) * T.β
    end

    if T[idx][1] < 0.0
        T[idx][1] += T.β
    elseif T[idx][1] > T.β
        T[idx][1] -= T.β
    end

    if T[idx][2] < 0.0
        T[idx][2] += T.β
    elseif T[idx][2] > T.β
        T[idx][2] -= T.β
    end

    return 1.0
    # return 0.0
end

@inline function shiftRollback!(T::TauPair, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[idx] .= T[end]
end

@inline function swap!(T::TauPair, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
    return 1.0
end

@inline function swapRollback!(T::TauPair, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
end

"""
    create!(theta::Angle, idx::Int, rng=GLOBAL_RNG)

Propose to generate new angle (uniformly) randomly in [0, 2π), return proposal probability

# Arguments
- `theta`:  angle variable
- `idx`: theta.t[idx] will be updated
"""
@inline function create!(theta::Angle, idx::Int, config)
    (idx >= length(theta.data) - 1) && error("$idx overflow!")
    theta[idx] = rand(config.rng) * 2π
    return 2π
end
@inline createRollback!(theta::Angle, idx::Int, config) = nothing

"""
    remove(theta::Angle, idx::Int, rng=GLOBAL_RNG)

Propose to remove old theta in [0, 2π), return proposal probability

# Arguments
    - `theta`:  Tau variable
- `idx`: theta.t[idx] will be updated
"""
@inline function remove!(theta::Angle, idx::Int, config)
    (idx >= length(theta.data) - 1) && error("$idx overflow!")
    return 1.0 / 2.0 / π
end
@inline removeRollback!(theta::Angle, idx::Int, config) = nothing

"""
    shift!(theta::Angle, idx::Int, rng=GLOBAL_RNG)

Propose to shift the old theta to new theta, both in [0, 2π), return proposal probability

# Arguments
- `theta`:  angle variable
- `idx`: theta.t[idx] will be updated
"""
@inline function shift!(theta::Angle, idx::Int, config)
    (idx >= length(theta.data) - 1) && error("$idx overflow!")
    theta[end] = theta[idx]
    rng = config.rng
    x = rand(rng)
    if x < 1.0 / 3
        theta[idx] = theta[idx] + 2 * theta.λ * (rand(rng) - 0.5)
    elseif x < 2.0 / 3
        theta[idx] = 2π - theta[idx]
    else
        theta[idx] = rand(rng) * 2π
    end

    if theta[idx] < 0.0
        theta[idx] += 2π
    elseif theta[idx] > 2π
        theta[idx] -= 2π
    end

    return 1.0
end

@inline function shiftRollback!(theta::Angle, idx::Int, config)
    (idx >= length(theta.data) - 1) && error("$idx overflow!")
    theta[idx] = theta[end]
end

@inline function RFK_Weight(K::RadialFermiK, k)
    # norm = atan(K.kF/K.δk)/K.δk+π/2.0/K.δk
    # return 1.0/((k-K.kF)^2+K.δk^2)/norm
    if 0.0 <= k < K.kF - K.δk
        return 1.0 / 3.0 / (K.kF - K.δk)
    elseif k < K.kF + K.δk
        return 2.0 * K.δk / (3π) / ((k - K.kF)^2 + K.δk^2)
    else
        return 1.0 / 3 / (K.kF + K.δk) * exp(1 - k / (K.kF + K.δk))
    end
    return 0.0
end

@inline function RFK_Propose(K::RadialFermiK, x, y)
    # norm = atan(K.kF/K.δk)/K.δk+π/2.0/K.δk
    # return K.kF-K.δk*tan(atan(K.kF/K.δk)-K.δk*x*norm )
    if x < 1.0 / 3
        return y * (K.kF - K.δk)
    elseif x < 2.0 / 3
        return K.kF + K.δk * tan(π * (2.0 * y - 1) / 4.0)
    else
        return (K.kF + K.δk) * (1 - log1p(-y))
    end
end

"""
    create!(K::RadialFermiK, idx::Int, rng=GLOBAL_RNG)

Propose to generate new k randomly in [0, +inf), return proposal probability
k is generated uniformly on [0, K.kF-K.δk), Lorentzianly on [K.kF-K.δk,K.kF+K.δk),
and exponentially on [K.kF-K.δk, +inf).

# Arguments
- `K`:  k variable
- `idx`: K.data[idx] will be updated
"""
@inline function create!(K::RadialFermiK, idx::Int, config)
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    rng = config.rng
    x = rand(rng)
    y = rand(rng)
    K[idx] = RFK_Propose(K, x, y)
    if isnan(K[idx])
        return 0.0
    end
    return 1.0 / RFK_Weight(K, K[idx])
end
@inline createRollback!(K::RadialFermiK, idx::Int, config) = nothing

"""
    remove(K::RadialFermiK, idx::Int, rng=GLOBAL_RNG)

Propose to remove old k in [0, +inf), return proposal probability

# Arguments
- `K`:  K variable
- `idx`: K.data[idx] will be updated
"""
@inline function remove!(K::RadialFermiK, idx::Int, config)
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    return RFK_Weight(K, K[idx])
end
@inline removeRollback!(K::RadialFermiK, idx::Int, config) = nothing

"""
    shift!(K::RadialFermiK, idx::Int, rng=GLOBAL_RNG)

Propose to shift an existing k to a new k, both in [0, +inf), return proposal probability

# Arguments
- `K`:  K variable
- `idx`: K.data[idx] will be updated
"""
@inline function shift!(K::RadialFermiK, idx::Int, config)
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    K[end] = K[idx]
    rng = config.rng
    # K[idx] = RFK_Propose(K, rand(rng), rand(rng))
    # return RFK_Weight(K, K[end])/RFK_Weight(K, K[idx])
    x = rand(rng)
    if x < 1.0 / 3
        K[idx] = K[idx] + 2 * K.δk * (rand(rng) - 0.5)
        if K[idx] < 0.0
            K[idx] = -K[idx]
        end
        return 1.0
    elseif x < 2.0 / 3
        if K[idx] < K.kF - K.δk
            y = (K.kF - K.δk - K[idx]) / (K.kF - K.δk)
            K[idx] = RFK_Propose(K, 1.0, y)
        elseif K[idx] < K.kF + K.δk
            K[idx] = 2 * K.kF - K[idx]
        else
            y = exp(1 - K[idx] / (K.kF + K.δk))
            K[idx] = RFK_Propose(K, 0.0, y)
        end
        if isnan(K[idx])
            return 0.0
        end
        return RFK_Weight(K, K[end]) / RFK_Weight(K, K[idx])
    else
        K[idx] = RFK_Propose(K, rand(rng), rand(rng))
        if isnan(K[idx])
            return 0.0
        end
        return RFK_Weight(K, K[end]) / RFK_Weight(K, K[idx])
    end

end

@inline function shiftRollback!(K::RadialFermiK, idx::Int, config)
    (idx >= length(K.data) - 1) && error("$idx overflow!")
    K[idx] = K[end]
end

@inline function swap!(K::RadialFermiK, idx1::Int, idx2::Int, config)
    ((idx1 >= length(K.data) - 1) || (idx2 >= length(K.data) - 1)) && error("$idx1 or $idx2 overflow!")
    K[idx1], K[idx2] = K[idx2], K[idx1]
    return 1.0
end

@inline function swapRollback!(K::RadialFermiK, idx1::Int, idx2::Int, config)
    ((idx1 >= length(K.data) - 1) || (idx2 >= length(K.data) - 1)) && error("$idx1 or $idx2 overflow!")
    K[idx1], K[idx2] = K[idx2], K[idx1]
end

"""
    create!(T::Continuous, idx::Int, rng=GLOBAL_RNG)

Propose to generate new (uniform) variable randomly in [T.lower, T.lower+T.range), return proposal probability

# Arguments
- `T`:  Continuous variable
- `idx`: T.data[idx] will be updated
"""
@inline function create!(T::Continuous, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[idx] = rand(config.rng) * T.range + T.lower
    return T.range
end
@inline createRollback!(T::Continuous, idx::Int, config) = nothing

"""
    remove(T::Continuous, idx::Int, rng=GLOBAL_RNG)

Propose to remove old variable in [T.lower, T.lower+T.range), return proposal probability

# Arguments
- `T`:  Continuous variable
- `idx`: T.data[idx] will be updated
"""
@inline function remove!(T::Continuous, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    return 1.0 / T.range
end
@inline removeRollback!(T::Continuous, idx::Int, config) = nothing

"""
    shift!(T::Continuous, idx::Int, rng=GLOBAL_RNG)

Propose to shift an existing variable to a new one, both in [T.lower, T.lower+T.range), return proposal probability

# Arguments
- `T`:  Continuous variable
- `idx`: T.data[idx] will be updated
"""
@inline function shift!(T::Continuous, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[end] = T[idx]
    rng = config.rng
    x = rand(rng)
    if x < 1.0 / 2
        T[idx] = T[idx] + 2 * T.λ * (rand(rng) - 0.5)
    else
        T[idx] = rand(rng) * T.range + T.lower
    end

    if T[idx] < T.lower
        T[idx] += T.range
    elseif T[idx] > T.lower + T.range
        T[idx] -= T.range
    end

    return 1.0
end

@inline function shiftRollback!(T::Continuous, idx::Int, config)
    (idx >= length(T.data) - 1) && error("$idx overflow!")
    T[idx] = T[end]
end

@inline function swap!(T::Continuous, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
    return 1.0
end

@inline function swapRollback!(T::Continuous, idx1::Int, idx2::Int, config)
    ((idx1 >= length(T.data) - 1) || (idx2 >= length(T.data) - 1)) && error("$idx1 or $idx2 overflow!")
    T[idx1], T[idx2] = T[idx2], T[idx1]
end

# @inline function sample(model::Uniform{Int,D}, data, idx) where {D}
#     p = 1.0
#     for di in 1:D
#         data[idx, di] = rand(config.rng, model.lower[di]:model.upper[di])
#         p *= Float64(d.upper - d.lower + 1) # lower:upper has upper-lower+1 elements!
#     end
#     return p
# end

# @inline function sample(model::Uniform{T,D}, data, idx) where {T<:Real,D}
#     p = T(1)
#     for di in 1:D
#         data[idx, di] = rand(config.rng, T) * (model.upper[di] - model.lower[di]) + model.lower[di]
#         p *= d.upper[di] - d.lower[di]
#     end
#     return p
# end

# """
#     create!(newIdx::Int, size::Int, rng=GLOBAL_RNG)

# Propose to generate new index (uniformly) randomly in [1, size]

# # Arguments
# - `newIdx`:  index ∈ [1, size]
# - `size` : up limit of the index
# - `rng=GLOBAL_RNG` : random number generator
# """
# @inline function create!(d::Var{T,D,M}, idx::Int, config) where {T,D,M}
#     (idx >= length(d.data) - 1) && error("$idx overflow!")
#     d[idx] = rand(config.rng, d.lower:d.upper)
#     return Float64(d.upper - d.lower + 1) # lower:upper has upper-lower+1 elements!
# end

# @inline createRollback!(d::Discrete, idx::Int, config) = nothing

# """
#     remove!(newIdx::Int, size::Int, rng=GLOBAL_RNG)

# Propose to remove the old index in [1, size]

# # Arguments
# - `oldIdx`:  index ∈ [1, size]
# - `size` : up limit of the index
# - `rng=GLOBAL_RNG` : random number generator
# """
# @inline function remove!(d::Discrete, idx::Int, config)
#     (idx >= length(d.data) - 1) && error("$idx overflow!")
#     return 1.0 / Float64(d.upper - d.lower + 1)
# end

# @inline removeRollback!(d::Discrete, idx::Int, config) = nothing

# """
#     shift!(d::Discrete, idx::Int, config)

# Propose to shift the old index in [1, size] to a new index

# # Arguments
# - `oldIdx`:  old index ∈ [1, size]
# - `newIdx`:  new index ∈ [1, size], will be modified!
# - `size` : up limit of the index
# - `rng=GLOBAL_RNG` : random number generator
# """
# @inline function shift!(d::Discrete, idx::Int, config)
#     (idx >= length(d.data) - 1) && error("$idx overflow!")
#     d[end] = d[idx] # save the current variable
#     d[idx] = rand(config.rng, d.lower:d.upper)
#     return 1.0
# end

# @inline function shiftRollback!(d::Discrete, idx::Int, config)
#     (idx >= length(d.data) - 1) && error("$idx overflow!")
#     d[idx] = d[end]
# end

# """
#     swap!(d::Discrete, idx1::Int, idx2::Int, config)

#  Swap the variables idx1 and idx2

# """
# @inline function swap!(d::Discrete, idx1::Int, idx2::Int, config)
#     ((idx1 >= length(d.data) - 1) || (idx2 >= length(d.data) - 1)) && error("$idx overflow!")
#     d[idx1], d[idx2] = d[idx2], d[idx1]
#     return 1.0
# end

# @inline function swapRollback!(d::Discrete, idx1::Int, idx2::Int, config)
#     ((idx1 >= length(d.data) - 1) || (idx2 >= length(d.data) - 1)) && error("$idx overflow!")
#     d[idx1], d[idx2] = d[idx2], d[idx1]
# end


text:
function goToDirect(r, side)
    sum = 0
    n = 0
    while !isborder(r, side)
        move!(r, side)
        if ismarker(r)
            sum+=temperature(r)
            n+=1
        end
    end
    return sum, n
end

function draw(r)
    sum = 0
    n = 0

    if ismarker(r)
        n+=1
        sum+=temperature(r)
    end

    k = 0
    directions = (Ost, West)

    while !isborder(r, Nord)
        sum1, n1 = goToDirect(r, directions[k%2+1])
        sum+=sum1
        n+=n1
        move!(r, Nord)
        if ismarker(r)
            n+=1
            sum+=temperature(r)
        end
        k += 1
    end
    sum1, n1 = goToDirect(r, directions[k%2+1])
    sum+=sum1
    n+=n1
    println(sum/7)
end

ext: jl

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
module TestCrossValidator

using Test
using Random
using AutoMLPipeline
using AutoMLPipeline.CrossValidators
using AutoMLPipeline.DecisionTreeLearners
using AutoMLPipeline.Pipelines
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.Utils

function test_crossvalidator()
  racc = 50.0
  Random.seed!(123)
  acc(X,Y) = score(:accuracy,X,Y)
  data=getiris()
  X=data[:,1:4]
  Y=data[:,5] |> Vector{String}
  rf = RandomForest()
  @test crossvalidate(rf,X,Y,acc,10,false).mean > racc
  Random.seed!(123)
  ppl1 = Pipeline([RandomForest()])
  @test crossvalidate(ppl1,X,Y,acc,10,false).mean > racc
  Random.seed!(123)
  ohe = OneHotEncoder()
  stdsc= SKPreprocessor("StandardScaler")
  ppl2 = Pipeline([ohe,stdsc,RandomForest()])
  @test crossvalidate(ppl2,X,Y,acc,10,false).mean > racc
  Random.seed!(123)
  mpca = SKPreprocessor("PCA")
  mppca = SKPreprocessor("KernelPCA")
  mfa = SKPreprocessor("FactorAnalysis")
  mica = SKPreprocessor("FastICA")
  mrb = SKPreprocessor("RobustScaler")
  ppl3 = Pipeline(Dict(:machines=>[mrb,mica,mpca,mppca,RandomForest()]))
  @test crossvalidate(ppl3,X,Y,acc,10,false).mean > racc
  Random.seed!(123)
  fit!(ppl3,X,Y)
  @test size(transform!(ppl3,X))[1] == length(Y)
  Random.seed!(123)
  ppl5 = Pipeline(Dict(:machines=>[mrb,mica,mppca,RandomForest()]))
  @test crossvalidate(ppl5,X,Y,acc,10,false).mean > racc
end
@testset "CrossValidator" begin
  test_crossvalidator()
end

end


text:
module Agents

using Requires
using Distributed
using DataStructures
using LightGraphs
using DataFrames
using Random
using OpenStreetMapX

import Base.length
import LinearAlgebra

# Core structures of Agents.jl
include("core/agents.jl")
include("core/model.jl")
include("core/schedule.jl")
include("core/space_interaction_API.jl")

# Existing spaces
include("spaces/nothing.jl")
include("spaces/graph.jl")
include("spaces/grid.jl")
include("spaces/discrete.jl")
include("spaces/continuous.jl")
include("spaces/openstreetmap.jl")
include("spaces/utilities.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")
include("simulations/ensemblerun.jl")

# Other advanced features
include("pathfinding/grid_pathfinder.jl")

function __init__()
    # Plot recipes
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("visualization/plot-recipes.jl")
    end
    # Workaround for Documenter.jl, so we don't need to include
    # heavy dependencies to build documentation
    @require Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" begin
        include("visualization/plot-recipes.jl")
    end
end

# 4.0 Depreciations
@deprecate space_neighbors nearby_ids
@deprecate node_neighbors nearby_positions
@deprecate get_node_contents ids_in_position
@deprecate get_node_agents agents_in_position
@deprecate pick_empty random_empty
@deprecate find_empty_nodes empty_positions
@deprecate has_empty_nodes has_empty_positions
@deprecate nodes positions

# Predefined models
include("models/Models.jl")
export Models

# Update message:
display_update = true
version_number = "4.2"
update_name = "update_v$(version_number)"

if display_update
    if !isfile(joinpath(@__DIR__, update_name))
        printstyled(
            stdout,
            """
            \nUpdate message: Agents v$(version_number)

            We now have full support for pathfinding, using the A* algorithm, in `GridSpace`!
            See the "Path-finding" section in the online docs!
            """;
            color = :light_magenta,
        )
        touch(joinpath(@__DIR__, update_name))
    end
end

end # module


text:
using SuperFastGraphs
using Base.Test
using LightGraphs

# g = loadgraph("./testdata/barabasi_example.lg", "graph")
g = loadgraph("./testdata/karate.lg", "graph")
g = loadgraph("./testdata/internet.lg", "graph")
# g = barabasi_albert(1000, 100)

# g = Graph(5)

# add_edge!(g, 1, 2)
# add_edge!(g, 1, 3)
# add_edge!(g, 1, 4)
# add_edge!(g, 1, 5)
# add_edge!(g, 2, 4)
# add_edge!(g, 2, 3)

tic()
t = triangles(g)
toc()

tic()
triangleNumberD = triangleCountingDegree(g)
toc()

tic()
nt = triangleForward(g)
toc()

println(nt == sum(t)/3)

max_stars_repo_name: mameli/SuperFastGraphs

max_issues_repo_name: mameli/SuperFastGraphs

max_forks_repo_name: mameli/SuperFastGraphs


text:
using BilevelExperiments
using Test

@testset "BilevelExperiments.jl" begin
    # Write your tests here.
end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
using BinDeps

@BinDeps.setup

enable_sensitivities = true

if enable_sensitivities
    cvodes = library_dependency("libsundials_cvodes")
    idas = library_dependency("libsundials_idas")
else
    cvode = library_dependency("libsundials_cvode")
    ida = library_dependency("libsundials_ida")
end
kinsol = library_dependency("libsundials_kinsol")
nvecserial = library_dependency("libsundials_nvecserial")
sundialslibs = enable_sensitivities ? [cvodes, idas, kinsol, nvecserial] : [cvode, ida, kinsol, nvecserial]

sundialsver = "sundials-2.5.0"
provides(Sources, URI("http://ftp.mcs.anl.gov/pub/petsc/externalpackages/$sundialsver.tar.gz"), sundialslibs)

provides(Binaries, URI("https://cache.julialang.org/https://bintray.com/artifact/download/tkelman/generic/$sundialsver.7z"),
    sundialslibs, unpacked_dir="usr$WORD_SIZE/bin", os = :Windows)

prefix = joinpath(BinDeps.depsdir(sundialslibs[1]),"usr")
srcdir = joinpath(BinDeps.depsdir(sundialslibs[1]),"src",sundialsver)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(sundialslibs[1])
        @build_steps begin
            ChangeDirectory(srcdir)
            `./configure --prefix=$prefix --enable-shared`
            `make install`
        end
    end), sundialslibs)

if enable_sensitivities
@BinDeps.install Dict(:libsundials_cvodes => :libsundials_cvodes,
                  :libsundials_idas => :libsundials_idas,
                  :libsundials_kinsol => :libsundials_kinsol,
                  :libsundials_nvecserial => :libsundials_nvecserial)
else
@BinDeps.install Dict(:libsundials_cvode => :libsundials_cvode,
                  :libsundials_ida => :libsundials_ida,
                  :libsundials_kinsol => :libsundials_kinsol,
                  :libsundials_nvecserial => :libsundials_nvecserial)
end


text:
include("testFunctions.jl")
include("polyFunctions.jl")
include("../src/algorithm.jl")

@testset "countZeros counts zeros" begin
    @test_throws ErrorException countZeros(-2 - 5im, 2 - 3im, poly10, 0.01, 0.3)
    @test countZeros(-2 + 2im, 2 - 2im, poly1, 0.1, 0.5) == 2
    @test countZeros(-2 + 2im, 2 - 2im, poly2, 0.1, 0.5) == 5
    @test countZeros(-2 + 2im, 2 - 2im, poly3, 0.01, 0.5) == 2
    @test countZeros(-2 + 2im, 3 - 2im, poly4, 0.01, 0.5) == 2
    @test countZeros(-2 + 2im, 2 - 2im, poly5, 0.01, 0.5) == 2
    @test countZeros(-2 + 5im, 2 - 2im, poly6, 0.01, 0.6) == 1
    @test countZeros(-2 + 2im, 2 - 2im, poly7, 0.01, 0.5) == 2
    @test countZeros(-2 + 2im, 2 - 2im, poly8, 0.01, 0.8) == 3
    @test countZeros(-4 + 5im, 4 - 2im, poly9, 0.01, 0.3) == 3
    @test countZeros(-2 + 4im, 2 - 3im, poly10, 0.01, 0.1) == 2
    @test countZeros(-2 + 4im, 5 - 3im, poly10, 0.01, 0.1) == 3
    @test countZeros(-2 + 5im, 3 - 3im, poly11, 0.01, 0.1) == 3
end;

# Randomised testing

pass = randomised_countzeros(100)
@testset "Randomised testing for countZeros" begin @test pass == 0 end


text:
function maximize_bellman(state_variables_in)
"""
Right hand side of the Bellman equation that gets maximized at each collocation
node
"""

    ## Pull variables out of stateVariables

    # Current states
    states = state_variables_in.states

    # Controls [consumption, abatement cost]
    controls = state_variables_in.controls

    # Exogenous variables
    exog_states = state_variables_in.exog_states

    # Parameters
    mparams = state_variables_in.mparams

    # Approximation space
    space = state_variables_in.space

    # Boolean options
    opts = state_variables_in.opts

    # Terminal approximation space
    space_terminal = state_variables_in.space_terminal

    # Current node
    node = state_variables_in.node

    # Compute flow utility + expected continuation value for learning model
    function value_func(controls, grad)

        cons = controls[1]
        cost = controls[2]

        # Abatement rate
        abatement = (cost ./ exog_states[:Psi]).^(1 ./ mparams[:a2])

        # Next period's states
        g = zeros(mparams[:nqnodes], mparams[:nqnodes], length(states))

        # Cumulative emissions
        g[:, :, 2] .= states[2] .+ exog_states[:emint].*(1 .- abatement[1]).*states[1].^mparams[:kappa] .+ exog_states[:Bs]

        # Damage quadrature
        dams, total_weight = qnwlogn(mparams[:nqnodes],
            states[3]*log.(mparams[:ccr]*g[1, 1, 2]),
            states[4]*log.(mparams[:ccr]*g[1, 1, 2])^2)

        # Fractional net output
        g[:, :, 5] .= 1 ./ (ones(size(dams, 1), size(mparams[:coeff], 1)) + dams*mparams[:coeff]')

        # d_2 mean
        g[:, :, 3] .= ((mparams[:coeff_scale])*states[3] .+
            states[4].*log.(mparams[:ccr]*g[1, 1, 2]).*(log.(1 ./ g[:, :, 5] .- 1) .- mparams[:coeff_loc]))./
            (mparams[:coeff_scale] .+ states[4]*log.(mparams[:ccr]*g[1, 1, 2])^2)

        # d_2 variance
        g[:, :, 4] .= states[4]*(mparams[:coeff_scale])./
            (mparams[:coeff_scale] .+ states[4]*log.(mparams[:ccr]*g[1, 1, 2])^2)

        # Net output
        exog_states[:net_output] = states[1].^(mparams[:kappa]).*states[5]

        # Capital
        g[:, :, 1] .= (states[1].*(1 .- mparams[:deltak]) .+ (exog_states[:net_output]*(1 .- cost[1]) .- cons[1])) .* exog_states[:effkgrowth]

        # Reshape next period states for evaluating continuation value
        g = reshape(g, mparams[:nqnodes]^2, size(g, 3))

        # Evaluate continuation value
        if mparams[:year] .== mparams[:time_horizon]
            vn,eval_out = funeval(space_terminal[:c_cont], space_terminal, [g[:, 1] g[:, 2] g[:, 5]])
            vn /= mparams[:scaling_factor]
        else
            vn,eval_out = funeval(space[:c_cont] ,space, g)
        end

        # Reshape for taking expectations
        vn = reshape(vn, mparams[:nqnodes], mparams[:nqnodes])

        if !opts[:robust]

            # Expected continuation value
            Ev = total_weight'*vn*total_weight

            # Flow utility
            f = (cons[1]^mparams[:rho])/mparams[:rho]/mparams[:scaling_factor]

            # Total utility
            v = (f .+ exog_states[:betaeff]*Ev)

        else

            # Center value for risk operator about zero
            rc_center = -1 .* total_weight'*vn*total_weight

            # Expectation over d_1,omega_{t+1} inside log, expectation over d_2 outside log
            # Undo centering: .- exog_states[:betaeff]*rc_center
            Ev = -mparams[:theta_robust]*total_weight'*
                log.((exp.(-exog_states[:betaeff]*(vn .+ rc_center)/mparams[:theta_robust]))*total_weight) .- exog_states[:betaeff]*rc_center

            # Flow utility
            f = (cons[1]^mparams[:rho])/mparams[:rho]/mparams[:scaling_factor]

            # Total utility
            v = (f .+ Ev)

        end

        return v[1]

    end

    # Compute flow utility .+ expected continuation value for non-learning model
    function value_func_unc(controls, grad)

        cons = controls[1]
        cost = controls[2]

        # Abatement rate
        abatement = (cost./exog_states[:Psi]).^(1 ./ mparams[:a2])

        # Next period's states
        g = zeros(mparams[:nqnodes], mparams[:nqnodes], length(states))

        # Cumulative emissions
        g[:, :, 2] .= states[2] .+ exog_states[:emint].*(1 .- abatement[1]).*states[1].^mparams[:kappa] .+ exog_states[:Bs]

        if ~opts[:perfect_info]

            # Damage quadrature
            dams, total_weight = qnwlogn(mparams[:nqnodes],
                mparams[:exp_mean]*log.(mparams[:ccr]*g[1, 1, 2]),
                mparams[:exp_var]*log.(mparams[:ccr]*g[1, 1, 2])^2)

            # Fractional net output
            g[:, :, 3] .= 1.0./(ones(size(dams,1),size(mparams[:coeff],1))+dams*mparams[:coeff]')

        else

            total_weight = ones(mparams[:nqnodes], 1)/mparams[:nqnodes]

            # Fractional net output
            g[:, :, 3] .= 1 ./ (1 + exp(mparams[:coeff_loc])*(mparams[:ccr]*g[1, 1, 2])^mparams[:exp_mean])

        end

        # Net output
        exog_states[:net_output] = states[1].^(mparams[:kappa]).*states[3]

        # Capital
        g[:, :, 1] .= (states[1].*(1 .- mparams[:deltak]) .+ (exog_states[:net_output]*(1 .- cost[1]) .- cons[1])) .* exog_states[:effkgrowth]

        # Reshape next period states for evaluating continuation value
        g = reshape(g, mparams[:nqnodes]^2, size(g, 3))

        # Evaluate continuation value
        if mparams[:year] .== mparams[:time_horizon]
            vn,eval_out = funeval(space_terminal[:c_cont], space_terminal, g)
            vn /= mparams[:scaling_factor]
        else
            vn,eval_out = funeval(space[:c_cont], space, g)
        end

        # Reshape for taking expectations
        vn = reshape(vn, mparams[:nqnodes], mparams[:nqnodes])

        if !opts[:robust]

            # Expected continuation value
            Ev = total_weight'*vn*total_weight

            # Flow utility
            f = (cons[1]^mparams[:rho])/mparams[:rho]/mparams[:scaling_factor]

            # Total utility
            v = (f .+ exog_states[:betaeff]*Ev)

        else

            # Center value for risk operator
            rc_center = -1 .* total_weight'*vn*total_weight

            # Expectation over d_1,omega_{t+1} inside log, expectation over d_2 outside log
            # Undo centering: .- exog_states[:betaeff]*rc_center
            Ev = -mparams[:theta_robust]*total_weight'*
                log.((exp.(-exog_states[:betaeff]*(vn .+ rc_center)/mparams[:theta_robust]))*total_weight) .- exog_states[:betaeff]*rc_center

            # Flow utility
            f = (cons[1]^mparams[:rho])/mparams[:rho]/mparams[:scaling_factor]

            # Total utility
            v = (f .+ Ev)
        end

        return v[1]

    end

    # Initialize optimizer
    opt = Opt(:LN_COBYLA,2)

    # Bounds on consumption (~0, net output) and abatement cost (0, cost of full abatement)
    lower_bounds!(opt, [1e-10, 0.])

    # Initial guess needs to be inside the bounds or else NLopt throws an error
    if opts[:learn_e]
        upper_bounds!(opt, [states[5]*states[1].^mparams[:kappa], exog_states[:Psi]])
        controls[1] = max(min(states[5]*states[1].^mparams[:kappa], controls[1]), 1e-5)
        controls[2] = max(min(exog_states[:Psi], controls[2]), 0.)
    else
        upper_bounds!(opt, [states[3]*states[1].^mparams[:kappa], exog_states[:Psi]])
        controls[1] = max(min(states[3]*states[1].^mparams[:kappa], controls[1]), 1e-5)
        controls[2] = max(min(exog_states[:Psi], controls[2]), 0.)
    end

    function resource_constraint(xx,grad)
    """
    Inequality constraint requiring that consumption and abatement are less than
    net output. Returns negative (required by NLopt) net output after consumption
    and abatement.
    """

        if opts[:learn_e]
            resource_left = xx[1] .+ (xx[2] - 1)*states[1].^mparams[:kappa].*states[5]
        else
            resource_left = xx[1] .+ (xx[2] - 1)*states[1].^mparams[:kappa].*states[3]
        end

        return resource_left

    end

    # Set resource constraint
    inequality_constraint!(opt, (xx,g) -> resource_constraint(xx,g), 1e-4)

    # Solution tolerances
    if opts[:sim]
        ftol_abs!(opt, 1e-13)
        ftol_rel!(opt, 1e-13)
        xtol_rel!(opt, 1e-13)
        xtol_abs!(opt, 1e-13)
    else
        ftol_abs!(opt, 1e-8)
        ftol_rel!(opt, 1e-8)
        xtol_rel!(opt, 1e-8)
        xtol_abs!(opt, 1e-8)
    end

    # Set max evaluations
    maxeval!(opt, 2000)

    # Declare maximization problem
    if opts[:learn_e]
        max_objective!(opt, value_func)
    else
        max_objective!(opt, value_func_unc)
    end

    # NLopt does a bad job handling initial guesses where the constraints are violated
    # so if the initial guess violates the constraint keep reducing consumption
    # until it doesnt
    while resource_constraint(controls, [0.]) > 0
        controls[1] /= 2
    end

    # Maximize right hand side of Bellman
    maxf, maxx, ret = optimize(opt, controls)

    # Return argmaxs and maximum value for polynomial fitting
    return maxx, maxf[1]

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
using Dolo

model = Model("/home/pablo/.julia/dev/Dolo/examples/models/rbc.yaml")

model

@time perturb(model)

fff = model.factories[:arbitrage]

fff

using OrderedCollections
using SymEngine
import Dolang
using StaticArrays

ff2 = Dolang.FlatFunctionFactory(OrderedCollections.OrderedDict{Symbol,Union{Expr, Number, Symbol}}(:_out_1_ => :(_chi_ * _n__0_ ^ _eta_ * _c__0_ ^ _sigma_ - ((1 - _alpha_) * exp(_z__0_) * _k__0_ ^ _alpha_ * _n__0_ ^ (1 - _alpha_)) / _n__0_),:_out_2_ => :(1 - _beta_ * (_c__0_ / _c__1_) ^ _sigma_ * ((1 - _delta_) + _rk__1_))), OrderedCollections.OrderedDict(:m => [:_z__0_],:s => [:_k__0_],:x => [:_n__0_, :_i__0_],:M => [:_z__1_],:S => [:_k__1_],:X => [:_n__1_, :_i__1_],:p => [:_beta_, :_sigma_, :_eta_, :_chi_, :_delta_, :_alpha_, :_rho_, :_zbar_, :_sig_z_]), [:_out_1_, :_out_2_], OrderedCollections.OrderedDict{Symbol,Union{Expr, Number, Symbol}}(:_y__1_ => :(exp(_z__1_) * _k__1_ ^ _alpha_ * _n__1_ ^ (1 - _alpha_)),:_y__0_ => :(exp(_z__0_) * _k__0_ ^ _alpha_ * _n__0_ ^ (1 - _alpha_)),:_c__0_ => :(_y__0_ - _i__0_),:_rk__1_ => :((_alpha_ * _y__1_) / _k__1_),:_c__1_ => :(_y__1_ - _i__1_)), :arbitrage)

l = ([0.0], [9.354978290145986], [0.33, 0.23387445725364966], [0.99, 5.0, 1.0, 23.95785990938192, 0.025, 0.33, 0.8, 0.0, 0.016])
m, s, x, p = [SVector(e...) for e in l]

code = Dolang.gen_kernel(fff, [0,3])
arb = eval(code)

res, f_x = arb(m,s,x,m,s,x,p)

# true result
# array([[ 34.95244649, -13.27057815],
#         [-13.27057912,   6.56871655]]),

# ff2.preamble
# ff2.equations

# all_eqs = OrderedDict{Symbol, Expr}()

# for (k,v) in ff2.preamble
#     all_eqs[k] = v
# end
# for (k,v) in ff2.equations
#     all_eqs[k] = v
# end

# all_eqs

# import Dolang: diff_symbol
# using Dolang: diff_symbol, deriv

# function add_gradients(all_eqs, gradient_components)

#     new_eqs = OrderedDict{Symbol, Union{Expr, Int64}}()
#     for g in gradient_components
#         k = diff_symbol(g,g)
#         new_eqs[k] = 1
#     end

#     for (sym, eq) in all_eqs

#         new_eqs[sym] = eq
#         # symbols the equation depends on
#         deps = Dolang.list_parameters(eq)

#         for θ in gradient_components

#             deq_terms = []

#             for a in deps

#                 # compute dy/da * da/dθ

#                 if a == θ
#                     push!(deq_terms, ( deriv(eq, a), 1) )
#                 elseif !(a in keys(new_eqs))  # I should keep a list of defined variables
#                     # nothing to do
#                 else
#                     da_dθ = diff_symbol(a, θ)
#                     if da_dθ in keys(new_eqs) # otherwise it is 0
#                         dy_da = diff_symbol(sym, a)
#                         if !(dy_da in keys(new_eqs))
#                             new_eqs[dy_da] = deriv(eq, a)
#                         end
#                         push!(deq_terms, (dy_da, da_dθ) )
#                     end
#                 end
#             end
#             if length(deq_terms)>0
#                 if length(deq_terms)==1
#                     rhs = Expr(:call, :*, deq_terms[1]...)
#                 else
#                     rhs = Expr(:call, :+, [Expr(:call, :*, a, b) for (a,b) in deq_terms]...)
#                 end
#                 new_eqs[diff_symbol(sym, θ)] = rhs
#             end
#         end
#     end

#     return new_eqs

# end

# new_eqs = add_gradients(all_eqs, [:_n__1_])

# for (k,v) in new_eqs

#     println(k, " : ", v)
# end

max_issues_repo_name: EconForge/Dolang


text:
using Documenter, IntervalArithmetic

makedocs(
    modules = [IntervalArithmetic],
    format = :html,
    sitename = "IntervalArithmetic",
    pages = [
        "Package" => "index.md",
        "Interval Arithmetic" => "intro.md",
        "Constructing intervals" => "construction.md",
        "Basic usage" => "usage.md",
        "Decorations" => "decorations.md",
        "Multi-dimensional boxes" => "multidim.md",
        "Rounding" => "rounding.md",
        "API" => "api.md"
    ]
)

deploydocs(
    repo   = "github.com/JuliaIntervals/IntervalArithmetic.jl.git",
    target = "build",
    julia  = "1.0",
    osname = "linux",
    deps   = nothing,
    make   = nothing
)


text:
using .Plots

function preprocess(name, plot::Plots.Plot, data)
    push!(data, name=>plot)
    return data
end

function preprocess(name, plots::AbstractArray{<:Plots.Plot}, data)
    for (i, plot)=enumerate(plots)
        preprocess(name*"/$i", plot, data)
    end
    return data
end

function log_plot(lg::CLogger, name::AbstractString, obj::Plots.Plot; step=nothing)
    mktempdir(;prefix="jlcom_") do dirname
        fname = joinpath(dirname, name, "plot.png")
        savefig(obj, fname)
        log_image(lg, name, fname; step=step, copy_to_tmp=false)
    end
end

process(lg::CLogger, name::AbstractString, obj::Plots.Plot, step::Int) = log_plot(lg, name, obj; step=step)


text:
# Adjoints needed in polar formulation
is_constraint(::typeof(active_power_generation)) = true
size_constraint(polar::PolarForm, ::typeof(active_power_generation)) = get(polar, PS.NumberOfGenerators())

function adjoint!(
    polar::PolarForm,
    pbm::AutoDiff.TapeMemory{F, S, I},
    pg, ∂pg,
    vm, ∂vm,
    va, ∂va,
    pinj, ∂pinj,
) where {F<:typeof(active_power_generation), S, I}
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    nref = PS.get(polar.network, PS.NumberOfSlackBuses())
    index_pv = polar.indexing.index_pv
    index_ref = polar.indexing.index_ref
    pv2gen = polar.indexing.index_pv_to_gen
    ref2gen = polar.indexing.index_ref_to_gen

    ngen = get(polar, PS.NumberOfGenerators())
    ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix())

    fill!(∂vm, 0.0)
    fill!(∂va, 0.0)
    fill!(∂pinj, 0.0)
    ev = adj_active_power_kernel!(polar.device)(
        ∂pg,
        vm, ∂vm,
        va, ∂va,
        ∂pinj,
        index_pv, index_ref, pv2gen, ref2gen,
        ybus_re.nzval, ybus_re.colptr, ybus_re.rowval, ybus_im.nzval,
        ndrange=ngen,
        dependencies=Event(polar.device)
    )
    wait(ev)
    return
end
#

function pullback_objective(polar::PolarForm)
    return AutoDiff.TapeMemory(
        active_power_generation,
        AdjointStackObjective(polar),
        nothing,
    )
end

function cost_production(polar::PolarForm, buffer::PolarNetworkState)
    pg = buffer.pg
    ngen = PS.get(polar, PS.NumberOfGenerators())
    coefs = polar.costs_coefficients
    c2 = @view coefs[:, 2]
    c3 = @view coefs[:, 3]
    c4 = @view coefs[:, 4]
    # Return quadratic cost
    # NB: this operation induces three allocations on the GPU,
    #     but is faster than writing the sum manually
    return sum(c2 .+ c3 .* pg .+ c4 .* pg.^2)
end

function adjoint_cost!(polar::PolarForm, ∂f, pg)
    coefs = polar.costs_coefficients
    c3 = @view coefs[:, 3]
    c4 = @view coefs[:, 4]
    # Return adjoint of quadratic cost
    ∂f .= c3 .+ 2.0 .* c4 .* pg
    return
end

function hessian_cost!(polar, ∂²f)
    coefs = polar.costs_coefficients
    c4 = @view coefs[:, 4]
    ∂²f .= 2.0 .* c4
    return
end

function put(
    polar::PolarForm{T, IT, VT, MT},
    ::PS.Generators,
    ::PS.ActivePower,
    pbm::AutoDiff.TapeMemory,
    buffer::PolarNetworkState
) where {T, IT, VT, MT}

    obj_autodiff = pbm.stack

    adj_pg = obj_autodiff.∂pg
    adj_x = obj_autodiff.∇fₓ
    adj_u = obj_autodiff.∇fᵤ
    adj_vmag = obj_autodiff.∂vm
    adj_vang = obj_autodiff.∂va
    adj_pinj = obj_autodiff.∂pinj

    # Adjoint of active power generation
    adjoint!(polar, pbm,
        buffer.pg, adj_pg,
        buffer.vmag, adj_vmag,
        buffer.vang, adj_vang,
        buffer.pinj, adj_pinj
    )

    # Adjoint w.r.t. x and u
    fill!(adj_x, 0.0)
    fill!(adj_u, 0.0)
    adjoint_transfer!(polar, adj_u, adj_x, adj_vmag, adj_vang, adj_pinj)

    return
end

function gradient_objective!(polar::PolarForm, ∂obj::AutoDiff.TapeMemory, buffer::PolarNetworkState)
    ∂pg = ∂obj.stack.∂pg
    adjoint_cost!(polar, ∂pg, buffer.pg)
    put(polar, PS.Generators(), PS.ActivePower(), ∂obj, buffer)
    return
end

function hessian_prod_objective!(
    polar::PolarForm,
    ∇²f::AutoDiff.Hessian, adj_obj::AutoDiff.TapeMemory,
    hv::AbstractVector,
    ∂²f::AbstractVector, ∂f::AbstractVector,
    buffer::PolarNetworkState,
    tgt::AbstractVector,
)
    nx = get(polar, NumberOfState())
    nu = get(polar, NumberOfControl())
    # Indexing of generators
    pv2gen = polar.indexing.index_pv_to_gen ; npv = length(pv2gen)
    ref2gen = polar.indexing.index_ref_to_gen ; nref = length(ref2gen)

    ∇f = adj_obj.stack

    # Remember that
    # ```math
    # ∂f = (∂f / ∂pg) * ∂pg
    # ```
    # Using the chain-rule, we get
    # ```math
    # ∂²f = (∂f / ∂pg) * ∂²pg + ∂pg' * (∂²f / ∂²pg) * ∂pg
    # ```

    ## Step 1: evaluate (∂f / ∂pg) * ∂²pg
    ∂pg_ref = @view ∂f[ref2gen]
    AutoDiff.adj_hessian_prod!(polar, ∇²f, hv, buffer, ∂pg_ref, tgt)

    ## Step 2: evaluate ∂pg' * (∂²f / ∂²pg) * ∂pg
    # Compute adjoint w.r.t. ∂pg_ref
    # copyto!(∇f.∂pg, ∂²f)
    ∇f.∂pg .= 0.0

    ∇f.∂pg[ref2gen] .= 1.0
    put(polar, PS.Generators(), PS.ActivePower(), adj_obj, buffer)

    tx = view(tgt, 1:nx)
    tu = view(tgt, 1+nx:nx+nu)
    tpg = view(tgt, nx+nu-npv+1:nx+nu)

    hvx = view(hv, 1:nx)
    hvu = view(hv, 1+nx:nx+nu)
    hvpg = view(hv, nx+nu-npv+1:nx+nu)

    ∇pgₓ = ∇f.∇fₓ
    ∇pgᵤ = ∇f.∇fᵤ

    scale_x = dot(∇pgₓ, ∂²f[ref2gen[1]], tx)
    scale_u = dot(∇pgᵤ, ∂²f[ref2gen[1]], tu)
    α = scale_x + scale_u
    # Contribution of slack node
    axpy!(α, ∇pgₓ, hvx)
    axpy!(α, ∇pgᵤ, hvu)
    # Contribution of PV nodes (only through power generation)
    axpy!(1.0, tpg .* ∂²f[pv2gen], hvpg)

    return nothing
end

#=
    Special function to compute Hessian of ProxAL's objective.
    This avoid to reuse intermediate results, for efficiency.

    This function collect the contribution of the state and
    the control (in `hv_xu`) and the contribution of the slack
    variable (`hv_s`).

    For ProxAL, we have:
    H = [ H_xx  H_ux  J_x' ]
        [ H_xu  H_uu  J_u' ]
        [ J_x   J_u   ρ I  ]

    so, if `v = [v_x; v_u; v_s]`, we get

    H * v = [ H_xx v_x  +   H_ux v_u  +  J_x' v_s ]
            [ H_xu v_x  +   H_uu v_u  +  J_u' v_s ]
            [  J_x v_x  +    J_u v_u  +   ρ I     ]

=#
function hessian_prod_objective_proxal!(
    polar::PolarForm,
    ∇²f::AutoDiff.Hessian, adj_obj::AutoDiff.TapeMemory,
    hv_xu::AbstractVector,
    hv_s::AbstractVector,
    ∂²f::AbstractVector, ∂f::AbstractVector,
    buffer::PolarNetworkState,
    tgt::AbstractVector,
    vs::AbstractVector,
    ρ::Float64,
    has_slack::Bool,
)
    nx = get(polar, NumberOfState())
    nu = get(polar, NumberOfControl())
    ngen = get(polar, PS.NumberOfGenerators())

    # Stack
    ∇f = adj_obj.stack
    # Indexing of generators
    pv2gen = polar.indexing.index_pv_to_gen ; npv = length(pv2gen)
    ref2gen = polar.indexing.index_ref_to_gen ; nref = length(ref2gen)

    # Split tangent wrt x part and u part
    tx = @view tgt[1:nx]
    tu = @view tgt[1+nx:nx+nu]

    #= BLOCK UU
        [ H_xx v_x  +   H_ux v_u ]
        [ H_xu v_x  +   H_uu v_u ]
    =#
    ## Step 1: evaluate (∂f / ∂pg) * ∂²pg
    ∂pg_ref = @view ∂f[ref2gen]
    AutoDiff.adj_hessian_prod!(polar, ∇²f, hv_xu, buffer, ∂pg_ref, tgt)

    ## Step 2: evaluate ∂pg' * (∂²f / ∂²pg) * ∂pg
    # Compute adjoint w.r.t. ∂pg_ref
    ∇f.∂pg .= 0.0
    ∇f.∂pg[ref2gen] .= 1.0
    put(polar, PS.Generators(), PS.ActivePower(), adj_obj, buffer)

    ∇pgₓ = ∇f.∇fₓ
    ∇pgᵤ = ∇f.∇fᵤ

    @views begin
        tpg = tgt[nx+nu-npv+1:end]

        scale_x = dot(∇pgₓ, ∂²f[ref2gen[1]], tx)
        scale_u = dot(∇pgᵤ, ∂²f[ref2gen[1]], tu)
        # Contribution of slack node
        hv_xu[1:nx]       .+= (scale_x + scale_u) .* ∇pgₓ
        hv_xu[1+nx:nx+nu] .+= (scale_x + scale_u) .* ∇pgᵤ
        # Contribution of PV nodes (only through power generation)
        hv_xu[nx+nu-npv+1:end] .+= ∂²f[pv2gen] .* tpg
    end

    if has_slack
        @views begin
            #= BLOCK SU
                [  J_x' v_s ]
                [  J_u' v_s ]
                [  J_x v_x + J_u v_u ]
            =#
            # 1. Contribution of slack node
            hv_s[ref2gen] .+= - ρ * (dot(∇pgᵤ, tu) + dot(∇pgₓ, tx))
            hv_xu[1:nx] .-= ρ * ∇pgₓ .* vs[ref2gen]
            hv_xu[1+nx:nx+nu] .-= ρ * ∇pgᵤ .* vs[ref2gen]

            # 2. Contribution of PV node
            hv_s[pv2gen] .-= ρ .* tu[nref+npv+1:end]
            hv_xu[nx+nref+npv+1:end] .-= ρ .* vs[pv2gen]

            #= BLOCK SS
                 ρ I
            =#
            # Hessian w.r.t. slack
            hv_s .+= ρ .* vs
        end
    end

    return nothing
end


text:
module Units

import Unitful
using Unitful: @unit, @derived_dimension,𝐋

# Re-export units
import Unitful: nm, m, cm, s, kg, g, eV, keV, MeV, c, u, b, C, J, A, H, Hz, rad, T, V, Ω, W, K, Pa, angstrom, Å, cal, atm, Torr, mol, lm
import Unitful: k, h, c0

# Define New Units
@unit Gauss "Gauss"    Gauss       (1//10_000)T   true
@unit kB     "kB"    Boltzmann          k         false
@unit hc     "hc"   PlanckLight        h*c0       false

@derived_dimension NumberDensity 𝐋^-3

export nm, m, cm, s, kg, g, eV, keV, MeV, c, u, b, C, J, A, Hz, H, rad, T, V, Ω, W, K, Pa, angstrom, Å, cal, atm, Torr, mol, lm
export Gauss, kB, hc, NumberDensity

#Register Units
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(Units)
end

end


text:
using ShaderToy, GLAbstraction, FileIO, ColorTypes
shadertoy("submerged.frag", Dict{Symbol, Any}(
	:iChannel0 => Texture(rand(Float32, 64,64), x_repeat=:repeat, minfilter=:linear),
))


text:
#!/usr/bin/env julia

using BinaryBuilder, LibGit2, SHA

function build_tarballs_path(dep_name::String, base::String = @__DIR__)
    # Strip _jll from the end of the dep name, if it has it
    if endswith(dep_name, "_jll")
        dep_name = dep_name[1:end-4]
    end
    return joinpath(base, uppercase(dep_name[1:1]), dep_name, "build_tarballs.jl")
end

build_tarballs = build_tarballs_path(get(Sys.ARGS, 1, ""))

if !isfile(build_tarballs)
    println(stderr, "usage: recursively_regenerate_jlls.jl <jll name>")
    exit(1)
end

function recursively_collect_dependencies(build_tarballs::String, dependencies = Set{String}())
    # Generate meta.json for this script, read it in
    @info("Parsing $(build_tarballs)")
    meta_json = String(read(`$(Base.julia_cmd()) --project=$(Base.active_project()) $(build_tarballs) --meta-json`))
    meta_obj = BinaryBuilder.JSON.parse(meta_json)
    for dep in meta_obj["dependencies"]
        if !(dep["name"] ∈ dependencies)
            dep_build_tarballs = build_tarballs_path(dep["name"])
            if !isfile(dep_build_tarballs)
                @warn("Unable to find build_tarballs.jl file for dependency \"$(dep["name"])\"")
            else
                push!(dependencies, dep["name"])
                recursively_collect_dependencies(dep_build_tarballs, dependencies)
            end
        end
    end
    return dependencies
end

# This borrows heavily from BinaryBuilder.yggdrasil_deploy()
function open_jll_bump_pr(dep_name::String)
    @info("Opening JLL-bumping PR for $(dep_name)")
    gh_auth = BinaryBuilder.github_auth(;allow_anonymous=false)
    fork = BinaryBuilder.GitHub.create_fork("JuliaPackaging/Yggdrasil"; auth=gh_auth)

    mktempdir() do tmp
        repo = LibGit2.clone(BinaryBuilder.get_yggdrasil(), tmp)
        build_tarballs = build_tarballs_path(dep_name, tmp)
        if !isfile(build_tarballs)
            throw(ArgumentError("Invalid dep_name \"$(dep_name)\", no file $(build_tarballs)"))
        end

        # We're going to need to alter the build_tarballs.jl file.
        # We will just add a newline onto the end if it doesn't already
        # have two, otherwise we remove one.
        build_tarballs_content = String(read(build_tarballs))
        if build_tarballs_content[end-1:end] != "\n\n"
            build_tarballs_content = string(build_tarballs_content, "\n")
        else
            build_tarballs_content = build_tarballs_content[1:end-1]
        end

        recipe_hash = bytes2hex(sha256(build_tarballs_content))[end-3:end]
        branch_name = "jllbump/$(dep_name)_$(recipe_hash)"
        LibGit2.branch!(repo, branch_name)

        # Write out the modified build_tarballs.jl file
        open(build_tarballs, "w") do io
            write(build_tarballs, build_tarballs_content)
        end

        # Commit modified build_tarballs.jl file
        LibGit2.add!(repo, build_tarballs_path(dep_name, ""))
        LibGit2.commit(repo, "Bump JLL wrapper: $(dep_name)\n\n[skip build]")
        creds = LibGit2.UserPasswordCredential(
            dirname(fork.full_name),
            deepcopy(gh_auth.token),
        )
        try
            LibGit2.push(
                repo,
                refspecs=["+HEAD:refs/heads/$(branch_name)"],
                remoteurl="https://github.com/$(fork.full_name).git",
                credentials=creds,
                # This doesn't work :rage: instead we use `+HEAD:` at the beginning of our
                # refspec: https://github.com/JuliaLang/julia/issues/23057
                #force=true,
            )
        finally
            Base.shred!(creds)
        end

		params = Dict(
			"base" => "master",
			"head" => "$(dirname(fork.full_name)):$(branch_name)",
			"maintainer_can_modify" => true,
			"title" => "JLL bump: $(dep_name)",
			"body" => """
			This pull request bumps the JLL version of $(dep_name).
			It was generated via the `recursively_regenerate_jlls.jl` script.
			"""
		)
		pr = BinaryBuilder.create_or_update_pull_request("JuliaPackaging/Yggdrasil", params, auth=gh_auth)
        println("https://github.com/$(fork.full_name)/pull/new/$(BinaryBuilder.HTTP.escapeuri(branch_name))?expand=1")
    end
end

deps = recursively_collect_dependencies(build_tarballs)
push!(deps, Sys.ARGS[1])
println("Discovered dependencies:")
println(deps)

if BinaryBuilder.yn_prompt(BinaryBuilder.WizardState(), "Open JLL-bumping PRs?", :y) == :y
	for dep in deps
		open_jll_bump_pr(dep)
	end
end


text:
struct UnimplementedFeaturization <: SeaPearl.AbstractFeaturization end
struct TestFeatState{F,TS} <: SeaPearl.FeaturizedStateRepresentation{F,TS} end

struct TestFeaturization <: SeaPearl.AbstractFeaturization end

function SeaPearl.featurize(::TestFeatState{TestFeaturization})
    return [1, 2, 3]
end

@testset "representation.jl" begin

    include("default/cp_layer/cp_layer.jl")
    include("default/defaulttrajectorystate.jl")
    include("default/defaultstaterepresentation.jl")
    include("tsptw/tsptwtrajectorystate.jl")
    include("tsptw/tsptwstaterepresentation.jl")

    @testset "Unimplemented featurization" begin
        @test_throws ErrorException SeaPearl.featurize(TestFeatState{UnimplementedFeaturization,SeaPearl.DefaultTrajectoryState}())
    end

    @testset "Custom featurization" begin
        @test SeaPearl.featurize(TestFeatState{TestFeaturization,SeaPearl.DefaultTrajectoryState}()) == [1, 2, 3]
    end

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: PHE.mso 250 2007-04-27 16:32:02Z bicca $
#*------------------------------------------------------------------
type PHE_PressureDrop
	PHE_PressureDrop()=begin
		new(
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Channel Pressure Drop",
				:Default=>0.01,
				:Lower=>1E-10,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P^{channel}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Ports Pressure Drop",
				:Default=>0.01,
				:Lower=>1E-10,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P^{ports}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Total Pressure Drop",
				:Default=>0.01,
				:Lower=>1E-10,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			fricfactor(Dict{Symbol,Any}(
				:Brief=>"Friction Factor",
				:Default=>0.05,
				:Lower=>1E-10,
				:Upper=>2000
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Stream Velocity in Channel",
				:Lower=>1E-8,
				:Symbol=>"V^{channel}"
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"Stream Velocity in Ports",
				:Lower=>1E-8,
				:Symbol=>"V^{ports}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Number of  Channels per Pass",
				:Symbol=>"N^{passage}"
			)),
			[:DPchannel,:DPports,:Pdrop,:fi,:Vchannel,:Vports,:Npassage,]
		)
	end
	DPchannel::press_delta
	DPports::press_delta
	Pdrop::press_delta
	fi::fricfactor
	Vchannel::velocity
	Vports::velocity
	Npassage::positive
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PHE_PressureDrop
function atributes(in::PHE_PressureDrop,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
PHE_PressureDrop(_::Dict{Symbol,Any})=begin
	newModel=PHE_PressureDrop()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
### Pedro Cavalcante
### pc_oliveira@id.uff.br

### Implementação da Função de Weierstrass
### Tópicos Especiais de Métodos Estatísticos

function weierstrass(x, a = 0.5, b = 7, eps = 9^(-5), n = 10000)

    if a*b <= (1 + 3*pi/2)
        print("a*b must be greater than 1 + 3*pi/2")
    end

    if a < 0 | a > 1
        print("a must be a positive number smaller than 1")
    end

    if isodd(b) == FALSE
        print("b must be an odd number")
    end

    soma[1] = 0

    for i = 2:n
        if a^n * cos(b^n * pi * x) > eps
            soma[i] = a^n * cos(b^n * pi * x)
        end
    end

    return sum(soma)

end


text:
"""
Filter contains middleware that is able to mutate the request map.
These functions are most commonly used to make accessing data from the request more simple and straight forward.
"""
module Filter

using HttpServer
import HttpCommon: Cookie

"""
# strtodict

converts
"""
function strtodict(str, delim)
    str = split(str, delim)
    str = map((kv) -> split(kv, "="), str)
    str = filter((l) -> length(l) > 1, str)
    [ strip(k) => v for (k, v) in str ]
end

"""
# parse_multipart
"""
function parse_multipart(raw, boundary)
    println(raw)
    raw = convert(ASCIIString, raw)
    raw = split(raw, "\r\n")
    raw = map(strip, raw)
    raw = map((kv) -> split(kv, boundary), raw)
    # str = filter((l) -> length(l) > 1, str)
    println(raw)
    [ strip(k) => v for (k, v) in raw ]
end

"""
# body_todict(app, ctx)

Converts the query from the request to a dictionary. Then stores in to locations
```julia
ctx.data[:request][:data] ## overwriting the raw version
ctx.data[:body] ## for convenience
```
"""
function body_todict(app, ctx)
    if !in(:data, ctx.data[:request] |> keys)
        ctx.data[:data] = Dict()
        return app(ctx)
    end

    if :type in keys(ctx.data[:contenttype])
        raw_data = convert(ASCIIString, ctx.data[:request][:data])

        ctx.data[:body] = if ctx.data[:contenttype][:type] != "multipart"
            strtodict(raw_data, "&")
        else
            parse_multipart(raw_data,
                            ctx.data[:contenttype][:parameters]["boundary"])
        end

        return app(ctx)
    end

    raw_data = convert(ASCIIString, ctx.data[:request][:data])
    data = strtodict(raw_data, "&")
    ctx.data[:body] = data
    app(ctx)
end

# hacky
function contenttype_todict(app, ctx)
    if !in("Content-Type", ctx.data[:request][:headers] |> keys)
        ctx.data[:contenttype] = Dict()
        return app(ctx)
    end

    ct = Dict()

    raw_data = convert(ASCIIString, ctx.data[:request][:headers]["Content-Type"])

    println(raw_data)

    parts = split(raw_data, ";")
    ct[:parameters] = length(parts) > 1 ?
        strtodict(parts[2], ";") :
        Dict()

    types                  = split(parts[1], "/")
    ct[:type]              = types[1]
    ct[:subtype]           = types[2]
    ctx.data[:contenttype] = ct
    app(ctx)
end

"""
# body_todict(app, ctx)

Converts the query from the request to a dictionary. Then stores in to locations
```julia
ctx.data[:request][:cookie] ## overwriting the raw version
ctx.data[:cookie] ## for convenience
```
"""
function cookie_todict(app::Function, ctx)
    if !in("Cookie", keys(ctx.data[:request][:headers]))
        ctx.data[:cookies] = Dict()
        resp = app(ctx)
        return resp
    end
    dough = ctx.data[:request][:headers]["Cookie"]
    cookie = strtodict(dough, "; ")
    ctx.data[:cookies] = cookie

    resp = app(ctx)

    # adding new cookies added from the handler to the response
    for (k, v) in ctx.data[:cookies]
        ctx.response.cookies[string(k)] = Cookie(k, v)
    end
    return resp
end

"""
# query_todict

Converts the query from the request to a dictionary. Then stores in to locations
```julia
ctx.data[:request][:query] ## overwriting the raw version
ctx.data[:query] ## for convenience
```
"""
function query_todict(app::Function, ctx)
    raw_query = ctx.data[:request][:query]
    marshalled = strtodict(raw_query, "&")
    ctx.data[:request][:query] = marshalled
    ctx.data[:query]           = marshalled
    app(ctx)
end

function mime!(ctx, value::AbstractString)
    ctx.response.headers["Content-Type"] = value * "; charset=utf-8"

end

end


text:
int_rules_3_2_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 3.2.1*(f+g*x)^m*(A+B*log(e*((a+b*x)*over*(c+d*x))^n))^p =#
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => ((~a + ~b * ~x) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ ~p) / ~b - ((~B * ~n * ~p * (~b * ~c - ~a * ~d)) / ~b) * Antiderivative((~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ (~p - 1) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => ((~a + ~b * ~x) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ ~p) / ~b - ((~B * ~n * ~p * (~b * ~c - ~a * ~d)) / ~b) * Antiderivative((~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ (~p - 1) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && IGtQ(~p, 0)))
    #= Antiderivative(((~A')+(~B')*log((~e')*(((~a')+(~b')*(~x))^(~n1')*((~c')+(~d')*(~x))^(~n2))^(~n')))^(~p') ,~x) := (a+b*x)*(A+B*log(e*((a+b*x)^n1/(c+d*x)^n1)^n))^p/b - B*n*p*(b*c-a*d)/b*Antiderivative((A+B*log(e*((a+b*x)^n1/(c+d*x)^n1)^n))^(p-1)/( c+d*x),x) <-- FreeQ([a,b,c,d,e,A,B,n],x) && EqQ(n1+n2,0) && GtQ(n1,0) &&  (EqQ(n1,1) || EqQ(n,1)) && NeQ(b*c-a*d,0) && IGtQ(p,0) =#
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~p, ~x) => Unintegrable((~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~p, ~x) => Unintegrable((~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~n, ~p], ~x) && EqQ(~n + ~mn, 0)
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~u / ~v) ^ ~(n'))) ^ ~(p'), ~x) => Antiderivative((~A + ~B * log(~e * (ExpandToSum(~u, ~x) / ExpandToSum(~v, ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~n, ~p], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~u) ^ ~(n') * (~v) ^ ~mn)) ^ ~(p'), ~x) => Antiderivative((~A + ~B * log((~e * ExpandToSum(~u, ~x) ^ ~n) / ExpandToSum(~v, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) / (~(f') + ~(g') * ~x), ~x) => (-(log(-((~b * ~c - ~a * ~d)) / (~d * (~a + ~b * ~x)))) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n))) / ~g + ((~B * ~n * (~b * ~c - ~a * ~d)) / ~g) * Antiderivative(log(-((~b * ~c - ~a * ~d)) / (~d * (~a + ~b * ~x))) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~b * ~f - ~a * ~g, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) / (~(f') + ~(g') * ~x), ~x) => (-(log(-((~b * ~c - ~a * ~d)) / (~d * (~a + ~b * ~x)))) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n))) / ~g + ((~B * ~n * (~b * ~c - ~a * ~d)) / ~g) * Antiderivative(log(-((~b * ~c - ~a * ~d)) / (~d * (~a + ~b * ~x))) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~b * ~f - ~a * ~g, 0)))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) / (~(f') + ~(g') * ~x), ~x) => (-(log((~b * ~c - ~a * ~d) / (~b * (~c + ~d * ~x)))) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n))) / ~g + ((~B * ~n * (~b * ~c - ~a * ~d)) / ~g) * Antiderivative(log((~b * ~c - ~a * ~d) / (~b * (~c + ~d * ~x))) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~d * ~f - ~c * ~g, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) / (~(f') + ~(g') * ~x), ~x) => (-(log((~b * ~c - ~a * ~d) / (~b * (~c + ~d * ~x)))) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n))) / ~g + ((~B * ~n * (~b * ~c - ~a * ~d)) / ~g) * Antiderivative(log((~b * ~c - ~a * ~d) / (~b * (~c + ~d * ~x))) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~d * ~f - ~c * ~g, 0)))
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) / (~(f') + ~(g') * ~x), ~x) => ((log(~f + ~g * ~x) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n))) / ~g - ((~b * ~B * ~n) / ~g) * Antiderivative(log(~f + ~g * ~x) / (~a + ~b * ~x), ~x)) + ((~B * ~d * ~n) / ~g) * Antiderivative(log(~f + ~g * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) / (~(f') + ~(g') * ~x), ~x) => ((log(~f + ~g * ~x) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n))) / ~g - ((~b * ~B * ~n) / ~g) * Antiderivative(log(~f + ~g * ~x) / (~a + ~b * ~x), ~x)) + ((~B * ~d * ~n) / ~g) * Antiderivative(log(~f + ~g * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && NeQ(~b * ~c - ~a * ~d, 0))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))), ~x) => ((~f + ~g * ~x) ^ (~m + 1) * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n))) / (~g * (~m + 1)) - ((~B * ~n * (~b * ~c - ~a * ~d)) / (~g * (~m + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m + 1) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ(~m, -1) && NeQ(~m, -2)))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)), ~x) => ((~f + ~g * ~x) ^ (~m + 1) * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n))) / (~g * (~m + 1)) - ((~B * ~n * (~b * ~c - ~a * ~d)) / (~g * (~m + 1))) * Antiderivative((~f + ~g * ~x) ^ (~m + 1) / ((~a + ~b * ~x) * (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~m, ~n], ~x) && (EqQ(~n + ~mn, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ(~m, -1) && Not(EqQ(~m, -2) && IntegerQ(~n)))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + 1) * (~g / ~b) ^ ~m * Subst(Antiderivative(((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~p) && (EqQ(~b * ~f - ~a * ~g, 0) && (GtQ(~p, 0) || LtQ(~m, -1)))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + 1) * (~g / ~b) ^ ~m * Subst(Antiderivative(((~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~p) && (EqQ(~b * ~f - ~a * ~g, 0) && (GtQ(~p, 0) || LtQ(~m, -1)))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + 1) * (~g / ~d) ^ ~m * Subst(Antiderivative((~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~p) && (EqQ(~d * ~f - ~c * ~g, 0) && (GtQ(~p, 0) || LtQ(~m, -1)))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) ^ (~m + 1) * (~g / ~d) ^ ~m * Subst(Antiderivative((~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~p) && (EqQ(~d * ~f - ~c * ~g, 0) && (GtQ(~p, 0) || LtQ(~m, -1)))))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && IGtQ(~p, 0)))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => (~b * ~c - ~a * ~d) * Subst(Antiderivative((((~b * ~f - ~a * ~g) - (~d * ~f - ~c * ~g) * ~x) ^ ~m * (~A + ~B * log(~e * (~x) ^ ~n)) ^ ~p) / (~b - ~d * ~x) ^ (~m + 2), ~x), ~x, (~a + ~b * ~x) / (~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~n], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && IGtQ(~p, 0)))))
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * ((~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x)) ^ ~(n'))) ^ ~(p'), ~x) => Unintegrable((~f + ~g * ~x) ^ ~m * (~A + ~B * log(~e * ((~a + ~b * ~x) / (~c + ~d * ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~(a') + ~(b') * ~x) ^ ~(n') * (~(c') + ~(d') * ~x) ^ ~mn)) ^ ~(p'), ~x) => Unintegrable((~f + ~g * ~x) ^ ~m * (~A + ~B * log((~e * (~a + ~b * ~x) ^ ~n) / (~c + ~d * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && IntegerQ(~n))
    @apply_utils Antiderivative((~w) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~u / ~v) ^ ~(n'))) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~w, ~x) ^ ~m * (~A + ~B * log(~e * (ExpandToSum(~u, ~x) / ExpandToSum(~v, ~x)) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~m, ~n, ~p], ~x) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x)))
    @apply_utils Antiderivative((~w) ^ ~(m') * (~(A') + ~(B') * log(~(e') * (~u) ^ ~(n') * (~v) ^ ~mn)) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~w, ~x) ^ ~m * (~A + ~B * log((~e * ExpandToSum(~u, ~x) ^ ~n) / ExpandToSum(~v, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~e, ~A, ~B, ~m, ~n, ~p], ~x) && (EqQ(~n + ~mn, 0) && (IGtQ(~n, 0) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x)))))
end

max_stars_repo_path: src/rules/3 Logarithms/3.2.1 (f+g x)^m (A+B log(e ((a+b x) over (c+d x))^n))^p.jl


text:
function eigenproblem(P, q::AbstractVector{T}, r::T, nev=1;
	tol=0.0, maxiter=300, v0=zeros((0,))) where {T}
	"""
	Calculates rightmost eigenvalues/vectors of

	|-P   qq'/r^2|  |v1|  =  λ|v1|
	| I        -P|  |v2|      |v2|

	with Arpack's eigs.
	"""
	n = length(q)
	if norm(q) < 1e-9
		# Solution is the minimum eigenvector of P
		(λ, V, nconv, niter, nmult, _) = eigs(P, nev=1, which=:SR, tol=tol, maxiter=maxiter)
		return -λ, [0*V; V], niter, nmult
	end

	function A(y::AbstractVector, x::AbstractVector)
		@inbounds y1 = view(y, 1:n); @inbounds y2 = view(y, n+1:2*n)
		@inbounds x1 = view(x, 1:n); @inbounds x2 = view(x, n+1:2*n)
		mul!(y1, P, x1)
		mul!(y2, P, x2)
		axpy!(-dot(q, x2)/r^2, q, y1)
		axpy!(-one(T), x1, y2)
	end
	D = LinearMap{T}(A, 2*n; ismutating=true)

	(λ, V, nconv, niter, nmult, _) = eigs(-D, nev=nev, which=:LR,
	tol=tol, maxiter=maxiter, v0=v0)
	@assert(nconv >= min(nev, 2), "Eigensolver has failed to converge.
		Try decreasing tolerance or changing parameters of eigs.")

	return λ, V, niter, 2*nmult
end

function gen_eigenproblem(P, q::AbstractVector{T}, r::T, C::AbstractMatrix, nev=1;
	tol=0.0, maxiter=500, v0=zeros((0,))) where {T}
	"""
	Calculates rightmost eigenvalues/vectors of

	|-P    qq'/r^2|  |v1|  =  λ |C   0| |v1|
	| C         -P|  |v2|       |0   C| |v2|

	with Arpack's eigs.
	"""
	n = length(q)
	function A(y::AbstractVector, x::AbstractVector)
		@inbounds y1 = view(y, 1:n); @inbounds y2 = view(y, n+1:2*n)
		@inbounds x1 = view(x, 1:n); @inbounds x2 = view(x, n+1:2*n)
		mul!(y1, C, x1)
		mul!(y2, P, x2)
		axpy!(-one(T), y1, y2)

		mul!(y1, P, x1)
		axpy!(-dot(q, x2)/r^2, q, y1)
	end
	D = LinearMap{T}(A, 2*n; ismutating=true)

	(λ, V, nconv, niter, nmult, _) = eigs(-D, [C 0*I; 0*I C], nev=nev, which=:LR, tol=tol)
	@assert(nconv >= min(nev, 2), "Eigensolver has failed to converge.
		Try decreasing tolerance or changing parameters of eigs.")

	return λ, V, niter, 2*nmult
end

function eigenproblem_small(P::AbstractMatrix, q::AbstractVector{T}, r::T, nev=1) where {T}
	"""
	Calculates eigenvalues/vectors of

	|-P    qq'/r^2|  |v1|  =  λ |v1|
	| I         -P|  |v2|       |v2|

	with eigen.
	"""
	n = length(q)
	A = zeros(T, 2*n, 2*n) # The matrix to perform eigendecomposition
	@inbounds for i in 1:n
		for j in 1:n
			A[i, j] = -P[i, j]
		end
		A[i + n, i] = one(T)
		for j in 1:n
			A[i + n, j + n] = -P[i, j]
		end
	end
	@inbounds for i in 1:n
		c = q[i]/r^2
		for j in n+1:2*n
			A[i, j] = c*q[j-n]
		end
	end
	λ, V = eigen!(A)
	return λ, V, 0, 0
end

function gen_eigenproblem_small(P::AbstractMatrix, q::AbstractVector{T}, r::T, C::AbstractMatrix, nev=1) where {T}
	"""
	Calculates rightmost eigenvalues/vectors of

	|-P    qq'/r^2|  |v1|  =  λ |C   0| |v1|
	| C         -P|  |v2|       |0   C| |v2|

	with Arpack's eigs.
	"""
	n = length(q)
	A = zeros(T, 2*n, 2*n) # Left matrix of eigenproblem
	@inbounds for i in 1:n
		for j in 1:n
			A[i, j] = -P[i, j]
		end
		for j in 1:n
			A[i + n, j] = C[i, j]
		end
		for j in 1:n
			A[i + n, j + n] = -P[i, j]
		end
	end
	@inbounds for i in 1:n
		c = q[i]/r^2
		for j in n+1:2*n
			A[i, j] = c*q[j-n]
		end
	end
	B = zeros(T, 2*n, 2*n) # Right matrix of eigenproblem
	@inbounds for i in 1:n
		for j = 1:n
			B[i, j] = C[i, j]
		end
		for j in 1:n
			B[i + n, j + n] = C[i, j]
		end
	end
	λ, V = eigen!(A, B)
	return λ, V, 0, 0
end


text:
module GraphRecipes

using LightGraphs
using NetworkLayout
using PlotUtils         # ColorGradient
using RecipesBase

using InteractiveUtils  # subtypes
using LinearAlgebra
using SparseArrays
using Statistics
using NaNMath
using GeometryTypes
using Interpolations

include("utils.jl")
include("graph_layouts.jl")
include("graphs.jl")
include("misc.jl")

end # module


text:
struct RegressionHyperParametersAuxiliary
    beta_names::Vector{String}
    function RegressionHyperParametersAuxiliary(num_states::Int)
        return new(beta_names(num_states))
    end
end

@doc raw"""
    LinearRegression(X::Matrix{Fl}, y::Vector{Fl}) where Fl

The linear regression state-space model is defined by:
```math
\begin{gather*}
    \begin{aligned}
        y_{t} &=  X_{1,t} \cdot \beta_{1,t} + \dots + X_{n,t} \cdot \beta_{n,t} + \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\
        \beta_{1,t+1} &= \beta_{1,t}\\
        \dots &= \dots\\
        \beta_{n,t+1} &= \beta_{n,t}\\
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = LinearRegression(rand(100, 2), rand(100))
LinearRegression model
```
"""
mutable struct LinearRegression <: StateSpaceModel
    hyperparameters::HyperParameters
    hyperparameters_auxiliary::RegressionHyperParametersAuxiliary
    system::LinearUnivariateTimeVariant
    results::Results
    exogenous::Matrix

    function LinearRegression(X::Matrix{Fl}, y::Vector{Fl}) where Fl

        # Treat possible input errors
        @assert length(y) == size(X, 1)

        num_observations = size(X, 1)
        num_exogenous = size(X, 2)

        Z = [X[t, :] for t in 1:num_observations]
        T = [Matrix{Fl}(I, num_exogenous, num_exogenous) for _ in 1:num_observations]
        R = [zeros(num_exogenous, 1) for _ in 1:num_observations]
        d = [zero(Fl) for _ in 1:num_observations]
        c = [zeros(num_exogenous) for _ in 1:num_observations]
        H = [one(Fl) for _ in 1:num_observations]
        Q = [zeros(Fl, 1, 1) for _ in 1:num_observations]

        system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q)

        names = [["β_$i" for i in 1:num_exogenous]; "sigma2_ε"]
        hyperparameters = HyperParameters{Fl}(names)

        hyperparameters_auxiliary = RegressionHyperParametersAuxiliary(num_states(system))

        return new(hyperparameters, hyperparameters_auxiliary, system, Results{Fl}(), X)
    end
end

function beta_names(num_states::Int)
    str = Vector{String}(undef, num_states)
    for i in eachindex(str)
        str[i] = "β_$i"
    end
    return str
end

function get_beta_name(model::LinearRegression, i::Int)
    return model.hyperparameters_auxiliary.beta_names[i]
end

function fill_H_in_time(model::LinearRegression, H::Fl) where Fl
    return fill_system_matrice_with_value_in_time(model.system.H, H)
end

# Obligatory functions
function default_filter(model::LinearRegression)
    Fl = typeof_model_elements(model)
    a1 = zeros(Fl, num_states(model))
    return RegressionKalmanFilter(a1)
end

function initial_hyperparameters!(model::LinearRegression)
    Fl = typeof_model_elements(model)
    initial_hyperparameters = Dict{String,Fl}("sigma2_ε" => var(model.system.y))
    # The optimal regressors are the result of X \ y
    betas = model.exogenous \ model.system.y
    for i in 1:num_states(model)
        initial_hyperparameters[get_beta_name(model, i)] = betas[i]
    end
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return model
end

function constrain_hyperparameters!(model::LinearRegression)
    for i in 1:num_states(model)
        constrain_identity!(model, get_beta_name(model, i))
    end
    constrain_variance!(model, "sigma2_ε")
    return model
end

function unconstrain_hyperparameters!(model::LinearRegression)
    for i in 1:num_states(model)
        unconstrain_identity!(model, get_beta_name(model, i))
    end
    unconstrain_variance!(model, "sigma2_ε")
    return model
end

function fill_model_system!(model::LinearRegression)
    # Fill the same H for every timestamp
    H = get_constrained_value(model, "sigma2_ε")
    fill_H_in_time(model, H)
    return nothing
end

function fill_model_filter!(filter::KalmanFilter, model::LinearRegression)
    for i in axes(filter.kalman_state.a, 1)
        filter.kalman_state.a[i] = get_constrained_value(model, get_beta_name(model, i))
    end
    return filter
end

function reinstantiate(::LinearRegression, y::Vector{Fl}, X::Matrix{Fl}) where Fl
    return LinearRegression(X, y)
end

has_exogenous(::LinearRegression) = true


text:
module IbnrClaims

export Triangle,
Incremental,
Accumulated,
IncrementalCompleted,
AccumulatedCompleted,
taccum,
tdeacum,
tfactors,
fillcl,
ibnrcl,
tquotas,
expositionbf,
fillbf,
ibnrbf,
randtri,
lossesbf,
incurred,
lossratiosbf

using Distributions

include("trianglestructs.jl")
include("trianglefunctions.jl")
include("prettyprint.jl")
include("chainladder.jl")
include("bornhuetterferguson.jl")

end


text:
using AbstractActuators
using Documenter

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

makedocs(;
    modules=[AbstractActuator],
    authors="Paulo Jabardo <pjabardo@ipt.br>",
    repo="https://github.com/pjsjipt/AbstractActuators.jl/blob/{commit}{path}#{line}",
    sitename="AbstractActuators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
module Distance

using Images

using DeformableRegistration
using DeformableRegistration.Transformation
using DeformableRegistration.Interpolation
using DeformableRegistration.ImageProcessing

export ssdDistance, maskedSsdDistance, ngfDistance, estimateNGFEpsilon
export calculateHessianStaggered

include("distances/ssdDistance.jl")
include("distances/maskedSsdDistance.jl")
include("distances/ngfDistance.jl")
include("distances/estimateNGFEpsilon.jl")

end