text:
# # [Lotka-Volterra Three Ways](@id lotka_volterra)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Lotka-Volterra.ipynb)

# We will give three examples of deriving the Lotka-Volterra equations $$ \dot r(t) = \alpha r(t) - \beta r(t)f(t), \dot f(t) = \gamma f(t) r(t) - \delta f(t)$$ as a composition of primitive systems. First, we will show it as the composition of resource sharers where the undirected composition pattern is an undirected wiring diagram. Second, we will show it as the composition of machines where the directed composition pattern is a wiring diagram. Lastly, we will show it as the composition of machines where the directed composition pattern is an open circular port graph (open CPG).

# Although these examples are quite small and you can easily work them out by hand,  this approach is based on the solid mathematical theory of operad algebras that allows you to scale to larger and more complex problems.

#-

# ## Undirected composition
# ### Undirected wiring diagrams

# A standard Lotka-Volterra predator-prey model is the composition of three resource sharers:

# 1. a model of rabbit growth: this resource sharer has dynamics $\dot r(t) = \alpha r(t)$ and one port which exposes the rabbit population.
# 2. a model of rabbit/fox predation: this resource sharer has dynamics $$\dot r(t) = -\beta r(t) f(t), \dot f(t) = \gamma r(t)f(t)$$ and two ports which expose the rabbit and fox populations respectively
# 3. a model of fox population decline: this resource sharer has dynamics $\dot f(t) = -\delta f(t)$ and one port which exposes the fox population.

# However, there are not two independent rabbit populations --- one that grows and one that gets eaten by foxes. Likewise, there are not two independent fox populations --- one that declines and one that feasts on rabbits. To capture these interactions between the trio of resource sharers, we compose them by identifying the exposed rabbit populations and identifying the exposed fox populations.

using AlgebraicDynamics
using AlgebraicDynamics.UWDDynam
using Catlab.WiringDiagrams, Catlab.Programs
using LabelledArrays
using OrdinaryDiffEq, Plots, Plots.PlotMeasures

const UWD = UndirectedWiringDiagram

## Define the primitive systems
dotr(u,p,t) = p.α*u
dotrf(u,p,t) = [-p.β*u[1]*u[2], p.γ*u[1]*u[2]]
dotf(u,p,t) = -p.δ*u

rabbit_growth = ContinuousResourceSharer{Float64}(1, dotr)
rabbitfox_predation = ContinuousResourceSharer{Float64}(2, dotrf)
fox_decline = ContinuousResourceSharer{Float64}(1, dotf)

## Define the composition pattern
rf = @relation (rabbits,foxes) begin
    growth(rabbits)
    predation(rabbits,foxes)
    decline(foxes)
end

## Compose
rabbitfox_system = oapply(rf, [rabbit_growth, rabbitfox_predation, fox_decline])

## Solve and plot
u0 = [10.0, 100.0]
params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7)
tspan = (0.0, 100.0)

prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, rabbitfox_system,
    lw=2, title = "Lotka-Volterra Predator-Prey Model",
    xlabel = "time", ylabel = "population size"
)

# ## Directed composition
# A standard Lotka-Volterra predator-prey model is the composition of two machines:

# 1. Evolution of a rabbit population: this machine has one input which represents a population of predators, $h(t)$, that hunt rabbits. This machine has one output which emits the rabbit population $r(t)$. The dynamics of this machine is the driven ODE $\dot r(t) = \alpha r(t) - \beta r(t) h(t).$

# 2. Evoluation of a fox population: this machine has one input which represents a population of prey, $e(t)$, that are eaten by foxes. This machine has one output which emits the fox population $f(t)$. The dynamics of this machine is the driven ODE $\dot f(t) =\gamma f(t)e(t) - \delta f(t).$

# ### Wiring diagrams
# Since foxes hunt rabbit, these machines compose by setting the fox population to be the input for rabbit evolution. Likewise, we set the rabbit population to be the input for fox evolution. We first implement this setting of inputs using a wiring diagram as our composition pattern.

using AlgebraicDynamics.DWDDynam

## Define the primitive systems
dotr(u, x, p, t) = [p.α*u[1] - p.β*u[1]*x[1]]
dotf(u, x, p, t) = [p.γ*u[1]*x[1] - p.δ*u[1]]

rabbit = ContinuousMachine{Float64}(1,1,1, dotr, (u, p, t) -> u)
fox    = ContinuousMachine{Float64}(1,1,1, dotf, (u, p, t) -> u)

## Define the composition pattern
rabbitfox_pattern = WiringDiagram([], [:rabbits, :foxes])
rabbit_box = add_box!(rabbitfox_pattern, Box(:rabbit, [:pop], [:pop]))
fox_box = add_box!(rabbitfox_pattern, Box(:fox, [:pop], [:pop]))

add_wires!(rabbitfox_pattern, Pair[
    (rabbit_box, 1) => (fox_box, 1),
    (fox_box, 1)    => (rabbit_box, 1),
    (rabbit_box, 1) => (output_id(rabbitfox_pattern), 1),
    (fox_box, 1)    => (output_id(rabbitfox_pattern), 2)
])

## Compose
rabbitfox_system = oapply(rabbitfox_pattern, [rabbit, fox])

## Solve and plot
u0 = [10.0, 100.0]
params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7)
tspan = (0.0, 100.0)

prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, rabbitfox_system, params,
    lw=2, title = "Lotka-Volterra Predator-Prey Model",
    xlabel = "time", ylabel = "population size"
)

# ### Open CPG
# We next implement the setting of inputs using an open CPG as our composition pattern. We will use a barbell CPG. A barbell has two boxes connected by $n$ wires. In this instance we will set $n$ to $1$ since each machine receives and emits exactly one piece of information.

using AlgebraicDynamics.CPortGraphDynam
using AlgebraicDynamics.CPortGraphDynam: barbell

## Define the composition pattern
rabbitfox_pattern = barbell(1)

## Compose
rabbitfox_system = oapply(rabbitfox_pattern, [rabbit, fox])

## Solve and plot
u0 = [10.0, 100.0]
params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7)
tspan = (0.0, 100.0)

prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, FRK65(0))

plot(sol, lw=2, bottom_margin=10mm, left_margin=10mm, title = "Lotka-Volterra Predator-Prey Model", label=["rabbits" "foxes"])
xlabel!("Time")
ylabel!("Population size")


text:
@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}, Message{Gaussian}),
                :name          => SPAdditionOutNGG)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Message{Gaussian}, Nothing),
                :name          => SPAdditionIn2GGN)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Nothing, Message{Gaussian}),
                :name          => SPAdditionIn1GNG)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}, Message{PointMass}),
                :name          => SPAdditionOutNGP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{PointMass}, Message{Gaussian}),
                :name          => SPAdditionOutNPG)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{PointMass}, Nothing, Message{Gaussian}),
                :name          => SPAdditionIn1PNG)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Nothing, Message{PointMass}),
                :name          => SPAdditionIn1GNP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{PointMass}, Message{Gaussian}, Nothing),
                :name          => SPAdditionIn2PGN)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Message{PointMass}, Nothing),
                :name          => SPAdditionIn2GPN)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{PointMass},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPAdditionOutNPP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{PointMass},
                :inbound_types => (Message{PointMass}, Message{PointMass}, Nothing),
                :name          => SPAdditionIn2PPN)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{PointMass},
                :inbound_types => (Message{PointMass}, Nothing, Message{PointMass}),
                :name          => SPAdditionIn1PNP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message{SampleList}, Message{PointMass}),
                :name          => SPAdditionOutNSP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message{PointMass}, Message{SampleList}),
                :name          => SPAdditionOutNPS)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{PointMass}, Nothing, Message{SampleList}),
                :name          => SPAdditionIn1PNS)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{SampleList}, Nothing, Message{PointMass}),
                :name          => SPAdditionIn1SNP)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{PointMass}, Message{SampleList}, Nothing),
                :name          => SPAdditionIn2PSN)

@sumProductRule(:node_type     => Addition,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{SampleList}, Message{PointMass}, Nothing),
                :name          => SPAdditionIn2SPN)

@marginalRule(:node_type     => Addition,
              :inbound_types => (Nothing, Message{Gaussian}, Message{Gaussian}), # "Nothing" indicates marginalization over outbound edge
              :name          => MAdditionNGG)

ext: jl


text:
include("downsample.jl")

export downsample1
export downsample2
export downsample3


text:
"""
    out = abc_mcmc(theta0, prior, s_obs, sim!, n; [tol=0.0, b=0.25n, adapt_cov=true,
             adapt_tol=true, cutoff=nothing, sc=1.0, tol_targetAcceptance=0.10,
             rng=MersenneTwister(Random.make.seed())], dist=ABCdefaultDist,
             record_summaries=true)

Sample from ABC posterior using tolerance adaptive ABC-MCMC.

# Arguments
- `theta_0`: Initial value of parameter vector
- `prior`: Function returning prior log density values.
- `s_obs`: Observed summary statistics vector
- `sim!`: Function simulating summaries; sim!(s, θ, rng) simulates observed
          summaries corresponding to parameters θ and records them to s.
- `n`: Length of MCMC chain.
- `tol`: Tolerance (initial tolerance if `adapt_tol=true`).
- `b`: Burn-in length. The total number of simulations is b+n.
- `adapt_cov`: Whether covariance adaptation is done.
- `adapt_tol`: Whether tolerance adaptation is done.
- `cutoff`: Cut-off function ϕ. If `nothing`, simple cut-off is used,
            "gaussian" => Gaussian, "epa" => Epanechnikov. For other cut-offs, specifies
            the desired function.
- `sc`: Scaling of proposal distribution (only used if `adapt_cov=false`)
- `tol_targetAcceptance`: Target acceptance rate of tolerance adaptation.
- `rng`: Random number generator
- `dist`: Distance function (default Euclidean distance).
- `record_summaries`: Whether the summaries should be recorded.
- `out`: Output named tuple with fields `Theta` (matrix of d×n, each column
         is simulated parameter), `Dist` (vector of corresponding tolerances),
         `Summary` (matrix of corresponding summaries),
         `cutoff` (the cut-off function w/ tolerance), `Stats` (various
         statistics)
"""
function abc_mcmc(theta0::thetaT, prior::PriT, s_obs::summaryT,
    sim!::SimT, n; tol=zero(FT),
    b=convert(Int64, ceil(0.25n)),
    adapt_cov::Bool=true, adapt_tol::Bool=true,
    cutoff::CutoffT=nothing, sc=one(FT), stepSize_eta=FT(0.66),
    tol_targetAcceptance=FT(0.10), rng=Random.GLOBAL_RNG,
    record_summaries::Bool=true, dist=ABCdefaultDist) where {
        FT <: AbstractFloat,
        thetaT <: AbstractVector{FT},
        summaryT <: AbstractVector{FT},
        PriT <: Function, SimT <: Function,
        CutoffT <: Union{Nothing,String,Function,ABCCutoff}}

    d = length(theta0)
    rwm = AdaptiveMCMC.RWMState(theta0, rng)

    d_obs = length(s_obs)
    s = deepcopy(s_obs); s_ = deepcopy(s_obs)

    p_theta = prior(theta0)
    y_dist = FT(Inf)

    # Ensure that initial distance is finite & postive:
    while !isfinite(y_dist) || y_dist < 0.0
        sim!(s, theta0, rng)
        y_dist = dist(s, s_obs)
    end

    if tol == zero(FT)
        tol = y_dist + (y_dist == zero(FT))
    end
    if adapt_tol
        step = PolynomialStepSize(stepSize_eta)
        tolerance = AdaptiveMCMC.AdaptiveScalingMetropolis(tol_targetAcceptance, one(FT), step)
        tolerance.sc[] = -log(tol)
        tols_ = zeros(FT,b)
    else
        tols_ = nothing
    end
    if typeof(cutoff) == ABCCutoff
        phi = cutoff
    elseif typeof(cutoff) == String
        if cutoff=="gaussian"
            phi = ABCGaussCutoff(tol)
        elseif cutoff=="epa"
            phi = ABCEpaCutoff(tol)
        else
            error("Unknown cutoff")
        end
    else
        phi = ABCCutoff(FT(tol), cutoff)
    end

    if adapt_cov
        if adapt_tol
            proposal = AdaptiveMCMC.AdaptiveMetropolis(theta0, 2.38/sqrt(d), step)
        else
            proposal = AdaptiveMCMC.AdaptiveMetropolis(theta0)
        end
    else
        rwm = AdaptiveMCMC.RWMState(theta0, rng)
        proposal = sc
    end

    if record_summaries
        Summary = Array{FT}(undef,d_obs,n)
    else
        Summary = missing
    end

    proposal_burnin = proposal

    y_phi = phi(y_dist)
    if y_phi == 0.0
        # If fail to accept the initial observation, enforce acceptance
        y_phi = max(phi.fun(1.0), eps(Float64)) # phi.fun(1.0)
        y_dist = phi.tol
    end
    Theta = Array{FT}(undef,d,n)
    Dist = Vector{FT}(undef,n)
    all_acc = 0; acc = 0
    z = zeros(d)
    for k in 1:(n+b)
        AdaptiveMCMC.draw!(rwm, proposal)
        p_theta_ = prior(rwm.y)
        sim!(s_, rwm.y, rng)
        y_dist_ = dist(s_, s_obs)
        y_phi_ = phi(y_dist_)
        α = min(one(FT), exp(p_theta_ - p_theta)*y_phi_/y_phi)
        accept = (rand(rng) < α)
        if accept
            accept!(rwm)
            p_theta = p_theta_;
            y_dist = y_dist_; y_phi = y_phi_;
            s, s_ = s_, s
            all_acc += 1
            acc += (k>b)
        end
        if adapt_cov
            AdaptiveMCMC.adapt!(proposal, rwm, α, k)
        end
        if adapt_tol && k <= b
            AdaptiveMCMC.adapt!(tolerance, rwm, α, k)
            phi.tol = tols_[k] = exp(-tolerance.sc[])
        end
        if (k == b)
            proposal_burnin = deepcopy(proposal)
        end
        if (k > b)
            Theta[:,k-b] .= rwm.x
            Dist[k-b] = y_dist
            if record_summaries
                Summary[:,k-b] .= s
            end
        end
    end
    stats = (all_acc = all_acc/(n+b), acc = acc/n,
       rwm = rwm, adapt = proposal, adapt_burnin = proposal_burnin, tol = tols_)

    (Theta = Theta, Dist = Dist, Summary = Summary, s_obs = s_obs, cutoff = phi,
    Stats = stats, theta0 = theta0)
end


text:
#this example is modified from https://numericalenvironmental.wordpress.com/2018/05/26/a-steady-state-variably-saturated-flow-model-in-vertical-cross-section-a-finite-difference-approach-using-julia/
using Test
import Flux
import Zygote

include("utilities.jl")
include("inputdeck.jl")
include("ex.jl")

@Zygote.adjoint function solveforpsi(Ks)
	psi = solveforpsi(Ks)
	back = delta->begin
		lambda = transpose(f_psi(psi, Ks)) \ delta
		return (-transpose(f_Ks(psi, Ks)) * lambda, )
	end
	return psi, back
end

predict() = solveforpsi(Ks)
loss() = predict()[obsnode]
theta = Flux.params([Ks])
grads = Flux.gradient(() -> loss(), theta)
@test grads[Ks] ≈ grad

predict() = solveforpsi(Ks)
loss() = predict()[obsnode] ^ 2
theta = Flux.params([Ks])
grads = Flux.gradient(() -> loss(), theta)
@test grads[Ks] ≈ 2 * psi[obsnode] * grad

predict() = solveforpsi(exp.(log.(Ks)))
loss() = predict()[obsnode]
theta = Flux.params([Ks])
grads = Flux.gradient(() -> loss(), theta)
@test grads[Ks] ≈ grad

logKs = log.(Ks)
predict() = solveforpsi(exp.(logKs))
loss() = predict()[obsnode]
theta = Flux.params([logKs])
grads = Flux.gradient(() -> loss(), theta)
@test grads[logKs] ≈ grad .* Ks

logKs = log.(Ks)
predict() = solveforpsi(exp.(logKs))
loss() = predict()[obsnode] ^ 2
theta = Flux.params([logKs])
grads = Flux.gradient(() -> loss(), theta)
@test grads[logKs] ≈ 2 * psi[obsnode] * grad .* Ks


text:
# ---
# title: 719. Find K-th Smallest Pair Distance
# id: problem719
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Binary Search, Heap
# link: <https://leetcode.com/problems/find-k-th-smallest-pair-distance/description/>
# hidden: true
# ---
#
# Given an integer array, return the k-th smallest **distance** among all the
# pairs. The distance of a pair (A, B) is defined as the absolute difference
# between A and B.
#
# **Example 1:**
#
#
#
#     Input:
#     nums = [1,3,1]
#     k = 1
#     **Output: 0**
#     Explanation:
#     Here are all the pairs:
#     (1,3) -> 2
#     (1,1) -> 0
#     (3,1) -> 2
#     Then the 1st smallest distance pair is (1,1), and its distance is 0.
#
#
# **Note:**
#
#   1. `2 <= len(nums) <= 10000`.
#   2. `0 <= nums[i] < 1000000`.
#   3. `1 <= k <= len(nums) * (len(nums) - 1) / 2`.
#
#
## @lc code=start
using LeetCode

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


text:
#!/usr/bin/env julia

# PolyAAnalysis.jl
# =======
#
# A julia package for the representation and manipulation of biological sequences.
#
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

__precompile__()

module PolyAAnalysis

export
    annotate_polya_sites,
    annotate_clusters,
    check_polyA_prefixes,
    clust!,
    cluster!,
    detect_polyA_in_a_string,
    enumeratenames!,
    extend_poly_A,
    get_polyA_prefixes,
    get_split_key,
    get_transcripts_from_dict,
    get_transcripts_from_gff,
    getintervals,
    parseGFF3,
    parserecord,
    readbam,
    rmdups,
    stats_poly_a,
    stats_clusters,
    trim_polyA_3end,
    trim_polyA_file_records,
    trim_polyA_from_fastq_pair,
    trim_polyA_from_fastq_record,
    trim_polyA_from_fastq_pair_pararell,
    merge_adj_clusters!

import BioAlignments: BAM
import BioSequences: BioSymbols, @dna_str, FASTA, reverse_complement!, sequence
import BioSequences: FASTQ
import BufferedStreams
import CodecZlib
import CSV
import DataFrames: DataFrame, DataFrameRow
import DataFrames: eachrow, deleterows!
import DataStructures
import GenomicFeatures: eachoverlap, isoverlapping, strand, metadata, seqname, first, last
import GenomicFeatures: GFF3
import GenomicFeatures: Interval
import GenomicFeatures: IntervalCollection
import GenomicFeatures: Strand
import StringDistances: Levenshtein, evaluate
import TranscodingStreams
import FMIndexes: FMIndexes
import StatsBase: nquantile

include("AnnotatePolyA.jl")
include("ParseGFF.jl")
include("TrimmPolyA.jl")

end  # module PolyAAnalysis

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
import GR
using Images, Plots, Triplot

include("common.jl")

function to_rgba(x::UInt32)
    a = ((x & 0xff000000)>>24)/255
    b = ((x & 0x00ff0000)>>16)/255
    g = ((x & 0x0000ff00)>>8)/255
    r = (x & 0x000000ff)/255
    RGBA(r, g, b, a)
end

x,y,t = cartmesh(100)
z = peaks_unit.(x,y)

GR.setcolormap(1)
img = to_rgba.(Triplot.rasterize(x,y,z,t)')
plot([extrema(x)...],[extrema(y)...],img,yflip=false)


text:
@testset "homemade_bridging.jl" begin
    @testset "get_cp_variable(::VariableIndex)" begin
        opt = SeaPearl.Optimizer()

        x = SeaPearl.IntVar(1, 2, "x", opt.cpmodel.trailer)
        SeaPearl.addVariable!(opt.cpmodel, x)

        push!(opt.moimodel.variables, SeaPearl.MOIVariable("x", 1, 2, MOI.VariableIndex(1)))

        @test SeaPearl.get_cp_variable(opt, MOI.VariableIndex(1)) == x
    end
    @testset "get_cp_variable(::AffineIndex)" begin
        opt = SeaPearl.Optimizer()

        x = SeaPearl.IntVar(1, 2, "x", opt.cpmodel.trailer)
        y = SeaPearl.IntVar(1, 2, "y", opt.cpmodel.trailer)
        SeaPearl.addVariable!(opt.cpmodel, x)
        SeaPearl.addVariable!(opt.cpmodel, y)

        push!(opt.moimodel.variables, SeaPearl.MOIVariable("x", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("y", 1, 2, MOI.VariableIndex(2)))

        aff = SeaPearl.MOIAffineFunction("aff", MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5, -3], [MOI.VariableIndex(1), MOI.VariableIndex(2)]), 3))
        push!(opt.moimodel.affines, aff)

        # Dummy variable
        cp_aff = SeaPearl.IntVar(1, 2, "aff", opt.cpmodel.trailer)
        SeaPearl.addVariable!(opt.cpmodel, cp_aff)

        @test SeaPearl.get_cp_variable(opt, SeaPearl.AffineIndex(1)) == cp_aff
    end
    @testset "bridge_constraints!()" begin
        opt = SeaPearl.Optimizer()

        x = SeaPearl.IntVar(1, 2, "x", opt.cpmodel.trailer)
        y = SeaPearl.IntVar(1, 2, "y", opt.cpmodel.trailer)
        z = SeaPearl.IntVar(1, 2, "z", opt.cpmodel.trailer)
        SeaPearl.addVariable!(opt.cpmodel, x)
        SeaPearl.addVariable!(opt.cpmodel, y)
        SeaPearl.addVariable!(opt.cpmodel, z)
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("x", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("y", 1, 2, MOI.VariableIndex(2)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("z", 1, 2, MOI.VariableIndex(3)))

        ci1 = MOI.add_constraint(opt, MOI.VectorOfVariables([MOI.VariableIndex(1), MOI.VariableIndex(2)]), SeaPearl.NotEqualSet())
        ci2 = MOI.add_constraint(opt, MOI.VectorOfVariables([MOI.VariableIndex(3), MOI.VariableIndex(2)]), SeaPearl.NotEqualSet())

        SeaPearl.bridge_constraints!(opt)

        @test length(opt.cpmodel.constraints) == 2
        @test isa(opt.cpmodel.constraints[1], SeaPearl.NotEqual)
        @test SeaPearl.variablesArray(opt.cpmodel.constraints[1]) == [x, y]
        @test isa(opt.cpmodel.constraints[2], SeaPearl.NotEqual)
        @test SeaPearl.variablesArray(opt.cpmodel.constraints[2]) == [z, y]
    end
    @testset "bridge_variables!()" begin
        opt = SeaPearl.Optimizer()

        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(2)))

        SeaPearl.bridge_variables!(opt)

        @test length(keys(opt.cpmodel.variables)) == 2
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVar)
        @test length(opt.cpmodel.variables["1"].domain) == 2

        opt = SeaPearl.Optimizer()

        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, nothing, MOI.VariableIndex(2)))

        @test_throws AssertionError SeaPearl.bridge_variables!(opt)
    end
    @testset "build_affine_term!()" begin
        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        SeaPearl.bridge_variables!(opt)
        new_var = SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), -5.)

        @test length(keys(opt.cpmodel.variables)) == 3
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVarViewMul)
        @test isa(opt.cpmodel.variables["3"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["3"] == new_var
        @test new_var.x == opt.cpmodel.variables["2"]
        @test opt.cpmodel.variables["2"].a == 5
        @test opt.cpmodel.variables["2"].x == opt.cpmodel.variables["1"]

        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        SeaPearl.bridge_variables!(opt)
        new_var = SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), 5.)

        @test length(keys(opt.cpmodel.variables)) == 2
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["2"] == new_var
        @test new_var.a == 5
        @test opt.cpmodel.variables["2"].x == opt.cpmodel.variables["1"]

        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        SeaPearl.bridge_variables!(opt)
        new_var = SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), -1.)

        @test length(keys(opt.cpmodel.variables)) == 2
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["2"] == new_var
        @test opt.cpmodel.variables["2"].x == opt.cpmodel.variables["1"]

        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        SeaPearl.bridge_variables!(opt)
        new_var = SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), 1.)

        @test length(keys(opt.cpmodel.variables)) == 1
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test opt.cpmodel.variables["1"] == new_var

        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        SeaPearl.bridge_variables!(opt)
        @test_throws AssertionError SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), 0.)
        @test_throws AssertionError SeaPearl.build_affine_term!(opt, MOI.VariableIndex(1), -.2)
    end

    @testset "build_affine!()" begin
        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(2)))
        SeaPearl.bridge_variables!(opt)
        aff = SeaPearl.MOIAffineFunction(nothing, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [MOI.VariableIndex(1), MOI.VariableIndex(2)]), 3.))
        push!(opt.moimodel.affines, aff)

        @test SeaPearl.build_affine!(opt, aff) == "7"
        @test SeaPearl.get_cp_variable(opt, SeaPearl.AffineIndex(1)) == opt.cpmodel.variables["7"]

        @test length(keys(opt.cpmodel.variables)) == 8
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["3"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["3"].a == 5
        @test opt.cpmodel.variables["3"].x == opt.cpmodel.variables["1"]

        @test isa(opt.cpmodel.variables["4"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["4"].a == 3
        @test opt.cpmodel.variables["4"].x == opt.cpmodel.variables["2"]

        @test isa(opt.cpmodel.variables["5"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["5"].x == opt.cpmodel.variables["4"]

        @test isa(opt.cpmodel.variables["6"], SeaPearl.IntVarViewOffset)
        @test opt.cpmodel.variables["6"].c == 3
        @test opt.cpmodel.variables["6"].x == opt.cpmodel.variables["5"]

        @test isa(opt.cpmodel.variables["7"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["8"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["8"].x == opt.cpmodel.variables["7"]

        @test length(opt.cpmodel.constraints) == 1
        @test isa(opt.cpmodel.constraints[1], SeaPearl.SumToZero)
        @test opt.cpmodel.constraints[1].x == [
            opt.cpmodel.variables["3"],
            opt.cpmodel.variables["6"],
            opt.cpmodel.variables["8"]
        ]
    end

    @testset "bridge_affines!()" begin
        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(2)))
        SeaPearl.bridge_variables!(opt)
        aff = SeaPearl.MOIAffineFunction(nothing, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [MOI.VariableIndex(1), MOI.VariableIndex(2)]), 3.))
        push!(opt.moimodel.affines, aff)

        SeaPearl.bridge_affines!(opt)

        @test SeaPearl.get_cp_variable(opt, SeaPearl.AffineIndex(1)) == opt.cpmodel.variables["7"]

        @test length(keys(opt.cpmodel.variables)) == 8
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["3"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["3"].a == 5
        @test opt.cpmodel.variables["3"].x == opt.cpmodel.variables["1"]

        @test isa(opt.cpmodel.variables["4"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["4"].a == 3
        @test opt.cpmodel.variables["4"].x == opt.cpmodel.variables["2"]

        @test isa(opt.cpmodel.variables["5"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["5"].x == opt.cpmodel.variables["4"]

        @test isa(opt.cpmodel.variables["6"], SeaPearl.IntVarViewOffset)
        @test opt.cpmodel.variables["6"].c == 3
        @test opt.cpmodel.variables["6"].x == opt.cpmodel.variables["5"]

        @test isa(opt.cpmodel.variables["7"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["8"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["8"].x == opt.cpmodel.variables["7"]

        @test length(opt.cpmodel.constraints) == 1
        @test isa(opt.cpmodel.constraints[1], SeaPearl.SumToZero)
        @test opt.cpmodel.constraints[1].x == [
            opt.cpmodel.variables["3"],
            opt.cpmodel.variables["6"],
            opt.cpmodel.variables["8"]
        ]
    end

    @testset "fill_cpmodel!()" begin
        opt = SeaPearl.Optimizer()
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(1)))
        push!(opt.moimodel.variables, SeaPearl.MOIVariable("", 1, 2, MOI.VariableIndex(2)))
        aff = SeaPearl.MOIAffineFunction(nothing, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([5., -3.], [MOI.VariableIndex(1), MOI.VariableIndex(2)]), 3.))
        push!(opt.moimodel.affines, aff)

        ci1 = MOI.add_constraint(opt, MOI.VectorOfVariables([MOI.VariableIndex(1), MOI.VariableIndex(2)]), SeaPearl.NotEqualSet())

        SeaPearl.fill_cpmodel!(opt)

        @test SeaPearl.get_cp_variable(opt, SeaPearl.AffineIndex(1)) == opt.cpmodel.variables["7"]

        @test length(keys(opt.cpmodel.variables)) == 8
        @test isa(opt.cpmodel.variables["1"], SeaPearl.IntVar)
        @test isa(opt.cpmodel.variables["2"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["3"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["3"].a == 5
        @test opt.cpmodel.variables["3"].x == opt.cpmodel.variables["1"]

        @test isa(opt.cpmodel.variables["4"], SeaPearl.IntVarViewMul)
        @test opt.cpmodel.variables["4"].a == 3
        @test opt.cpmodel.variables["4"].x == opt.cpmodel.variables["2"]

        @test isa(opt.cpmodel.variables["5"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["5"].x == opt.cpmodel.variables["4"]

        @test isa(opt.cpmodel.variables["6"], SeaPearl.IntVarViewOffset)
        @test opt.cpmodel.variables["6"].c == 3
        @test opt.cpmodel.variables["6"].x == opt.cpmodel.variables["5"]

        @test isa(opt.cpmodel.variables["7"], SeaPearl.IntVar)

        @test isa(opt.cpmodel.variables["8"], SeaPearl.IntVarViewOpposite)
        @test opt.cpmodel.variables["8"].x == opt.cpmodel.variables["7"]

        @test length(opt.cpmodel.constraints) == 2
        @test isa(opt.cpmodel.constraints[2], SeaPearl.NotEqual)
        @test SeaPearl.variablesArray(opt.cpmodel.constraints[2]) == [opt.cpmodel.variables["1"], opt.cpmodel.variables["2"]]
        @test isa(opt.cpmodel.constraints[1], SeaPearl.SumToZero)
        @test opt.cpmodel.constraints[1].x == [
            opt.cpmodel.variables["3"],
            opt.cpmodel.variables["6"],
            opt.cpmodel.variables["8"]
        ]
    end

end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ d9097a02-51d0-11eb-304e-cb6b5940b4c7
begin
	using ReinforcementLearning
	using Flux

	const NX = 7
	const NY = 10
	const Wind = [CartesianIndex(w, 0) for w in [0, 0, 0, -1, -1, -1, -2, -2, -1, 0]]
	const StartPosition = CartesianIndex(4, 1)
	const Goal = CartesianIndex(4, 8)
	const ACTION = [
		CartesianIndex(0, -1),  # left
		CartesianIndex(0, 1),   # right
		CartesianIndex(-1, 0),  # up
		CartesianIndex(1, 0),   # down
	]

	const LinearInds = LinearIndices((NX, NY))

	Base.@kwdef mutable struct WindyGridWorldEnv <: AbstractEnv
		position::CartesianIndex{2} = StartPosition
	end

	RLBase.state_space(env::WindyGridWorldEnv) = Base.OneTo(length(LinearInds))
	RLBase.action_space(env::WindyGridWorldEnv) = Base.OneTo(length(ACTION))

	function (env::WindyGridWorldEnv)(a::Int)
		p = env.position + Wind[env.position[2]] + ACTION[a]
		p = CartesianIndex(min(max(p[1], 1), NX), min(max(p[2], 1), NY))
		env.position = p
		nothing
	end

	RLBase.state(env::WindyGridWorldEnv) = LinearInds[env.position]
	RLBase.is_terminated(env::WindyGridWorldEnv) = env.position == Goal
	RLBase.reward(env::WindyGridWorldEnv) = env.position == Goal ? 0.0 : -1.0

	RLBase.reset!(env::WindyGridWorldEnv) = env.position = StartPosition
end

# ╔═╡ 807e29b8-51d1-11eb-05bd-87574d1e6154
using Plots

# ╔═╡ feaaa9bc-51cf-11eb-1664-9f23f5cdc759
md"""
# Example 6.5: Windy Gridworld

First, let's define this environment by implementing the interfaces defined in `RLBase`.
"""

# ╔═╡ 412e68f4-51d1-11eb-2fe9-590a4f9d0cd3
world = WindyGridWorldEnv()

# ╔═╡ 63aec388-51d1-11eb-0d45-413471e91ea5
agent = Agent(
    policy=QBasedPolicy(
        learner=TDLearner(
            approximator=TabularQApproximator(
                ;n_state=length(state_space(world)),
                n_action=length(action_space(world)),
                opt=Descent(0.5)
                ),
            method=:SARSA
        ),
        explorer=EpsilonGreedyExplorer(0.1)
    ),
    trajectory=VectorSARTTrajectory()
)

# ╔═╡ 6e6a7f7e-51d1-11eb-296a-bb55d6d3a66c
hook = StepsPerEpisode()

# ╔═╡ 748e0e98-51d1-11eb-059f-5392e547fbfe
run(agent, world, StopAfterStep(8000),hook)

# ╔═╡ 840462b2-51d1-11eb-3bf4-f33179eed650
plot([i for (i, x) in enumerate(hook.steps) for _ in 1:x])

# ╔═╡ Cell order:
# ╟─feaaa9bc-51cf-11eb-1664-9f23f5cdc759
# ╠═d9097a02-51d0-11eb-304e-cb6b5940b4c7
# ╠═412e68f4-51d1-11eb-2fe9-590a4f9d0cd3
# ╠═63aec388-51d1-11eb-0d45-413471e91ea5
# ╠═6e6a7f7e-51d1-11eb-296a-bb55d6d3a66c
# ╠═748e0e98-51d1-11eb-059f-5392e547fbfe
# ╠═807e29b8-51d1-11eb-05bd-87574d1e6154
# ╠═840462b2-51d1-11eb-3bf4-f33179eed650