text:
"""
    ContourFiber2D(f; w1_fixed=true)

The fiber for the contour of `f` with the first axis fixed if `w1_fixed == true`.
"""
function ContourFiber2D(p::MP.AbstractPolynomialLike, w_i=0.0; w1_fixed=true)
    @assert (length(MP.variables(p)) == 2) "Expected a bivariate polynomial, but got $(p)."
    # MP.@similarvariable p s
    z1, z2 = MP.variables(p)
    realvars, imagvars = realifyvariables(p)

    fre, fim = realify(p, realvars, imagvars)
    g1re , g1im = realify(z1 * MP.differentiate(p, z1), realvars, imagvars)
    g2re , g2im = realify(z2 * MP.differentiate(p, z2), realvars, imagvars)
    g = g1re * g2im - g1im * g2re

    append!(realvars, imagvars)
    vars = realvars

    f_re = SP.Polynomial(1.0 * fre, vars)
    f_im = SP.Polynomial(1.0 * fim, vars)
    g = SP.Polynomial(1.0 * g, vars)

    v_i = exp(w_i)
    exponents = SP.exponents.((f_re, f_im, g))
    original_coefficients = copy.(SP.coefficients.((f_re, f_im, g)))
    ContourFiber2D(f_re, f_im, g, exponents, v_i, SVector(0.0, 0.0, 0.0, 0.0), zeros(3, 4), w1_fixed, original_coefficients)
end

function condition_contour_fiber!(f::SP.Polynomial, E, w1_fixed::Bool, v_i::Float64)
    if w1_fixed
        v1 = v_i
        m = maximum(map(j -> E[1,j]+E[3,j], 1:size(E, 2)))
        SP.scale_coefficients!(f, clamp(inv(v1^m), 1e-12, 1e12))
    else
        v2 = v_i
        m = maximum(map(j -> E[2,j]+E[4,j], 1:size(E, 2)))
        SP.scale_coefficients!(f, clamp(inv(v2^m), 1e-12, 1e12))
    end
end
function condition!(fiber::ContourFiber2D)
    fiber.f_re.coefficients .= fiber.original_coefficients[1]
    fiber.f_im.coefficients .= fiber.original_coefficients[2]
    fiber.g.coefficients .= fiber.original_coefficients[3]

    condition_contour_fiber!(fiber.f_re, fiber.exponents[1], fiber.w1_fixed, fiber.v_i)
    condition_contour_fiber!(fiber.f_im, fiber.exponents[2], fiber.w1_fixed, fiber.v_i)
    condition_contour_fiber!(fiber.g, fiber.exponents[3], fiber.w1_fixed, fiber.v_i)
end

function update_fiber!(c::ContourFiber2D, w_i::Float64)
    c.v_i = exp(w_i)
    condition!(c)
    c
end

function fix_axis!(c::ContourFiber2D, axis::Symbol)
    if axis == :w1
        c.w1_fixed = true
    elseif axis == :w2
        c.w1_fixed = false
    else
        throw(DomainError(":$axis is not a valid argument. Expected :w1 or :w2"))
    end
end

function phi!(fiber::ContourFiber2D, u)
    if fiber.w1_fixed
        v1 = fiber.v_i
        v2 = exp(u[3])
    else
        v1 = exp(u[3])
        v2 = fiber.v_i
    end
    x = SVector(v1 * cos(u[1]), v2 * cos(u[2]), v1 * sin(u[1]), v2 * sin(u[2]))
    fiber.x = x
    x
end

function evaluate(fiber::ContourFiber2D{T}, u::SVector{3,T}) where T
    x = phi!(fiber, u)

    SVector{3, T}(
        SP.evaluate(fiber.f_re, x),
        SP.evaluate(fiber.f_im, x),
        SP.evaluate(fiber.g, x))
end

function jacobian(fiber::ContourFiber2D, u::AbstractVector, precomputed=false)
    x = precomputed ? fiber.x : phi!(fiber, u)
    U = fiber.U
    @inbounds begin
        ∇f_re = SP.gradient(fiber.f_re, x)
        ∇f_im = SP.gradient(fiber.f_im, x)
        ∇g = SP.gradient(fiber.g, x)
        for i=1:4
            U[1, i] = ∇f_re[i]
            U[2, i] = ∇f_im[i]
            U[3, i] = ∇g[i]
        end
    end

    @inbounds out = begin
        mx3, mx4 = -x[3], -x[4]
        ∂f_re∂θ1 = muladd(U[1, 3], x[1], U[1, 1] * mx3)
        ∂f_re∂θ2 = muladd(U[1, 4], x[2], U[1, 2] * mx4)

        ∂f_im∂θ1 = muladd(U[2, 3], x[1], U[2, 1] * mx3)
        ∂f_im∂θ2 = muladd(U[2, 4], x[2], U[2, 2] * mx4)

        ∂g∂θ1 = muladd(U[3, 3], x[1], U[3, 1] * mx3)
        ∂g∂θ2 = muladd(U[3, 4], x[2], U[3, 2] * mx4)

        if fiber.w1_fixed
            ∂f_re∂w_i = muladd(U[1, 2], x[2], U[1, 4] * x[4])
            ∂f_im∂w_i = muladd(U[2, 2], x[2], U[2, 4] * x[4])
            ∂g∂w_i = muladd(U[3, 2], x[2], U[3, 4] * x[4])
        else
            ∂f_re∂w_i = muladd(U[1, 1], x[1], U[1, 3] * x[3])
            ∂f_im∂w_i = muladd(U[2, 1], x[1], U[2, 3] * x[3])
            ∂g∂w_i = muladd(U[3, 1], x[1], U[3, 4] * x[3])
        end

        # SMatrix has column-major storage
        SMatrix{3,3}(
            ∂f_re∂θ1, ∂f_im∂θ1, ∂g∂θ1,
            ∂f_re∂θ2, ∂f_im∂θ2, ∂g∂θ2,
            ∂f_re∂w_i, ∂f_im∂w_i, ∂g∂w_i)
    end
    out
end


text:
#Define type AMFMmodel as a vector of type AMFMcomp
mutable struct AMFMmodel
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  comps::Vector{AMFMcomp}
end

#Method for type AMFMmodel when called with a time index
function (z::AMFMmodel)(t::Vector{Float64})
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  temp = zeros(size(t))
  for k = 1:length(z.comps)
    temp += z.comps[k](t)
  end
  return temp
end

#Construction of type AMFMmodel from a component set
function AMFMmodel(triplets::Vector{Tuple{Function, Function, Real}})
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
  components = []
  for k = 1:length(triplets)
      push!(components,AMFMcomp(triplets[k]))
  end
  return AMFMmodel(Vector(components))
end

function (z::AMFMmodel)(t::Float64)
  return (z::AMFMmodel)([t::Float64])[1]
end

function (z::AMFMmodel)(t::StepRangeLen)
  return (z::AMFMmodel)(collect(t))
end

#Construct fourierSeries as a type AMFMmodel with a sepecific form
function fourierSeries(T::Real, aₖ::Function, kInds::Array{Int,1}=Vector(-1000:1000))
  components = []
  for k in kInds
    ia(t) = aₖ(k)
    ω(t) = k / T
    push!(components, AMFMcomp(ia, ω, 0))
  end
  return AMFMmodel(components)
end

ext: jl

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
# ---
# title: 884. Uncommon Words from Two Sentences
# id: problem884
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/uncommon-words-from-two-sentences/description/>
# hidden: true
# ---
#
# We are given two sentences `A` and `B`.  (A _sentence_  is a string of space
# separated words.  Each _word_ consists only of lowercase letters.)
#
# A word is _uncommon_  if it appears exactly once in one of the sentences, and
# does not appear in the other sentence.
#
# Return a list of all uncommon words.
#
# You may return the list in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = "this apple is sweet", B = "this apple is sour"
#     Output: ["sweet","sour"]
#
#
# **Example 2:**
#
#
#
#     Input: A = "apple apple", B = "banana"
#     Output: ["banana"]
#
#
#
#
# **Note:**
#
#   1. `0 <= A.length <= 200`
#   2. `0 <= B.length <= 200`
#   3. `A` and `B` both contain only spaces and lowercase letters.
#
#
## @lc code=start
using LeetCode

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


text:
using  NonLinearReactionAdvectionDiffusionWithFrontData;
using  NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity;
using  NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using  NonLinearReactionAdvectionDiffusionWithFrontData: apply_on_dynamic_mesh;

@testset "Градиент на точных данных и статической сетке    " begin

    # Создадим априорную информацию и точное решение для точного q(x)
    a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = NonLinearReactionAdvectionDiffusionWithFrontData.dparams();
    #
    u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
    #' ## Генерация априорной информации
    #########################################################################################
    ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                      # Левый вырожденный корень
    ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true);     # Правый вырожденный корень
    ϕ       = Φ(ϕl, ϕr, N, M);                                              # Серединный корень
    f1_data = f1(ϕ, u, Xₙ, N, M);                                           # Положение переходного слоя
    f2_data = f2(f1_data, u, Xₙ, N, M);                                     # Значение функции на переходном слое
    #########################################################################################
    nothing #hide

    # Функция `heterogeneity` — для внутреннего использования,
    # она принимает аргументы без граничных точек по x.
    Uₙₘ = u[2:N, :];
    X = XX[2:N, :];
    # Создадим замыкания, для удобства
    hmap(z) = [  - heterogeneity(n, m, X[:, m], N, Tₘ, M, Uₙₘ, f1_data, f2_data, z) for n in 1:N-1, m in 1:M+1]
    nnz(arr) = length( findall( x -> ! isapprox(x, 0), arr ))

    # Количество ненулевых элементов — 0
    @test       nnz(hmap(0.0001)) < M/4        # 0, Не рекомендовано
    @test 0 < nnz(hmap(0.0004)) < M/10    # 5, норм!

    # Для визуального контроля, можно использовать следующий код дающий анимацию
    #=
    # {{{
    using Plots, LaTeXStrings;
    pyplot();

    heatmap(h(0.0004)', label=L"- 2 \delta(x - f_1(t)) ( u^s(x,t) - f_2(t))")

    d = [ deltaw( n, f1[m], Xₙ[2:N], N) for n in 1:N-1, m in 1:M+1];    # Значения первого множителя
    uf = [ ( Uₙₘ[n+1, m] - f2[m] ) for n in 1:N-1, m in 1:M+1];         # Значения второго множителя

    yl = extrema( [H; d; uf]  );
    a = Animation()
    for m in 1:M
        nnd = findfirst( x-> x!= 0, d[:,m]);                # Первый ненулевой элемент на m-ом шаге по T
        nnh = findfirst( x-> x!= 0, H[:,m]);                # Первый ненулевой элемент на m-ом шаге по T
        ym = d[nnd, m] > H[nnh, m] ? d[nnd,m] : H[nnh, m]   # Выбор наибольшего ненулевого элемента
                                                            # чтобы нарисовать красивую вертикальную лиинию
        scatter(X,  H[:, m],    label=L"- 2 \delta(x - f_1(t)) ( u^s(x,t) - f_2(t))", ylims = yl)
        scatter!(X, d[:, m],    label=L"\delta(x - f_1(t))")
        scatter!(X, uf[:, m],   label=L"( u^s(x,t) - f_2(t))")
        plot!([Xₙ[nnd+1], Xₙ[nnd+1]], [ym, 0], line=:dash, label="") # Вертикальная линия
        frame(a)
    end
    g = gif(a, "heterogenyry.gif")

    # Карта поверхности неоднородности
    heatmap(X, Tₘ, H',  xlabel = "x", ylabel="t", title=L"- 2 \delta(x - f_1(t)) ( u^s(x,t) - f_2(t))")
    savefig("heterogeneity.png")
    # }}}
    =#

end

@testset "Градиент на точных данных и динамической сетке   " begin

    using NonLinearReactionAdvectionDiffusionWithFrontData;
    using NonLinearReactionAdvectionDiffusionWithFrontData: Φ, apply_on_dynamic_mesh;
    using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;

    a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀, mshfrm = NonLinearReactionAdvectionDiffusionWithFrontData.dparams_nonuniform();
    #
    u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, create_mesh = mshfrm);
    #' ## Генерация априорной информации
    #########################################################################################
    ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                          # Левый вырожденный корень
    ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true);         # Правый вырожденный корень
    ϕ       = Φ(ϕl, ϕr, N, M);                                                  # Полуразность вырожденных корней
    ϕ       = apply_on_dynamic_mesh(ϕ, XX, N, M);                               # Аппроксимация на переменную сетку
    ϕl      = apply_on_dynamic_mesh(ϕl, XX, N, M);                              # Аппроксимация на переменную сетку
    ϕr      = apply_on_dynamic_mesh(ϕr, XX, N, M);                              # Аппроксимация на переменную сетку
    f1_data = f1(ϕ, u, XX, N, M);                                               # Положение переходного слоя
    f2_data = f2(f1_data, u, XX, N, M);                                         # Значение функции на переходном слое
    #########################################################################################

    using  NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity
    # Функция `heterogeneity` — для внутреннего использования,
    # она принимает аргументы без граничных точек по x.
    Uₙₘ = u[2:N, :];
    X = XX[2:N, :];
    hmap(z) = [  - heterogeneity(n, m, X[:, m], N, Tₘ, M, Uₙₘ, f1_data, f2_data, z) for n in 1:N-1, m in 1:M+1]
    nnz(arr) = length( findall( x -> ! isapprox(x, 0), arr ));

    @test nnz(hmap(0.00002)) > M * 0.05         # Много ненулевых элементов,
                                                # значит итерационный процесс продолжится и уйдет от решения
    @test nnz(hmap(0.0000001)) < M * 0.05        # 0, Не рекомендовано
    @test 0 < nnz(hmap(0.000003)) < M * 0.01    # 4, Должно быть норм!

end


text:
###
# This script generates julia bindings for different Vulkan versions

using Clang.wrap_c
using Clang.cindex

# Which variable contains the library later onr
header_library(x) = "libslurm"

slurm_headers = ["/usr/include/slurm/slurm.h"]

const wc = wrap_c.init(;
                        headers = slurm_headers,
                        output_file = "api/slurm.jl",
                        common_file = "api/slurm_common.jl",
                        header_library = header_library,
                        header_wrapped=(x,y)->(println(y);contains(y,"slurm"))
                        )

  # wrap_structs, immutable_structs
  wc.options = wrap_c.InternalOptions(true, true)
  run(wc)


text:
# Xpress model attributes

############################################
#
#   Low level attribute getters
#
############################################

function get_intattr(model::Model, ipar::Int)
    a = Array(Cint, 1)
    ipar = convert(Cint,ipar)
    ret = @xprs_ccall(getintattrib, Cint,
        (Ptr{Void}, Cint, Ptr{Cint}),
        model.ptr_model, ipar, a);
    if ret != 0
        throw( XpressError(model,ret  ) )
    end
    convert(Int, a[1])
end

function get_dblattr(model::Model, ipar::Int)
    a = Array(Float64, 1)
    ipar = convert(Cint,ipar)
    ret = @xprs_ccall(getdblattrib, Cint,
        (Ptr{Void}, Cint, Ptr{Float64}),
        model.ptr_model, ipar, a);
    if ret != 0
        throw(XpressError(model,convert(Int, ret[1])))
    end
    a[1]::Float64
end

function get_strattr(model::Model, ipar::Int)
    a = Array(Cchar, 256)
    ipar = convert(Cint,ipar)
    ret = @xprs_ccall(getstrattrib, Cint,
        (Ptr{Void}, Cint, Ptr{Cchar}),
        model.ptr_model, ipar, a)
    if ret != 0
        throw(XpressError(model,convert(Int, ret[1])))
    end
    String(pointer(a))
end

############################################
#
#   Macros for array definition
#
############################################

macro xprs_int_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_intattr(model, $(attrname))
end

macro xprs_dbl_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_dblattr(model, $(attrname))
end

macro xprs_str_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_strattr(model, $(attrname))
end

############################################
#
#   Model attributes
#
############################################

# basic attributes

#@xprs_str_attr model_name  XPRS_PROBNAME

@xprs_int_attr num_vars     XPRS_COLS
@xprs_int_attr num_constrs  XPRS_ROWS
@xprs_int_attr num_sos      XPRS_SETS
@xprs_int_attr num_qconstrs XPRS_QCONSTRAINTS
@xprs_int_attr num_cnzs     XPRS_ELEMS
@xprs_int_attr num_qnzs     XPRS_QELEMS
@xprs_int_attr num_qcnzs    XPRS_QCELEMS

@xprs_int_attr num_intents  XPRS_MIPENTS

@xprs_int_attr num_setmembers SETMEMBERS

@xprs_dbl_attr obj_sense    XPRS_OBJSENSE
#@xprs_int_attr num_binvars  "NumBinVars"

@xprs_int_attr num_originalvars     XPRS_ORIGINALCOLS
@xprs_int_attr num_originalconstrs  XPRS_ORIGINALROWS

# derived attribute functions

num_linconstrs(model::Model) = num_constrs(model) - num_qconstrs(model)

model_sense(model::Model) = obj_sense(model) == XPRS_OBJ_MINIMIZE ? (:minimize) : (:maximize)

is_qcp(model::Model) = num_qconstrs(model) > 0
is_mip(model::Model) = num_intents(model)+num_sos(model) > 0
is_qp(model::Model) = num_qnzs(model)>0

function model_type(model::Model)
    is_qp(model)  ? (:QP)  :
    is_qcp(model) ? (:QCP) : (:LP)
end

############################################
#
#   The show method for model
#
#   - Based on attributes
#
############################################

function show(io::IO, model::Model)
    if model.ptr_model != C_NULL
        println(io, "Xpress Model:"     )# $(model_name(model))")
        if is_mip(model)
            println(io, "    type   : $(model_type(model)) (MIP)")
        else
            println(io, "    type   : $(model_type(model))")
        end
        println(io, "    sense  : $(model_sense(model))")
        println(io, "    number of variables             = $(num_vars(model))")
        println(io, "    number of linear constraints    = $(num_constrs(model))")
        println(io, "    number of quadratic constraints = $(num_qconstrs(model))")
        println(io, "    number of sos constraints       = $(num_sos(model))")
        println(io, "    number of non-zero coeffs       = $(num_cnzs(model))")
        println(io, "    number of non-zero qp objective terms  = $(num_qnzs(model))")
        println(io, "    number of non-zero qp constraint terms = $(num_qcnzs(model))")
        println(io, "    number of integer entities = $(num_intents(model))")
    else
        println(io, "Xpress Model: NULL")
    end
end

############################################
#
#   more setters and getters
#
############################################
function set_sense!(model::Model, sense::Symbol)
    v = sense == :maximize ? XPRS_OBJ_MAXIMIZE :
        sense == :minimize ? XPRS_OBJ_MINIMIZE :
        throw(ArgumentError("Invalid model sense."))

    ret = @xprs_ccall(chgobjsense, Cint, (
            Ptr{Void},    # model
            Cint          # sense
            ),
            model.ptr_model, v)

    if ret != 0
        throw(XpressError(model))
    end

end

# note: this takes effect only after update_model! is called:
function set_objcoeffs!(model::Model, inds::Vector{Int}, c::Vector)
    n = num_vars(model)
    length(c) == length(inds) || error("Inconsistent argument dimensions.")
    n >= maximum(inds) || error("Inconsistent argument dimensions.")

        ret = @xprs_ccall(chgobj, Cint, (
            Ptr{Void},    # model
            Cint,          # sense
            Ptr{Cint},
            Ptr{Float64}
            ),
            model.ptr_model, Cint( length(c) ), ivec(inds)-1, fvec(c) )

        if ret != 0
            throw(XpressError(model))
        end

end
function set_objcoeffs!(model::Model,c::Vector)
    n = num_vars(model)
    length(c) == n || error("Inconsistent argument dimensions.")

    set_objcoeffs!(model,  collect(1:(n)) , c)
end
set_obj!(model::Model,c::Vector) = set_objcoeffs!(model,c)

function get_lb(model::Model)

    cols = num_vars(model)

    out = Array(Float64,cols)

    ret = @xprs_ccall(getlb, Cint, (
        Ptr{Void},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, cols-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end
lowerbounds(model::Model) = get_lb(model)

function get_ub(model::Model)

    cols = num_vars(model)

    out = Array(Float64,cols)

    ret = @xprs_ccall(getub, Cint, (
        Ptr{Void},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, cols-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end
upperbounds(model::Model) = get_ub(model)

function get_obj(model::Model)

    cols = num_vars(model)

    out = Array(Float64,cols)

    ret = @xprs_ccall(getobj, Cint, (
        Ptr{Void},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, cols-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end
objcoeffs(model::Model) = get_obj(model)

function get_rhs(model::Model)

    rows = num_constrs(model)

    out = Array(Float64,rows)

    ret = @xprs_ccall(getrhs, Cint, (
        Ptr{Void},    # model
        Ptr{Float64},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, rows-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end
function get_rowtype(model::Model)

    rows = num_constrs(model)

    out = Array(Cchar,rows)

    ret = @xprs_ccall(getrowtype, Cint, (
        Ptr{Void},    # model
        Ptr{Cchar},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, rows-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end

function get_coltype(model::Model)

    cols = num_vars(model)

    out = Array(Cchar,cols)

    ret = @xprs_ccall(getcoltype, Cint, (
        Ptr{Void},    # model
        Ptr{Cchar},
        Cint,
        Cint
        ),
        model.ptr_model, out, 0, cols-1)

    if ret != 0
        throw(XpressError(model))
    end

    return out
end

function set_lb!(model::Model,ind::Vector{Int},lb::Vector)

    nbnds = length(ind)

    ret = @xprs_ccall(chgbounds, Cint, (
        Ptr{Void},    # model
        Cint,
        Ptr{Cint},
        Ptr{Cchar},
        Ptr{Float64}
        ),
        model.ptr_model, Cint(nbnds) , ivec(ind)-1, cvecx('L',nbnds), fvec(lb) )

    if ret != 0
        throw(XpressError(model))
    end
end
function set_lb!(model::Model,lb::Vector)

    cols = num_vars(model)
    ( cols == length(lb) ) || error("wrong size of LB vector")

    ind = collect(1:cols)

    set_lb!(model,ind,lb)
end
function set_ub!(model::Model,ind::Vector{Int},ub::Vector)

    nbnds = length(ind)

    ret = @xprs_ccall(chgbounds, Cint, (
        Ptr{Void},    # model
        Cint,
        Ptr{Cint},
        Ptr{Cchar},
        Ptr{Float64}
        ),
        model.ptr_model, Cint(nbnds) , ivec(ind)-1, cvecx('U',nbnds), fvec(ub) )

    if ret != 0
        throw(XpressError(model))
    end
end
function set_ub!(model::Model,ub::Vector)

    cols = num_vars(model)
    ( cols == length(ub) ) || error("wrong size of UB vector")

    ind = collect(1:cols)

    set_ub!(model,ind,ub)
end

function set_rhs!(model::Model,ind::Vector{Int},rhs::Vector)

    nels = length(ind)

    ret = @xprs_ccall(chgrhs, Cint, (
        Ptr{Void},    # model
        Cint,
        Ptr{Cint},
        Ptr{Float64}
        ),
        model.ptr_model, Cint(nels) , ivec(ind)-1, fvec(rhs) )

    if ret != 0
        throw(XpressError(model))
    end
end
function set_rhs!(model::Model, rhs::Vector)
    rows = num_constrs(model)
    set_rhs!(model, collect(1:rows) ,fvec(rhs) )
end

function set_rowtype!(model::Model,senses::Vector)

    rows = num_constrs(model)
    ind = collect(1:rows)
    ret = @xprs_ccall(chgrowtype, Cint, (
        Ptr{Void},    # model
        Cint,
        Ptr{Cint},
        Ptr{Cchar}
        ),
        model.ptr_model, Cint(rows) , ivec(ind)-1, cvec(senses) )

    if ret != 0
        throw(XpressError(model))
    end
end

function set_constrLB!(model::Model, lb)
    senses = get_rowtype(model)
    rhs    = get_rhs(model)
    sense_changed = false
    for i = 1:num_constrs(model)
        if senses[i] == XPRS_GEQ || senses[i] == XPRS_EQ
            # Do nothing
        elseif senses[i] == XPRS_LEQ
            if lb[i] != -Inf
                # LEQ constraint with non-NegInf LB implies a range
                if isapprox(lb[i], rhs[i])
                    # seems to be an equality
                    senses[i] = XPRS_EQ
                    sense_changed = true
                else
                    error("Tried to set LB != -Inf on a LEQ constraint (index $i)")
                end
            else
                lb[i] = rhs[i]
            end
        end
    end
    if sense_changed
        set_rowtype!(model, senses)
    end
    set_rhs!(model, lb)
end

function set_constrUB!(model::Model, ub)
    senses = get_rowtype(model)
    rhs    = get_rhs(model)
    sense_changed = false
    for i = 1:num_constrs(model)
        if senses[i] == XPRS_LEQ || senses[i] == XPRS_EQ
            # Do nothing
        elseif senses[i] == XPRS_GEQ
            if ub[i] != Inf
                # GEQ constraint with non-PosInf UB implies a range
                if isapprox(ub[i], rhs[i])
                    # seems to be an equality
                    senses[i] = XPRS_EQ
                    sense_changed = true
                else
                    error("Tried to set UB != +Inf on a GEQ constraint (index $i)")
                end
            else
              ub[i] = rhs[i]
            end
        end
    end
    if sense_changed
        set_rowtype!(model, senses)
    end
    set_rhs!(model, ub)
end


text:
# Cost functions
"""
    mse(ŷ, y)

Return the mean squared error `sum((ŷ .- y).^2) / length(y)`.
"""
mse(ŷ, y) = sum((ŷ .- y).^2) * 1 // length(y)

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Nothing)
  return -sum(y .* log.(ŷ)) * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Number)
  return -sum(y .* log.(ŷ)) .* weight * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::AbstractVector)
  return -sum(y .* log.(ŷ) .* weight) * 1 // size(y, 2)
end

"""
    crossentropy(ŷ, y; weight=1)

Return the crossentropy computed as `-sum(y .* log.(ŷ) .* weight) / size(y, 2)`.

See also [`logitcrossentropy`](@ref), [`binarycrossentropy`](@ref).
"""
crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight=nothing) = _crossentropy(ŷ, y, weight)

"""
    logitcrossentropy(ŷ, y; weight=1)

Return the crossentropy computed after a [softmax](@ref) operation:

  -sum(y .* logsoftmax(ŷ) .* weight) / size(y, 2)

See also [`crossentropy`](@ref), [`binarycrossentropy`](@ref).
"""
function logitcrossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight = 1)
  return -sum(y .* logsoftmax(ŷ) .* weight) * 1 // size(y, 2)
end

"""
    binarycrossentropy(ŷ, y; ϵ=eps(ŷ))

Return `-y*log(ŷ + ϵ) - (1-y)*log(1-ŷ + ϵ)`. The ϵ term provides numerical stability.

Typically, the prediction `ŷ` is given by the output of a [`sigmoid`](@ref) activation.
"""
binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

"""
    logitbinarycrossentropy(ŷ, y)

`logitbinarycrossentropy(ŷ, y)` is mathematically equivalent to `binarycrossentropy(σ(ŷ), y)`
but it is more numerically stable.

See also [`binarycrossentropy`](@ref), [`sigmoid`](@ref), [`logsigmoid`](@ref).
"""
logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

"""
    normalise(x; dims=1)

Normalises `x` to mean 0 and standard deviation 1, across the dimensions given by `dims`. Defaults to normalising over columns.

```julia-repl
julia> a = reshape(collect(1:9), 3, 3)
3×3 Array{Int64,2}:
  1  4  7
  2  5  8
  3  6  9

julia> normalise(a)
3×3 Array{Float64,2}:
  -1.22474  -1.22474  -1.22474
  0.0       0.0       0.0
  1.22474   1.22474   1.22474

julia> normalise(a, dims=2)
3×3 Array{Float64,2}:
  -1.22474  0.0  1.22474
  -1.22474  0.0  1.22474
  -1.22474  0.0  1.22474
```
"""
function normalise(x::AbstractArray; dims=1)
  μ′ = mean(x, dims = dims)
  σ′ = std(x, dims = dims, mean = μ′, corrected=false)
  return (x .- μ′) ./ σ′
end

"""
    kldivergence(ŷ, y)

KLDivergence is a measure of how much one probability distribution is different from the other.
It is always non-negative and zero only when both the distributions are equal everywhere.
[KL Divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence).
"""
function kldivergence(ŷ, y)
  entropy = sum(y .* log.(y)) *1 //size(y,2)
  cross_entropy = crossentropy(ŷ, y)
  return entropy + cross_entropy
end

"""
    poisson(ŷ, y)

Poisson loss function is a measure of how the predicted distribution diverges from the expected distribution.
[Poisson Loss](https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/poisson).
"""
poisson(ŷ, y) = sum(ŷ .- y .* log.(ŷ)) *1 // size(y,2)

"""
    hinge(ŷ, y)

Measures the loss given the prediction `ŷ` and true labels `y` (containing 1 or -1).
[Hinge Loss](https://en.wikipedia.org/wiki/Hinge_loss).
"""
hinge(ŷ, y) = sum(max.(0, 1 .-  ŷ .* y)) *1 // size(y,2)


text:
# This file contains the implementation of various operators.
# Tests for them is at test/runtests.jl.

using Base
using Statistics
# TODO: we need kwarg support for many of these

# Generic
get_tuple(x) = (x...,)
get_tuple() = nothing
convert_type(x) = Base.convert(Array{Float32, 1}, x)

ops[:Reshape] = function(params, tensor1, shape...)
  if haskey(params, :shape)
    return vcall(:reshape, tensor1, vcall(:broadcast, Int64, vcall(:Tuple, params[:shape])))
  end
  vcall(:reshape, tensor1, vcall(:broadcast, Int64, vcall(:Tuple, vcall(:reverse, shape[1]))))
end

ops[:Conv] = function(params, x, w, b...)
  if (haskey(params, Symbol("auto_pad")))
    if (String(params[:auto_pad]) == "SAME_UPPER" || String(params[:auto_pad] == "SAME_LOWER"))
      temp = Base.convert(Array{Int64,1}, (params[:kernel_shape] .- 1)./2) # Only for strides = [1,1]
      params[:pads] = vcat(temp, temp)                                    # To Do: Add support for other stride values.
    elseif String(params[:auto_pad]) == "VALID"
      params[:pads] = [0,0,0,0]
    end
  end
  if length(params[:pads]) == 4
    params[:pads] = params[:pads][1:2]
  end
    return vcall(:conv4, w, x, Symbol("padding=$(params[:pads])"), Symbol("mode = 1"))
end

ops[:Add] = function(params, A, B)
  s1 = vcall(:size, A)
  s2 = vcall(:size, B)
  if (s1==s2)
    return vcall(:Add, params[:axis], A, B)
  else
    return vcall(:.+, A, B)
  end
end

ops[:Relu] = function(params, x)
  vcall(broadcast, :relu, x)
end

ops[:MaxPool] = function(params, x)
  return vcall(:pool, x, Symbol("mode=0"), Symbol("window=$(params[:kernel_shape])"))
end

ops[:MatMul] = function(params, A, B)
  #tempa = vcall(:permutedims, A, vcall(:reverse, vcall(:range, 1, vcall(:ndims, A))))
  #tempb = vcall(:permutedims, B, vcall(:reverse, vcall(:range, 1, vcall(:ndims, B))))
  vcall(:*, B, A)
end


text:
## make LightGraphs functions work on PowerModels graph directly? i.e.
# function dijkstra_shortest_paths(PMG::PowerModelsGraph)
#     return LightGraphs.dijkstra_shortest_paths(PMG.graph)
# end


text:
import DiffEqBase: init
"""
    init(prob::LEProblem; <keyword arguments>) :: LESolver
    init(prob::CLVProblem; <keyword arguments>) :: CLVSolver

These are simply the aliases of `LESolver(prob; <keyword arguments>)`
and `CLVSolver(prob; <keyword arguments>)`.  See [`LESolver`](@ref)
and [`CLVSolver`](@ref) for supported keyword arguments.
"""
function init end

import DiffEqBase: solve
"""
    solve(prob::LEProblem; <keyword arguments>) :: LESolution
    solve(prob::CLVProblem; <keyword arguments>) :: CLVSolution

Equivalent to:

```julia
solver = init(prob; <keyword arguments except progress>)
solve!(solver; progress=progress)
solver.sol
```
"""
function solve end
using DiffEqBase: solve

import DiffEqBase: solve!
"""
    solve!(solver::LESolver; <keyword arguments>) :: LESolver
    solve!(solver::CLVSolver; <keyword arguments>) :: CLVSolver

Solve pre-initialized problem.
"""
function solve! end
using DiffEqBase: solve!


text:
function Attenuation{T<:AbstractFloat}(nz::Int,nx::Int,ext::Int,atten_max::T)

# This function builds the complex-valued diagonal attenuation operator/matrix,
# where the main diagonal contains the attenuation factor for each point in
# the modeling region AND the surrounding boundary region. It is assumed
# that no attenuation occurs in the modeling regionl. It is also assumed
# that all 4 sides of the modeling region are surrounded by an attenuating
# boundary material to prevent boundary reflections. The imaginary component
# in the attenuation factor is zero in the modeling region and increases
# parabolically outward in the boundary region.
#
# INPUTS:     nz        - Number of grid points in z-direction INCLUDING the absorbing boundary region
#             nx        - Number of grid points in x-direction INCLUDING the absorbing boundary region
#             ext       - Thickness (number of grid points) of attenuating boundary region
#             atten_max - Maximum complex amplitude in the attenuating boundary layer
#
# OUTPUTS:    A         - Diagonal matrix containing the attenuation factors

    atten = (ext:-1:1).^2/(ext.^2)*atten_max
    a_boundary1 = kron(ones(nz,1),atten')
    a_boundary2 = kron(atten,ones(1,nx))
    a = zeros(T,nz,nx)
    a[:,1:ext] += a_boundary1
    a[:,nx-ext+1:nx] += flipdim(a_boundary1,2)
    a[1:ext,:] += a_boundary2
    a[nz-ext+1:nz,:] += flipdim(a_boundary2,1)

    A = 1 - im*a

    return spdiagm((A[:]),(0))

end


text:
@testset "Testing Scheduler with PriorityTime" begin
    s = Scheduler(EmptyState(), PriorityTime{Int, Int})
    ptv = []
    for i in 1:100
        a = rand(1:5)
        b = rand(1:5)
        register!(s, x -> push!(ptv, s.now), PriorityTime(a,b))
    end
    go!(s, 10)

    for i in 1:length(ptv)-1
        a = ptv[i]
        b = ptv[i+1]
        @test a<=b
        @test (a.time < b.time) || (a.time == b.time && a.priority <= b.priority)
    end
end

@testset "Testing +, -, promotion and conversion" begin
    pt1 = PriorityTime(1)
    pt2 = PriorityTime(2.0, 1)
    pt3 = pt1 + pt2
    pt4 = pt1 - pt2
    @test pt3.time == pt1.time + pt2.time
    @test pt3.priority == pt1.priority + pt2.priority
    @test pt4.time == pt1.time - pt2.time
    @test pt4.priority == pt1.priority - pt2.priority

    pt5 = pt1 + 1.0
    @test pt5.time == pt1.time + 1.0
    @test pt5.priority == pt1.priority
    pt6 = 1.0 + pt1
    @test pt6.time == pt1.time + 1.0
    @test pt6.priority == pt1.priority

    pt_if = PriorityTime(1, 1.0)
    pt_bb = PriorityTime(big(1.0), big(1))
    @test pt_if+pt_bb == PriorityTime(big(2.0), big(2.0))
end

@testset "Testing zero" begin
    z1 = zero(PriorityTime{Int, Int16})
    z2 = zero(z1)
    @test z1 == z2
    @test typeof(z1) == typeof(z2)
    @test isa(z1, PriorityTime{Int, Int16})
end

@testset "Testing comparison" begin
    ptc1 = PriorityTime(1,1)
    ptc2 = PriorityTime(1.0,1.0)
    @test isequal(ptc1, ptc2)
    @test hash(ptc1) == hash(ptc2)
    @test ptc1 == ptc2
    @test !(ptc1 != ptc2)
    @test ptc1 <= ptc2
    @test !(ptc1 < ptc2)
    @test ptc1 >= ptc2
    @test !(ptc1 > ptc2)

    ptc3 = PriorityTime(1.0,1.1)
    @test !isequal(ptc1, ptc3)
    @test hash(ptc1) != hash(ptc3)
    @test !(ptc1 == ptc3)
    @test ptc1 != ptc3
    @test ptc1 <= ptc3
    @test ptc1 < ptc3
    @test !(ptc1 >= ptc3)
    @test !(ptc1 > ptc3)

    ptc4 = PriorityTime(1.1,0.1)
    @test !isequal(ptc1, ptc4)
    @test hash(ptc1) != hash(ptc4)
    @test !(ptc1 == ptc4)
    @test ptc1 != ptc4
    @test ptc1 <= ptc4
    @test ptc1 < ptc4
    @test !(ptc1 >= ptc4)
    @test !(ptc1 > ptc4)

    tc4 = Int16(2)
    @test !isequal(ptc1, tc4)
    @test hash(ptc1) != hash(tc4)
    @test !(ptc1 == tc4)
    @test ptc1 != tc4
    @test ptc1 <= tc4
    @test ptc1 < tc4
    @test !(ptc1 >= tc4)
    @test !(ptc1 > tc4)
end


text:
using RunTests
using Base.Test

@testmodule SecondSetOfTests begin

  function test_hello()
    @test true
  end

end


text:
using GLMakie
using ColorTypes

scene = GLMakie.Scene()
width = 640
height = 360
#im = convert.(ColorTypes.RGBA, rand(Float32, width, height))
@time img = [ColorTypes.RGBA(rand(Float32, 3)...) for x in 1:width, y in 1:height]
glim = GLMakie.image!(scene, img, show_axis=false)
makimg = scene.plots[1]

@time makimg[:image][] = img
@run scene.plots[1][:image][] = rand(Float32, width, height)
#glim = scene.plots[1]

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
function test_power_layer(backend::Backend, scale, shift, power, n_input, T, eps)
  tensor_dim = abs(rand(Int)) % 6 + 1
  println("    > scale=$scale, shift=$shift, power=$power, tensor_dim=$tensor_dim")

  dims = [rand(5:9, tensor_dim) for i = 1:n_input]
  input = [rand(T, dims[i]...) for i = 1:n_input]
  input_blob = Blob[make_blob(backend, x) for x in input]
  grad_blob = Blob[make_blob(backend, x) for x in input]

  layer = PowerLayer(tops=Array{Symbol}(n_input), bottoms=Array{Symbol}(n_input),
      scale=scale, shift=shift, power=power)
  state = setup(backend, layer, input_blob, grad_blob)

  forward(backend, state, input_blob)

  for i = 1:n_input
    output = (scale * input[i] + shift) .^ power
    got_output = zeros(T, size(output))
    copy!(got_output, state.blobs[i])

    @test all(-eps .< output - got_output .< eps)
  end

  top_diff = [rand(T, size(input[i])) for i = 1:n_input]
  for i = 1:n_input
    copy!(state.blobs_diff[i], top_diff[i])
  end

  backward(backend, state, input_blob, grad_blob)

  for i = 1:n_input
    grad = power * scale * (scale * input[i] + shift) .^ (power - 1) .* top_diff[i]
    got_grad = zeros(T, size(grad))
    copy!(got_grad, grad_blob[i])
    @test all(-eps .< got_grad - grad .< eps)
  end

  shutdown(backend, state)
end

function simple_rand()
  round(Int, 100*rand())/100
end
function test_power_layer(backend::Backend, n_input, T, eps)
  println("-- Testing PowerLayer on $(typeof(backend)){$T}...")
  test_power_layer(backend, simple_rand(), simple_rand(), 2, n_input, T, eps)
  test_power_layer(backend, 0, simple_rand(), abs(rand(Int)) % 5 + 2, n_input, T, eps)
  test_power_layer(backend, simple_rand(), simple_rand(), 2, n_input, T, eps)
  test_power_layer(backend, simple_rand(), 0, 3, n_input, T, eps)
  test_power_layer(backend, simple_rand(), simple_rand(), 4, n_input, T, eps)

  test_power_layer(backend, simple_rand(), simple_rand(), 0, n_input, T, eps)
  test_power_layer(backend, simple_rand(), simple_rand(), 1, n_input, T, eps)
  test_power_layer(backend, simple_rand(), simple_rand(), -1, n_input, T, eps)
end

function test_power_layer(backend::Backend)
  test_power_layer(backend, 3, Float32, 8e-2)
  test_power_layer(backend, 3, Float64, 1e-8)
end

if test_cpu
  test_power_layer(backend_cpu)
end
if test_gpu
  test_power_layer(backend_gpu)
end


text:
function sumpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CUDA.zeros(T, dims.us_dims[1], c)
    ScatterNNlib.scatter_add!(Y, X, cluster)
    Y
end

sumpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    sumpool(CuArray{Int64}(cluster), X, c)

function subpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CUDA.zeros(T, dims.us_dims[1], c)
    ScatterNNlib.scatter_sub!(Y, X, cluster)
    Y
end

subpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    subpool(CuArray{Int64}(cluster), X, c)

function prodpool(cluster::CuArray{Int}, X::CuArray{T},
                  c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CUDA.ones(T, dims.us_dims[1], c)
    ScatterNNlib.scatter_mul!(Y, X, cluster)
    Y
end

prodpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    prodpool(CuArray{Int64}(cluster), X, c)

function divpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    FT = float(T)
    Y = CUDA.ones(FT, dims.us_dims[1], c)
    ScatterNNlib.scatter_div!(Y, FT.(X), cluster)
    Y
end

divpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    divpool(CuArray{Int64}(cluster), X, c)

function maxpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CUDA.fill(typemin(T), dims.us_dims[1], c)
    ScatterNNlib.scatter_max!(Y, X, cluster)
    Y
end

maxpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    maxpool(CuArray{Int64}(cluster), X, c)

function minpool(cluster::CuArray{Int}, X::CuArray{T},
                 c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    Y = CUDA.fill(typemax(T), dims.us_dims[1], c)
    ScatterNNlib.scatter_min!(Y, X, cluster)
    Y
end

minpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    minpool(CuArray{Int64}(cluster), X, c)

function meanpool(cluster::CuArray{Int}, X::CuArray{T},
                  c::Integer=length(Set(cluster))) where {T<:Real}
    dims = Dims(cluster, X)
    FT = float(T)
    Y = CUDA.zeros(FT, dims.us_dims[1], c)
    ScatterNNlib.scatter_mean!(Y, FT.(X), cluster)
    Y
end

meanpool(cluster::Array{Int}, X::CuArray{T}, c::Integer=length(Set(cluster))) where {T<:Real} =
    meanpool(CuArray{Int64}(cluster), X, c)

@adjoint function prodpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    prodpool(cluster, X), function (Δ)
        rev_cluster = ScatterNNlib.gather_indices(cluster)
        ∇X = ScatterNNlib.gather(zero(Δ)+Δ, cluster)
        @inbounds for ind = CartesianIndices(cluster)
            ind = Tuple(ind)
            inds = [x for x in rev_cluster[cluster[ind...]] if x != ind]
            for i = 1:size(X, 1)
                multiplier = one(T)
                for j = inds
                    multiplier *= X[i, j...]
                end
                ∇X[i, ind...] *= multiplier
            end
        end
        (nothing, ∇X)
    end
end

@adjoint function divpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    divpool(cluster, X), function (Δ)
        rev_cluster = ScatterNNlib.gather_indices(cluster)
        ∇X = -ScatterNNlib.gather(zero(Δ)+Δ, cluster)
        ∇X ./= X.^2
        @inbounds for ind = CartesianIndices(cluster)
            ind = Tuple(ind)
            inds = [x for x in rev_cluster[cluster[ind...]] if x != ind]
            for i = 1:size(X, 1)
                denom = one(T)
                for j = inds
                    denom *= X[i, j...]
                end
                ∇X[i, ind...] /= denom
            end
        end
        (nothing, ∇X)
    end
end

@adjoint function maxpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    max = maxpool(cluster, X)
    max, function (Δ)
       Δu = (X .== ScatterNNlib.gather(max, cluster))
       Δu .*= ScatterNNlib.gather(zero(Δ)+Δ, cluster)
       (nothing, Δu)
    end
end

@adjoint function minpool(cluster::CuArray{Int}, X::CuArray{T}) where {T<:Real}
    min = minpool(cluster, X)
    min, function (Δ)
       Δu = (X .== ScatterNNlib.gather(min, cluster))
       Δu .*= ScatterNNlib.gather(zero(Δ)+Δ, cluster)
       (nothing, Δu)
    end
end


text:
module BOMoD

    using Kronecker
    using StatsBase
    using Stheno
    using Optim
    using Random
    using Statistics
    using Distributions
    using DataFrames
    using Plots
    #using BayesianLinearRegressors: BayesianLinearRegressor, rand, posterior, marginals

    import Base: +, *, getindex,length, eltype,isequal,isless,push!,in
    import Stheno: ew,pw,Kernel
    import LinearAlgebra: norm, eigvals!, Diagonal
    import Optim: minimizer
    import StatsBase: sample

    export Mod, Group_Mod , group_mod ,getspace
    export construct_design
    export No_Constrain, Ordered_Constrain, UnOrdered_Constrain, Possition_Constrain
    export Ordered_Construct, Unordered_Construct
    export Frame_Space,Computed_Space,Full_Ordered_space
    export GP_optimised,gp_optimised
    export thompson_sampling, save_thompson_sampling
    export LevStehnoexp

    # include construction of space
    include(joinpath("design", "Mod.jl"))
    include(joinpath("design", "Combinatroics.jl"))
    include(joinpath("design", "Constrains.jl"))
    include(joinpath("design", "Space.jl"))
    include(joinpath("design", "Construct.jl"))
    include(joinpath("design", "Design_space.jl"))
    include(joinpath("design", "Sample.jl"))
    #BO
    include(joinpath("BO", "Hyper_opt.jl"))
    include(joinpath("BO", "Kernels.jl"))
    include(joinpath("BO", "TS_sampeling.jl"))
    #linear_model
    include(joinpath("linear_model", "linear_model.jl"))

end # module


text:
#!/usr/bin/env julia

module QJuliaSUNtest

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "../core"))
push!(LOAD_PATH, joinpath(@__DIR__, "../main/fields"))
push!(LOAD_PATH, @__DIR__)

import QJuliaRegisters
import QJuliaEnums
import QJuliaGrid
import QJuliaFields
import QJuliaCUDAFields

using CUDAdrv
using CUDAnative
using CuArrays
using BenchmarkTools
using StaticArrays

#Float2 version
@inline function axpy_(a, x::Complex{T}, y::Complex{T})::Complex{T} where T <: AbstractFloat

   rres = CUDAnative.fma(a, real(x), real(y))
   ires = CUDAnative.fma(a, imag(x), imag(y))

   return Complex{T}(rres, ires)
 end

 @inline function caxpy_(a::Complex{T}, x::Complex{T}, y::Complex{T})::Complex{T} where T <: AbstractFloat

    rres = CUDAnative.fma(+real(a), real(x), real(y))
    rres = CUDAnative.fma(-imag(a), imag(x), rres   )

    ires = CUDAnative.fma(+imag(a), real(x), imag(y))
    ires = CUDAnative.fma(+real(a), imag(x), ires   )

    return Complex{T}(rres, ires)
  end

 @inline function axpy_(a::T, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat

   res1 = CUDAnative.fma(a, x[1].value, y[1].value)
   res2 = CUDAnative.fma(a, x[2].value, y[2].value)

   return NTuple{2, VecElement{T}}((res1, res2))
 end

 @inline function caxpy_(a::Complex{T}, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat

   rres = CUDAnative.fma(+real(a), x[1].value, y[1].value)
   rres = CUDAnative.fma(-imag(a), x[2].value, rres      )

   ires = CUDAnative.fma(+real(a), x[2].value, y[2].value)
   ires = CUDAnative.fma(+imag(a), x[1].value, ires      )

   return NTuple{2, VecElement{T}}((res1, res2))
 end

 @inline function axpy_(a::NTuple{2, VecElement{T}}, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat
   return (VecElement(CUDAnative.fma(a[1].value, x[1].value, y[1].value)), VecElement(CUDAnative.fma(a[2].value, x[2].value, y[2].value)))
   #return (VecElement(a[1].value*x[1].value + y[1].value), VecElement(a[2].value*x[1].value + y[2].value))
 end

 @inline function caxpy_(a::NTuple{2, VecElement{T}}, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat
   res  = (VecElement(CUDAnative.fma(+a[1].value, x[1].value, y[1].value)),
           VecElement(CUDAnative.fma(+a[1].value, x[2].value, y[2].value)))
   return (VecElement(CUDAnative.fma(-a[2].value, x[2].value, res[1].value)),
           VecElement(CUDAnative.fma(+a[2].value, x[1].value, res[2].value)))
 end

 @inline function caxpy_(a::NTuple{4, VecElement{T}}, x::NTuple{4, VecElement{T}}, y::NTuple{4, VecElement{T}})::NTuple{4, VecElement{T}} where T <: AbstractFloat
   res  = (VecElement(CUDAnative.fma(+a[1].value, x[1].value, y[1].value)),
           VecElement(CUDAnative.fma(+a[1].value, x[2].value, y[2].value)),
           VecElement(CUDAnative.fma(+a[3].value, x[3].value, y[3].value)),
           VecElement(CUDAnative.fma(+a[3].value, x[4].value, y[4].value)))
   return (VecElement(CUDAnative.fma(-a[2].value, x[2].value, res[1].value)),
           VecElement(CUDAnative.fma(+a[2].value, x[1].value, res[2].value)),
           VecElement(CUDAnative.fma(-a[4].value, x[4].value, res[3].value)),
           VecElement(CUDAnative.fma(+a[4].value, x[3].value, res[4].value)))
 end

 # a few helper methods:
 @inline function get_regN(meta::NTuple{N, VecElement{T}}) where T <: AbstractFloat where N; return N; end
 @inline function get_regT(meta::NTuple{N, VecElement{T}}) where T <: AbstractFloat where N; return T; end

function axpy(a, xaccessor, yaccessor)
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x
    # get data arrays from the accessors
    xdata = xaccessor[1]; ydata = yaccessor[1]
    args  = xaccessor[2]
    #
    if xcb > args.length; return; end

#    if (Int64(blockIdx().x) == 1 && threadIdx().x == 1)
#      for s in 1:args.nSpin
#        for c in 1:args.nColor
#          @cuprintf("color %d  ", c)
#          @cuprintf("spin  %d\n", s)
#        end
#      end
#      @cuprintf("Info block %ld, thread %ld : %d  (length), %d \n", Int64(blockIdx().x), Int64(threadIdx().x), args.length, args.length)
#      @cuprintf("Info block %ld, thread %ld : %d  (colors) \n", Int64(blockIdx().x), Int64(threadIdx().x), args.nColor)
#    end

    y_site_view = view(ydata,xcb,:,:,1)
    x_site_view = view(xdata,xcb,:,:,1)

    for s in 1:args.nSpin
      for c in 1:args.nColor
        y_site_view[c, s] = axpy_(a, x_site_view[c, s], y_site_view[c, s])
      end
    end

  #  if (Int64(blockIdx().x) == 1 && threadIdx().x == 4)
  #    for s in 1:args.nSpin
  #      for c in 1:args.nColor
  #        @cuprintf("Info block %ld, thread %ld : %le, %le!\n", Int64(blockIdx().x), Int64(threadIdx().x), x_site_view[c,s][1].value, x_site_view[c,s][2].value)
  #      end
  #    end
  #  end
    return
end

function set_value(x, a)
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x

    (xlength, xcolors, xspins, xvecs) = size(x)

    if xcb > xlength; return; end

    xview = view(x,:,:,:,1)

    for s in 1:xspins
      for c in 1:xcolors
        xview[xcb, c, s] =  NTuple{2, Float32}((Float32(s), Float32(c)))
      end
    end

    return
end

function sun_test(x, y)
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x
    #    wid, lane = fldmod1(threadIdx().x, CUDAnative.warpsize())
    lid = Int64(threadIdx().x)

    (nElems, nCols, nRows, nDirs) = size(x)

    if xcb > nElems; return; end

    dir = 1

    y_site_view = view(y,xcb,:,:,dir)
    x_site_view = view(x,xcb,:,:,dir)

    #shared_dims = blockDim().x*nCols*nRows*get_regN(x_site_view[1,1])
    #tmp = @cuDynamicSharedMem(get_regT(x_site_view[1,1]), shared_dims)
    msize = nCols*nRows*4 #get_regN(x_site_view[1,1])

    shared_dims = blockDim().x*msize
    tmp = @cuDynamicSharedMem(Float32, shared_dims)

    #if (Int64(blockIdx().x) == 1)
    #  @cuprintf("Info block %ld, thread %ld : %d  (wid), %d (lane), %d (warp size) \n", Int64(blockIdx().x), Int64(threadIdx().x), wid, lane, CUDAnative.warpsize())
    #end

    for r in 1:nRows
      for c in 1:nCols
        temp = NTuple{4, VecElement{Float32}}((0.0f0, 0.0f0, 0.0f0, 0.0f0))
        for j in 1:nCols
          temp = caxpy_(x_site_view[r, j], y_site_view[j, c], temp)
        end
        @inbounds tmp[msize*(lid-1)+(r+c*3)*2+0] = temp[1].value
        @inbounds tmp[msize*(lid-1)+(r+c*3)*2+1] = temp[2].value
        @inbounds tmp[msize*(lid-1)+(r+c*3)*2+2] = temp[3].value
        @inbounds tmp[msize*(lid-1)+(r+c*3)*2+3] = temp[4].value
      end
    end

    for r in 1:nRows
      for c in 1:nCols
        temp = NTuple{4, VecElement{Float32}}((tmp[msize*(lid-1)+(r+c*3)*2+0], tmp[msize*(lid-1)+(r+c*3)*2+1], tmp[msize*(lid-1)+(r+c*3)*2+2], tmp[msize*(lid-1)+(r+c*3)*2+3]))
        y_site_view[r, c] = temp
      end
    end

#    if (Int64(blockIdx().x) == 1 && threadIdx().x == 4)
#      for c in 1:xrows
#        for r in 1:xcols
#          @cuprintf("Info block %ld, thread %ld : %le, %le!\n", Int64(blockIdx().x), Int64(threadIdx().x), y_site_view[r,c][1].value, y_site_view[r,c][2].value)
#        end
#      end
#    end

    return
end

function sun_test_pure_arg_templ(x::CuDeviceArray{NTuple{N,VecElement{T}},4}, y::CuDeviceArray{NTuple{N,VecElement{T}},4}) where T <: AbstractFloat where N
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x
    #    wid, lane = fldmod1(threadIdx().x, CUDAnative.warpsize())
    tid  = Int64(threadIdx().x)
    bdim = Int64(blockDim().x)

    (nElems, nCols, nRows, nDirs) = size(x)

    if xcb > nElems; return; end

    dir = 1

    y_site_view = view(y,xcb,:,:,dir)
    x_site_view = view(x,xcb,:,:,dir)

    #shared_dims = blockDim().x*nCols*nRows*get_regN(x_site_view[1,1])
    #tmp = @cuDynamicSharedMem(get_regT(x_site_view[1,1]), shared_dims)

    shared_dims = bdim*nCols*nRows*N
    tmp = @cuDynamicSharedMem(T, shared_dims)

    #if (Int64(blockIdx().x) == 1)
    #  @cuprintf("Info block %ld, thread %ld : %d  (wid), %d (lane), %d (warp size) \n", Int64(blockIdx().x), Int64(threadIdx().x), wid, lane, CUDAnative.warpsize())
    #end

    for r in 1:nRows
      for c in 1:nCols
        temp = NTuple{N, VecElement{T}}(ntuple(i->T(0.0), N))
        for j in 1:nCols
          temp = caxpy_(x_site_view[r, j], y_site_view[j, c], temp)
        end
        offset = tid + ((r-1)+(c-1)*3)*bdim
        for l in 1:N
          @inbounds tmp[offset+(l-1)*9*bdim] = temp[l].value
        end
      end
    end

    for r in 1:nRows
      for c in 1:nCols
        offset = tid + ((r-1)+(c-1)*3)*bdim
        temp = N == 2 ? NTuple{N, VecElement{T}}( (tmp[offset+0*9*bdim], tmp[offset+1*9*bdim])) : NTuple{N, VecElement{T}}( (tmp[offset+0*9*bdim], tmp[offset+1*9*bdim],tmp[offset+2*9*bdim], tmp[offset+3*9*bdim]) )
        @inbounds y_site_view[r, c] = temp
      end
    end

#    if (Int64(blockIdx().x) == 1 && threadIdx().x == 4)
#      for c in 1:xrows
#        for r in 1:xcols
#          @cuprintf("Info block %ld, thread %ld : %le, %le!\n", Int64(blockIdx().x), Int64(threadIdx().x), y_site_view[r,c][1].value, y_site_view[r,c][2].value)
#        end
#      end
#    end

    return
end

function sun_extended(z::CuDeviceArray{NTuple{N, VecElement{T}},4},  x::CuDeviceArray{NTuple{N,VecElement{T}},4}, y::CuDeviceArray{NTuple{N,VecElement{T}},4}) where T <: AbstractFloat where N
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x

    (nElems, nCols, nRows, nDirs) = size(x)

    if xcb > nElems; return; end

    dir = 1

    z_site_view = view(z,xcb,:,:,dir)
    y_site_view = view(y,xcb,:,:,dir)
    x_site_view = view(x,xcb,:,:,dir)

    #if (Int64(blockIdx().x) == 1)
    #  @cuprintf("Info block %ld, thread %ld : %d  (wid), %d (lane), %d (warp size) \n", Int64(blockIdx().x), Int64(threadIdx().x), wid, lane, CUDAnative.warpsize())
    #end

    for r in 1:nRows
      for c in 1:nCols
        temp = NTuple{N, VecElement{T}}(ntuple(i->T(0.0), N))
        for j in 1:nCols
          temp = caxpy_(x_site_view[r, j], y_site_view[j, c], temp)
        end
        z_site_view[r,c] = temp
      end
    end

    return
 end

function sun_test_templ(x::Tuple{CuDeviceArray{NTuple{N,VecElement{T}},4},QJuliaFields.AccessorArgs}, y::Tuple{CuDeviceArray{NTuple{N,VecElement{T}},4},QJuliaFields.AccessorArgs}) where T <: AbstractFloat where N
    xcb = (blockIdx().x-1) * blockDim().x + threadIdx().x
    #    wid, lane = fldmod1(threadIdx().x, CUDAnative.warpsize())
    tid  = Int64(threadIdx().x)
    bdim = Int64(blockDim().x)

    (nElems, nCols, nRows, nDirs) = size(x[1])

    if xcb > nElems; return; end

    dir = 1

    y_site_view = view(y[1],xcb,:,:,dir)
    x_site_view = view(x[1],xcb,:,:,dir)

    #shared_dims = blockDim().x*nCols*nRows*get_regN(x_site_view[1,1])
    #tmp = @cuDynamicSharedMem(get_regT(x_site_view[1,1]), shared_dims)

    shared_dims = bdim*nCols*nRows*N
    tmp = @cuDynamicSharedMem(T, shared_dims)

    #if (Int64(blockIdx().x) == 1)
    #  @cuprintf("Info block %ld, thread %ld : %d  (wid), %d (lane), %d (warp size) \n", Int64(blockIdx().x), Int64(threadIdx().x), wid, lane, CUDAnative.warpsize())
    #end

    for r in 1:nRows
      for c in 1:nCols
        temp = NTuple{N, VecElement{T}}(ntuple(i->T(0.0), N))
        for j in 1:nCols
          temp = caxpy_(x_site_view[r, j], y_site_view[j, c], temp)
        end
        offset = tid + ((r-1)+(c-1)*3)*bdim
        for l in 1:N
          @inbounds tmp[offset+(l-1)*9*bdim] = temp[l].value
        end
      end
    end

    for r in 1:nRows
      for c in 1:nCols
        offset = tid + ((r-1)+(c-1)*3)*bdim
        temp = N == 2 ? NTuple{N, VecElement{T}}( (tmp[offset+0*9*bdim], tmp[offset+1*9*bdim])) : NTuple{N, VecElement{T}}( (tmp[offset+0*9*bdim], tmp[offset+1*9*bdim],tmp[offset+2*9*bdim], tmp[offset+3*9*bdim]) )
        y_site_view[r, c] = temp
      end
    end

#    if (Int64(blockIdx().x) == 1 && threadIdx().x == 4)
#      for c in 1:xrows
#        for r in 1:xcols
#          @cuprintf("Info block %ld, thread %ld : %le, %le!\n", Int64(blockIdx().x), Int64(threadIdx().x), y_site_view[r,c][1].value, y_site_view[r,c][2].value)
#        end
#      end
#    end

    return
end

  const  nthreads = 128

  #data_type = ComplexF32
  #data_type = QJuliaRegisters.float2
  data_type = QJuliaRegisters.float4

  const N = 32
  println("Test lattice size ::" , N)

  csGrid   = QJuliaGrid.QJuliaGridDescr_qj{data_type}(QJuliaEnums.QJULIA_CUDA_FIELD_LOCATION, 0, (N,N,N,N,))

  # compilation time
  gaugeParam  = QJuliaGrid.CreateGaugeParams( csGrid )

  cuda_sun_m1 = QJuliaCUDAFields.CreateGenericField( gaugeParam )
  cuda_sun_m2 = QJuliaCUDAFields.CreateGenericField( gaugeParam )
  cuda_sun_m3 = QJuliaCUDAFields.CreateGenericField( gaugeParam )

  println("\n=====GAUGE FIELD INFO=====")
  QJuliaFields.field_info(cuda_sun_m1)

  println("====================")

  cuda_sun_m1_accessor = QJuliaFields.create_field_accessor(cuda_sun_m1)
  cuda_sun_m2_accessor = QJuliaFields.create_field_accessor(cuda_sun_m2)
  cuda_sun_m3_accessor = QJuliaFields.create_field_accessor(cuda_sun_m3)

  accessor_dims = size(cuda_sun_m1_accessor[1])

  println("  ::: ", typeof(cuda_sun_m1_accessor))

  len = accessor_dims[1]

  nblocks= ceil.(Int, len ./ nthreads)
  # calculate size of dynamic shared memory
  regN    = QJuliaRegisters.register_size(typeof(cuda_sun_m1_accessor[1][1]))
  regT    = QJuliaRegisters.register_type(typeof(cuda_sun_m1_accessor[1][1]))

  ncolors = accessor_dims[2]

  shmem_bytes = ncolors*ncolors*regN* nthreads * sizeof(regT)

  println("Execution config for length ", len, " : ", nthreads , " threads, ", nblocks, " blocks , shared mem ", shmem_bytes )

  # compilation options
#code inspections:
#@device_code_typed
#@device_code_llvm
#@device_code_ptx
#@device_code_sass
  @device_code_ptx @cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_test_pure_arg_templ(cuda_sun_m1_accessor[1], cuda_sun_m2_accessor[1])
  #@time @cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_test_templ(cuda_sun_m1_accessor, cuda_sun_m2_accessor)
  #@device_code_ptx @cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_extended(cuda_sun_m3_accessor[1],  cuda_sun_m1_accessor[1], cuda_sun_m2_accessor[1])

  # execution time
  niters = 10000
  secs = CUDAdrv.@elapsed begin # begin profile
  for i in 1:niters
    @cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_test_pure_arg_templ(cuda_sun_m1_accessor[1], cuda_sun_m2_accessor[1])
    #@cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_test_templ(cuda_sun_m1_accessor, cuda_sun_m2_accessor)
    #@cuda blocks=nblocks threads=nthreads shmem=shmem_bytes sun_extended(cuda_sun_m3_accessor[1], cuda_sun_m1_accessor[1], cuda_sun_m2_accessor[1])
  end
  end # end profile

  vol  = regN*(len / 2.0) #complex array length
  ncol = 3
  bytes_per_link_dir = 3.0*vol*length(cuda_sun_m1.v) / 4.0
  footprint_per_link_dir = 2.0*vol*length(cuda_sun_m1.v) / 4.0
  flops = (8.0*ncol - 2.0)*ncol*ncol*vol

  secs_per_iter = (secs / Float64(niters))

  println("Bytes :: ", bytes_per_link_dir, " footprint :: ", footprint_per_link_dir, " flops ",  flops / secs_per_iter)

end # module

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
using CARMA
using Test

# Ideas for tests to implement Any Day Now (TM)
#
# 1. Check against full-matrix likelihood.
# 2. Check that init returns always within prior.
# 3. Check to_params/to_array are inverses.

@testset "CARMA Tests" begin
    include("TestCARMAKepler.jl")
    include("TestCelerite.jl")
    include("TestKalman.jl")
    include("TestLowPass.jl")
end


text:
using Test
using ParametricInversion
using IRTools

function bubblesort(arr)
  n = length(arr)
  for i=1:n, j=1:n-i
      if (arr[j] > arr[j+1])
          arr[j], arr[j+1] = arr[j+1], arr[j]
      end
  end
  arr
end

function test()
  ir = IRTools.explicitbranch!(@code_ir bubblesort(rand(5)))
  ParametricInversion.passvars!(ir)
end


text:
module EnergyExpressions

using AtomicLevels
using LinearAlgebra
using SparseArrays
using Combinatorics
using UnicodeFun
using Formatting
using ProgressMeter

if VERSION < v"1.1-DEV"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
    Base.:(:)(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
        CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

include("conjugate_orbitals.jl")
include("slater_determinants.jl")
include("show_coefficients.jl")
include("nbody_operators.jl")
include("loop_macros.jl")
include("minors.jl")
include("nbody_matrix_elements.jl")
include("bit_configurations.jl")
include("common_operators.jl")
include("equations.jl")
include("variations.jl")
include("multi_configurational_equations.jl")
include("invariant_sets.jl")

end # module


text:
module RPOMDPModels

using RPOMDPs
using RPOMDPToolbox
using Distributions, SimpleProbabilitySets
using JuMP, Clp
using StaticArrays
using AutoHashEquals
using StatsBase
using RCall

importall RPOMDPs

import Base.rand!
import Base.rand
import Base.==
import Base.hash

import RPOMDPs: initial_state, generate_s, generate_o, generate_sor
import RPOMDPs: observation

include("Baby2POMDP.jl")
export
    Baby2POMDP,
    Baby2RPOMDP,
    Baby2IPOMDP,
    Baby2RIPOMDP

include("BabyRobust.jl")
export
    BabyPOMDP,
    BabyRPOMDP,
    BabyIPOMDP,
    BabyRIPOMDP

include("BabyRobustInfo.jl")
export
    BabyInfoPOMDP,
    BabyInfoRPOMDP

include("SimpleBaby2Robust.jl")
export
    SimpleBaby2RPOMDP

include("SimpleTigerRobust.jl")
export
    SimpleTigerRPOMDP

include("TigerRobust.jl")
export
    TigerPOMDP,
    TigerRPOMDP,
    TigerIPOMDP,
    TigerRIPOMDP

include("TigerRobustInfo.jl")
export
    TigerInfoPOMDP,
    TigerInfoRPOMDP

include("CyberRobust.jl")
export
    CyberPOMDP,
    CyberRPOMDP,
    CyberIPOMDP,
    CyberRIPOMDP

include("CyberRobustMisspec.jl")
export
    CyberTestIPOMDP

include("SimpleRobust.jl")
export
    SimpleIPOMDP,
    SimpleRIPOMDP

include("RockDiagnosis.jl")
export
    RockIPOMDP,
    RockRIPOMDP

export
    n_states,
    n_actions,
    n_observations,
    state_index,
    action_index,
    observation_index,
    observation,
    dynamics,
    reward,
    rewardalpha,
    transition,
    observation,
    action,
    initial_belief,
    generate_sor,
    generate_sor_worst,
    initial_state_distribution,
    initial_belief_distribution

end # module


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    delete_text(
        text::Cstring
    )::Nothing

Free a string allocated by the Tesseract library.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | text |         | The string to call Tesseract to free.

__Details:__

This method is meant for internal use.  This library only exports normal Julia objects that can be
garbage collected like normal.
"""
function delete_text(
            text::Cstring
        )::Nothing
    ccall((:TessDeleteText, TESSERACT), Cstring, (Cstring,), text)
    nothing
end

# =================================================================================================
"""
    delete_text(
        text::Ptr{CUInt8}
    )::Nothing

Free a byte array allocated by the Tesseract library.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | text |         | The string to call Tesseract to free.

__Details:__

This method is meant for internal use.  This library only exports normal Julia objects that can be
garbage collected like normal.
"""
function delete_text(
            text::Ptr{UInt8}
        )::Nothing
    ccall((:TessDeleteText, TESSERACT), Cstring, (Ptr{UInt8},), text)
    nothing
end

# =================================================================================================
"""
    delete_array(
        array::Ptr{Ptr{UInt8}}
    )::Nothing

Free an array of strings allocated by the Tesseract library.

__Arguments:__

| T | Name  | Default | Description
|---| :---- | :------ | :----------
| R | array |         | The array of strings to free.

__Details:__

This method is meant for internal use.  This library only exports normal Julia objects that can be
garbage collected like normal.
"""
function delete_array(
            array::Ptr{Ptr{UInt8}}
        )::Nothing
    ccall((:TessDeleteTextArray, TESSERACT), Cvoid, (Ptr{Ptr{UInt8}},), array)
    nothing
end

# =================================================================================================
"""
    delete_array(
        array::Ptr{Cint}
    )::Nothing

Free an array of integers allocated by the Tesseract library.

__Arguments:__

| T | Name  | Default | Description
|---| :---- | :------ | :----------
| R | array |         | The array to ask Tesseract to free.

__Details:__

This method is meant for internal use.  This library only exports normal Julia objects that can be
garbage collected like normal.
"""
function delete_array(
            array::Ptr{Cint}
        )::Nothing
    ccall((:TessDeleteIntArray, TESSERACT), Cvoid, (Ptr{Cint},), array)
    nothing
end


text:
module AngularMomentumAlgebra

using AtomicLevels
using EnergyExpressions
import EnergyExpressions: QuantumOperator, NBodyTermFactor, NBodyTerm,
    OrbitalMatrixElement, NBodyMatrixElement
using LinearAlgebra
using SparseArrays
using WignerSymbols
using HalfIntegers
using UnicodeFun

include("common.jl")
include("clebsch_gordan.jl")
include("linear_combinations.jl")
include("orbitals.jl")
include("tensors.jl")
include("tensor_dsl.jl")
include("tensor_matrix_elements.jl")
include("cartesian.jl")
include("angular_momenta.jl")
include("spherical_tensors.jl")
include("gradients.jl")
include("coulomb.jl")
include("multipole_expansion.jl")
include("energy_expressions.jl")

export # Re-exports from EnergyExpressions.jl
    Conjugate, OrbitalOverlap, EnergyExpression,
    FieldFreeOneBodyHamiltonian, OneBodyHamiltonian, CoulombInteraction, CoulombPotential,
    # Re-exports from LinearAlgebra
    ⋅, dot,
    # Own types
    CoulombInteractionMultipole, CoulombPotentialMultipole

end # module


text:
using Test
import Thermodynamics
const TD = Thermodynamics
const TP = TD.Parameters

import UnPack
import BenchmarkTools

import RootSolvers
const RS = RootSolvers

import CLIMAParameters
const CP = CLIMAParameters
const FT = Float64
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
aliases = string.(fieldnames(TP.ThermodynamicsParameters))
param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
const param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)

ArrayType = Array{FT}
profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
UnPack.@unpack e_int, T, ρ, p, θ_liq_ice, q_tot = profiles

function thermo_state_ρeq()
    ts = TD.PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot)
end

function thermo_state_pθq()
    ts = TD.PhaseEquil_pθq.(param_set, p, θ_liq_ice, q_tot)
end

function thermo_state_pTq()
    ts = TD.PhaseEquil_pTq.(param_set, p, T, q_tot)
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra

# ## system
gravity=-9.81
timestep = 0.1
env = get_environment(:cartpole,
    representation=:minimal,
    timestep=timestep,
    gravity=gravity);

## mujoco_inertia!(env.mechanism)

# ## visualizer
open(env.vis)

# ## dimensions
n = env.num_states
m = env.num_inputs

# ## states
z1 = maximal_to_minimal(env.mechanism, Dojo.cartpole_nominal_max())
zT = maximal_to_minimal(env.mechanism, Dojo.cartpole_goal_max())

# ## horizon
T = 26

# ## model
dyn = IterativeLQR.Dynamics(
    (y, x, u, w) -> dynamics(y, env, x, u, w),
    (dx, x, u, w) -> dynamics_jacobian_state(dx, env, x, u, w),
    (du, x, u, w) -> dynamics_jacobian_input(du, env, x, u, w),
    n, n, m)
model = [dyn for t = 1:T-1]

# ## rollout
ū = [t < 5 ? 1.0 * rand(m) : (t < 10 ? -1.0 * rand(m) : zeros(m)) for t = 1:T-1]
x̄ = IterativeLQR.rollout(model, z1, ū)
visualize(env, x̄)

# ## objective
ot = (x, u, w) -> transpose(x - zT) * Diagonal(1.0e-1 * ones(n)) * (x - zT) + transpose(u) * Diagonal(1.0e-3 * ones(m)) * u
oT = (x, u, w) -> transpose(x - zT) * Diagonal(1.0e-1 * ones(n)) * (x - zT)

ct = IterativeLQR.Cost(ot, n, m)
cT = IterativeLQR.Cost(oT, n, 0)
obj = [[ct for t = 1:T-1]..., cT]

# ## constraints
function goal(x, u, w)
    x - zT
end

cont = IterativeLQR.Constraint()
conT = IterativeLQR.Constraint(goal, n, 0)
cons = [[cont for t = 1:T-1]..., conT]

# ## solver
s = IterativeLQR.solver(model, obj, cons,
    opts=IterativeLQR.Options(
        max_al_iter=10,
        verbose=true))
IterativeLQR.initialize_controls!(s, ū)
IterativeLQR.initialize_states!(s, x̄)

# ## solve
@time IterativeLQR.solve!(s)

# ## solution
z_sol, u_sol = IterativeLQR.get_trajectory(s)
@show IterativeLQR.eval_obj(s.m_data.obj.costs, s.m_data.x, s.m_data.u, s.m_data.w)
@show s.s_data.iter[1]
@show norm(goal(s.m_data.x[T], zeros(0), zeros(0)), Inf)

# ## visualize
visualize(env, z_sol)


text:
using AbstractGPs
using BayesianLinearRegressors
using Distributions
using FiniteDifferences
using LinearAlgebra
using Random
using Test
using Zygote

using BayesianLinearRegressors: BayesianLinearRegressor, posterior, marginals, cov, mean
using FiniteDifferences: j′vp

@testset "BayesianLinearRegressors" begin
    include("bayesian_linear_regression.jl")
end


text:
# Module BLAST
# ============
#
# A module for running command line BLAST and parsing BLAST output files.

module BLAST

export
    blastn,
    blastp,
    readblastXML,
    BLASTResult

import BioAlignments: AlignedSequence
import BioSequences: BioSequence, DNASequence, AminoAcidSequence
import EzXML

include("blastcommandline.jl")

end # module BLAST


text:
using SubglacialPlumes
using Test

#params = Params( )
@testset "ExampleTests" begin
    #2x + 3y
    @test SubglacialPlumes.f(2,1) == 7
    @test SubglacialPlumes.f(2,3) == 13
end

max_stars_repo_name: alextbradley/SubglacialPlumes

max_issues_repo_name: alextbradley/SubglacialPlumes

max_forks_repo_name: alextbradley/SubglacialPlumes


text:
using LinearAlgebra

@testset "Test QXZoo Circ" begin

    @testset "Test circuit creation" begin
        @testset "Test empty circuit" begin
            @test begin
                cct = QXZoo.Circuit.Circ()
                cct.num_qubits == 0
                cct.circ_ops.len == 0
                length(cct.gate_set) == 0
            end
        end

        @testset "Test gate n-qubit circuit" begin
            cct = QXZoo.Circuit.Circ(5)
            cct << QXZoo.DefaultGates.x(1)
            @test length(cct.gate_set) == 1
            @test cct.circ_ops.len == 1
            @test cct.num_qubits == 5
            @test_throws ErrorException cct << QXZoo.DefaultGates.x(6)
            @test_throws ErrorException cct << QXZoo.DefaultGates.c_x(1,6)
            @test_throws ErrorException cct << QXZoo.DefaultGates.c_x(6,1)
        end

        @testset "Test circuit iteration" begin
            cct1 = QXZoo.Circuit.Circ(5)

            for i in 1:5
                cct1 << QXZoo.DefaultGates.x(i)
            end

            for i in cct1.circ_ops
                @test i.gate_symbol === QXZoo.DefaultGates.GateSymbols.x
            end
        end

        @testset "Test circuit IO" begin
            cct = QXZoo.Circuit.Circ(5)

            cct << QXZoo.DefaultGates.x(1)
            cct << QXZoo.DefaultGates.x(2)
            cct << QXZoo.DefaultGates.h(3)
            cct << QXZoo.DefaultGates.c_z(3,1)

            QXZoo.Circuit.export_circuit(cct, "test_circ")
            @test isfile("test_circ.jld2") === true

            data = QXZoo.Circuit.load_circuit("test_circ")
            cct2 = data["circuit"]

            @test cct.circ_ops.len == cct2.circ_ops.len
            @test cct.gate_set == cct2.gate_set

            h1 = cct.circ_ops.head
            h2 = cct2.circ_ops.head

            for i in 1:cct.circ_ops.len
                @test h1.data === h2.data
                h1 = h1.next
                h2 = h2.next
            end
        end

        @testset "Test defined gateset" begin
            gates = Set([QXZoo.DefaultGates.GateSymbols.x, QXZoo.DefaultGates.GateSymbols.h, QXZoo.DefaultGates.GateSymbols.r_phase(pi/5)])
            cct = QXZoo.Circuit.Circ(5, gates)
            @test cct.num_qubits == 5
            @test cct.gate_set == gates
        end

        @testset "Test circuit combine" begin
            cct1 = QXZoo.Circuit.Circ(5)
            cct2 = QXZoo.Circuit.Circ(5)
            cct3 = QXZoo.Circuit.Circ(5)
            cct4 = QXZoo.Circuit.Circ(6)

            for i in 1:5
                cct1 << QXZoo.DefaultGates.x(i)
                cct3 << QXZoo.DefaultGates.x(i)
                cct4 << QXZoo.DefaultGates.x(i)
            end
            for i in 1:4
                cct2 << QXZoo.DefaultGates.c_x(i, i+1)
                cct3 << QXZoo.DefaultGates.c_x(i, i+1)
                cct4 << QXZoo.DefaultGates.c_x(i, i+1)
            end

            @test cct3.circ_ops.len == cct1.circ_ops.len + cct2.circ_ops.len
            @test cct3.gate_set == union(cct1.gate_set, cct2.gate_set)

            cct1 << cct2

            @test cct3.circ_ops.len == cct1.circ_ops.len
            @test cct2.circ_ops.len == 0
            @test cct3.gate_set == cct1.gate_set

            h1 = cct1.circ_ops.head
            h3 = cct3.circ_ops.head

            for i in 1:cct1.circ_ops.len
                @test h1.data === h3.data
                h1 = h1.next
                h3 = h3.next
            end

            @test_throws ErrorException cct3 << cct4

        end
    end
end


text:
# example-4.jl
# Example nested simulation models with data collection

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

run_number = 0
const GENERATOR_RNG = 1
const CUSTOMER_RNG = 2

"Process to generate n customers arriving into the system."
@process generator(n::Integer) begin
    rng = MersenneTwister(seed(run_number, GENERATOR_RNG))
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    for i = 1:n
        work(rand(rng, dist))
        @schedule now customer(i)
    end
end

"The ith customer into the system."
@process customer(i::Integer) begin
    rng = MersenneTwister(seed(run_number, CUSTOMER_RNG, i))
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    @with_resource tellers begin
        work(rand(rng, dist))
    end
end

"Run the simulation for n customers."
function run_simulation(n1::Integer, n2::Integer)
    println("SimLynx Repeatable Open Loop Processing Example")
    println("Number of runs = $n1, number of customers = $n2")
    println("Single Queue, $N_TELLERS Teller Bank Model")
    @simulation begin
        avg_wait = Variable{Float64}(data=:tally, history=true)
        for i = 1:n1
            global run_number = i
            @simulation begin
                global tellers = Resource(N_TELLERS, "tellers")
                @schedule at 0.0 generator(n2)
                start_simulation()
                avg_wait.value = tellers.wait.stats.mean
            end
        end
        print_stats(avg_wait, title="Average Wait Statistics")
        plot_history(avg_wait, title="Average Wait History")
    end
end

run_simulation(10_000, 100)


text:
mutable struct Clock{T<:Number}
  time::T
  iteration::Int
end


text:
module ReadQNM

export qnm

include("../src/CustomTypes.jl")
using .CustomTypes

export qnm

function eq(x,y)
   return abs(x-y)<1e-14 ? true : false
end

"""
    qnm(n::myI,s::myI,m::myI,l::myI,a::myF)

Returns omega, lambda
"""
function qnm(n::myI,s::myI,m::myI,l::myI,a::myF)

   for line in eachline("qnmvals.txt")
      line = split(line,",")
      nl = parse(myF,line[1])
      sl = parse(myF,line[2])
      ml = parse(myF,line[3])
      ll = parse(myF,line[4])
      al = parse(myF,line[5])

      if eq(n,nl) && eq(s,sl) && eq(m,ml) && eq(l,ll) && eq(a,al)
         return parse(myF,line[6])+im*parse(myF,line[7]), parse(myF,line[8])+ im*parse(myF,line[9])
      end
   end
   return NaN, NaN
end

end


text:
function exponential_fg!(∇f, x)

    if !(∇f==nothing)
        ∇f[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
        ∇f[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
    end

    fx = exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)

    return fx
end

function exponential_hessian!(storage, x)
    storage[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    storage[1, 2] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
end

tp_fletch_powell_fg! = TestProblem(exponential_fg!, [0.0, 0.0], I, NearlyNewton.OptOptions())
tp_fletch_powell_fg!_alt = TestProblem(exponential_fg!, [0.0, 0.0], I, NearlyNewton.OptOptions())


text:
# ======================================
# Construction of the linearized system
# ======================================

"""
    build_matrix(F₁, F₂, N)

Compute the Carleman linearization matrix associated to the quadratic
system ``x' = F₁x + F₂(x⊗x)``, truncated at order ``N``.

### Input

- `F₁` -- sparse matrix of size ``n × n``
- `F₂` -- sparse matrix of size ``n × n^2``
- `N`  -- integer representing the truncation order, should be at least two

### Output

Sparse matrix `A`.

### Algorithm

See references [1] and [2] of CARLIN.md.
"""
function build_matrix(F₁, F₂, N)
    if N < 1
        throw(ArgumentError("expected the truncation order to be at least 1, got N=$N"))
    elseif N == 1
        _build_matrix_N1(F₁, F₂)
    elseif N == 2
        _build_matrix_N2(F₁, F₂)
    elseif N == 3
        _build_matrix_N3(F₁, F₂)
    elseif N == 4
        _build_matrix_N4(F₁, F₂)
    else
        _build_matrix_N(F₁, F₂, N) # general case
    end
end

function _build_matrix_N1(F₁, F₂)
    return F₁
end

function _build_matrix_N2(F₁, F₂)
    n = size(F₁, 1)
    a = hcat(kron_sum(F₁, 1), kron_sum(F₂, 1))
    b = hcat(spzeros(n^2, n), kron_sum(F₁, 2))
    return vcat(a, b)
end

function _build_matrix_N3(F₁, F₂)
    n = size(F₁, 1)
    a = hcat(kron_sum(F₁, 1), kron_sum(F₂, 1), spzeros(n, n^3))
    b = hcat(spzeros(n^2, n), kron_sum(F₁, 2), kron_sum(F₂, 2))
    c = hcat(spzeros(n^3, n), spzeros(n^3, n^2), kron_sum(F₁, 3))
    return vcat(a, b, c)
end

function _build_matrix_N4(F₁, F₂)
    n = size(F₁, 1)
    a = hcat(kron_sum(F₁, 1), kron_sum(F₂, 1), spzeros(n, n^3), spzeros(n, n^4))
    b = hcat(spzeros(n^2, n), kron_sum(F₁, 2), kron_sum(F₂, 2), spzeros(n^2, n^4))
    c = hcat(spzeros(n^3, n), spzeros(n^3, n^2), kron_sum(F₁, 3), kron_sum(F₂, 3))
    d = hcat(spzeros(n^4, n), spzeros(n^4, n^2), spzeros(n^4, n^3), kron_sum(F₁, 4))
    return vcat(a, b, c, d)
end

function _build_matrix_N(F₁, F₂, N)
    @assert N >= 3 "expected N to be at least 3, got N=$N"
    n = size(F₁, 1)

    out = Vector{typeof(F₁)}()
    for j in 1:N
        if j == N
            F12_j = kron_sum(F₁, j)
            Zleft_j = spzeros(n^j, sum(n^i for i in 1:j-1))
            push!(out, hcat(Zleft_j, F12_j))
        else
            F12_j = hcat(kron_sum(F₁, j), kron_sum(F₂, j))
            if j == 1
                Zright_j = spzeros(n^j, sum(n^i for i in j+2:N))
                push!(out, hcat(F12_j, Zright_j))
            elseif j == N-1
                Zleft_j = spzeros(n^j, sum(n^i for i in 1:j-1))
                push!(out, hcat(Zleft_j, F12_j))
            else # general case
                Zleft_j = spzeros(n^j, sum(n^i for i in 1:j-1))
                Zright_j = spzeros(n^j, sum(n^i for i in j+2:N))
                push!(out, hcat(Zleft_j, F12_j, Zright_j))
            end
        end
    end
    return reduce(vcat, out)
end


text:
module MultiFEOperators

using Gridap
import Gridap.FEOperators: LinearFEOperator
import Gridap.FEOperators: NonLinearFEOperator

function LinearFEOperator(
  biform::Function,
  liform::Function,
  testfesp::Vector{<:FESpaceWithDirichletData},
  trialfesp::Vector{<:FESpaceWithDirichletData},
  assem::MultiAssembler,
  trian::Triangulation{Z},
  quad::CellQuadrature{Z}) where Z
  V = MultiFESpace(testfesp)
  U = MultiFESpace(trialfesp)
  LinearFEOperator(biform,liform,V,U,assem,trian,quad)
end

function LinearFEOperator(
  testfesp::Vector{<:FESpaceWithDirichletData},
  trialfesp::Vector{<:FESpaceWithDirichletData},
  assem::MultiAssembler,
  terms::Vararg{<:AffineFETerm})

  V = MultiFESpace(testfesp)
  U = MultiFESpace(trialfesp)
  LinearFEOperator(V,U,assem,terms...)
end

function NonLinearFEOperator(
  res::Function,
  jac::Function,
  testfesp::Vector{<:FESpaceWithDirichletData},
  trialfesp::Vector{<:FESpaceWithDirichletData},
  assem::MultiAssembler,
  trian::Triangulation{Z},
  quad::CellQuadrature{Z}) where Z
  V = MultiFESpace(testfesp)
  U = MultiFESpace(trialfesp)
  NonLinearFEOperator(res,jac,V,U,assem,trian,quad)
end

function NonLinearFEOperator(
  testfesp::Vector{<:FESpaceWithDirichletData},
  trialfesp::Vector{<:FESpaceWithDirichletData},
  assem::MultiAssembler,
  terms::Vararg{<:FETerm})

  V = MultiFESpace(testfesp)
  U = MultiFESpace(trialfesp)
  NonLinearFEOperator(V,U,assem,terms...)
end

function NonLinearFEOperator(
  testfesp::Vector{<:FESpaceWithDirichletData},
  trialfesp::Vector{<:FESpaceWithDirichletData},
  terms::Vararg{<:FETerm})

  assem = SparseMatrixAssembler(testfesp,trialfesp)

  V = MultiFESpace(testfesp)
  U = MultiFESpace(trialfesp)
  NonLinearFEOperator(V,U,assem,terms...)
end

end # module MultiFEOperators


text:
using Documenter
using Comrade
using ComradeBase

using Literate

using ComradeAHMC
using ComradeGalactic
using GalacticBBO
using Glob
using Plots

# Make the examples using Literate
GENERATED = joinpath(@__DIR__, "src", "examples")
SOURCE_FILES = Glob.glob("*.jl", GENERATED)
foreach(fn -> Literate.markdown(fn, GENERATED), SOURCE_FILES)

format = Documenter.HTML(edit_link = "source",
                         prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = String[])

ENV["GKSwstype"] = "nul" # needed for the GR backend on headless servers
gr()

makedocs(;
    modules=[ComradeBase, Comrade],
    authors="Paul Tiede",
    repo="https://github.com/ptiede/Comrade.jl/blob/{commit}{path}#L{line}",
    sitename="Comrade.jl",
    format=format,
    pages=Any[
        "Home" => "index.md",
        "vlbi_imaging_problem.md",
        "Tutorials" => [
                       "examples/data.md",
                       "examples/black_hole_image.md",
                       "examples/nonanalytic.md"
                       ],
        "interface.md",
        "base_api.md",
        "api.md"
    ],
)

deploydocs(;
    repo="github.com/ptiede/Comrade.jl",
    push_preview=true,
    forcepush=true,
    devbranch = "main",
)


text:
module BaselineModels

using Random
using DataFrames
using StatsBase: mode

using ..Utils
using ..AbsTypes

import ..AbsTypes: fit!, transform!

export fit!,transform!
export Baseline, Identity

"""
    Baseline(
       default_args = Dict(
	       :name => "baseline",
          :output => :class,
          :strat => mode
       )
    )

Baseline model that returns the mode during classification.
"""
mutable struct Baseline <: Learner
   name::String
   model::Dict{Symbol,Any}

   function Baseline(args=Dict())
      default_args = Dict{Symbol,Any}(
         :name      => "baseline",
         :output    => :class,
         :strat     => mode
      )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
    Baseline(name::String,opt...)

 Helper function
"""
function Baseline(name::String;opt...)
   Baseline(Dict(:name=>name,Dict(pairs(opt))...))
end

"""
    fit!(bsl::Baseline,x::DataFrame,y::Vector)

Get the mode of the training data.
"""
function fit!(bsl::Baseline,x::DataFrame,y::Vector)
   bsl.model[:choice] = bsl.model[:strat](y)
end

"""
    transform!(bsl::Baseline,x::DataFrame)

Return the mode in classification.
"""
function transform!(bsl::Baseline,x::DataFrame)
  fill(bsl.model[:choice],size(x,1))
end

"""
    Identity(args=Dict())

Returns the input as output.
"""
mutable struct Identity <: Transformer
  name::String
  model::Dict{Symbol,Any}

  function Identity(args=Dict())
	 default_args = Dict{Symbol,Any}(
				:name => "identity"
			)
	 cargs = nested_dict_merge(default_args, args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
	 new(cargs[:name],cargs)
  end
end

"""
    Baseline(name::String,opt...)

 Helper function
"""
function Identity(name::String)
   Baseline(Dict(:name=>name))
end

"""
    fit!(idy::Identity,x::DataFrame,y::Vector)

Does nothing.
"""
function fit!(idy::Identity,x::DataFrame,y::Vector)
    nothing
end

"""
    transform!(idy::Identity,x::DataFrame)

Return the input as output.
"""
function transform!(idy::Identity,x::DataFrame)
    return x
end

end