text:
module Telegrambot

import HTTP, Test, JSON, UUIDs
export InlineQueryResultArticle
export startBot, getUpdates, sendText

struct InlineQueryResultArticle
    id::String
    title::String
    message::String
end

function startBot(botApi=""; textHandle=Dict(), inlineQueryHandle=Dict(), customHandle=_->nothing)
    #!isempty(textHandle) || error("You need to pass repond function as parameter to startBot")
    # in case people put in what botfather spits
    botApi = botApi[1:3]=="bot" ? botApi : "bot" * botApi
    # to be used to clear msg que
    offset = 0

    msgDict = Dict()
    while true
        msgQuery = getUpdates(botApi,offset)
        if length(msgQuery)==0
            #this repeast every timeout seconds, now 30
            continue
        end
        offset = maximum([ i["update_id"] for i in msgQuery ]) + 1 #update to clear the msg query next loop

        # HTTP respond will have this field if a user @'ed bot in a chat
        # find match cmd to pass to correspond function to handle
        for rawCmd in msgQuery
            # if this is a message
            if haskey(rawCmd, "message")
                msg = rawCmd["message"]
                cmdName=" "
                cmdPara = " "
                try
                    cmdName = string(match(r"/([^@\s]+)", msg["text"])[1]) #match till first @ or space
                catch
                    #= @warn "Not a command start with /" =#
                    continue
                end

                try
                    cmdPara = string(match(r"\s(.*)$", msg["text"])[1]) #match from first space to end
                catch
                    cmdPara = " "
                    #= @warn "Command may got passed empty parameter" =#
                end
                if haskey(textHandle, cmdName)
                    reply = textHandle[cmdName](cmdPara)  #encode for GET purpose
                    reply_id = string(msg["chat"]["id"])  #encode for GET purpose
                    # all space msg is not allower either
                    !isempty(strip(reply)) || (@warn " message must be non-empty, also can't be all spaces";
                                               reply = "Using the command incorrectly or command is bad")
                    sendText(botApi, reply_id, reply)
                else
                    #reply_id = string(msg["chat"]["id"])  #encode for GET purpose
                    #no_cmd_prompt= "The command $cmdName is not found"
                    #sendText(botApi, reply_id, no_cmd_prompt)
                    #= @warn backtrace() =#
                end
            end

            # HTTP respond will have this field if a user started an inline request
            if haskey(rawCmd, "inline_query")
                qry = rawCmd["inline_query"]
                if qry["query"]≠"" #multiple inline query fills up as user type so initially there is an empty string
                    articleList = InlineQueryResultArticle[] #make an array of Articles waiting to be converted
                    for (key, inlineOpt) in inlineQueryHandle
                        articleEntry = InlineQueryResultArticle(string(UUIDs.uuid4()), key, inlineOpt(qry["query"]))
                        push!(articleList, articleEntry)
                    end
                    results = articleList |> ArticleListtoJSON  #encode according to https://core.telegram.org/bots/api#inlinequeryresult
                    query_id = qry["id"]  #encode for GET escape
                    answerInlineQuery(botApi, query_id, results)
                end
            end

            customHandle(rawCmd)
        end
    end
end

# GET request sendig to telegram for text repond
function sendText(botApi, id, text)
    # can't pass empty text, results 400
    text = text |> HTTP.URIs.escapeuri #encode for GET purpose
    id = id |> HTTP.URIs.escapeuri #encode for GET purpose
    tQuery="""chat_id=$id&text=$text"""
    try
        updates = HTTP.request("GET","https://api.telegram.org/$botApi/sendMessage";query="$tQuery")
    catch e
        errmsg = JSON.parse(String(e.response.body))
        @warn "$(errmsg["description"]): $id"
    end
    sleep(0.01)
end

# GET request sendig to telegram for inline respond
function answerInlineQuery(botApi, query_id, results::String)
    results = results |> HTTP.URIs.escapeuri #encode for GET purpose
    query_id = query_id |> HTTP.URIs.escapeuri #encode for GET purpose
    tQuery="""inline_query_id=$query_id&results=$results"""
    updates = HTTP.request("GET","https://api.telegram.org/$botApi/answerInlineQuery";query="$tQuery")
    sleep(0.01)
end

#encode a list of InlineQueryResultArticle according to telegram api JSON format
function ArticleListtoJSON(articles::Array{InlineQueryResultArticle})
    return JSON.json([ Dict(
                            "type"=>"article","id"=>article.id,
                            "title"=>article.title,
                            "input_message_content"=>
                            Dict{String,String}("message_text"=>article.message)
                           )
                      for article in articles])
end

function getUpdates(botApi="", offset=0)
    # this is already a long-poll handled on telegram's side
    # telling telegram how long we want to timeout once
    tQuery="""timeout=30&offset=$offset"""
    updates = JSON.parse(String(HTTP.request("GET","https://api.telegram.org/$botApi/getUpdates";query="$tQuery").body))
    result = updates["result"]
    return result
end

end # module


text:
widened{S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C}) = Flex{S,Q,C}(prevFloat(a.lo), nextFloat(a.hi))

function narrowed{S<:Sculpt,Q<:Qualia,C<:Clay}(a::Flex{S,Q,C})
   lo, hi = nextFloat(a.lo), prevFloat(a.hi)
   if (lo > hi)
       a
   else
       Flex{S,Q,C}(lo,hi)
   end
end

# from "Complete Interval Arithmetic and its Implementation on the Computer" by Ulrich W. Kulisch

function glb{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}, y::Flex{S,Q,C})
    lo = min(x.lo,y.lo)
    hi = min(x.hi,y.hi)
    Flex{S,Q,C}(minmax(lo,hi)...)
end

function lub{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}, y::Flex{S,Q,C})
    lo = max(x.lo,y.lo)
    hi = max(x.hi,y.hi)
    Flex{S,Q,C}(minmax(lo,hi)...)
end

# supremum (join, interval [convex] hull) and an infimum (meet, intersection).
function supremum{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}, y::Flex{S,Q,C})
    lo = min(x.lo,y.lo)
    hi = max(x.hi,y.hi)
    Flex{S,Q,C}(minmax(lo,hi)...)
end

function infimum{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}, y::Flex{S,Q,C})
    lo = max(x.lo,y.lo)
    hi = min(x.hi,y.hi)
    Flex{S,Q,C}(minmax(lo,hi)...)
end

# from "Standardized notation in interval analysis", R. B. Kearfott, et. al.
# width(diameter), radius, midpoint, mignitude, magnitude, deviation, absolutevalue

dia{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = (x.hi - x.lo)
rad{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = (x.hi - x.lo)/2
mid{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = (x.lo + x.hi)/2
mig{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) =
   (signbit(x.lo)==signbit(x.hi)) ? min(abs(x.lo), abs(x.hi)) : 0
mag{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = max(abs(x.lo), abs(x.hi))
dev{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = (abs(x.lo) >= abs(x.hi)) ? x.lo : x.hi
abs{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}) = Flex{S,Q,C}(mig(x),mag(x))
dist{S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C}, y::Flex{S,Q,C}) = max(abs(x.lo-y.lo), abs(x.hi-y.hi))

for fn in (:widened, :narrowed, :opened, :closed, :glb, :lub, :supremum, :infimum,
           :dia, :rad, :mid, :mig, :mag, :dev, :abs, :dist, :value, :isexact, :isinexact,
           :isclosedclosed, :isclosedopen, :isopenclosed, :isopenopen)
    @eval begin
        @vectorize_1arg Flex ($fn)
    end
end

ext: jl


text:
# === MODEL FUNCTIONS === #

num_neurons(model::SeqModel) = model.priors.bkgd_proportions.dim
num_sequence_types(model::SeqModel) = model.priors.seq_type_proportions.dim
num_warp_values(model::SeqModel) = length(model.priors.warp_values)
num_sequence_events(model::SeqModel) = length(model.sequence_events)

"""Helper function to initialize warp probabilities."""
function default_warps(num_warp_values::Int64,
                       τ_max::Float64,
                       warp_variance::Float64)
    @assert num_warp_values >= 1
    @assert τ_max >= 1 "max warp value must be >= 1 by convention."
    @assert warp_variance >= 0

    # Initialize the warp parameters.
    if num_warp_values == 1
        warp_values = ones(1)
        warp_log_proportions = zeros(1)
    else
        # Log-spaced values between 1/τ_max and τ_max
        warp_values = τ_max .^ range(-1, 1, length=num_warp_values)
        # Set a mean-zero Gaussian prior on the log warp values
        warp_log_proportions = -0.5 / warp_variance * range(-1, 1, length=num_warp_values) .^ 2
        warp_log_proportions .-= logsumexp(warp_log_proportions)
    end
    warp_values, warp_log_proportions
end

"""User-Friendly Constructor of SeqModel model."""
function SeqModel(
        # constants
        max_time::Float64,
        max_sequence_length::Float64,

        # warp parameters
        num_warp_values::Int64,
        max_warp::Float64,
        warp_variance::Float64,

        # priors
        seq_event_rate::Float64,
        seq_type_proportions::SymmetricDirichlet,
        seq_event_amplitude::RateGamma,
        neuron_response_proportions::SymmetricDirichlet,
        neuron_response_profile::NormalInvChisq,
        bkgd_amplitude::RateGamma,
        bkgd_proportions::SymmetricDirichlet
    )

    # Number of neurons (N), sequence types (R), and warp values (W).
    N = neuron_response_proportions.dim
    R = seq_type_proportions.dim
    W = num_warp_values

    # Specify discrete distribution over warp values.
    warp_values, warp_log_proportions =
        default_warps(num_warp_values, max_warp, warp_variance)
    @show warp_values, warp_log_proportions

    # Create prior distributions.
    priors = SeqPriors(
        seq_event_rate,
        seq_type_proportions,
        seq_event_amplitude,
        neuron_response_proportions,
        neuron_response_profile,
        bkgd_amplitude,
        bkgd_proportions,
        warp_values,
        warp_log_proportions
    )

    # Initialize global variables.
    globals = sample(priors)

    # Initialize sequence events.
    sequence_events = SeqEventList(R, W)

    # Compute probability of introducing a new cluster.
    α = seq_event_amplitude.α
    β = seq_event_amplitude.β
    λ = seq_event_rate

    # Compute probability of introducing a background spike.
    bkgd_log_prob = (
        log(globals.bkgd_amplitude)
        + log(max_time)
        + log(1 + β)
    )
    new_cluster_log_prob = (
        log(α)
        + log(λ)
        + log(max_time)
        + α * (log(β) - log(1 + β))
    )

    # Package it all together into the model object.
    return SeqModel(
        max_time,
        max_sequence_length,
        priors,
        globals,
        sequence_events,
        new_cluster_log_prob,
        bkgd_log_prob,
        zeros(R),      # _R_buffer
        zeros(R, W),   # _RW_buffer
        Float64[],     # _K_buffer
    )
end

"""
By default, assume no warping.  (Or equivalently, warp_values=[1].)
"""
function SeqModel(
        # constants
        max_time::Float64,
        max_sequence_length::Float64,

        # priors
        seq_event_rate::Float64,
        seq_type_proportions::SymmetricDirichlet,
        seq_event_amplitude::RateGamma,
        neuron_response_proportions::SymmetricDirichlet,
        neuron_response_profile::NormalInvChisq,
        bkgd_amplitude::RateGamma,
        bkgd_proportions::SymmetricDirichlet
    )

    SeqModel(max_time,
          max_sequence_length,
          1, # number of warps
          1.0, # max warp value
          1.0, # warp variance
          seq_event_rate,
          seq_type_proportions,
          seq_event_amplitude,
          neuron_response_proportions,
          neuron_response_profile,
          bkgd_amplitude,
          bkgd_proportions)
end

# ======================================================================== #
# ===                                                                  === #
# === Methods to draw samples from a model (or the prior distribution) === #
# ===                                                                  === #
# ======================================================================== #

function sample(priors::SeqPriors)

    N = priors.neuron_response_proportions.dim
    R = priors.seq_type_proportions.dim

    # Initialize length-R probability vector over sequence types.
    seq_type_log_proportions = log.(rand(priors.seq_type_proportions))

    # Draw N x R matrix holding neuron response amplitudes.
    neuron_response_log_proportions =
        log.(rand(priors.neuron_response_proportions, R))

    # Draw N x R matrix holding neuron response offsets and widths.
    neuron_response_offsets = zeros(N, R)
    neuron_response_widths = zeros(N, R)
    for n = 1:N
        for r = 1:R
            μ, σ2 = rand(priors.neuron_response_profile)
            neuron_response_offsets[n, r] = μ
            neuron_response_widths[n, r] = σ2
        end
    end

    # Draw background rate parameter.
    bkgd_amplitude = rand(priors.bkgd_amplitude)

    # Draw length-N probability vector holding normalized background firing rates.
    bkgd_log_proportions = log.(rand(priors.bkgd_proportions))

    return SeqGlobals(
        seq_type_log_proportions,
        neuron_response_log_proportions,
        neuron_response_offsets,
        neuron_response_widths,
        bkgd_amplitude,
        bkgd_log_proportions
    )
end

function sample(
        model::SeqModel;
        resample_latents::Bool=false,
        resample_globals::Bool=false
    )

    # =========== (i) SAMPLE GLOBAL PARAMS =========== #

    globals =
        resample_globals ? sample(model.priors) : deepcopy(model.globals)

    # =========== (ii) SAMPLE LATENT EVENTS =========== #

    if resample_latents

        K = rand(Poisson(model.priors.seq_event_rate * model.max_time))

        seq_type_dist = Categorical(exp.(globals.seq_type_log_proportions))
        warp_dist = Categorical(exp.(model.priors.warp_log_proportions))
        events = EventSummaryInfo[]

        for k = 1:K
            τ = rand() * model.max_time
            r = rand(seq_type_dist)
            ω = model.priors.warp_values[rand(warp_dist)]
            A = rand(model.priors.seq_event_amplitude)
            push!(events, EventSummaryInfo(k, τ, r, ω, A))
        end

    else
        events = event_list_summary(model)
        # TODO -- need to sample empty latent events too...
    end

    # =========== (iii) SAMPLE SPIKES =========== #

    spikes = Spike[]
    assignments = Int64[]

    # Sample background spikes.
    S_bkgd = rand(Poisson(globals.bkgd_amplitude * model.max_time))
    bkgd_dist = Categorical(exp.(model.globals.bkgd_log_proportions))

    n_bkgd = rand(bkgd_dist, S_bkgd)
    t_bkgd = rand(S_bkgd) * model.max_time

    for (n, t) in zip(n_bkgd, t_bkgd)
        push!(spikes, Spike(n, t))
        push!(assignments, -1)
    end

    # Compute neuron probabilities.
    neuron_rel_amps =
        exp.(globals.neuron_response_log_proportions)

    # Sample sequence-evoked spikes.
    for e in events

        # Num spikes evoked by latent event.
        S = rand(Poisson(e.amplitude))

        # Neuron response proportions.
        nrn_dist = Categorical(neuron_rel_amps[:, e.seq_type])

        # Sample neuron, then spike time, for each spike.
        for n in rand(nrn_dist, S)

            # Sample spike time.
            μ = globals.neuron_response_offsets[n, e.seq_type]
            σ = sqrt(globals.neuron_response_widths[n, e.seq_type])
            t = e.seq_warp * (σ * randn() + μ) + e.timestamp

            # Exclude spikes outside of time window.
            if (0 < t < model.max_time)
                push!(spikes, Spike(n, t))
                push!(assignments, e.assignment_id)
            end

        end
    end

    return spikes, assignments, events

end

"""
Return firing (num_neurons x num_timebins) firing rate matrix.
"""
function firing_rates(
        globals::SeqGlobals,
        events::Vector{EventSummaryInfo},
        timebins::Vector{Float64};
        resample_latents::Bool=false,
        resample_globals::Bool=false
    )

    # Discretized time bins
    num_neurons = length(globals.bkgd_log_proportions)

    # Compute background intensities for each neuron.
    bkgd_firing_rates = globals.bkgd_amplitude .* exp.(globals.bkgd_log_proportions)

    # Initialize matrix of firing rates with background rates.
    F = repeat(bkgd_firing_rates, outer=(1, length(timebins)))

    # Compute neuron probabilities.
    neuron_rel_amps =
        exp.(globals.neuron_response_log_proportions)

    # Sample sequence-evoked spikes.
    for e in events
        for n = 1:num_neurons
            μ = e.timestamp + globals.neuron_response_offsets[n, e.seq_type]
            σ = sqrt(globals.neuron_response_widths[n, e.seq_type])
            F[n, :] += e.amplitude * neuron_rel_amps[n] * normpdf.(μ, σ, timebins)
        end
    end

    return F

end

priors(model::SeqModel) = model.priors

function set_priors!(model::SeqModel, priors::SeqPriors)
    model.priors = priors
end

function set_globals!(model::SeqModel, globals::SeqGlobals)
    model.globals = globals
end

function set_new_cluster_log_prob!(model::SeqModel, prob::Float64)
    model.new_cluster_log_prob = prob
end

get_max_time(model::SeqModel) = model.max_time

get_bkgd_amplitude(model::SeqModel) = model.globals.bkgd_amplitude

function set_bkgd_log_prob!(model::SeqModel, prob::Float64)
    model.bkgd_log_prob = prob
end


text:
@testset "scaling" begin

    @testset "IdentityScaling" begin
        scaling = IdentityScaling()
        @test scaling isa Transform
        @test cardinality(scaling) == OneToOne()

        @testset "Arguments do nothing" begin
            @test IdentityScaling(123) == IdentityScaling()
            @test IdentityScaling([1, 2, 3]) == IdentityScaling()
        end
    end

    @testset "MeanStdScaling" begin
        @testset "Constructor" begin
            M = [0.0 -0.5 0.5; 0.0 1.0 2.0]
            nt = (a = [0.0, -0.5, 0.5], b = [1.0, 0.0, 2.0])

            @testset "simple" for x in (M, nt)
                x_copy = deepcopy(x)
                scaling = MeanStdScaling(x)
                @test cardinality(scaling) == OneToOne()
                @test scaling isa Transform
                @test x == x_copy  # data is not mutated
                # constructor uses all data by default
                @test scaling.μ == 0.5
                @test scaling.σ ≈ 0.89443 atol=1e-5
            end

            @testset "use certain slices to compute statistics" begin
                @testset "Array" begin
                    scaling = MeanStdScaling(M; dims=1, inds=[2])
                    @test scaling.μ == 1.0
                    @test scaling.σ == 1.0
                end

                @testset "Table" begin
                    scaling = MeanStdScaling(nt; cols=:a)
                    @test scaling.μ == 0.0
                    @test scaling.σ == 0.5
                end
            end
        end

        @testset "Re-apply" begin
            M = [0.0 -0.5 0.5; 0.0 1.0 2.0]
            scaling = MeanStdScaling(M; dims=2)
            new_M = [1.0 -2.0 -1.0; 0.5 0.0 0.5]
            @test M !== new_M
            # Expect scaling parameters to be fixed to the first data applied to
            expected_reapply = [0.559017 -2.79508 -1.67705; 0.0 -0.55901 0.0]
            @test FeatureTransforms.apply(new_M, scaling; dims=2) ≈ expected_reapply atol=1e-5
        end

        @testset "Inverse" begin
            M = [0.0 -0.5 0.5; 0.0 1.0 2.0]
            M_expected = [-0.559017 -1.11803 0.0; -0.559017 0.559017 1.67705]
            scaling = MeanStdScaling(M)
            transformed = FeatureTransforms.apply(M, scaling)

            @test transformed ≈ M_expected atol=1e-5
            @test FeatureTransforms.apply(transformed, scaling; inverse=true) ≈ M atol=1e-5
        end

        @testset "Zero std" begin
            x = ones(Float64, 3)
            expected = zeros(Float64, 3)

            scaling = MeanStdScaling(x)

            @test FeatureTransforms.apply(x, scaling) == expected  # default `eps`
            @test FeatureTransforms.apply(x, scaling; eps=1) == expected
            @test all(isnan.(FeatureTransforms.apply(x, scaling; eps=0)))  # 0/0
        end
    end
end


text:
using Distributions

# config
Ks = [3, 10, 30, 100, 1000, 10000]
n = 10^5

function run_sample(sampler, n::Int)
	for i = 1:n
		rand(sampler)
	end
end

macro bench_sampler(construct_expr)
	quote
		s = $(esc(construct_expr)) # construct sampler
		tname = typeof(s)

		# warm up
		run_sample(s, 100)

		et = @elapsed run_sample(s, n)

		mps = n * 1.0e-6 / et
		@printf("  %-28s :  elapsed = %12.6f sec  |  %8.4f MPS\n",
			tname, et, mps)
	end
end

for K in Ks
	println("K = $K")
	println("-----------------------")

	p = fill(1/K, K)
	pu = fill(uint64(1), K)

	@bench_sampler DiscreteUniform(1, K)   # samples generated by rand(1:K)
	@bench_sampler Categorical(p)
	# @bench_sampler Distributions.DiscreteDistributionTable(p)
	# @bench_sampler Distributions.huffman(1:K, pu)
	@bench_sampler Distributions.AliasTable(p)

	println()
end


text:
@testset "Inner products" begin
    rmax = 1.0
    ρ = 0.1
    N = ceil(Int, rmax/ρ)
    k = 7

    @testset "$(grid_type)" for (grid_type,MS,MatElMetric) = [
        (:fd, Diagonal, ρ*I),
        (:sfd, Diagonal, ρ*I),
        (:implicit_fd, Diagonal, ρ*I),
        (:fedvr, Diagonal, I),
        (:bsplines, BandedMatrix, I)
    ]
        R = if grid_type == :fd
            FiniteDifferences(N, ρ)
        elseif grid_type == :sfd
            StaggeredFiniteDifferences(N, ρ)
        elseif grid_type == :implicit_fd
            ImplicitFiniteDifferences(N, ρ)
        elseif grid_type == :fedvr
            t = range(0, stop=rmax, length=3)
            FEDVR(t, k)
        elseif grid_type == :bsplines
            t = LinearKnotSet(k, 0, rmax, 3)
            BSpline(t)
        end

        r = axes(R,1)

        c = R \ sin.(2π*r)
        f = R*c
        d = R \ cos.(2π*r)
        g = R*d

        F = R*hcat(c,d)
        FF = F'F
        @test size(FF) == (2,2)
        @test isapprox(FF, 0.5I, rtol=1e-3)

        @test f isa CompactBases.FuncVector
        @test isapprox(norm(f), 1/√2, rtol=1e-3)
        f′ = deepcopy(f)
        normalize!(f′)
        @test isapprox(norm(f′), 1, rtol=1e-3)

        S = R'R

        @test apply(*, c', R', R, c) isa Real
        @test isapprox(apply(*, c', R', R, c), 0.5, rtol=1e-3)
        @test dot(c, S, c) isa Real
        @test isapprox(dot(c, S, c), 0.5, rtol=1e-3)
        @test dot(c, S, d) isa Real
        @test isapprox(dot(c, S, d), 0.0, atol=1e-14)
        @test f'f isa Real
        @test isapprox(f'f, 0.5, rtol=1e-3)
        @test f'g isa Real
        @test isapprox(f'g, 0.0, atol=1e-14)
        @test g'g isa Real
        @test isapprox(g'g, 0.5, rtol=1e-3)

        @test metric(R) == S
        @test metric_shape(R) <: MS
        @test matrix_element_metric(R) == MatElMetric
    end

    @warn "Need to test inner products with restricted bases as well"
end


text:
# This file was generated, do not modify it. # hide
rms(y_wind[train], MLJ.predict(m_linReg, rows=train))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
## ------------------------------------------------------------------
function _setproperty!(callinfo::LineNumberNode, ex)

    # obj
    obj = _check_currobj()

    # check
    _check_call_inbook(callinfo, obj)
    _check_bang_call(callinfo, obj)

    # collect
    meta = _extract_expr_pairs(ex)

    # quote
    exout = :()
    for (key, val) in meta
        exout = quote
            $(exout)
            ResearchBooks.setproperty!($(obj), $(Meta.quot(key)), $(esc(val)))
        end
    end
    return quote
        $(exout); $(esc(obj))
    end
end

function _invoke_setter(setter::Function, callinfo::LineNumberNode, ex...)
    obj = _check_currobj()
    _check_call_inbook(callinfo, obj)
    _check_bang_call(callinfo, obj)

    setter(obj, ex...)

    return obj
end

## ------------------------------------------------------------------
# Meta setters

macro setproperty!(ex...)
    _setproperty!(__source__, ex)
end

macro set_title!(ex::String)
    _invoke_setter(set_title!, __source__, ex)
end

macro set_doi!(ex::String)
    _invoke_setter(set_doi!, __source__, ex)
end

macro set_bibkey!(ex::String)
    _invoke_setter(set_bibkey!, __source__, ex)
end

macro set_abstract!(ex::String)
    _invoke_setter(set_abstract!, __source__, ex)
end

macro set_text!(ex::String)
    _invoke_setter(set_text!, __source__, ex)
end

macro set_year!(ex::String)
    _invoke_setter(set_year!, __source__, string(ex))
end

macro set_author!(ex::Union{String, Symbol}...)
    _invoke_setter(set_author!, __source__, string.(ex)...)
end

macro set_ctime!(ex::String)
    _invoke_setter(set_ctime!, __source__, ex)
end

macro set_similars!(ex::Union{String, Symbol}...)
    _invoke_setter(set_similars!, __source__, string.(ex)...)
end

macro set_antagonists!(ex::Union{String, Symbol}...)
    _invoke_setter(set_antagonists!, __source__, string.(ex)...)
end

macro set_description!(ex::String)
    _invoke_setter(set_description!, __source__, ex)
end

macro set_supertag!(ex::String)
    _invoke_setter(set_supertag!, __source__, ex)
end

## ------------------------------------------------------------------
# Meta adders

macro add_tag!(ex::Union{String, Symbol}...)
    _invoke_setter(add_tag!, __source__, string.(ex)...)
end

macro add_todo!(ex::Union{String, Symbol}...)
    _invoke_setter(add_todo!, __source__, string.(ex)...)
end


text:
@testset "Search" begin

    @testset "is_key" begin
        @test AAindex.is_key("CORJ870107")
        @test !AAindex.is_key("701078JROC")
        @test !AAindex.is_key("COJ80107")
    end

    @testset "search" begin
        id, description, index_type = only(search("PHAT"))

        @test id == "NGPC000101"
        @test index_type === AAindex.AMatrix

        result = aaindex_by_id(id)

        @test result isa AAindex.AMatrix
        @test result.metadata.title == "PHAT: a transmembrane-specific substitution matrix"
    end

end


text:
module autodiff

using PyCall
include("utils.jl")

export Variable,
    Matrix,
    Constant,
    Empty,
    add,
    mul,
    sub,
    add_byconst,
    mul_byconst,
    sub_byconst,
    ones,
    oneslike,
    zeroslike,
    negative,
    power,
    einsum,
    norm,
    identity,
    tensorinv,
    scalar,
    transpose,
    sum,
    tensordot,
    Executor,
    jacobians,
    jvps,
    jtjvps,
    transposed_vjps,
    gradients,
    hvp,
    hessian

const ad = PyNULL()

function __init__()
    copy!(ad, pyimport_conda("autohoot.autodiff", "autohoot"))
end

@func_loader ad Variable Matrix Constant Empty add mul sub add_byconst mul_byconst sub_byconst ones oneslike zeroslike negative power einsum norm identity tensorinv scalar transpose sum tensordot
@func_loader ad Executor jacobians jvps jtjvps transposed_vjps gradients hvp hessian

end


text:
function task_to_array(task)
    arr = []
    for e=task
        push!(arr, e)
    end
    arr
end

# These keys only give access to the JULIA_TESTING table. They are included on purpose.
# See the bottom of this file for more information.
aws = aws_config(region="eu-west-1")

# could also be done using a julia Type or Immutable if you prefer
const table = dynamo_table(Dict, "JULIA_TESTING", :id, :order)

# unique key for testing purposes
const id_key = random_key()

# let's load some data
for i=1:26:260
    # the POST api takes a max of 25 items at once, DynamoDB.jl multiplexes client-side
    batch_put_item(aws, table,
                   Dict("id" => id_key, "order" => i + 0),
                   Dict("id" => id_key, "order" => i + 1),
                   Dict("id" => id_key, "order" => i + 2),
                   Dict("id" => id_key, "order" => i + 3),
                   Dict("id" => id_key, "order" => i + 4),
                   Dict("id" => id_key, "order" => i + 5),
                   Dict("id" => id_key, "order" => i + 6),
                   Dict("id" => id_key, "order" => i + 7),
                   Dict("id" => id_key, "order" => i + 8),
                   Dict("id" => id_key, "order" => i + 9),
                   Dict("id" => id_key, "order" => i + 10),
                   Dict("id" => id_key, "order" => i + 11),
                   Dict("id" => id_key, "order" => i + 12),
                   Dict("id" => id_key, "order" => i + 13),
                   Dict("id" => id_key, "order" => i + 14),
                   Dict("id" => id_key, "order" => i + 15),
                   Dict("id" => id_key, "order" => i + 16),
                   Dict("id" => id_key, "order" => i + 17),
                   Dict("id" => id_key, "order" => i + 18),
                   Dict("id" => id_key, "order" => i + 19),
                   Dict("id" => id_key, "order" => i + 20),
                   Dict("id" => id_key, "order" => i + 21),
                   Dict("id" => id_key, "order" => i + 22),
                   Dict("id" => id_key, "order" => i + 23),
                   Dict("id" => id_key, "order" => i + 24),
                   Dict("id" => id_key, "order" => i + 25))
end

itr = query(aws, table, id_key, between(attr("order"), 14, 16); consistant_read=true)
items = task_to_array(itr)
@test length(items) == 3 # 14, 15, and 16
items[2]["order"] = 200
items[2]["some_new_field"] = "another value we might want"
items[2]["sub_doc"] = Dict("a" => 1, "b" => 2)

put_item(aws, table, items[2]) # writes the value as a *new* item.

items[2]["some_new_field"] = "secondary value"
@test put_item(aws, table, items[2]; return_old=true)["some_new_field"] == "another value we might want"

itr = query(aws, table, id_key, between(attr("order"), 14, 16); consistant_read=true)
items = task_to_array(itr)
@test length(items) == 3

item = get_item(aws, table, id_key, 200)
@test item["some_new_field"] == "secondary value"
@test item["sub_doc"] == Dict("a" => 1, "b" => 2)

update_item(aws, table, id_key, 200, assign(attr("sub_doc", "c"), 3))

item = get_item(aws, table, id_key, 200)
@test item["some_new_field"] == "secondary value"
@test item["sub_doc"] == Dict("a" => 1, "b" => 2, "c" => 3)

delete_item(aws, table, id_key, 200)
@test get_item(aws, table, id_key, 200) == nothing

items = batch_get_item(aws, table, (id_key, 1), (id_key, 2), (id_key, 3), (id_key, 200))
@assert length(items) == 3 # since 200 was deleted

items = batch_get_item(aws, table, (id_key, 1), (id_key, 2), (id_key, 3), (id_key, 200);
                       only_returning=[attr("id")])
@assert length(items) == 3 # since 200 was deleted

# look at some random items in the table. without a limit this goes on forever
# ... hopefully 6k is enough to exercise the rate limiting and retrying code :)
for e=scan(aws, table; limit=60_000)
    @assert e["id"] != nothing
    @assert e["order"] != nothing
end

# cleanup. the projection asks DynamoDB to only return the order attribute to us.
for e = query(aws, table, id_key; projection=[attr("order")])
    delete_item(aws, table, id_key, e["order"])
end

# Pssst... wanna set something like this up for yourself?
# Here's the policy file:

# {
#     "Version": "2012-10-17",
#     "Statement": [
#         {
#             "Effect": "Allow",
#             "Action": [
#                 "dynamodb:GetItem",
#                 "dynamodb:PutItem",
#                 "dynamodb:UpdateItem",
#                 "dynamodb:DeleteItem",
#                 "dynamodb:BatchGetItem",
#                 "dynamodb:BatchWriteItem",
#                 "dynamodb:Query",
#                 "dynamodb:Scan"
#             ],
#             "Resource": [
#                 "arn:aws:dynamodb:us-east-1:300203084218:table/JULIA_TESTING",
#                 "arn:aws:dynamodb:us-east-1:300203084218:table/JULIA_TESTING/index/*"
#             ]
#         }
#     ]
# }

# You can obviously add operations, and tables to suit your needs :)
# ... also note you can also provide row-level permissions using
# https://docs.aws.amazon.com/amazondynamodb/latest/developerguide/FGAC_DDB.Examples.html


text:
# ps3000_get_times_and_values.jl

export ps3000_get_times_and_values!, TimesAndValuesBuffer

immutable TimesAndValuesBuffer
  times  :: Array{Int32,1}
  a      :: Array{Int16,1}
  b      :: Array{Int16,1}
  c      :: Array{Int16,1}
  d      :: Array{Int16,1}
  TimesAndValuesBuffer(length) =
    new(Array(Int32,length), Array(Int16,length), Array(Int16,length),
    Array(Int16,length), Array(Int16,length))
end

function ps3000_get_times_and_values!(handle       :: Int16,
                                     buffer       :: TimesAndValuesBuffer,
                                     time_units,
                                     no_of_values)
  overflow = Ref{Int16}()
  number_of_data_values_per_channel = ccall((:ps3000_get_times_and_values, ps3000driver),Int32,
    (Int16,
     Ptr{Int32}, Ptr{Int16}, Ptr{Int16}, Ptr{Int16}, Ptr{Int16},
     Ref{Int16}, Int16, Int32),
    handle,
    buffer.times, buffer.a, buffer.b, buffer.c, buffer.d,
    overflow, time_units, no_of_values)
  if number_of_data_values_per_channel == 0
    error("ps3000 Error: one or more of the parameters are out of range or the times will overflow with the time_units requested")
  end
  return (number_of_data_values_per_channel, overflow[])
end


text:
export rrefx

"""
`rref(A)`  returns the row-reduced echelon form of the exact
matrix `A`.
"""
function rrefx(A::AbstractMatrix{T}) where {T<:IntegerX}
    AA = big.(A // 1)
    return rrefx(AA)
end

function rrefx(A::AbstractMatrix{T}) where {T}
    AA = copy(A)
    rrefx!(AA)
    return AA
end

function _pivot!(A::AbstractMatrix{T}, i::Int, j::Int) where {T}
    r, c = size(A)
    s = A[i, j]
    row_scale!(A, i, invx(s))

    for k = 1:r
        if k != i
            s = -A[k, j]
            row_add_mult!(A, i, s, k)
        end
    end
end

function rrefx!(A::AbstractMatrix{T}) where {T}
    rn, nc = size(A)

    r = 1
    for j = 1:nc
        # see if column j's first one is at level r or below
        col = A[r:end, j]
        if all(col .== 0)
            continue
        end
        k = findfirst(col .!= 0)
        row_swap!(A, r, k + r - 1)
        _pivot!(A, r, j)
        r += 1
    end

end


text:
export SampleAvg,
       CachedSampleAvg,
       CachedSum,
       discount_rewards,
       discount_rewards!

import StatsBase: countmap

"extend the countmap in StatsBase to support general iterator"
function countmap(iter)
    res = Dict{eltype(iter),Int}()
    for x in iter
        if haskey(res, x)
            res[x] += 1
        else
            res[x] = 1
        end
    end
    res
end

Base.@kwdef mutable struct SampleAvg
    t::Int = 0
    avg::Float64 = 0.0
end

function (s::SampleAvg)(x)
    s.t += 1
    s.avg += (x - s.avg) / s.t
    s.avg
end

Base.@kwdef struct CachedSampleAvg{T}
    cache::Dict{T,SampleAvg} = Dict{T,SampleAvg}()
end

function (c::CachedSampleAvg)(k, x)
    if !haskey(c.cache, k)
        c.cache[k] = SampleAvg()
    end
    c.cache[k](x)
end

Base.Base.@kwdef struct CachedSum{T}
    cache::Dict{T,Float64} = Dict{T,Float64}()
end

function (c::CachedSum)(k, x)
    c.cache[k] = get!(c.cache, k, 0.0) + x
end

function discount_rewards!(new_rewards, rewards, γ)
    new_rewards[end] = rewards[end]
    for i = (length(rewards) - 1):-1:1
        new_rewards[i] = rewards[i] + new_rewards[i+1] * γ
    end
    new_rewards
end

discount_rewards(rewards, γ) = discount_rewards!(similar(rewards), rewards, γ)

discount_rewards_reduced(rewards, γ) = foldr((r, g) -> r + γ * g, rewards)


text:
cd("LAP") do
    print(pwd())
    run(`make`)
end

cd("UMFLP") do
    print(pwd())
    run(`make`)
end

cd("Combo") do
    run(`make`)
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
#=
VT stability algorithm
=#
struct VTStabilityBuffer{B}
    solver_buffer::B
end

@inline function Base.getproperty(b::VTStabilityBuffer, p::Symbol)
    sol_buf = getfield(b, :solver_buffer)
    return getproperty(sol_buf, p)
end

@inline function Base.getindex(b::VTStabilityBuffer, p::Symbol)
    sol_buf = getfield(b, :solver_buffer)
    return getproperty(sol_buf, p)
end

@inline function Base.propertynames(b::VTStabilityBuffer)
    sol_buf = getfield(b, :solver_buffer)
    return propertynames(sol_buf)
end

function VTStabilityBuffer(mix::BrusilovskyEoSMixture)
    thermo_buf = BrusilovskyThermoBuffer(mix)
    loga_parent = similar(thermo_buf.vec1)
    loga_test = similar(loga_parent)
    p_sat = similar(loga_parent)
    bi = similar(loga_parent)
    nmol_test = similar(loga_parent)
    jacobian = similar(thermo_buf.matr)

    solver_buffer = (
        ;
        loga_parent = loga_parent,
        loga_test = loga_test,
        p_sat = p_sat,
        bi = bi,
        nmol_test = nmol_test,
        jacobian = jacobian,
        thermo_buf = thermo_buf,
    )
    return VTStabilityBuffer(solver_buffer)
end

function VTStabilityBuffer(
    thermo_buf::BrusilovskyThermoBuffer
)
    loga_parent = similar(thermo_buf.vec1)
    loga_test = similar(loga_parent)
    p_sat = similar(loga_parent)
    bi = similar(loga_parent)
    nmol_test = similar(loga_parent)
    jacobian = similar(thermo_buf.matr)

    solver_buffer = (
        ;
        loga_parent = loga_parent,
        loga_test = loga_test,
        p_sat = p_sat,
        bi = bi,
        nmol_test = nmol_test,
        jacobian = jacobian,
        thermo_buf = thermo_buf,
    )
    return VTStabilityBuffer(solver_buffer)
end

"""
    vt_stability_buffer(mix)

Create a buffer for intermediate calculations of one-phase stability of a mixture.

For more info see ?vt_stability(mixture).

See also: [`vt_stability`](@ref), [`thermo_buffer`](@ref)
"""
vt_stability_buffer(x) = VTStabilityBuffer(x)

struct VTStabilityResult{T}
    # converged::Bool  # for future
    isstable::Bool
    energy_density::T
    concentration::Vector{T}

    function VTStabilityResult{T}(isstable::Bool, energy_density, concentration) where {T}
        return new{T}(isstable, energy_density, copy(concentration))
    end
end

"""
Оптимизационная процедура для алгоритма проверки стабильности.
Минимизирует D(η').

Возвращает минимум функции D.
Обновляет аргумент функции D - концентрацию η_.

-> η_, D_thresh, Dfunc!, Dgrad!, isvalid_η, buffer
<- D_min (NaN if not converged)
"""
function vt_stability_optim_try!(
    optmethod,
    η_::AbstractVector,  # вектор концентраций, должен быть инициализирован начальными значениями
    grad_::AbstractVector,
    Dfunc!::Function,  # функция D вида D!(η', ∇D_) -> D::Number
    maxstep::Function,
)
    is_converged = false

    try
        result = Downhill.optimize!(Dfunc!, optmethod, η_;
            gtol=1e-3,
            maxiter=1000,
            constrain_step=maxstep,
            reset=false,
        )
        is_converged = result.converged
    catch e
        @warn e
        is_converged = false
    end
    D_min, _ = Dfunc!(optmethod.x, grad_)
    return D_min
end

function vt_stability(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume,
    RT;
    vtstab_buf::VTStabilityBuffer = VTStabilityBuffer(mix),
) where {T}

    nc = length(mix)

    loga_parent = vtstab_buf.loga_parent
    loga_test = vtstab_buf.loga_test
    jacobian = vtstab_buf.jacobian
    thermo_buf = vtstab_buf.thermo_buf
    bi = vtstab_buf.bi
    p_sat = vtstab_buf.p_sat
    nmol_test = vtstab_buf.nmol_test

    log_c_activity!(loga_parent, mix, nmol, volume, RT; buf = thermo_buf)
    loga_parent .+= log.(nmol ./ volume)
    p_parent = pressure(mix, nmol, volume, RT; buf = thermo_buf)

    thresh = -1e-5

    comp = components(mix)
    map!(subst -> subst.b, bi, comp)
    p_sat .= wilson_saturation_pressure.(comp, RT)
    function stabilitytest!(nmol_test, grad)
        log_c_activity!(grad, mix, nmol_test, 1, RT; buf = thermo_buf)
        p_test = pressure(mix, nmol_test, 1, RT; buf = thermo_buf)
        grad .+= log.(nmol_test)
        grad .-= loga_parent
        D = dot(grad, nmol_test) - (p_test - p_parent) / RT
        return D, grad
    end

    function stability_step(nmol_test, dir)
        αm = convert(eltype(nmol_test), Inf)
        @inbounds for i in eachindex(nmol_test, dir)
            α = -nmol_test[i] / dir[i]
            if 0 < α < αm
                αm = α
            end
        end
        α = (1 - nmol_test' * bi) / (dir' * bi)  # Σ(xᵢ+αdᵢ)bᵢ < 1, covolume constrain
        if 0 < α < αm
            αm = α
        end
        return αm
    end

    function D_min(nmol_test, optmethod)
        _, jacobian = log_c_activity_wj!(
            loga_test,
            jacobian,
            mix,
            nmol_test,
            1,
            RT;
            buf = thermo_buf,
        )
        for i in 1:nc
            jacobian[i,i] += 1 / nmol_test[i]
        end
        Downhill.reset!(optmethod, nmol_test, jacobian)
        return vt_stability_optim_try!(
            optmethod,
            nmol_test,
            loga_test,
            stabilitytest!,
            stability_step,
        )
    end

    optmethod = Downhill.CholBFGS(nmol)
    results = Vector{VTStabilityResult{T}}(undef, 4)

    # Initial - gas
    p_init = dot(p_sat, nmol) / sum(nmol)
    nmol_test .= nmol .* p_sat ./ p_init

    # Test - gas
    z_gg = compressibility(mix, nmol_test, p_init, RT, 'g')
    nmol_test .*= p_init / (z_gg * RT * sum(nmol_test))

    D_gg = D_min(nmol_test, optmethod)
    results[1] = VTStabilityResult{T}(D_gg ≥ thresh, D_gg, optmethod.x)

    # Test - liquid
    z_gl = compressibility(mix, nmol_test, p_init, RT, 'l')
    nmol_test .*= z_gg / z_gl

    D_gl = D_min(nmol_test, optmethod)
    results[2] = VTStabilityResult{T}(D_gl ≥ thresh, D_gl, optmethod.x)

    # Initial - liquid
    nmol_test .= nmol ./ p_sat ./ sum(nmol[i] / p_sat[i] for i in 1:nc)
    p_init = dot(p_sat, nmol_test)

    # Test - gas

    z_lg = compressibility(mix, nmol_test, p_init, RT, 'g')
    nmol_test .*= p_init / (z_lg * RT * sum(nmol_test))

    D_lg = D_min(nmol_test, optmethod)
    results[3] = VTStabilityResult{T}(D_lg ≥ thresh, D_lg, optmethod.x)

    # Test - liquid
    z_ll = compressibility(mix, nmol_test, p_init, RT, 'l')
    nmol_test .*= z_lg / z_ll

    D_ll = D_min(nmol_test, optmethod)
    results[4] = VTStabilityResult{T}(D_ll ≥ thresh, D_ll, optmethod.x)

    if isnan(D_gg) && isnan(D_gl) && isnan(D_lg) && isnan(D_ll)  # все 4 попытки провалились
        error("VTStability: all tries have failed")
    end

    isstable = results[1].isstable && results[2].isstable && results[3].isstable && results[4].isstable

    return isstable, results
end


text:
export PutBlock, put, Swap, swap, PSwap, pswap

"""
    PutBlock <: AbstractContainer

Type for putting a block at given locations.
"""
struct PutBlock{N,C,GT<:AbstractBlock} <: AbstractContainer{GT,N}
    content::GT
    locs::NTuple{C,Int}

    function PutBlock{N}(block::GT, locs::NTuple{C,Int}) where {N,M,C,GT<:AbstractBlock{M}}
        @assert_locs_safe N locs
        @assert nqubits(block) == C "number of locations doesn't match the size of block"
        return new{N,C,GT}(block, locs)
    end
end

"""
    put(total::Int, pair)

Create a [`PutBlock`](@ref) with total number of active qubits, and a pair of
location and block to put on.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> put(4, 1=>X)
nqubits: 4
put on (1)
└─ X gate
```

If you want to put a multi-qubit gate on specific locations, you need to write down all possible locations.

```jldoctest; setup=:(using YaoBlocks)
julia> put(4, (1, 3)=>kron(X, Y))
nqubits: 4
put on (1, 3)
└─ kron
   ├─ 1=>X gate
   └─ 2=>Y gate
```

The outter locations creates a scope which make it seems to be a contiguous two qubits for the block inside `PutBlock`.

!!! tips
    It is better to use [`concentrate`](@ref) instead of `put` for large blocks, since put will use the matrix of its contents
    directly instead of making use of what's in it. `put` is more efficient for small blocks.
"""
put(total::Int, pa::Pair{NTuple{M,Int},<:AbstractBlock}) where {M} =
    PutBlock{total}(pa.second, pa.first)
put(total::Int, pa::Pair{Int,<:AbstractBlock}) = PutBlock{total}(pa.second, (pa.first,))
put(total::Int, pa::Pair{<:Any,<:AbstractBlock}) = PutBlock{total}(pa.second, Tuple(pa.first))

"""
    put(pair) -> f(n)

Lazy curried version of [`put`](@ref).

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> put(1=>X)
(n -> put(n, 1 => X gate))
```
"""
put(pa::Pair) = @λ(n -> put(n, pa))

occupied_locs(x::PutBlock) = map(i -> x.locs[i], x.content |> occupied_locs)
chsubblocks(x::PutBlock{N}, b::AbstractBlock) where {N} = PutBlock{N}(b, x.locs)
PropertyTrait(::PutBlock) = PreserveAll()
cache_key(pb::PutBlock) = cache_key(pb.content)

mat(::Type{T}, pb::PutBlock{N,1}) where {T,N} = u1mat(N, mat(T, pb.content), pb.locs...)
mat(::Type{T}, pb::PutBlock{N,C}) where {T,N,C} = unmat(N, mat(T, pb.content), pb.locs)

function apply!(r::AbstractRegister, pb::PutBlock{N}) where {N}
    _check_size(r, pb)
    instruct!(r, mat_matchreg(r, pb.content), pb.locs)
    return r
end

# specialization
for G in [:X, :Y, :Z, :T, :S, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    @eval function apply!(r::AbstractRegister, pb::PutBlock{N,C,<:$GT}) where {N,C}
        _check_size(r, pb)
        instruct!(r, Val($(QuoteNode(G))), pb.locs)
        return r
    end
end

Base.adjoint(x::PutBlock{N}) where {N} = PutBlock{N}(adjoint(content(x)), x.locs)
Base.copy(x::PutBlock{N}) where {N} = PutBlock{N}(x.content, x.locs)
function Base.:(==)(lhs::PutBlock{N,C,GT}, rhs::PutBlock{N,C,GT}) where {N,C,GT}
    return (lhs.content == rhs.content) && (lhs.locs == rhs.locs)
end

function YaoBase.iscommute(x::PutBlock{N}, y::PutBlock{N}) where {N}
    if x.locs == y.locs
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end

const Swap{N} = PutBlock{N,2,G} where {G<:ConstGate.SWAPGate}
const PSwap{N,T} = PutBlock{N,2,RotationGate{2,T,G}} where {G<:ConstGate.SWAPGate}
Swap{N}(locs::Tuple{Int,Int}) where {N} = PutBlock{N}(ConstGate.SWAPGate(), locs)
PSwap{N}(locs::Tuple{Int,Int}, θ::Real) where {N} = PutBlock{N}(rot(ConstGate.SWAPGate(), θ), locs)

"""
    swap(n, loc1, loc2)

Create a `n`-qubit [`Swap`](@ref) gate which swap `loc1` and `loc2`.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> swap(4, 1, 2)
nqubits: 4
put on (1, 2)
└─ SWAP gate
```
"""
swap(n::Int, loc1::Int, loc2::Int) = Swap{n}((loc1, loc2))

"""
    swap(loc1, loc2) -> f(n)

Create a lambda that takes the total number of active qubits as input. Lazy curried
version of `swap(n, loc1, loc2)`. See also [`Swap`](@ref).

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> swap(1, 2)
(n -> swap(n, 1, 2))
```
"""
swap(loc1::Int, loc2::Int) = @λ(n -> swap(n, loc1, loc2))

function mat(::Type{T}, g::Swap{N}) where {T,N}
    mask = bmask(g.locs[1], g.locs[2])
    orders = map(b -> swapbits(b, mask) + 1, basis(N))
    return PermMatrix(orders, ones(T, 1 << N))
end

apply!(r::AbstractRegister, g::Swap) = (instruct!(r, Val(:SWAP), g.locs); r)
occupied_locs(g::Swap) = g.locs

"""
    pswap(n::Int, i::Int, j::Int, α::Real)
    pswap(i::Int, j::Int, α::Real) -> f(n)

parametrized swap gate.
"""
pswap(n::Int, i::Int, j::Int, α::Real) = PSwap{n}((i, j), α)
pswap(i::Int, j::Int, α::Real) = n -> pswap(n, i, j, α)

for (G, GT) in [
    (:Rx, :(PutBlock{N,1,RotationGate{1,T,XGate}} where {N,T})),
    (:Ry, :(PutBlock{N,1,RotationGate{1,T,YGate}} where {N,T})),
    (:Rz, :(PutBlock{N,1,RotationGate{1,T,ZGate}} where {N,T})),
    (:PSWAP, :(PSwap)),
]
    @eval function apply!(reg::AbstractRegister, g::$GT)
        instruct!(reg, Val($(QuoteNode(G))), g.locs, g.content.theta)
        return reg
    end
end


text:
export denseblocks

mutable struct DenseBlockIterator
    readers::Array{Any}
    blockSize::Int64
    blockWidth::Int64
    block::Array{Float64,2}
    offset::Int64
    done::Bool
    constantColumn::Bool
    loop::Bool
end

function denseblocks(readers, blockSize::Int64; constantColumn=false, loop=false)
    blockWidth = constantColumn ? length(readers) + 1 : length(readers)
    block = ones(Float64, blockSize, blockWidth)
    if constantColumn
        block[:,end] .= 1.0
    end
    DenseBlockIterator(readers, blockSize, blockWidth, block, 0, false, constantColumn, loop)
end

#Depricated for Julia 1.0
#Base.start(it::DenseBlockIterator) = 0
#Base.done(it::DenseBlockIterator, nil) = it.done && !it.loop # never done with a constant column

#function Base.next(it::DenseBlockIterator, nil)
function next(it::DenseBlockIterator, nil)

    if it.constantColumn
        it.block[:,1:end-1] .= 0.0
    else
        it.block[:,:] .= 0.0
    end

    # Fill in the block
    if !it.done

        foundRead = false
        for i in 1:length(it.readers)
            reader = it.readers[i]

            while !eof(reader) && position(reader) <= it.offset + it.blockSize
                it.block[position(reader) - it.offset, i] += value(reader)
                advance!(reader)
                foundRead = true
            end
        end

        # See if we are really done or just found a blank block
        if !foundRead
            it.done = true
            for i in 1:length(it.readers)
                it.done = it.done && eof(it.readers[i])
            end
        end
    end

    # update the offset
    it.offset += it.blockSize

    it.block, 0
end

function Base.iterate(it::DenseBlockIterator, state=0)
    if it.done && !it.loop
        return nothing
    else
        return next(it, state)
    end
end


text:
test01_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0499989727508645
        1.0397264813951748
    ],
    "V_I" => [
        2.380952380932277e-6
        0.023811904761703702
    ],
    "generator-102-1" => [
        0.16852719659039553
        1.0
    ],
)

test02_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031176878
        1.013581093277158
        1.0058147016829673
    ],
    "V_I" => [
        -8.347022011372903e-6
        -0.035395838898434936
        0.013070165151318975
    ],
    "generator-102-1" => [
        1.0767922333048168
        0.42832288470316265
        0.5202660714969694
        1.0
        2.417740387578741
        0.42897866591580786
        -0.43519326976417344
        1.0141989440244172
    ],
    "generator-103-1" => [
        0.8114924277843656
        0.6033931227898625
        0.8590280121051105
        1.0
        1.7138556549957245
        0.11317883762504757
        -0.3084940178992305
        1.0058996189176528
    ],
)

test03_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031176878
        1.0135810932772933
        1.005814701683217
    ],
    "V_I" => [
        -8.347022011402048e-6
        -0.03539583889890502
        0.013070165151052577
    ],
    "generator-102-1" => [
        1.1646938632048018
        0.4165011730835588
        1.027834348960626
        0.45807343235867964
        0.5202660714937425
        1.0
        2.4177403875765755
        0.4289786659140603
        -0.43519326976378364
        1.014198944024569
    ],
    "generator-103-1" => [
        0.8706440051085687
        0.5867394726000086
        0.7785472357034733
        0.6453037385810354
        0.8590280121029226
        1.0
        1.7138556549981903
        0.11317883762557865
        -0.3084940178996743
        1.0058996189178993
    ],
)

test04_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031176878
        1.0135810932778297
        1.0058147016833365
    ],
    "V_I" => [
        -8.34702201124049e-6
        -0.035395838896950645
        0.013070165152948175
    ],
    "generator-102-1" => [
        -0.5345748406265437
        0.861875418973418
        1.16469386320582
        0.4165011730882875
        1.0278343489597315
        0.45807343236388043
        0.5202660715024426
        1.0
        2.4177403875951065
        0.4289786659290129
        -0.43519326976711925
        1.0141989440250365
    ],
    "generator-103-1" => [
        -0.7530738934810763
        0.6668686184674737
        0.8706440051057203
        0.5867394726048358
        0.7785472356988481
        0.6453037385863445
        0.8590280121139345
        1.0
        1.7138556550116086
        0.11317883762846868
        -0.3084940179020896
        1.0058996189180431
    ],
)

test05_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031793623
        1.0135821475828384
        1.0058150887248778
    ],
    "V_I" => [
        -8.34704149711476e-6
        -0.03539621844488235
        0.013069953660280445
    ],
    "generator-102-1" => [
        1.265927354629895
        0.0
        1.1892309624086803
        0.1492441538451619
        0.3618503246775007
        1.0
        0.605565899879176
        0.012111599195329345
        -0.1090018619782517
        1.0142000109341534
    ],
    "generator-103-1" => [
        1.0350129715205239
        -7.90345313694394e-27
        0.9943720398921774
        0.19876653759719853
        0.5585603876595846
        1.0
        0.6850916263924302
        0.014604005256708285
        -0.12331649275063744
        1.0059000031789025
    ],
)

test06_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031793616
        1.013582147582662
        1.0058150887129174
    ],
    "V_I" => [
        -8.347041497038558e-6
        -0.035396218444544984
        0.013069953662795964
    ],
    "generator-102-1" => [
        -0.3919175747163956
        0.9354155636971591
        1.2659273546594862
        0.0
        1.1892309624328359
        0.14924415383937045
        0.3618503246616462
        1.0
        0.6055658998619685
        0.012111599194823144
        -0.10900186197515435
        1.0142000109339653
    ],
    "generator-103-1" => [
        -0.5219641597591145
        0.8598768704287161
        1.035012971442613
        1.786839370084875e-26
        0.9943720398216922
        0.1987665376172479
        0.5585603877306666
        1.0
        0.6850916263784529
        0.014604005256241823
        -0.12331649274812152
        1.0059000031669756
    ],
)

test07_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199954445260693
        1.0121032812279551
        1.0057639923670012
    ],
    "V_I" => [
        -1.0808535692176882e-5
        -0.06516685391537277
        -0.016523127201602223
    ],
    "generator-102-1" => [
        1.1131398694301942
        0.3370004397979819
        0.36332527065271414
        1.0
        2.30028045920761
        0.3438431783682635
        -0.4140504826573698
        1.0141990784464463
    ],
    "generator-103-1" => [
        0.8754905367233644
        0.5136808045551609
        0.6746234081353792
        1.0
        0.7403215561452087
        1.0
        0.7176423907384956
        1.0
        0.6914094421550748
        1.0
        0.6746234081353792
        1.0
        1.5017208842419558
        0.07552452190539284
        -0.2703097591635521
        1.0058997077614296
    ],
)

test08_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0000099992980975
        0.9999994707586406
    ],
    "V_I" => [
        6.874931250857116e-8
        0.0010313084369410759
    ],
    "generator-102-1" => [
        1.0
        0.1978641793142158
        0.04541511706491953
        0.0006674495903061723
        -0.00012674142328973991
        0.0703325301924954
        -0.006874509029034731
        1.0042596418725966
        -0.09824857969306856
        1.0090541173324878
        -1.3552527156068803e-23
        0.0
        0.1003426665815884
        0.49251497972474095
        0.07939915329276713
        1.0039784714785935
        0.10108135591344172
        0.49999500006233566
        0.005104746403351219
    ],
)

test09_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953029915412
        1.0135805859171578
        1.0058122208606277
    ],
    "V_I" => [
        -8.347002209520271e-6
        -0.03539559882532268
        0.013070652636616562
    ],
    "generator-102-1" => [
        1.076794967921721
        0.4283206474479955
        0.5202633660985048
        1.0
        2.4177531630463567
        0.4289889745060745
        -0.43519556934834425
        1.014198428594904
    ],
    "generator-103-1" => [
        1.0
        0.3749605370771409
        0.2800031488177231
        0.0012244421556582116
        -0.0005806861125950668
        0.07258390820648088
        -0.012677638073392193
        1.0352063062624437
        -0.18023788531288878
        1.0507796114445667
        3.30872245013387e-27
        -5.847986998874822e-29
        0.202580321610166
        0.9794556814296821
        0.008524280226319543
        1.0292918881467439
        0.21141428717781932
        0.9951003386808408
        -0.0676433194965395
    ],
)

test10_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199995546547758
        1.0138894577714104
        1.0046423611323474
    ],
    "V_I" => [
        -1.4828810811188414e-6
        -0.02509557882440288
        0.050277280031934914
    ],
    "generator-102-1" => [
        0.6478478687865317
        0.6672388799376429
        0.9385852366586044
        1.0
        1.0780870419914552
        0.033260063213494444
        -0.19405566755846196
        1.0141999904636838
    ],
    "generator-103-1" => [
        1.0
        0.4573111218123713
        -0.4453409115740306
        0.0013438428599890165
        0.00038704202037608665
        0.06150607054133235
        -0.013759754650798525
        0.8765256890716812
        -0.19781366899038325
        0.8985697141802758
        -1.0339757656912844e-28
        0.0
        0.23535024402517596
        0.7461517229537096
        0.7569627550747495
        0.8737986000941161
        0.20953171052512662
        0.7616570695325691
        0.6923016586677849
    ],
)

test12_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0200000009709098
        1.014199890107407
        1.0170999455391585
    ],
    "V_I" => [
        2.480607756117223e-8
        0.00048060307478791004
        0.00024031394043273134
    ],
    "generator-102-1" => [
        0.0004738860549890352
        0.9999999991950035
        8.049965538940007e-9
    ],
    "generator-101-1" => [
        9.889760280537507e-6
        0.9999999991950035
    ],
)

test13_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0199953031703135
        1.0135819409151956
        1.0058150851669077
    ],
    "V_I" => [
        -8.347038015388167e-6
        -0.03539614593111441
        0.013069982811242055
    ],
    "generator-102-1" => [
        1.076796543168192
        0.42832051528716575
        0.5202618423978124
        1.0
        2.4177538458928827
        1.8133153844196617
        -0.1087989230651797
        1.014199801861641
        1.0000000000007443
        1.0000000000007443
        0.7600000000005657
    ],
    "generator-103-1" => [
        0.8114925633188212
        0.6033933448331982
        0.8590278136367715
        1.0
        1.7138550022515144
    ],
)

test14_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0200000013609518
        1.0141998903792513
        1.0170999458701016
    ],
    "V_I" => [
        -2.3863302996949183e-8
        0.00048055378548599335
        0.00024026496109149398
    ],
    "generator-102-1" => [
        0.0004738340333812817
        0.9999999994250409
        5.749590553056069e-9
    ],
    "generator-101-1" => [
        0.9999999994250409
        -0.0006780551644918766
        0.07499970515180326
        1.134256731497628e-7
        -9.850337915313427e-5
        0.07234269603262689
        -3.169688095409349e-7
        1.0345002991161087
        -1.6696259078046452e-5
        1.034500299250843
        -8.195468443874853e-21
        6.645113843347989e-23
        -0.0006941946079418535
        8.74659973432602e-5
        0.004054476659144544
        1.0345000499848231
        -0.000718144471974579
        3.4323306447695294e-5
        -0.07249852699571752
    ],
)

test_psse_genrou_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992131
        1.0198615747329156
        0.9817964396228305
    ],
    "V_I" => [
        -1.9156999446242012e-12
        -0.016803815334368564
        -0.15145841258537568
    ],
    "generator-102-1" => [
        0.808414719298974
        0.5302614330069606
        0.7288775916772273
        0.7311890437210803
        0.9615872654389475
        1.0
    ],
)

test_psse_genrou_no_sat_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992116
        1.0198615747190438
        0.981796439616456
    ],
    "V_I" => [
        -1.9153833490517204e-12
        -0.016803815331282758
        -0.15145841258294843
    ],
    "generator-102-1" => [
        0.761127552111316
        0.5969928185692338
        0.6778994251511813
        0.778686285090305
        1.0290879496303451
        1.0
    ],
)

test_psse_genrou_high_sat_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992136
        1.0198615747361248
        0.9817964396243052
    ],
    "V_I" => [
        -1.9157731855354334e-12
        -0.01680381533508244
        -0.15145841258593723
    ],
    "generator-102-1" => [
        0.8848138955743352
        0.4029730032295409
        0.8129023592015162
        0.6364665865027268
        0.8388673542776348
        1.0
    ],
)

test_psse_genroe_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992122
        1.0198615747241417
        0.9817964396187986
    ],
    "V_I" => [
        -1.9154996983794634e-12
        -0.0168038153324168
        -0.15145841258384046
    ],
    "generator-102-1" => [
        0.7765987734123033
        0.5759873205986621
        0.6945076285567474
        0.7639104104022088
        1.007555994095225
        1.0
    ],
)

test_psse_genroe_high_sat_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992125
        1.0198615747274278
        0.9817964396203087
    ],
    "V_I" => [
        -1.9155746976060157e-12
        -0.0168038153331478
        -0.15145841258441545
    ],
    "generator-102-1" => [
        0.780985312007146
        0.5698915612533341
        0.6992284922760335
        0.7595916513872448
        1.0013586298474484
        1.0
    ],
)

test_psse_genrou_avr_init = Dict{String, Any}(
    "V_R" => [
        1.049999999999611
        1.0198615748018152
        0.9817964396548093
    ],
    "V_I" => [
        -9.586411659696706e-13
        -0.016803815348740852
        -0.15145841259655002
    ],
    "generator-102-1" => [
        0.7611275529260925
        0.5969928183151659
        0.6778994259655987
        0.7786862847589121
        1.0290879487397262
        1.0
        2.009549457333499
        0.1984434602979826
        -0.36171890232002984
        1.0199999999840748
    ],
)

test_psse_gensal_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992136
        1.0198615747359734
        0.9817964396242356
    ],
    "V_I" => [
        -1.91576973049762e-12
        -0.01680381533504876
        -0.15145841258591072
    ],
    "generator-102-1" => [
        1.0432937733330578
        0.8763054685303016
        -0.40384020950991173
        0.6194010895374469
        1.0
    ],
)

test_psse_gensae_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992138
        1.0198615747367235
        0.9817964396245803
    ],
    "V_I" => [
        -1.9157868491132568e-12
        -0.016803815335215615
        -0.15145841258604198
    ],
    "generator-102-1" => [
        1.0525147160721513
        0.8940751738672431
        0.37013399906022626
        0.5845318226206091
        1.0
    ],
)

test_psse_ac1a_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992227
        1.0198615748099735
        0.9817964396582409
    ],
    "V_I" => [
        -1.917458635215167e-12
        -0.01680381535151026
        -0.15145841259885914
    ],
    "generator-102-1" => [
        0.8084147201810992
        0.5302614325066757
        0.7288775925696391
        0.7311890432071789
        0.961587264395278
        1.0
        1.0199999999922775
        0.010765577988970956
        2.153115597794191
        2.153115597794191
        -0.05382788994485477
    ],
)

test_psse_ac1a_sat_init = Dict{String, Any}(
    "V_R" => [
        1.0499999999992224
        1.0198615748070743
        0.9817964396569085
    ],
    "V_I" => [
        -1.9173924649791205e-12
        -0.01680381535086531
        -0.15145841259835183
    ],
    "generator-102-1" => [
        0.8084147201479094
        0.5302614325254987
        0.7288775925360624
        0.7311890432265143
        0.9615872644345458
        1.0
        1.019999999989368
        0.013083531377287442
        2.6167062754574886
        2.1531155977360585
        -0.0654176568864372
    ],
)

test_psse_tgov1_init = Dict{String, Any}(
    "V_R" => [
        1.049999999999223
        1.0198615748129898
        0.9817964396596418
    ],
    "V_I" => [
        -1.9175011705717703e-12
        -0.01680381535196654
        -0.15145841259928494
    ],
    "generator-102-1" => [
        0.7611275530397932
        0.596992818279838
        0.6778994260789002
        0.7786862847128322
        1.0290879486156426
        1.0
        1.0199999999953009
        0.010047747287265943
        2.0095494574531885
        2.0095494574531885
        -0.0502387364363297
        1.0000000000027502
        0.6250000000017188
    ],
)

test_psse_gast_init = Dict{String, Any}(
    "V_R" => [
        1.049999999999223
        1.0198615748129447
        0.9817964396596628
    ],
    "V_I" => [
        -1.9174264624943186e-12
        -0.016803815351355156
        -0.15145841259899145
    ],
    "generator-102-1" => [
        0.7611275530385
        0.596992818281269
        0.6778994260767844
        0.7786862847146988
        1.029087948620471
        1.0
        1.0199999999952458
        0.010047747287321161
        2.009549457464232
        2.009549457464232
        -0.0502387364366058
        1.0000000000037501
        1.000000000009531
        1.000000000009531
    ],
)

test23_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0000100000000947
        0.9999994715292945
    ],
    "V_I" => [
        2.151752000488142e-18
        0.0010312396876354022
    ],
    "generator-102-1" => [
        0.19786411039169022
        0.5025002105975629
        0.04541511590746641
        0.0006674495901551267
        -0.00012674142163426048
        0.07033253022570093
        -0.006874509027216917
        1.0042596423477739
        -0.09824857967083465
        0.49251498484008666
        0.07939912055435007
        1.0039784789069033
        0.10108128683191751
        0.4999950000656485
        0.005104713115239256
    ],
)

test24_x0_init = Dict{String, Any}(
    "V_R" => [
        1.0000100000000953
        0.999999471539245
    ],
    "V_I" => [
        6.658282939054623e-19
        0.0010312396876131118
    ],
    "generator-102-1" => [
        0.0835937101719736
        0.500750063179269
        0.016666812918242877
        -0.004089048789580191
        0.011428671715366634
        0.05732140126078192
        8.170271546498093e-24
        -8.170271546498093e-25
        0.002043497690829951
        0.4948786405233513
        2.5088278260982038
        1.0014892471228272
        0.002046543812596204
        0.4999950000548418
        0.005104708291136156
    ],
)

test25_x0_init = Dict{String, Any}(
    "V_R" => [
        1.02
        1.0141688016598351
        0.980914345234313
    ],
    "V_I" => [
        5.5237548171800665e-8
        -0.007954982077548196
        -0.10878650410837654
    ],
    "generator-101-1" => [
        -0.7019759099414052
        0.740020148680559
        0.9598377510862162
        0.5469277035702602
        0.8604906642063641
        0.6015182348027084
        0.7590212970938108
        1.0000000006770182
        1.8694312317349822
        -6.770182514300934e-9
    ],
    "generator-102-1" => [
        -0.6633219887358436
        0.7672063463415573
        0.9194793489539811
        0.5168114274138487
        0.8506592041865252
        0.5683959607723741
        0.7050621597099493
        1.0000000006770182
        1.5495768886381311
    ],
    "Line BUS 1-BUS 3-i_1" => [
        0.9272581840805404
        -0.2484422743740139
    ],
    "Line BUS 1-BUS 2-i_2" => [
        0.06985630766987676
        -0.04277196052888417
    ],
    "Line BUS 2-BUS 3-i_3" => [
        0.8574018764106637
        -0.2056703138451297
    ],
)


text:
export minkowski_difference, pontryagin_difference

"""
    minkowski_difference(P::LazySet{N}, Q::LazySet{N}) where {N<:Real}

Concrete Minkowski difference (geometric difference) for a pair of
convex sets.

### Input

- `P` -- polytopic set
- `Q` -- compact convex set that is subtracted from `P`

### Output

An `HPolytope` that corresponds to the Minkowski difference of `P` minus `Q` if
`P` is bounded, and an `HPolyhedron` if `P` is unbounded.

### Notes

This function requires that the list of constraints of the set `P` is
available and that the set `Q` is bounded.

### Algorithm

This function implements Theorem 2.3 in [1], which we state next.

Suppose ``P`` is a polyhedron
```math
P = \\{z ∈ ℝ^n: sᵢᵀz ≤ rᵢ,~i = 1, …, N\\}.
```
where ``sᵢ ∈ ℝ^n, sᵢ ≠ 0``, and ``rᵢ ∈ ℝ``.
Assume ``ρ(sᵢ,Q)`` is defined for ``i = 1, …, N``. Then,

```math
P ⊖ Q = \\{z ∈ ℝ^n: sᵢᵀz ≤ rᵢ - ρ(sᵢ,Q),~i = 1, …, N\\}.
```

where ``⊖`` is defined as ``P ⊖ Q = \\{z ∈ ℝ^n: z + v ∈ P  ~∀~v ∈ Q\\}`` and is called
the *Minkowski difference* (also referenced as *Pontryagin difference*, or geometric difference).
It is denoted in [1] as the operation `P ~ Q`.

[1] Ilya Kolmanovsky and Elmer G. Gilbert (1997). *Theory and computation
of disturbance invariant sets for discrete-time linear systems.*
[Mathematical Problems in Engineering Volume 4, Issue 4, Pages
317-367.](http://dx.doi.org/10.1155/S1024123X98000866)
"""
function minkowski_difference(P::LazySet{N}, Q::LazySet{N}) where {N<:Real}

    @assert applicable(constraints_list, P)  "this function " *
        "requires that the list of constraints of its first argument is applicable, but it is not; " *
        "if it is bounded, try overapproximating with an `HPolytope` first"
    @assert isbounded(Q) "this function requires that its second argument is bounded, but it is not"

    A, b = tosimplehrep(P)
    g_PminusQ = [b[i] - ρ(A[i, :], Q) for i in eachindex(b)]
    if isbounded(P)
        return HPolytope(A, g_PminusQ)
    else
        return HPolyhedron(A, g_PminusQ)
    end
end

"""
    pontryagin_difference(P::LazySet{N}, Q::LazySet{N}) where {N<:Real}

An alias for the function `minkowski_difference`.

### Notes

Due to inconsistent naming conventions, both the name *Minkowski difference* and
*Pontryagin difference* are used to refer to the geometric difference of two sets.

"""
const pontryagin_difference = minkowski_difference

# concrete minkowski difference with singleton
minkowski_difference(X::LazySet, S::AbstractSingleton) = translate(X, -element(S))
minkowski_difference(X::LazySet, ::ZeroSet) = X


text:
include("defs.jl")

# Zooming into standard map
using DynamicalSystems
figure = Figure(resolution = (800,800))
ax = figure[1,1] = Axis(figure)
deactivate_interaction!(ax, :rectanglezoom)
rect = select_rectangle(ax.scene)
sm = Systems.standardmap(; k = 1.0)
g = 10 # grid density
inrect(u, r) = (r.origin[1] ≤ u[1] ≤ r.origin[1]+r.widths[1]) && (r.origin[2] ≤ u[2] ≤ r.origin[2] + r.widths[2])
integ = integrator(sm)
marker = InteractiveDynamics.MARKER

N = 100 # how many points to plot inside the rect
Nslider = labelslider!(
    figure, "N =", round.(Int, range(100, 100000; length = 100))
)
N = Nslider.slider.value

recompute = Button(figure; label = "recompute")

figure[2,1][1,1] = recompute
figure[2,1][1,2] = Nslider.layout

sco = Observable([Point2f0(sm.u0)])

scatter!(ax, sco; markersize = 2*px, strokewidth = 0, color = :black)

on(rect) do r # new subgrid selected
    xs = range(r.origin[1], r.origin[1]+r.widths[1]; length = g)
    ys = range(r.origin[2], r.origin[2]+r.widths[2]; length = g)
    s = sco[]
    n = N[]
    empty!(s)
    for x in xs
        for y in ys
            DynamicalSystems.reinit!(integ, SVector(x, y))
            push!(s, integ.u)
            i = 0
            for i in 1:n
                DynamicalSystems.step!(integ)
                if inrect(integ.u, r)
                    push!(s, integ.u)
                end
            end
        end
    end
    sco[] = s
    xlims!(ax, xs[1], xs[end])
    ylims!(ax, ys[1], ys[end])
end

on(recompute.clicks) do c
    rect[] = rect[]
end

rect[] = FRect2D([0,0],[2π,2π])
hidexdecorations!(ax)
hideydecorations!(ax)

display(figure)

# record_interaction(figure, joinpath(@__DIR__, "sm_zoom.mp4"); total_time = 20)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
macro hierarchical_dict()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, InnerNode, nodeshow, children
                   NodeType(::Type{<:Dict}) = InnerNode()
                   children(d::Dict) = d
                   nodeshow(io::IO, d::Dict) = print(io, isempty(d) ? "Empty Dict" : "Dict of")
                   @specialize
                   nothing
               end
              )
end

macro hierarchical_vector()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, InnerNode, nodeshow, children
                   NodeType(::Type{<:Vector}) = InnerNode()
                   children(v::Vector) = v
                   nodeshow(io::IO, v::Vector) = print(io, isempty(v) ? "[]" : "Vector of")
                   @specialize
                   nothing
               end
              )
end

macro hierarchical_tuple()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, InnerNode, nodeshow, children
                   NodeType(::Type{<:Tuple}) = InnerNode()
                   children(v::Tuple) = v
                   nodeshow(io::IO, v::Tuple) = print(io, isempty(v) ? "()" : "Tuple of")
                   @specialize
                   nothing
               end
              )
end

macro hierarchical_namedtuple()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, InnerNode, nodeshow, children
                   NodeType(::Type{<:NamedTuple}) = InnerNode()
                   children(v::NamedTuple) = v
                   nodeshow(io::IO, v::NamedTuple) = print(io, isempty(v) ? "()" : "NamedTuple of")
                   @specialize
                   nothing
               end
              )
end

macro hierarchical_pairvector()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, InnerNode, nodeshow, children
                   NodeType(::Type{<:HierarchicalUtils.PairVec}) = InnerNode()
                   children(v::HierarchicalUtils.PairVec) = v
                   function nodeshow(io::IO, v::HierarchicalUtils.PairVec)
                       print(io, isempty(v) ? "Empty Vector" : "Vector of pairs of")
                   end
                   @specialize
                   nothing
               end
              )
end

macro primitives()
    return esc(
               quote
                   @nospecialize
                   import HierarchicalUtils: NodeType, LeafNode
                   NodeType(::Type{<:Number}) = LeafNode()
                   NodeType(::Type{<:AbstractString}) = LeafNode()
                   NodeType(::Type{<:AbstractChar}) = LeafNode()
                   NodeType(::Type{Symbol}) = LeafNode()
                   NodeType(::Type{Missing}) = LeafNode()
                   NodeType(::Type{Nothing}) = LeafNode()
                   @specialize
                   nothing
               end
              )
end


text:
module JuliaBlackBoard

using CommonMark
using Mustache
using Base64

# for LaTeX
using Tectonic
using ImageMagick_jll
using TtH_jll

# Exports
export question
export MC, MA, TF, ESS, ORD, MAT, FIB, FIB_PLUS, FIL, NUM, SR, OP, JUMBLED_SENTENCE #, QUIZ_BOWL
export Plot, File, Tikz, LaTeX, preview, LaTeX′
export lquestion, mdquestion # tentative
export @ltx_str, @mdltx_str, @MT_str, @mt_str # reexoprt

include("mustache-additions.jl")
include("questions.jl")

end


text:
include("../AWSCorePrototypeServices.jl")
using .Services: xray

"""
    ListTagsForResource()

Required Parameters
{
  "ResourceARN": ""
}

Optional Parameters
{
  "NextToken": ""
}
"""
ListTagsForResource(args) = xray("POST", "/ListTagsForResource", args)

"""
    DeleteGroup()

Deletes a group resource.

Optional Parameters
{
  "GroupARN": "The ARN of the group that was generated on creation.",
  "GroupName": "The case-sensitive name of the group."
}
"""

DeleteGroup() = xray("POST", "/DeleteGroup")
DeleteGroup(args) = xray("POST", "/DeleteGroup", args)

"""
    PutTraceSegments()

Uploads segment documents to AWS X-Ray. The X-Ray SDK generates segment documents and sends them to the X-Ray daemon, which uploads them in batches. A segment document can be a completed segment, an in-progress segment, or an array of subsegments. Segments must include the following fields. For the full segment document schema, see AWS X-Ray Segment Documents in the AWS X-Ray Developer Guide.  Required Segment Document Fields     name - The name of the service that handled the request.    id - A 64-bit identifier for the segment, unique among segments in the same trace, in 16 hexadecimal digits.    trace_id - A unique identifier that connects all segments and subsegments originating from a single client request.    start_time - Time the segment or subsegment was created, in floating point seconds in epoch time, accurate to milliseconds. For example, 1480615200.010 or 1.480615200010E9.    end_time - Time the segment or subsegment was closed. For example, 1480615200.090 or 1.480615200090E9. Specify either an end_time or in_progress.    in_progress - Set to true instead of specifying an end_time to record that a segment has been started, but is not complete. Send an in progress segment when your application receives a request that will take a long time to serve, to trace the fact that the request was received. When the response is sent, send the complete segment to overwrite the in-progress segment.   A trace_id consists of three numbers separated by hyphens. For example, 1-58406520-a006649127e371903a2de979. This includes:  Trace ID Format    The version number, i.e. 1.   The time of the original request, in Unix epoch time, in 8 hexadecimal digits. For example, 10:00AM December 2nd, 2016 PST in epoch time is 1480615200 seconds, or 58406520 in hexadecimal.   A 96-bit identifier for the trace, globally unique, in 24 hexadecimal digits.

Required Parameters
{
  "TraceSegmentDocuments": "A string containing a JSON document defining one or more segments or subsegments."
}
"""
PutTraceSegments(args) = xray("POST", "/TraceSegments", args)

"""
    GetGroup()

Retrieves group resource details.

Optional Parameters
{
  "GroupARN": "The ARN of the group that was generated on creation.",
  "GroupName": "The case-sensitive name of the group."
}
"""

GetGroup() = xray("POST", "/GetGroup")
GetGroup(args) = xray("POST", "/GetGroup", args)

"""
    GetEncryptionConfig()

Retrieves the current encryption configuration for X-Ray data.
"""

GetEncryptionConfig() = xray("POST", "/EncryptionConfig")
GetEncryptionConfig(args) = xray("POST", "/EncryptionConfig", args)

"""
    GetTraceGraph()

Retrieves a service graph for one or more specific trace IDs.

Required Parameters
{
  "TraceIds": "Trace IDs of requests for which to generate a service graph."
}

Optional Parameters
{
  "NextToken": "Pagination token."
}
"""
GetTraceGraph(args) = xray("POST", "/TraceGraph", args)

"""
    CreateSamplingRule()

Creates a rule to control sampling behavior for instrumented applications. Services retrieve rules with GetSamplingRules, and evaluate each rule in ascending order of priority for each request. If a rule matches, the service records a trace, borrowing it from the reservoir size. After 10 seconds, the service reports back to X-Ray with GetSamplingTargets to get updated versions of each in-use rule. The updated rule contains a trace quota that the service can use instead of borrowing from the reservoir.

Required Parameters
{
  "SamplingRule": "The rule definition."
}

Optional Parameters
{
  "Tags": ""
}
"""
CreateSamplingRule(args) = xray("POST", "/CreateSamplingRule", args)

"""
    PutEncryptionConfig()

Updates the encryption configuration for X-Ray data.

Required Parameters
{
  "Type": "The type of encryption. Set to KMS to use your own key for encryption. Set to NONE for default encryption."
}

Optional Parameters
{
  "KeyId": "An AWS KMS customer master key (CMK) in one of the following formats:    Alias - The name of the key. For example, alias/MyKey.    Key ID - The KMS key ID of the key. For example, ae4aa6d49-a4d8-9df9-a475-4ff6d7898456. AWS X-Ray does not support asymmetric CMKs.    ARN - The full Amazon Resource Name of the key ID or alias. For example, arn:aws:kms:us-east-2:123456789012:key/ae4aa6d49-a4d8-9df9-a475-4ff6d7898456. Use this format to specify a key in a different account.   Omit this key if you set Type to NONE."
}
"""
PutEncryptionConfig(args) = xray("POST", "/PutEncryptionConfig", args)

"""
    PutTelemetryRecords()

Used by the AWS X-Ray daemon to upload telemetry.

Required Parameters
{
  "TelemetryRecords": ""
}

Optional Parameters
{
  "Hostname": "",
  "ResourceARN": "",
  "EC2InstanceId": ""
}
"""
PutTelemetryRecords(args) = xray("POST", "/TelemetryRecords", args)

"""
    TagResource()

Required Parameters
{
  "ResourceARN": "",
  "Tags": ""
}
"""
TagResource(args) = xray("POST", "/TagResource", args)

"""
    UntagResource()

Required Parameters
{
  "ResourceARN": "",
  "TagKeys": ""
}
"""
UntagResource(args) = xray("POST", "/UntagResource", args)

"""
    GetTimeSeriesServiceStatistics()

Get an aggregation of service statistics defined by a specific time range.

Required Parameters
{
  "StartTime": "The start of the time frame for which to aggregate statistics.",
  "EndTime": "The end of the time frame for which to aggregate statistics."
}

Optional Parameters
{
  "Period": "Aggregation period in seconds.",
  "GroupARN": "The ARN of the group for which to pull statistics from.",
  "NextToken": "Pagination token.",
  "EntitySelectorExpression": "A filter expression defining entities that will be aggregated for statistics. Supports ID, service, and edge functions. If no selector expression is specified, edge statistics are returned. ",
  "GroupName": "The case-sensitive name of the group for which to pull statistics from."
}
"""
GetTimeSeriesServiceStatistics(args) = xray("POST", "/TimeSeriesServiceStatistics", args)

"""
    GetSamplingTargets()

Requests a sampling quota for rules that the service is using to sample requests.

Required Parameters
{
  "SamplingStatisticsDocuments": "Information about rules that the service is using to sample requests."
}
"""
GetSamplingTargets(args) = xray("POST", "/SamplingTargets", args)

"""
    UpdateSamplingRule()

Modifies a sampling rule's configuration.

Required Parameters
{
  "SamplingRuleUpdate": "The rule and fields to change."
}
"""
UpdateSamplingRule(args) = xray("POST", "/UpdateSamplingRule", args)

"""
    BatchGetTraces()

Retrieves a list of traces specified by ID. Each trace is a collection of segment documents that originates from a single request. Use GetTraceSummaries to get a list of trace IDs.

Required Parameters
{
  "TraceIds": "Specify the trace IDs of requests for which to retrieve segments."
}

Optional Parameters
{
  "NextToken": "Pagination token."
}
"""
BatchGetTraces(args) = xray("POST", "/Traces", args)

"""
    CreateGroup()

Creates a group resource with a name and a filter expression.

Required Parameters
{
  "GroupName": "The case-sensitive name of the new group. Default is a reserved name and names must be unique."
}

Optional Parameters
{
  "Tags": "",
  "FilterExpression": "The filter expression defining criteria by which to group traces."
}
"""
CreateGroup(args) = xray("POST", "/CreateGroup", args)

"""
    GetSamplingRules()

Retrieves all sampling rules.

Optional Parameters
{
  "NextToken": "Pagination token."
}
"""

GetSamplingRules() = xray("POST", "/GetSamplingRules")
GetSamplingRules(args) = xray("POST", "/GetSamplingRules", args)

"""
    DeleteSamplingRule()

Deletes a sampling rule.

Optional Parameters
{
  "RuleARN": "The ARN of the sampling rule. Specify a rule by either name or ARN, but not both.",
  "RuleName": "The name of the sampling rule. Specify a rule by either name or ARN, but not both."
}
"""

DeleteSamplingRule() = xray("POST", "/DeleteSamplingRule")
DeleteSamplingRule(args) = xray("POST", "/DeleteSamplingRule", args)

"""
    GetGroups()

Retrieves all active group details.

Optional Parameters
{
  "NextToken": "Pagination token."
}
"""

GetGroups() = xray("POST", "/Groups")
GetGroups(args) = xray("POST", "/Groups", args)

"""
    GetSamplingStatisticSummaries()

Retrieves information about recent sampling results for all sampling rules.

Optional Parameters
{
  "NextToken": "Pagination token."
}
"""

GetSamplingStatisticSummaries() = xray("POST", "/SamplingStatisticSummaries")
GetSamplingStatisticSummaries(args) = xray("POST", "/SamplingStatisticSummaries", args)

"""
    GetServiceGraph()

Retrieves a document that describes services that process incoming requests, and downstream services that they call as a result. Root services process incoming requests and make calls to downstream services. Root services are applications that use the AWS X-Ray SDK. Downstream services can be other applications, AWS resources, HTTP web APIs, or SQL databases.

Required Parameters
{
  "StartTime": "The start of the time frame for which to generate a graph.",
  "EndTime": "The end of the timeframe for which to generate a graph."
}

Optional Parameters
{
  "GroupARN": "The ARN of a group to generate a graph based on.",
  "NextToken": "Pagination token.",
  "GroupName": "The name of a group to generate a graph based on."
}
"""
GetServiceGraph(args) = xray("POST", "/ServiceGraph", args)

"""
    GetTraceSummaries()

Retrieves IDs and annotations for traces available for a specified time frame using an optional filter. To get the full traces, pass the trace IDs to BatchGetTraces. A filter expression can target traced requests that hit specific service nodes or edges, have errors, or come from a known user. For example, the following filter expression targets traces that pass through api.example.com:  service("api.example.com")  This filter expression finds traces that have an annotation named account with the value 12345:  annotation.account = "12345"  For a full list of indexed fields and keywords that you can use in filter expressions, see Using Filter Expressions in the AWS X-Ray Developer Guide.

Required Parameters
{
  "StartTime": "The start of the time frame for which to retrieve traces.",
  "EndTime": "The end of the time frame for which to retrieve traces."
}

Optional Parameters
{
  "SamplingStrategy": "A paramater to indicate whether to enable sampling on trace summaries. Input parameters are Name and Value.",
  "NextToken": "Specify the pagination token returned by a previous request to retrieve the next page of results.",
  "TimeRangeType": "A parameter to indicate whether to query trace summaries by TraceId or Event time.",
  "FilterExpression": "Specify a filter expression to retrieve trace summaries for services or requests that meet certain requirements.",
  "Sampling": "Set to true to get summaries for only a subset of available traces."
}
"""
GetTraceSummaries(args) = xray("POST", "/TraceSummaries", args)

"""
    UpdateGroup()

Updates a group resource.

Optional Parameters
{
  "GroupARN": "The ARN that was generated upon creation.",
  "FilterExpression": "The updated filter expression defining criteria by which to group traces.",
  "GroupName": "The case-sensitive name of the group."
}
"""

UpdateGroup() = xray("POST", "/UpdateGroup")
UpdateGroup(args) = xray("POST", "/UpdateGroup", args)

max_stars_repo_name: mattBrzezinski/Julia-AWS-Status

max_issues_repo_name: mattBrzezinski/Julia-AWS-Status

max_forks_repo_name: mattBrzezinski/Julia-AWS-Status


text:
#Accuracy -> efficiency
for acc in 0:2:100
    noToken = mean([diversityWithChallenge(1000, 50, acc, 10, [Benchmarks.registryMean]) for i in 1:20])[1]
    token = mean([binaryTokenChallenge(1000, 50, acc, 10, [Benchmarks.registryMean]) for i in 1:20])[1]
    println("$acc, $noToken, $token")
end

#Agents of steps -> efficiency
for agents in 10:10:200
    noToken = mean([diversityWithChallenge(1000, agents, 50, 20, [Benchmarks.registryMean]) for i in 1:20])[1]
    token = mean([binaryTokenChallenge(1000, agents, 50, 20, [Benchmarks.registryMean]) for i in 1:20])[1]
    println("$agents, $noToken, $token")
end

#Number of steps -> efficiency
for steps in 2000:100:3000
    noToken = mean([diversityWithChallenge(steps, 50, 50, 20, [Benchmarks.registryMean]) for i in 1:20])[1]
    token = mean([binaryTokenChallenge(steps, 50, 50, 20, [Benchmarks.registryMean]) for i in 1:20])[1]
    println("$steps, $noToken, $token")
end

for acc in 0:10:100
    token = mean([binaryTokenChallenge(1000, 100, acc, 20, [Benchmarks.registryMean]) for i in 1:20])[1]
    println("$acc, $token")
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
function simpleDFT(hk::Union{Array{T}, SubArray{T}}) where T
    N = length(hk)
    if N == 2
        return [Complex(hk[1] + hk[2]), Complex(hk[1] - hk[2])]
    end
    ev = simpleDFT(@view hk[1:2:end]) #even
    od = simpleDFT(@view hk[2:2:end]) #odd

    twiddle = exp(-2 * π * im / N)
    w = one(Complex{T}) #twiddle^0

    Hk = Array{Complex{T}}(undef, N)
    NHalf = Int(N / 2)
    @inbounds @simd for i = 1:NHalf
        Hk[i] = ev[i] + w * od[i]
        Hk[i + NHalf] = ev[i] - w * od[i] #took me longer than i care to admit - had to look it up
        w = w * twiddle
    end

    return Hk
end

max_stars_repo_name: JulienKluge/ComputationalPhotonics

max_issues_repo_name: JulienKluge/ComputationalPhotonics

max_forks_repo_name: JulienKluge/ComputationalPhotonics


text:
# # Basic usage of `CoreModel` and `CoreModelCoupled`

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# In this tutorial we will introduce `COBREXA`'s `CoreModel` and
# `CoreModelCoupled`. We will use *E. coli*'s toy model to start with.

!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")

using COBREXA

# ## Loading a `CoreModel`

model = load_model(CoreModel, "e_coli_core.xml") # we specifically want to load a CoreModel from the model file

# ## Basic analysis on `CoreModel`

# As before, for optimization based analysis we need to load an optimizer. Here we
# will use [`Tulip.jl`](https://github.com/ds4dm/Tulip.jl) to optimize the linear
# programs of this tutorial. Refer to the constraint-based analysis basics
# tutorial if you are confused by any functions in this section.

# All the normal analysis functions work on `CoreModel`, due to it also having
# the same generic accessor interface as all the other model types.

using Tulip

dict_sol = flux_balance_analysis_dict(
    model,
    Tulip.Optimizer;
    modifications = [
        change_objective("R_BIOMASS_Ecoli_core_w_GAM"),
        change_constraint("R_EX_glc__D_e", -12, -12),
        change_constraint("R_EX_o2_e", 0, 0),
    ],
)

# ## Structure of `CoreModel`

# `CoreModel` is a special `COBREXA` type that is optimized for large scale
# analysis of large models. It stores data in a sparse format where possible.

fieldnames(CoreModel)
#
model.S

# ## `CoreModelCoupled` adds coupling constraints to `CoreModel`

# `CoreModelCoupled` extends `CoreModel` by adding coupling constraints.

fieldnames(CoreModelCoupled)

# In short, coupling constraints can be used to ensure that fluxes scale with
# the growth rate (μ) of a model. This reduces the impact of biologically
# infeasible cycles from occurring. Here we will model coupling constraints by
# assuming that they have the form: `-γ ≤ vᵢ/μ  ≤ γ`, where `γ` is the ratio
# between each individual flux (vᵢ) in the model and the growth rate.

gamma = 40 # arbitrary

nr = n_reactions(model) # number of reactions
biomass_index = first(indexin(["R_BIOMASS_Ecoli_core_w_GAM"], reactions(model)))

using LinearAlgebra, SparseArrays

Cf = sparse(1.0I, nr, nr)
Cf[:, biomass_index] .= -gamma
Cb = sparse(1.0I, nr, nr)
Cb[:, biomass_index] .= gamma
C = [Cf; Cb] # coupling constraint matrix

clb = spzeros(2 * nr)
clb[1:nr] .= -1000.0
cub = spzeros(2 * nr)
cub[nr+1:end] .= 1000

cmodel = CoreModelCoupled(model, C, clb, cub)
#
d = flux_balance_analysis_dict(model, Tulip.Optimizer)
d["R_BIOMASS_Ecoli_core_w_GAM"]
#
dc = flux_balance_analysis_dict(cmodel, Tulip.Optimizer)
dc["R_BIOMASS_Ecoli_core_w_GAM"]


text:
"Compute etendue of the instrument"
function etendue(instrument::GratingNoiseModel)
    @unpack detector_size, Fnumber = instrument
    "Compute detector area"
    A = detector_size^2;
    "Compute solid angle Ω"
    Ω = uconvert(u"sr",π/(4Fnumber^2))
    "Compute etendue (same unit as area A)"
    A * Ω
end

"Compute photons converted to electrons at the FPA"
function photons_at_fpa(ins, λ, radiance)
    @unpack t_int, Δλ, Qₑ, η = ins
    # Common unit:
    #luminosity = uconvert.(u"J/s/μm^2/sr/μm",radiance)
    energy = radiance * etendue(ins) * t_int * Δλ * Qₑ * η
    uconvert.(u"nm^-1", energy, Spectral()) .* λ
end

function noise_at_fpa(ins, photons)
    @unpack t_int, dark_current, σ_read = ins
    total_electrons = photons .+ dark_current .* t_int
    # Sum variances; Can add more noise terms here later if needed!
    total_variance = total_electrons .+ σ_read^2
    # Return 1-σ noise level
    sqrt.(total_variance)
end

"Compute noise equivalent spectral radiance"
function noise_equivalent_radiance(ins, λ, radiance)
    pfpa  = photons_at_fpa(ins, λ, radiance)
    noise = noise_at_fpa(ins,pfpa)
    SNR = pfpa./noise
    return radiance./SNR
end

function noise_components(ins, λ, radiance)
    @unpack t_int, dark_current, σ_read = ins
    photons  = photons_at_fpa(ins, λ, radiance)
    σ_shot           = sqrt.(photons)
    σ_shot_dark      = sqrt.(photons .+ dark_current .* t_int)
    σ_shot_dark_read = sqrt.(photons .+ dark_current .* t_int .+ σ_read^2)
    return [σ_shot σ_shot_dark σ_shot_dark_read]
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    html_p(id="tcid", children = "Hello Disable UI"),
    dcc_graph(id="broken", animate = 3)
end

run_server(app, debug = true,
    dev_tools_hot_reload = false, dev_tools_ui = false
)