text:
@pointwise_diagnostic(
    AtmosGCMConfigType,
    u,
    "m s^-1",
    "zonal wind",
    "eastward_wind",
    true,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[1] / states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    v,
    "m s^-1",
    "meridional wind",
    "northward_wind",
    true,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[2] / states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    w,
    "m s^-1",
    "vertical wind",
    "upward_air_velocity",
    true,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3] / states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    rho,
    "kg m^-3",
    "air density",
    "air_density",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    et,
    "J kg^-1",
    "total specific energy",
    "specific_dry_energy_of_air",
) do (
    energy::TotalEnergyModel,
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    states.prognostic.energy.ρe / states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    temp,
    "K",
    "air temperature",
    "air_temperature",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    air_temperature(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    pres,
    "Pa",
    "air pressure",
    "air_pressure",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    air_pressure(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    thd,
    "K",
    "dry potential temperature",
    "air_potential_temperature",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    dry_pottemp(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    ei,
    "J kg^-1",
    "specific internal energy",
    "internal_energy",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    internal_energy(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    ht,
    "J kg^-1",
    "specific enthalpy based on total energy",
    "",
) do (
    energy::TotalEnergyModel,
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    e_tot = states.prognostic.energy.ρe / states.prognostic.ρ
    total_specific_enthalpy(ts, e_tot)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    hi,
    "J kg^-1",
    "specific enthalpy based on internal energy",
    "atmosphere_enthalpy_content",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    specific_enthalpy(ts)
end

#= TODO
@XXX_diagnostic(
    "vort",
    AtmosGCMConfigType,
    GridInterpolated,
    "s^-1",
    "vertical component of relative velocity",
    "atmosphere_relative_velocity",
) do (atmos::AtmosModel, states::States, curr_time, cache)
end
=#

@pointwise_diagnostic(
    AtmosGCMConfigType,
    qt,
    "kg kg^-1",
    "mass fraction of total water in air (qv+ql+qi)",
    "mass_fraction_of_water_in_air",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    states.prognostic.moisture.ρq_tot / states.prognostic.ρ
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    ql,
    "kg kg^-1",
    "mass fraction of liquid water in air",
    "mass_fraction_of_cloud_liquid_water_in_air",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    liquid_specific_humidity(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    qv,
    "kg kg^-1",
    "mass fraction of water vapor in air",
    "specific_humidity",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    vapor_specific_humidity(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    qi,
    "kg kg^-1",
    "mass fraction of ice in air",
    "mass_fraction_of_cloud_ice_in_air",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    ice_specific_humidity(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    thv,
    "K",
    "virtual potential temperature",
    "virtual_potential_temperature",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    virtual_pottemp(ts)
end

@pointwise_diagnostic(
    AtmosGCMConfigType,
    thl,
    "K",
    "liquid-ice potential temperature",
    "",
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    liquid_ice_pottemp(ts)
end


text:
@testset "ridge.jl" begin
    A = rand(3, 3)
    b = rand(3)
    x = ridge(A, b)
    @assert x ≈ A \ b
end

ext: jl


text:
module TestGMMMath

using Compat
using Statistics
using Test
import TreeParzen: GMM

# GMM1 Math
weights = [0.1, 0.3, 0.4, 0.2]
mus = [1.0, 2.0, 3.0, 4.0]
sigmas = [0.1, 0.4, 0.8, 2.0]

# No low or high
samples = GMM.GMM1(weights, mus, sigmas, 10_001)
samples = sort(samples)
edges = samples[1:500:end]
pdf = exp.(GMM.GMM1_lpdf(edges[1:end - 1], weights, mus, sigmas))
dx = edges[2:end] .- edges[1:end - 1]
y = 1 ./ dx ./ length(dx)
err = (pdf .- y) .^ 2
@test maximum(err) < .1
@test mean(err) < .01
@test median(err) < .01

# Low and high
samples = GMM.GMM1(weights, mus, sigmas, 2.5, 3.5, 10_001)
samples = sort(samples)
edges = samples[1:500:end]
pdf = exp.(GMM.GMM1_lpdf(edges[1:end - 1], weights, mus, sigmas, 2.5, 3.5))
dx = edges[2:end] .- edges[1:end - 1]
y = 1 ./ dx ./ length(dx)
err = (pdf .- y) .^ 2
@test maximum(err) < 0.1
@test mean(err) < 0.01
@test median(err) < 0.01

# Low is >= high
@test_throws ArgumentError GMM.GMM1(weights, mus, sigmas, 3.5, 3.5, 10_001)

# QGMM1 Math
weights_t = [0.1, 0.3, 0.4, 0.2]
mus_t = [1.0, 2.0, 3.0, 4.0]
sigmas_t = [0.1, 0.4, 0.8, 2.0]
n_samples_t = 1_001

function test_samples(samples, c)
    @test samples == Int.(samples)
    samples_min = minimum(samples)
    samples_max = maximum(samples)
    bincount = samples .- samples_min
    counts = [count(x -> x == i, bincount) for i in 0:maximum(bincount)]
    @test sum(counts) == c.n_samples
    xcoords = @compat range(samples_min, samples_max; length = length(counts)) * c.q
    prob = if :low in propertynames(c)
        exp.(GMM.GMM1_lpdf(xcoords |> collect, c.weights, c.mus, c.sigmas, c.low, c.high, c.q))
    else
        exp.(GMM.GMM1_lpdf(xcoords |> collect, c.weights, c.mus, c.sigmas, c.q))
    end
    y = counts ./ c.n_samples
    err = (prob .- y) .^ 2

    @test maximum(err) < 0.1
    @test mean(err) < 0.01
    @test median(err) < 0.01
end

for c in (
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 1.0, n_samples = n_samples_t),
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 2.0, n_samples = n_samples_t),
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 0.5, n_samples = n_samples_t),
)
    samples = GMM.GMM1(c.weights, c.mus, c.sigmas, c.q, c.n_samples) / c.q
    test_samples(samples, c)
end

for c in (
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 1.0, low = 2.0, high = 4.0,
    n_samples = n_samples_t),
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 2.0, low = 2.0, high = 4.0,
    n_samples = n_samples_t),
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 1.0, low = 1.0, high = 4.1,
    n_samples = n_samples_t),
    (weights = weights_t, mus = mus_t, sigmas = sigmas_t, q = 2.0, low = 1.0, high = 4.1,
    n_samples = n_samples_t),
    (weights = [0.14285714, 0.28571429, 0.28571429, 0.28571429], mus = [5.505, 7., 2., 10.],
    sigmas = [8.99, 5., 8., 8.], q = 1.0, low = 1.01, high = 10.0, n_samples = 10_000),
    (weights = [0.33333333, 0.66666667], mus = [5.505, 5.], sigmas = [8.99, 5.19], q = 1.0,
    low = 1.01, high = 10.0, n_samples = 10_000),
)
    samples = GMM.GMM1(c.weights, c.mus, c.sigmas, c.low, c.high, c.q, c.n_samples) / c.q
    test_samples(samples, c)
end

end
true


text:
export LSTM

mutable struct LSTM
    insize::Int
    hsize::Int
    nlayers::Int
    droprate::Float64
    bidir::Bool
    weights::Vector
end

doc"""
    LSTM(::Type{T}, insize::Int, hsize::Int, nlayers::Int, droprate::Float64, bidir::Bool,
        [init_w=Orthogonal(), init_u=Orthogonal(), init_b=Fill(0)])

Long Short-Term Memory network.

```math
\begin{align*}
\mathbf{f}_{t} & =\sigma_{g}(W_{f}\mathbf{x}_{t}+U_{f}\mathbf{h}_{t-1}+\mathbf{b}_{f})\\
\mathbf{i}_{t} & =\sigma_{g}(W_{i}\mathbf{x}_{t}+U_{i}\mathbf{h}_{t-1}+\mathbf{b}_{i})\\
\mathbf{o}_{t} & =\sigma_{g}(W_{o}\mathbf{x}_{t}+U_{o}\mathbf{h}_{t-1}+\mathbf{b}_{o})\\
\mathbf{c}_{t} & =\mathbf{f}_{t}\odot\mathbf{c}_{t-1}+\mathbf{i}_{t}\odot\sigma_{c}(W_{c}\mathbf{x}_{t}+U_{c}\mathbf{h}_{t-1}+\mathbf{b}_{c})\\
\mathbf{h}_{t} & =\mathbf{o}_{t}\odot\sigma_{h}(\mathbf{c}_{t})
\end{align*}
```

* ``x_t \in R^{d}``: input vector to the LSTM block
* ``f_t \in R^{h}``: forget gate's activation vector
* ``i_t \in R^{h}``: input gate's activation vector
* ``o_t \in R^{h}``: output gate's activation vector
* ``h_t \in R^{h}``: output vector of the LSTM block
* ``c_t \in R^{h}``: cell state vector
* ``W \in R^{h \times d}``, ``U \in R^{h \times h}`` and ``b \in R^{h}``: weight matrices and bias vectors
* ``\sigma_g``: sigmoid function
* ``\sigma_c``: hyperbolic tangent function
* ``\sigma_h``: hyperbolic tangent function

```julia
T = Float32
x = Var(rand(T,100,10))
f = LSTM(T, 100, 100)
h = f(x)
```
"""
function LSTM(::Type{T}, insize::Int, hsize::Int, nlayers::Int, droprate::Float64, bidir::Bool;
    init_W=Normal(0,0.001), init_U=Orthogonal(), init_b=Fill(0), init_h=Fill(0), init_c=Fill(0)) where T

    weights = []
    coef = bidir ? 2 : 1
    for l = 1:nlayers
        for _ = 1:coef
            Ws = []
            Us = []
            for i = 1:4
                s = l == 1 ? insize : hsize*coef
                push!(Ws, init_W(T,s,hsize))
                push!(Us, init_U(T,hsize,hsize))
            end
            W = parameter(cat(Ws...,dims=2))
            U = parameter(cat(Us...,dims=2))
            b = parameter(init_b(T,4hsize))
            h = parameter(init_h(T,hsize))
            c = parameter(init_c(T,hsize))
            push!(weights, W, U, b, h, c)
        end
    end
    LSTM(insize, hsize, nlayers, droprate, bidir, weights)
end

function (f::LSTM)(x, dims)
    p = (insize=f.insize,hsize=f.hsize,nlayers=f.nlayers,droprate=f.droprate,bidir=f.bidir)
    lstm(x, dims, p, f.weights...)
end

function lstm(x::Var, dims, p, weights::Var...)
    @assert sum(dims) == size(x,2)
    @assert issorted(dims, rev=true)
    if isa(x.data, Array)
        lstm_cpu(x, dims, p, weights)
    elseif isa(x.data, CuArray)
        lstm_cuda(x, dims, p, weights)
    else
        throw("Invalid backend.")
    end
end
lstm(x::Node, args...) = Node(lstm, (x,args...))

function lstm_cpu(x::Var, dims, p, weights)
    h = x
    #coef = p.bidir ? 2 : 1
    i = 0
    for l = 1:p.nlayers
        #i = (l-1) * coef + 1
        #p = lstm.params[i]
        h1 = lstm_tstep(h, dims, weights[i+1:i+5]..., false)
        i += 5
        if p.bidir
            #p = lstm.params[i+1]
            h2 = lstm_tstep(h, dims, weights[i+1:i+5]..., true)
            h = concat(1, h1, h2)
            i += 5
        else
            h = h1
        end
    end
    h
end

function lstm_tstep(x::Var, dims, W::Var, U::Var, b::Var, h::Var, c::Var, rev::Bool)
    WU = concat(1, W, U)
    cumdims = Array{Int}(undef, length(dims)+1)
    cumdims[1] = 1
    for i = 1:length(dims)
        cumdims[i+1] = cumdims[i] + dims[i]
    end

    hsize = length(h)
    ht = concat(2, [h for i=1:length(dims)]...)
    ct = concat(2, [c for i=1:length(dims)]...)
    hts = Array{Var}(undef, size(x,2))
    cts = Array{Var}(undef, size(x,2))
    for t = 1:dims[1]
        xts = Var[]
        for j = 1:length(dims)
            d = dims[j]
            t > d && break
            k = cumdims[j]
            k += rev ? d-t : t-1
            push!(xts, x[:,k:k])
        end
        xt = concat(2, xts...)
        if size(ht,2) > size(xt,2)
            ht = ht[:,1:size(xt,2)]
            ct = ct[:,1:size(xt,2)]
        end
        ht, ct = lstm_onestep(xt, WU, b, ht, ct)
        for j = 1:length(dims)
            d = dims[j]
            t > d && break
            k = cumdims[j]
            k += rev ? d-t : t-1
            hts[k] = ht[:,j:j]
            cts[k] = ct[:,j:j]
        end
    end
    concat(2, hts...)
end

function lstm_onestep(xt::Var, WU::Var, b::Var, ht::Var, ct::Var)
    a = linear(concat(1,xt,ht), WU, b)
    n = size(a,1) ÷ 4
    i = sigmoid(a[1:n,:])
    f = sigmoid(a[n+1:2n,:])
    ct = (f .* ct) + (i .* tanh(a[2n+1:3n,:]))
    o = sigmoid(a[3n+1:4n,:])
    ht = o .* tanh(ct)
    ht, ct
end

function lstm_cuda(x::Var, dims, p, weights)
    Ws, hs, cs = [], [], []
    for i = 1:5:length(weights)
        W, U, b, h, c = weights[i:i+4]
        push!(Ws, vec(W.data), vec(U.data))
        for _ = 1:length(dims)
            push!(hs, h.data)
            push!(cs, c.data)
        end
    end
    for i = 1:5:length(weights)
        W, U, b, h, c = weights[i:i+4]
        push!(Ws, b.data, fill!(similar(b.data),0)) # b
    end

    #for (W,U,b,h,c) in lstm.params
    #    push!(Ws, vec(W.data), vec(U.data))
    #end
    #for (W,U,b,h,c) in lstm.params
    #    push!(Ws, b.data, fill!(similar(b.data),0))
    #    for i = 1:length(dims)
    #        push!(hs, h.data)
    #        push!(cs, c.data)
    #    end
    #end
    W = cat(Ws..., dims=1)
    h = cat(hs..., dims=1)
    c = cat(cs..., dims=1)

    t_xdata, t_dims = transpose_batch(x.data, dims)
    dir = p.bidir ? CUDNN.CUDNN_BIDIRECTIONAL : CUDNN.CUDNN_UNIDIRECTIONAL
    mode = CUDNN.CUDNN_LSTM
    t_ydata, work = CUDNN.rnn(p.insize, p.hsize, p.nlayers, p.droprate, dir, mode,
        t_xdata, t_dims, h, c, W, istrain())
    ydata, _ = transpose_batch(t_ydata, t_dims)
    Var(ydata, ∇lstm_cuda!, (x,dims,p,weights,W,work))
end

function ∇lstm_cuda!(y::Var, x::Var, dims, p, weights, dW, work)
    @assert !isnothing(x.grad)
    t_gy, t_dims = transpose_batch(y.grad, dims)
    t_gx, gh, gc = CUDNN.∇rnn_data(t_gy, work) # this call is required for ∇rnn_weights!
    gx, _ = transpose_batch(t_gx, t_dims)
    addto!(x.grad, gx)

    gW = fill!(similar(dW), 0)
    CUDNN.∇rnn_weights!(gW, work)
    Wi = 0
    hi = 0
    for i = 1:5:length(weights)
        W, U, b, h, c = weights[i:i+4]
        addto!(W.grad, 1, gW, Wi+1, length(W))
        Wi += length(W)
        addto!(U.grad, 1, gW, Wi+1, length(U))
        Wi += length(U)
        for _ = 1:length(dims)
            addto!(h.grad, 1, gh, hi+1, length(h))
            addto!(c.grad, 1, gc, hi+1, length(c))
            hi += length(h)
        end
    end
    for i = 1:5:length(weights)
        W, U, b, h, c = weights[i:i+4]
        addto!(b.grad, 1, gW, Wi+1, length(b))
        Wi += 2length(b)
    end
end

@generated function transpose_batch(x::CuMatrix{T}, dims::Vector{Int}) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void transpose_batch(int n, $Ct *t_x, int *t_cumdims, $Ct *x, int *cumdims) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= n*t_cumdims[1]*cumdims[1]) return;

        int vj = idx / n;
        int vi = idx - vj * n;
        int j = vj / t_cumdims[1];
        int i = vj - j * t_cumdims[1];
        if (t_cumdims[j] + i < t_cumdims[j+1]) {
            int t_idx = (t_cumdims[j] + i) * n + vi;
            idx = (cumdims[i] + j) * n + vi;
            t_x[t_idx] = x[idx];
        }
    }""")
    quote
        t_x = similar(x)
        t_dims = transpose_dims(dims)
        cumdims = Array{Cint}(undef, length(dims)+1)
        t_cumdims = Array{Cint}(undef, length(t_dims)+1)
        cumdims[1] = 0
        t_cumdims[1] = 0
        for i = 2:length(cumdims)
            cumdims[i] = cumdims[i-1] + dims[i-1]
        end
        for i = 2:length(t_cumdims)
            t_cumdims[i] = t_cumdims[i-1] + t_dims[i-1]
        end
        cumdims = CuArray(cumdims)
        t_cumdims = CuArray(t_cumdims)

        gdims, bdims = cudims(size(x,1)*dims[1]*t_dims[1])
        $k(gdims, bdims, size(x,1), pointer(t_x), pointer(t_cumdims), pointer(x), pointer(cumdims))
        t_x, t_dims
    end
end

function transpose_dims(dims::Vector{Int})
    k = length(dims)
    t_dims = Int[]
    for t = 1:dims[1]
        while dims[k] < t
            k -= 1
        end
        push!(t_dims, k)
    end
    t_dims
end


text:
"""
    MetropolisHastings{D}

`MetropolisHastings` has one field, `proposal`.
`proposal` is a `Proposal`, `NamedTuple` of `Proposal`, or `Array{Proposal}` in the shape of your data.
For example, if you wanted the sampler to return a `NamedTuple` with shape

```julia
x = (a = 1.0, b=3.8)
```

The proposal would be

```julia
proposal = (a=Proposal(Static(), Normal(0,1)), b=Proposal(Static(), Normal(0,1)))
````

Other allowed proposal styles are

```
p1 = Proposal(Static(), Normal(0,1))
p2 = Proposal(Static(), [Normal(0,1), InverseGamma(2,3)])
p3 = Proposal(Static(), (a=Normal(0,1), b=InverseGamma(2,3)))
p4 = Proposal(Static(), (x=1.0) -> Normal(x, 1))
```

The sampler is constructed using

```julia
spl = MetropolisHastings(proposal)
```

When using `MetropolisHastings` with the function `sample`, the following keyword
arguments are allowed:

- `init_params` defines the initial parameterization for your model. If
none is given, the initial parameters will be drawn from the sampler's proposals.
- `param_names` is a vector of strings to be assigned to parameters. This is only
used if `chain_type=Chains`.
- `chain_type` is the type of chain you would like returned to you. Supported
types are `chain_type=Chains` if `MCMCChains` is imported, or
`chain_type=StructArray` if `StructArrays` is imported.
"""
mutable struct MetropolisHastings{D} <: Metropolis
    proposal :: D
end

StaticMH(d) = MetropolisHastings(Proposal(Static(), d))
RWMH(d) = MetropolisHastings(Proposal(RandomWalk(), d))

# Propose from a vector of proposals
function propose(
    spl::MetropolisHastings{<:AbstractArray},
    model::DensityModel
)
    proposal = map(i -> propose(spl.proposal[i], model), 1:length(spl.proposal))
    return Transition(model, proposal)
end

function propose(
    spl::MetropolisHastings{<:AbstractArray},
    model::DensityModel,
    params_prev::Transition
)
    proposal = map(i -> propose(spl.proposal[i], model, params_prev.params[i]), 1:length(spl.proposal))
    return Transition(model, proposal)
end

# Make a proposal from one Proposal struct.
function propose(
    spl::MetropolisHastings{<:Proposal},
    model::DensityModel
)
    proposal = propose(spl.proposal, model)
    return Transition(model, proposal)
end

function propose(
    spl::MetropolisHastings{<:Proposal},
    model::DensityModel,
    params_prev::Transition
)
    proposal = propose(spl.proposal, model, params_prev.params)
    return Transition(model, proposal)
end

# Make a proposal from a NamedTuple of Proposal.
function propose(
    spl::MetropolisHastings{<:NamedTuple},
    model::DensityModel
)
    proposal = _propose(spl.proposal, model)
    return Transition(model, proposal)
end

@generated function _propose(
    proposals::NamedTuple{names},
    model::DensityModel
) where {names}
    expr = Expr(:tuple)
    map(names) do f
        push!(expr.args, Expr(:(=), f, :(propose(proposals.$f, model)) ))
    end
    return expr
end

function propose(
    spl::MetropolisHastings{<:NamedTuple},
    model::DensityModel,
    params_prev::Transition
)
    proposal = _propose(spl.proposal, model, params_prev.params)
    return Transition(model, proposal)
end

@generated function _propose(
    proposals::NamedTuple{names},
    model::DensityModel,
    params_prev::NamedTuple
) where {names}
    expr = Expr(:tuple)
    map(names) do f
        push!(expr.args, Expr(:(=), f, :(propose(proposals.$f, model, params_prev.$f)) ))
    end
    return expr
end

# Evaluate the likelihood of t conditional on t_cond.
function q(
    spl::MetropolisHastings{<:AbstractArray},
    t::Transition,
    t_cond::Transition
)
    return sum(map(i -> q(spl.proposal[i], t.params[i], t_cond.params[i]), 1:length(spl.proposal)))
end

function q(
    spl::MetropolisHastings{<:Proposal},
    t::Transition,
    t_cond::Transition
)
    return q(spl.proposal, t.params, t_cond.params)
end

function q(
    spl::MetropolisHastings{<:NamedTuple},
    t::Transition,
    t_cond::Transition
)
    ks = keys(t.params)
    total = 0.0

    for k in ks
        total += q(spl.proposal[k], t.params[k], t_cond.params[k])
    end

    return total
end

# Define the first step! function, which is called at the
# beginning of sampling. Return the initial parameter used
# to define the sampler.
function step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    N::Integer;
    init_params=nothing,
    kwargs...
)
    if init_params === nothing
        return propose(spl, model)
    else
        return Transition(model, init_params)
    end
end

# Define the other step functions. Returns a Transition containing
# either a new proposal (if accepted) or the previous proposal
# (if not accepted).
function step!(
    rng::AbstractRNG,
    model::DensityModel,
    spl::MetropolisHastings,
    ::Integer,
    params_prev::Transition;
    kwargs...
)
    # Generate a new proposal.
    params = propose(spl, model, params_prev)

    # Calculate the log acceptance probability.
    α = logdensity(model, params) - logdensity(model, params_prev) +
        q(spl, params_prev, params) - q(spl, params, params_prev)

    # Decide whether to return the previous params or the new one.
    if log(rand(rng)) < min(0.0, α)
        return params
    else
        return params_prev
    end
end


text:
include("lib/utils.jl")
include("lib/json.jl")

function is_ident_char(c)
  ('a' <= c && c <= 'z') ||
  ('0' <= c && c <= '9') ||
  c == '_'
end

function match_ident(str)
  bi = 1 # byte index

  while bi <= lastindex(str)
    if is_ident_char(str[bi])
      bi = nextind(str, bi)
    else
      break
    end
  end

  bi - 1
end

function match_sym(str)
  if str[1:2] in ["==", "!="]
    return 2
  end

  if occursin(str[1], "(){};=,+*")
    1
  else
    0
  end
end

function match_comment(str)
  if str[1:2] != "//"
    return 0
  end

  find_index(str, '\n')
end

function is_kw(str)
  str in [
    "func",
    "set",
    "var",
    "call",
    "call_set",
    "return",
    "while",
    "case"
  ]
end

function print_token(lineno, kind, str)
  json_print([lineno, kind, str], false)
  print("\n")
end

function lex()
  src = read_stdin_all()

  bi = 1
  lineno = 1

  while bi <= lastindex(src)
    rest = src[bi:end]

    if rest[1] == '\n'
      bi += 1
      lineno += 1
    elseif rest[1] == ' '
      bi += 1
    elseif 0 < match_sym(rest)
      # byte size
      bsize = match_sym(rest)
      str = rest[1:bsize]
      print_token(lineno, "sym", str)
      bi += bsize
    elseif 0 < match_comment(rest)
      bsize = match_comment(rest)
      str = rest[1:bsize]
      bi += bsize
    elseif 0 < match_int(rest)
      bsize = match_int(rest)
      str = rest[1:bsize]
      print_token(lineno, "int", str)
      bi += bsize
    elseif 0 <= match_str(rest)
      bsize = match_str(rest)
      str = rest[2:prevind(rest, bsize + 2)]
      print_token(lineno, "str", str)
      bi += bsize + 2
    elseif 0 < match_ident(rest)
      bsize = match_ident(rest)
      str = rest[1:bsize]
      if is_kw(str)
        print_token(lineno, "kw", str)
      else
        print_token(lineno, "ident", str)
      end
      bi += bsize
    else
      error("unexpected pattern line($(lineno)) rest($(rest))")
    end
  end
end

# --------------------------------

lex()


text:
function maprange(s, a, b)
    a₁, a₂ = minimum(a), maximum(a)
    b₁, b₂ = minimum(b), maximum(b)
    return b₁ + (s - a₁) * (b₂ - b₁) / (a₂ - a₁)
end

@show maprange(6, 1:10, -1:0)
@show maprange(0:10, 0:10, -1:0)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type Medium end

struct MediumFlat<:Medium

    minBoundaryType::String
    maxBoundaryType::String
    N::Int

end

function MediumFlat(t::String,N::Int)

    if t == "Periodic"
        return MediumFlat("Periodic","Periodic",N)
    elseif t == "Newmann"
        return MediumFlat("Newmann","Newmann",N)
    elseif t == "Dirichlet"
        return MediumFlat("Dirichlet","Dirichlet",N)
    else
        t = split(t,"-")
        for i in t
            if !(i in ["Periodic","Newmann","Dirichlet"])
                error("N boundary type implemented called ", i, " in ", t, ".")
            end
        end
        if length(t) != 2
            error("Boundary has to be declared as 'nameOfBoundary' or 'nameOfLowerBoundary-nameOfUpperBoundary'")
        end
        if "Periodic" in t && t[1] != t[2]
            error("If boundary is Periodic, upper and lower limits have to be the same.")
        end

        return MediumFlat(t[1],t[2],N)
    end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

BenchmarkTools.DEFAULT_PARAMETERS.samples = 5

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime gb = Groebner.groebner($system, reduced=false)
    # gb = Groebner.groebner(system, reduced=false)
end

function benchmark_system_my_certify(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime gb = Groebner.groebner($system, reduced=false, linalg=:prob)
    # gb = Groebner.groebner(system, reduced=false)
end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
        ("siwr", read_SIWR()),
        ("sear", read_SEAIJRC()),
        ("root 11", Groebner.rootn(11, ground=ground)),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("eco 10"    ,Groebner.eco10(ground=ground)),
        ("katsura 8"    ,Groebner.katsuran(8, ground=ground))
    ]

    println()
    for (name, system) in systems
        println("$name")
        # benchmark_system_my(system)
        benchmark_system_my_certify(system)
    end
end

ground = AbstractAlgebra.QQ
run_f4_ff_degrevlex_benchmarks(ground)

#=
root 11
  55.078 ms (183751 allocations: 23.39 MiB)
noon 7
  805.825 ms (1773797 allocations: 276.06 MiB)
eco 10
  911.046 ms (1324278 allocations: 294.65 MiB)
katsura 8
  3.124 s (1794620 allocations: 312.78 MiB)
=#


text:
α = 0.4
β = 0.96
μ = 0
s = 0.1

c1 = log(1 - α * β) / (1 - β)
c2 = (μ + α * log(α * β)) / (1 - α)
c3 = 1 / (1 - β)
c4 = 1 / (1 - α * β)

# Utility
u(c) = log(c)

u_prime(c) = 1 / c

# Deterministic part of production function
f(k) = k^α

f_prime(k) = α * k^(α - 1)

# True optimal policy
c_star(y) = (1 - α * β) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y)


text:
# https://docs.julialang.org/en/v1/manual/asynchronous-programming/

using FileWatching

function learn_basic_task_operations()
    # 3 ways to async

    func = function()       # function
        sleep(5)
        println("Done-1")
    end
    task1 = Task(func)      # create task from a function
    schedule(task1)         # run task

    task2 = @task begin     # create task from code block
        sleep(4)
        println("Done-2")
    end
    schedule(task2)         # run task

    task3 = @async begin    # task & run in 1 step
        sleep(3)
        println("Done-3")
    end

    #wait(task1)
end

function learn_Channel()
    order_chn = Channel{Int}(32)
    ship_chn = Channel{Tuple}(32)

    n = 10

    take_order = function()
        for order_id in 1:n
            println("Received order $order_id")
            put!(order_chn, order_id)
        end
    end

    manufacture = function()
        for order_id in order_chn

            manufacture_time = rand()
            sleep(manufacture_time)
            println("Order $order_id manufactured in $manufacture_time")
            put!(ship_chn, (order_id, manufacture_time))
        end
    end

    ship = function()
        for (order_id,manufacture_time) in ship_chn
            println("Shipped $order_id")
        end
    end

    @async ship()
    @async manufacture()
    @async take_order()
end

"""

@async does not seem to work, see same function in `learn_Threads.jl`
"""
function learn_write_read_file()
    n = 5000000
    done = Threads.Atomic{Bool}(false)

    writefile = function()
        open("output.txt", "w") do io
            for i in 1:n
                #sleep(0.1)
                println(io, "Hello! Number $i")
                #flush(io)
                # must use yield() for this to work
                yield()
            end
        end

        Threads.atomic_or!(done, true)
        println("Done writing")
    end

    readfile = function()
        while !isfile("output.txt")
            sleep(0.1)
        end

        pos = 0
        first_line = ""
        final_line = ""

        while !done[]
            sleep(0.2) # won't work without sleep
            open("output.txt", "r") do io
                skip(io, pos)

                first_line = readline(io)
                for line in eachline(io)
                    final_line = line
                end

                pos = position(io)
            end

            # When first_line is blank meaning no update on the file between
            # opening.
            println("First line: $first_line")
            println("Final line: $final_line")

            # no need for yield(), guessing because sleep() above
            #yield()
        end

        # do one more time here?
    end

    write_task = Task(writefile)
    read_task = Task(readfile)

    schedule(read_task)
    schedule(write_task)

    wait(read_task)
    wait(write_task)
end

#learn_basic_task_operations()
#learn_Channel()
learn_write_read_file()

nothing


text:
module Adapters # Octo

include(joinpath("Adapters", "SQLite.jl"))
include(joinpath("Adapters", "MySQL.jl"))
include(joinpath("Adapters", "PostgreSQL.jl"))
include(joinpath("Adapters", "SQL.jl"))
include(joinpath("Adapters", "ODBC.jl"))
include(joinpath("Adapters", "JDBC.jl"))
include(joinpath("Adapters", "Hive.jl"))

end # module Octo.Adapters


text:
function _make_discrete_gui()

    discrete_grid = Grid()
    discrete_space_button = SpinButton(2:100)
    discrete_grid[1,1] = discrete_space_button
    discrete_grid[2,1] = Label("Space Between Points")

    discrete_max_points_button = SpinButton(4:20)
    discrete_grid[1,2] = discrete_max_points_button
    discrete_grid[2,2] = Label("Max number of Points")

    discrete_auto_calc = CheckButton("Auto Calculate")
    setproperty!(discrete_auto_calc,:active,true)
    discrete_grid[1,3] = discrete_auto_calc

    discrete_add_button = CheckButton("Add Discrete Point")
    discrete_grid[1,4] = discrete_add_button

    discrete_delete_button = Button("Delete All Points in Frame")
    discrete_grid[1,5] = discrete_delete_button

    discrete_win=Window(discrete_grid)
    Gtk.showall(discrete_win)
    visible(discrete_win,false)

    d_widgets=discrete_widgets(discrete_win,discrete_space_button,discrete_max_points_button,discrete_auto_calc,
    discrete_add_button,discrete_delete_button)
end

function add_discrete_callbacks(w::discrete_widgets,handles::Tracker_Handles)

    signal_connect(discrete_distance_cb,w.space_button,"value-changed",Void,(),false,(handles,))
    signal_connect(discrete_points_cb,w.points_button,"value-changed",Void,(),false,(handles,))
    signal_connect(discrete_auto_cb,w.calc_button,"clicked",Void,(),false,(handles,))
    signal_connect(discrete_add_cb,w.add_button,"clicked",Void,(),false,(handles,))
    signal_connect(discrete_delete_cb,w.delete_button,"clicked",Void,(),false,(handles,))

    nothing
end

function discrete_distance_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    num_dist=getproperty(han.d_widgets.space_button,:value,Int)

    han.d_spacing = num_dist

    make_discrete_woi(han.wt,han.woi,han.tracked,num_dist)

    redraw_all(han)

    nothing
end

function discrete_points_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    num_points=getproperty(han.d_widgets.points_button,:value,Int)
    num_dist=getproperty(han.d_widgets.space_button,:value,Int)

    change_discrete_size(han.wt,num_points)
    make_discrete_woi(han.wt,han.woi,han.tracked,num_dist)

    redraw_all(han)

    nothing

end

function discrete_auto_cb(w::Ptr, user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.discrete_auto_calc=getproperty(han.d_widgets.calc_button,:active,Bool)

    nothing
end

function discrete_add_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    if getproperty(han.d_widgets.add_button,:active,Bool)
        han.selection_mode = 13
    else
        han.selection_mode = 1
        determine_viewers(han)
    end

    nothing
end

function discrete_delete_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    #han.wt.w_p[:,han.frame] .= 0.0
    redraw_all(han)

    nothing
end

#=
function add_discrete_point(han::Tracker_Handles,x,y)

    num_points = div(size(han.wt.w_p,1),2)
    ss=han.d_spacing

    last_point=1

    loc_x=zeros(Float64,num_points)
    loc_y=zeros(Float64,num_points)
    point_exists=falses(num_points)
    for i=1:num_points

        if han.wt.w_p[i*2-1,han.frame] != 0.0
            loc_x[i]=han.wt.w_p[i*2-1,han.frame]
            loc_y[i]=han.wt.w_p[i*2,han.frame]
            point_exists[i] = true
            last_point = i
        end
    end

    if last_point < num_points

        total_dist=0.0
        for i=2:last_point
            total_dist += sqrt((loc_x[i]-loc_x[i-1])^2 + (loc_y[i]-loc_y[i-1])^2)
        end

        true_dist=sqrt((x - loc_x[last_point])^2 + (y-loc_y[last_point])^2)

        pos_vals=((last_point+1):(num_points))*ss

        j=argmin(abs.(pos_vals .- true_dist))

        han.wt.w_p[(last_point+j)*2-1,han.frame] = x
        han.wt.w_p[(last_point+j)*2,han.frame] = y
    end

    nothing
end
=#


text:
function mkdge2dgf(facecon,entmax)

#[~,nf] = size(facecon);
ndgf = length(facecon);

ent = sort(unique(facecon[:]));
ent = ent[ent .>0 ];
ndof = length(ent);
#entmax = ent(end); # largest entity

# entity-to-index mapping
ent2ind = zeros(Int,entmax,1);
ent2ind[ent] = 1:ndof;
 # store number of nfighboring elements for each entity
rowdge2dgf = zeros(Int,ndof,1);
for i = 1:ndgf  # for each face entity
    k = facecon[i];     # entities on element i
    if k>0
        ind = ent2ind[k];   # get entity indices
        rowdge2dgf[ind] = rowdge2dgf[ind] + 1;
    end
end
rowdge2dgf=[0; cumsum(rowdge2dgf[:])];

 # store nfighboring-element indices for each entity
coldge2dgf = zeros(Int,rowdge2dgf[end],1);
inc = ones(Int,ndof,1);
for i = 1:ndgf
    k = facecon[i];     # entities on element i
    if k>0
        ind = ent2ind[k];   # get entity indices
        # rowdge2dgf(ind): pointer to the list of entities k
        coldge2dgf[rowdge2dgf[ind]+inc[ind]] = i;
        inc[ind] = inc[ind] + 1;
    end
end

return rowdge2dgf,coldge2dgf,ent2ind,ent

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Ported https://github.com/bpalmer4 Python implementation
# Reference:
# - https://www.abs.gov.au/websitedbs/d3310114.nsf/4a256353001af3ed4b2562bb00121564/5fc845406def2c3dca256ce100188f8e!OpenDocument

# Henderson symmetric weights
function hmaSymmetricWeights(n::Int)
    m = (n-1)÷2
    m1 = (m+1)^2
    m2 = (m+2)^2
    m3 = (m+3)^2
    d = BigInt(8*(m+2)*(m2-1)*(4*m2-1)*(4*m2-9)*(4*m2-25))

    w = map(x -> round(Float16(BigFloat(315*(m1-(x^2))*(m2-(x^2))*(m3-(x^2))*(3*m2-11*(x^2)-16))/d), digits = 3), range(0, m+1, step = 1))[end-1:-1:1]
    u = vcat(w, w[end-1:-1:1])
    return if mod(n, 2) != 0
        u
    else
        vcat(u, missing)
    end
end

# Henderson asymmetric weights for dataset edges
function hmaAsymmetricWeights(m::Int, w::AbstractArray{<:Number})
    n = length(w)

    @assert m <= n "The m argument must be less than w"
    @assert m >= (n-1)÷2 "The m argument must be greater than (n-1)/2"

    sumResidual = sum(w[range(m + 1, n, step = 1)]) / m
    sumEnd = sum(map(x -> (x - ((m+1.0)/2.0)) * w[x], range(m+1, n, step = 1)))

    ic = if n >= 13 && n < 15
        3.5
    elseif n >= 15
        4.5
    else 1.0
    end
    b2s2 = BigFloat(4/pi/ic^2)
    denominator = 1.0 + ((m*(m-1.0)*(m+1.0) / 12.0 ) * b2s2)

    u = map(r -> round(Float16(w[r] + sumResidual + (((r-1+1.0) - (m+1.0)/2.0) * b2s2) / denominator * sumEnd), digits = 3),
        range(1, m, step = 1))
    return (u)
end

"""
Package: Forecast

    hma(s, n)

Applies the Henderson moving average filter to dataset `s` with `n`-term.

Information about involved processes and application can be found at the following:

("Time Series Analysis: The Process of Seasonal Adjustment")[https://www.abs.gov.au/websitedbs/d3310114.nsf/4a256353001af3ed4b2562bb00121564/5fc845406def2c3dca256ce100188f8e!OpenDocument]
Australian Bureau of Statistics

# Arguments
- `s`: Observations' support.
- `n`: Observation values.

# Returns
An array of Henderson filter smoothed values provided in `s`.

# Examples
```julia-repl
julia> hma(rand(24), 13)
24-element Array{Float64,1}:
[...]
```
"""
function hma(s::AbstractArray{<:Number}, n::Int)
    @assert isodd(n) "n must be odd"
    @assert n >= 5 "n must be >= 5"
    @assert length(s) >= n "dataset must be >= than n"

    w = hmaSymmetricWeights(n)
    m = (n-1) ÷ 2
    l = length(s)
    map(i -> if i - 1 < m
            u = hmaAsymmetricWeights(m + i, w)[end:-1:1]
            sum(s[1:i + m] .* u)
        elseif i - 1 + m >= l
            u = hmaAsymmetricWeights(m + l - i + 1, w)
            sum(s[i-m:l] .* u)
        else
            sum(s[i-m:i+m] .* w)
        end, range(1, l, step = 1))
end


text:
function eject(Q::AbstractVector)
    if isempty(Q)
        return Q
    end

    pop = Q[1]
    deleteat!(Q, 1)
    return [pop]
end


text:
export plot_realtime_corr0t,plot_realtime_corr0t!, plot_full_corr0t,
       plot_Gpm_Gmp, plot_eucledian_corr0t, plot_xi_param_corr0t,
       plot_avg_avg2, plot_avg_avg2_nonEquil,
       plot_corrections, plot_avg2_with_corrections

####################################################

using Plots,LaTeXStrings
using StatsPlots
using Measures

M_SIZE = 1.8
SOL_LW = 0.8

plot_setup = Dict(:legend     => :topright, #(0.1, 0.8), #(1.08, 0.6),
                  :grid       => false,
                  :foreground_color_legend=> nothing,
                  :background_color_legend=> nothing,
                  :framestyle => :box,
                  :thickness_scaling => 1.5)

markers_dict(color,markershape) = Dict(
                        :color => color,
                        :markershape => markershape,
                        :markersize  => M_SIZE,
                        :markerstrokecolor => color,
                        :lw => false)

markers_dict(color) = markers_dict(color,:square)

solution_line_dict = Dict(:color => "black",:lw => SOL_LW)

color_1 = 1
color_2 = 2
color_3 = "black"
color_4 = 3

## PLOT Real-Time CORR0T
function plot_realtime_corr0t(corr0t_Re, corr0t_Im, CC;
        true_solution=0,true_solution2=0,label="",theta=0.1,
        save=false,fig_name="realtime_corrx0xt")

    x = real(CC)
    tpmax = argmax(x)

    fig = plot(xlabel=L"$x_0$",ylabel= true ? "" : L"$\langle \phi(0) \phi(t) \rangle$";plot_setup...)#,ylim=[-0.3,0.32])#,xlim=[0.3,0.501])#,ylim=[0.15,0.25],xlim=[0.3,0.501])

    if true_solution != 0
        tmax_inx = argmin(abs.(true_solution["x"] .- maximum(x[1:tpmax])))
        plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tRe"][1:tmax_inx],label=L"$\textrm{Solution}$", color=color_3, lw=SOL_LW)
        plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tIm"][1:tmax_inx],label=false, color=color_3, lw=SOL_LW)
    end

    if true_solution2 != 0
        tmax_inx = argmin(abs.(true_solution2["x"] .- maximum(x[1:tpmax])))
        plot!(fig,true_solution2["x"][1:tmax_inx],true_solution2["corr0tRe"][1:tmax_inx],label=L"$\textrm{Solution Tilted}$",color=color_4)
        plot!(fig,true_solution2["x"][1:tmax_inx],true_solution2["corr0tIm"][1:tmax_inx],label=false,color=color_4)
    end

    fig = plot_realtime_corr0t!(fig, corr0t_Re, corr0t_Im, CC; label=label,theta=theta)
    #lens!([0.39,0.505], [-0.205,-0.165], inset = (1,bbox(0.65, 0.4, 0.3, 0.3) ),xticks= nothing,yticks=nothing,framestyle = :box,subplot = 2)
    #lens!([0.39,0.505], [-0.205,-0.165], inset = (1,bbox(0.68, 0.4, 0.3, 0.3) ),xticks= nothing,yticks=nothing,framestyle = :box,subplot = 2)
    #lens!([0.41,0.505], [0.155,0.21], inset = (1, bbox(0.20, 0.13, 0.3, 0.3)), xticks= nothing,yticks=nothing, framestyle = :box, subplot=3)
    #lens!([0.41,0.505], [0.155,0.21], inset = (1, bbox(0.37, 0.22, 0.3, 0.3)), xticks= nothing,yticks=nothing, framestyle = :box, subplot=3)

    if save
        savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    return fig
end

function plot_realtime_corr0t!(fig,corr0t_Re, corr0t_Im, CC; label="",theta=0.1)
    x = real(CC)
    tpmax = argmax(x)

    scatter!(fig,x[1:tpmax],corr0t_Re[1:tpmax],label= (label=="" ? L"$\textrm{Re} \; G_{++}(x_0)$" : label);markers_dict(color_1,:square)...)#color=color_1, markershape=:square,markersize=M_SIZE,markerstrokecolor = color_1,lw=0.2)
    scatter!(fig,x[1:tpmax],corr0t_Im[1:tpmax],label= (label=="" ? L"$\textrm{Im} \; G_{++}(x_0)$" : label);markers_dict(color_2,:utriangle)...)#color=color_2, markershape=:square,markersize=M_SIZE,markerstrokecolor = color_2,lw=0.2)
    return fig
end

## PLOT Eucledian CORR0T
function plot_full_corr0t(corr0t_Re, corr0t_Im, CC;
    true_solution=0,label="",
    save=false,fig_name="realtime_corrx0xt0t_full")
    x = real(CC)
    tpmax = argmax(x)

    fig = plot(xlabel = L"$\xi$",
    ylabel = true ? "" : L"$ \langle \phi(0) \phi(\xi) \rangle - \langle \phi(\xi) \rangle\langle \phi(0) \rangle $";plot_setup...)#, ylim=[-0.25,0.32])
    tmax_inx = argmin(abs.(true_solution["x"] .- maximum(x[1:tpmax])))
    plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tRe"][1:tmax_inx],label=L"$\textrm{Solution}$",color=color_3, lw=SOL_LW)
    plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tIm"][1:tmax_inx],label=false,color=color_3, lw=SOL_LW)

    plot!(fig,true_solution["x"][tmax_inx] .+ true_solution["x"][1:tmax_inx],
    reverse(true_solution["corr0tRe"][1:tmax_inx]),label=false,color=color_3,lw=SOL_LW)
    plot!(fig,true_solution["x"][tmax_inx] .+ true_solution["x"][1:tmax_inx],
    reverse(true_solution["corr0tIm"][1:tmax_inx]),label=false,color=color_3,lw=SOL_LW)

    # Forward
    fw_interval = 1:tpmax
    scatter!(fig,x[fw_interval],corr0t_Re[fw_interval],
    label= (label=="" ? L"$\textrm{Re}\; G_{++}(\xi)$" : label); markers_dict(color_1)...)
    scatter!(fig,x[fw_interval],corr0t_Im[fw_interval],
    label= (label=="" ? L"$\textrm{Im}\; G_{++}(\xi)$" : label); markers_dict(color_2,:utriangle)...)

    # Backward
    bw_interval = tpmax+1:2*tpmax-1
    scatter!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),corr0t_Re[bw_interval],
    label= L"$\textrm{Re}\; G_{+-}(\xi)$";markers_dict(3,:dtriangle)...)
    scatter!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),corr0t_Im[bw_interval],
    label= L"$\textrm{Im}\; G_{+-}(\xi)$";markers_dict(4,:star)...)

    lens!([0.4,0.505], [-0.21,-0.17], inset = (1, bbox(0.55, 0.48, 0.25, 0.25)),xticks= nothing,yticks=nothing,framestyle = :box,subplot = 2)
    lens!([0.4,0.505], [0.15,0.21], inset = (1, bbox(0.7, 0.15, 0.25, 0.25)), xticks= nothing,yticks=nothing, framestyle = :box, subplot=3)

    if save
        savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    display(fig)
end

## PLOT Eucledian CORR0T
function plot_Gpm_Gmp(Gpm_Re, Gpm_Im, Gmp_Re, Gmp_Im, CC;
    true_solution=0,label="",
    save=false,fig_name="realtime_corrx0xt0t_full")

    x = real(CC)
    tpmax = argmax(x)

    fig = plot(xlabel = L"$x_0$",
    ylabel = true ? "" : L"$ \langle \phi(0) \phi(t) \rangle - \langle \phi(x_0) \rangle\langle \phi(0) \rangle $";plot_setup...)#, ylim=[-0.25,0.32])
    tmax_inx = argmin(abs.(true_solution["x"] .- maximum(x[1:tpmax])))
    plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tRe"][1:tmax_inx],label=L"$\textrm{Solution } G_{+-}$",color=color_3, lw=SOL_LW)
    plot!(fig,true_solution["x"][1:tmax_inx],true_solution["corr0tIm"][1:tmax_inx],label=false,color="black", lw=SOL_LW)
    plot!(fig,true_solution["x"][1:tmax_inx],(-1) .* true_solution["corr0tIm"][1:tmax_inx],label=L"$\textrm{Solution } G_{-+}$",color="gray", lw=SOL_LW)

    # Forward
    fw_interval = 1:tpmax
    scatter!(fig,x[fw_interval],Gpm_Re[fw_interval],
    label= (label=="" ? L"$\textrm{Re} \;G_{+-}(x_0)$" : label);markers_dict(1,:square)...)
    scatter!(fig,x[fw_interval],Gpm_Im[fw_interval],
    label= (label=="" ? L"$\textrm{Im} \;G_{+-}(x_0)$" : label);markers_dict(3,:square)...)

    # Forward
    fw_interval = 1:tpmax
    scatter!(fig,x[fw_interval],Gmp_Re[fw_interval],
    label= (label=="" ? L"$\textrm{Re} \; G_{-+}(x_0)$" : label);markers_dict(2,:dtriangle)...)
    scatter!(fig,x[fw_interval],Gmp_Im[fw_interval],
    label= (label=="" ? L"$\textrm{Im} \; G_{-+}(x_0)$" : label);markers_dict(4,:dtriangle)...)

    if save
        savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    display(fig)
end

## PLOT Eucledian CORR0T
function plot_eucledian_corr0t(corr0t_Re, corr0t_Im, CC; true_solution=0,
            save=false,fig_name="eucledian_corrx0xt")

    x = -1 .* imag(CC)
    fig = plot(xlabel=L"$\xi$", ylabel=true ? "" : L"$\langle \phi(0)\phi(\xi) \rangle - \langle \phi(0) \rangle \langle\phi(\xi) \rangle$";
    ylim=[-Inf,Inf],plot_setup...)
    plot!(1 .+ true_solution["x"],true_solution["corr0tRe"],label=L"$\textrm{Solution}$", color=color_3,lw=SOL_LW)
    plot!(fig,1 .+ true_solution["x"],true_solution["corr0tIm"],label=false, color=color_3,lw=SOL_LW)
    scatter!(fig,1 .+ vcat(x[ceil(Int64,length(x)/2)+1:end],[1]),vcat(corr0t_Re[ceil(Int64,length(x)/2)+1:end],corr0t_Re[1]),label=L"$\textrm{Re}\; G_E(\xi)$",color=color_1, markershape=:square,markersize=M_SIZE,markerstrokecolor = color_1,lw=0)
    scatter!(fig,1 .+ vcat(x[ceil(Int64,length(x)/2)+1:end],[1]),vcat(corr0t_Im[ceil(Int64,length(x)/2)+1:end],corr0t_Im[1]),label=L"$\textrm{Im}\; G_E(\xi)$",color=color_2, markershape=:utriangle,markersize=M_SIZE,markerstrokecolor = color_2,lw=0)

    if save
    savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    display(fig)
end

## PLOT corr0t along the xi contour parmeter
function plot_xi_param_corr0t(corr0t_Re, corr0t_Im, tp; true_solution=0)
    fig = plot(legend=:bottomright)#,ylim=[-0.3,0.4])
    #plot!(true_solution["x"],true_solution["corr0tRe"],label="ScroSol Re")
    #plot!(fig,true_solution["x"],true_solution["corr0tIm"],label="ScroSol Im")
    plot!(fig,tp,corr0t_Re,label=L"$Re(\phi(0)\phi(t))$",color=1, markershape=:square,markerstrokecolor = 1,lw=0.1)
    plot!(fig,tp,corr0t_Im,label=L"$Im(\phi(0)\phi(t))$",color=2, markershape=:square,markerstrokecolor = 2,lw=0.1)
end

## PLOT AVG and AVG2
function plot_avg_avg2(avg_Re, avg_Im, avg2_Re, avg2_Im, tp; true_solution=0,save=false,fig_name="param_corrx")

    fig = plot(xlabel=L"$\xi$"; plot_setup...)#,ylim=[-0.02,0.325])

    tpmax = argmax(tp)
    tmax_inx = argmin(abs.(true_solution["x"] .- maximum(tp[1:tpmax])))
    plot!(fig,2*true_solution["x"][1:tmax_inx],true_solution["avgRe"][1:tmax_inx],label=L"$\textrm{Solution}$"; solution_line_dict...)
    plot!(fig,2*true_solution["x"][1:tmax_inx],true_solution["avg2Re"][1:tmax_inx],label=false; solution_line_dict...)

    scatter!(fig,2*tp,avg_Re,label=L"$\textrm{Re}\langle \phi \rangle $"; markers_dict(1, :square)...)
    scatter!(fig,2*tp,avg_Im,label=L"$\textrm{Im}\langle \phi \rangle $"; markers_dict(2, :utriangle)...)
    scatter!(fig,2*tp,avg2_Re,label=L"$\textrm{Re}\langle \phi^2\rangle $"; markers_dict(3, :star)...)
    scatter!(fig,2*tp,avg2_Im,label=L"$\textrm{Im}\langle \phi^2 \rangle $"; markers_dict(4, :xcross)...)

    if save
        savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    display(fig)
end

## PLOT AVG and AVG
function plot_avg_avg2_nonEquil(avg_Re, avg_Im, avg2_Re, avg2_Im, CC; true_solution=0,save=false,fig_name="param_corrx")

    x = real(CC)
    tpmax = argmax(x)

    fig = plot(xlabel=L"$\xi$"; plot_setup...)

    x_true = true_solution["x"]
    if true_solution != 0
        tmax_inx = argmin(abs.(true_solution["x"] .- maximum(x[1:tpmax])))

        plot!(fig,x_true[1:tmax_inx],true_solution["avgRe"][1:tmax_inx],
        label=L"\textrm{Solution}"; solution_line_dict...)
        plot!(fig,x_true[1:tmax_inx],true_solution["avgIm"][1:tmax_inx],
        label=false;solution_line_dict...)

        plot!(fig,x_true[tmax_inx] .+ x_true[1:tmax_inx],
        reverse(true_solution["avgRe"][1:tmax_inx]),
        label=false;solution_line_dict...)
        plot!(fig,x_true[tmax_inx] .+ x_true[1:tmax_inx],
        reverse(true_solution["avgIm"][1:tmax_inx]),
        label=false;solution_line_dict...)

        plot!(fig,x_true[1:tmax_inx],true_solution["avg2Re"][1:tmax_inx],
        label=false; solution_line_dict...)
        plot!(fig,x_true[1:tmax_inx],true_solution["avg2Im"][1:tmax_inx],
        label=false;solution_line_dict...)

        plot!(fig,x_true[tmax_inx] .+ x_true[1:tmax_inx],
        reverse(true_solution["avg2Re"][1:tmax_inx]),
        label=false;solution_line_dict...)
        plot!(fig,x_true[tmax_inx] .+ x_true[1:tmax_inx],
        reverse(true_solution["avg2Im"][1:tmax_inx]),
        label=false;solution_line_dict...)
    end

    # Forward
    fw_interval = 1:tpmax
    scatter!(fig,x[fw_interval],avg_Re[fw_interval],
    label= L"$\textrm{Re}\langle \phi \rangle$"; markers_dict(1,:square)...)
    scatter!(fig,x[fw_interval],avg_Im[fw_interval],
    label= L"$\textrm{Im}\langle \phi \rangle$"; markers_dict(2,:utriangle)...)

    # Forward
    fw_interval = 1:tpmax
    scatter!(fig,x[fw_interval],avg2_Re[fw_interval],
    label= L"$\textrm{Re}\langle \phi^2 \rangle$"; markers_dict(3,:star)...)
    scatter!(fig,x[fw_interval],avg2_Im[fw_interval],
    label= L"$\textrm{Im}\langle \phi^2 \rangle$"; markers_dict(4,:cross)...)

    # Backward
    bw_interval = tpmax+1:2*tpmax-1
    plot!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),avg_Re[bw_interval],
    label= false; markers_dict(1,:square)...)
    plot!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),avg_Im[bw_interval],
    label= false; markers_dict(2,:utriangle)...)

    # Backward
    bw_interval = tpmax+1:2*tpmax-1
    plot!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),avg2_Re[bw_interval],
    label= false; markers_dict(3,:star)...)
    plot!(fig,x[tpmax] .+ (x[tpmax] .- x[bw_interval]),avg2_Im[bw_interval],
    label= false; markers_dict(4,:cross)...)

    if save
        savefig(fig, string("Imgs/",fig_name,".pdf"))
    end
    display(fig)
end

function plot_corrections(SS_Re, SS2_Re, SS_Im, SS2_Im, avg2_Re, avg2_Im, tp;true_solution=0, save=false)

    fig_Re = plot(xlabel="";plot_setup...,xticks=nothing,bottom_margin=-3mm)#,ylim=[-0.15,0.015])

    if true_solution != 0
        scatter!(fig_Re,2*tp[1:end-1],(avg2_Re .- true_solution[1] ); markers_dict(2,:square)...,
                 label=L"$\textrm{Re}\langle \Delta \phi^2 \rangle $")
    end

    scatter!(fig_Re,2*tp[1:end-1], SS2_Re; markers_dict(1,:circle)...,
             label=L"$\textrm{Re} \; \Sigma$")
    scatter!(fig_Re,2*tp[1:end-1],SS_Re .+ SS2_Re; markers_dict(5,:utriangle)...,
             label=L"$\textrm{Re}\left\langle \left(\tilde\phi - \phi\right)^2 \right\rangle $")

    fig_Im = plot(xlabel=L"$\xi$"; plot_setup...,top_margin=0mm)#,ylim=[-0.015,0.009])

    if true_solution != 0
        scatter!(fig_Im,2*tp[1:end-1],avg2_Im; markers_dict(4,:square)...,
                 label=L"$\textrm{Im}\langle \Delta \phi^2 \rangle $")
    end
    scatter!(fig_Im,2*tp[1:end-1],SS2_Im; markers_dict(3,:circle)...,
             label=L"$\textrm{Im} \; \Sigma$")
    scatter!(fig_Im,2*tp[1:end-1], SS_Im .+ SS2_Im; markers_dict(6,:utriangle)...,
             label=L"$\textrm{Im}\left\langle \left(\tilde\phi - \phi\right)^2 \right\rangle $")

    fig = plot(fig_Re, fig_Im, layout = (2, 1),size=(750, 400),right_margin=25mm)
    if save
        savefig(fig, string("Imgs/Correcting_vs_difference_HO_Avg2.pdf"))
    end
    display(fig)
end

function plot_avg2_with_corrections(true_solution=0)
    fig = plot(xlabel=L"$\xi$";plot_setup...)

    plot!(fig,2*tp[1:end-1],SS_Re .+ SS2_Re;
                    markers_dict(1,:square)...,label=L"$\textrm{Re}\langle \left(\tilde\phi - \phi\right)^2 \rangle $")
    plot!(fig,2*tp[1:end-1], SS2_Re;
                    markers_dict(1,:square)...,alpha=0.3,label=L"$\textrm{Re}\left\langle \phi (iS) \right\rangle$")
    plot!(fig,2*tp[1:end-1],(avg2_Re .- dfScroSolMinkHO["avg2Re"][1] );
                    markers_dict(2,:utriangle)...,label=L"$\textrm{Re}\langle \Delta \phi^2 \rangle $")

    plot!(fig,2*tp[1:end-1], SS_Im .+ SS2_Im;
                    markers_dict(3,:star)...,label=L"$\textrm{Im}\langle \left(\tilde\phi - \phi\right)^2 \rangle $")
    plot!(fig,2*tp[1:end-1],SS2_Im;
                    markers_dict(3,:star)...,alpha=0.3,label=L"$\textrm{Im}\left\langle \phi (iS) \right\rangle$")
    plot!(fig,2*tp[1:end-1],avg2_Im;
                    markers_dict(4,:utriangle)...,label=L"$\textrm{Im}\langle \Delta \phi^2 \rangle $")

    if save
        savefig(fig, string("Imgs/Correcting_vs_difference_AHO_Avg2.pdf"))
    end

    display(fig)
end


text:
function insert_lines!(js_scene, points)
    evaljs(js_scene, @js begin
        @var scene = window.scene
        @var THREE = window.THREE
        console.log(THREE)
        @var material = @new THREE.LineBasicMaterial(d(
            color = "#001100"
        ));
        @var geometry = @new THREE.Geometry();
        @var points = $(points[]);
        for i=0:points.length-1
            # TODO can we directly serialize this as a Vector3?
            geometry.vertices.push(@new(THREE.Vector3(points[i][0], points[i][1], points[i][2])))
        end
        scene.add(@new THREE.Line(geometry, material))
    end)
end

w = Window()
js_scene = ThreeScene();
tools(w)
points = tuple.(rand(10).* 100, rand(10).* 100, rand(10).* 100)
body!(w, js_scene(dom"div#container"()))
points = Observable(js_scene, "points", tuple.(rand(10).* 10, rand(10).* 10, rand(10).* 10))
insert_lines!(js_scene, points)


text:
module MatrixProductStates

export mpsgroundstate, mpsproductstate, mpsrandom
export makempo, makemps, addsiteMPS, steady_state_search!, steady_state_search_two_site!
export applyMPOtoMPO, apply_site_MPOtoMPO, applyMPOtoMPS, build_env
export mps_to_mpo, mpo_to_mps, mps_to_mpo_dens, mpo_to_mps_dens, mpo_to_mps_site
export leftorthonormalizationQR!, leftorthogonalizeQR!, compressionSVD!
export compress_var, compress_var_apply_H, compress_sum_var, RK4stp_apply_H_half, RK4stp_apply_H_half!
export compress_sum_var_apply_H, RK4stp_apply_H, ground_state_search!, compress_sum_var_apply_H!
export Arnoldi_expm!, normalize_lo, compress_var_apply_H!, compress_mpo!, compress_var!
export compress_var_two_site, compress_var_apply_H_two_site, compress_sum_var_two_site
export compress_sum_var_apply_H_two_site, RK4stp_apply_H_half_two_site, RK4stp_apply_H_half_full_two_site
export ground_state_search_full, expandMPS, mpsdims, normMPS, local_exp
export normalize_lo!, norm_lo, conj_site_mpo, conj_mpo, scal_prod_no_conj, scal_op_prod
export apply_site_operator!, update_renv, prod_LR

using TensorOperations
using LinearMaps
#using CuArrays

import Base.randn
#import CuArrays.CUSOLVER.qrq!

Base.@irrational SQRT_HALF 0.7071067811865475244008  sqrt(big(0.5))

randn(rng::AbstractRNG, ::Type{Complex{T}}) where {T<:AbstractFloat} =
    Complex{T}(SQRT_HALF*randn(rng, T), SQRT_HALF*randn(rng, T))

include("mps_base.jl")
include("constructors.jl")
include("compression.jl")
include("steadystate.jl")
include("groundstate.jl")
include("timeevolution.jl")
include("orthonormalisation.jl")

"""
return a vec of the bond dimensions
"""

function mpsdims(mps)

    na = length(mps)
    dims = ones(Int64, na + 1)
    for jj = 1:na
        dims[jj] = size(mps[jj], 1)
    end

    return dims

end

function diffMP(mps1,mps2)

    n = length(mps2);
    diffm = 0
    for i = 1:n
        diffm += sum(abs2.(mps1[i][:] - mps2[i][:]))
    end
    return diffm
end

"""
make operator product mpo2*mpo1
"""

function applyMPOtoMPO(mpo1,mpo2)

        n = length(mpo2)
        mpoout = similar(mpo1)

        for i = 1:n
            mpoout[i] = apply_site_MPOtoMPO(mpo1[i],mpo2[i])
        end
        return mpoout
end

function applyMPOtoMPO_to(mpo1,mpo2)

        n = length(mpo2)
        mpoout = similar(mpo1)

        for i = 1:n
            mpoout[i] = apply_site_MPOtoMPO_to(mpo1[i],mpo2[i])
        end
        return mpoout
end

function apply_site_MPOtoMPO(mpo1,mpo2)

        (dwl1,d,dwr1,_) = size(mpo1)
        (dwl2,_,dwr2,_) = size(mpo2)
        ai = reshape(mpo2, dwl2*d*dwr2, d)*reshape(permutedims(mpo1, [2, 1, 3, 4]), d, dwl1*dwr1*d)
        ai = reshape(ai, dwl2, d, dwr2, dwl1, dwr1, d)
        ai = permutedims(ai, [4, 1, 2, 5, 3, 6])
        reshape(ai,(dwl1*dwl2,d,dwr1*dwr2,d))

end

function apply_site_MPOtoMPO_to(mpo1,mpo2)

        (dwl1,d,dwr1,_) = size(mpo1)
        (dwl2,_,dwr2,_) = size(mpo2)
        @tensor ai[-1,-2,-3,-5,-4,-6] := mpo1[-1,1,-5,-6]*mpo2[-2,-3,-4,1]
        mpoout = reshape(ai,(dwl1*dwl2,d,dwr1*dwr2,d))

end

"""
applyMPOtoMPS(mpsin,mpo);

Multiplies the MPS A and the MPO U
"""

function applyMPOtoMPS(mpsin,mpo)

    n = length(mpsin)
    mpsout = similar(mpsin)

    for i = 1:n
        mpsout[i] = apply_site_MPOtoMPS(mpsin[i], mpo[i])
    end
    return mpsout
end

function apply_site_MPOtoMPS(ms,mo)

    (dl,_,dr) = size(ms)
    (dwl,d,dwr,_) = size(mo)
    @tensor ai[-1,-2,-3,-5,-4] := ms[-1,1,-5]*mo[-2,-3,-4,1]
    mpsout = reshape(ai,(dl*dwl,d,dr*dwr))

end

function apply_site_operator!(mps_site, op)

    m1, m2, m3 = size(mps_site)
    temp_site = reshape(permutedims(mps_site, (2, 1, 3)), m2, m1*m3)
    temp_site = op*temp_site
    temp_site = permutedims(reshape(temp_site, m2, m1, m3), (2, 1, 3))
    copy!(mps_site, temp_site)

end

function apply_site_operator_to!(mps_site, op)

    @tensor mps_site[-1, -2, -3] = mps_site[-1, 1, -3]*op[-2, 1]

end

"""
product of state and local operator with left and right environments operators
"""

function O_prod_LR(L,O,R)

    @tensor A[-1,-2,-3,-4,-5,-6] := L[-4,1,-1]*O[1,-2,2,-5]*R[-6,2,-3]

end

"""
product of two sites states and local operators with left and right environments operators
"""

function A_prod_LOR_two_site(L,A1,O1,A2,O2,R)

    @tensor A[-1,-2,-3,-4] := L[1,2,-1]*A1[1,3,4]*O1[2,-2,5,3]*A2[4,6,7]*O2[5,-3,8,6]*R[7,8,-4]

end

function prod_AOR_direct(A,O,R)

    a1,a2,a3 = size(A)
    _,r2,r3 = size(R)
    o1 = size(O, 1)
    AOR = complex(zeros(a1, o1, a2, r3))

    @inbounds for r3i in 1:r3, o2i in 1:a2, o1i in 1:o1, a1i in 1:a1, a3i in 1:a3, a2i in 1:a2, r2i in 1:r2
        @fastmath AOR[a1i, o1i, o2i, r3i] += A[a1i, a2i, a3i]*R[a3i, r2i, r3i]*
                                                O[o1i, o2i, r2i, a2i]
    end

    return AOR

end

"""
update right environment for operator
"""

function update_re_large(A0,O,A,FR)

    @tensor FR[-1,-2,-3,-4] := FR[1,2,3,4]*conj(A0[-4,3,4])*A[-1,5,1]*O[-2,-3,2,5]

end

function update_le_large(A0,O,A,FL)

    @tensor FL[-1,-2,-3,-4] := FL[1,2,3,4]*conj(A0[4,3,-4])*A[1,5,-1]*O[2,-3,-2,5]

end

function A_prod_LOR(L,A,O,R)

	@tensor A[-1,-2,-3] := L[1,2,-1]*A[1,3,4]*O[2,-2,5,3]*R[4,5,-3]

end

function A_prod_LOR_test(L,A,O,R)

    al,d,ar = size(A)
    rl,dwr,rr = size(R)
    ll,dwl,lr = size(L)

    A = reshape(A,al*d,ar)*reshape(R,ar,dwr*rr)
    A = permutedims(reshape(A,al,d,dwr,rr),[3,2,1,4])
    A = reshape(O,dwl*d,d*dwr)*reshape(A,d*dwr,rr*al)
    A = permutedims(reshape(A,dwl,d,al,rr),[3,1,2,4])
    A = reshape(L,ll*dwl,lr).'*reshape(A,ll*dwl,rr*d)
    A = reshape(A,lr,d,rr)

end

function A_prod_LOR_test_2(L,A,O,R)

    al,d,ar = size(A)
    rl,dwr,rr = size(R)
    ll,dwl,lr = size(L)

    A = reshape(A,al*d,ar)*reshape(R,ar,dwr*rr)
    A = reshape(L,ll,dwl*lr).'*reshape(A,ll,d*dwr*rr)
    A = permutedims(reshape(A,dwl,lr,d*dwr,rr),[1,3,2,4])
    O = reshape(permutedims(O,[2,1,4,3]),d,dwr*d*dwl)
    A = O*reshape(A,dwr*d*dwl,lr*rr)
    A = permutedims(reshape(A,d,lr,rr),[2,1,3])

end

"""
product of state and local operator with left and right environments operators
designed for use with eigs
"""

function A_prod_LOR_eigs(L,Ain,O,R)

	A = reshape(Ain[:],size(L,1),size(O,4),size(R,1))
	@tensor A[-1,-2,-3] := L[1,2,-1]*A[1,3,4]*O[2,-2,5,3]*R[4,5,-3]
	A[:]

end

"""
product of state and two local operators with left and right environments operators
designed for use with eigs
"""

function A_prod_L_OO_R_eigs(L,Ain,O1,O2,R)

    A = reshape(Ain[:],size(L,1),size(O1,4),size(R,1))
    @tensor A[-1,-2,-3] := L[1,2,3,-1]*A[1,4,6]*O1[2,5,7,4]*O2[3,-2,8,5]*R[6,7,8,-3]
    A[:]

end

function neighbour_contract(A1,A2)

    @tensor A[-1,-2,-3,-4] := A1[-1,-2,1]*A2[1,-3,-4]
    A[:]

end

function A_prod_L_OO_R_two_site(L,Ain,L1,cL1,L2,cL2,R)

    A = reshape(Ain[:],size(L,1),size(L1,4),size(L2,4),size(R,1))
    @tensor A[-1,-2,-3,-4] := L[1,2,3,-1]*A[1,4,5,10]*L1[2,7,6,4]*L2[6,8,11,5]*
                                    cL1[3,-2,9,7]*cL2[9,-3,12,8]*R[10,11,12,-4]
    A[:]

end

"""
product of two local A tensors with left and right environments
"""

function A_prod_LR_two_site(L,A1,A2,R)

    @tensor A[-1,-2,-3,-4] := L[1,-1]*A1[1,-2,2]*A2[2,-3,3]*R[3,-4]

end

"""
update left environment for two operators
"""

function update_lenv_two_op(A0,O1,O2,A,FL)

    @tensor FL[-1,-2,-3,-4] := FL[1,2,3,4]*A[1,5,-1]*O1[2,6,-2,5]*O2[3,7,-3,6]*conj(A0[4,7,-4])

end

"""
update right environment for two operators
"""

function update_renv_two_op(A0,O1,O2,A,FR)

    @tensor FR[-1,-2,-3,-4] := FR[1,2,3,4]*A[-1,5,1]*O1[-2,6,2,5]*O2[-3,7,3,6]*conj(A0[-4,7,4])

end

"""
Compute norm of arbitrary MPS A, no left or right orthogonality assumed
"""
function normMPS(mps)

    abs(scal_prod(mps, mps))

end

"""
Compute sandwcich of operator by mps, <mps1|O|mps2>
"""
function scal_op_prod(mps1, op, mps2)

    n = length(mps1)
    lenv = similar(mps1[1], 1, 1, 1)
    lenv[1] = 1.0

    for i = 1:n
        #Update lenv from A[i], conj(A[i]) and the previous lenv
        lenv = update_le(mps1[i], op[i], mps2[i], lenv)
    end

    return lenv[1]

end

"""
Compute scalar product of two mps, <mps1|mps2>
"""
function scal_prod(mps1, mps2)

    na = length(mps1)
    lenv = similar(mps1[1], 1, 1)
    lenv[1] = 1.0

    for n = 1:na
        #Update lenv from A[i], conj(A[i]) and the previous lenv
        lenv = update_lenv(mps1[n], mps2[n], lenv)
    end

    return lenv[1]

end

"""
scal_prod(mps1,mps2)

Compute scalar product of two mps, <mps1|mps2> without conjugating mps1
"""

function scal_prod_no_conj(mps1, mps2)

    na = length(mps1)
    lenv = similar(mps1[1], 1, 1)
    lenv[1] = 1.0

    for n = 1:na
        #Update lenv from A[i], conj(A[i]) and the previous lenv
        lenv = update_lenv_nc(mps1[n], mps2[n], lenv)
    end

    return lenv[1]

end

"""
diff_norm(mps1,mps2)

Compute norm of difference between two mps,<mps1|mps1> + <mps2|mps2> - <mps1|mps2> - <mps2|mps1>
"""

function diff_norm(mps1,mps2)

    normMPS(mps1) + normMPS(mps2) - 2*real(scal_prod(mps1,mps2))

end

function conj_mpo(mpo)

	na = length(mpo)
	mpoconj = similar(mpo)

	for jj = 1:na
    	mpoconj[jj] = conj_site_mpo(mpo[jj])
	end

	return mpoconj

end

function conj_site_mpo(mpo_site)

    mpoconj = permutedims(conj(mpo_site),[1,4,3,2])

end

end


text:
# Alias sampling
struct FastSampler{Tv,Ti}
    F::Array{Tv,1}    # array of p[s.A[i]] / (1 / n)
    A::Array{Ti,1}    # array of indices of elements < (1 / n)
    V::Array{Ti,1}    # arary of indices of elements > (1 / n)
    n::Ti
    rng::AbstractRNG  # the pseudo random number generator
end

# sample a random number
function sample(s::FastSampler{Tv,Ti}) where {Tv,Ti}
    #i = rand(1:s.n)
    i = ceil(Ti,rand(s.rng)*s.n)
    f = rand(s.rng)
    if f < s.F[i]
        return s.A[i]
    else
        return s.V[i]
    end
end

# sample a random number
function sampleMany(s::FastSampler{Tv,Ti},sampCount::Ti) where {Tv,Ti}
    #i = rand(1:s.n)
    samples = Array{Ti,1}(sampCount)
    for j = 1:sampCount
        i = ceil(Ti,rand(s.rng)*s.n)
        f = rand(s.rng)
        if f < s.F[i]
            samples[j] = s.A[i]
        else
            samples[j] = s.V[i]
        end
    end
    return samples
end

"""
    s = FastSampler(p)
"""
function FastSampler(p::Array{Tv,1}; residual::Bool = false, rng::AbstractRNG=Random.GLOBAL_RNG) where Tv

    @assert(minimum(p) > 0, "The probability vector has a negative entry")

    n = length(p)

    # generate the Sampler
    F = Array{Tv}(undef, n)
    A = Array{Int64}(undef, n)
    V = Array{Int64}(undef, n)

    # generate a normalized p, so that it sums up to n
    newp = n * copy(p) / sum(p)

    # fill up the arays A and V
    pos = 0                      # number of filled or half filled buckets
    posFilled = 0                # number of filled buckets

    for i in 1:n
        if newp[i] < 1 #&& !epsequal(newp[i], 1)
            pos = pos + 1
            A[pos] = i
            F[pos] = newp[i]
        end
    end

    err = 0
    for i in 1:n
        if newp[i] >= 1
            val = newp[i]

            while val >= 1
                if posFilled < pos
                    # fill up a bucket already containing a small element, and subtract that value from the current element
                    posFilled = posFilled + 1
                    V[posFilled] = i
                    val = val - (1 - F[posFilled])
                else
                    # create a bucket with only one large element. note posFilled = pos
                    pos = pos + 1
                    posFilled = posFilled + 1
                    A[pos] = V[pos] = i
                    F[pos] = 1
                    val = val - 1
                end
            end

            if val > 0 && pos == n
            	err = err + val
            end

            # treat the case in which val becomes smaller than 1. Have to do the pos < n check because of precision errors
            if val > 0 && pos < n
                pos = pos + 1
                A[pos] = i
                F[pos] = val
            end
        end
    end

    while posFilled < n
        # if we are in this case, posFilled should be equal to pos - 1
        posFilled = posFilled + 1
        V[posFilled] = A[posFilled]
        err += 1-F[posFilled]
    end

    if residual
        return FastSampler(F, A, V, n, rng), err
    else
        return FastSampler(F, A, V, n, rng)
    end
end

"""
    s = blockSample(p; k = length(p))

Compute numbers sampled with probability proportional to the vector p.
They are returned in order.
If need be, they can be permuted with randperm.
"""
function blockSample(p; k = length(p))
    samp = zeros(Int,k)

    r = sort(rand(k))

    cs = cumsum(p)
    cs = cs / cs[end]

    j = 1
    for i in 1:k
        while r[i] > cs[j]
            j += 1
        end

        samp[i] = j
    end

    return samp
end


text:
using OrdinaryDiffEq
using Trixi

# Define new structs inside a module to allow re-evaluating the file.
# This module name needs to be unique among all examples, otherwise Julia will throw warnings
# if multiple test cases using the same module name are run in the same session.
module TrixiExtensionEulerRotated

using Trixi

# initial_condition_convergence_test transformed to the rotated rectangle
struct InitialConditionSourceTermsRotated
  sin_alpha::Float64
  cos_alpha::Float64
end

function InitialConditionSourceTermsRotated(alpha)
  sin_alpha, cos_alpha = sincos(alpha)

  InitialConditionSourceTermsRotated(sin_alpha, cos_alpha)
end

function (initial_condition::InitialConditionSourceTermsRotated)(x, t, equations::CompressibleEulerEquations2D)
  sin_ = initial_condition.sin_alpha
  cos_ = initial_condition.cos_alpha

  # Rotate back to unit square and translate from [-1, 1]^2 to [0, 2]^2

  # Clockwise rotation by α and translation by 1
  # Multiply with [  cos(α)  sin(α);
  #                 -sin(α)  cos(α)]
  x1 =  cos_ * x[1] + sin_ * x[2] + 1
  x2 = -sin_ * x[1] + cos_ * x[2] + 1

  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  ini = c + A * sin(ω * (x1 + x2 - t))

  rho = ini
  rho_v1 = ini
  rho_v2 = ini
  rho_e = ini^2

  # Rotate velocity vector counterclockwise
  # Multiply with [ cos(α)  -sin(α);
  #                 sin(α)   cos(α)]
  rho_v1_rot = cos_ * rho_v1 - sin_ * rho_v2
  rho_v2_rot = sin_ * rho_v1 + cos_ * rho_v2

  return SVector(rho, rho_v1_rot, rho_v2_rot, rho_e)
end

@inline function (source_terms::InitialConditionSourceTermsRotated)(u, x, t, equations::CompressibleEulerEquations2D)
  sin_ = source_terms.sin_alpha
  cos_ = source_terms.cos_alpha

  # Rotate back to unit square and translate from [-1, 1]^2 to [0, 2]^2

  # Clockwise rotation by α and translation by 1
  # Multiply with [  cos(α)  sin(α);
  #                 -sin(α)  cos(α)]
  x1 =  cos_ * x[1] + sin_ * x[2] + 1
  x2 = -sin_ * x[1] + cos_ * x[2] + 1

  # Same settings as in `initial_condition`
  c = 2
  A = 0.1
  L = 2
  f = 1/L
  ω = 2 * pi * f
  γ = equations.gamma

  si, co = sincos((x1 + x2 - t)*ω)
  tmp1 = co * A * ω
  tmp2 = si * A
  tmp3 = γ - 1
  tmp4 = (2*c - 1)*tmp3
  tmp5 = (2*tmp2*γ - 2*tmp2 + tmp4 + 1)*tmp1
  tmp6 = tmp2 + c

  du1 = tmp1
  du2 = tmp5
  du3 = tmp5
  du4 = 2*((tmp6 - 1)*tmp3 + tmp6*γ)*tmp1

  # Original terms (without performanc enhancements)
  # du1 = cos((x1 + x2 - t)*ω)*A*ω
  # du2 = (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
  #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
  # du3 = (2*sin((x1 + x2 - t)*ω)*A*γ - 2*sin((x1 + x2 - t)*ω)*A +
  #                             2*c*γ - 2*c - γ + 2)*cos((x1 + x2 - t)*ω)*A*ω
  # du3 = 2*((c - 1 + sin((x1 + x2 - t)*ω)*A)*(γ - 1) +
  #                             (sin((x1 + x2 - t)*ω)*A + c)*γ)*cos((x1 + x2 - t)*ω)*A*ω

  # Rotate velocity vector counterclockwise
  # Multiply with [ cos(α)  -sin(α);
  #                 sin(α)   cos(α)]
  du2_rotated = cos_ * du2 - sin_ * du3
  du3_rotated = sin_ * du2 + cos_ * du3

  return SVector(du1, du2_rotated, du3_rotated, du4)
end

end # module TrixiExtensionEulerRotated

import .TrixiExtensionEulerRotated

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

alpha = 0.1
initial_condition_source_terms = TrixiExtensionEulerRotated.InitialConditionSourceTermsRotated(alpha)
sin_ = initial_condition_source_terms.sin_alpha
cos_ = initial_condition_source_terms.cos_alpha
T = [cos_ -sin_; sin_ cos_]

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

mapping(xi, eta) = T * SVector(xi, eta)

cells_per_dimension = (16, 16)

mesh = CurvedMesh(cells_per_dimension, mapping)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_source_terms, solver,
                                    source_terms=initial_condition_source_terms)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
struct ServeFiles <: ServerMode end

function handle(server, ::ServeFiles, event::DocUpdated)
    addsource!(server, event.name, event.doc)
    addrewrite!(server, event.name)
    addbuild!(server, event.name)
end

function geteventsource(::ServeFiles, server, ch)
    return FileServer(server.builder.dir)
end

function initialize(::ServeFiles, server)
    @info "Starting initial build..."
    fullbuild(server.project, server.builder)
    @info "Done."
end

struct ServeFilesLazy <: ServerMode end

function initialize(::ServeFilesLazy, server)
    build(server.builder, server.project)
end

function handle(server, ::ServeFilesLazy, event::DocUpdated)
    addsource!(server, event.name, event.doc)
    if event.name in keys(server.project.outputs)
        addrewrite!(server, event.name)
        addbuild!(server, event.name)
    end
end

function handle(server, ::ServeFilesLazy, event::DocRequested)
    if (event.name ∈ keys(server.project.sources)) && !(event.name ∈ keys(server.project.outputs))
        @info "Building $(event.name) for the first time..."
        addrewrite!(server, event.name)
        addbuild!(server, event.name)
    end
end

function geteventsource(::ServeFilesLazy, server, ch)
    builddir = server.builder.dir
    return FileServer(builddir, preprocess_request = req -> _lazyservecallback(req, ch, builddir))
end

function _lazyservecallback(req, ch, builddir)
    if !endswith(req.target, ".html")
        return req
    else
        buildpath = joinpath(builddir, req.target[2:end])
        try mkpath(parent(buildpath)) catch end
        touch(buildpath)
        #=
        open(joinpath(builddir, req.target[2:end]), "w") do f
            write(f, "Building...")
        end
        =#
        sourcepath = Path(req.target[2:end-5])  # cut off / and .html
        put!(ch, DocRequested(sourcepath))
        # Give time to build so file server doesn't instantly return a 404
        sleep(0.05)
        return req
    end
end


text:
export gemm_batch

doc"""
    gemm_batch(tA::Char, tB::Char, alpha, A::Var, B::Var)
"""
function gemm_batch(tA::Char, tB::Char, alpha, A::Var, dimsA, B::Var, dimsB)
    T = eltype(A)
    ydata = gemm_batch(tA, tB, T(alpha), A.data, B.data)
    Var(ydata, ∇gemm_batch!, (tA,tB,alpha,A,B))
end
gemm_batch(tA::Char, tb::Char, alpha, A::Node, B) = Node(gemm_batch, (tA,tB,alpha,A,B))
gemm_batch(tA::Char, tb::Char, alpha, A, B::Node) = Node(gemm_batch, (tA,tB,alpha,A,B))
gemm_batch(A, B; tA='N', tB='N', alpha=1) = gemm_batch(tA, tB, alpha, A, B)

function gemm_batch(tA::Char, tB::Char, alpha, A::Array{T,3}, B::Array{T,3}) where T
    @assert size(A,3) == size(B,3)
    m = size(A, tA == 'N' ? 1 : 2)
    n = size(B, tB == 'N' ? 2 : 1)
    C = Array{T}(m, n, size(A,3))
    for i = 1:size(A,3)
        gemm!(tA, tB, T(alpha), view(A,:,:,i), view(B,:,:,i), T(0), view(C,:,:,i))
    end
    C
end

function gemm_batch(tA::Char, tB::Char, alpha, A::CuArray{T,3}, B::CuArray{T,3}) where T
    CUBLAS.gemm_batched(tA, tB, T(alpha), A, B)
end

function ∇gemm_batch!(C::Var, tA::Char, tB::Char, alpha, A::Var, B::Var)
    isnothing(A.grad) || ∇gemm_batch_A!(C.grad, tA, tB, alpha, A.grad, B.data)
    isnothing(B.grad) || ∇gemm_batch_B!(C.grad, tA, tB, alpha, A.data, B.grad)
end

function ∇gemm_batch_A!(gC::Array{T,3}, tA::Char, tB::Char, alpha, gA::Array{T,3}, B::Array{T,3}) where T
    for i = 1:size(gC,3)
        ∇gemm_A!(view(gC,:,:,i), tA, tB, alpha, view(gA,:,:,i), view(B,:,:,i))
    end
end
function ∇gemm_batch_B!(gC::Array{T,3}, tA::Char, tB::Char, alpha, A::Array{T,3}, gB::Array{T,3}) where T
    for i = 1:size(gC,3)
        ∇gemm_B!(view(gC,:,:,i), tA, tB, alpha, view(A,:,:,i), view(gB,:,:,i))
    end
end

function ∇gemm_batch_A!(gC::CuArray{T,3}, tA::Char, tB::Char, alpha, gA::CuArray{T,3}, B::CuArray{T,3}) where T
    if tA == 'N'
        CUBLAS.gemm_batched!('N', tB=='N' ? 'T' : 'N', T(alpha), gC, B, T(1), gA)
    else
        CUBLAS.gemm_batched!(tB, 'T', T(alpha), B, gC, T(1), gA)
    end
end
function ∇gemm_batch_B!(gC::CuArray{T,3}, tA::Char, tB::Char, alpha, A::CuArray{T,3}, gB::CuArray{T,3}) where T
    if tB == 'N'
        CUBLAS.gemm_batched!(tA=='N' ? 'T' : 'N', 'N', T(alpha), A, gC, T(1), gB)
    else
        CUBLAS.gemm_batched!('T', tA, T(alpha), gC, A, T(1), gB)
    end
end


text:
int_rules_4_1_0_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.0.2*(a*trg)^m*(b*tan)^n =#
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~b) * (~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~m) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && EqQ((~m + ~n) - 1, 0)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * tan(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => (-1 / ~f) * Subst(Antiderivative((1 - (~x) ^ 2) ^ (((~m + ~n) - 1) / 2) / (~x) ^ ~n, ~x), ~x, cos(~e + ~f * ~x)) <-- FreeQ([~e, ~f], ~x) && IntegersQ(~m, ~n, ((~m + ~n) - 1) / 2)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~b * ff) / ~f) * Subst(Antiderivative((ff * ~x) ^ (~m + ~n) / ((~b) ^ 2 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, (~b * tan(~e + ~f * ~x)) / ff)) <-- FreeQ([~b, ~e, ~f, ~n], ~x) && IntegerQ(~m / 2)
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * tan(~(e') + ~(f') * ~x) ^ ~(n'), ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((ff * ~x) ^ (~m + ~n) / ((~a) ^ 2 - ff ^ 2 * (~x) ^ 2) ^ ((~n + 1) / 2), ~x), ~x, (~a * sin(~e + ~f * ~x)) / ff)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && IntegerQ((~n + 1) / 2)
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / ((~a) ^ 2 * ~f * (~n - 1)) - (((~b) ^ 2 * (~m + 2)) / ((~a) ^ 2 * (~n - 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~n, 1) && ((LtQ(~m, -1) || EqQ(~m, -1) && EqQ(~n, 3 / 2)) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~n - 1)) - (((~b) ^ 2 * ((~m + ~n) - 1)) / (~n - 1)) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (GtQ(~n, 1) && (IntegersQ(2 * ~m, 2 * ~n) && Not(GtQ(~m, 1) && Not(IntegerQ((~m - 1) / 2)))))
    @apply_utils Antiderivative(sqrt(~(a') * sin(~(e') + ~(f') * ~x)) / (~(b') * tan(~(e') + ~(f') * ~x)) ^ (3 / 2), ~x) => (2 * sqrt(~a * sin(~e + ~f * ~x))) / (~b * ~f * sqrt(~b * tan(~e + ~f * ~x))) + ((~a) ^ 2 / (~b) ^ 2) * Antiderivative(sqrt(~b * tan(~e + ~f * ~x)) / (~a * sin(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x)
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * ~m) - (((~a) ^ 2 * (~n + 1)) / ((~b) ^ 2 * ~m)) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (LtQ(~n, -1) && (GtQ(~m, 1) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~m + ~n + 1)) - ((~n + 1) / ((~b) ^ 2 * (~m + ~n + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (LtQ(~n, -1) && (NeQ(~m + ~n + 1, 0) && (IntegersQ(2 * ~m, 2 * ~n) && Not(EqQ(~n, -3 / 2) && EqQ(~m, 1)))))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-(~b) * (~a * sin(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~m) + (((~a) ^ 2 * ((~m + ~n) - 1)) / ~m) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && ((GtQ(~m, 1) || EqQ(~m, 1) && EqQ(~n, 1 / 2)) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~b * (~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / ((~a) ^ 2 * ~f * (~m + ~n + 1)) + ((~m + 2) / ((~a) ^ 2 * (~m + ~n + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (LtQ(~m, -1) && (NeQ(~m + ~n + 1, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~m * tan(~(e') + ~(f') * ~x) ^ ~n, ~x) => (1 / (~a) ^ ~n) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + ~n) / cos(~e + ~f * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ(~n) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((cos(~e + ~f * ~x) ^ ~n * (~b * tan(~e + ~f * ~x)) ^ ~n) / (~a * sin(~e + ~f * ~x)) ^ ~n) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + ~n) / cos(~e + ~f * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~n)) && (ILtQ(~m, 0) || (EqQ(~m, 1) && EqQ(~n, -1 / 2) || IntegersQ(~m - 1 / 2, ~n - 1 / 2))))
    @apply_utils Antiderivative((~(a') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * cos(~e + ~f * ~x) ^ (~n + 1) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * (~a * sin(~e + ~f * ~x)) ^ (~n + 1))) * Antiderivative((~a * sin(~e + ~f * ~x)) ^ (~m + ~n) / cos(~e + ~f * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~(a') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * cos(~e + ~f * ~x)) ^ FracPart(~m) * (sec(~e + ~f * ~x) / ~a) ^ FracPart(~m) * Antiderivative((~b * tan(~e + ~f * ~x)) ^ ~n / (sec(~e + ~f * ~x) / ~a) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~(a') * cot(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * cot(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ ~m * Antiderivative((~b * tan(~e + ~f * ~x)) ^ (~n - ~m), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-((~a * sec(~e + ~f * ~x)) ^ ~m) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * ~m) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && EqQ(~m + ~n + 1, 0)
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~a / ~f) * Subst(Antiderivative((~a * ~x) ^ (~m - 1) * (-1 + (~x) ^ 2) ^ ((~n - 1) / 2), ~x), ~x, sec(~e + ~f * ~x)) <-- FreeQ([~a, ~e, ~f, ~m], ~x) && (IntegerQ((~n - 1) / 2) && Not(IntegerQ(~m / 2) && LtQ(0, ~m, ~n + 1)))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (1 / ~f) * Subst(Antiderivative((~b * ~x) ^ ~n * (1 + (~x) ^ 2) ^ (~m / 2 - 1), ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~b, ~e, ~f, ~n], ~x) && (IntegerQ(~m / 2) && Not(IntegerQ((~n - 1) / 2) && LtQ(0, ~n, ~m - 1)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a) ^ 2 * (~a * sec(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) - (((~a) ^ 2 * (~m - 2)) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (LtQ(~n, -1) && ((GtQ(~m, 1) || EqQ(~m, 1) && EqQ(~n, -3 / 2)) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * (~n + 1)) - ((~m + ~n + 1) / ((~b) ^ 2 * (~n + 1))) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (LtQ(~n, -1) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ~m) - (((~b) ^ 2 * (~n - 1)) / ((~a) ^ 2 * ~m)) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (GtQ(~n, 1) && ((LtQ(~m, -1) || EqQ(~m, -1) && EqQ(~n, 3 / 2)) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * (~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 1)) / (~f * ((~m + ~n) - 1)) - (((~b) ^ 2 * (~n - 1)) / ((~m + ~n) - 1)) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (GtQ(~n, 1) && (NeQ((~m + ~n) - 1, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-((~a * sec(~e + ~f * ~x)) ^ ~m) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * ~m) + ((~m + ~n + 1) / ((~a) ^ 2 * ~m)) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ (~m + 2) * (~b * tan(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && ((LtQ(~m, -1) || EqQ(~m, -1) && EqQ(~n, -1 / 2)) && IntegersQ(2 * ~m, 2 * ~n))
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a) ^ 2 * (~a * sec(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~b * ~f * ((~m + ~n) - 1)) + (((~a) ^ 2 * (~m - 2)) / ((~m + ~n) - 1)) * Antiderivative((~a * sec(~e + ~f * ~x)) ^ (~m - 2) * (~b * tan(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && ((GtQ(~m, 1) || EqQ(~m, 1) && EqQ(~n, 1 / 2)) && (NeQ((~m + ~n) - 1, 0) && IntegersQ(2 * ~m, 2 * ~n)))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) / sqrt(~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (sqrt(sin(~e + ~f * ~x)) / (sqrt(cos(~e + ~f * ~x)) * sqrt(~b * tan(~e + ~f * ~x)))) * Antiderivative(1 / (sqrt(cos(~e + ~f * ~x)) * sqrt(sin(~e + ~f * ~x))), ~x) <-- FreeQ([~b, ~e, ~f], ~x)
    @apply_utils Antiderivative(sqrt(~(b') * tan(~(e') + ~(f') * ~x)) / sec(~(e') + ~(f') * ~x), ~x) => ((sqrt(cos(~e + ~f * ~x)) * sqrt(~b * tan(~e + ~f * ~x))) / sqrt(sin(~e + ~f * ~x))) * Antiderivative(sqrt(cos(~e + ~f * ~x)) * sqrt(sin(~e + ~f * ~x)), ~x) <-- FreeQ([~b, ~e, ~f], ~x)
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~a) ^ (~m + ~n) * (~b * tan(~e + ~f * ~x)) ^ ~n) / ((~a * sec(~e + ~f * ~x)) ^ ~n * (~b * sin(~e + ~f * ~x)) ^ ~n)) * Antiderivative((~b * sin(~e + ~f * ~x)) ^ ~n / cos(~e + ~f * ~x) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (IntegerQ(~n + 1 / 2) && IntegerQ(~m + 1 / 2))
    #= Antiderivative(((~a')*sec((~e')+(~f')*(~x)))^(~m')*((~b')*tan((~e')+(~f')*(~x)))^(~n),~x):=   (a*sec(e+f*x))^m*(b*tan(e+f*x))^(n+1)*(cos(e+f*x)^2)^((m+n+1)/2)/(b* f*(b*sin(e+f*x))^(n+1))* Subst(Antiderivative(x^n/(1-x^2/b^2)^((m+n+1)/2),x),x,b*sin(e+f*x)) <--  FreeQ([a,b,e,f,m,n],x) && Not(IntegerQ((n-1)/2)) &&  Not(IntegerQ(m/2)) =#
    @apply_utils Antiderivative((~(a') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~a * sec(~e + ~f * ~x)) ^ ~m * (~b * tan(~e + ~f * ~x)) ^ (~n + 1) * (cos(~e + ~f * ~x) ^ 2) ^ ((~m + ~n + 1) / 2)) / (~b * ~f * (~n + 1))) * HypergeometricFunctions._₂F₁((~n + 1) / 2, (~m + ~n + 1) / 2, (~n + 3) / 2, sin(~e + ~f * ~x) ^ 2) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ((~n - 1) / 2)) && Not(IntegerQ(~m / 2)))
    @apply_utils Antiderivative((~(a') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * tan(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~a * csc(~e + ~f * ~x)) ^ FracPart(~m) * (sin(~e + ~f * ~x) / ~a) ^ FracPart(~m) * Antiderivative((~b * tan(~e + ~f * ~x)) ^ ~n / (sin(~e + ~f * ~x) / ~a) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.0.2 (a trg)^m (b tan)^n.jl


text:
#: ============================================

@inline lookup_code(P::Vector{TPattern},c::Int) = any(x->any(y->y==c,x),P)
@inline lookup_code(P::Vector{TPattern},cs::Vector{Int}) = any(x->any(c->any(y->y==c,x),cs),P)
@inline lookup_code(P::Vector{TPattern},cs::NTuple{N,Int}) where N = any(x->any(c->any(y->y==c,x),cs),P)

lookup_codes(P::Vector{TPattern},cs::Vector) = all(lookup_code(P,c) for c in cs)

#: ============================================

function get_DATA_by_codes(
    b::Block{RT},
    codes::Vector
    )::Vector{DATATYPE} where RT
    if is_multi(b)
        @assert length(b.R[2])>1  # R[2] = [[1,0,2],[23,24,25,1],...]
        return vcat([get_DATA_by_codes(c,codes)
                        for c in children(b)
                            if lookup_codes(c.R[2],codes) ]...)
    elseif is_single(b) && lookup_codes(b.R[2],codes)
        return DATATYPE[b.DATA]
    else
        return DATATYPE[]
    end
end

get_DATA_by_codes(b::Block{RT}, code::Int) where RT = get_DATA_by_codes(b, [code])

#: ============================================

function get_n_blocks_by_codes(b::Block, codes::Vector; n=1)
    status = 0
    function f(x)
        if is_single(x)
            @assert length(x.R[2])==1
            if lookup_codes(x.R[2],codes)
                status = 1
            end
        end
        return
    end
    return DFS( b,
                f,
                x->nothing,
                Vb->( status==0
                        ? vcat(Vb[1]...)    # hit; reset status and return
                        : (status>=1+n  ? (status=0 ; vcat(Vb[1]...))    # enough recorded; reset status
                                        : (status+=1; vcat(Vb[1]..., Any[(status-1,Vb[2])]))) )  # not yet
            )
                #Vb->( status>=1+n
                #        ? (status=0 ; vcat(Vb[1]...))    # enough recorded; reset status
                #        : (status>0 ? (status+=1; vcat(Vb[1]..., Any[(status-1,Vb[2])]))
                #                    : vcat(Vb[1]...)) ) ) # not yet
end

get_n_blocks_by_code(b::Block, code::Int; n=1) =
    get_n_blocks_by_codes(b, [code]; n=n)

#: ============================================

function get_blocks_max_by_codes(
    b::Block,
    codes::Vector,
    stop_codes::Vector;
    n=800
    )
    status = 0
    function f(x)
        if is_single(x)
            @assert length(x.R[2])==1
            if lookup_codes(x.R[2],codes)
                status = 1
            elseif lookup_codes(x.R[2],stop_codes)
                status = 0
            end
        end
        return
    end
    return DFS( b,
                f,
                x->nothing,
                Vb->(status==0 ? vcat(Vb[1]...)
                               : (status+=1; vcat(Vb[1]..., Any[(status-1,Vb[2])])))  # not yet
            )
end

get_blocks_max_by_code(b::Block, code::Int, stop_code::Int; n=300) =
    get_blocks_max_by_codes(b, [code], [stop_code]; n=n)

#: ============================================

function next_block_by_codes(b::Block, codes::Vector; delay=1)
    status = 0
    function f(x)
        if is_single(x)
            @assert length(x.R[2])==1
            if lookup_codes(x.R[2],codes)
                status = 1
            end
        end
        return
    end
    return DFS( b,
                f,
                x->nothing,
                Vb->( status==0
                        ? vcat(Vb[1]...)    # hit; reset status and return
                        : (status>=1+delay  ? (status=0;  Any[Vb[2],])
                                            : (status+=1; vcat(Vb[1]...)) ) )  # not yet
            )
                #Vb->( status>=1+delay
                #        ? (status=0; Any[Vb[2],])    # hit; reset status and return
                #        : (status>0 ? (status+=1; vcat(Vb[1]...))
                #                    : vcat(Vb[1]...)) ) ) # not yet
end

next_block_by_code(b::Block, code::Int; delay=1) = next_block_by_codes(b, [code]; delay=delay)

#: ============================================

comp(p1::TPattern,p2::TPattern) = (length(p1)==length(p2) && all(i->(p1[i]==p2[i]), 1:length(p1)))

lookup_patt(P::Vector{TPattern},patt::TPattern) = any(x->comp(patt,x),P)

function select_block_by_patt(
    b::Block{TR},
    patt::TPattern
    )::Vector{Block{TR}} where TR
    if lookup_patt(x.R[2],patt)  # R[2] = [[1,0,2],[23,24,25,1],...]
        return vcat([ select_block_by_patt(c,patt)
                        for c in children(b)
                            if lookup_patt(c.R[2],patt) ]...)
    else
        return Block{TR}[b,]
    end
end

function get_DATA_by_patt(
    b::Block{TR},
    patt::TPattern
    )::Vector{DATATYPE} where TR
    if  is_multi(b)  # R[2] = [[1,0,2],[23,24,25,1],...]
        return vcat([ get_DATA_by_patt(c,patt)
                        for c in children(b)
                            if lookup_patt(c.R[2],patt) ]...)
    elseif is_single(b) && lookup_patt(b.R[2],patt)
        return  DATATYPE[ b.DATA, ]
    else
        return  DATATYPE[  ]
    end
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser