text:
# Stochastic Lanczos Quadrature method
# Based on Ubaru, Shashanka, Jie Chen, and Yousef Saad.
# "Fast estimation of tr(F(A)) via stochastic lanczos quadrature, 2016."
# URL: http://www-users.cs.umn.edu/~saad/PDF/ys-2016-04.pdf.

export SLQWorkspace, slq

using LinearAlgebra
using Parameters

struct SLQWorkspace{
    elt, TM<:AbstractMatrix{elt}, FN<:Function, FN2<:Function,
    I<:Integer, TV<:AbstractVector{elt}, AV<:AbstractVector{elt},
    TS<:SymTridiagonal, TM2<:AbstractMatrix{elt}, R<:Real,
}
    A::TM
    fn::FN
    dfn::FN2
    m::I
    v::AV
    nᵥ::I
    ctol::R
    T::TS
    α::TV
    β::TV
    ω::TV
    Y::TM2
    Θ::TV
    v₀::AV
end

# Using v = Vector{elt}(undef, n) will increase the performance but due to
# support issue for CuArray, it must be kept similar(A, n).
# This is a good use-case for full (which was removed in Julia 1.0) where
# we could just get a dense vector from a sparse one solving the performance
# issue
"""
    SLQWorkspace(A::AbstractMatrix; fn::Function, dfn::Function, ctol, m, nv)

Create an SLQWorkspace for supplied SPD Matrix A.
Use it to calculate tr(fn(A)).

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `fn` : Function to apply. By default uses inverse function
- `dfn` : Distribution function for v (random dist. with norm(v) = 1). By default uses rademacherDistribution!(v::Vector, t::Type)
- `ctol` : SLQ Convergence Tolerance value. By default ctol = 0.1
- `m` : Specify value for lanczos steps. By default m = 15
- `nv` : Specify value for SLQ iterations. By default nb = 10
"""
function SLQWorkspace(A; fn::Function=invfun, dfn::Function=rademacherDistribution!, ctol=0.1, m=15, nv=10)
    elt = eltype(A)
    Atype = typeof(A)
    n = size(A, 1)
    #v = Vector{elt}(undef, n)
    v = similar(A, n)
    v₀ = similar(v)
    α = Vector{elt}(undef, m)
    β = Vector{elt}(undef, m-1)
    ω = Vector{elt}(undef, n)
    Y = similar(A, m, m)
    Θ = similar(α)
    T = SymTridiagonal(α, β)
    return SLQWorkspace(A, fn, dfn, m, v, nv, ctol, T, α, β, ω, Y, Θ, v₀)
end

function lcz(w::SLQWorkspace)
    α₀ = zero(eltype(w.A))
    β₀ = zero(eltype(w.A))
    fill!(w.v₀, 0)
    # Following loop executes lanczos steps
    @unpack A, v, ω, v₀, α, β, m, T = w
    for i in 1:w.m
        mul!(ω, A, v)
        α₀ = dot(ω, v)
        ω .= ω .- (α₀ .* v) .- (β₀ .* v₀)
        β₀ = norm(ω)

        α[i] = α₀
        if i < m
            β[i] = β₀
        end

        # Sparse Vectors do not support copy!
        #copy!(v₀, v)
        v₀ .= v
        v .= ω ./ β₀
    end
end

"""
    slq(w::SLQWorkspace; skipverify = false)
    slq(A::AbstractMatrix; skipverify = false, fn::Function = invfun,
    dfn::Function = Base.rand, ctol = 0.1, eps = ϵ, mtol = tol)

SLQ method to calculate tr(fn(A)) for a Symmetric Positive Definite matrix
A and an analytic function fn.

# Arguments
- `A` : Symmetric Positive Definite Matrix
- `skipverify` : Skip isposdef(A) verification. By default skipverify = false
- `fn` : Function to apply on A before trace calculation. fn must be analytic λₘ and λ₁ of A. By default fn = inv
- `dfn` : Distribution function for v (random dist. with norm(v) = 1). By default uses rademacherDistribution!(v::Vector, t::Type)
- `ctol` : SLQ Convergence Tolerance value. Decrease this for higher precision. By default ctol = 0.1
- `eps` : Error bound for lanczos steps calculation. Decrease this for higher accuracy. By default eps = 0.05
- `mtol` : Tolerance for eigenvalue Convergence. Decrease this for precision. By default mtol = 0.01
"""
function slq(w::SLQWorkspace; skipverify = false)
    @unpack A, v, T, Y, Θ, m, nᵥ, ctol, fn, dfn = w
    tr = zero(eltype(w.A))
    if skipverify || isposdef(w.A)
        actual_nᵥ = nᵥ
        for i in 1:w.nᵥ
            prev_tr = tr
            # Create a uniform random distribution with norm(v) = 1
            dfn(v)
            v .= v ./ norm(v)
            # Run lanczos algorithm to find estimate Ritz SymTridiagonal
            lcz(w)
            Y .= eigvecs(T)
            Θ .= eigvals(T)
            for j in 1:m
                τ = Y[1,j]
                tr = tr + τ^2 * fn(Θ[j])
            end
            if isapprox(prev_tr, tr, rtol = ctol)
                actual_nᵥ = i
                break
            end
        end
        tr = size(w.A, 1)/actual_nᵥ * tr
    else
        throw("Given Matrix is NOT Symmetric Positive Definite")
    end

    return tr
end

function slq(A::AbstractMatrix; skipverify = false, fn::Function = invfun, dfn::Function = rademacherDistribution!, ctol = 0.1, eps = ϵ, mtol = ξ)

    # Estimate eigmax and eigmin for SLQ bounds
    mval = Int64(ceil(log(eps/(1.648 * sqrt(size(A, 1))))/(-2 * sqrt(mtol))))
    w = SLQWorkspace(A, fn = fn, dfn = dfn, m = mval)
    #rademacherDistribution!(w.v)
    w.dfn(w.v)
    w.v .= w.v ./ norm(w.v)
    lcz(w)
    λₘ = eigmax(w.T)
    λ₁ = eigmin(w.T)

    if λ₁ < 0 && λₘ > 0
        @warn "Eigenvalues cross zero. Functions like log may not give correct results. Try scaling the input matrix."
    end

    # SLQ bounds
    # Todo: Research and create better bounds for λ₁ < 1 && λₘ > 1 case
    κ = λₘ/λ₁
    Mₚ = fn(λₘ)
    mₚ = fn(λ₁)
    ρ = (sqrt(κ) + 1)/(sqrt(κ) - 1)
    K = ((λₘ - λ₁) * (sqrt(κ) - 1)^2 * Mₚ)/(sqrt(κ) * mₚ)
    mval = Int64(ceil((sqrt(κ)/4) * log(K/eps)))
    if mval < 10
        @warn "Low lanczos step value. Try decreasing eps and mtol for better accuracy."
        mval = 5
    end
    nval = Int64(ceil((24/ϵ^2) * log(2/mtol)))

    # Re-construct SLQWorkspace
    w = SLQWorkspace(A, fn = fn, dfn = dfn, m = mval, nv = nval, ctol = ctol)
    slq(w, skipverify = skipverify)
end


text:
"""
    threshold(w::AbstractVector, q::Number, method=:knockoff)

Chooses a threshold τ > 0 by setting
τ = min{ t > 0 : {#j: w[j] ≤ -t} / {#j: w[j] ≥ t} ≤ q }.

# Inputs
+ `w`: Vector of feature important statistics
+ `q`: target FDR (between 0 and 1)
+ `method`: either `:knockoff` or `:knockoff_plus`

# Reference:
Equation 3.10 (`method=:knockoff`) or 3.11 (`method=:knockoff_plus`) of
"Panning for Gold: Model-X Knockoffs for High-dimensional
Controlled Variable Selection" by Candes, Fan, Janson, and Lv (2018)
"""
function threshold(w::AbstractVector{T}, q::Number, method=:knockoff) where T <: AbstractFloat
    0 ≤ q ≤ 1 || error("Target FDR should be between 0 and 1 but got $q")
    offset = method == :knockoff ? 0 : method == :knockoff_plus ? 1 :
        error("method should be :knockoff or :knockoff_plus but was $method.")
    τ = typemax(T)
    for t in sort!(abs.(w), rev=true) # t starts from largest |W|
        ratio = (offset + count(x -> x ≤ -t, w)) / count(x -> x ≥ t, w)
        ratio ≤ q && (τ = t)
    end
    return τ
end

ext: jl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	energy_share_requirement(EP::Model, inputs::Dict, setup::Dict)
This function establishes constraints that can be flexibily applied to define alternative forms of policies that require generation of a minimum quantity of megawatt-hours from a set of qualifying resources, such as renewable portfolio standard (RPS) or clean electricity standard (CES) policies prevalent in different jurisdictions. These policies usually require that the annual MWh generation from a subset of qualifying generators has to be higher than a pre-specified percentage of load from qualifying zones.
The implementation allows for user to define one or multiple RPS/CES style minimum energy share constraints, where each constraint can cover different combination of model zones to mimic real-world policy implementation (e.g. multiple state policies, multiple RPS tiers or overlapping RPS and CES policies). The number of energy share requirement constraints is specified by the user by the value of the GenX settings parameter ```EnergyShareRequirement``` (this value should be an integer >=0). For each constraint $p \in \mathcal{P}^{ESR}$, we define a subset of zones $z \in \mathcal{Z}^{ESR}_{p} \subset \mathcal{Z}$ that are eligible for trading renewable/clean energy credits to meet the corresponding renewable/clean energy requirement.
For each energy share requirement constraint $p \in \mathcal{P}^{ESR}$, we specify the share of total demand in each eligible model zone, $z \in \mathcal{Z}^{ESR}_{p}$, that must be served by qualifying resources, $\mathcal{G}_{p}^{ESR} \subset \mathcal{G}$:
```math
\begin{aligned}
&\sum_{z \in \mathcal{Z}_{p}^{ESR}} \sum_{y \in \mathcal{G}_{p}^{ESR}} \sum_{t \in \mathcal{T}} (\omega_{t} \times  \Theta_{y,z,t}) \geq  \sum_{z \in \mathcal{Z}^{ESR}_{p}} \sum_{t \in \mathcal{T}} (\mu_{p,z}^{ESR} \times \omega_{t} \times D_{z,t}) + \\
&\sum_{y \in \mathcal{O}} \sum_{z \in \mathcal{Z}^{ESR}_{p}} \sum_{t \in \mathcal{T}} \left(\mu_{p,z}^{ESR} \times \omega_{t} \times (\Pi_{y,z,t} - \Theta_{y,z,t}) \right) \hspace{1 cm}  \forall p \in \mathcal{P}^{ESR} \\
\end{aligned}
```
The final term in the summation above adds roundtrip storage losses to the total load to which the energy share obligation applies. This term is included in the constraint if the GenX setup parameter ```StorageLosses=1```. If ```StorageLosses=0```, this term is removed from the constraint. In practice, most existing renewable portfolio standard policies do not account for storage losses when determining energy share requirements. However, with 100% RPS or CES policies enacted in several jurisdictions, policy makers may wish to include storage losses in the minimum energy share, as otherwise there will be a difference between total generation and total load that will permit continued use of non-qualifying resources (e.g. emitting generators).
"""
function energy_share_requirement(EP::Model, inputs::Dict, setup::Dict)

	println("Energy Share Requirement Policies Module")

	## Energy Share Requirements (minimum energy share from qualifying renewable resources) constraint
	@constraint(EP, cESRShare[ESR=1:inputs["nESR"]], EP[:eESR][ESR] >= 0)

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
function demoCars()

    wife,husband,income=1,2,3 # give variables a number

    pi=PotArray(income,[0.8 0.2]) # income distribution
    phgi=PotArray([husband income],[0.2 0; 0.8 0; 0 0.3; 0 0.7]) # p(husband_car|income)
    pwgi=PotArray([wife,income],[0.7 0.2; 0.3 0.1; 0 0.4; 0 0.3]) # p(wife_car|income)
    joint=pwgi*phgi*pi # p(wife_car,husband_car,income)=p(wife_car|income)p(husband_car|income)p(income)
    pwh=sum(joint,income)

    #
    println("Check if w is independent of h.  mean(|p(w,h)-p(w)p(h)|) :")
    println(mean(abs.((pwh-sum(pwh,husband)*sum(pwh,wife)).content)))
    println("This is not zero, so w and h are dependent")

end


text:
using Revise
include("execute.jl")
using .execute
import LPDM, D3Trees, LightDarkLPDM

# execute.batch_execute(dims=2)
# Juno.@run execute.batch_execute(dims=1, debugger=true)
execute.batch_execute(dims=2, debugger=false)
# execute.batch_execute(dims=1, debugger=true)

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
abstract type AbstractDiffusion end
abstract type AbstractStaticDiffusion <: AbstractDiffusion end
abstract type AbstractDynamicDiffusion <: AbstractDiffusion end
isstatic(diffusion::AbstractStaticDiffusion) = true
isdynamic(diffusion::AbstractStaticDiffusion) = false
isstatic(diffusion::AbstractDynamicDiffusion) = false
isdynamic(diffusion::AbstractDynamicDiffusion) = true
initial_diffusion(diffusion::AbstractDiffusion, d, q, Eltype) = one(Eltype)

struct FixedDiffusion <: AbstractStaticDiffusion end
function estimate_diffusion(rule::FixedDiffusion, integ)
    @unpack d, measurement, m_tmp = integ.cache
    sol_diffusions = integ.sol.diffusions

    v, S = measurement.μ, measurement.Σ
    e, _ = m_tmp.μ, m_tmp.Σ.mat
    _S = copy!(m_tmp.Σ.mat, S.mat)

    if iszero(v)
        return zero(integ.cache.global_diffusion)
    end
    if iszero(S)
        return Inf
    end

    # diffusion_t = v' * inv(S) * v / d
    if _S isa Diagonal
        e .= v ./ _S.diag
    else
        S_chol = cholesky!(_S)
        ldiv!(e, S_chol, v)
    end
    diffusion_t = dot(v, e) / d

    if integ.success_iter == 0
        @assert length(sol_diffusions) == 0
        return diffusion_t, diffusion_t
    else
        @assert length(sol_diffusions) == integ.success_iter
        diffusion_prev = sol_diffusions[end]
        global_diffusion =
            diffusion_prev + (diffusion_t - diffusion_prev) / integ.success_iter
        return diffusion_t, global_diffusion
    end
end

"""Maximum a-posteriori Diffusion estimate when using an InverseGamma(1/2,1/2) prior

The mode of an InverseGamma(α,β) distribution is given by β/(α+1)
To compute this in an on-line basis from the previous Diffusion, we reverse the computation to
get the previous sum of residuals from Diffusion, and then modify that sum and compute the new
Diffusion.
"""
struct MAPFixedDiffusion <: AbstractStaticDiffusion end
function estimate_diffusion(rule::MAPFixedDiffusion, integ)
    @unpack d, measurement = integ.cache
    diffusions = integ.sol.diffusions

    N = integ.success_iter + 1
    v, S = measurement.μ, measurement.Σ
    res_t = v' * inv(S) * v / d

    α, β = 1 / 2, 1 / 2
    if integ.success_iter == 0
        @assert length(diffusions) == 0
        diffusion_t = (β + 1 / 2 * res_t) / (α + N * d / 2 + 1)
        return res_t, diffusion_t
    else
        @assert length(diffusions) == integ.success_iter
        diffusion_prev = diffusions[end]
        res_prev = (diffusion_prev * (α + (N - 1) * d / 2 + 1) - β) * 2
        res_sum_t = res_prev + res_t
        diffusion = (β + 1 / 2 * res_sum_t) / (α + N * d / 2 + 1)
        return res_t, diffusion
    end
end

struct DynamicDiffusion <: AbstractDynamicDiffusion end
function estimate_diffusion(kind::DynamicDiffusion, integ)
    @unpack d, R, H, Qh, measurement, m_tmp = integ.cache
    # @assert all(R .== 0) "The dynamic-diffusion assumes R==0!"
    z = measurement.μ
    e, HQH = m_tmp.μ, m_tmp.Σ
    X_A_Xt!(HQH, Qh, H)
    if HQH.mat isa Diagonal
        e .= z ./ HQH.mat.diag
        σ² = dot(e, z) / d
        return σ², σ²
    else
        C = cholesky!(HQH.mat)
        ldiv!(e, C, z)
        σ² = dot(z, e) / d
        return σ², σ²
    end
end

struct MVDynamicDiffusion <: AbstractDynamicDiffusion end
initial_diffusion(diffusion::MVDynamicDiffusion, d, q, Eltype) =
    kron(Diagonal(ones(Eltype, d)), Diagonal(ones(Eltype, q + 1)))
function estimate_diffusion(kind::MVDynamicDiffusion, integ)
    @unpack d, q, R, P, PI, E1 = integ.cache
    @unpack H, Qh, measurement, m_tmp = integ.cache
    z = measurement.μ

    # @assert all(R .== 0) "The dynamic-diffusion assumes R==0!"

    # Assert EK0
    @assert all(H .== E1)

    # More safety checks
    HQH = X_A_Xt!(m_tmp.Σ, Qh, H)
    @assert isdiag(HQH)
    @assert length(unique(diag(HQH))) == 1
    Q0_11 = diag(HQH)[1]

    Σ_ii = z .^ 2 ./ Q0_11
    Σ_ii .= max.(Σ_ii, eps(eltype(Σ_ii)))
    Σ = Diagonal(Σ_ii)

    Σ_out = kron(Σ, I(q + 1))
    @assert isdiag(Σ_out)
    # @info "MVDynamic diffusion" Σ Σ_out
    return Σ_out, Σ_out
end

struct MVFixedDiffusion <: AbstractStaticDiffusion end
initial_diffusion(diffusion::MVFixedDiffusion, d, q, Eltype) =
    kron(Diagonal(ones(Eltype, d)), Diagonal(ones(Eltype, q + 1)))
function estimate_diffusion(kind::MVFixedDiffusion, integ)
    @unpack d, q, R, P, PI, E1 = integ.cache
    @unpack measurement, H = integ.cache
    @unpack d, measurement = integ.cache
    diffusions = integ.sol.diffusions

    # Assert EK0
    @assert all(H .== E1)

    @unpack measurement = integ.cache
    v, S = measurement.μ, measurement.Σ

    # More safety checks
    @assert isdiag(S)
    # @assert length(unique(diag(S))) == 1
    S_11 = diag(S)[1]

    Σ_ii = v .^ 2 ./ S_11
    Σ = Diagonal(Σ_ii)
    Σ_out = kron(Σ, I(q + 1))
    @assert isdiag(Σ_out)
    # @info "MV-MLE-Diffusion" v S Σ Σ_out

    if integ.success_iter == 0
        @assert length(diffusions) == 0
        return Σ_out, Σ_out
    else
        @assert length(diffusions) == integ.success_iter
        diffusion_prev = diffusions[end]
        diffusion = diffusion_prev + (Σ_out - diffusion_prev) / integ.success_iter
        return Σ_out, diffusion
    end
end


text:
@testset "Utils" begin
    J = 8
    n = 2^J
    x = testfunction(n, "Bumps")
    y = dwt(x, wavelet(WT.cdf97, WT.Lifting))
    d,l = wplotdots(y, 0.1, n)
    A = wplotim(y)

    x = randn(64,64)
    L = 2
    xts = wplotim(x, L, wavelet(WT.db3, WT.Filter))
end


text:
using Bio.Seq
using BioBenchmarks

chrom = first(open(datafile("chr1.fa"), FASTA)).seq
@benchmark reverse_complement(reverse_complement(chrom))


text:
"Convert a gridded drainage direction to a directed graph"
function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcr_dir::AbstractVector)
    # prepare a directed graph to be filled
    n = length(inds)
    graph = DiGraph(n)

    # loop over ldd, adding the edge to the downstream node
    for (from_node, from_index) in enumerate(inds)
        ldd_val = ldd[from_node]
        # skip pits to prevent cycles
        ldd_val == 5 && continue
        to_index = from_index + pcr_dir[ldd_val]
        # find the node id of the downstream cell
        to_node = searchsortedfirst(inds, to_index)
        add_edge!(graph, from_node, to_node)
    end
    if is_cyclic(graph)
        error("""One or more cycles detected in flow graph.
            The provided local drainage direction map may be unsound.
            Verify that each active flow cell flows towards a pit.
            """)
    end
    return graph
end

"Kinematic wave flow rate for a single cell and timestep"
function kinematic_wave(Qin, Qold, q, α, β, Δt, Δx)
    ϵ = 1.0e-12
    max_iters = 3000

    if Qin + Qold + q ≈ 0.0
        return 0.0
    else
        # common terms
        ab_pQ = α * β * pow(((Qold + Qin) / 2.0), (β - 1.0))
        Δtx = Δt / Δx
        C = Δtx * Qin + α * pow(Qold, β) + Δt * q

        Qkx = (Δtx * Qin + Qold * ab_pQ + Δt * q) / (Δtx + ab_pQ)
        if isnan(Qkx)
            Qkx = 0.0
        end
        Qkx = max(Qkx, 1.0e-30)
        count = 1

        while true
            fQkx = Δtx * Qkx + α * pow(Qkx, β) - C
            dfQkx = Δtx + α * β * pow(Qkx, (β - 1.0))
            Qkx = Qkx - fQkx / dfQkx
            Qkx = max(Qkx, 1.0e-30)
            if (abs(fQkx) <= ϵ) || (count >= max_iters)
                break
            end
            count += 1
        end

        return Qkx
    end
end

"Kinematic wave flow rate over the whole network for a single timestep"
function kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt)
    for v in toposort
        upstream_nodes = inneighbors(graph, v)
        Qin = sum_at(Q, upstream_nodes)
        Q[v] = kinematic_wave(Qin, Qold[v], q, α[v], β, Δt, DCL[v])
    end
    return Q
end

"Kinematic wave for lateral subsurface flow for a single cell and timestep"
function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, dw, ssfmax)

    ϵ = 1.0e-12
    max_iters = 3000

    if ssfin + ssfₜ₋₁ ≈ 0.0 && r <= 0.0
        return 0.0, d, 0.0
    else
        # initial estimate
        ssf = (ssfₜ₋₁ + ssfin) / 2.0
        count = 1

        # Estimate zi on the basis of the relation between subsurfacel flow and zi
        zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f
        # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation)
        Cn = (kh₀ * exp(-f * zi) * β) / θₑ
        # Term of the continuity equation for Newton-Raphson iteration for iteration 1
        # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
        # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
        c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssf + Δt * (r - (ziₜ₋₁ - zi) * θₑ * dw)

        # Continuity equation of which solution should be zero
        fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
        # Derivative of the continuity equation w.r.t. Q_out for iteration 1
        dfQ = (Δt / Δx) + 1.0 / Cn
        # Update lateral outflow estimate ssf (Q_out) for iteration 1
        ssf = ssf - (fQ / dfQ)
        if isnan(ssf)
            ssf = 0.0
        end
        ssf = max(ssf, 1.0e-30)

        # Start while loop of Newton-Raphson iteration m until continuity equation approaches zero
        while true
            # Estimate zi on the basis of the relation between lateral flow rate and groundwater level
            zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f
            # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation
            Cn = (kh₀ * exp(-f * zi) * β) / θₑ

            # Term of the continuity equation for given Newton-Raphson iteration m
            # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
            # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
            c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssf + Δt * (r - (ziₜ₋₁ - zi) * θₑ * dw)

            # Continuity equation of which solution should be zero
            fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
            # Derivative of the continuity equation w.r.t. Q_out for iteration m+1
            dfQ = (Δt / Δx) + 1.0 / Cn
            # Update lateral outflow estimate ssf (Q_out) for iteration m+1
            ssf = ssf - (fQ / dfQ)
            if isnan(ssf)
                ssf = 0.0
            end
            ssf = max(ssf, 1.0e-30)
            if (abs(fQ) <= ϵ) || (count >= max_iters)
                break
            end
            count += 1
        end

        # Constrain the lateral flow rate ssf
        ssf = min(ssf, (ssfmax * dw))
        # On the basis of the lateral flow rate, estimate the amount of groundwater level above surface (saturation excess conditions), then rest = negative
        rest = ziₜ₋₁ - (ssfin + r * Δx - ssf) / (dw * Δx) / θₑ
        # In case the groundwater level lies above surface (saturation excess conditions, rest = negative), calculate the exfiltration rate and set groundwater back to zero.
        exfilt = min(rest, 0.0) * -θₑ
        zi = clamp(zi, 0.0, d)

        return ssf, zi, exfilt

    end
end

"""
    accucapacitystate!(material, network, capacity) -> material

Transport of material downstream with a limited transport capacity over a directed graph.
Mutates the material input. The network is expected to hold a graph and order field, where
the graph implements the Graphs interface, and the order is a valid topological ordering
such as that returned by `Graphs.topological_sort_by_dfs`.

Returns the material state after transport.
"""
function accucapacitystate!(material, network, capacity)
    @unpack graph, order = network
    for v in order
        downstream_nodes = outneighbors(graph, v)
        n = length(downstream_nodes)
        flux_val = min(material[v], capacity[v])
        material[v] -= flux_val
        if n == 0
            # pit: material is transported out of the map if a capacity is set,
            # cannot add the material anywhere
        elseif n == 1
            material[only(downstream_nodes)] += flux_val
        else
            error("bifurcations not supported")
        end
    end
    return material
end

"""
    accucapacitystate!(material, network, capacity) -> material

Non mutating version of `accucapacitystate!`.
"""
function accucapacitystate(material, network, capacity)
    accucapacitystate!(copy(material), network, capacity)
end

"""
    accucapacityflux!(flux, material, network, capacity) -> flux, material

Transport of material downstream with a limited transport capacity over a directed graph.
Updates the material input, and overwrites the flux input, not using existing values. The
network is expected to hold a graph and order field, where the graph implements the Graphs
interface, and the order is a valid topological ordering such as that returned by
`Graphs.topological_sort_by_dfs`.

Returns the flux (material leaving each cell), and material (left after transport).
"""
function accucapacityflux!(flux, material, network, capacity)
    @unpack graph, order = network
    for v in order
        downstream_nodes = outneighbors(graph, v)
        n = length(downstream_nodes)
        flux_val = min(material[v], capacity[v])
        material[v] -= flux_val
        flux[v] = flux_val
        if n == 0
            # pit: material is transported out of the map if a capacity is set,
            # cannot add the material anywhere
        elseif n == 1
            material[only(downstream_nodes)] += flux_val
        else
            error("bifurcations not supported")
        end
    end
    return flux, material
end

"""
    accucapacityflux!(material, network, capacity) -> flux, material

Non mutating version of `accucapacityflux!`.
"""
function accucapacityflux(material, network, capacity)
    accucapacityflux!(zero(material), copy(material), network, capacity)
end

"""
    lateral_snow_transport!(snow, snowwater, slope, network)

Lateral snow transport. Transports snow downhill. Mutates `snow` and `snowwater`.
"""
function lateral_snow_transport!(snow, snowwater, slope, network)
    snowflux_frac = min.(0.5, slope ./ 5.67) .* min.(1.0, snow ./ 10000.0)
    maxflux = snowflux_frac .* snow
    snow = accucapacitystate!(snow, network, maxflux)
    snowwater = accucapacitystate!(snowwater, network, snowwater .* snowflux_frac)
    return snow, snowwater
end


text:
    function catch_phrase()
        resultado =""
        for word_list in data["en"]["faker"]["company"]["buzzwords"]
            resultado *= executor(word_list)*" "
        end
        resultado
    end

    function bs()
        resultado =""
        for word_list in data["en"]["faker"]["company"]["bs"]
            resultado *= executor(word_list)*" "
        end
        resultado
    end

    company_suffix()= executor(data["en"]["faker"]["company"]["suffix"])

    company_prefix() = executor(data["en"]["faker"]["company"]["prefix"])
    company()=executor(data["en"]["faker"]["company"]["name"])


text:
module OptimizationLH

using Distributions, Formatting, Parameters, Random
using NLopt

export Model, calibrate, make_endowments, decision_probabilities
export make_random_parameters

include("discretize.jl")
include("types.jl")
include("model.jl")

"""
    calibrate()

Calibrate the model
"""
function calibrate(mdl :: Model, algo :: Symbol)
    # algo = :LN_COBYLA;
    # algo = :GN_CRS2_LM;
    # algo = :LN_SBPLX;

    # Solve with random parameters to make attainable calibration targets
    mpTrue = make_random_parameters(mdl, 392);
    solnV = mp_to_guess(mdl, mpTrue);

    # Calibration targets = model solution
    tgStats = solve(mdl, mpTrue);

    # Objective should recover tgStats
    dev = objective(mdl, tgStats, mpTrue);
    @assert isa(dev, Float64);
    @assert dev ≈ 0.0

    # Objective wrapper to be called by solver
    function dev_fct(guessInV :: Vector{Float64}, gradV)
        @assert length(gradV) == 0
        mpDev = guess_to_mp(mdl, guessInV);
        return objective(mdl, tgStats, mpDev);
    end

    # Check that wrapper returns the same as calling objective directly
    dev2 = dev_fct(solnV, []);
    @assert dev2 ≈ dev

    ## ------  Optimization

    nParams = n_params(mdl);
    println("\n--------------------")
    println("Starting $algo with $nParams parameters");
    nTgMoments = n_tg_moments(tgStats);
    println("  No of target moments: $nTgMoments")

    mpGuess = make_random_parameters(mdl, 493);
    guessV = mp_to_guess(mdl, mpGuess);
    dev0 = dev_fct(guessV, []);
    println("Initial deviation: $dev0")

    optS = Opt(algo,  nParams);
    optS.min_objective = dev_fct;
    optS.lower_bounds = 1.0;
    optS.upper_bounds = 2.0;
    optS.stopval = 0.001;
    optS.stopval = 0.001;
    optS.ftol_rel = 0.001;
    optS.maxeval = 4_000;
    mdl.currentIter = 0;
    optf, optx, ret = optimize(optS, guessV);

    optfr = round(optf, digits = 3);
    nIter = mdl.currentIter;
    println("Optimization complete.   $nIter iterations.   f = $optfr");
    println("  Return value: $ret");
    mp = guess_to_mp(mdl, optx);
    devSol = objective(mdl, tgStats, mp);
    @assert abs(devSol - optf) < 0.001

    meanAbsDev = round(mean(abs.(optx .- solnV)), digits = 3);
    println("  Mean abs deviation vs true solution x: $meanAbsDev");
    println("----- Initial parameters:")
    show_params(mpTrue);
    println("----- Solution parameters:")
    show_params(mp);

    return nothing
end

end # module


text:
function montecarlo(
    f::Function,
    a::Real,
    b::Real,
    iterations::Integer = 10_000,
)
    width = (b - a)
    estimate = 0.0
    for i in 1:iterations
        x = width * rand() + a
        estimate += f(x) * width
    end
    return estimate / iterations
end


text:
import Base: +, -, /, *, <=, >=

# Random Normal Variable
(+)(lhs::FlexibilityVariable,rhs::FlexibilityVariable) = FlexibilityExpr([lhs,rhs],[1.0,1.0],0.0)
(-)(lhs::FlexibilityVariable,rhs::FlexibilityVariable) = FlexibilityExpr([lhs,rhs],[1.0,-1.0],0.0)
(/)(lhs::FlexibilityVariable,rhs::Number) = FlexibilityExpr([lhs],[1/rhs],0.0)
(-)(rhs::FlexibilityVariable) = FlexibilityExpr([rhs],[-1.0],0.0)

for op in (:+, :-, :*)
    @eval begin
        ($op)(lhs::Number, rhs::FlexibilityVariable) = ($op)(lhs, FlexibilityExpr([rhs],[1.0],0.0))
        ($op)(lhs::Variable, rhs::FlexibilityVariable) = ($op)(convert(AffExpr,lhs),rhs)
        ($op)(lhs::Variable, rhs::FlexibilityExpr) = ($op)(convert(AffExpr,lhs),rhs)
    end
    if op == :-
        @eval ($op)(lhs::FlexibilityExpr, rhs::Variable) = (+)(lhs,-rhs)
        @eval ($op)(lhs::FlexibilityVariable, rhs::Union{Variable,Number}) = (+)(lhs,-rhs)
    else
        @eval ($op)(lhs::FlexibilityVariable, rhs::Union{Variable,Number}) = ($op)(rhs,lhs)
        @eval ($op)(lhs::FlexibilityExpr, rhs::Variable) = ($op)(rhs,lhs)
    end
end

# AffExpr
# AffExpr--FlexibilityVariable
(+)(lhs::AffExpr, rhs::FlexibilityVariable) = FlexibilityExpr([rhs],[convert(AffExpr,1.0)],lhs)
(-)(lhs::AffExpr, rhs::FlexibilityVariable) = FlexibilityExpr([rhs],[convert(AffExpr,-1.0)],lhs)
(*)(lhs::AffExpr, rhs::FlexibilityVariable) = FlexibilityExpr([rhs],[lhs],AffExpr())
(+)(lhs::FlexibilityVariable, rhs::AffExpr) = rhs+lhs
(-)(lhs::FlexibilityVariable, rhs::AffExpr) = (+)(lhs,-rhs)
(*)(lhs::FlexibilityVariable, rhs::AffExpr) = rhs*lhs

# AffExpr--FlexibilityExpr
(+)(lhs::AffExpr,rhs::FlexibilityExpr) = convert(FlexibilityExpr,lhs)+rhs
(+)(lhs::FlexibilityExpr,rhs::AffExpr) = rhs+lhs
(-)(lhs::AffExpr,rhs::FlexibilityExpr) = convert(FlexibilityExpr,lhs)-rhs
(-)(lhs::FlexibilityExpr,rhs::AffExpr) = lhs-convert(FlexibilityExpr,rhs)
Base.promote_rule(::Type{AffExpr},::Type{FlexibilityExpr}) = FlexibilityExpr
Base.promote_rule(::Type{Float64},::Type{FlexibilityExpr}) = FlexibilityExpr
Base.convert(::Type{FlexibilityExpr},a::AffExpr) = FlexibilityExpr(FlexibilityVariable[],AffExpr[],a)

# AffExpr--FlexibilityExpr
(*)(lhs::AffExpr,rhs::FlexibilityExpr) = FlexibilityExpr(rhs.vars, [lhs*c for c in rhs.coeffs], lhs*rhs.constant)
(*)(lhs::FlexibilityExpr,rhs::AffExpr) = rhs*lhs


text:
# High-level plotting functions for use with `BlockTask`s

"""
    showsample([backend], task, sample)

Show an unprocessed `sample` for `LearningTask` `task` to
`backend::`[`ShowBackend`](#).

## Examples

```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
task = ImageClassificationSingle(data)
sample = getobs(data, 1)
showsample(task, sample)  # select backend automatically
showsample(ShowText(), task, sample)
```
"""
function showsample(backend::ShowBackend, task::AbstractBlockTask, sample)
    blocks = ("Input" => getblocks(task)[1], "Target" => getblocks(task)[2])
    showblock(backend, blocks, sample)
end
showsample(task::AbstractBlockTask, sample) =
    showsample(default_showbackend(), task, sample)

"""
    showsample([backend], task, sample)

Show a vector of unprocessed `samples` for `LearningTask` `task` to
`backend::`[`ShowBackend`](#).

## Examples

```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
task = ImageClassificationSingle(data)
samples = [getobs(data, i) for i in 1:4]
showsamples(task, samples)  # select backend automatically
showsamples(ShowText(), task, samples)
```
"""
function showsamples(backend::ShowBackend, task::AbstractBlockTask, samples)
    showblocks(backend, "Sample" => getblocks(task).sample, samples)
end
showsamples(task::AbstractBlockTask, samples) =
    showsamples(default_showbackend(), task, sample)

"""
    showencodedsample([backend], task, encsample)

Show an encoded sample `encsample` to `backend`.
"""
function showencodedsample(backend::ShowBackend, task::AbstractBlockTask, encsample)
    showblockinterpretable(
        backend,
        getencodings(task),
        getblocks(task).encodedsample,
        encsample,
    )
end
showencodedsample(task, encsample) =
    showencodedsample(default_showbackend(), task, encsample)

"""
    showencodedsamples([backend], task, encsamples)

Show a vector of encoded samples `encsamples` to `backend`.
"""
function showencodedsamples(
    backend::ShowBackend,
    task::AbstractBlockTask,
    encsamples::AbstractVector,
)
    xblock, yblock = encodedblockfilled(getencodings(task), getblocks(task).encodedsample)
    showblocksinterpretable(
        backend,
        getencodings(task),
        ("x" => xblock, "y" => yblock),
        encsamples,
    )
end

"""
    showbatch([backend], task, batch)

Show a collated batch of encoded samples to `backend`.
"""
function showbatch(backend::ShowBackend, task::AbstractBlockTask, batch)
    encsamples = collect(DataLoaders.obsslices(batch))
    showencodedsamples(backend, task, encsamples)
end
showbatch(task, batch) = showbatch(default_showbackend(), task, batch)

"""
    showprediction([backend], task, pred)
    showprediction([backend], task, sample, pred)

Show a prediction `pred`. If a `sample` is also given, show it next to
the prediction. ŷ
"""
function showprediction(backend::ShowBackend, task::AbstractBlockTask, pred)
    showblock(backend, "Prediction" => getblocks(task).pred, pred)
end

function showprediction(backend::ShowBackend, task::AbstractBlockTask, sample, pred)
    blocks = getblocks(task)
    showblock(
        backend,
        ("Sample" => blocks.sample, "Prediction" => blocks.pred),
        (sample, pred),
    )
end

showprediction(task::AbstractBlockTask, args...) =
    showprediction(default_showbackend(), task, args...)

"""
    showpredictions([backend], task, preds)
    showpredictions([backend], task, samples, preds)

Show predictions `pred`. If `samples` are also given, show them next to
the prediction.
"""
function showpredictions(backend::ShowBackend, task::AbstractBlockTask, preds)
    predblock = decodedblockfilled(getencodings(task), getblocks(task).ŷ)
    showblocks(backend, "Prediction" => predblock, preds)
end

function showpredictions(backend::ShowBackend, task::AbstractBlockTask, samples, preds)
    predblock = decodedblockfilled(getencodings(task), getblocks(task).ŷ)
    showblocks(
        backend,
        ("Sample" => getblocks(task), "Prediction" => predblock),
        collect(zip(samples, preds)),
    )
end

showpredictions(task::AbstractBlockTask, args...) =
    showpredictions(default_showbackend(), task, args...)

"""
    showoutput([backend], task, output)
    showoutput([backend], task, encsample, output)

Show a model output to `backend`. If an encoded sample `encsample` is also
given, show it next to the output.
"""
function showoutput(backend::ShowBackend, task::AbstractBlockTask, output)
    showblockinterpretable(
        backend,
        getencodings(task),
        "Output" => getblocks(task).ŷ,
        output,
    )
end
function showoutput(backend::ShowBackend, task::AbstractBlockTask, encsample, output)
    blocks = getblocks(task)
    showblockinterpretable(
        backend,
        getencodings(task),
        ("Encoded sample" => blocks.encodedsample, "Output" => blocks.ŷ),
        (encsample, output),
    )
end
showoutput(task::AbstractBlockTask, args...) =
    showoutput(default_showbackend(), task, args...)

"""
    showoutputs([backend], task, outputs)
    showoutputs([backend], task, encsamples, outputs)

Show model outputs to `backend`. If a vector of encoded samples `encsamples` is also
given, show them next to the outputs. Use [`showoutputbatch`](#) to show collated
batches of outputs.
"""
function showoutputs(backend::ShowBackend, task::AbstractBlockTask, outputs)
    showblocksinterpretable(
        backend,
        getencodings(task),
        "Output" => getblocks(task).ŷ,
        outputs,
    )
end
function showoutputs(backend::ShowBackend, task::AbstractBlockTask, encsamples, outputs)
    blocks = getblocks(task)
    showblocksinterpretable(
        backend,
        getencodings(task),
        ("Encoded sample" => blocks.encodedsample, "Output" => blocks.ŷ),
        collect(zip(encsamples, outputs)),
    )
end

showoutputs(task::AbstractBlockTask, args...) =
    showoutputs(default_showbackend(), task, args...)

"""
    showoutputbatch([backend], task, outputbatch)
    showoutputbatch([backend], task, batch, outputbatch)

Show collated batch of outputs to `backend`. If a collated batch of encoded samples
`batch` is also given, show them next to the outputs. See [`showoutputs`](#) if you
have vectors of outputs and not collated batches.
"""
function showoutputbatch(backend::ShowBackend, task::AbstractBlockTask, outputbatch)
    outputs = collect(DataLoaders.obsslices(outputbatch))
    return showoutputs(backend, task, outputs)
end
function showoutputbatch(
    backend::ShowBackend,
    task::AbstractBlockTask,
    batch,
    outputbatch,
)
    encsamples = collect(DataLoaders.obsslices(batch))
    outputs = collect(DataLoaders.obsslices(outputbatch))
    return showoutputs(backend, task, encsamples, outputs)
end

showoutputbatch(task::AbstractBlockTask, args...) =
    showoutputbatch(default_showbackend(), task, args...)

# Testing helper

"""
    test_task_show(task, backend::ShowBackend)

Test suite that tests that all learning task-related `show*` functions
work for `backend`

## Keyword arguments

- `sample = mockblock(getblocks(task))`: Sample data to use for tests.
- `output = mockblock(getblocks(task).ŷ)`: Model output data to use for tests.
"""
function test_task_show(
    task::LearningTask,
    backend::ShowBackend;
    sample = mockblock(getblocks(task).sample),
    output = mockblock(getblocks(task).ŷ),
    context = Training(),
)

    encsample = encodesample(task, context, sample)
    pred = decodeypred(task, context, output)

    Test.@testset "`show*` test suite for learning task $task and backend $(typeof(backend))" begin
        @test (showsample(backend, task, sample); true)
        @test (showencodedsample(backend, task, encsample); true)
        @test (showoutput(backend, task, output); true)
        @test (showoutput(backend, task, encsample, output); true)
        @test (showprediction(backend, task, pred); true)
        @test (showprediction(backend, task, sample, pred); true)
    end
end

# Deprecations


text:
function endinvariants!(W1end::Vector{Float64},Q::Float64,A::Float64,Qm1::Float64,Am1::Float64,
    c0::Float64,beta::Float64,rho::Float64,mu::Float64,diff::Float64,k::Float64,h::Float64)
    lf = Float64;
    W1 = Vector{Float64}[];
    dW1 = Float64;
    dQ = Float64;
    dA = Float64;
    Qint = Float64;
    Aint = Float64;

    # forward wave speed
    lf = (Q./A .+ sqrt.(0.5.*beta./rho).*A.^0.25);
    # right-running invariants at current time
    W1 = [Qm1./Am1 .+ 4.*(sqrt.(0.5.*beta./rho).*Am1.^0.25 .- c0)];
    push!(W1,Q./A .+ 4.*(sqrt.(0.5.*beta./rho).*A.^0.25 .- c0));
    dW1 = (W1[2] .- W1[1])./k;
    # change in solution variables
    dQ = (Qm1 .- Q)./k;
    dA = (Am1 .- A)./k;
    # interpolated solution variables
    Qint = Qm1 .- dQ.*lf.*h;
    Aint = Am1 .- dA.*lf.*h;
    # update right-running invariant by linear interp.
    W1end[1] = W1[2] .- dW1.*lf.*h .- diff.*mu./rho.*h.*(Qint./Aint.^2);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
@testset "1552.magnetic-force-between-two-balls.jl" begin
    @test max_distance([1,2,3,4,7], 3) == 3
    @test max_distance([5,4,3,2,1,1000000000], 2) == 999999999
    @test max_distance([153, 66, 174, 97, 76, 101, 189, 47, 94, 141, 30, 135, 112, 15, 59, 134, 52, 108, 10, 31, 166, 160, 99, 146, 28, 58, 100, 18, 109, 9, 11, 54, 132, 113, 80, 88, 90, 46, 7, 145, 56, 33, 6, 79, 82, 157, 67, 199, 25, 2, 4, 178, 23, 57, 32, 151, 156, 73, 142, 116, 187, 77, 163, 38, 179, 154, 26, 122, 62, 136, 186, 87, 162, 159, 98, 117, 42, 161, 14], 15) == 12
    @test max_distance([153, 66, 174, 97, 76, 101, 189, 47, 94, 141, 30, 135, 112, 15, 59, 134, 52, 108, 10, 31, 166, 160, 99, 146, 28, 58, 100, 18, 109, 9, 11, 54, 132, 113, 80, 88, 90, 46, 7, 145, 56, 33, 6, 79, 82, 157, 67, 199, 25, 2, 4, 178, 23, 57, 32, 151, 156, 73, 142, 116, 187, 77, 163, 38, 179, 154, 26, 122, 62, 136, 186, 87, 162, 159, 98, 117, 42, 161, 14], 20) == 8
end


text:
using Test
using RadiiPolynomial, LinearAlgebra

include("vector_spaces.jl")

include("sequence.jl")

include("linear_operator.jl")

include("norm.jl")

include("arithmetic.jl")

include("special_operators.jl")


text:
function solve()
    n = parse(Int, readline())
    ret = 10^9 + 1
    for _=1:n
        a, p, x = [parse(Int, x) for x in split(readline())]
        if a < x
            ret = min(ret, p)
        end
    end
    ret == 10^9+1 ? -1 : ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module CovidModels

using SemanticModels
using OrdinaryDiffEq

include("odemodel.jl")
include("covid.jl")
include("data.jl")

end


text:
doc"""
    Erlang(α,θ)

The *Erlang distribution* is a special case of a [`Gamma`](:func:`Gamma`) distribution with integer shape parameter.

```julia
Erlang()       # Erlang distribution with unit shape and unit scale, i.e. Erlang(1.0, 1.0)
Erlang(a)      # Erlang distribution with shape parameter a and unit scale, i.e. Erlang(a, 1.0)
Erlang(a, s)   # Erlang distribution with shape parameter a and scale b
```

External links

* [Erlang distribution on Wikipedia](http://en.wikipedia.org/wiki/Erlang_distribution)

"""
immutable Erlang <: ContinuousUnivariateDistribution
    α::Int
    θ::Float64

    function Erlang(α::Real, θ::Real)
        @check_args(Erlang, isinteger(α) && α >= zero(α))
        new(α, θ)
    end
    Erlang(α::Real) = Erlang(α, 1.0)
    Erlang() = new(1, 1.0)
end

@distr_support Erlang 0.0 Inf

#### Parameters

shape(d::Erlang) = d.α
scale(d::Erlang) = d.θ
rate(d::Erlang) = inv(d.θ)
params(d::Erlang) = (d.α, d.θ)

#### Statistics

mean(d::Erlang) = d.α * d.θ
var(d::Erlang) = d.α * d.θ^2
skewness(d::Erlang) = 2.0 / sqrt(d.α)
kurtosis(d::Erlang) = 6.0 / d.α

function mode(d::Erlang)
    (α, θ) = params(d)
    α >= 1.0 ? θ * (α - 1.0) : error("Erlang has no mode when α < 1.0")
end

function entropy(d::Erlang)
    (α, θ) = params(d)
    α + lgamma(α) + (1.0 - α) * digamma(α) + log(θ)
end

mgf(d::Erlang, t::Real) = (1.0 - t * d.θ)^(-d.α)
cf(d::Erlang, t::Real)  = (1.0 - im * t * d.θ)^(-d.α)

#### Evaluation & Sampling

@_delegate_statsfuns Erlang gamma α θ

rand(d::Erlang) = StatsFuns.Rmath.gammarand(d.α, d.θ)


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libcurl"], :libcurl),
]

# Install BinaryBuilder dependencies
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaWeb/MbedTLSBuilder/releases/download/v0.11/build_MbedTLS.v1.0.0.jl"
]

for url in dependencies
    build_file = joinpath(@__DIR__, basename(url))
    if !isfile(build_file)
        download(url, build_file)
    end
end

# Execute the build scripts for the dependencies in an isolated module to avoid overwriting
# any variables/constants here
for url in dependencies
    build_file = joinpath(@__DIR__, basename(url))
    m = @eval module $(gensym()); include($build_file); end
    append!(products, m.products)
end

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaWeb/LibCURLBuilder/releases/download/v0.3.1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/LibCURL.v7.61.0.aarch64-linux-gnu.tar.gz", "69e02dcc8f2290551593ed68eac37615daf1d1bf185b940c105a93f5d95d6546"),
    Linux(:aarch64, :musl) => ("$bin_prefix/LibCURL.v7.61.0.aarch64-linux-musl.tar.gz", "eb562b8ebd0ef7ab3ddb5726049c12fd7666a40e93d76ef66bb53bbd5ab3b051"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/LibCURL.v7.61.0.arm-linux-gnueabihf.tar.gz", "4a08bab5a6d85ee65bf59c05d1fdf3617552a1c0989f3eac4d535ac23a306f62"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/LibCURL.v7.61.0.arm-linux-musleabihf.tar.gz", "154285daba1e01de10e38742fcf562b704df0126a369c1b78f829db06f08ab77"),
    Linux(:i686, :glibc) => ("$bin_prefix/LibCURL.v7.61.0.i686-linux-gnu.tar.gz", "7573fb35ff104c9badfa6f928e1c30dfe3c01eda72bd3a28773a87da42297249"),
    Linux(:i686, :musl) => ("$bin_prefix/LibCURL.v7.61.0.i686-linux-musl.tar.gz", "8c5cbadbbe39f9c6fa192edac58cdb4af8a29f1e260a2466a25fd5a7e87ae786"),
    Windows(:i686) => ("$bin_prefix/LibCURL.v7.61.0.i686-w64-mingw32.tar.gz", "fe41ca4f12ef8e720047a079076848b335579a006f6ea7e4a30b758484dcbb4d"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/LibCURL.v7.61.0.powerpc64le-linux-gnu.tar.gz", "3288401b57b46e6fc2f307ce6aa1c93b4d449764155ff33fcdbaeeed58b218c9"),
    MacOS(:x86_64) => ("$bin_prefix/LibCURL.v7.61.0.x86_64-apple-darwin14.tar.gz", "c1169390ff0dad5bf104a0645e8638c636b6f694d3a5ce10ad31384b0eeb2c3b"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/LibCURL.v7.61.0.x86_64-linux-gnu.tar.gz", "e1fb2a82877b77363630b5fde1097a26cc84ec7ac497f80908c50df21d3a5e2d"),
    Linux(:x86_64, :musl) => ("$bin_prefix/LibCURL.v7.61.0.x86_64-linux-musl.tar.gz", "7f45b7f955f80b68ae9430c7ae081626f916515b6ccaf25550655149b43c4382"),
    FreeBSD(:x86_64) => ("$bin_prefix/LibCURL.v7.61.0.x86_64-unknown-freebsd11.1.tar.gz", "c3945f2871f621d825203da419e274e490fab0b6e7ed0b1fb34908a173681a42"),
    Windows(:x86_64) => ("$bin_prefix/LibCURL.v7.61.0.x86_64-w64-mingw32.tar.gz", "6a3fb40bea8f2c7987ef46af8a423abfdf60287f69397f0f1df8973ea66e17fa"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
sir_model = @reaction_network rn begin
    0.1/1000, s + i --> 2i
    0.01, i --> r
end

sir_prob = DiscreteProblem([999,1,0],(0.0,250.0))
sir_jump_prob = JumpProblem(sir_prob,Direct(),sir_model)
sir_sol = solve(sir_jump_prob,FunctionMap())

prob = DiscreteProblem([999,1,0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),sir_model)

Random.seed!(100)
sol = solve(jump_prob,FunctionMap())

nums = Int[]
@time for i in 1:1000
  jump_prob = JumpProblem(prob,Direct(),sir_model)
  sol = solve(jump_prob,FunctionMap())
  push!(nums,sol[end][3])
end
@test mean(nums) - 740 < 20

f = function (du,u,p,t)
  du[4] = u[2]*u[3]/100000 - u[1]*u[2]/100000
end

prob = ODEProblem(f,[999.0,1.0,0.0,100.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),sir_model)

sol = solve(jump_prob,Tsit5())


text:
export ParticleMeshCoupling1D

"""
    ParticleMeshCoupling1D( mesh, no_particles, spline_degree,
                          smoothing_type )

Kernel smoother with splines of arbitrary degree placed on a uniform mesh.
Spline with index i starts at point i

- `delta_x` : Value of grid spacing along both directions.
- `n_grid` : Array containing number ofpoints along each direction
- `no_particles` : Number of particles of underlying PIC method
- `spline_degree` : Degree of smoothing kernel spline
- `n_span` : Number of intervals where spline non zero (spline_degree + 1)
- `scaling` : Scaling factor depending on whether :galerkin or :collocation
- `n_quad_points` : Number of quadrature points
- `spline_val`: scratch data for spline evaluation
- `spline_val_more` : more scratch data for spline evaluation
- `quad_x, quad_w` : quadrature weights and points

!!! note

    Only 1D version is implemented for now

"""
struct ParticleMeshCoupling1D <: AbstractParticleMeshCoupling
    dims::Int
    delta_x::Float64
    n_grid::Int
    n_dofs::Int
    no_particles::Int
    spline_degree::Int
    n_span::Int
    scaling::Float64
    n_quad_points::Int
    spline_val::Vector{Float64}
    spline_val_more::Vector{Float64}
    quad_x::Vector{Float64}
    quad_w::Vector{Float64}
    spline_pp::SplinePP
    xmin::Float64
    Lx::Float64

    function ParticleMeshCoupling1D(
        mesh::OneDGrid, no_particles::Int, spline_degree::Int, smoothing_type::Symbol
    )
        dims = 1
        n_grid = mesh.nx
        n_dofs = n_grid
        xmin = mesh.xmin
        xmax = mesh.xmax
        Lx = xmax - xmin
        delta_x = (xmax - xmin) / n_grid
        n_span = spline_degree + 1

        if smoothing_type == :collocation
            scaling = 1.0 / delta_x
        elseif smoothing_type == :galerkin
            scaling = 1.0
        else
            throw(ArgumentError("""
                 Smoothing Type $smoothing_type not implemented
                 for kernel_smoother_spline_1d.
           """))
        end

        n_quad_points = (spline_degree + 2) ÷ 2

        spline_val = zeros(Float64, n_span)
        spline_val_more = zeros(Float64, n_span)

        quad_x, quad_w = gausslegendre(n_quad_points)

        spline_pp = SplinePP(spline_degree, n_grid[1])

        return new(
            dims,
            delta_x,
            n_grid,
            n_dofs,
            no_particles,
            spline_degree,
            n_span,
            scaling,
            n_quad_points,
            spline_val,
            spline_val_more,
            quad_x,
            quad_w,
            spline_pp,
            xmin,
            Lx,
        )
    end
end

"""
    add_charge_pp!(rho_dofs, p, position, marker_charge)

Add charge of one particle
- `p`             : kernel smoother object
- `position`      : Position of the particle
- `marker_charge` : Particle weights time charge
- `rho_dofs`      : Coefficient vector of the charge distribution
"""
function add_charge_pp!(
    rho_dofs::Vector{Float64}, p::ParticleMeshCoupling1D, position::Float64, marker_charge
)
    xi = (position - p.xmin) / p.delta_x
    index = trunc(Int, xi) + 1
    xi = xi - (index - 1)
    index = index - p.spline_degree

    @inbounds for i in eachindex(p.spline_val)
        p.spline_val[i] = horner_1d(p.spline_degree, p.spline_pp.poly_coeffs, xi, i)
    end

    @inbounds for i in 1:(p.n_span)
        index1d = mod(index + i - 2, p.n_grid[1]) + 1
        rho_dofs[index1d] += marker_charge * p.spline_val[i] * p.scaling
    end
end

"""
    add_current_update_v_pp!( j_dofs, p, position_old, position_new,
                              marker_charge, qoverm, bfield_dofs, vi)

Add current for one particle and update v
(according to `H_p1` part in Hamiltonian splitting)
"""
function add_current_update_v_pp!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    position_old,
    position_new,
    marker_charge::Float64,
    qoverm::Float64,
    bfield_dofs::Vector{Float64},
    vi::Vector{Float64},
)

    # Read out particle position and velocity
    # Compute index_old, the index of the last DoF on the grid the particle
    # contributes to, and r_old, its position (normalized to cell size one).

    xi = (position_old[1] - p.xmin) / p.delta_x[1]
    index_old = trunc(Int, xi)
    r_old = xi - index_old

    # Compute the new box index index_new and normalized position r_old.
    xi = (position_new[1] - p.xmin) / p.delta_x[1]
    index_new = trunc(Int, xi)
    r_new = xi - index_new

    if index_old == index_new
        vi[2] = update_jv_pp!(
            j_dofs, p, r_old, r_new, index_old, marker_charge, qoverm, vi[2], bfield_dofs
        )

    elseif index_old < index_new
        vi[2] = update_jv_pp!(
            j_dofs, p, r_old, 1.0, index_old, marker_charge, qoverm, vi[2], bfield_dofs
        )

        vi[2] = update_jv_pp!(
            j_dofs, p, 0.0, r_new, index_new, marker_charge, qoverm, vi[2], bfield_dofs
        )

        for ind in (index_old + 1):(index_new - 1)
            vi[2] = update_jv_pp!(
                j_dofs, p, 0.0, 1.0, ind, marker_charge, qoverm, vi[2], bfield_dofs
            )
        end

    else
        vi[2] = update_jv_pp!(
            j_dofs, p, 1.0, r_new, index_new, marker_charge, qoverm, vi[2], bfield_dofs
        )

        vi[2] = update_jv_pp!(
            j_dofs, p, r_old, 0.0, index_old, marker_charge, qoverm, vi[2], bfield_dofs
        )

        for ind in (index_new + 1):(index_old - 1)
            vi[2] = update_jv_pp!(
                j_dofs, p, 1.0, 0.0, ind, marker_charge, qoverm, vi[2], bfield_dofs
            )
        end
    end
end

"""
    update_jv_pp!( j_dofs, p, lower, upper, index, marker_charge,
                   qoverm, vi, bfield_dofs)

Helper function for `add_current_update_v`.
"""
function update_jv_pp!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    lower::Float64,
    upper::Float64,
    index::Int64,
    marker_charge::Float64,
    qoverm::Float64,
    vi::Float64,
    bfield_dofs::Vector{Float64},
)
    n_cells = p.n_grid[1]

    # Evaluation of the primitive integral at the lower
    # and upper bound of the gridcell

    horner_primitive_1d(p.spline_val, p.spline_degree, p.spline_pp.poly_coeffs_fp, lower)
    horner_primitive_1d(
        p.spline_val_more, p.spline_degree, p.spline_pp.poly_coeffs_fp, upper
    )

    p.spline_val .= (p.spline_val_more .- p.spline_val) .* (p.delta_x[1])

    ind = 1
    @inbounds for i_grid in (index - p.spline_degree):index
        i_mod = mod(i_grid, n_cells) + 1
        j_dofs[i_mod] = j_dofs[i_mod] + (marker_charge * p.spline_val[ind] * p.scaling)
        vi = vi - qoverm * p.spline_val[ind] * bfield_dofs[i_mod]
        ind = ind + 1
    end

    return vi
end

"""
    evaluate_pp(p, position, field_dofs_pp)

Evaluate field at `position` using horner scheme
- `p` : Kernel smoother object
- `position` : Position of the particle
- `field_dofs_pp` : Degrees of freedom in kernel representation.
- `field_value` : Value(s) of the electric fields at given position

"""
function evaluate_pp(
    p::ParticleMeshCoupling1D, position::Float64, field_dofs_pp::Array{Float64,2}
)
    xi = (position - p.xmin) / p.delta_x
    index = trunc(Int, xi) + 1
    xi = xi - (index - 1)

    return horner_1d(p.spline_degree, field_dofs_pp, xi, index)
end

"""
    add_charge!( rho, p, position, marker_charge)

Add charge of one particle
- `p`             : kernel smoother object
- `position`      : Position of the particle
- `marker_charge` : Particle weights time charge
- `rho_dofs`      : Coefficient vector of the charge distribution
"""
function add_charge!(
    rho_dofs::Vector{Float64},
    p::ParticleMeshCoupling1D,
    position::Float64,
    marker_charge::Float64,
)
    xi = (position - p.xmin) / p.delta_x
    index = trunc(Int, xi)
    xi = xi - index
    index = index - p.spline_degree

    uniform_bsplines_eval_basis!(p.spline_val, p.spline_degree, xi)

    nx = p.n_grid[1]

    @inbounds for i in 1:(p.n_span)
        index1d = mod1(index + i, nx)
        rho_dofs[index1d] += marker_charge * p.spline_val[i] * p.scaling
    end
end

"""
    add_current_update_v!( j_dofs, p,
                           position_old, position_new,
                           marker_charge, qoverm,
                           bfield_dofs, vi)

Add current for one particle and update v (according to ``H_{p1}``
part in Hamiltonian splitting)

- Read out particle position and velocity
- Compute index_old, the index of the last DoF on the grid the
particle contributes to, and `r_old`, its position (normalized to cell size one).

"""
function add_current_update_v!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    position_old::Float64,
    position_new::Float64,
    marker_charge::Float64,
    qoverm::Float64,
    bfield_dofs::Vector{Float64},
    vi::Float64,
)
    xi = (position_old - p.xmin) / p.delta_x
    index_old = trunc(Int, xi)
    r_old = xi - index_old

    # Compute the new box index index_new and normalized position r_old.

    xi = (position_new - p.xmin) / p.delta_x
    index_new = trunc(Int, xi)
    r_new = xi - index_new

    if index_old == index_new
        if r_old < r_new
            vi = update_jv!(
                j_dofs,
                p,
                r_old,
                r_new,
                index_old,
                marker_charge,
                qoverm,
                1.0,
                vi,
                bfield_dofs,
            )
        else
            vi = update_jv!(
                j_dofs,
                p,
                r_new,
                r_old,
                index_old,
                marker_charge,
                qoverm,
                -1.0,
                vi,
                bfield_dofs,
            )
        end

    elseif index_old < index_new
        vi = update_jv!(
            j_dofs, p, r_old, 1.0, index_old, marker_charge, qoverm, 1.0, vi, bfield_dofs
        )

        vi = update_jv!(
            j_dofs, p, 0.0, r_new, index_new, marker_charge, qoverm, 1.0, vi, bfield_dofs
        )

        for ind in (index_old + 1):(index_new - 1)
            vi = update_jv!(
                j_dofs, p, 0.0, 1.0, ind, marker_charge, qoverm, 1.0, vi, bfield_dofs
            )
        end

    else
        vi = update_jv!(
            j_dofs, p, r_new, 1.0, index_new, marker_charge, qoverm, -1.0, vi, bfield_dofs
        )
        vi = update_jv!(
            j_dofs, p, 0.0, r_old, index_old, marker_charge, qoverm, -1.0, vi, bfield_dofs
        )

        for ind in (index_new + 1):(index_old - 1)
            vi = update_jv!(
                j_dofs, p, 0.0, 1.0, ind, marker_charge, qoverm, -1.0, vi, bfield_dofs
            )
        end
    end

    return vi
end

"""
    update_jv!(j_dofs, p,
               lower, upper, index, marker_charge,
               qoverm, sign, vi, bfield_dofs)

Helper function for `add_current_update_v`.
"""
function update_jv!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    lower::Float64,
    upper::Float64,
    index::Int64,
    marker_charge::Float64,
    qoverm::Float64,
    sign::Float64,
    vi::Float64,
    bfield_dofs::Vector{Float64},
)
    n_cells = p.n_grid[1]

    c1 = 0.5 * (upper - lower)
    c2 = 0.5 * (upper + lower)

    uniform_bsplines_eval_basis!(p.spline_val, p.spline_degree, c1 * p.quad_x[1] + c2)

    p.spline_val .*= p.quad_w[1] * c1

    @inbounds for j in 2:(p.n_quad_points)
        uniform_bsplines_eval_basis!(
            p.spline_val_more, p.spline_degree, c1 * p.quad_x[j] + c2
        )

        p.spline_val .+= p.spline_val_more .* p.quad_w[j] .* c1
    end

    p.spline_val .*= sign * p.delta_x[1]

    ind = 1
    @inbounds for i_grid in (index - p.spline_degree):index
        i_mod = mod(i_grid, n_cells) + 1
        j_dofs[i_mod] += marker_charge * p.spline_val[ind] * p.scaling
        vi = vi - qoverm * p.spline_val[ind] * bfield_dofs[i_mod]
        ind = ind + 1
    end

    return vi
end

"""
    evaluate(p, position, field_dofs)

Evaluate field at `position`

- `p` : Kernel smoother object
- `position` : Position of the particle
- `field_dofs` : Coefficient vector for the field DoFs
- `field_value` : Value(s) of the electric fields at given position

"""
function evaluate(p::ParticleMeshCoupling1D, position::Float64, field_dofs::Vector{Float64})
    xi = (position - p.xmin) / p.delta_x
    index = trunc(Int, xi)
    dxi = xi - index
    index = index - p.spline_degree
    uniform_bsplines_eval_basis!(p.spline_val, p.spline_degree, dxi)

    nx = p.n_grid[1]
    field_value = 0.0
    @inbounds for i in 1:(p.n_span)
        index1d = mod1(index + i, nx)
        field_value += field_dofs[index1d] * p.spline_val[i]
    end

    return field_value
end


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

mutable struct BillingProfile <: SwaggerModel
    maxPrice::Any # spec type: Union{ Nothing, Float64 } # spec name: maxPrice

    function BillingProfile(;maxPrice=nothing)
        o = new()
        validate_property(BillingProfile, Symbol("maxPrice"), maxPrice)
        setfield!(o, Symbol("maxPrice"), maxPrice)
        o
    end
end # type BillingProfile

const _property_map_BillingProfile = Dict{Symbol,Symbol}(Symbol("maxPrice")=>Symbol("maxPrice"))
const _property_types_BillingProfile = Dict{Symbol,String}(Symbol("maxPrice")=>"Float64")
Base.propertynames(::Type{ BillingProfile }) = collect(keys(_property_map_BillingProfile))
Swagger.property_type(::Type{ BillingProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BillingProfile[name]))}
Swagger.field_name(::Type{ BillingProfile }, property_name::Symbol) =  _property_map_BillingProfile[property_name]

function check_required(o::BillingProfile)
    true
end

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


text:
include("../benchmark.jl")
include("create_report.jl")

results = get_CUTEst_results()
#["one_phase"] = load("../results/one_phase/Dec10/summary.jld", "summary")

using DataFrames,CSV

df_results = Dict()
for (method_name,method_results) in results
  df_results[method_name] = DataFrame(name=Array{String,1}(),it=Array{Int64,1}(),time=Array{Float64,1}(),fval=Array{Float64,1}(),con=Array{Float64,1}(),status=[])
  for (name,info) in method_results
    push!(df_results[method_name],[name,info.it_count,info.total_time,info.fval,info.con_vio,string(info.status)])
  end
  CSV.write("$folder/table_$(folder)_$(method_name).csv", df_results[method_name])
end

###
### ITER/TIMING BREAK DOWN
###

df_summary = DataFrame(name=[],median_it=[],median_time=[],total_time=[])
for (method_name,method_results) in results
    push!(df_summary,[method_name,median(df_results[method_name][:it]),median(df_results[method_name][:time]),sum(df_results[method_name][:time])])
end

df_summary

CSV.write("$folder/runtime_summary.csv", df_summary)

using FreqTables

##
## FAILURE REASONS
##
freqtable(df_results["ipopt"][:status])
freqtable(df_results["one phase"][:status])


text:
using ClearswiApp
using Test

@testset "ClearswiApp.jl" begin

niread = ClearswiApp.niread
savenii = ClearswiApp.savenii

p = joinpath("data", "small")
phasefile_me = joinpath(p, "Phase.nii")
phasefile_me_nan = joinpath(p, "phase_with_nan.nii")
magfile_me = joinpath(p, "Mag.nii")
tmpdir = mktempdir()
phasefile_1eco = joinpath(tmpdir, "Phase.nii")
magfile_1eco = joinpath(tmpdir, "Mag.nii")
phasefile_1arreco = joinpath(tmpdir, "Phase.nii")
magfile_1arreco = joinpath(tmpdir, "Mag.nii")
magfile_me_nan_size = joinpath(tmpdir, "mag_nan_size.nii")
maskfile = joinpath(tmpdir, "Mask.nii")
savenii(niread(magfile_me)[:,:,:,1] |> I -> I .> ClearswiApp.MriResearchTools.median(I), maskfile)
savenii(niread(phasefile_me)[:,:,:,1], phasefile_1eco)
savenii(niread(magfile_me)[:,:,:,1], magfile_1eco)
savenii(niread(phasefile_me)[:,:,:,[1]], phasefile_1arreco)
savenii(niread(magfile_me)[:,:,:,[1]], magfile_1arreco)
savenii(ones(size(niread(phasefile_me_nan))), magfile_me_nan_size)

function test_clearswi(args)
    folder = tempname()
    args = [args..., "-o", folder]
    try
        msg = clearswi_main(args)
        @test msg == 0
        @test isfile(joinpath(folder, "clearswi.nii"))
    catch e
        println(args)
        println(sprint(showerror, e, catch_backtrace()))
        @test "test failed" == "with error" # signal a failed test
    end
end

configurations_se(pf, mf) = configurations_se(["-p", pf, "-m", mf])
configurations_se(pm) = [
    [pm...],
    [pm..., "--mag-combine", "SNR"],
    [pm..., "--mag-combine", "average"],
    [pm..., "--mag-combine", "echo", "2"],
    [pm..., "--mag-combine", "SE", "2.4"],
    [pm..., "--mag-sensitivity-correction", "off"],
    [pm..., "--mag-softplus-scaling", "off"],
    [pm..., "--unwrapping-algorithm", "romeo"],
    [pm..., "--unwrapping-algorithm", "laplacianslice"],
    [pm..., "--filter-size", "[2,2,3]"],
    [pm..., "--phase-scaling-type", "negativetanh"],
    [pm..., "--phase-scaling-type", "positive"],
    [pm..., "--phase-scaling-type", "negative"],
    [pm..., "--phase-scaling-type", "triangular"],
    [pm..., "--phase-scaling-strength", "1"],
    [pm..., "--phase-scaling-strength", "10"],
    [pm..., "-N"],
    [pm..., "--no-phase-rescale"],
    [pm..., "--writesteps", tmpdir],
]
configurations_me(phasefile_me, magfile_me) = configurations_me(["-p", phasefile_me, "-m", magfile_me])
configurations_me(pm) = [
    [pm..., "-e", "1:2", "-t", "[2,4]"], # giving two echo times for two echoes used out of three
    [pm..., "-e", "[1,3]", "-t", "[2,4,6]"], # giving three echo times for two echoes used out of three
    [pm..., "-e", "[1", "3]", "-t", "[2,4,6]"],
    [pm..., "-t", "[2,4,6]"],
    [pm..., "-t", "2:2:6"],
    [pm..., "-t", "[2.1,4.2,6.3]"],
]

files = [(phasefile_1eco, magfile_1eco), (phasefile_1arreco, magfile_1arreco), (phasefile_1eco, magfile_1arreco), (phasefile_1arreco, magfile_1eco)]
for (pf, mf) in files, args in configurations_se(pf, mf)
    test_clearswi(args)
end
for args in configurations_me(phasefile_me, magfile_me)
    test_clearswi(args)
end
for args in configurations_se(["-p", phasefile_me, "-m", magfile_me, "-t", "[2,4,6]"])
    test_clearswi(args)
end

test_clearswi(["-p", phasefile_1eco, "-m", magfile_1eco, "-t", "5"])
#test_clearswi(["-p", phasefile_me_nan, "-m", magfile_me_nan_size, "-t", "[2,4]"])

end

## print version to verify
println()
clearswi_main(["--version"])
#=
## Test error and warning messages
m = "multi-echo data is used, but no echo times are given. Please specify the echo times using the -t option."
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v"])

m = "masking option '0.8' is undefined (Maybe '-k qualitymask 0.8' was meant?)"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-v", "-k", "0.8"])

m = "masking option 'blub' is undefined"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-v", "-k", "blub"])

m = "Phase offset determination requires all echo times!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me_5D, "-o", tmpdir, "-v", "-t", "[1,2]", "-e", "[1,2]", "--phase-offset-correction"])

m = "5D phase is given but no coil combination is selected"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me_5D, "-o", tmpdir, "-v", "-t", "[1,2,3]"])

m = "echoes=[1,5]: specified echo out of range! Number of echoes is 3"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-e", "[1,5]"])

m = "echoes=[1,5} wrongly formatted!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-e", "[1,5}"])

m = "Number of chosen echoes is 2 (3 in .nii data), but 5 TEs were specified!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3,4,5]", "-e", "[1,2]"])

m = "size of magnitude and phase does not match!"
@test_throws ErrorException(m) unwrapping_main(["-p", phasefile_me, "-o", tmpdir, "-v", "-t", "[1,2,3]", "-m", magfile_1eco])

m = "robustmask was chosen but no magnitude is available. No mask is used!"
@test_logs (:warn, m) match_mode=:any unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir])

@test_logs unwrapping_main(["-p", phasefile_1eco, "-o", tmpdir, "-m", magfile_1eco]) # test that no warning appears

end

=#


text:
using Yao
using BenchmarkTools
using Statistics, DelimitedFiles

const NL = 10:3:25

function bgate(ops, filename)
    println("Benchmarking $filename ...")
    res_list = zeros(6, ops |> length)
    for (ib, nbit) in enumerate(NL)
        reg = rand_state(nbit)
        for (ig, G) in enumerate(ops)
            println(G(nbit))
            res = @benchmark $reg |> $(G(nbit))
            res_list[ib, ig] = median(res).time
        end
    end
    println(res_list)
    writedlm(filename, res_list/1e3)
end

function bench_xyz()
    #bgate([repeat(G, [3]) for G in [X,Y,Z]], "xyz-report.dat")
    bgate([put(3=>G) for G in [X,Y,Z]], "xyz-report.dat")
end

function bench_cxyz()
    bgate([control((7,), 3=>G) for G in [X, Y, Z]], "cxyz-report.dat")
end

function bench_repeatxyz()
    bgate([repeat(G, collect(2:7)) for G in [X, Y, Z]], "repeatxyz-report.dat")
end

function bench_hgate()
   bgate([repeat(H, [3]), control((7,), 3=>H), repeat(H, collect(2:7))], "h-report.dat")
end

function bench_hgate()
    bgate([put(3=>H), control((7,), 3=>H), repeat(H, 2:7)], "h-report.dat")
end

function bench_rot()
    gates = [repeat(rot(G, 0.5), [3]) for G in [X, Y, Z]]
    bgate(gates, "rot-report.dat")
end

function bench_crot()
    gates = [control((7,), 3=>rot(G, 0.5)) for G in [X, Y, Z]]
    bgate(gates, "crot-report.dat")
end

function bench_toffoli()
    bgate([control((2,3), 5=>X)], "toffoli-report.dat")
end

function bench_all()
    bench_xyz()
    bench_cxyz()
    bench_repeatxyz()
    bench_hgate()
    bench_rot()
    bench_crot()
    bench_toffoli()
end

#bench_all()
bench_rot()
bench_crot()
bench_toffoli()
#bench_hgate()


text:
using Images
using FileIO

using Plots
pyplot()

using AkazeFeatures

# imagename = "images/dirac.png"
# imagename = "images/square.png"
# imagename = "images/concave.png"
# imagename = "images/wiggly.png"
# imagename = "images/wiggly-blur.png"
# imagename = "images/pirate.png"
imagename = "/home/user/src/data/mars/NLB_629553760EDR_F0780612CCAM05613M_.jpeg"
# imagename = "../vggaffine/boat/img1.pgm"
# imagename = "images/cameraman.png"
# imagename = "images/zezao.png"
# imagename = "images/polly.png"
# imagename = "images/polly-small.png"

oriimg = load(imagename)
# img = rawview(channelview(oriimg)) / 255

img = load_image_as_grayscale(imagename)

img_height, img_width = size(img)

opt = AKAZEOptions(
    omin = 3,
    omax = 4,
    nsublevels=2,
    img_width = img_width,
    img_height = img_height,
    # diffusivity = PM_G1,
    # diffusivity = PM_G2,
    # diffusivity = WEICKERT,
    diffusivity = CHARBONNIER,
    # dthreshold = 16e-3,
    # dthreshold = 8e-3,
    # dthreshold = 4e-3,
    dthreshold = 2e-3,
)
akaze = AKAZE(opt)

Create_Nonlinear_Scale_Space(akaze, img)
kpts = Feature_Detection(akaze)
Compute_Main_Orientation.([akaze], kpts)
println("Extracted $(length(kpts)) points")
desc = Compute_Descriptors(akaze, kpts)

origpt, origdesc = original_akaze_features(imagename, diffusivity=Int(opt.diffusivity), nsublevels=opt.nsublevels, omax=opt.omax, dthreshold=opt.dthreshold)
println("Original: $(length(origpt)) points")

plot(size=(800,800))
plot!([-0.5, size(img)[2]-0.5], [-0.5,size(img)[1]-0.5], RGB.(oriimg), yflip = true)
plot_features(origpt, "Blues", :solid);
plot_features(kpts[1:length(kpts)], "Oranges");
plot!(xticks = :native, yticks = :native)
plot!()
plot!(yflip = false)


text: triangle(n::Int) = n >= 0 ? [[binomial(i, j) for j in 0:i] for i in 0:n-1] : throw(DomainError())


text:
module Router

using Dates, HTTP, JSON3
using ..DataModel, ..Service, ..Auth, ..Workers # ..Contexts

const ROUTER = HTTP.Router()

# ------------------------------------------------------------ #
# comment sample code
# ------------------------------------------------------------ #
# createAlbum(req) = Service.createAlbum(JSON3.read(req.body))::Album
# HTTP.@register(ROUTER, "POST", "/album", createAlbum)

# getAlbum(req) = Service.getAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]))::Album
# HTTP.@register(ROUTER, "GET", "/album/*", getAlbum)

# updateAlbum(req) = Service.updateAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]), JSON3.read(req.body, Album))::Album
# HTTP.@register(ROUTER, "PUT", "/album/*", updateAlbum)

# deleteAlbum(req) = Service.deleteAlbum(parse(Int, HTTP.URIs.splitpath(req.target)[2]))
# HTTP.@register(ROUTER, "DELETE", "/album/*", deleteAlbum)

# pickAlbumToListen(req) = fetch(Workers.@async(Service.pickAlbumToListen()::Album))
# HTTP.@register(ROUTER, "GET", "/", pickAlbumToListen)

# ------------------------------------------------------------ #
# RestAPI routers
# ------------------------------------------------------------ #
getSystemAlerts(req) = fetch(Workers.@async(Service.getSystemAlerts()::SystemAlert))
HTTP.@register(ROUTER, "GET", "/sysai", getSystemAlerts)

# ------------------------------------------------------------ #
# this is demo app, call simple router instead of AUTH_ROUTER
# authentication is skipped
# ------------------------------------------------------------ #

# function contextHandler(req)
#     withcontext(User(req)) do
#         HTTP.Response(200, JSON3.write(HTTP.handle(ROUTER, req)))
#     end
# end

# const AUTH_ROUTER = HTTP.Router(contextHandler)

# function authenticate(user::User)
#     resp = HTTP.Response(200, JSON3.write(user))
#     return Auth.addtoken!(resp, user)
# end

# createUser(req) = authenticate(Service.createUser(JSON3.read(req.body))::User)
# HTTP.@register(AUTH_ROUTER, "POST", "/user", createUser)

# loginUser(req) = authenticate(Service.loginUser(JSON3.read(req.body, User))::User)
# HTTP.@register(AUTH_ROUTER, "POST", "/user/login", loginUser)

function requestHandler(req)
    start = Dates.now(Dates.UTC)
    @info (timestamp = start, event = "ServiceRequestBegin", tid = Threads.threadid(), method = req.method, target = req.target)
    local resp
    try
        # ------------------------------------------------------------ #
        # this is demo app, call simple router instead of AUTH_ROUTER
        # authentication is skipped
        # ------------------------------------------------------------ #

        # resp = HTTP.handle(AUTH_ROUTER, req)
        resp = HTTP.handle(ROUTER, req)
    catch e
        if e isa Auth.Unauthenticated
            resp = HTTP.Response(401)
        else
            s = IOBuffer()
            showerror(s, e, catch_backtrace(); backtrace = true)
            errormsg = String(resize!(s.data, s.size))
            @error errormsg
            resp = HTTP.Response(500, errormsg)
        end
    end
    stop = Dates.now(Dates.UTC)
    @info (timestamp = stop, event = "ServiceRequestEnd", tid = Threads.threadid(), method = req.method, target = req.target, duration = Dates.value(stop - start), status = resp.status, bodysize = length(resp.body))
    return resp
end

function run()
    HTTP.serve(requestHandler, "0.0.0.0", 8080)
end

end # module


text:
NIScope_lib_found = Libdl.find_library(libniscope)
if ~isempty(NIScope_lib_found)
#const lib = Libdl.dlopen(libniscope) # Open the library explicitly. not working with dlsym in packages...
else
@warn "NI PXI libraries not found!\nPXI functionality broken\nPlease install NI drivers!"
end
#Libdl.dlclose(lib) # close dll
#sym = Libdl.dlsym(lib, :niScope_init)
#ccall(sym, ...) # Use the pointer `sym` instead of the (symbol, library) tuple
# ccall basic:
# ccall(:func, "dllname", return_type,(input types) ,input arguments)

include("wrapper/IviScope_common.jl")
include("wrapper/IviScope_api.jl")
include("wrapper/niScope_common.jl")
include("wrapper/niScope_api.jl")
include("wrapper/niScopeCal_api.jl")
include("wrapper/niScope_api_dispatcher.jl")

include("NIScope_Acquisition.jl")
include("NIScope_high_level_built_in.jl")
include("NIScope_Horizontal.jl")
include("NIScope_Math.jl")
include("NIScope_Measurements.jl")
include("NIScope_Misc.jl")
include("NIScope_Synchronization.jl")
include("NIScope_Triggering.jl")
include("NIScope_Vertical.jl")


text:
"""
$(TYPEDEF)

A scalar equation for optimization.

# Fields
$(FIELDS)

# Examples

```
@variables x y z
@parameters σ ρ β

op = σ*(y-x) + x*(ρ-z)-y + x*y - β*z
os = OptimizationSystem(eqs, [x,y,z],[σ,ρ,β])
```
"""
struct OptimizationSystem <: AbstractSystem
    """Vector of equations defining the system."""
    op::Any
    """Unknown variables."""
    states::Vector
    """Parameters."""
    ps::Vector
    observed::Vector{Equation}
    equality_constraints::Vector{Equation}
    inequality_constraints::Vector
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{OptimizationSystem}
    """
    default_u0: The default initial conditions to use when initial conditions
    are not supplied in `ODEProblem`.
    """
    default_u0::Dict
    """
    default_p: The default parameters to use when parameters are not supplied
    in `ODEProblem`.
    """
    default_p::Dict
end

function OptimizationSystem(op, states, ps;
                            observed = [],
                            equality_constraints = Equation[],
                            inequality_constraints = [],
                            default_u0=Dict(),
                            default_p=Dict(),
                            name = gensym(:OptimizationSystem),
                            systems = OptimizationSystem[])

    default_u0 isa Dict || (default_u0 = Dict(default_u0))
    default_p isa Dict || (default_p = Dict(default_p))
    default_u0 = Dict(value(k) => value(default_u0[k]) for k in keys(default_u0))
    default_p = Dict(value(k) => value(default_p[k]) for k in keys(default_p))

    OptimizationSystem(
                       value(op), value.(states), value.(ps),
                       observed,
                       equality_constraints, inequality_constraints,
                       name, systems, default_u0, default_p
                      )
end

function calculate_gradient(sys::OptimizationSystem)
    expand_derivatives.(gradient(equations(sys), states(sys)))
end

function generate_gradient(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys); kwargs...)
    grad = calculate_gradient(sys)
    return build_function(grad, vs, ps;
                          conv = AbstractSysToExpr(sys),kwargs...)
end

function calculate_hessian(sys::OptimizationSystem)
    expand_derivatives.(hessian(equations(sys), states(sys)))
end

function generate_hessian(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys);
                          sparse = false, kwargs...)
    if sparse
        hess = sparsehessian(equations(sys),states(sys))
    else
        hess = calculate_hessian(sys)
    end
    return build_function(hess, vs, ps;
                          conv = AbstractSysToExpr(sys),kwargs...)
end

function generate_function(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys); kwargs...)
    return build_function(equations(sys), vs, ps;
                          conv = AbstractSysToExpr(sys),kwargs...)
end

equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_expr.(sys.systems))
namespace_expr(sys::OptimizationSystem) = namespace_expr(sys.op,sys.name,nothing)

hessian_sparsity(sys::OptimizationSystem) =
    hessian_sparsity(sys.op, states(sys))

struct AutoModelingToolkit <: DiffEqBase.AbstractADType end

DiffEqBase.OptimizationProblem(sys::OptimizationSystem,args...;kwargs...) =
    DiffEqBase.OptimizationProblem{true}(sys::OptimizationSystem,args...;kwargs...)

"""
```julia
function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          grad = false,
                                          hess = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates an OptimizationProblem from an OptimizationSystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0,
                                          parammap=DiffEqBase.NullParameters();
                                          lb=nothing, ub=nothing,
                                          grad = false,
                                          hess = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)

    f = generate_function(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                              expression=Val{false})

    if grad
        grad_oop,grad_iip = generate_gradient(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                                  parallel=parallel,expression=Val{false})
        _grad(u,p) = grad_oop(u,p)
        _grad(J,u,p) = (grad_iip(J,u,p); J)
    else
        _grad = nothing
    end

    if hess
        hess_oop,hess_iip = generate_hessian(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                                 sparse=sparse,parallel=parallel,expression=Val{false})
       _hess(u,p) = hess_oop(u,p)
       _hess(J,u,p) = (hess_iip(J,u,p); J)
    else
        _hess = nothing
    end

    _f = DiffEqBase.OptimizationFunction{iip,AutoModelingToolkit,typeof(f),typeof(_grad),typeof(_hess),Nothing,Nothing,Nothing,Nothing}(f,AutoModelingToolkit(),_grad,_hess,nothing,nothing,nothing,nothing)

    u0 = varmap_to_vars(u0,dvs; defaults=default_u0(sys))
    p = varmap_to_vars(parammap,ps; defaults=default_p(sys))
    lb = varmap_to_vars(lb,dvs)
    ub = varmap_to_vars(ub,dvs)
    OptimizationProblem{iip}(_f,u0,p;lb=lb,ub=ub,kwargs...)
end

"""
```julia
function DiffEqBase.OptimizationProblemExpr{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          grad = false,
                                          hes = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates a Julia expression for an OptimizationProblem from an
OptimizationSystem and allows for automatically symbolically
calculating numerical enhancements.
"""
struct OptimizationProblemExpr{iip} end

OptimizationProblemExpr(sys::OptimizationSystem,args...;kwargs...) =
    OptimizationProblemExpr{true}(sys::OptimizationSystem,args...;kwargs...)

function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
                                          parammap=DiffEqBase.NullParameters();
                                          lb=nothing, ub=nothing,
                                          grad = false,
                                          hess = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = false, parallel=SerialForm(),
                                          kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)
    idx = iip ? 2 : 1
    f = generate_function(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                              expression=Val{true})
    if grad
        _grad = generate_gradient(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                             parallel=parallel,expression=Val{false})[idx]
    else
        _grad = :nothing
    end

    if hess
        _hess = generate_hessian(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                                         sparse=sparse,parallel=parallel,expression=Val{false})[idx]
    else
        _hess = :nothing
    end

    u0 = varmap_to_vars(u0,dvs; defaults=default_u0(sys))
    p = varmap_to_vars(parammap,ps; defaults=default_p(sys))
    lb = varmap_to_vars(lb,dvs)
    ub = varmap_to_vars(ub,dvs)
    quote
        f = $f
        p = $p
        u0 = $u0
        grad = $_grad
        hess = $_hess
        lb = $lb
        ub = $ub
        _f = OptimizationFunction{$iip,typeof(f),typeof(grad),typeof(hess),Nothing,Nothing,Nothing,Nothing}(f,grad,hess,nothing,AutoModelingToolkit(),nothing,nothing,nothing,0)
        OptimizationProblem{$iip}(_f,u0,p;lb=lb,ub=ub,kwargs...)
    end
end

function DiffEqBase.OptimizationFunction{iip}(f, ::AutoModelingToolkit, x, p = DiffEqBase.NullParameters();
                              grad=false, hess=false, cons = nothing, cons_j = nothing, cons_h = nothing,
                              num_cons = 0, chunksize = 1, hv = nothing) where iip

    sys = modelingtoolkitize(OptimizationProblem(f,x,p))
    u0map = states(sys) .=> x
    if p == DiffEqBase.NullParameters()
        parammap = DiffEqBase.NullParameters()
    else
        parammap = parameters(sys) .=> p
    end
    OptimizationProblem(sys,u0map,parammap,grad=grad,hess=hess).f
end


text:
module IBMQ

export @ibmq

using YaoCompiler
using IBMQClient
using YaoAPI
using REPL.TerminalMenus
using Crayons.Box
using IBMQClient: AccountInfo, JobInfo, DataInfo, IBMQDevice, DeviceMenu

# RL: qiskit stores it on disk anyway, so I think it's fine for us to
# store it temporarily in a Julia instance just for convenience

"account cache to speed up register creation and queries"
const account_cache = Dict{String,AccountInfo}()

"""
    IBMQReg <: AbstractRegister{1}

A virtual IBM Q register for uploading jobs to IBM Q REST API.
"""
struct IBMQReg <: AbstractRegister{1}
    account::AccountInfo
    device::IBMQDevice

    nshots::Int
    memory_slots::Int
    nqubits::Int

    options::NamedTuple
    # max_credits::Int
    # seed::Int
    # schema_version::String
    # type::String
    # memory::Bool
    # init_qubits::Bool
    # parameter_binds::Vector
    # parametric_pulses
end

function Base.show(io::IO, m::MIME"text/plain", r::IBMQReg)
    println(io, "IBM Q (virtual) register:")
    show(IOContext(io, :indent => 2), m, r.account)
    println(io)
    println(io)
    show(IOContext(io, :indent => 2), m, r.device)
end

function IBMQReg(;
    token::Union{Nothing,String} = nothing,
    device::Union{Nothing,String} = nothing,
    nshots::Int = 1024,
    memory_slots::Int = -1,
    nqubits::Int = -1, # NOTE: we use nqubits here for consistency
    kw...,
    # max_credits::Int = 3,
    # seed::Int = 1,
    # schema_version::String="1.3.0",
    # type::String="QASM",
    # memory::Bool=false,
    # init_qubits::Bool=true,
    # parameter_binds::Vector=[],
    # parametric_pulses=[],
)

    if token === nothing
        if haskey(ENV, "IBMQ_TOKEN")
            token = ENV["IBMQ_TOKEN"]
        else
            buf = Base.getpass("IBM API token (https://quantum-computing.ibm.com/account)")
            token = read(buf, String)
            Base.shred!(buf)
        end
    end

    global account_cache
    if haskey(account_cache, token)
        account = account_cache[token]
    else
        account = AccountInfo(token)
        account_cache[token] = account
    end

    devices = IBMQClient.devices(account)

    if device === nothing
        menu = DeviceMenu(devices, pagesize = 6)
        choice = request("choose a device:", menu)
        if choice < 0
            return
        end
    else
        dev = findfirst(x -> x.name == device, devices)
        dev === nothing || error("device \"$device\" is not available.")
    end

    return IBMQReg(account, devices[choice], nshots, memory_slots, nqubits, NamedTuple(kw))
end

function create_main_qobj(r::IBMQReg, specs::RoutineSpec...)
    target = YaoCompiler.TargetQobjQASM(;
        nshots = r.nshots,
        max_credits = get(r.options, :max_credits, 3),
        seed = get(r.options, :seed, 1),
    )

    experiments = map(specs) do spec
        # TODO: support parametric routine spec, mark parameters as Const in typeinf
        ci, _ = YaoCompiler.code_yao(YaoCompiler.Semantic.main, spec; optimize = true, passes = [:julia])
        return YaoCompiler.codegen(target, ci)
    end

    qobj = IBMQClient.create_qobj(
        r.device,
        experiments...;
        memory_slots = r.memory_slots,
        shots = r.nshots,
        n_qubits = r.nqubits,
        r.options...,
    )

    return qobj
end

function submit(r::IBMQReg, specs::RoutineSpec...)
    qobj = create_main_qobj(r, specs...)
    return IBMQClient.submit(r.account, r.device, qobj)
end

# syntax sugar
# TODO: spawn a task in background to check status
#       so we don't need to do that manually
# function (spec::RoutineSpec)(r::IBMQReg)
#     return submit(r, spec)
# end

# disable YaoLang entry
function YaoCompiler.execute(::typeof(YaoCompiler.Semantic.main), ::IBMQReg, ::RoutineSpec)
    error(
        "IBM Q devices is not fully compatible with YaoLang, please use @ibmq to define and submit jobs with a subset of YaoLang",
    )
end

macro ibmq(ex...)
    return esc(ibmq_m(ex...))
end

function ibmq_m(ex...)
    options = ex[1:end-1]
    kwargs = []
    for opt in options
        opt isa Expr && opt.head === :(=) || error("invalid syntax $opt")
        opt.args[1] in [:nshots, :max_credits, :seed, :schema_version] ||
            error("invalid option: $opt")
        push!(kwargs, Expr(:kw, opt.args[1], opt.args[2]))
    end

    experiment = last(ex)
    experiment isa Expr || error("invalid syntax: $experiment")

    seq = gensym(:seq)
    return quote
        $seq = Any[]
        $(collect_experiments(seq, experiment))
        r = $IBMQReg(; $(kwargs...))
        $submit(r, $seq...)
    end
end

function collect_experiments(seq::Symbol, @nospecialize(ex))
    ex isa Expr || return ex # skip other things
    ex.head == :call && return :(push!($seq, $ex))

    if ex.head === :block
        return Expr(:block, map(x -> collect_experiments(seq, x), ex.args)...)
    elseif ex.head === :for
        return Expr(:for, ex.args[1], collect_experiments(seq, ex.args[2]))
    else
        error("invalid syntax: $ex")
    end
end

end


text:
module Config

using JSON
using Dates
using Mmap
using Logging
using HDF5
using Statistics

# utilities
include("Utils.jl")

# data managers
include("Save.jl")

# Base Class
include("Manager.jl")
export ConfigManager, iterator, total_combinations, add!, parse!,
    get_sweep_params, param_setting_from_id, params_with, param_values,
    get_datafile, get_stepslog, get_expdir, get_logdir,
    log_config, get_subconfig, get_attrs,
    getindex, setindex!, haskey

# Saving/Loading/Parsing data
include("Data.jl")
export save, save!, load, update!, clear!,
     get_run_data, get_data, get_summary, get_best,
    sensitivity, get_max_runs,
    Summary

# Checklist to prevent overwriting results or redundantly re-running
# already-ran experiments
include("Checklist.jl")
export dump_config, mark_config, check_config, unmarked_configs

end # module


text:
#equation of states

#gamma law
function eos_press(rho, eps, gamma)
    const pmin = 1.0e-10
    const pmax = 1.0e20

    press = (gamma - 1.0) .* rho .* eps
    press = clamp(press, pmin, pmax)
    return press
end

#soundspeed
function eos_cs2(rho, eps, gamma)
    prs = (gamma - 1.0) .* rho .* eps
    dpde = (gamma - 1.0) .* rho
    dpdrho = (gamma - 1.0) .* eps
    cs2 = dpdrho .+ dpde .* prs ./ (rho .+ 1.0e-30).^2.0

    cs2 = clamp(cs2, 1.0e-10, 1.0e10)

    #for i = 1:length(cs2)
    #    if cs2[i] < 0.0
    #        println("cs2 = $(cs2[i]) i=$i")
    #        println(eps[i])
            #println(dpdrho[i])
            #println(dpde[i])
            #println(prs[i])
    #    end
    #end

    return cs2
end

max_stars_repo_name: natj/hydro

max_issues_repo_name: natj/hydro

max_forks_repo_name: natj/hydro