text:
rules = Dict()
open("day07/day07.in") do io
    global rules
    for line in readlines(io)
        line = replace(line, ("." => ""))
        line = split(line, "contain ")
        key_bag = split(line[1], " ")[1] * " " * split(line[1], " ")[2]
        bags = split(line[2], ", ")
        tmpDict = Dict()
        for bag in bags
            if bag != "no other bags"
                bag_type = split(bag, " ")[2] * " " * split(bag, " ")[3]
                tmpDict[bag_type] = parse(Int64, split(bag, " ")[1])
            end
        end
        rules[key_bag] = tmpDict
    end
end

function num_of_bags(bag_type::String)::Int64
    global rules
    num = 1
    if get(rules, bag_type, 0) != 0
        bags = rules[bag_type]
        for (bag, mult) in bags
            num += mult * num_of_bags(bag)
        end
    end
    return num
end

println(num_of_bags("shiny gold")-1)

ext: jl

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using Yao, YaoExtensions
using YaoBlocks.Optimise: replace_block
using YaoBlocks.ConstGate
using YaoBlocks
using Test
using LinearAlgebra

"""
    TODO:
    * optimize to optimal value,
    * if this optimal value's loss is close to `θ = 0.0`, kill it.
"""
performance(params, gradient) = abs.(params .* gradient)

nbit = 5
c = variational_circuit(nbit, 50)

dispatch!(c, :random)
h = heisenberg(nbit)
cc = replace_block(x->x isa RotationGate ? Bag(x) : x, c)

function train(h, c; kill_rate=0.1, kill_step=10, niter=100)
    bags = blockfilter(b->b isa Bag && isenabled(b), c)
    for i=1:niter
        regδ, paramsδ = expect'(h, zero_state(nbit) => c)
        dispatch!(-, c, 0.01 .* paramsδ)
        loss = expect(h, zero_state(nbit)=>c) |> real

        Nk = round(Int, kill_rate*length(paramsδ))

        if i%kill_step == 0 && kill_rate != 0
            ps = performance(parameters(c), paramsδ)
            cut = sort(ps)[Nk]
            disable_block!.(bags[ps .<= cut])
            bags = bags[ps .> cut]
        end
        @show i,nparameters(c), loss
    end
end

train(h, cc; kill_rate=0.01, niter=200)

cr = variational_circuit(nbit, 5)
dispatch!(cr, :random)
train(h, cr; kill_rate=0.0, niter=200)


text:
using Test

@testset "ITensors.jl" begin
  @testset "$filename" for filename in [
    "tagset.jl",
    "smallstring.jl",
    "symmetrystyle.jl",
    "index.jl",
    "indexset.jl",
    "not.jl",
    "itensor.jl",
    "itensor_slice.jl",
    "itensor_scalar_contract.jl",
    "itensor_combine_contract.jl",
    "itensor_contraction_sequence_optimization.jl",
    "broadcast.jl",
    "diagitensor.jl",
    "contract.jl",
    "combiner.jl",
    "debug_checks.jl",
    "trg.jl",
    "ctmrg.jl",
    "iterativesolvers.jl",
    "dmrg.jl",
    "sitetype.jl",
    "phys_site_types.jl",
    "decomp.jl",
    "lattices.jl",
    "mps.jl",
    "mpo.jl",
    "sweeps.jl",
    "sweepnext.jl",
    "autompo.jl",
    "svd.jl",
    "qn.jl",
    "qnindex.jl",
    "qnitensor.jl",
    "qncombiner.jl",
    "qndiagitensor.jl",
    "empty.jl",
    "qnmpo.jl",
    "readwrite.jl",
    "readme.jl",
    "examples.jl",
  ]
    println("Running $filename")
    include(filename)
  end
end


text:
# mtd_inference.jl

export ParamsMTD, PriorMTD, ModMTD,
  sim_mtd, symmetricDirPrior_mtd, transTensor_mtd,
  counttrans_mtd,
  mcmc!, timemod!, etr, llik_MTD;

mutable struct ParamsMTD
  lλ::Vector{Float64}
  ζ::Vector{Int} # will be length TT - L
  lQ::Matrix{Float64} # organized so first index is now and lag 1 is the next index

  ParamsMTD(lλ, ζ, lQ) = new(deepcopy(lλ), deepcopy(ζ), deepcopy(lQ))
end

mutable struct PriorMTD
  λ::Union{Vector{Float64}, SparseDirMix, SBMprior}
  Q::Union{Matrix{Float64}, Vector{SparseDirMix}, Vector{SBMprior}}

  PriorMTD(λ, Q) = new(deepcopy(λ), deepcopy(Q))
end

mutable struct ModMTD
  L::Int # maximal order
  K::Int # number of states
  TT::Int
  S::Vector{Int}
  prior::PriorMTD
  state::ParamsMTD
  iter::Int

  ModMTD(L, K, TT, S, prior, state) = new(L, K, TT, deepcopy(S), deepcopy(prior), deepcopy(state), 0)
end

"""
    sim_mtd(TT, nburn, L, K, λ, Q)
"""
function sim_mtd(TT::Int, nburn::Int, L::Int, K::Int,
  λ::Vector{Float64}, Q::Matrix{Float64})

  Nsim = nburn + TT

  ζ = [ StatsBase.sample(Weights(λ)) for i in 1:(Nsim-L) ]

  S = Vector{Int}(undef, Nsim)
  S[1:L] = StatsBase.sample(1:K, L)

  for tt in (L+1):(Nsim)
    i = tt - L
    Slagrev_now = S[range(tt-1, step=-1, length=L)]
    pvec = deepcopy( Q[:, Slagrev_now[ζ[i]] ] )
    S[tt] = StatsBase.sample(Weights( pvec ))
  end

  S[(nburn+1):(Nsim)], ζ[(nburn+1):(Nsim-L)]
end

"""
    symmetricPrior_mtd(size_λ, size_Q, L, K)
"""
function symmetricDirPrior_mtd(size_λ::Float64, size_Q::Float64,
  L::Int, K::Int)

  α0_λ = fill( size_λ / float(L), L )
  a0_Q = size_Q / float(K)
  α0_Q = fill( a0_Q, (K, K) )

  (α0_λ, α0_Q)
end

"""
    transTensor_mtd(L, K, λ, Q)

Calculate full transition tensor from λ and Q.
"""
function transTensor_mtd(L::Int, K::Int, λ::Vector{Float64}, Q::Matrix{Float64})

    froms, nfroms = create_froms(K, L) # in this case, ordered by now, lag1, lag2, etc.
    Ωmat = zeros(Float64, (K, nfroms))
    for i in 1:nfroms
        for k in 1:K
            for ℓ in 1:L
                Ωmat[k,i] += λ[ℓ] .* Q[k,froms[i][ℓ]]
            end
        end
    end
    reshape(Ωmat, fill(K,L+1)...)
end

"""
    rpost_lλ_mtd(α0_λ, ζ, L)
"""
function rpost_lλ_mtd(α0_λ::Vector{Float64}, ζ::Vector{Int}, L::Int)

    Nζ = StatsBase.counts(ζ, 1:L)
    α1_λ = α0_λ .+ Nζ
    lλ_out = SparseProbVec.rDirichlet(α1_λ, logout=true)

    lλ_out
end
function rpost_lλ_mtd(prior::SparseDirMix, ζ::Vector{Int}, L::Int)

    Nζ = StatsBase.counts(ζ, 1:L)
    α1_λ = prior.α .+ Nζ
    d = SparseProbVec.SparseDirMix(α1_λ, prior.β)
    lλ_out = SparseProbVec.rand(d, logout=true)

    lλ_out
end
function rpost_lλ_mtd(prior::SBMprior, ζ::Vector{Int}, L::Int)

    Nζ = StatsBase.counts(ζ, 1:L)
    post_lλ = SparseProbVec.SBM_multinom_post(prior, Nζ)
    lλ_out = SparseProbVec.rand(post_lλ, logout=true)

    lλ_out
end

"""
    counttrans_mtd(S, TT, ζ, L, K)

    ### Example
    ```julia
    L = 2
    K = 3
    TT = 12
    S = [1,2,1,3,3,1,2,1,3,2,1,1]
    ζ = [1,2,1,2,1,2,1,2,1,2]
      counttrans_mtd(S, TT, ζ, L, K)
    ```
"""
function counttrans_mtd(S::Vector{Int}, TT::Int, ζ::Vector{Int},
  L::Int, K::Int)

  ## initialize
  N_out = zeros(Int, (K,K))

  ## pass through data and add counts
  for tt in (L+1):(TT)
    Slagrev_now = deepcopy( S[range(tt-1, step=-1, length=L)] )
    from = deepcopy( Slagrev_now[ ζ[tt-L] ] )
    N_out[ S[tt], from ] += 1
  end

  N_out
end

"""
    rpost_lQ_mtd(S, TT, prior, ζ, L, K)
"""
function rpost_lQ_mtd(S::Vector{Int}, TT::Int, prior::Matrix{Float64},
    ζ::Vector{Int}, L::Int, K::Int)

    ## initialize
    α0_Q = deepcopy(prior)
    lQ_out = Matrix{Float64}(undef, K, K)

    N = counttrans_mtd(S, TT, ζ, L, K)

    α1_Q = α0_Q .+ N
    for j in 1:K
        lQ_out[:,j] = SparseProbVec.rDirichlet(α1_Q[:,j], logout=true)
    end

    lQ_out
end
function rpost_lQ_mtd(S::Vector{Int}, TT::Int,
    prior::Vector{SparseDirMix},
    ζ::Vector{Int}, L::Int, K::Int)

    ## initialize
    lQ_out = Matrix{Float64}(undef, K, K)

    N = counttrans_mtd(S, TT, ζ, L, K)

    for j in 1:K
        α1 = prior[j].α .+ N[:,j]
        d = SparseProbVec.SparseDirMix(α1, prior[j].β)
        lQ_out[:,j] = SparseProbVec.rand(d, logout=true)
    end

    lQ_out
end
function rpost_lQ_mtd(S::Vector{Int}, TT::Int,
    prior::Vector{SBMprior}, ζ::Vector{Int},
    L::Int, K::Int)

    ## initialize
    lQ_out = Matrix{Float64}(undef, K, K)

    N = counttrans_mtd(S, TT, ζ, L, K)

    for j in 1:K
        d = SparseProbVec.SBM_multinom_post(prior[j], N[:,j])
        lQ_out[:,j] = SparseProbVec.rand(d, logout=true)
    end

    lQ_out
end

"""
    rpost_ζ_mtd(S, TT, lλ, lQ, L, K)
"""
function rpost_ζ_mtd(S::Vector{Int}, TT::Int,
  lλ::Vector{Float64}, lQ::Array{Float64,2},
  L::Int, K::Int)

  ζ_out = Vector{Int}(undef, TT-L)

  for i in 1:(TT-L)
      tt = i + L
      Slagrev_now = S[range(tt-1, step=-1, length=L)]
      lp = Vector{Float64}(undef, L)
      for j in 1:L
          lp[j] = lλ[j] + lQ[ append!([deepcopy(S[tt])], deepcopy(Slagrev_now[j]))... ]
      end

      w = exp.( lp .- maximum(lp) )
      ζ_out[i] = StatsBase.sample(Weights(w))
  end

  ζ_out
end

"""
    rpost_ζ_mtd_marg(S, ζ_old, lλ, prior_Q, TT, L, K)

Collapsed conditinal updates for ζ marginalizing over Q
"""
function rpost_ζ_mtd_marg(S::Vector{Int}, ζ_old::Vector{Int},
    prior_Q::Matrix{Float64},
    lλ::Vector{Float64},
    TT::Int, L::Int, K::Int)

  α0_Q = deepcopy(prior_Q)
  ζ_out = deepcopy(ζ_old)
  N_now = counttrans_mtd(S, TT, ζ_old, L, K) # rows are tos, cols are froms

  for i in 1:(TT-L)  # i indexes ζ, tt indexes S
    tt = i + L
    Slagrev_now = S[range(tt-1, step=-1, length=L)]
    N0 = deepcopy(N_now)
    N0[ S[tt], Slagrev_now[ ζ_out[i] ] ] -= 1
    α1_Q = α0_Q .+ N0
    eSt = [1.0*(ii==S[tt]) for ii in 1:K]

    kuse = unique(Slagrev_now)
    nkuse = length(kuse)

    lmvbn0 = [ SparseProbVec.lmvbeta( α1_Q[:,kk] ) for kk in kuse ]
    lmvbn1 = [ SparseProbVec.lmvbeta( α1_Q[:,kk] + eSt ) for kk in kuse ]

    lw = zeros(Float64, L)

    for ℓ in 1:L
        lw[ℓ] = deepcopy(lλ[ℓ])
        for kk in 1:nkuse
            if Slagrev_now[ℓ] == kuse[kk]
                lw[ℓ] += lmvbn1[kk]
            else
                lw[ℓ] += lmvbn0[kk]
            end
        end
    end

    w = exp.( lw .- maximum(lw) )
    ζ_out[i] = StatsBase.sample(Weights( w ))
    N_now = deepcopy(N0)
    N_now[ S[tt], Slagrev_now[ ζ_out[i] ] ] += 1

  end

  ζ_out
end
function rpost_ζ_mtd_marg(S::Vector{Int}, ζ_old::Vector{Int},
    prior_Q::Union{Vector{SparseDirMix}, Vector{SBMprior}},
    lλ::Vector{Float64},
    TT::Int, L::Int, K::Int)

  ζ_out = deepcopy(ζ_old)
  N_now = counttrans_mtd(S, TT, ζ_old, L, K) # rows are tos, cols are froms

  for i in 1:(TT-L)  # i indexes ζ, tt indexes S
    tt = i + L
    Slagrev_now = S[range(tt-1, step=-1, length=L)]
    N0 = deepcopy(N_now)
    N0[ S[tt], Slagrev_now[ ζ_out[i] ] ] -= 1
    eSt = [1*(ii==S[tt]) for ii in 1:K]

    kuse = unique(Slagrev_now)
    nkuse = length(kuse)

    lmarg0 = [ SparseProbVec.logMarginal(prior_Q[kk], N0[:,kk]) for kk in kuse ]
    lmarg1 = [ SparseProbVec.logMarginal(prior_Q[kk], N0[:,kk] + eSt) for kk in kuse ]

    lw = zeros(Float64, L)

    for ℓ in 1:L
        lw[ℓ] = deepcopy(lλ[ℓ])
        for kk in 1:nkuse
            if Slagrev_now[ℓ] == kuse[kk]
                lw[ℓ] += lmarg1[kk]
            else
                lw[ℓ] += lmarg0[kk]
            end
        end
    end

    w = exp.( lw .- maximum(lw) )
    ζ_out[i] = StatsBase.sample(Weights( w ))
    N_now = deepcopy(N0)
    N_now[ S[tt], Slagrev_now[ ζ_out[i] ] ] += 1

  end

  ζ_out
end

"""
    MetropIndep_λζ(S::Vector{Int}, lλ_old::Vector{Float64}, ζ_old::,
        prior_λ, prior_Q,
        TT::Int, L::Int, K::Int)

    Independence Metropolis step for λ and ζ.
"""
function MetropIndep_λζ(S::Vector{Int}, lλ_old::Vector{Float64}, ζ_old::Vector{Int},
    prior_λ::Union{Vector{Float64}, SparseDirMix, SBMprior},
    prior_Q::Union{Matrix{Float64}, Vector{SparseDirMix}, Vector{SBMprior}},
    TT::Int, L::Int, K::Int)

  if typeof(prior_λ) == Vector{Float64}
      lλ_cand = SparseProbVec.rDirichlet(prior_λ, logout=true)
  elseif typeof(prior_λ) == SBMprior
      lλ_cand = SparseProbVec.rand(prior_λ, logout=true)
  elseif typeof(prior_λ) == SparseDirMix
      lλ_cand = SparseProbVec.rand(prior_λ, logout=true)
  end

  λ_cand = exp.(lλ_cand)
  ζ_cand = [ StatsBase.sample( Weights(λ_cand) ) for i in 1:(TT-L) ]

  N_cand = counttrans_mtd(S, TT, ζ_cand, L, K) # rows are tos, cols are froms
  N_old = counttrans_mtd(S, TT, ζ_old, L, K) # rows are tos, cols are froms

  if typeof(prior_Q) == Matrix{Float64}
      lmargQ_cand = [ SparseProbVec.lmvbeta(N_cand[:,kk] .+ prior_Q[:,kk]) for kk in 1:K ] # can ignore denominator
      lmargQ_old = [ SparseProbVec.lmvbeta(N_old[:,kk] .+ prior_Q[:,kk]) for kk in 1:K ]
  else
      lmargQ_cand = [ SparseProbVec.logMarginal(prior_Q[kk], N_cand[:,kk]) for kk = 1:K ]
      lmargQ_old = [ SparseProbVec.logMarginal(prior_Q[kk], N_old[:,kk]) for kk = 1:K ]
  end

  ll_cand = sum( lmargQ_cand )
  ll_old = sum( lmargQ_old )

  lu = log(rand())
  if lu < (ll_cand - ll_old)
      lλ_out = lλ_cand
      ζ_out = ζ_cand
  else
      lλ_out = lλ_old
      ζ_out = ζ_old
  end

  (lλ_out, ζ_out)
end

"""
mcmc!(model::ModMTD, n_keep::Int, save::Bool=true,
    report_filename::String="out_progress.txt", thin::Int=1, jmpstart_iter::Int=25,
    report_freq::Int=1000;
    monitor::Vector{Symbol}=[:lλ, :lQ])
"""
function mcmc!(model::ModMTD, n_keep::Int; save::Bool=true,
    report_filename::String="out_progress.txt", thin::Int=1, jmpstart_iter::Int=25,
    report_freq::Int=1000,
    monitor::Vector{Symbol}=[:lλ, :lQ])

    ## output files
    report_file = open(report_filename, "a+")
    write(report_file, "Commencing MCMC at $(Dates.now()) at iteration $(model.iter) for $(n_keep * thin) iterations.\n")
    close(report_file)

    if save
        sims = postSimsInit(monitor, n_keep, model.state)
    end

    ## sampling
    for i in 1:n_keep
        for j in 1:thin

            jmpstart = (model.iter % jmpstart_iter == 0)

            if jmpstart

                model.state.lλ, model.state.ζ = MetropIndep_λζ(model.S,
                    model.state.lλ, model.state.ζ, model.prior.λ,
                    model.prior.Q,
                    model.TT, model.L, model.K)

            else

                model.state.ζ = rpost_ζ_mtd_marg(model.S, model.state.ζ,
                    model.prior.Q, model.state.lλ,
                    model.TT, model.L, model.K)

                model.state.lλ = rpost_lλ_mtd(model.prior.λ, model.state.ζ, model.L)

            end

            model.state.lQ = rpost_lQ_mtd(model.S, model.TT, model.prior.Q,
                    model.state.ζ, model.L, model.K)

            model.iter += 1
            if model.iter % report_freq == 0
                report_file = open(report_filename, "a+")
                write(report_file, "Iter $(model.iter) at $(Dates.now())\n")
                llik_now = llik_MTD(model.S, model.state.lλ, model.state.lQ)
                write(report_file, "Log-likelihood $(llik_now)\n")
                write(report_file, "Current lag weights: $( round.(exp.(model.state.lλ), digits=3) )\n\n")
                close(report_file)
            end
        end

        if save
            for field in monitor
              sims[i][field] = deepcopy(getfield(model.state, field))
            end
            sims[i][:llik] = llik_MTD(model.S, model.state.lλ, model.state.lQ)
        end
    end

    if save
        return sims
    else
        return model.iter
    end

end

"""
    llik_MTD(S::Vector{Int}, lλ::Vector{Vector{Float64}}, lQ::Matrix{Float64})

Computes the (conditional on intial values) log-likelihood for a time series.

### Example
```julia
    lλ = log.([0.0, 0.0, 1.0])
    lQ = log.(reshape([0.7, 0.3, 0.25, 0.75], (2,2)))
    S = [1,1,2,1,1,2,2,1,2,2,1,2,1,2,1,1,1,1,2,1,2]
    llik_MTD(S, lλ, lQ)
```
"""
function llik_MTD(S::Vector{Int}, lλ::Vector{Float64}, lQ::Matrix{Float64})

    TT = length(S)
    L = length(lλ)

    n = TT - L

    lpt = Matrix{Float64}(undef, n, L)

    for i in 1:n
        tt = i + L
        Slagrev_now = S[range(tt-1, step=-1, length=L)]

        for ℓ in 1:L
            lpt[i,ℓ] = lλ[ℓ] + lQ[S[tt], Slagrev_now[ℓ]]
        end
    end

    lptvec = [ SparseProbVec.logsumexp(lpt[i,:]) for i in 1:n ]
    llik = sum(lptvec)

    return llik
end

max_stars_repo_name: mheiner/MTD

max_issues_repo_name: mheiner/MTD

max_forks_repo_name: mheiner/MTD


text:
import FastTransforms: chebyshevpoints#, transform, itransform
import OffsetArrays: OffsetArray

@inline function _chebyshev(::Val{1}, j::Int, x::T) where {T}
    if j < 0
        return zero(T)
    elseif j == 0
        return one(T)
    elseif j == 1
        return x
    else
        return _chebyshev(Val(1), j-1, x) * 2x - _chebyshev(Val(1), j-2, x)
    end
end

@inline function _chebyshev(::Val{2}, j::Int, x::T) where {T}
    if j < 0
        return zero(T)
    elseif j == 0
        return one(T)
    elseif j == 1
        return 2x
    else
        return _chebyshev(Val(2), j-1, x) * 2x - _chebyshev(Val(2), j-2, x)
    end
end

"""
Chebyshev basis on the interval [-1..+1].
"""
struct Chebyshev{kind, T, BT, XT <: AbstractVector{T}} <: Basis{T}
    b::BT
    x::XT

    function Chebyshev{kind, T}(n::Int) where {kind, T}
        p = n-1
        x = chebyshevpoints(T, n, Val(kind))
        b = OffsetArray([y -> _chebyshev(Val(kind), i, y) for i in 0:p], 0:p)
        new{kind, T, typeof(b), typeof(x)}(b, x)
    end

    Chebyshev{kind}(::Type{T}, n::Int) where {kind, T} = Chebyshev{kind, T}(n)
    Chebyshev{kind}(n::Int) where {kind} = Chebyshev{kind, Float64}(n)
end

Chebyshev(::Type{T}, n::Int, ::Val{kind}) where {kind, T} = Chebyshev{kind, T}(n)
Chebyshev(n::Int, ::Val{kind}) where {kind} = Chebyshev(Float64, n, Val(kind))

const ChebyshevT = Chebyshev{1}
const ChebyshevU = Chebyshev{2}

(C::Chebyshev)(x::Number, j::Int) = C.b[j](x)

basis(C::Chebyshev) = C.b
nodes(C::Chebyshev) = C.x
nbasis(C::Chebyshev) = length(basis(C))
nnodes(C::Chebyshev) = length(nodes(C))
order(C::Chebyshev)  = nnodes(C)
degree(C::Chebyshev) = nnodes(C) - 1

Base.eltype(::Chebyshev{kind,T}) where {kind,T} = T
Base.eachindex(C::Chebyshev) = eachindex(C.b)
Base.axes(C::Chebyshev) = (Inclusion(0..1), eachindex(C))
ContinuumArrays.grid(C::Chebyshev) = nodes(C)

Base.hash(C::Chebyshev{kind}, h::UInt) where {kind} = hash(C.x, hash(kind, h))
Base.:(==)(C1::Chebyshev{kind1}, C2::Chebyshev{kind2}) where {kind1,kind2} = (C1.x == C2.x && kind1 == kind2)
Base.isequal(C1::Chebyshev{kind1,T1}, C2::Chebyshev{kind2,T2}) where {kind1,T1,kind2,T2} = (T1 == T2 && C1 == C2)
Base.isapprox(C1::Chebyshev, C2::Chebyshev; kwargs...) = isapprox(C1.x, C2.x; kwargs...)

Base.getindex(C::Chebyshev, x::Number, j::Int) = C(x,j)
Base.getindex(C::Chebyshev, x::Number,  ::Colon) = [b(x) for b in C.b]
Base.getindex(C::Chebyshev, X::AbstractVector, j::Int) = C.(X,j)
Base.getindex(C::Chebyshev, X::AbstractVector,  ::Colon) = [b(x) for x in X, b in C.b]

## Derivative

@simplify *(D::Derivative, C::Chebyshev) = Mul(D,C)

const ChebyshevDerivative  = QMul2{<:Derivative,<:Chebyshev}
const ChebyshevTDerivative = QMul2{<:Derivative,<:ChebyshevT}
const ChebyshevUDerivative = QMul2{<:Derivative,<:ChebyshevU}

function _eval(D::ChebyshevTDerivative, x::DT, i::Int) where {DT}
    local C = D.B
    local x̃ = promote_type(eltype(C), DT)(x)
    @assert i ≥ 0 && i < nbasis(C)
    _chebyshev(Val(2), i-1, x̃) * i
end

function _eval(D::ChebyshevUDerivative, x::DT, i::Int) where {DT}
    local C = D.B
    local x̃ = promote_type(eltype(C), DT)(x)
    @assert i ≥ 0 && i < nbasis(C)
    ( _chebyshev(Val(1), i+1, x̃) * (i+1) - _chebyshev(Val(2), i, x̃) * x̃ ) / (x̃^2 - 1)
end

Base.getindex(D::ChebyshevDerivative, x::Number, j::Integer) = _eval(D, x, j)
Base.getindex(D::ChebyshevDerivative, x::Number,  ::Colon) = [_eval(D, x, j) for j in eachindex(D.B)]
Base.getindex(D::ChebyshevDerivative, X::AbstractVector, j::Integer) = [_eval(D, x, j) for x in X]
Base.getindex(D::ChebyshevDerivative, X::AbstractVector,  ::Colon) = [_eval(D, x, j) for x in X, j in eachindex(D.B)]

Base.adjoint(C::Chebyshev) = Derivative(axes(C,1)) * C


text:
using Documenter, DiffusionDefinition

makedocs(;
    modules=[DiffusionDefinition],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                )
            )
        )
    ),
    pages=[
        "Home" => "index.md",
        "Overview" => Any[
            "Defining diffusions" => joinpath("overview", "defining_diffusions.md"),
            "Loading diffusions" => joinpath("overview", "load_diff.md"),
            "Sampling trajectories" => joinpath("overview", "sampling.md"),
        ],
        "Examples" => Any[
            "Lotka-Volterra system" => joinpath("examples", "lotka_volterra.md"),
            "Sine diffusion" => joinpath("examples", "sine.md"),
            "SIR model" => joinpath("examples", "sir.md"),
            "Lorenz63 system" => joinpath("examples", "lorenz63.md"),
            "Lorenz96 system" => joinpath("examples", "lorenz96.md"),
            "Favetto-Samson model" => joinpath("examples", "favetto_samson.md"),
            "Prokaryotic autoregulatory gene network" => joinpath("examples", "prokaryote.md"),
            "FitzHugh-Nagumo model" => joinpath("examples", "fitzhugh_nagumo.md"),
            "Jansen-Rit model" => joinpath("examples", "jansen_rit.md"),
        ],
        "Index" => "post_index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/DiffusionDefinition.jl/blob/{commit}{path}#L{line}",
    sitename="DiffusionDefinition.jl",
    authors="Marcin Mider",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/DiffusionDefinition.jl",
)


text:
# For testing Tangent
struct Foo
    x
    y::Float64
end

# For testing Primal + Tangent performance
struct Bar
    x::Float64
end

# For testing Tangent: it is an invarient of the type that x2 = 2x
# so simple addition can not be defined
struct StructWithInvariant
    x
    x2

    StructWithInvariant(x) = new(x, 2x)
end

@testset "Tangent" begin
    @testset "empty types" begin
        @test typeof(Tangent{Tuple{}}()) == Tangent{Tuple{}, Tuple{}}
    end
    @testset "convert" begin
        @test convert(NamedTuple, Tangent{Foo}(x=2.5)) == (; x=2.5)
        @test convert(Tuple, Tangent{Tuple{Float64,}}(2.0)) == (2.0,)
        @test convert(Dict, Tangent{Dict}(Dict(4 => 3))) == Dict(4 => 3)
    end

    @testset "==" begin
        @test Tangent{Foo}(x=0.1, y=2.5) == Tangent{Foo}(x=0.1, y=2.5)
        @test Tangent{Foo}(x=0.1, y=2.5) == Tangent{Foo}(y=2.5, x=0.1)
        @test Tangent{Foo}(y=2.5, x=ZeroTangent()) == Tangent{Foo}(y=2.5)

        @test Tangent{Tuple{Float64,}}(2.0) == Tangent{Tuple{Float64,}}(2.0)
        @test Tangent{Dict}(Dict(4 => 3)) == Tangent{Dict}(Dict(4 => 3))

        tup = (1.0, 2.0)
        @test Tangent{typeof(tup)}(1.0, 2.0) == Tangent{typeof(tup)}(1.0, @thunk(2*1.0))
        @test Tangent{typeof(tup)}(1.0, 2.0) == Tangent{typeof(tup)}(1.0, 2)

        @test Tangent{Foo}(;y=2.0,) == Tangent{Foo}(;x=ZeroTangent(), y=Float32(2.0),)
    end

    @testset "hash" begin
        @test hash(Tangent{Foo}(x=0.1, y=2.5)) == hash(Tangent{Foo}(y=2.5, x=0.1))
        @test hash(Tangent{Foo}(y=2.5, x=ZeroTangent())) == hash(Tangent{Foo}(y=2.5))
    end

    @testset "indexing, iterating, and properties" begin
        @test keys(Tangent{Foo}(x=2.5)) == (:x,)
        @test propertynames(Tangent{Foo}(x=2.5)) == (:x,)
        @test haskey(Tangent{Foo}(x=2.5), :x) == true
        if isdefined(Base, :hasproperty)
            @test hasproperty(Tangent{Foo}(x=2.5), :y) == false
        end
        @test Tangent{Foo}(x=2.5).x == 2.5

        @test keys(Tangent{Tuple{Float64,}}(2.0)) == Base.OneTo(1)
        @test propertynames(Tangent{Tuple{Float64,}}(2.0)) == (1,)
        @test getproperty(Tangent{Tuple{Float64,}}(2.0), 1) == 2.0
        @test getproperty(Tangent{Tuple{Float64,}}(@thunk 2.0^2), 1) == 4.0
        @test getproperty(Tangent{Tuple{Float64,}}(a=(@thunk 2.0^2),), :a) == 4.0

        # TODO: uncomment this once https://github.com/JuliaLang/julia/issues/35516
        @test_broken haskey(Tangent{Tuple{Float64}}(2.0), 1) == true
        @test_broken hasproperty(Tangent{Tuple{Float64}}(2.0), 2) == false

        @test length(Tangent{Foo}(x=2.5)) == 1
        @test length(Tangent{Tuple{Float64,}}(2.0)) == 1

        @test eltype(Tangent{Foo}(x=2.5)) == Float64
        @test eltype(Tangent{Tuple{Float64,}}(2.0)) == Float64

        # Testing iterate via collect
        @test collect(Tangent{Foo}(x=2.5)) == [2.5]
        @test collect(Tangent{Tuple{Float64,}}(2.0)) == [2.0]

        # Test indexed_iterate
        ctup = Tangent{Tuple{Float64,Int64}}(2.0, 3)
        _unpack2tuple = function(comp)
            a, b = comp
            return (a, b)
        end
        @inferred _unpack2tuple(ctup)
        @test _unpack2tuple(ctup) === (2.0, 3)

        # Test getproperty is inferrable
        _unpacknamedtuple = comp -> (comp.x, comp.y)
        if VERSION ≥ v"1.2"
            @inferred _unpacknamedtuple(Tangent{Foo}(x=2, y=3.0))
            @inferred _unpacknamedtuple(Tangent{Foo}(y=3.0))
        end
    end

    @testset "reverse" begin
        c = Tangent{Tuple{Int, Int, String}}(1, 2, "something")
        cr = Tangent{Tuple{String, Int, Int}}("something", 2, 1)
        @test reverse(c) === cr

        # can't reverse a named tuple or a dict
        @test_throws MethodError reverse(Tangent{Foo}(;x=1.0, y=2.0))

        d = Dict(:x => 1, :y => 2.0)
        cdict = Tangent{Foo, typeof(d)}(d)
        @test_throws MethodError reverse(Tangent{Foo}())
    end

    @testset "unset properties" begin
        @test Tangent{Foo}(; x=1.4).y === ZeroTangent()
    end

    @testset "conj" begin
        @test conj(Tangent{Foo}(x=2.0+3.0im)) == Tangent{Foo}(x=2.0-3.0im)
        @test ==(
            conj(Tangent{Tuple{Float64,}}(2.0+3.0im)),
            Tangent{Tuple{Float64,}}(2.0-3.0im)
        )
        @test ==(
            conj(Tangent{Dict}(Dict(4 => 2.0 + 3.0im))),
            Tangent{Dict}(Dict(4 => 2.0 + -3.0im)),
        )
    end

    @testset "canonicalize" begin
        # Testing iterate via collect
        @test ==(
            canonicalize(Tangent{Tuple{Float64,}}(2.0)),
            Tangent{Tuple{Float64,}}(2.0)
        )

        @test ==(
            canonicalize(Tangent{Dict}(Dict(4 => 3))),
            Tangent{Dict}(Dict(4 => 3)),
        )

        # For structure it needs to match order and ZeroTangent() fill to match primal
        CFoo = Tangent{Foo}
        @test canonicalize(CFoo(x=2.5, y=10)) == CFoo(x=2.5, y=10)
        @test canonicalize(CFoo(y=10, x=2.5)) == CFoo(x=2.5, y=10)
        @test canonicalize(CFoo(y=10)) == CFoo(x=ZeroTangent(), y=10)

        @test_throws ArgumentError canonicalize(CFoo(q=99.0, x=2.5))

        @testset "unspecified primal type" begin
            c1 = Tangent{Any}(;a=1, b=2)
            c2 = Tangent{Any}(1, 2)
            c3 = Tangent{Any}(Dict(4 => 3))

            @test c1 == canonicalize(c1)
            @test c2 == canonicalize(c2)
            @test c3 == canonicalize(c3)
        end
    end

    @testset "+ with other composites" begin
        @testset "Structs" begin
            CFoo = Tangent{Foo}
            @test  CFoo(x=1.5) + CFoo(x=2.5) == CFoo(x=4.0)
            @test CFoo(y=1.5) + CFoo(x=2.5) == CFoo(y=1.5, x=2.5)
            @test CFoo(y=1.5, x=1.5) + CFoo(x=2.5) == CFoo(y=1.5, x=4.0)
        end

        @testset "Tuples" begin
            @test ==(
                typeof(Tangent{Tuple{}}() + Tangent{Tuple{}}()),
                Tangent{Tuple{}, Tuple{}}
            )
            @test (
                Tangent{Tuple{Float64, Float64}}(1.0, 2.0) +
                Tangent{Tuple{Float64, Float64}}(1.0, 1.0)
            ) == Tangent{Tuple{Float64, Float64}}(2.0, 3.0)
        end

        @testset "NamedTuples" begin
            nt1 = (;a=1.5, b=0.0)
            nt2 = (;a=0.0, b=2.5)
            nt_sum = (a=1.5, b=2.5)
            @test (
                Tangent{typeof(nt1)}(; nt1...) +
                Tangent{typeof(nt2)}(; nt2...)
            ) == Tangent{typeof(nt_sum)}(; nt_sum...)
        end

        @testset "Dicts" begin
            d1 = Tangent{Dict}(Dict(4 => 3.0, 3 => 2.0))
            d2 = Tangent{Dict}(Dict(4 => 3.0, 2 => 2.0))
            d_sum = Tangent{Dict}(Dict(4 => 3.0 + 3.0, 3 => 2.0, 2 => 2.0))
            @test d1 + d2 == d_sum
        end

        @testset "Fields of type NotImplemented" begin
            CFoo = Tangent{Foo}
            a = CFoo(x=1.5)
            b = CFoo(x=@not_implemented(""))
            for (x, y) in ((a, b), (b, a), (b, b))
                z = x + y
                @test z isa CFoo
                @test z.x isa ChainRulesCore.NotImplemented
            end

            a = Tangent{Tuple}(1.5)
            b = Tangent{Tuple}(@not_implemented(""))
            for (x, y) in ((a, b), (b, a), (b, b))
                z = x + y
                @test z isa Tangent{Tuple}
                @test first(z) isa ChainRulesCore.NotImplemented
            end

            a = Tangent{NamedTuple{(:x,)}}(x=1.5)
            b = Tangent{NamedTuple{(:x,)}}(x=@not_implemented(""))
            for (x, y) in ((a, b), (b, a), (b, b))
                z = x + y
                @test z isa Tangent{NamedTuple{(:x,)}}
                @test z.x isa ChainRulesCore.NotImplemented
            end

            a = Tangent{Dict}(Dict(:x => 1.5))
            b = Tangent{Dict}(Dict(:x => @not_implemented("")))
            for (x, y) in ((a, b), (b, a), (b, b))
                z = x + y
                @test z isa Tangent{Dict}
                @test z[:x] isa ChainRulesCore.NotImplemented
            end
        end
    end

    @testset "+ with Primals" begin
        @testset "Structs" begin
            @test Foo(3.5, 1.5) + Tangent{Foo}(x=2.5) == Foo(6.0, 1.5)
            @test Tangent{Foo}(x=2.5) + Foo(3.5, 1.5) == Foo(6.0, 1.5)
            @test (@ballocated Bar(0.5) + Tangent{Bar}(; x=0.5)) == 0
        end

        @testset "Tuples" begin
            @test Tangent{Tuple{}}() + () == ()
            @test ((1.0, 2.0) + Tangent{Tuple{Float64, Float64}}(1.0, 1.0)) == (2.0, 3.0)
            @test (Tangent{Tuple{Float64, Float64}}(1.0, 1.0)) + (1.0, 2.0) == (2.0, 3.0)
        end

        @testset "NamedTuple" begin
            ntx = (; a=1.5)
            @test Tangent{typeof(ntx)}(; ntx...) + ntx == (; a=3.0)

            nty = (; a=1.5, b=0.5)
            @test Tangent{typeof(nty)}(; nty...) + nty == (; a=3.0, b=1.0)
        end

        @testset "Dicts" begin
            d_primal = Dict(4 => 3.0, 3 => 2.0)
            d_tangent = Tangent{typeof(d_primal)}(Dict(4 =>5.0))
            @test d_primal + d_tangent == Dict(4 => 3.0 + 5.0, 3 => 2.0)
        end
    end

    @testset "+ with Primals, with inner constructor" begin
        value = StructWithInvariant(10.0)
        diff = Tangent{StructWithInvariant}(x=2.0, x2=6.0)

        @testset "with and without debug mode" begin
            @assert ChainRulesCore.debug_mode() == false
            @test_throws MethodError (value + diff)
            @test_throws MethodError (diff + value)

            ChainRulesCore.debug_mode() = true  # enable debug mode
            @test_throws ChainRulesCore.PrimalAdditionFailedException (value + diff)
            @test_throws ChainRulesCore.PrimalAdditionFailedException (diff + value)
            ChainRulesCore.debug_mode() = false  # disable it again
        end

        # Now we define constuction for ChainRulesCore.jl's purposes:
        # It is going to determine the root quanity of the invarient
        function ChainRulesCore.construct(::Type{StructWithInvariant}, nt::NamedTuple)
            x = (nt.x + nt.x2/2)/2
            return StructWithInvariant(x)
        end
        @test value + diff == StructWithInvariant(12.5)
        @test diff + value == StructWithInvariant(12.5)
    end

    @testset "differential arithmetic" begin
        c = Tangent{Foo}(y=1.5, x=2.5)

        @test NoTangent() * c == NoTangent()
        @test c * NoTangent() == NoTangent()
        @test dot(NoTangent(), c) == NoTangent()
        @test dot(c, NoTangent()) == NoTangent()

        @test ZeroTangent() * c == ZeroTangent()
        @test c * ZeroTangent() == ZeroTangent()
        @test dot(ZeroTangent(), c) == ZeroTangent()
        @test dot(c, ZeroTangent()) == ZeroTangent()

        @test true * c === c
        @test c * true === c

        t = @thunk 2
        @test t * c == 2 * c
        @test c * t == c * 2
    end

    @testset "scaling" begin
        @test (
            2 *  Tangent{Foo}(y=1.5, x=2.5)
            == Tangent{Foo}(y=3.0, x=5.0)
            == Tangent{Foo}(y=1.5, x=2.5) * 2
        )
        @test (
            2 * Tangent{Tuple{Float64, Float64}}(2.0, 4.0)
            == Tangent{Tuple{Float64, Float64}}(4.0, 8.0)
            == Tangent{Tuple{Float64, Float64}}(2.0, 4.0) * 2
        )
        d = Tangent{Dict}(Dict(4 => 3.0))
        two_d = Tangent{Dict}(Dict(4 => 2 * 3.0))
        @test 2 * d == two_d == d * 2
    end

    @testset "show" begin
        @test repr(Tangent{Foo}(x=1,)) == "Tangent{Foo}(x = 1,)"
        # check for exact regex match not occurence( `^...$`)
        # and allowing optional whitespace (`\s?`)
        @test occursin(
            r"^Tangent{Tuple{Int64,\s?Int64}}\(1,\s?2\)$",
            repr(Tangent{Tuple{Int64,Int64}}(1, 2)),
        )
    end

    @testset "internals" begin
        @testset "Can't do backing on primative type" begin
            @test_throws Exception ChainRulesCore.backing(1.4)
        end

        @testset "Internals don't allocate a ton" begin
            bk = (; x=1.0, y=2.0)
            VERSION >= v"1.5" && @test (@ballocated(ChainRulesCore.construct($Foo, $bk))) <= 32

            # weaker version of the above (which should pass on all versions)
            @test (@ballocated(ChainRulesCore.construct($Foo, $bk))) <= 48
            @test (@ballocated ChainRulesCore.elementwise_add($bk, $bk)) <= 48
        end
    end

    @testset "non-same-typed differential arithmetic" begin
        nt = (; a=1, b=2.0)
        c = Tangent{typeof(nt)}(; a=NoTangent(), b=0.1)
        @test nt + c == (; a=1, b=2.1);
    end
end


text:
using CUDA
include("./types")

mutable struct exec_config
    pad::Int
    single_step_offset::Int
    has_multi_step::Bool
    m_step::Int

end

function config(st_inst, data)
    radius = st_inst.radius
    if st_inst.combined_time_step != false
        ctp = st_inst.combined_time_step
        radius = radius*(2^(ctp-1))

end

function apply_single_step()

    conf = config()
    i = 0

    while  i < t_steps
        ## Assign the single time step kernel
        foo = st_inst.kernel
        offset = conf.single_step_offset
        time_inc = 1
        ## If multistep feasible, assign the m_step kernel
        if st_inst.combined_time_step != false
            if i + st_inst.combined_time_step <= t_steps
                foo = st_inst.m_kernel
                offset = 0
                time_inc = conf.m_step
            end
        end
        ## Construct arguments for kenerl
        args = (dev_data, dev_out, 0, 0, offset, offset)
        if st_inst.uses_vsq
           args = (args..., dev_vsq)
        end
        ## Kernel call
        @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
                   shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
                   foo(args...))
        dev_data,dev_out = dev_out,dev_data
        println("t = $(i), (bx,by) = $((bx,by))")
        i += 1
    end
end

function apply_multi_step()

end


text:
"""
    SystemWithOutput{ST<:AbstractSystem, MT<:AbstractMap} <: AbstractSystem

Parametric composite type for systems with outputs. It is parameterized in the
system's type (`ST`) and in the map's type (`MT`).

### Fields

- `s`         -- system of type `ST`
- `outputmap` -- output map of type `MT`
"""
struct SystemWithOutput{ST<:AbstractSystem, MT<:AbstractMap} <: AbstractSystem
    s::ST
    outputmap::MT
end
statedim(swo::SystemWithOutput) = statedim(swo.s)
inputdim(swo::SystemWithOutput) = inputdim(swo.s)
inputset(swo::SystemWithOutput) = inputset(swo.s)
outputdim(swo::SystemWithOutput) = outputdim(swo.outputmap)
outputmap(swo::SystemWithOutput) = swo.outputmap
islinear(swo::SystemWithOutput) = islinear(swo.s) && islinear(swo.outputmap)
isaffine(swo::SystemWithOutput) = isaffine(swo.s) && isaffine(swo.outputmap)

"""
    LinearTimeInvariantSystem(A, B, C, D)

A linear time-invariant system with of the form

```math
x' = Ax + Bu,\\\\
y = Cx + Du.
```

### Input

- `A` -- matrix
- `B` -- matrix
- `C` -- matrix
- `D` -- matrix

### Output

A system with output such that the system is a linear control
continuous system and the output map is a linear control map.
"""
function LinearTimeInvariantSystem(A, B, C, D)
    system = LinearControlContinuousSystem(A, B)
    outputmap = LinearControlMap(C, D)
    return SystemWithOutput(system, outputmap)
end

"""
    LinearTimeInvariantSystem(A, B, C, D, X, U)

A linear time-invariant system with state and input constraints
of the form

```math
x' = Ax + Bu,\\\\
y = Cx + Du,
```
where ``x(t) ∈ X`` and ``u(t) ∈ U`` for all ``t``.

### Input

- `A` -- matrix
- `B` -- matrix
- `C` -- matrix
- `D` -- matrix
- `X` -- state constraints
- `U` -- input constraints

### Output

A system with output such that the system is a constrained linear control
continuous system and the output map is a constrained linear control map.
"""
function LinearTimeInvariantSystem(A, B, C, D, X, U)
    system = ConstrainedLinearControlContinuousSystem(A, B, X, U)
    outputmap = ConstrainedLinearControlMap(C, D, X, U)
    return SystemWithOutput(system, outputmap)
end

"""
    LTISystem

`LTISystem` is an alias for `LinearTimeInvariantSystem`.
"""
const LTISystem = LinearTimeInvariantSystem


text:
using ANema
using ATools
using DataFrames
using Test

dir = "../data/"
conf = "n3-w-20mm-phot"
path = joinpath(dir, conf)

@testset "ANema.jl" begin
    ndf = readdf(path)
    @test nrow(ndf) == 180079
    @test length(names(ndf)) ==  46

    ndfq =selectinterval(ndf, "q1", "q2", 1700.0, 2100.0);
    @test nrow(ndfq) == 42875

    ndfz =selectinterval(ndfq, "zstd1", "zstd2", 1.0, 16.0)
    @test nrow(ndfz) == 42806

    fz, pz = fit_profile(ndfz.zstd1, ndfz.r1, "σz", "r", "pol2")
    fpars = fz.fitpar
    fstds = fz.fitstd
    @test isapprox(fpars[1], 372.96, atol=0.1)

    r1z = fz.g.(ndfz.zstd1)
    fg,p = fitg1(r1z - ndfz.r1, "r1 - r1z", 200, -5.0, 5.0, xgmin=-1.5, xgmax=1.5)
    @test isapprox(fg.std[1], 0.4, atol=0.1)

end

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
using DynamicalSystems
using DataFrames
using StatsPlots
using BenchmarkTools

@inline function rate_to_proportion(r,t)
    1-exp(-r*t)
end;

function sir_map!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    infection = rate_to_proportion(β*c*I/N,δt)*S
    recovery = rate_to_proportion(γ,δt)*I
    @inbounds begin
        du[1] = S-infection
        du[2] = I+infection-recovery
        du[3] = R+recovery
    end
    nothing
end;

δt = 0.1
nsteps = 400
tmax = nsteps*δt
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0];

p = [0.05,10.0,0.25,δt]; # β,c,γ,δt

ds = DiscreteDynamicalSystem(sir_map!, u0, p, t0 = 0)

sol = trajectory(ds,nsteps)

df = DataFrame(Matrix(sol))
df[!,:t] = t;

@df df plot(:t,
    [:x1, :x2, :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark trajectory(ds,nsteps)


text:
const _Variable = VariableRef
const _Decision = Union{DecisionRef, KnownRef}
const _Scalar = Union{_Variable, _Decision, _Constant}
const _DecisionAffOrQuadExpr{C} = Union{DecisionAffExpr{C}, DecisionQuadExpr{C}}

function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               ::Type{<:_Constant}, S::Type{<:_Decision})
    return DAE
end
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               S::Type{<:_Decision}, ::Type{<:_Constant})
    return DAE
end
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               ::Type{<:_Constant}, S::Type{<:_DecisionAffOrQuadExpr})
    return S
end
function MA.promote_operation(::Union{typeof(+), typeof(-), typeof(*)},
                               S::Type{<:_DecisionAffOrQuadExpr}, ::Type{<:_Constant})
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, ::Type{<:Union{_Variable, _Decision}},
    S::Type{<:_DecisionAffOrQuadExpr})
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, S::Type{<:_DecisionAffOrQuadExpr},
    ::Type{<:Union{_Variable, _Decision}})
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, ::Type{GenericAffExpr{C,V}},
    S::Type{_DecisionAffOrQuadExpr{C}}) where {C, V <: Union{_Variable, _Decision}}
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, S::Type{_DecisionAffOrQuadExpr{C}},
    ::Type{GenericAffExpr{C,V}}) where {C, V <: Union{_Variable, _Decision}}
    return S
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, ::Type{GenericQuadExpr{C,V}},
    S::Type{_DecisionAffOrQuadExpr{C}}) where {C, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(
    ::Union{typeof(+), typeof(-)}, S::Type{_DecisionAffOrQuadExpr{C}},
    ::Type{GenericQuadExpr{C,V}}) where {C, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{C}
end
function MA.promote_operation(::Union{typeof(+), typeof(-)}, ::Type{A},
                               ::Type{A}) where {A <: _DecisionAffOrQuadExpr}
    return A
end
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               ::Type{<:DecisionAffExpr},
                               S::Type{<:DecisionQuadExpr})
    return S
end
function MA.promote_operation(::Union{typeof(+), typeof(-)},
                               S::Type{<:DecisionQuadExpr},
                               ::Type{<:DecisionAffExpr})
    return S
end

function MA.promote_operation(::typeof(*), ::Type{<:_Variable}, ::Type{<:_Decision})
    return DQE
end
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{<:_Variable})
    return DQE
end
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{<:_Decision})
    return DQE
end
function MA.promote_operation(::typeof(*), ::Type{<:_Variable}, ::Type{GenericAffExpr{T, D}}) where {T, D <: _Decision}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{T, D}}, ::Type{<:_Variable}) where {T, D <: _Decision}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{<:_Decision}, ::Type{GenericAffExpr{T, V}}) where {T, V <: _Variable}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{T, V}}, ::Type{<:_Decision}) where {T, V <: _Variable}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{D}, ::Type{GenericAffExpr{T, D}}) where {T, D <: _Decision}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{GenericAffExpr{T, D}}, ::Type{D}) where {T, D <: _Decision}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{V}, ::Type{DecisionAffExpr{T}}) where {T, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{T}
end
function MA.promote_operation(::typeof(*), ::Type{DecisionAffExpr{T}}, ::Type{V}) where {T, V <: Union{_Variable, _Decision}}
    return DecisionQuadExpr{T}
end

function MA.scaling(aff::DecisionAffExpr{C}) where C
    if !isempty(aff.variables.terms) || !isempty(aff.decisions.terms) || !isempty(aff.knowns.terms)
        throw(InexactError("Cannot convert `$aff` to `$C`."))
    end
    return MA.scaling(aff.variables.constant)
end
function MA.scaling(quad::DecisionQuadExpr{C}) where C
    if !isempty(quad.variables.terms) || !isempty(quad.decisions.terms) || !isempty(quad.decisions.quad.terms) || !isempty(quad.knowns.terms) || !isempty(quad.knowns.aff.terms) || !isempty(quad.cross_terms) || !isempty(quad.known_variable_terms.terms) || !isempty(quad.known_decisions_terms.terms)
        throw(InexactError("Cannot convert `$quad` to `$C`."))
    end
    return MA.scaling(quad.variables.aff)
end

MA.mutability(::Type{<:_DecisionAffOrQuadExpr}) = MA.IsMutable()
function MA.mutable_copy(expr::_DecisionAffOrQuadExpr)
    return map_coefficients(MA.copy_if_mutable, expr)
end

function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, aff::DecisionAffExpr)
    MA.mutable_operate!(op, aff.variables)
    MA.mutable_operate!(zerl, aff.decisions)
    MA.mutable_operate!(zero, aff.knowns)
    return aff
end
function MA.mutable_operate!(op::Union{typeof(zero), typeof(one)}, quad::DecisionQuadExpr)
    MA.mutable_operate!(op, quad.variables)
    MA.mutable_operate!(zero, quad.decisions)
    MA.mutable_operate!(zero, quad.knowns)
    empty!(quad.cross_terms)
    MA.mutable_operate!(zero, quad.known_variable_terms)
    MA.mutable_operate!(zero, quad.known_decisions_terms)
    return quad
end

function MA.mutable_operate!(::typeof(*), expr::_DecisionAffOrQuadExpr, α::_Constant)
    if iszero(α)
        return MA.mutable_operate!(zero, expr)
    else
        return map_coefficients_inplace!(x -> MA.mul!(x, α), expr)
    end
end

function MA.mutable_operate!(::typeof(+), expr::_DecisionAffOrQuadExpr, x)
    return JuMP.add_to_expression!(expr, x)
end
function MA.mutable_operate!(::typeof(-), expr::_DecisionAffOrQuadExpr, x)
    return JuMP.add_to_expression!(expr, -1, x)
end

function MA.mutable_operate!(::typeof(MA.add_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar)
    return JuMP.add_to_expression!(expr, x)
end
function MA.mutable_operate!(::typeof(MA.add_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar, y::_Scalar)
    return JuMP.add_to_expression!(expr, x, y)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar)
    return JuMP.add_to_expression!(expr, -1.0, x)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::_Scalar, y::_Scalar)
    return JuMP.add_to_expression!(expr, -x, y)
end
function MA.mutable_operate!(::typeof(MA.sub_mul), expr::_DecisionAffOrQuadExpr, x::Union{_Variable, _Decision}, y::_Constant)
    return JuMP.add_to_expression!(expr, x, -y)
end
function MA.mutable_operate!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, x, y)
    return MA.mutable_operate!(op, expr, x * y)
end
@generated function _add_sub_mul_reorder!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, args::Vararg{Any, N}) where N
    n = length(args)
    @assert n ≥ 3
    varidx = findall(t -> (t <: Union{_Variable, _Decision}), collect(args))
    allscalar = all(t -> (t <: _Constant), args[setdiff(1:n, varidx)])
    idx = (allscalar && length(varidx) == 1) ? varidx[1] : n
    coef = Expr(:call, :*, [:(args[$i]) for i in setdiff(1:n, idx)]...)
    return :(MA.mutable_operate!(op, expr, $coef, args[$idx]))
end
function MA.mutable_operate!(op::MA.AddSubMul, expr::_DecisionAffOrQuadExpr, x, y, z, other_args::Vararg{Any, N}) where N
    return _add_sub_mul_reorder!(op, expr, x, y, z, other_args...)
end
function JuMP.add_to_expression!(expr::_DecisionAffOrQuadExpr, α::_Constant, β::_Constant)
    return JuMP.add_to_expression!(expr, *(α, β))
end


text:
"""
Supertype of all dictionary solver operators. A solver operator typically implements an
efficient algorithm to apply the inverse of an operator. Examples include
a solver based on QR or SVD factorizations.

A solver operator stores the operator it has inverted.
"""
abstract type AbstractSolverOperator <: AbstractOperator
end

operator(op::AbstractSolverOperator) = op.op

src_space(op::AbstractSolverOperator) = dest_space(operator(op))

dest_space(op::AbstractSolverOperator) = src_space(operator(op))
dest(op::AbstractSolverOperator) = src(operator(op))

inv(op::AbstractSolverOperator) = operator(op)

"""
Supertype of all dictionary solver operators. A solver operator typically implements an
efficient algorithm to apply the inverse of an operator. Examples include
a solver based on QR or SVD factorizations.

A solver operator stores the operator it has inverted.
"""
abstract type DictionarySolverOperator{T} <: DictionaryOperator{T}
end

operator(op::DictionarySolverOperator) = op.op

src(op::DictionarySolverOperator) = dest(operator(op))

dest(op::DictionarySolverOperator) = src(operator(op))

inv(op::DictionarySolverOperator) = operator(op)

abstract type VectorizingSolverOperator{T} <: DictionarySolverOperator{T} end

apply!(op::VectorizingSolverOperator, coef_dest::Vector, coef_src::Vector) =
    linearized_apply!(op, coef_dest, coef_src)

function apply!(op::VectorizingSolverOperator, coef_dest, coef_src)
    copyto!(op.src_linear, coef_src)
    linearized_apply!(op, op.dest_linear, op.src_linear)
    copyto!(coef_dest, op.dest_linear)
end

function apply!(op::VectorizingSolverOperator, coef_dest, coef_src::Vector)
    linearized_apply!(op, op.dest_linear, coef_src)
    copyto!(coef_dest, op.dest_linear)
end

function apply!(op::VectorizingSolverOperator, coef_dest::Vector, coef_src)
    copyto!(op.src_linear, coef_src)
    linearized_apply!(op, coef_dest, op.src_linear)
end

"A GenericSolverOperator wraps around a generic solver type."
struct GenericSolverOperator{Q,T} <: VectorizingSolverOperator{T}
    op      ::  DictionaryOperator
    solver  ::  Q
    # In case the operator does not map between vectors, we allocate space
    # for two vectors so that we can convert between representations.
    src_linear  ::  Vector{T}
    dest_linear ::  Vector{T}

    GenericSolverOperator{Q,T}(op, solver) where {Q,T} =
        new{Q,T}(op, solver, zeros(T, length(dest(op))), zeros(T, length(src(op))))
end

linearized_apply!(op::GenericSolverOperator, dest_vc::Vector, src_vc::Vector) =
    linearized_apply!(op, dest_vc, src_vc, op.solver)

# The solver should be the inverse of the given operator
GenericSolverOperator(op::DictionaryOperator{T}, solver) where T =
    GenericSolverOperator{typeof(solver),T}(op, solver)

similar_operator(op::GenericSolverOperator, src, dest) =
    GenericSolverOperator(similar_operator(operator(op), dest, src), op.solver)

# This is the generic case
function linearized_apply!(op::GenericSolverOperator, coef_dest::Vector, coef_src::Vector, solver)
    coef_dest[:] = solver \ coef_src
    coef_dest
end

# # More efficient version for vectors and factorization solvers
# # TODO: these don't work for rectangular matrices
#     function _apply!(op::GenericSolverOperator, coef_dest::Vector, coef_src::Vector, solver::Factorization)
#         copyto!(coef_dest, coef_src)
#         ldiv!(solver, coef_dest)
#     end
#     function _apply!(op::GenericSolverOperator, coef_dest::Vector, coef_src::Vector, solver::LinearAlgebra.SVD)
#         copyto!(coef_dest, coef_src)
#         coef_dest[:] = ldiv!(solver, coef_dest)
#     end

adjoint(op::GenericSolverOperator) = @warn("not implemented")

if VERSION ≥ v"1.7-"
    qr_factorization(op::DictionaryOperator) = qr(matrix(op), ColumnNorm())
else
    qr_factorization(op::DictionaryOperator) = qr(matrix(op), Val(true))
end

function svd_factorization(op::DictionaryOperator)
    local F
    A = Matrix(op)
    try
        F = svd(A, full=false)
    catch
        # Sometimes svd does not work because A is a structured matrix.
        # In that case, we convert to a dense matrix.
        F = svd(collect(A), full=false)
    end
    return F
end

function regularized_svd_factorization(op::DictionaryOperator;
            verbose = false, threshold = 1000*eps(real(eltype(op))), options...)
    F = svd_factorization(op)
    U = F.U
    S = F.S
    Vt = F.Vt
    I = findfirst(F.S .< threshold)
    if I == nothing
        verbose && println("SVD: no regularization needed (threshold: $threshold)")
        return F
    else
        verbose && println("SVD truncated at rank $I of $(minimum(size(op))) (threshold: $threshold)")
        U_reg = U[:,1:I]
        Vt_reg = Vt[1:I,:]
        S_reg = S[1:I]
        return SVD(U_reg, S_reg, Vt_reg)
    end
end

QR_solver(op::DictionaryOperator; options...) =
    GenericSolverOperator(op, qr_factorization(op))
SVD_solver(op::DictionaryOperator; options...) =
    GenericSolverOperator(op, svd_factorization(op))

regularized_SVD_solver(op::DictionaryOperator; options...) =
    GenericSolverOperator(op, regularized_svd_factorization(op; options...))

for (damp, LS) in ((:λ, :LSMR), (:damp, :LSQR))
    OPTIONS = Symbol(LS,"Options")
    @eval begin
        export $OPTIONS
        mutable struct $OPTIONS{T}
            atol::T
            btol::T
            conlim::T
            $damp::T
            maxiter::Int
            verbose::Bool

            function $OPTIONS{T}(;atol=T(1e-6), btol=T(1e-6), conlim=T(Inf), λ=T(0), maxiter=-1, verbose=false) where T
                new{T}(atol, btol, conlim, λ, maxiter, verbose)
            end
            function $OPTIONS(;atol=1e-6, btol=1e-6, conlim=Inf, λ=0, maxiter=-1, verbose=false)
                T = promote_type(map(typeof, (atol, btol, conlim, λ))...)
                new{T}(T(atol), T(btol), T(conlim), T(λ), maxiter, verbose)
            end
            $OPTIONS{T}(itoptions::$OPTIONS{T}) where T =
                itoptions
            $OPTIONS{T}(itoptions::$OPTIONS) where T =
                $OPTIONS{T}(;convert(NamedTuple, itoptions)...)
        end

        convert(::Type{<:NamedTuple}, opts::$OPTIONS) =
            (atol=real(opts.atol), btol=real(opts.btol), conlim=real(opts.conlim), $damp=real(opts.$damp), verbose=real(opts.verbose), maxiter=real(opts.maxiter))
    end
end

export ItOptions
const ItOptions = LSMROptions
LSQROptions{T}(options::LSMROptions) where T =
    LSQROptions{T}(;convert(NamedTuple,options)...)

struct LSQR_solver{T} <: VectorizingSolverOperator{T}
    op      ::  DictionaryOperator{T}
    options ::  LSQROptions{T}

    src_linear  ::  Vector{T}
    dest_linear ::  Vector{T}
    LSQR_solver(op::DictionaryOperator{T}; itoptions=LSQROptions{T}(), options...) where T =
        LSQR_solver(op, LSQROptions{T}(itoptions); options...)

    function LSQR_solver(op::DictionaryOperator{T}, opts::LSQROptions{T}; verbose=false, options...) where T
        if opts.maxiter == -1
             opts.maxiter=100*max(size(op)...)
        end
        opts.verbose = verbose
        new{T}(op, opts,
            Vector{T}(undef, length(dest(op))), Vector{T}(undef, length(src(op)))
        )
    end
end

function linearized_apply!(op::LSQR_solver, coef_dest::Vector, coef_src::Vector)
    keys = convert(NamedTuple, op.options)
    op.options.verbose && @info "LSQR with $keys"
    copy!(coef_dest, lsqr(op.op, coef_src; keys...))
end

struct LSMR_solver{T} <: VectorizingSolverOperator{T}
    op      ::  DictionaryOperator{T}
    options ::  LSMROptions{T}

    src_linear  ::  Vector{T}
    dest_linear ::  Vector{T}

    LSMR_solver(op::DictionaryOperator{T}; itoptions=LSMROptions{T}(), options...) where T =
        LSMR_solver(op, LSMROptions{T}(itoptions); options...)

    function LSMR_solver(op::DictionaryOperator{T}, opts::LSMROptions{T}; verbose=false, options...) where T
        if opts.maxiter == -1
             opts.maxiter=100*max(size(op)...)
        end
        opts.verbose = verbose
        new{T}(op, opts,
            Vector{T}(undef, length(dest(op))), Vector{T}(undef, length(src(op)))
        )
    end
end

function linearized_apply!(op::LSMR_solver, coef_dest::Vector, coef_src::Vector)
    keys = convert(NamedTuple, op.options)
    op.options.verbose && @info "LSMR with $keys"
    copy!(coef_dest, lsmr(op.op, coef_src; keys...))
end


text:
# Minimal server using the 'listen' syntax, starting with the 'inner' functions
using Sockets
using WebSockets
using DelimitedFiles, LuxurySparse
import WebSockets.handle

export run_server, naive_handler

include("webparse.jl")

function coroutine(ws)
    state = YaoWSState()
    while isopen(ws)
        data, = readguarded(ws)
        s = String(data)
        println("Received: ", s)
        if strip(s) == "q"
            writeguarded(ws, "Goodbye!")
            break
        end
        try
            msg = yaoweb_handler(s, state)
            writeguarded(ws, string(msg))
        catch err
            writeguarded(ws, string(err))
        end
    end
end

function gatekeeper(req, ws)
    orig = WebSockets.origin(req)
    println("\nOrigin:", orig, "    Target:", target(req), "    subprotocol:", subprotocol(req))
    if occursin(WEBCONFIG[:LOCALIP], orig)
        coroutine(ws)
    else
        @info "Non-browser clients don't send Origin. We liberally accept the update request in this case."
        coroutine(ws)
    end
end

function run_server(handler)
    LOCALIP = WEBCONFIG[:LOCALIP]
    PORT = WEBCONFIG[:PORT]
    @info("Browser > $LOCALIP:$PORT , F12> console > ws = new WebSocket(\"ws://$LOCALIP:$PORT\") ")

    handler_wrap = WebSockets.RequestHandlerFunction(req->handler(req))
    SERVER = Sockets.listen(Sockets.InetAddr(parse(IPAddr, LOCALIP), PORT))
    @async try
        WebSockets.HTTP.listen(LOCALIP, PORT, server = SERVER, readtimeout = 0 ) do http
            if WebSockets.is_upgrade(http.message)
                WebSockets.upgrade((req, ws)->gatekeeper(req, ws), http)
            else
                handle(handler_wrap, http)
            end
        end
    catch err
        # Add your own error handling code; HTTP.jl sends error code to the client.
        @info err
        @info stacktrace(catch_backtrace())
    end
    @info("To stop serving: close(SERVER)")
    return SERVER
end

function naive_handler(req)
    BAREHTML = "<meta http-equiv=\"Content-Type\" content=\"text/html; charset=UTF-8\">
    <title>Empty.html</title>"

    LOCALIP = WEBCONFIG[:LOCALIP]
    PORT = WEBCONFIG[:PORT]
    BODY =  "
Press F12

ws = new WebSocket(\"ws://$LOCALIP:$PORT\")

ws.onmessage = function(e){console.log(e.data)}

ws.send(\"hello!\")
                "

    replace(BAREHTML, "" => BODY) |> WebSockets.Response
end


text:
using Debug

@debug function asianOpt(N=10000,T=100; S0=100.0,K= 100.0,r=0.05,v=0.2,tma=0.25)
  dt = tma/T;
  S = zeros(Float64,T);
  A = zeros(Float64,N);
  for n = 1:N
    S[1] = S0
    dW = randn(T)*sqrt(dt);
    for t = 2:T
      z0 = (r - 0.5*v*v)*S[t-1]*dt;
      z1 = v*S[t-1]*dW[t];
      z2 = 0.5*v*v*S[t-1]*dW[t]*dW[t];
      S[t] = S[t-1] + z0 + z1 + z2;
    end
    A[n] = mean(S);
    @bp
  end
  P = zeros(Float64,N);
  [ P[n] = max(A[n] - K, 0) for n = 1:N ];
  @bp
  price = exp(-r*tma)*mean(P);
end

max_stars_repo_path: Module 3/Chapter10/asian-debug.jl

max_issues_repo_path: Module 3/Chapter10/asian-debug.jl

max_forks_repo_path: Module 3/Chapter10/asian-debug.jl


text:
# Object-based wrapper for level-3 BLIS routines.
#

macro blis_object_api_level3_gemm(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α,     xObj,
                           a,     xObj,
                           b,     xObj,
                           β,     xObj,
                           c,     xObj)
    end
end

@blis_object_api_level3_gemm gemm
@blis_object_api_level3_gemm gemmt
@blis_object_api_level3_gemm her2k
@blis_object_api_level3_gemm syr2k

macro blis_object_api_level3_hemm(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           sidea, BliSide,
                           α,     xObj,
                           a,     xObj,
                           b,     xObj,
                           β,     xObj,
                           c,     xObj)
    end
end

@blis_object_api_level3_hemm hemm
@blis_object_api_level3_hemm symm
@blis_object_api_level3_hemm trmm3

macro blis_object_api_level3_herk(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α,     xObj,
                           a,     xObj,
                           β,     xObj,
                           c,     xObj)
    end
end

@blis_object_api_level3_herk herk
@blis_object_api_level3_herk syrk

macro blis_object_api_level3_trmm(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           sidea, BliSide,
                           α,     xObj,
                           a,     xObj,
                           b,     xObj)
    end
end

@blis_object_api_level3_trmm trmm
@blis_object_api_level3_trmm trsm


text:
# Functions that take as first input a named tuple or dictionary with different values
# and as second argument

# Note: I copy the Dict -> Tuple transfrom from NamedTupleTools.jl as it was not dispatching to the correct method...
function initialize_dual(values::NamedTuple, names::NTuple)
    d = initialize_dual(NamedTupleTools.dictionary(values), names)
    parts = (d...,)
    names = first.(parts)
    vals  = last.(parts)
    return NamedTuple{names}(vals)
end

# Note: This function is type unstable
function initialize_dual(values::Dict, names)
    # Check the inputs
    @assert length(values) >= length(names)
    all_names = Symbol.(keys(values))
    for name in names
        @assert name in all_names
    end
    # Convert to UnitDual with partials
    ref_values = [values[name] for name in names]
    new_values = [initialize_dual(values, names, name, ref_values) for name in all_names]
    # Return as named tuple
    Dict(zip(all_names, new_values))
end

function initialize_dual(values, names, name, ref_values)
    val = values[name]
    partials_values = [zero(val/refval) for refval in ref_values]
    if name in names
        partials_values[findfirst(x -> x == name, names)] = 1.0
    end
    partials = UnitPartials(NamedTupleTools.namedtuple(names, partials_values))
    UnitDual(val, partials)
end


text:
module ReverseDiff

using LinearAlgebra

export track, accum!, σ
export TrackedArray, TrackedMatrix, TrackedVector, TrackedReal

mutable struct TrackedArray{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::A
    grad::Union{Nothing,A}
    children::Dict
end

const TrackedVector{T} = TrackedArray{T,1}
const TrackedMatrix{T} = TrackedArray{T,2}

mutable struct TrackedReal{T<:Real}
    data::T
    grad::Union{Nothing,T}
    children::Dict
end
Base.show(io::IO, x::TrackedReal) = print(io, "$(x.data) (tracked)")

track(x::Real) = TrackedReal(x, nothing, Dict())
data(x::TrackedReal) = x.data

track(x::Array) = TrackedArray(x, nothing, Dict())
data(x::TrackedArray) = x.data
for f in [:size, :length, :eltype]
	eval(:(Base.$(f)(x::TrackedArray, args...) = $(f)(x.data, args...)))
end
Base.getindex(x::TrackedArray, args...) = track(getindex(x.data,args...))
Base.show(io::IO, x::TrackedArray) = print(io, "Tracked $(x.data)")
Base.print_array(io::IO, x::TrackedArray) = Base.print_array(io, x.data)

track(x::Union{TrackedArray,TrackedReal}) = x
function reset!(x::Union{TrackedReal,TrackedArray})
    x.grad = nothing
    x.children = Dict()
    nothing
end

function accum!(x::Union{TrackedReal,TrackedArray})
    if isnothing(x.grad)
        x.grad = sum(λ(accum!(Δ)) for (Δ,λ) in x.children)
    end
    x.grad
end

function gradient(f, args::TrackedArray...)
    y = f(args...)
    (y isa TrackedReal) || error("Output of `f` must be a scalar.")
    y.grad = 1.0
    accum!.(args)
    Tuple(a.grad for a in args)
end

##########  RULES  #############################################################

function Base.sum(x::TrackedArray)
    z = track(sum(x.data))
    x.children[z] = Δ -> Δ*ones(eltype(x), size(x)...)
    z
end

function Base.:*(X::TrackedArray, Y::TrackedArray)
    Z = track(X.data * Y.data)
    X.children[Z] = Δ -> Δ * Y.data'
    Y.children[Z] = Δ -> X.data' * Δ
    Z
end

function Base.:+(X::TrackedArray, Y::TrackedArray)
    Z = track(X.data + Y.data)
    X.children[Z] = Δ -> Δ
    Y.children[Z] = Δ -> Δ
    Z
end

function Base.:-(X::TrackedArray, Y::TrackedArray)
    Z = track(X.data - Y.data)
    X.children[Z] = Δ -> Δ
    Y.children[Z] = Δ -> -Δ
    Z
end

σ(x::Real) = 1/(1+exp(-x))
σ(x::AbstractArray) = σ.(x)
function σ(x::TrackedArray)
    z = track(σ(x.data))
    d = z.data
    x.children[z] = Δ -> Δ .* d .* (1 .- d)
    z
end

function Base.abs2(x::TrackedArray)
    y = track(abs2.(x.data))
    x.children[y] = Δ -> Δ .* 2x.data
    y
end

function Base.hcat(xs::TrackedArray...)
    y  = track(hcat(data.(xs)...))
    stops  = cumsum([size(x,2) for x in xs])
    starts = vcat([1], stops[1:end-1] .+ 1)
    for (start,stop,x) in zip(starts,stops,xs)
        x.children[y] = function (Δ)
            δ = if ndims(x) == 1
                Δ[:,start]
            else
                ds = map(_ -> :, size(x)) |> Base.tail |> Base.tail
                Δ[:, start:stop, ds...]
            end
            δ
        end
    end
    y
end

end


text:
module CLIntrinsics

using MacroTools

import ..Transpiler: AbstractCLIO, EmptyCLIO
import ..Transpiler: ints, floats, numbers, Numbers, Floats, int, Ints, uchar
import ..Transpiler: fixed_array_length, is_fixedsize_array
import ..Transpiler: ret, vecs, Vecs, vector_lengths, functions

using StaticArrays, Sugar
import Sugar: typename, vecname
using SpecialFunctions: erf, erfc

# TODO, these are rather global pointers and this should be represented in the type
immutable GlobalPointer{T} end
immutable LocalPointer{T} end
immutable LocalArray{T} end
Base.eltype(::Type{GlobalPointer{T}}) where T = T
Base.eltype(::Type{LocalPointer{T}}) where T = T
Base.eltype(::Type{LocalArray{T}}) where T = T

const DevicePointer = Union{GlobalPointer, LocalPointer}
const Types = Union{vecs..., numbers..., GlobalPointer, LocalPointer}

#########
# GLOBALS
const CLK_LOCAL_MEM_FENCE = Cuint(0)
const CLK_GLOBAL_MEM_FENCE = Cuint(0)

const intrinsic_signatures = Dict{Function, Any}()

macro cl_intrinsic(expr)
    matched = @capture(
        expr,
        (func_(args__) where {T__} = body_) | (func_(args__) = body_)
    )
    @assert matched "internal error: intrinsic wasn't matched: $expr"
    ret_expr = Expr(:block)

    # it's possible to define methods in base as intrinsic.
    # if they're not in base, we need to define a function stub
    if !isdefined(Base, func)
        push!(ret_expr.args, esc(expr))
    end
    types = map(args) do arg
        @assert isa(arg, Expr) && arg.head == :(::) "wrong type declaration"
        arg.args[end]
    end
    tuple_typ = :(Tuple{$(types...)})
    if T != nothing
        tuple_typ = :($tuple_typ where {$(T...)})
    end
    push!(ret_expr.args, :(intrinsic_signatures[$func] = $tuple_typ))
    ret_expr
end

@cl_intrinsic get_num_groups(dim::Integer) = ret(Cuint)
@cl_intrinsic get_global_id(dim::Integer) = ret(Cuint)
@cl_intrinsic get_local_id(::Integer) = ret(Cuint)
@cl_intrinsic get_group_id(::Integer) = ret(Cuint)
@cl_intrinsic get_local_size(::Integer) = ret(Cuint)
@cl_intrinsic get_global_size(::Integer) = ret(Cuint)
@cl_intrinsic get_work_dim() = ret(Cuint)
@cl_intrinsic sizeof(x::Any) = ret(Cuint)

@cl_intrinsic intrinsic_select(a::T, b::T, c::Any) where T = ret(T)

for (a, b) in (
        Float32 => UInt32,
        Float64 => UInt64,
        Int32 => UInt32,
        Int64 => UInt64,
        Bool => Bool
    )
    @eval cl_select(a::$a, b::$a, c::Bool) = intrinsic_select(b, a, $b(c))
end
cl_select(a, b, c::Bool) = c ? a : b

# @cl_intrinsic clt(::T, ::T, ::Bool) where {T} = ret(T)

@cl_intrinsic barrier(::Cuint) = nothing
@cl_intrinsic mem_fence(::Cuint) = nothing

@cl_intrinsic erfc(::T) where T <: Floats = ret(T)
@cl_intrinsic erf(::T) where T <: Floats = ret(T)
@cl_intrinsic remainder(::T, ::T) where T <: Floats = ret(T)
@cl_intrinsic exp(::T) where T <: Floats = ret(T)

for N in vector_lengths
    fload = Symbol(string("vload", N))
    fstore = Symbol(string("vstore", N))
    @eval begin
        @cl_intrinsic $(fload)(i::Integer, a::GlobalPointer{T}) where {T <: Numbers} = ret(NTuple{$N, T})
        @cl_intrinsic $(fstore)(x::NTuple{$N, T}, i::Integer, a::GlobalPointer{T}) where {T <: Numbers} = nothing
    end
end

end # end CLIntrinsics

using .CLIntrinsics

const cli = CLIntrinsics

import .cli: GlobalPointer, DevicePointer

import Sugar: typename, isintrinsic

function is_native_type(m::CLMethod, T)
    T <: cli.Types || is_fixedsize_array(m, T) || T <: Tuple{T} where T <: cli.Numbers
end

function isintrinsic(m::CLMethod, func::ANY, sig_tuple::ANY)
    # constructors are intrinsic. TODO more thorow lookup to match actual inbuild constructor
    isa(func, DataType) && is_native_type(m, func) && return true
    func == tuple && return true # TODO match against all Base intrinsics?
    func == getfield && sig_tuple <: (Tuple{X, Symbol} where X) && return true
    func == getfield && sig_tuple <: (Tuple{X, Integer} where X <: Tuple) && return true
    func == muladd && sig_tuple <: (Tuple{T, T, T} where T <: AbstractFloat) && return true
    Sugar.isintrinsic(func) && return true
    #TODO better julia/cl intrinsic matching
    func == Base.select_value && return true
    # Symbol(func) == Symbol("GPUArrays.LocalMemory") && return true
    # shared intrinsic functions should all work on all native types.
    # TODO, find exceptions where this isn't true
    func in functions && all(x-> is_native_type(m, x), Sugar.to_tuple(sig_tuple)) && return true
    haskey(cli.intrinsic_signatures, func) || return false
    sig = cli.intrinsic_signatures[func]
    sig_tuple <: sig
end

Base.getindex{T}(a::cli.LocalPointer{T}, i::Integer) = cli.ret(T)
Base.getindex{T}(a::GlobalPointer{T}, i::Integer) = cli.ret(T)

Base.setindex!{T}(::cli.LocalPointer{T}, ::T, ::Integer) = nothing
Base.setindex!{T}(a::GlobalPointer{T}, value::T, i::Integer) = nothing
function Base.setindex!(a::GlobalPointer{T}, value::T2, i::Integer) where {T, T2}
    setindex!(a, T(value), i)
    nothing
end

# TODO overload SIMD.vload, so that code can run seamlessly on the CPU as well.
for VecType in (NTuple, SVector)
    for N in cli.vector_lengths
        fload = Symbol(string("vload", N))
        fstore = Symbol(string("vstore", N))
        VType = :($VecType{$N, T})
        @eval begin
            function vload{T <: cli.Numbers, IT <: Integer}(::Type{$VType}, a::GlobalPointer{$VType}, i::IT)
                $VType(cli.$(fload)(i - IT(1), GlobalPointer{T}(a)))
            end
            function vstore{T <: cli.Numbers, IT <: Integer}(x::$VType, a::GlobalPointer{$VType}, i::IT)
                cli.$(fstore)(Tuple(x), i - IT(1), GlobalPointer{T}(a))
            end
        end
    end
end

Base.getindex{T <: Vecs}(a::GlobalPointer{T}, i::Integer) = vload(T, a, i)
function Base.setindex!{T <: Vecs}(a::GlobalPointer{T}, value::T, i::Integer)
    vstore(value, a, i)
end

function supports_indexing(m::LazyMethod, ::Type{T}) where T
    is_fixedsize_array(T) || T <: DevicePointer
end

function supports_indices(m::LazyMethod, ::Type{<: GlobalPointer{T}}, index_types) where T
    is_fixedsize_array(m, T) && return false # fixed size arrays are implemented via vstore/load
    length(index_types) == 1 && index_types[1] <: Integer
end
function supports_indices(m::LazyMethod, ::Type{<: DevicePointer}, index_types)
    length(index_types) == 1 && index_types[1] <: Integer
end
function supports_indices(m::LazyMethod, ::Type{<: Tuple}, index_types)
    length(index_types) == 1 && index_types[1] <: Integer
end

function typename(io::AbstractCLIO, x::Type{GlobalPointer{T}}) where T
    tname = typename(io, T)
    # restrict should be fine for now, since we haven't implemented views yet!
    "__global $tname * "
end
function typename(io::AbstractCLIO, x::Type{cli.LocalPointer{T}}) where T
    tname = typename(io, T)
    "__local $tname * "
end

function Sugar.vecname(io::AbstractCLIO, t::Type{T}) where T
    N = fixed_array_length(T)
    ET = eltype(T)
    etname = typename(io, ET)
    return string(etname, N)
end


text:
println("===============================================")
println("Testing kelvinvoigt.jl")
println("===============================================")

function Fract_KelvinVoigt_G_reduce()
    dt = 0.001
    # fractional modulus is NaN not inf at 0.0 so start at dt
    t = collect(dt:dt:10.0)

    fm_model = relaxmod(Fract_KelvinVoigt, t, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = relaxmod(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(t))
end
@test Fract_KelvinVoigt_G_reduce()

function FractS_KelvinVoigt_G_reduce()
    dt = 0.001
    # fractional modulus is NaN not inf at 0.0 so start at dt
    t = collect(dt:dt:10.0)

    fms_model = relaxmod(FractS_KelvinVoigt, t, [1.0, 1.0, 1.0])
    kelvinvoigt = relaxmod(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(t))
end
@test FractS_KelvinVoigt_G_reduce()

function FractD_KelvinVoigt_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = relaxmod(FractD_KelvinVoigt, t, [1.0, 1.0, 0.0])
    kelvinvoigt = relaxmod(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(t))
end
@test FractD_KelvinVoigt_G_reduce()

function Fract_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = creepcomp(Fract_KelvinVoigt, t, [1.0, 0.99999, 1.0, 0.00001])
    kelvinvoigt = creepcomp(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test Fract_KelvinVoigt_J_reduce(tol)

function FractS_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = creepcomp(FractS_KelvinVoigt, t, [1.0, 0.99999, 1.0])
    kelvinvoigt = creepcomp(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test FractS_KelvinVoigt_J_reduce(tol)

function FractD_KelvinVoigt_J_reduce(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = creepcomp(FractD_KelvinVoigt, t, [1.0, 1.0, 0.00001])
    kelvinvoigt = creepcomp(KelvinVoigt, t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i], atol=tol), eachindex(t))
end
@test FractD_KelvinVoigt_J_reduce(tol)

function Fract_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = storagemod(Fract_KelvinVoigt, chirp, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = storagemod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test Fract_KelvinVoigt_Gp_reduce()

function FractS_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = storagemod(FractS_KelvinVoigt, chirp, [1.0, 1.0, 1.0])
    kelvinvoigt = storagemod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractS_KelvinVoigt_Gp_reduce()

function FractD_KelvinVoigt_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = storagemod(FractD_KelvinVoigt, chirp, [1.0, 1.0, 0.0])
    kelvinvoigt = storagemod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractD_KelvinVoigt_Gp_reduce()

function Fract_KelvinVoigt_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = lossmod(Fract_KelvinVoigt, chirp, [1.0, 1.0, 1.0, 0.0])
    kelvinvoigt = lossmod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test Fract_KelvinVoigt_Gpp_reduce()

function FractS_KelvinVoigt_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = lossmod(FractS_KelvinVoigt, chirp, [1.0, 1.0, 1.0])
    kelvinvoigt = lossmod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractS_KelvinVoigt_Gpp_reduce()

function FractD_KelvinVoigt_Gpp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = lossmod(FractD_KelvinVoigt, chirp, [1.0, 1.0, 0.0])
    kelvinvoigt = lossmod(KelvinVoigt, chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], kelvinvoigt[i]), eachindex(chirp))
end
@test FractD_KelvinVoigt_Gpp_reduce()


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

mutable struct TroubleshootingResult <: SwaggerModel
    startTime::Any # spec type: Union{ Nothing, DateTime } # spec name: startTime
    endTime::Any # spec type: Union{ Nothing, DateTime } # spec name: endTime
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    results::Any # spec type: Union{ Nothing, Vector{TroubleshootingDetails} } # spec name: results

    function TroubleshootingResult(;startTime=nothing, endTime=nothing, code=nothing, results=nothing)
        o = new()
        validate_property(TroubleshootingResult, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        validate_property(TroubleshootingResult, Symbol("endTime"), endTime)
        setfield!(o, Symbol("endTime"), endTime)
        validate_property(TroubleshootingResult, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(TroubleshootingResult, Symbol("results"), results)
        setfield!(o, Symbol("results"), results)
        o
    end
end # type TroubleshootingResult

const _property_map_TroubleshootingResult = Dict{Symbol,Symbol}(Symbol("startTime")=>Symbol("startTime"), Symbol("endTime")=>Symbol("endTime"), Symbol("code")=>Symbol("code"), Symbol("results")=>Symbol("results"))
const _property_types_TroubleshootingResult = Dict{Symbol,String}(Symbol("startTime")=>"DateTime", Symbol("endTime")=>"DateTime", Symbol("code")=>"String", Symbol("results")=>"Vector{TroubleshootingDetails}")
Base.propertynames(::Type{ TroubleshootingResult }) = collect(keys(_property_map_TroubleshootingResult))
Swagger.property_type(::Type{ TroubleshootingResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TroubleshootingResult[name]))}
Swagger.field_name(::Type{ TroubleshootingResult }, property_name::Symbol) =  _property_map_TroubleshootingResult[property_name]

function check_required(o::TroubleshootingResult)
    true
end

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


text:
import Base.+, Base.-, Base.*, LinearAlgebra.mul!

function mul!(res::AbstractVector, op::AbstractLinearOperator{T}, v::AbstractVector, α, β) where T
  (size(v, 1) == size(op, 2) && size(res, 1) == size(op, 1)) || throw(LinearOperatorException("shape mismatch"))
  increase_nprod(op)
  op.prod!(res, v, α, β)
end

function mul!(res::AbstractVector, op::AbstractLinearOperator{T}, v::AbstractVector) where T
  mul!(res, op, v, one(T), zero(T))
end

# Apply an operator to a vector.
function *(op::AbstractLinearOperator{T}, v::AbstractVector{S}) where {T, S}
  nrow, ncol = size(op)
  res = Vector{promote_type(T, S)}(undef, nrow)
  res .= 0 # in case v has some NaN
  mul!(res, op, v)
  return res
end

# Unary operations.
+(op::AbstractLinearOperator) = op

function -(op::AbstractLinearOperator{T}) where T
  prod! = @closure (res, v, α, β) -> mul!(res, op, v, -α, β)
  tprod! = @closure (res, u, α, β) -> mul!(res, transpose(op), u, -α, β)
  ctprod! = @closure (res, w, α, β) -> mul!(res, adjoint(op), w, -α, β)
  LinearOperator{T}(op.nrow, op.ncol, op.symmetric, op.hermitian, prod!, tprod!, ctprod!, op.Mv, op.Mtu, op.Maw)
end

function prod_op!(res::AbstractVector, op1::AbstractLinearOperator, op2::AbstractLinearOperator,
                  vtmp::AbstractVector, v::AbstractVector, α, β)
  mul!(vtmp, op2, v)
  mul!(res, op1, vtmp, α, β)
end

function tprod_op!(res::AbstractVector, op1::AbstractLinearOperator, op2::AbstractLinearOperator,
                   utmp::AbstractVector, u::AbstractVector, α, β)
  mul!(utmp, transpose(op1), u)
  mul!(res, transpose(op2), utmp, α, β)
end

function ctprod_op!(res::AbstractVector, op1::AbstractLinearOperator, op2::AbstractLinearOperator,
                    wtmp::AbstractVector, w::AbstractVector, α, β)
  mul!(wtmp, adjoint(op1), w)
  mul!(res, adjoint(op2), wtmp, α, β)
end

## Operator times operator.
function *(op1::AbstractLinearOperator, op2::AbstractLinearOperator)
  T = promote_type(eltype(op1), eltype(op2))
  (m1, n1) = size(op1)
  (m2, n2) = size(op2)
  if m2 != n1
    throw(LinearOperatorException("shape mismatch"))
  end
  #tmp vector for products
  if typeof(op2) <: AdjointLinearOperator || typeof(op2) <: TransposeLinearOperator
    vtmp = op2.parent.Mtu
  elseif typeof(op2) <: ConjugateLinearOperator
    vtmp = op2.parent.Mv
  else
    vtmp = op2.Mv
  end
  if typeof(op1) <: AdjointLinearOperator || typeof(op1) <: TransposeLinearOperator
    utmp = op1.parent.Mv
    wtmp = op1.parent.Mv
  elseif typeof(op1) <: ConjugateLinearOperator
    utmp = op1.parent.Mtu
    wtmp = op1.parent.Maw
  else
    utmp = op1.Mtu
    wtmp = op1.Maw
  end
  prod! = @closure (res, v, α, β) -> prod_op!(res, op1, op2, vtmp, v, α, β)
  tprod! = @closure (res, u, α, β) -> tprod_op!(res, op1, op2, utmp, u, α, β)
  ctprod! = @closure (res, w, α, β) -> ctprod_op!(res, op1, op2, wtmp, w, α, β)
  Mv = Vector{T}(undef, m1)
  Mtu = Vector{T}(undef, n2)
  Maw = Vector{T}(undef, n2)
  LinearOperator{T}(m1, n2, false, false, prod!, tprod!, ctprod!, Mv, Mtu, Maw)
end

## Matrix times operator.
*(M::AbstractMatrix, op::AbstractLinearOperator) = LinearOperator(M) * op
*(op::AbstractLinearOperator, M::AbstractMatrix) = op * LinearOperator(M)

## Scalar times operator. (# commutation α*v ???)
function *(op::AbstractLinearOperator, x::Number)
  S = promote_type(eltype(op), typeof(x))
  prod! = @closure (res, v, α, β) -> mul!(res, op, v, x * α, β)
  tprod! = @closure (res, u, α, β) -> mul!(res, transpose(op), u, x * α, β)
  ctprod! = @closure (res, w, α, β) -> mul!(res, adjoint(op), w, x' * α, β)
  LinearOperator{S}(op.nrow, op.ncol, op.symmetric, op.hermitian && isreal(x), prod!, tprod!, ctprod!, op.Mv, op.Mtu, op.Maw)
end

function *(x::Number, op::AbstractLinearOperator)
  return op * x
end

# Operator + operator.

function sum_prod!(res::AbstractVector, op1::AbstractLinearOperator, op2::AbstractLinearOperator{T}, v::AbstractVector, α, β) where T
  mul!(res, op1, v, α, β)
  mul!(res, op2, v, α, one(T))
end

function +(op1::AbstractLinearOperator, op2::AbstractLinearOperator)
  (m1, n1) = size(op1)
  (m2, n2) = size(op2)
  if (m1 != m2) || (n1 != n2)
    throw(LinearOperatorException("shape mismatch"))
  end
  S = promote_type(eltype(op1), eltype(op2))
  prod! = @closure (res, v, α, β) -> sum_prod!(res, op1, op2, v, α, β)
  tprod! = @closure (res, u, α, β) -> sum_prod!(res, transpose(op1), transpose(op2), u, α, β)
  ctprod! = @closure (res, w, α, β) -> sum_prod!(res, adjoint(op1), adjoint(op2), w, α, β)
  symm = (symmetric(op1) && symmetric(op2))
  herm = (hermitian(op1) && hermitian(op2))
  Mv = Vector{S}(undef, m1)
  Mtu = symm ? Mv : Vector{S}(undef, n1)
  Maw = herm ? Mv : Vector{S}(undef, n1)
  return LinearOperator{S}(
    m1,
    n1,
    symm,
    herm,
    prod!,
    tprod!,
    ctprod!,
    Mv,
    Mtu,
    Maw
  )
end

# Operator + matrix.
+(M::AbstractMatrix, op::AbstractLinearOperator) = LinearOperator(M) + op
+(op::AbstractLinearOperator, M::AbstractMatrix) = op + LinearOperator(M)

# # Operator .+ scalar.
+(op::AbstractLinearOperator, x::Number) = op + x * opOnes(op.nrow, op.ncol)
+(x::Number, op::AbstractLinearOperator) = x * opOnes(op.nrow, op.ncol) + op

# Operator - operator
-(op1::AbstractLinearOperator, op2::AbstractLinearOperator) = op1 + (-op2)

# Operator - matrix.
-(M::AbstractMatrix, op::AbstractLinearOperator) = LinearOperator(M) - op
-(op::AbstractLinearOperator, M::AbstractMatrix) = op - LinearOperator(M)

# # Operator - scalar.
-(op::AbstractLinearOperator, x::Number) = op + (-x)
-(x::Number, op::AbstractLinearOperator) = x + (-op)


text:
using OrdinaryDiffEq, DiffEqBase, DiffEqCallbacks, Test
using Random, DiffEqOperators
using ElasticArrays
Random.seed!(213)
CACHE_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK22(),SSPRK33(),SSPRK43(),SSPRK104(),
  CarpenterKennedy2N54(), SHLDDRK64(), ORK256(), DGLDDRK73_C(),
  CFRLDDRK64(), TSLDDRK74(),
  CKLLSRK43_2(),
  ParsaniKetchesonDeconinck3S32(),
  BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),Feagin14(),TanYam7(),
  Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9(),OwrenZen3(),OwrenZen4(),OwrenZen5(),
  AutoTsit5(Rosenbrock23()),TRBDF2(),KenCarp4(),ABDF2()]
broken_CACHE_TEST_ALGS = [QNDF(),ExtrapolationMidpointHairerWanner(),ImplicitEulerExtrapolation(),ImplicitDeuflhardExtrapolation()]
# AitkenNeville(threading=false) fails Elastic but not normal case

using InteractiveUtils

NON_IMPLICIT_ALGS = filter((x)->isconcretetype(x) && !OrdinaryDiffEq.isimplicit(x()),union(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAlgorithm),subtypes(OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm)))

f = function (du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.3/length(u))*u[i]
  end
end

condition = function (u,t,integrator)
  1-maximum(u)
end

affect! = function (integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()/5 + 0.25
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition,affect!)

u0 = [0.2]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)

println("Check for stochastic errors")
for i in 1:10
  @test_nowarn sol = solve(prob,Tsit5(),callback=callback)
end

println("Check some other integrators")
sol = solve(prob,Rosenbrock23(),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,Rosenbrock32(),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,KenCarp4(),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,TRBDF2(),callback=callback,dt=1/2)
@test length(sol[end]) > 1

Jv = JacVecOperator(f,u0,nothing,0.0)
ff2 = ODEFunction(f;jac_prototype=Jv)
prob2 = ODEProblem(ff2,u0,tspan)
sol = solve(prob2,TRBDF2(linsolve=LinSolveGMRES()),callback=callback)

for alg in CACHE_TEST_ALGS
  @show alg
  local sol = solve(prob,alg,callback=callback,dt=1/2)
  @test length(sol[end]) > 1
end

for alg in broken_CACHE_TEST_ALGS
  @show alg
  @test_broken length(solve(prob,alg,callback=callback,dt=1/2)[end]) > 1
end

sol = solve(prob,Rodas4(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,Rodas5(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1

# cache tests resizing multidimensional arrays
println("Check resizing multidimensional arrays")
u0_matrix = ElasticArray(ones(2, 2))
f_matrix = (du, u, p, t) -> du .= u
prob_matrix = ODEProblem(f_matrix, u0_matrix, (0.0, 2.0))
condition_matrix = (u, t, integrator) -> t - 1
affect_matrix! = function (integrator)
  resize!(integrator, (2,3))
  integrator.u .= 1
  nothing
end
callback_matrix = ContinuousCallback(condition_matrix, affect_matrix!)

for alg in CACHE_TEST_ALGS
  OrdinaryDiffEq.isimplicit(alg) && continue # this restriction should be removed in the future
  @show alg
  local sol = solve(prob_matrix, alg,callback=callback_matrix, dt=1/2)
  @test size(sol[end]) == (2,3)
end

# additional cache tests to find more bugs
println("Additional resize! checks")
u0resize3 = ones(2)
fresize3 = (du, u, p, t) -> du .= u
prob_resize3 = ODEProblem(fresize3, u0resize3, (0.0, 2.0))
condition_resize3 = (u, t, integrator) -> t - 1
affect!_resize3 = function (integrator)
  resize!(integrator, 3)
  integrator.u .= 1
  nothing
end
callback_resize3 = ContinuousCallback(condition_resize3, affect!_resize3)

for alg in CACHE_TEST_ALGS
  (OrdinaryDiffEq.isimplicit(alg) || OrdinaryDiffEq.alg_order(alg) < 2) && continue
  @show alg
  local sol = solve(prob_resize3, alg, callback=callback_resize3, dt=0.125)
  @test size(sol[end]) == (3,)
  @test all(sol[end] .== sol[end][1])
  @test sol[end][1] ≈ exp(1) atol=1.0e-2
end

# additional cache tests to find more bugs
println("Enforced adaptive dt checks")
u0_adapt = zeros(2)
f_adapt = (du, u, p, t) -> du .= t
prob_adapt = ODEProblem(f_adapt, u0_adapt, (0.0, 2.0))
condition_adapt = (u, t, integrator) -> true
affect!_adapt = function (integrator)
  dt = rand((0.0625, 0.125, 0.25))
  set_proposed_dt!(integrator, dt)
  integrator.opts.dtmax = dt
  integrator.dtcache = dt
  u_modified!(integrator, false)
  nothing
end
callback_adapt = DiscreteCallback(condition_adapt, affect!_adapt, save_positions=(false,false))

for alg in CACHE_TEST_ALGS
  (OrdinaryDiffEq.isimplicit(alg) || OrdinaryDiffEq.alg_order(alg) < 2) && continue
  @show alg
  local sol = solve(prob_adapt, alg, callback=callback_adapt, dt=0.125)
  @test all(idx -> all(isapprox.(sol.u[idx], 0.5*sol.t[idx]^2, atol=1.0e-6)), eachindex(sol.t))
end

# Force switching

function f2(du,u,p,t)
  @assert length(u) == length(du) "length(u) = $(length(u)), length(du) = $(length(du)) at time $(t)"
  for i in 1:length(u)
    if t > 10
      du[i] = -10000*u[i]
    else
      du[i] = 0.3*u[i]
    end
  end
  return du
end

function condition2(u, t, integrator)
  1-maximum(u)
end

function affect2!(integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition2,affect2!)
u0 = [0.2]
tspan = (0.0,20.0)
prob = ODEProblem(f2,u0,tspan)
sol = solve(prob, AutoTsit5(Rosenbrock23()), callback=callback)
@test length(sol[end]) > 1