text:
import OmegaCore
import OmegaCore: propose_and_logratio
export MH, mh, mh!

using InferenceBase: keepall, tonothing

"Metropolis Hastings Samplng Algorithm"
struct MHAlg end
const MH = MHAlg()

function propose_and_logratio end

"""
`mh(rng, ΩT, logdensity, f, n; proposal, state_init)`

Metropolis Hastings Sampler

Initialised at `state_init`, yields `n` samples `::ΩT` using Metropolis Hastings algorithm.

# Arguments
- `rng`: AbstractRng used to sample proposals in MH loop
- `ΩT`: Type of samples
- `n`: number of samples
- `logdensity`:  Density that `mh` will sample from. Function of `ω`.
- `f`: The random variable that we are taking samples from
- `proposal`: a proposal distributioion.   Has the form:
  `(ω_, log_p_q) = proposal(rng, ω)` where:
    `ω_` is the proposal
    `log_pqqp` is `log(g(p|q)/g(q|p))` the transition probability of moving from q to p
- `state_init`: point to initialise from
- `keep`: function from `i -> b::Bool` which says whether to keep the sample at ith iteration or not
Useful for defining burn in or thinning.
  Example `keep = i -> ignorefirstn(100)(i) & thin(m)(i)`
"""
function mh!(rng,
             logdensity,
             n,
             state_init,
             propose_and_logratio,
             samples;
             keep = keepall,
             prestore = identity,
             cb = tonothing)
  state = state_init
  plast = logdensity(state)
  qlast = 1.0
  accepted = 0
  i = 1
  s = 1
  while s <= n
    propstate, logtransitionp = propose_and_logratio(rng, state)
    p_ = logdensity(propstate)
    ratio = p_ - plast + logtransitionp
    if log(rand(rng)) < ratio
      state = propstate
      plast = p_
      accepted += 1
    end
    if keep(i)
      @inbounds samples[s] = prestore(state)
      s += 1
    end
    cb((i = i, p = p_, accepted = accepted))
    i += 1
  end
  samples
end

function mh(rng,
            logdensity,
            n,
            state_init::X,
            propose_and_logratio; keep = keepall,
            prestore = identity,
   cb = tonothing) where X
  mh!(rng,
      logdensity,
      n,
      state_init,
      propose_and_logratio,
      Vector{X}(undef, n);
      keep = keep,
      prestore = prestore, cb = cb)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# TODO: Move somewhere else
function loglikelihood(base::H3M, reduced::H3M, z::AbstractMatrix, lhmm::AbstractMatrix, N)
    sum(OneTo(length(base))) do i
        sum(OneTo(length(reduced))) do j
            z[i,j] * (log(reduced.ω[j] / (z[i,j] + eps())) + N*lhmm[i,j])
        end
    end
end

function cluster(base::H3M, reduced::H3M, τ::Integer, N::Integer; maxiter = 100, tol = 1e0)
    @argcheck maxiter >= 0

    z = zeros(length(base), length(reduced))
    history = EMHistory(false, 0, [])

    logtot = -Inf
    logtotp = logtot

    for it in 1:maxiter
        try
            reduced, lhmm, z = vhem_step(base, reduced, τ, N)
            logtotp = loglikelihood(base, reduced, z, lhmm, N)
        catch e
            @error e
            break
        end

        push!(history.logtots, logtotp)
        history.iterations += 1

        println("Iteration $it: logtot = $logtotp, diff = $(logtotp - logtot)")

        if logtotp - logtot < tol
            history.converged = true
            break
        end

        logtot = logtotp
    end

    history, z, reduced
end

mutable struct EMHistory
    converged::Bool
    iterations::Int
    logtots::Vector{Float64}
end

ext: jl


text:
interpolation_keywords = ["displacement","epsilon","pressure_interpolated","sigma","temperature_interpolated","\"displacement\"","\"epsilon\"","\"pressure_interpolated\"","\"sigma\"","\"temperature_interpolated\""]
uncompress_keywords = ["connectivity","offsets","bulk_node_ids","bulk_element_ids","Points","MaterialIDs","\"connectivity\"","\"offsets\"","\"bulk_node_ids\"","\"bulk_element_ids\"","\"Points\"","\"MaterialIDs\"","types","\"types\""]

mutable struct VTUKeyWords
	interpolation_keywords::Vector{String}
	uncompress_keywords::Vector{String}
end

const vtukeywords = VTUKeyWords(interpolation_keywords,uncompress_keywords)

function set_uncompress_keywords(uk::Vector{String})
	vtukeywords.uncompress_keywords = uk
	for k in copy(uk)
		push!(vtukeywords.uncompress_keywords, "\""*k*"\"")
	end
end

function set_interpolation_keywords(ik::Vector{String})
	vtukeywords.interpolation_keywords = ik
	for k in copy(ik)
		push!(vtukeywords.interpolation_keywords, "\""*k*"\"")
	end
end

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
using DataStructures: OrderedDict
using PrettyTables

# define `Subjects` as uppercase strings.

"""
    SubjectCounts

Type for counting the number of times each subject appears in a timetable.

Alias for `OrderedDict{String, Int}`.
"""
SubjectCounts = OrderedDict{String, Int}

"""
    Division

Type for representing a division of a timetable.

# Constructors
- `Division(grade::Int, section::Int)`
- `Division(grade::Int, section::Int, section_str::String)`

# Fields
- `grade::Int`: The grade of the division.
- `section::Int`: The section of the division.
- `section_str::String`: The section of the division as a string.
- `__str__::String`: The precomputed string representation of the division.
"""
struct Division
    grade::Int
    section::Int
    section_str::String
    __str__::String  # precomputed string representation of the division
end

"""
    Division(grade::Int, section::Int)

Return a division with the given grade and section.
`section_str` is inferred as the nth letter of the alphabet, where n=`section`.

Call the second constructor with `section_str` as the third argument.
"""
function Division(grade::Int, section::Int)
    if (section < 1) || (section > 26)
        error("`section::Int` not in bounds [1, 26]. Provide the `section_str::String` parameter.")
    end
    # uppercase letters start from char 65
    section_str = string(Char(64 + section))
    return Division(grade, section, section_str)
end

"""
    Division(grade::Int, section::Int, section_str::String)

Return a division with the given grade, section, section string and precomputed string representation.
"""
function Division(grade::Int, section::Int, section_str::String)
    # precompute __str__
    __str__ = "$(grade)-$(section_str)"
    return Division(grade, section, section_str, __str__)
end
Base.string(d::Division) = d.__str__
Base.show(io::IO, d::Division) = print(io, string(d))

"""
    Teacher

Type for representing a teacher in the timetable.

# Constructors
- `Teacher(name::String, id::String, subjects::Vector{String}, grades::Vector{Int})`

# Fields
- `name::String`: The name of the teacher.
- `id::String`: The unique id of the teacher.
- `subjects::Vector{String}`: The subjects taught by the teacher.
- `grades::Vector{Int}`: The grades taught by the teacher.
- `__str__::String`: The precomputed string representation of the teacher.
"""
struct Teacher
    name::String
    id::String
    subjects::Vector{String}
    grades::Vector{Int}
    __str__::String  # precomputed string representation of the teacher
end
"""
    Teacher(name::String, id::String, subjects::Vector{String}, grades::Vector{Int})

Return a teacher with the given name, id, subjects, grades and precomputed string representation.
"""
function Teacher(name::String, id::String, subjects::Vector{String}, grades::Vector{Int})
    # precompute __str__
    __str__ = id
    return Teacher(name, id, subjects, grades, __str__)
end
Base.string(t::Teacher) = t.__str__
Base.show(io::IO, t::Teacher) = print(io, string(t))

"""
    Period

Mutable type for representing a period in the timetable.
Empty period is initialized with `nothing` as both arguments.

# Fields
- `subject::Union{String,Nothing}`: The subject of the period.
- `teacher::Union{Teacher,Nothing}`: The teacher of the period.
"""
mutable struct Period
    subject::Union{String,Nothing}
    teacher::Union{Teacher,Nothing}
end
Base.string(p::Period) = "$(p.subject) ($(string(p.teacher)))"
Base.show(io::IO, p::Period) = print(io, string(p))

"""
    Timetable

Mutable type for representing a timetable.

# Constructors
- `Timetable(numperiods::Vector{Int}, subjectcounts::SubjectCounts, teachers::Vector{Teacher}, division::Division)`

# Fields
- `numperiods::Vector{Int}`: The number of periods in each row.
- `subjectcounts::SubjectCounts`: The number of times subjects must occur in the timetable, in total.
- `subjects::Vector{String}`: The subjects in the timetable.
- `teachers::Vector{Teacher}`: The teachers in the timetable.
- `subjectteachers::OrderedDict{String, Vector{Teacher}}`: The teachers teaching each subject.
- `teacher_strs::OrderedDict{String, Teacher}`: The teachers' string representations mapped to objects.
- `division::Division`: The division of the timetable.
- `data::Vector{Vector{Period}}`: The timetable data, as a vector of vectors of periods.

# Notes
- Only fields `numperiods`, `subjectcounts`, `teachers` and `division` are passed. The rest are inferred.
- See the constructor for the same.
"""
mutable struct Timetable
    numperiods::Vector{Int}
    subjectcounts::SubjectCounts
    subjects::Vector{String}
    teachers::Vector{Teacher}
    subjectteachers::OrderedDict{String,Vector{String}}
    teacher_strs::OrderedDict{String,Teacher}
    division::Division
    data::Vector{Vector{Period}}
end
"""
    Timetable(numperiods::Vector{Int}, subjectcounts::SubjectCounts, teachers::Vector{Teacher}, division::Division)

Return a timetable with the given number of periods, subject counts, teachers and division.

# Notes
- Fields `subjects`, `subjectteachers`, `teacher_strs` and `data` are inferred.
- `data` stores the actual vector of vectors of periods, initialized with (`nothing`, `nothing`) as arguments.
"""
function Timetable(numperiods::Vector{Int}, subjectcounts::SubjectCounts,
    teachers::Vector{Teacher}, division::Division)
    subjects = collect(String, keys(subjectcounts))
    subjectteachers = get_subjectteachers(subjects,teachers)
    teacher_strs = OrderedDict(zip(string.(teachers), teachers))
    # fill data
    data = [[] for _ in numperiods]
    for (i,n) in enumerate(numperiods)
        data[i] = [Period(nothing, nothing) for _ in 1:n]
    end
    return Timetable(numperiods, subjectcounts, subjects, teachers,
    subjectteachers, teacher_strs, division, data)
end

"""
    get_subjectteachers(subjects::Vector{String}, teachers::Vector{Teacher})

Return an OrderedDict mapping subjects to valid teachers, based on `subjects` and `teachers`.
"""
function get_subjectteachers(subjects::Vector{String},teachers::Vector{Teacher})
    # initialize dict with empty list mapped to each subject
    subjectteachers = OrderedDict{String,Vector{String}}(s=>[] for s in subjects)
    for t in teachers, s in subjects
        if s in t.subjects
            # push the teacher to the subject's list
            push!(subjectteachers[s], string(t))
        end
    end
    # filter out empty lists
    subjectteachers = filter(x->length(x)>0, subjectteachers)
    return subjectteachers
end
"""
    Base.string(tt::Timetable)::String

`Base.string` method for type Timetable.

Return a pretty table representation of the timetable.

# Notes
- Use the [`PrettyTables`](https://github.com/ronisbr/PrettyTables.jl) library to return a table-like string.
- Rows are numbered from 1 to `length(data)` and columns are numbered P1, P2 ... to longest row in data (`max(length.(data)...)`).
- Column 1 header is the division name.
- Call division and teacher's `Base.string` method for string representation.
"""
function Base.string(tt::Timetable)
    data = tt.data
    m, n = length(data), max(length.(data)...)
    header = vcat([string(tt.division)], ["P"*string(i) for i in 1:n])
    data_matrix = fill("", (m, n+1))
    for (i,row) in enumerate(data)
        matrix_row = string.(row)
        data_matrix[i,1] = string(i)
        data_matrix[i,2:lastindex(matrix_row)+1] = matrix_row
    end
    s = pretty_table(String, data_matrix, header=header, alignment=:L)
    return s
end
Base.show(io::IO, tt::Timetable) = print(io, string(tt))

"""
    Schedule

Mutable type for representing a schedule, i.e., a collection of timetables.

# Constructors
- `Schedule(timetables...)`

# Fields
- `data::OrderedDict{String,Timetable}`: The timetables in the schedule mapped to their division's string representation.
"""
mutable struct Schedule
    data::OrderedDict{String,Timetable}
end
"""
    Schedule(timetables...)

Return a schedule instance using the given timetables.

# Notes
`data` is an OrderedDict mapping division's string representation to the timetable.
"""
function Schedule(timetables...)
    # mapping div's string representation to timetable
    data = OrderedDict{String,Timetable}()
    for tt in timetables
        div = string(tt.division)
        data[div] = tt
    end
    return Schedule(data)
end

"""
    Base.string(s::Schedule)::String

`Base.string` method for type Schedule.

Return pretty table representation of the timetables.

# Notes
- Call the `Base.string` method for each timetable in the schedule and fill newlines in between.
"""
Base.string(s::Schedule) = join(string.(values(s.data)),"\n")
Base.show(io::IO, s::Schedule) = print(io, string(s))

"""
    modify!(s::Schedule, var::String, val::String)

Modify schedule `s` by replacing the value of the field `var` with `val`.

# Notes
- Modify the schedule `s` in place.
- Both `var` and `val` are strings.
- Call `modify!()` on the correct timetable in `s` with the same arguments.
"""
modify!(s::Schedule, var::String, val::String) = modify!(s.data[split(var, "_")[1]], var, val)

"""
    modify!(tt::Timetable, var::String, val::String)

Modify timetable `tt` by replacing the value of the field `var` with `val`.

# Notes
- Modify the timetable `tt` in place.
- Both `var` and `val` are strings.
- Row, period and type of variable are inferred by splitting `var` at `"_"`, and taking last 3 elements.
- If `type == "subject"` then the subject of the corresponding row and period is set to `val`.
- If `type == "teacher"` then the teacher of the corresponding row and period is set to `val` (after getting teacher object from teacher string using `tt.teacher_strs`).
"""
function modify!(tt::Timetable, var::String, val::String)
    # split var
    _, row, period, type = split(var, "_")
    row, period = parse(Int, row), parse(Int, period)

    # if else based on type
    if type == "subject"
        # change subject to val
        tt.data[row][period].subject = val
    elseif type == "teacher"
        # change teacher to Teacher instance matching val
        tt.data[row][period].teacher = tt.teacher_strs[val]
    else
        error("Invalid type $(type)")
    end
end


text:
facts("Trading metrics from transactions blotter") do
  context("Transaction blotter transformations") do
    blotter = TradingLogic.emptyblotter()

    # long than short with exit
    p1, p2, p3 = 108.0, 85.0, 74.0
    along, ashort = 50, 70
    tlong = Dates.DateTime(2015,1,1,3,8,17)
    tshort = Dates.DateTime(2015,1,1,5,2,4)
    texit = Dates.DateTime(2015,1,3,15,12,48)

    # enter transactions not in chronological order
    blotter[tshort] = (-(along + ashort), p2)
    blotter[tlong] = (along, p1)
    blotter[texit] = (ashort, p3)
    TradingLogic.printblotter(STDOUT, blotter)

    # sorted vectors
    @fact TradingLogic.vtblotter(blotter) --> [tlong, tshort, texit]
    va, vp = TradingLogic.vapblotter(blotter)
    @fact va --> [along, -(along + ashort), ashort]
    @fact vp --> roughly([p1, p2, p3])

    # print to file
    dtfmt = "yyyy-mm-dd HH:MM:SS" # not default
    fnm = "writetest.csv"
    TradingLogic.writeblotter(fnm, blotter, dtformat = dtfmt)
    rcsv = readcsv(fnm)
    run(`rm $fnm`)
    @fact size(rcsv) --> (4,3)
    @fact rcsv[1,2] --> "Amount"
    @fact rcsv[2,1] --> Dates.format(tlong, dtfmt)
    @fact rcsv[4,1] --> Dates.format(texit, dtfmt)
    ##@fact int(rcsv[3,2]) --> -(along + ashort)
    @fact round(Int64, rcsv[3,2]) --> -(along + ashort)
    @fact float(rcsv[4,3]) --> roughly(p3)
  end
  context("Long enter with partial fill") do
    blotter = TradingLogic.emptyblotter()
    # fill prices and exit price
    p1, p2, pexit = 120.0, 135.0, 150.0
    # amounts
    a1, a2 = 20, 30

    blotter[DateTime(2015,1,1,0)] = (a1, p1)
    blotter[DateTime(2015,1,1,1)] = (a2, p2)
    blotter[DateTime(2015,1,2)] = (-(a1 + a2), pexit)
    pnlfin = a1*(p2 - p1) + (a1 + a2)*(pexit - p2)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    ntrans = length(blotter)
    @fact length(vt) --> ntrans

    @fact length(perfm[:PnL]) --> ntrans

    @fact perfm[:Qty] --> [a1, a2, -(a1 + a2)]
    @fact perfm[:FillPrice] --> roughly([p1, p2, pexit])

    @fact perfm[:PnL] --> roughly([0.0, a1*(p2 - p1), pnlfin])
    @fact TradingLogic.tradepnlfinal(blotter) --> roughly(pnlfin)
  end
  context("Short enter and exit") do
    blotter = TradingLogic.emptyblotter()
    # short enter and exit price
    pent, pexit = 120.0, 115.2
    # short position held
    a = 100

    blotter[DateTime(2015,1,1,0)] = (-a, pent)
    blotter[DateTime(2015,1,1,1)] = (a, pexit)
    pnlfin = -a*(pexit - pent)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    ntrans = length(blotter)
    @fact length(vt) --> ntrans

    @fact length(perfm[:PnL]) --> ntrans

    @fact perfm[:Qty] --> [-a, a]
    @fact perfm[:FillPrice] --> roughly([pent, pexit])

    @fact perfm[:PnL] --> roughly([0.0, pnlfin])
    @fact TradingLogic.tradepnlfinal(blotter) --> roughly(pnlfin)
  end
  context("Long then short") do
    blotter = TradingLogic.emptyblotter()
    # prices and amounts
    p1, p2, p3 = 120.0, 85.0, 70.0
    along, ashort = 50, 70

    losslong = along*(p2 - p1)
    pnlfin = losslong + ashort*(p2 - p3)

    blotter[DateTime(2015,1,1,0)] = (along, p1)
    blotter[DateTime(2015,1,1,1)] = (-(along + ashort), p2)

    # not yet final PnL with incomplete blotter
    @fact abs(TradingLogic.tradepnlfinal(blotter) - pnlfin) > 10.0 --> true
    # test tradepnlfinal method with pnow: final PnL
    #  doing this before adding the final (exit) transaction to blotter
    @fact TradingLogic.tradepnlfinal(blotter, p3) --> roughly(pnlfin)

    # final transaction (e.g. stop backtest point)
    blotter[DateTime(2015,1,2)] = (0, p3)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    #println(vt, perfm)
    ntrans = length(blotter)
    @fact length(vt) --> ntrans
    @fact length(perfm[:PnL]) --> ntrans

    @fact perfm[:Qty] --> [along, -(along + ashort), 0]
    @fact perfm[:FillPrice] --> roughly([p1, p2, p3])

    @fact perfm[:PnL] --> roughly([0.0, losslong, pnlfin])
    @fact TradingLogic.tradepnlfinal(blotter) --> roughly(pnlfin)

    # return-based drawdown
    ddmax = abs(losslong)
    @fact perfm[:DDown] --> roughly([0.0, -ddmax, -abs(pnlfin)])
  end
end

facts("Running performance metrics") do
  context("Current maximum PnL and drawdown over trading history") do
    pnl_dd_max = TradingLogic.tradeperffold((10.0, 0.0), (true, 15.0))
    @fact pnl_dd_max[1] --> roughly(15.0)
    @fact pnl_dd_max[2] --> roughly(0.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 20.0), (true, 15.0))
    @fact pnl_dd_max[1] --> roughly(15.0)
    @fact pnl_dd_max[2] --> roughly(20.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 0.0), (true, 5.0))
    @fact pnl_dd_max[1] --> roughly(10.0)
    @fact pnl_dd_max[2] --> roughly(5.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 50.0), (true, 5.0))
    @fact pnl_dd_max[1] --> roughly(10.0)
    @fact pnl_dd_max[2] --> roughly(50.0)
  end
end

facts("Trades and corresponding metrics from transactions blotter") do
  blotter = TradingLogic.emptyblotter()
  dthour(h::Int64) = Dates.DateTime(2015,6,1,h)
  blotter[dthour(1)] = (25, 85.0)
  blotter[dthour(2)] = (25, 95.0)
  blotter[dthour(3)] = (-60, 100.0)
  blotter[dthour(4)] = (20, 80.0)
  blotter[dthour(5)] = (-10, 75.0)
  blotter[dthour(6)] = (10, 85.0)
  blotter[dthour(7)] = (-10, 70.0)
  blotter[dthour(8)] = (-50, 65.0)
  blotter[dthour(11)] = (50, 25.0)
  # open trade
  blotter[dthour(12)] = (20, 35.0)
  #TradingLogic.printblotter(STDOUT, blotter)
  @fact TradingLogic.tradepnlfinal(blotter) --> roughly(2500.0)

  context("Trades PnL vector with selected metrics from blotter") do
    vtrpnl, ntrpos, avwin, ntrneg, avloss = TradingLogic.vtradespnl(blotter)
    @fact length(vtrpnl) --> 5
    @fact sum(vtrpnl) --> roughly(2500.0)
    @fact ntrpos --> 3
    @fact avwin --> roughly(900.0)
    @fact ntrneg --> 2
    @fact avloss --> roughly(100.0)
  end

  context("Trades-based performance metrics: PROM") do
    marg = 5e3

    # no transactions
    @fact TradingLogic.perf_prom(
      TradingLogic.emptyblotter()) --> roughly(0.0)

    # from blotter, using all completed trades
    prom0 = ((3.0 - sqrt(3.0))*900.0 - (2.0 + sqrt(2.0))*100.0)/marg
    @fact TradingLogic.perf_prom(blotter, marg=marg) --> roughly(prom0)

    # removing specific number of best wins
    prom1 = ((2.0 - sqrt(2.0))*350.0 - (2.0 + sqrt(2.0))*100.0)/marg
    @fact TradingLogic.perf_prom(blotter, marg=marg,
                                 nbest_remove = 1) --> roughly(prom1)

    # from trades PnL vector
    vtrpnl, ntrpos, avwin, ntrneg, avloss = TradingLogic.vtradespnl(blotter)
    @fact TradingLogic.perf_prom(vtrpnl, marg=marg,
                                 nbest_remove = 1) --> roughly(prom1)
    prom2 = -(2.0 + sqrt(2.0))*100.0/marg
    @fact TradingLogic.perf_prom(vtrpnl, marg=marg,
                                 nbest_remove = 2) --> roughly(prom2)
    # if asking to remove more best wins than the total number of wins
    prom = TradingLogic.perf_prom(vtrpnl, marg=marg, nbest_remove = 5)
    # -Inf
    @fact prom < 0.0 --> true
    @fact isfinite(prom) --> false
    prom = TradingLogic.perf_prom(vtrpnl, marg=marg, nbest_remove = 6)
    # -Inf
    @fact prom < 0.0 --> true
    @fact isfinite(prom) --> false
  end

  context("Trades-based performance metrics: PROR") do
    pror0 = (3.0 - sqrt(3.0))*900.0 / ((2.0 + sqrt(2.0))*100.0)
    @fact TradingLogic.perf_prom(blotter, pror=true) --> roughly(pror0)
    pror1 = (2.0 - sqrt(2.0))*350.0 / ((2.0 + sqrt(2.0))*100.0)
    @fact TradingLogic.perf_prom(blotter, pror=true,
                                 nbest_remove = 1) --> roughly(pror1)

    # PROR without profits
    @fact TradingLogic.perf_prom(TradingLogic.emptyblotter(),
                                 pror=true) --> roughly(0.0)
    @fact TradingLogic.perf_prom(Array(Float64,0), pror=true) --> roughly(0.0)
    @fact TradingLogic.perf_prom([-2.0], pror=true) --> roughly(0.0)
    pror = TradingLogic.perf_prom([-2.0], pror=true, nbest_remove = 1)
    # -Inf
    @fact pror < 0.0 --> true
    @fact isfinite(pror) --> false

    # PROR without profits or losses: 1 - sqrt(1)
    @fact TradingLogic.perf_prom([5.0, -2.0], pror=true) --> roughly(0.0)
    @fact TradingLogic.perf_prom([5.0, -2.0], pror=true,
                                 nbest_remove = 1) --> roughly(0.0)
    # PROR without losses
    pror = TradingLogic.perf_prom([5.0, 2.0], pror=true, nbest_remove = 0)
    # +Inf
    @fact pror < 0.0 --> false
    @fact isfinite(pror) --> false
    @fact TradingLogic.perf_prom([5.0, 2.0], pror=true,
                                 nbest_remove = 1) --> roughly(0.0)
    pror = TradingLogic.perf_prom([5.0, 2.0], pror=true, nbest_remove = 2)
    # -Inf
    @fact pror < 0.0 --> true
    @fact isfinite(pror) --> false
  end
end


text:
using BQCESubroutine
using BQCESubroutine.Utilities
using LinearAlgebra
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using Test
using BenchmarkTools
using YaoArrayRegister
using LoopVectorization
using PaddedMatrices
using MatchCore
using MLStyle
using BQCESubroutine: BitContext, BitRoutine
using BQCESubroutine: diag_kernel, perm_kernel, generic_kernel, codegen_basic_broutine
using BQCESubroutine: codegen_broutine2x2_multi, codegen_broutine, hoist_entries!

BQCESubroutine.ENABLE_THREADS[]
BQCESubroutine.disable_threads()
U = rand(ComplexF64, 4, 4);
@benchmark broutine!(st, $U, Locations((3, 24))) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark instruct!(st, $U, (3, 24)) setup=(st=rand(ComplexF64, 1 << 25))

BQCESubroutine.disable_threads()
@benchmark broutine!(st, Val(:Z), Locations((3, 24))) setup=(st=rand(ComplexF64, 1 << 25))

BQCESubroutine.enable_threads()
@benchmark broutine!(st, Val(:Y), Locations((3, 24))) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark instruct!(st, Val(:Y), (3, 24)) setup=(st=rand(ComplexF64, 1 << 25))


text:
using Test, Pkg

@testset "Atmos" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false
    for submodule in ["Parameterizations", "TemperatureProfiles"]
        if all_tests || "$submodule" in ARGS || "Atmos" in ARGS
            include_test(submodule)
        end
    end

end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(y), 0.8, shuffle=true, rng=RANDOM_SEED);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# TODO: add more target distributions and make them iteratable

# Dimension of testing distribution
const D = 5
# Tolerance ratio
const TRATIO = Int == Int64 ? 1 : 2
# Deterministic tolerance
const DETATOL = 1e-3 * D * TRATIO
# Random tolerance
const RNDATOL = 5e-2 * D * TRATIO

# Hand-coded multivaite Gaussain

const gaussian_m = zeros(D)
const gaussian_s = ones(D)

function ℓπ(m, s, x::AbstractVecOrMat{T}) where {T}
    diff = x .- m
    v = s.^2
    return -(log(2 * T(pi)) .+ log.(v) .+ diff .* diff ./ v) / 2
end

function ℓπ(θ::AbstractVector)
    return sum(ℓπ(gaussian_m, gaussian_s, θ))
end

function ℓπ(θ::AbstractMatrix)
    return dropdims(sum(ℓπ(gaussian_m, gaussian_s, θ); dims=1); dims=1)
end

function ∂ℓπ∂θ(m, s, x::AbstractVecOrMat{T}) where {T}
    diff = x .- m
    v = s.^2
    v = -(log(2 * T(pi)) .+ log.(v) .+ diff .* diff ./ v) / 2
    g = -diff
    return v, g
end

function ∂ℓπ∂θ(θ::AbstractVector)
    v, g = ∂ℓπ∂θ(gaussian_m, gaussian_s, θ)
    return sum(v), g
end

function ∂ℓπ∂θ(θ::AbstractMatrix)
    v, g = ∂ℓπ∂θ(gaussian_m, gaussian_s, θ)
    return dropdims(sum(v; dims=1); dims=1), g
end

# For the Turing model
# @model gdemo() = begin
#     s ~ InverseGamma(2, 3)
#     m ~ Normal(0, sqrt(s))
#     1.5 ~ Normal(m, sqrt(s))
#     2.0 ~ Normal(m, sqrt(s))
#     return s, m
# end

using Distributions: logpdf, InverseGamma, Normal
using Bijectors: invlink, logpdf_with_trans

function invlink_gdemo(θ)
    s = invlink(InverseGamma(2, 3), θ[1])
    m = θ[2]
    return [s, m]
end

function ℓπ_gdemo(θ)
    s, m = invlink_gdemo(θ)
    logprior = logpdf_with_trans(InverseGamma(2, 3), s, true) + logpdf(Normal(0, sqrt(s)), m)
    loglikelihood = logpdf(Normal(m, sqrt(s)), 1.5) + logpdf(Normal(m, sqrt(s)), 2.0)
    return logprior + loglikelihood
end


text:
"""

Update the data with the data accumulated in a frame

"""
function update_counters_frame!(R::Result,
                                rdf_count_random_frame::Vector{Float64},
                                md_count_random_frame::Vector{Float64},
                                volume_frame::Volume,
                                solute::Selection, solvent::Selection,
                                n_solvent_in_bulk::Float64)

  # Volume of each bin shell and of the solute domain
  @. volume_frame.shell =
        volume_frame.total * (rdf_count_random_frame/(R.options.n_random_samples*solvent.nmols))
  volume_frame.domain = sum(volume_frame.shell)

  # Volume and density of the bulk region in this frame
  if ! R.options.usecutoff
    volume_frame.bulk = volume_frame.total - volume_frame.domain
  else
    ibulk = setbin(R.dbulk+0.5*R.options.binstep,R.options.binstep)
    for i in ibulk:R.nbins
      volume_frame.bulk = volume_frame.bulk + volume_frame.shell[i]
    end
  end

  bulk_density_at_frame = n_solvent_in_bulk / volume_frame.bulk

  # The number of "random" molecules generated must be that to fill the complete
  # volume with the density of the bulk. Now that we know the density of the bulk
  # and the volume of the frame, we can adjust that random count
  fillup_factor = bulk_density_at_frame * volume_frame.total / solvent.nmols
  @. R.rdf_count_random =
        R.rdf_count_random + rdf_count_random_frame*fillup_factor
  @. R.md_count_random =
        R.md_count_random + md_count_random_frame*fillup_factor

  # Update final counters, which have to be normalized by the number of frames
  # at the end
  @. R.volume.shell = R.volume.shell + volume_frame.shell
  R.volume.bulk = R.volume.bulk + volume_frame.bulk
  R.volume.domain = R.volume.domain + volume_frame.domain
  R.density.solvent_bulk = R.density.solvent_bulk + bulk_density_at_frame

  return nothing
end


text:
# taken from https://github.com/JuliaLang/julia/pull/31376
module MD
disable_nonforced() = (Symbol("llvm.loop.disable_nonforced"),)
interleave(n) = (Symbol("llvm.loop.interleave.count"), convert(Int, n))
vectorize_enable(flag) = (Symbol("llvm.loop.vectorize.enable"), convert(Bool, flag))
vectorize_width(n) = (Symbol("llvm.loop.vectorize.width"), convert(Int, n))
# ‘llvm.loop.vectorize.followup_vectorized’
# ‘llvm.loop.vectorize.followup_epilogue’
# ‘llvm.loop.vectorize.followup_all’
unroll_count(n) = (Symbol("llvm.loop.unroll.count"), convert(Int, n))
unroll_disable() = (Symbol("llvm.loop.unroll.disable"),)
unroll_enable() = (Symbol("llvm.loop.unroll.enable"),)
unroll_full() = (Symbol("llvm.loop.unroll.full"),)
# ‘llvm.loop.unroll.followup’
# ‘llvm.loop.unroll.followup_remainder’
jam_count(n) = (Symbol("llvm.loop.unroll_and_jam.count"), convert(Int, n))
jam_disable() = (Symbol("llvm.loop.unroll_and_jam.disable"),)
jam_enable() = (Symbol("llvm.loop.unroll_and_jam.enable"),)
# ‘llvm.loop.unroll_and_jam.followup_outer’
# ‘llvm.loop.unroll_and_jam.followup_inner’
# ‘llvm.loop.unroll_and_jam.followup_remainder_outer’
# ‘llvm.loop.unroll_and_jam.followup_remainder_inner’
# ‘llvm.loop.unroll_and_jam.followup_all’
# ‘llvm.loop.licm_versioning.disable’
# ‘llvm.loop.distribute.enable’
# ‘llvm.loop.distribute.followup_coincident’
# ‘llvm.loop.distribute.followup_sequential’
# ‘llvm.loop.distribute.followup_fallback’
# ‘llvm.loop.distribute.followup_all’
end

function loopinfo(name, expr, nodes...)
    if expr.head != :for
        error("Syntax error: pragma $name needs a for loop")
    end
    push!(expr.args[2].args, Expr(:loopinfo, nodes...))
    return expr
end

"""
   @unroll expr
Takes a for loop as `expr` and informs the LLVM unroller to fully unroll it, if
it is safe to do so and the loop count is known.
"""
macro unroll(expr)
    expr = loopinfo("@unroll", expr, MD.unroll_full())
    return esc(expr)
end

"""
    @unroll N expr
Takes a for loop as `expr` and informs the LLVM unroller to unroll it `N` times,
if it is safe to do so.
"""
macro unroll(N, expr)
    if !(N isa Integer)
        error("Syntax error: `@unroll N expr` needs a constant integer N")
    end
    expr = loopinfo("@unroll", expr, MD.unroll_count(N))
    return esc(expr)
end


text:
module AttachDirichletTests

using Gridap.Arrays
using Gridap.CellData

ncells = 10
ndofs = 3
cellmat = [ rand(ndofs,ndofs) for cell in 1:ncells]
cellvec = [ rand(ndofs) for cell in 1:ncells]
cellmatvec = pair_arrays(cellmat,cellvec)
cellvals = [ rand(ndofs) for cell in 1:ncells]

cellmatvec_with_dbc = attach_dirichlet(cellmatvec,cellvals)
r = collect(cellmatvec_with_dbc)
test_array(cellmatvec_with_dbc,r)

r = map((mat,vec,vals)-> (mat,vec-mat*vals),cellmat,cellvec,cellvals)
a,b = unpair_arrays(cellmatvec_with_dbc)
ra,rb = unpair_arrays(r)
test_array(a,ra,≈)
test_array(b,rb,≈)

cellmatvec_with_dbc = attach_dirichlet(cellmat,cellvals)
r = collect(cellmatvec_with_dbc)
test_array(cellmatvec_with_dbc,r)

r = map((mat,vals)-> (mat,-mat*vals),cellmat,cellvals)
a,b = unpair_arrays(cellmatvec_with_dbc)
ra,rb = unpair_arrays(r)
test_array(a,ra,≈)
test_array(b,rb,≈)

#a = cellmatvec_with_dbc
#cache = array_cache(a)
#using BenchmarkTools
#@btime getindex!($cache,$a,2)

end # module


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# Generates a realization of an cosine sequence

"""
	gen_cos_sequence(n = 100, f = 1.0, fs = 1.0, phase = 0.0)

Generates a realization of an cosine sequence

# Arguments
- `n`: Length of sequence
- `f`: Frequency
- `fs`: Sampling frequency
- `phase`: Phase shift

# Outputs
- `x`: vector of cosine process realization
- `t`: Time vector
"""
function gen_cos_sequence(n = 100, f = 1.0, fs = 1.0, phase = 0.0)

    n = round(abs(n[1])) # Samples
    n = (1 <= n) ? n : 1

    f = abs(f[1]) # Frequency
    fs = abs(fs[1]) # Sampling Frequency

    phase = phase[1] # Phase shift

    t = collect(0.0:(1.0 / fs):((n - 1) / fs)) # Time vector

    x = cos.(f * t .+ phase) # Cosines: time signature

    return x, t
end


text:
# Automatically generated using Clang.jl

const TIFF_VERSION_CLASSIC = 42
const TIFF_VERSION_BIG = 43
const TIFF_BIGENDIAN = 0x4d4d
const TIFF_LITTLEENDIAN = 0x4949
const MDI_LITTLEENDIAN = 0x5045
const MDI_BIGENDIAN = 0x4550
const TIFFTAG_SUBFILETYPE = 254
const FILETYPE_REDUCEDIMAGE = 0x01
const FILETYPE_PAGE = 0x02
const FILETYPE_MASK = 0x04
const TIFFTAG_OSUBFILETYPE = 255
const OFILETYPE_IMAGE = 1
const OFILETYPE_REDUCEDIMAGE = 2
const OFILETYPE_PAGE = 3
const TIFFTAG_IMAGEWIDTH = 256
const TIFFTAG_IMAGELENGTH = 257
const TIFFTAG_BITSPERSAMPLE = 258
const TIFFTAG_COMPRESSION = 259
const COMPRESSION_NONE = 1
const COMPRESSION_CCITTRLE = 2
const COMPRESSION_CCITTFAX3 = 3
const COMPRESSION_CCITT_T4 = 3
const COMPRESSION_CCITTFAX4 = 4
const COMPRESSION_CCITT_T6 = 4
const COMPRESSION_LZW = 5
const COMPRESSION_OJPEG = 6
const COMPRESSION_JPEG = 7
const COMPRESSION_T85 = 9
const COMPRESSION_T43 = 10
const COMPRESSION_NEXT = 32766
const COMPRESSION_CCITTRLEW = 32771
const COMPRESSION_PACKBITS = 32773
const COMPRESSION_THUNDERSCAN = 32809
const COMPRESSION_IT8CTPAD = 32895
const COMPRESSION_IT8LW = 32896
const COMPRESSION_IT8MP = 32897
const COMPRESSION_IT8BL = 32898
const COMPRESSION_PIXARFILM = 32908
const COMPRESSION_PIXARLOG = 32909
const COMPRESSION_DEFLATE = 32946
const COMPRESSION_ADOBE_DEFLATE = 8
const COMPRESSION_DCS = 32947
const COMPRESSION_JBIG = 34661
const COMPRESSION_SGILOG = 34676
const COMPRESSION_SGILOG24 = 34677
const COMPRESSION_JP2000 = 34712
const COMPRESSION_LERC = 34887
const COMPRESSION_LZMA = 34925
const COMPRESSION_ZSTD = 50000
const COMPRESSION_WEBP = 50001
const TIFFTAG_PHOTOMETRIC = 262
const PHOTOMETRIC_MINISWHITE = 0
const PHOTOMETRIC_MINISBLACK = 1
const PHOTOMETRIC_RGB = 2
const PHOTOMETRIC_PALETTE = 3
const PHOTOMETRIC_MASK = 4
const PHOTOMETRIC_SEPARATED = 5
const PHOTOMETRIC_YCBCR = 6
const PHOTOMETRIC_CIELAB = 8
const PHOTOMETRIC_ICCLAB = 9
const PHOTOMETRIC_ITULAB = 10
const PHOTOMETRIC_CFA = 32803
const PHOTOMETRIC_LOGL = 32844
const PHOTOMETRIC_LOGLUV = 32845
const TIFFTAG_THRESHHOLDING = 263
const THRESHHOLD_BILEVEL = 1
const THRESHHOLD_HALFTONE = 2
const THRESHHOLD_ERRORDIFFUSE = 3
const TIFFTAG_CELLWIDTH = 264
const TIFFTAG_CELLLENGTH = 265
const TIFFTAG_FILLORDER = 266
const FILLORDER_MSB2LSB = 1
const FILLORDER_LSB2MSB = 2
const TIFFTAG_DOCUMENTNAME = 269
const TIFFTAG_IMAGEDESCRIPTION = 270
const TIFFTAG_MAKE = 271
const TIFFTAG_MODEL = 272
const TIFFTAG_STRIPOFFSETS = 273
const TIFFTAG_ORIENTATION = 274
const ORIENTATION_TOPLEFT = 1
const ORIENTATION_TOPRIGHT = 2
const ORIENTATION_BOTRIGHT = 3
const ORIENTATION_BOTLEFT = 4
const ORIENTATION_LEFTTOP = 5
const ORIENTATION_RIGHTTOP = 6
const ORIENTATION_RIGHTBOT = 7
const ORIENTATION_LEFTBOT = 8
const TIFFTAG_SAMPLESPERPIXEL = 277
const TIFFTAG_ROWSPERSTRIP = 278
const TIFFTAG_STRIPBYTECOUNTS = 279
const TIFFTAG_MINSAMPLEVALUE = 280
const TIFFTAG_MAXSAMPLEVALUE = 281
const TIFFTAG_XRESOLUTION = 282
const TIFFTAG_YRESOLUTION = 283
const TIFFTAG_PLANARCONFIG = 284
const PLANARCONFIG_CONTIG = 1
const PLANARCONFIG_SEPARATE = 2
const TIFFTAG_PAGENAME = 285
const TIFFTAG_XPOSITION = 286
const TIFFTAG_YPOSITION = 287
const TIFFTAG_FREEOFFSETS = 288
const TIFFTAG_FREEBYTECOUNTS = 289
const TIFFTAG_GRAYRESPONSEUNIT = 290
const GRAYRESPONSEUNIT_10S = 1
const GRAYRESPONSEUNIT_100S = 2
const GRAYRESPONSEUNIT_1000S = 3
const GRAYRESPONSEUNIT_10000S = 4
const GRAYRESPONSEUNIT_100000S = 5
const TIFFTAG_GRAYRESPONSECURVE = 291
const TIFFTAG_GROUP3OPTIONS = 292
const TIFFTAG_T4OPTIONS = 292
const GROUP3OPT_2DENCODING = 0x01
const GROUP3OPT_UNCOMPRESSED = 0x02
const GROUP3OPT_FILLBITS = 0x04
const TIFFTAG_GROUP4OPTIONS = 293
const TIFFTAG_T6OPTIONS = 293
const GROUP4OPT_UNCOMPRESSED = 0x02
const TIFFTAG_RESOLUTIONUNIT = 296
const RESUNIT_NONE = 1
const RESUNIT_INCH = 2
const RESUNIT_CENTIMETER = 3
const TIFFTAG_PAGENUMBER = 297
const TIFFTAG_COLORRESPONSEUNIT = 300
const COLORRESPONSEUNIT_10S = 1
const COLORRESPONSEUNIT_100S = 2
const COLORRESPONSEUNIT_1000S = 3
const COLORRESPONSEUNIT_10000S = 4
const COLORRESPONSEUNIT_100000S = 5
const TIFFTAG_TRANSFERFUNCTION = 301
const TIFFTAG_SOFTWARE = 305
const TIFFTAG_DATETIME = 306
const TIFFTAG_ARTIST = 315
const TIFFTAG_HOSTCOMPUTER = 316
const TIFFTAG_PREDICTOR = 317
const PREDICTOR_NONE = 1
const PREDICTOR_HORIZONTAL = 2
const PREDICTOR_FLOATINGPOINT = 3
const TIFFTAG_WHITEPOINT = 318
const TIFFTAG_PRIMARYCHROMATICITIES = 319
const TIFFTAG_COLORMAP = 320
const TIFFTAG_HALFTONEHINTS = 321
const TIFFTAG_TILEWIDTH = 322
const TIFFTAG_TILELENGTH = 323
const TIFFTAG_TILEOFFSETS = 324
const TIFFTAG_TILEBYTECOUNTS = 325
const TIFFTAG_BADFAXLINES = 326
const TIFFTAG_CLEANFAXDATA = 327
const CLEANFAXDATA_CLEAN = 0
const CLEANFAXDATA_REGENERATED = 1
const CLEANFAXDATA_UNCLEAN = 2
const TIFFTAG_CONSECUTIVEBADFAXLINES = 328
const TIFFTAG_SUBIFD = 330
const TIFFTAG_INKSET = 332
const INKSET_CMYK = 1
const INKSET_MULTIINK = 2
const TIFFTAG_INKNAMES = 333
const TIFFTAG_NUMBEROFINKS = 334
const TIFFTAG_DOTRANGE = 336
const TIFFTAG_TARGETPRINTER = 337
const TIFFTAG_EXTRASAMPLES = 338
const EXTRASAMPLE_UNSPECIFIED = 0
const EXTRASAMPLE_ASSOCALPHA = 1
const EXTRASAMPLE_UNASSALPHA = 2
const TIFFTAG_SAMPLEFORMAT = 339
const SAMPLEFORMAT_UINT = 1
const SAMPLEFORMAT_INT = 2
const SAMPLEFORMAT_IEEEFP = 3
const SAMPLEFORMAT_VOID = 4
const SAMPLEFORMAT_COMPLEXINT = 5
const SAMPLEFORMAT_COMPLEXIEEEFP = 6
const TIFFTAG_SMINSAMPLEVALUE = 340
const TIFFTAG_SMAXSAMPLEVALUE = 341
const TIFFTAG_CLIPPATH = 343
const TIFFTAG_XCLIPPATHUNITS = 344
const TIFFTAG_YCLIPPATHUNITS = 345
const TIFFTAG_INDEXED = 346
const TIFFTAG_JPEGTABLES = 347
const TIFFTAG_OPIPROXY = 351
const TIFFTAG_GLOBALPARAMETERSIFD = 400
const TIFFTAG_PROFILETYPE = 401
const PROFILETYPE_UNSPECIFIED = 0
const PROFILETYPE_G3_FAX = 1
const TIFFTAG_FAXPROFILE = 402
const FAXPROFILE_S = 1
const FAXPROFILE_F = 2
const FAXPROFILE_J = 3
const FAXPROFILE_C = 4
const FAXPROFILE_L = 5
const FAXPROFILE_M = 6
const TIFFTAG_CODINGMETHODS = 403
const CODINGMETHODS_T4_1D = 1 << 1
const CODINGMETHODS_T4_2D = 1 << 2
const CODINGMETHODS_T6 = 1 << 3
const CODINGMETHODS_T85 = 1 << 4
const CODINGMETHODS_T42 = 1 << 5
const CODINGMETHODS_T43 = 1 << 6
const TIFFTAG_VERSIONYEAR = 404
const TIFFTAG_MODENUMBER = 405
const TIFFTAG_DECODE = 433
const TIFFTAG_IMAGEBASECOLOR = 434
const TIFFTAG_T82OPTIONS = 435
const TIFFTAG_JPEGPROC = 512
const JPEGPROC_BASELINE = 1
const JPEGPROC_LOSSLESS = 14
const TIFFTAG_JPEGIFOFFSET = 513
const TIFFTAG_JPEGIFBYTECOUNT = 514
const TIFFTAG_JPEGRESTARTINTERVAL = 515
const TIFFTAG_JPEGLOSSLESSPREDICTORS = 517
const TIFFTAG_JPEGPOINTTRANSFORM = 518
const TIFFTAG_JPEGQTABLES = 519
const TIFFTAG_JPEGDCTABLES = 520
const TIFFTAG_JPEGACTABLES = 521
const TIFFTAG_YCBCRCOEFFICIENTS = 529
const TIFFTAG_YCBCRSUBSAMPLING = 530
const TIFFTAG_YCBCRPOSITIONING = 531
const YCBCRPOSITION_CENTERED = 1
const YCBCRPOSITION_COSITED = 2
const TIFFTAG_REFERENCEBLACKWHITE = 532
const TIFFTAG_STRIPROWCOUNTS = 559
const TIFFTAG_XMLPACKET = 700
const TIFFTAG_OPIIMAGEID = 32781
const TIFFTAG_REFPTS = 32953
const TIFFTAG_REGIONTACKPOINT = 32954
const TIFFTAG_REGIONWARPCORNERS = 32955
const TIFFTAG_REGIONAFFINE = 32956
const TIFFTAG_MATTEING = 32995
const TIFFTAG_DATATYPE = 32996
const TIFFTAG_IMAGEDEPTH = 32997
const TIFFTAG_TILEDEPTH = 32998
const TIFFTAG_PIXAR_IMAGEFULLWIDTH = 33300
const TIFFTAG_PIXAR_IMAGEFULLLENGTH = 33301
const TIFFTAG_PIXAR_TEXTUREFORMAT = 33302
const TIFFTAG_PIXAR_WRAPMODES = 33303
const TIFFTAG_PIXAR_FOVCOT = 33304
const TIFFTAG_PIXAR_MATRIX_WORLDTOSCREEN = 33305
const TIFFTAG_PIXAR_MATRIX_WORLDTOCAMERA = 33306
const TIFFTAG_WRITERSERIALNUMBER = 33405
const TIFFTAG_CFAREPEATPATTERNDIM = 33421
const TIFFTAG_CFAPATTERN = 33422
const TIFFTAG_COPYRIGHT = 33432
const TIFFTAG_RICHTIFFIPTC = 33723
const TIFFTAG_IT8SITE = 34016
const TIFFTAG_IT8COLORSEQUENCE = 34017
const TIFFTAG_IT8HEADER = 34018
const TIFFTAG_IT8RASTERPADDING = 34019
const TIFFTAG_IT8BITSPERRUNLENGTH = 34020
const TIFFTAG_IT8BITSPEREXTENDEDRUNLENGTH = 34021
const TIFFTAG_IT8COLORTABLE = 34022
const TIFFTAG_IT8IMAGECOLORINDICATOR = 34023
const TIFFTAG_IT8BKGCOLORINDICATOR = 34024
const TIFFTAG_IT8IMAGECOLORVALUE = 34025
const TIFFTAG_IT8BKGCOLORVALUE = 34026
const TIFFTAG_IT8PIXELINTENSITYRANGE = 34027
const TIFFTAG_IT8TRANSPARENCYINDICATOR = 34028
const TIFFTAG_IT8COLORCHARACTERIZATION = 34029
const TIFFTAG_IT8HCUSAGE = 34030
const TIFFTAG_IT8TRAPINDICATOR = 34031
const TIFFTAG_IT8CMYKEQUIVALENT = 34032
const TIFFTAG_FRAMECOUNT = 34232
const TIFFTAG_PHOTOSHOP = 34377
const TIFFTAG_EXIFIFD = 34665
const TIFFTAG_ICCPROFILE = 34675
const TIFFTAG_IMAGELAYER = 34732
const TIFFTAG_JBIGOPTIONS = 34750
const TIFFTAG_GPSIFD = 34853
const TIFFTAG_FAXRECVPARAMS = 34908
const TIFFTAG_FAXSUBADDRESS = 34909
const TIFFTAG_FAXRECVTIME = 34910
const TIFFTAG_FAXDCS = 34911
const TIFFTAG_STONITS = 37439
const TIFFTAG_FEDEX_EDR = 34929
const TIFFTAG_INTEROPERABILITYIFD = 40965
const TIFFTAG_LERC_PARAMETERS = 50674
const TIFFTAG_DNGVERSION = 50706
const TIFFTAG_DNGBACKWARDVERSION = 50707
const TIFFTAG_UNIQUECAMERAMODEL = 50708
const TIFFTAG_LOCALIZEDCAMERAMODEL = 50709
const TIFFTAG_CFAPLANECOLOR = 50710
const TIFFTAG_CFALAYOUT = 50711
const TIFFTAG_LINEARIZATIONTABLE = 50712
const TIFFTAG_BLACKLEVELREPEATDIM = 50713
const TIFFTAG_BLACKLEVEL = 50714
const TIFFTAG_BLACKLEVELDELTAH = 50715
const TIFFTAG_BLACKLEVELDELTAV = 50716
const TIFFTAG_WHITELEVEL = 50717
const TIFFTAG_DEFAULTSCALE = 50718
const TIFFTAG_DEFAULTCROPORIGIN = 50719
const TIFFTAG_DEFAULTCROPSIZE = 50720
const TIFFTAG_COLORMATRIX1 = 50721
const TIFFTAG_COLORMATRIX2 = 50722
const TIFFTAG_CAMERACALIBRATION1 = 50723
const TIFFTAG_CAMERACALIBRATION2 = 50724
const TIFFTAG_REDUCTIONMATRIX1 = 50725
const TIFFTAG_REDUCTIONMATRIX2 = 50726
const TIFFTAG_ANALOGBALANCE = 50727
const TIFFTAG_ASSHOTNEUTRAL = 50728
const TIFFTAG_ASSHOTWHITEXY = 50729
const TIFFTAG_BASELINEEXPOSURE = 50730
const TIFFTAG_BASELINENOISE = 50731
const TIFFTAG_BASELINESHARPNESS = 50732
const TIFFTAG_BAYERGREENSPLIT = 50733
const TIFFTAG_LINEARRESPONSELIMIT = 50734
const TIFFTAG_CAMERASERIALNUMBER = 50735
const TIFFTAG_LENSINFO = 50736
const TIFFTAG_CHROMABLURRADIUS = 50737
const TIFFTAG_ANTIALIASSTRENGTH = 50738
const TIFFTAG_SHADOWSCALE = 50739
const TIFFTAG_DNGPRIVATEDATA = 50740
const TIFFTAG_MAKERNOTESAFETY = 50741
const TIFFTAG_CALIBRATIONILLUMINANT1 = 50778
const TIFFTAG_CALIBRATIONILLUMINANT2 = 50779
const TIFFTAG_BESTQUALITYSCALE = 50780
const TIFFTAG_RAWDATAUNIQUEID = 50781
const TIFFTAG_ORIGINALRAWFILENAME = 50827
const TIFFTAG_ORIGINALRAWFILEDATA = 50828
const TIFFTAG_ACTIVEAREA = 50829
const TIFFTAG_MASKEDAREAS = 50830
const TIFFTAG_ASSHOTICCPROFILE = 50831
const TIFFTAG_ASSHOTPREPROFILEMATRIX = 50832
const TIFFTAG_CURRENTICCPROFILE = 50833
const TIFFTAG_CURRENTPREPROFILEMATRIX = 50834
const TIFFTAG_DCSHUESHIFTVALUES = 65535
const TIFFTAG_FAXMODE = 65536
const FAXMODE_CLASSIC = 0x0000
const FAXMODE_NORTC = 0x0001
const FAXMODE_NOEOL = 0x0002
const FAXMODE_BYTEALIGN = 0x0004
const FAXMODE_WORDALIGN = 0x0008
const FAXMODE_CLASSF = FAXMODE_NORTC
const TIFFTAG_JPEGQUALITY = 65537
const TIFFTAG_JPEGCOLORMODE = 65538
const JPEGCOLORMODE_RAW = 0x0000
const JPEGCOLORMODE_RGB = 0x0001
const TIFFTAG_JPEGTABLESMODE = 65539
const JPEGTABLESMODE_QUANT = 0x0001
const JPEGTABLESMODE_HUFF = 0x0002
const TIFFTAG_FAXFILLFUNC = 65540
const TIFFTAG_PIXARLOGDATAFMT = 65549
const PIXARLOGDATAFMT_8BIT = 0
const PIXARLOGDATAFMT_8BITABGR = 1
const PIXARLOGDATAFMT_11BITLOG = 2
const PIXARLOGDATAFMT_12BITPICIO = 3
const PIXARLOGDATAFMT_16BIT = 4
const PIXARLOGDATAFMT_FLOAT = 5
const TIFFTAG_DCSIMAGERTYPE = 65550
const DCSIMAGERMODEL_M3 = 0
const DCSIMAGERMODEL_M5 = 1
const DCSIMAGERMODEL_M6 = 2
const DCSIMAGERFILTER_IR = 0
const DCSIMAGERFILTER_MONO = 1
const DCSIMAGERFILTER_CFA = 2
const DCSIMAGERFILTER_OTHER = 3
const TIFFTAG_DCSINTERPMODE = 65551
const DCSINTERPMODE_NORMAL = 0x00
const DCSINTERPMODE_PREVIEW = 0x01
const TIFFTAG_DCSBALANCEARRAY = 65552
const TIFFTAG_DCSCORRECTMATRIX = 65553
const TIFFTAG_DCSGAMMA = 65554
const TIFFTAG_DCSTOESHOULDERPTS = 65555
const TIFFTAG_DCSCALIBRATIONFD = 65556
const TIFFTAG_ZIPQUALITY = 65557
const TIFFTAG_PIXARLOGQUALITY = 65558
const TIFFTAG_DCSCLIPRECTANGLE = 65559
const TIFFTAG_SGILOGDATAFMT = 65560
const SGILOGDATAFMT_FLOAT = 0
const SGILOGDATAFMT_16BIT = 1
const SGILOGDATAFMT_RAW = 2
const SGILOGDATAFMT_8BIT = 3
const TIFFTAG_SGILOGENCODE = 65561
const SGILOGENCODE_NODITHER = 0
const SGILOGENCODE_RANDITHER = 1
const TIFFTAG_LZMAPRESET = 65562
const TIFFTAG_PERSAMPLE = 65563
const PERSAMPLE_MERGED = 0
const PERSAMPLE_MULTI = 1
const TIFFTAG_ZSTD_LEVEL = 65564
const TIFFTAG_LERC_VERSION = 65565
const LERC_VERSION_2_4 = 4
const TIFFTAG_LERC_ADD_COMPRESSION = 65566
const LERC_ADD_COMPRESSION_NONE = 0
const LERC_ADD_COMPRESSION_DEFLATE = 1
const LERC_ADD_COMPRESSION_ZSTD = 2
const TIFFTAG_LERC_MAXZERROR = 65567
const TIFFTAG_WEBP_LEVEL = 65568
const TIFFTAG_WEBP_LOSSLESS = 65569
const EXIFTAG_EXPOSURETIME = 33434
const EXIFTAG_FNUMBER = 33437
const EXIFTAG_EXPOSUREPROGRAM = 34850
const EXIFTAG_SPECTRALSENSITIVITY = 34852
const EXIFTAG_ISOSPEEDRATINGS = 34855
const EXIFTAG_OECF = 34856
const EXIFTAG_EXIFVERSION = 36864
const EXIFTAG_DATETIMEORIGINAL = 36867
const EXIFTAG_DATETIMEDIGITIZED = 36868
const EXIFTAG_COMPONENTSCONFIGURATION = 37121
const EXIFTAG_COMPRESSEDBITSPERPIXEL = 37122
const EXIFTAG_SHUTTERSPEEDVALUE = 37377
const EXIFTAG_APERTUREVALUE = 37378
const EXIFTAG_BRIGHTNESSVALUE = 37379
const EXIFTAG_EXPOSUREBIASVALUE = 37380
const EXIFTAG_MAXAPERTUREVALUE = 37381
const EXIFTAG_SUBJECTDISTANCE = 37382
const EXIFTAG_METERINGMODE = 37383
const EXIFTAG_LIGHTSOURCE = 37384
const EXIFTAG_FLASH = 37385
const EXIFTAG_FOCALLENGTH = 37386
const EXIFTAG_SUBJECTAREA = 37396
const EXIFTAG_MAKERNOTE = 37500
const EXIFTAG_USERCOMMENT = 37510
const EXIFTAG_SUBSECTIME = 37520
const EXIFTAG_SUBSECTIMEORIGINAL = 37521
const EXIFTAG_SUBSECTIMEDIGITIZED = 37522
const EXIFTAG_FLASHPIXVERSION = 40960
const EXIFTAG_COLORSPACE = 40961
const EXIFTAG_PIXELXDIMENSION = 40962
const EXIFTAG_PIXELYDIMENSION = 40963
const EXIFTAG_RELATEDSOUNDFILE = 40964
const EXIFTAG_FLASHENERGY = 41483
const EXIFTAG_SPATIALFREQUENCYRESPONSE = 41484
const EXIFTAG_FOCALPLANEXRESOLUTION = 41486
const EXIFTAG_FOCALPLANEYRESOLUTION = 41487
const EXIFTAG_FOCALPLANERESOLUTIONUNIT = 41488
const EXIFTAG_SUBJECTLOCATION = 41492
const EXIFTAG_EXPOSUREINDEX = 41493
const EXIFTAG_SENSINGMETHOD = 41495
const EXIFTAG_FILESOURCE = 41728
const EXIFTAG_SCENETYPE = 41729
const EXIFTAG_CFAPATTERN = 41730
const EXIFTAG_CUSTOMRENDERED = 41985
const EXIFTAG_EXPOSUREMODE = 41986
const EXIFTAG_WHITEBALANCE = 41987
const EXIFTAG_DIGITALZOOMRATIO = 41988
const EXIFTAG_FOCALLENGTHIN35MMFILM = 41989
const EXIFTAG_SCENECAPTURETYPE = 41990
const EXIFTAG_GAINCONTROL = 41991
const EXIFTAG_CONTRAST = 41992
const EXIFTAG_SATURATION = 41993
const EXIFTAG_SHARPNESS = 41994
const EXIFTAG_DEVICESETTINGDESCRIPTION = 41995
const EXIFTAG_SUBJECTDISTANCERANGE = 41996
const EXIFTAG_IMAGEUNIQUEID = 42016
const int8 = UInt8
const uint8 = Cuchar
const int16 = Int16
const uint16 = UInt16
const int32 = Cint
const uint32 = UInt32
const int64 = Clong
const uint64 = Culong
const uint16_vap = Cint

struct TIFFHeaderCommon
    tiff_magic::uint16
    tiff_version::uint16
end

struct TIFFHeaderClassic
    tiff_magic::uint16
    tiff_version::uint16
    tiff_diroff::uint32
end

struct TIFFHeaderBig
    tiff_magic::uint16
    tiff_version::uint16
    tiff_offsetsize::uint16
    tiff_unused::uint16
    tiff_diroff::uint64
end

@cenum TIFFDataType::UInt32 begin
    TIFF_NOTYPE = 0
    TIFF_BYTE = 1
    TIFF_ASCII = 2
    TIFF_SHORT = 3
    TIFF_LONG = 4
    TIFF_RATIONAL = 5
    TIFF_SBYTE = 6
    TIFF_UNDEFINED = 7
    TIFF_SSHORT = 8
    TIFF_SLONG = 9
    TIFF_SRATIONAL = 10
    TIFF_FLOAT = 11
    TIFF_DOUBLE = 12
    TIFF_IFD = 13
    TIFF_LONG8 = 16
    TIFF_SLONG8 = 17
    TIFF_IFD8 = 18
end

# Skipping MacroDefinition: TIFF_INT16_T signed short
# Skipping MacroDefinition: TIFF_INT32_T signed int
# Skipping MacroDefinition: TIFF_INT64_T signed long
# Skipping MacroDefinition: TIFF_INT8_T signed char
# Skipping MacroDefinition: TIFF_UINT16_T unsigned short
# Skipping MacroDefinition: TIFF_UINT32_T unsigned int
# Skipping MacroDefinition: TIFF_UINT64_T unsigned long
# Skipping MacroDefinition: TIFF_UINT8_T unsigned char
# Skipping MacroDefinition: TIFF_SSIZE_T signed long

const TIFF_PTRDIFF_T = ptrdiff_t
const HAVE_IEEEFP = 1
const HOST_FILLORDER = FILLORDER_LSB2MSB
const HOST_BIGENDIAN = 0
const CCITT_SUPPORT = 1
const LOGLUV_SUPPORT = 1
const LZW_SUPPORT = 1
const NEXT_SUPPORT = 1
const PACKBITS_SUPPORT = 1
const PIXARLOG_SUPPORT = 1
const THUNDER_SUPPORT = 1
const ZIP_SUPPORT = 1
const STRIPCHOP_DEFAULT = TIFF_STRIPCHOP
const SUBIFD_SUPPORT = 1
const DEFAULT_EXTRASAMPLE_AS_ALPHA = 1
const CHECK_JPEG_YCBCR_SUBSAMPLING = 1
const MDI_SUPPORT = 1
const TIFFPRINT_NONE = 0x00
const TIFFPRINT_STRIPS = 0x01
const TIFFPRINT_CURVES = 0x02
const TIFFPRINT_COLORMAP = 0x04
const TIFFPRINT_JPEGQTABLES = 0x0100
const TIFFPRINT_JPEGACTABLES = 0x0200
const TIFFPRINT_JPEGDCTABLES = 0x0200
const D65_X0 = Float32(95.047)
const D65_Y0 = Float32(100.0)
const D65_Z0 = Float32(108.8827)
const D50_X0 = Float32(96.425)
const D50_Y0 = Float32(100.0)
const D50_Z0 = Float32(82.468)
const CIELABTORGB_TABLE_RANGE = 1500

# Skipping MacroDefinition: TIFFGetR ( abgr ) ( ( abgr ) & 0xff )
# Skipping MacroDefinition: TIFFGetG ( abgr ) ( ( ( abgr ) >> 8 ) & 0xff )
# Skipping MacroDefinition: TIFFGetB ( abgr ) ( ( ( abgr ) >> 16 ) & 0xff )
# Skipping MacroDefinition: TIFFGetA ( abgr ) ( ( ( abgr ) >> 24 ) & 0xff )

const LOGLUV_PUBLIC = 1
const TIFF_ANY = TIFF_NOTYPE
const TIFF_VARIABLE = -1
const TIFF_SPP = -2
const TIFF_VARIABLE2 = -3
const FIELD_CUSTOM = 65
const U_NEU = 0.210526316
const V_NEU = 0.473684211
const UVSCALE = 410.0
const tiff = Cvoid
const TIFF = tiff
const tmsize_t = Clong
const toff_t = uint64
const ttag_t = uint32
const tdir_t = uint16
const tsample_t = uint16
const tstrile_t = uint32
const tstrip_t = tstrile_t
const ttile_t = tstrile_t
const tsize_t = tmsize_t
const tdata_t = Ptr{Cvoid}
const thandle_t = Ptr{Cvoid}
const TIFFRGBValue = Cuchar

struct TIFFDisplay
    d_mat::NTuple{3, NTuple{3, Cfloat}}
    d_YCR::Cfloat
    d_YCG::Cfloat
    d_YCB::Cfloat
    d_Vrwr::uint32
    d_Vrwg::uint32
    d_Vrwb::uint32
    d_Y0R::Cfloat
    d_Y0G::Cfloat
    d_Y0B::Cfloat
    d_gammaR::Cfloat
    d_gammaG::Cfloat
    d_gammaB::Cfloat
end

struct TIFFYCbCrToRGB
    clamptab::Ptr{TIFFRGBValue}
    Cr_r_tab::Ptr{Cint}
    Cb_b_tab::Ptr{Cint}
    Cr_g_tab::Ptr{int32}
    Cb_g_tab::Ptr{int32}
    Y_tab::Ptr{int32}
end

struct TIFFCIELabToRGB
    range::Cint
    rstep::Cfloat
    gstep::Cfloat
    bstep::Cfloat
    X0::Cfloat
    Y0::Cfloat
    Z0::Cfloat
    display::TIFFDisplay
    Yr2r::NTuple{1501, Cfloat}
    Yg2g::NTuple{1501, Cfloat}
    Yb2b::NTuple{1501, Cfloat}
end

struct ANONYMOUS1_put
    any::Ptr{Cvoid}
end

struct _TIFFRGBAImage
    tif::Ptr{TIFF}
    stoponerr::Cint
    isContig::Cint
    alpha::Cint
    width::uint32
    height::uint32
    bitspersample::uint16
    samplesperpixel::uint16
    orientation::uint16
    req_orientation::uint16
    photometric::uint16
    redcmap::Ptr{uint16}
    greencmap::Ptr{uint16}
    bluecmap::Ptr{uint16}
    get::Ptr{Cvoid}
    put::ANONYMOUS1_put
    Map::Ptr{TIFFRGBValue}
    BWmap::Ptr{Ptr{uint32}}
    PALmap::Ptr{Ptr{uint32}}
    ycbcr::Ptr{TIFFYCbCrToRGB}
    cielab::Ptr{TIFFCIELabToRGB}
    UaToAa::Ptr{uint8}
    Bitdepth16To8::Ptr{uint8}
    row_offset::Cint
    col_offset::Cint
end

const TIFFRGBAImage = _TIFFRGBAImage
const tileContigRoutine = Ptr{Cvoid}
const tileSeparateRoutine = Ptr{Cvoid}

struct ANONYMOUS2_put
    any::Ptr{Cvoid}
end

const TIFFInitMethod = Ptr{Cvoid}

struct TIFFCodec
    name::Cstring
    scheme::uint16
    init::TIFFInitMethod
end

const TIFFErrorHandler = Ptr{Cvoid}
const TIFFErrorHandlerExt = Ptr{Cvoid}
const TIFFReadWriteProc = Ptr{Cvoid}
const TIFFSeekProc = Ptr{Cvoid}
const TIFFCloseProc = Ptr{Cvoid}
const TIFFSizeProc = Ptr{Cvoid}
const TIFFMapFileProc = Ptr{Cvoid}
const TIFFUnmapFileProc = Ptr{Cvoid}
const TIFFExtendProc = Ptr{Cvoid}
const _TIFFField = Cvoid
const TIFFField = _TIFFField
const _TIFFFieldArray = Cvoid
const TIFFFieldArray = _TIFFFieldArray
const TIFFVSetMethod = Ptr{Cvoid}
const TIFFVGetMethod = Ptr{Cvoid}
const TIFFPrintMethod = Ptr{Cvoid}

struct TIFFTagMethods
    vsetfield::TIFFVSetMethod
    vgetfield::TIFFVGetMethod
    printdir::TIFFPrintMethod
end

struct TIFFFieldInfo
    field_tag::ttag_t
    field_readcount::Int16
    field_writecount::Int16
    field_type::TIFFDataType
    field_bit::UInt16
    field_oktochange::Cuchar
    field_passcount::Cuchar
    field_name::Cstring
end

const TIFFLIB_VERSION_STR = "LIBTIFF, Version 4.0.10\nCopyright (c) 1988-1996 Sam Leffler\nCopyright (c) 1991-1996 Silicon Graphics, Inc."
const TIFFLIB_VERSION = 20181110


text:
export IPEPS, SquareIPEPS
"""
    IPEPS{LT<:AbstractLattice, T, N}

Infinite projected entangled pair of states.
`LT` is the type of lattice, `T` and `N` are bulk tensor element type and order.
"""
struct IPEPS{LT<:AbstractLattice, T, N, AT<:AbstractArray{T, N}}
    bulk::AT
    # TODO: check input size in constructor
end
IPEPS{LT}(bulk::AT) where {LT,T,N,AT<:AbstractArray{T,N}} = IPEPS{LT,T,N,AT}(bulk)

############### IPEPS on square lattice ###################
# size of bulk is `d × d × d × d × s`
const SquareIPEPS{T} = IPEPS{SquareLattice, T, 5}
function SquareIPEPS(bulk::AT) where {T,AT<:AbstractArray{T, 5}}
    # NOTE: from here, wrapping `bulk` with a `IPEPS` type can prevent programing from illegal input with incorrect size.
    size(bulk,1) == size(bulk,2) == size(bulk,3) == size(bulk,4) || throw(DimensionMismatch(
        "size of tensor error, should be `(d, d, d, d, s)`, got $(size(bulk))."))
    IPEPS{SquareLattice,T,5,AT}(bulk)
end
getd(ipeps::SquareIPEPS) = size(ipeps.bulk, 1)
gets(ipeps::SquareIPEPS) = size(ipeps.bulk, 5)

"""
    indexperm_symmetrize(ipeps::SquareIPEPS)

return a `SquareIPEPS` based on `ipeps` that is symmetric under
permutation of its virtual indices.
```
        4
        │
 1 ── ipeps ── 3
        │
        2
```
"""
function indexperm_symmetrize(ipeps::SquareIPEPS)
    x = ipeps.bulk
    x += permutedims(x, (1,4,3,2,5)) # up-down
    x += permutedims(x, (3,2,1,4,5)) # left-right
    x += permutedims(x, (2,1,4,3,5)) # diagonal
    x += permutedims(x, (4,3,2,1,5)) # rotation
    return SquareIPEPS(x / norm(x))
end


text:
int_rules_4_5_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.7*(d*trig)^m*(a+b*(c*sec)^n)^p =#
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => (~b) ^ ~p * Antiderivative(ActivateTrig(~u * tan(~e + ~f * ~x) ^ (2 * ~p)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (EqQ(~a + ~b, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => Antiderivative(ActivateTrig(~u * (~b * tan(~e + ~f * ~x) ^ 2) ^ ~p), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && EqQ(~a + ~b, 0)
    @apply_utils Antiderivative((~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~b * ff) / ~f) * Subst(Antiderivative((~b + ~b * ff ^ 2 * (~x) ^ 2) ^ (~p - 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~b, ~e, ~f, ~p], ~x) && Not(IntegerQ(~p))
    @apply_utils Antiderivative((~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (((~b) ^ IntPart(~p) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ FracPart(~p)) / (~c * sec(~e + ~f * ~x)) ^ (~n * FracPart(~p))) * Antiderivative((~c * sec(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~b, ~c, ~e, ~f, ~n, ~p], ~x) && Not(IntegerQ(~p))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => (~b / (2 * ~f)) * Subst(Antiderivative((-1 + ~x) ^ ((~m - 1) / 2) * (~b * ~x) ^ (~p - 1), ~x), ~x, sec(~e + ~f * ~x) ^ 2) <-- FreeQ([~b, ~e, ~f, ~p], ~x) && (Not(IntegerQ(~p)) && IntegerQ((~m - 1) / 2))
    @apply_utils Antiderivative(~(u') * (~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(sec(~e + ~f * ~x), ~x)], (((~b * ff ^ ~n) ^ IntPart(~p) * (~b * sec(~e + ~f * ~x) ^ ~n) ^ FracPart(~p)) / (sec(~e + ~f * ~x) / ff) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (sec(~e + ~f * ~x) / ff) ^ (~n * ~p), ~x)) <-- FreeQ([~b, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (IntegerQ(~n) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    @apply_utils Antiderivative(~(u') * (~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (((~b) ^ IntPart(~p) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ FracPart(~p)) / (~c * sec(~e + ~f * ~x)) ^ (~n * FracPart(~p))) * Antiderivative(ActivateTrig(~u) * (~c * sec(~e + ~f * ~x)) ^ (~n * ~p), ~x) <-- FreeQ([~b, ~c, ~e, ~f, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~n)) && (EqQ(~u, 1) || MatchQ(~u, (~(d') * (~trig)(~e + ~f * ~x)) ^ ~(m') <-- FreeQ([d, m], ~x) && MemberQ([sin, cos, tan, cot, sec, csc], trig)))))
    @apply_utils Antiderivative(1 / (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ 2), ~x) => ~x / ~a - (~b / ~a) * Antiderivative(1 / (~b + ~a * cos(~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && NeQ(~a + ~b, 0)
    @apply_utils Antiderivative((~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + ~b + ~b * ff ^ 2 * (~x) ^ 2) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (NeQ(~a + ~b, 0) && NeQ(~p, -1))
    @apply_utils Antiderivative((~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ 4) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + ~b + 2 * ~b * ff ^ 2 * (~x) ^ 2 + ~b * ff ^ 4 * (~x) ^ 4) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && IntegerQ(2 * ~p)
    @apply_utils Antiderivative((~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + ~b * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~n / 2) && IGtQ(~p, -2))
    @apply_utils Antiderivative((~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => Unintegrable((~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff ^ (~m + 1) / ~f) * Subst(Antiderivative(((~x) ^ ~m * ExpandToSum(~a + ~b * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 + 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-ff / ~f) * Subst(Antiderivative(((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~b + ~a * (ff * ~x) ^ ~n) ^ ~p) / (ff * ~x) ^ (~n * ~p), ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n) && IntegerQ(~p)))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (1 / (~f * ff ^ ~m)) * Subst(Antiderivative(((-1 + ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * (~c * ff * ~x) ^ ~n) ^ ~p) / (~x) ^ (~m + 1), ~x), ~x, sec(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && (IntegerQ((~m - 1) / 2) && (GtQ(~m, 0) || (EqQ(~n, 2) || EqQ(~n, 4))))
    @apply_utils Antiderivative((~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * sin(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~(n')) ^ ~(p'), ~x) => (~d) ^ (~n * ~p) * Antiderivative((~d * cos(~e + ~f * ~x)) ^ (~m - ~n * ~p) * (~b + ~a * cos(~e + ~f * ~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~m)) && IntegersQ(~n, ~p))
    @apply_utils Antiderivative((~(d') * cos(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * cos(~e + ~f * ~x)) ^ FracPart(~m) * (sec(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((sec(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => Module([ff = FreeFactors(cos(~e + ~f * ~x), ~x)], (-1 / (~f * ff ^ ((~m + ~n * ~p) - 1))) * Subst(Antiderivative(((1 - ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~b + ~a * (ff * ~x) ^ ~n) ^ ~p) / (~x) ^ (~m + ~n * ~p), ~x), ~x, cos(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~n], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n) && IntegerQ(~p)))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(sec(~e + ~f * ~x), ~x)], (1 / ~f) * Subst(Antiderivative(((-1 + ff ^ 2 * (~x) ^ 2) ^ ((~m - 1) / 2) * (~a + ~b * (~c * ff * ~x) ^ ~n) ^ ~p) / ~x, ~x), ~x, sec(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~n, ~p], ~x) && (IntegerQ((~m - 1) / 2) && (GtQ(~m, 0) || (EqQ(~n, 2) || (EqQ(~n, 4) || (IGtQ(~p, 0) || IntegersQ(2 * ~n, ~p))))))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * sec(~(e') + ~(f') * ~x) ^ 2) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], ((~b * ff) / ~f) * Subst(Antiderivative((~d * ff * ~x) ^ ~m * (~b + ~b * ff ^ 2 * (~x) ^ 2) ^ (~p - 1), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~b, ~d, ~e, ~f, ~m, ~p], ~x)
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~(p'), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(((~d * ff * ~x) ^ ~m * (~a + ~b * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p) / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~d, ~e, ~f, ~m, ~p], ~x) && (IntegerQ(~n / 2) && (IntegerQ(~m / 2) || EqQ(~n, 2)))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => (~d * (~d * tan(~e + ~f * ~x)) ^ (~m - 1) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p) / (~f * ((~p * ~n + ~m) - 1)) - (((~d) ^ 2 * (~m - 1)) / ((~p * ~n + ~m) - 1)) * Antiderivative((~d * tan(~e + ~f * ~x)) ^ (~m - 2) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~p, ~n], ~x) && (GtQ(~m, 1) && (NeQ((~p * ~n + ~m) - 1, 0) && IntegersQ(2 * ~p * ~n, 2 * ~m)))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => ((~d * tan(~e + ~f * ~x)) ^ (~m + 1) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p) / (~d * ~f * (~m + 1)) - ((~p * ~n + ~m + 1) / ((~d) ^ 2 * (~m + 1))) * Antiderivative((~d * tan(~e + ~f * ~x)) ^ (~m + 2) * (~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~f, ~p, ~n], ~x) && (LtQ(~m, -1) && (NeQ(~p * ~n + ~m + 1, 0) && IntegersQ(2 * ~p * ~n, 2 * ~m)))
    @apply_utils Antiderivative((~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * tan(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * cot(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * cot(~e + ~f * ~x)) ^ FracPart(~m) * (tan(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((tan(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~m * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((1 + ff ^ 2 * (~x) ^ 2) ^ (~m / 2 - 1) * ExpandToSum(~a + ~b * (1 + ff ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) ^ ~p, ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ(~m / 2) && IntegerQ(~n / 2))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(ExpandToSum(~b + ~a * (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) ^ ~p / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + ~n * ~p + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && IntegerQ(~p)))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => With([ff = FreeFactors(sin(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative((~a + ~b / (1 - ff ^ 2 * (~x) ^ 2) ^ (~n / 2)) ^ ~p / (1 - ff ^ 2 * (~x) ^ 2) ^ ((~m + 1) / 2), ~x), ~x, sin(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~e, ~f, ~p], ~x) && (IntegerQ((~m - 1) / 2) && (IntegerQ(~n / 2) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(sec(~(e') + ~(f') * ~x) ^ ~(m') * (~a + ~(b') * sec(~(e') + ~(f') * ~x) ^ ~n) ^ ~p, ~x) => Antiderivative(ExpandTrig(sec(~e + ~f * ~x) ^ ~m * (~a + ~b * sec(~e + ~f * ~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IntegersQ(~m, ~n, ~p)
    @apply_utils Antiderivative((~(d') * sec(~(e') + ~(f') * ~x)) ^ ~(m') * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~(p'), ~x) => Unintegrable((~d * sec(~e + ~f * ~x)) ^ ~m * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~a + ~(b') * (~(c') * sec(~(e') + ~(f') * ~x)) ^ ~n) ^ ~p, ~x) => (~d * csc(~e + ~f * ~x)) ^ FracPart(~m) * (sin(~e + ~f * ~x) / ~d) ^ FracPart(~m) * Antiderivative((sin(~e + ~f * ~x) / ~d) ^ -(~m) * (~a + ~b * (~c * sec(~e + ~f * ~x)) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.7 (d trig)^m (a+b (c sec)^n)^p.jl


text:
using Bonsai
using Bonsai: Params
using CodeInfoTools
using HTTP: Stream
using Bonsai: handler_writes
using StructTypes
using Test

struct A2
    data
end

StructTypes.StructType(::Type{A2}) = StructTypes.Struct()

@testset "handler_writes" begin
    app = Bonsai.App()
    app.get("/path") do stream
        Bonsai.write(stream, Headers(header="value"))
        Bonsai.write(stream, Body(x=10))
        Bonsai.write(stream, Body(a=A2(10)))
    end
    h = match(app.paths, Params(), "GET", ["path"], 1)
    # for debuging
    # code_inferred(h.fn, Tuple{Stream})
    # Bonsai.handler_writes(h.fn)
    @test length(Bonsai.handler_writes(h.fn)) == 4
end

@testset "handler_reads" begin
    app = Bonsai.App()
    app.get("/path") do stream
        Bonsai.read(stream, Headers(A2))
        Bonsai.read(stream, Query(A2))
        Bonsai.read(stream, Body(A2))
    end
    h = match(app.paths, Params(), "GET", ["path"], 1)
    @test length(Bonsai.handler_reads(h.fn)) == 3
end

# using Bonsai, HTTP

# app = Bonsai.App()

# app.get("/data") do stream
#     @info "SERVER"
#     Bonsai.write(stream, Body(rand(10^6)))
# end

# start(app, port=7864)

# stop(app)

# res = HTTP.get("http://localhost:7864/data")

# try
#     start(app, port=7864)
#     sleep(1)
#     #http://localhost:7864/data
#     res = HTTP.get("http://localhost:7864/data")
#     @info res
# catch e
#     @error e
# finally
#     stop(app)
# end

# how to get this to pass
# app = Bonsai.App()

# "some function documenation"
# app.get("/test-docs") do stream
# end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
type FakeDriverModel <: DriverModel{FakeDriveAction, FakeActionContext} end

let
    trajdata = get_test_trajdata()
    roadway = trajdata.roadway
    veh = get_vehicle(trajdata, 1, 1)

    model = FakeDriverModel()
    reset_hidden_state!(model)
    observe!(model, Scene(), roadway, 1)
    prime_with_history!(model, trajdata, roadway, 1, 2, 1)

    @test get_name(model) == "???"
    @test action_type(model) <: FakeDriveAction
    @test_throws ErrorException action_context(model)
    @test_throws ErrorException rand(model)
    @test_throws ErrorException pdf(model, FakeDriveAction())
    @test_throws ErrorException logpdf(model, FakeDriveAction())
end


text:
import .CUDA
CUDA.libcufft()
CUDA.libcudnn()
CUDA.libcurand()
path_cublas = CUDA.libcublas()
using Libdl
# HACK
# CUDA does not yet? expose `libcudart`:
# https://discourse.julialang.org/t/shared-c-dependencies-and-artifacts/68525/4?u=jw3126
# But as an implementation detail libcudart does ship with the artifact
dir = splitdir(path_cublas)[1]
path_cudart = joinpath(dir, "libcudart.so.11.0")
dlopen(path_cudart)


text:
@testset "SSL/TLS" begin
    ssl_config = get_ssl_config(
        ssl_certfile=joinpath(@__DIR__, "..", "docker/ssl/redis.crt"),
        ssl_keyfile=joinpath(@__DIR__, "..", "docker/ssl/redis.key"),
        ssl_ca_certs=joinpath(@__DIR__, "..", "docker/ssl/ca.crt")
    )

    # Port 6380 corresponds to port defined for redis-ssl container in docker/docker-compose.yml
    client = Client(port=6380, ssl_config=ssl_config)
    set_global_client(client)

    set("test", "ssl success")
    @test get("test") == "ssl success"

    flushall()
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -4.116064765345723 <= q <= -0.12537054293327043)

                     add_NL_constraint(m, :(swish(-0.6620590692147021 + -0.9580027962881332*swish(-0.38090481240177176 + -0.02939531016329644*swish(0.2067810723024852 + 0.7469716803367321*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.9724733578803808*swish(-0.4274287128450034 + -0.22403315015562253*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))) + -0.699289524220764*swish(0.4498619401955626 + 0.5463276179102481*swish(0.2067810723024852 + 0.7469716803367321*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.02415027804893688*swish(-0.4274287128450034 + -0.22403315015562253*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) + swish(-0.8107501452143855 + 0.008470632510253484*swish(-0.38090481240177176 + -0.02939531016329644*swish(0.2067810723024852 + 0.7469716803367321*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.9724733578803808*swish(-0.4274287128450034 + -0.22403315015562253*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))) + 0.39739497887500974*swish(0.4498619401955626 + 0.5463276179102481*swish(0.2067810723024852 + 0.7469716803367321*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.3209180228281623*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))) + -0.02415027804893688*swish(-0.4274287128450034 + -0.22403315015562253*swish(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])) + -0.6136110647248842*swish(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export
    # From libcd.jl
    cdCanvas,
    cdCreateCanvas,
    cdKillCanvas,
    cdCanvasGetContext,
    cdCanvasActivate,
    cdCanvasDeactivate,
    cdUseContextPlus,
    cdInitContextPlus,
    cdFinishContextPlus,
    cdContextRegisterCallback,
    cdContextCaps,
    cdContextIsPlus,
    cdContextType,
    cdContext,
    cdContextNativeWindow,
    cdCanvasSimulate,
    cdCanvasFlush,
    cdCanvasClear,
    cdCanvasSaveState,
    cdCanvasRestoreState,
    cdReleaseState,
    cdCanvasSetAttribute,
    cdCanvasGetAttribute,
    cdCanvasPlay,
    cdCanvasGetSize,
    cdCanvasUpdateYAxis,
    cdfCanvasUpdateYAxis,
    cdCanvasInvertYAxis,
    cdfCanvasInvertYAxis,
    cdCanvasMM2Pixel,
    cdCanvasPixel2MM,
    cdfCanvasMM2Pixel,
    cdfCanvasPixel2MM,
    cdCanvasOrigin,
    cdfCanvasOrigin,
    cdCanvasGetOrigin,
    cdfCanvasGetOrigin,
    cdCanvasTransform,
    cdCanvasGetTransform,
    cdCanvasTransformMultiply,
    cdCanvasTransformRotate,
    cdCanvasTransformScale,
    cdCanvasTransformTranslate,
    cdCanvasTransformPoint,
    cdfCanvasTransformPoint,
    cdCanvasClip,
    cdCanvasClipArea,
    cdCanvasGetClipArea,
    cdfCanvasClipArea,
    cdfCanvasGetClipArea,
    cdCanvasIsPointInRegion,
    cdCanvasOffsetRegion,
    cdCanvasGetRegionBox,
    cdCanvasRegionCombineMode,
    cdCanvasPixel,
    cdCanvasMark,
    cdCanvasBegin,
    cdCanvasPathSet,
    cdCanvasEnd,
    cdCanvasLine,
    cdCanvasVertex,
    cdCanvasRect,
    cdCanvasBox,
    cdCanvasArc,
    cdCanvasSector,
    cdCanvasChord,
    cdCanvasText,
    cdfCanvasLine,
    cdfCanvasVertex,
    cdfCanvasRect,
    cdfCanvasBox,
    cdfCanvasArc,
    cdfCanvasSector,
    cdfCanvasChord,
    cdfCanvasText,
    cdCanvasSetBackground,
    cdCanvasSetForeground,
    cdCanvasBackground,
    cdCanvasForeground,
    cdCanvasBackOpacity,
    cdCanvasWriteMode,
    cdCanvasLineStyle,
    cdCanvasLineStyleDashes,
    cdCanvasLineWidth,
    cdCanvasLineJoin,
    cdCanvasLineCap,
    cdCanvasInteriorStyle,
    cdCanvasHatch,
    cdCanvasStipple,
    cdCanvasGetStipple,
    cdCanvasPattern,
    cdCanvasGetPattern,
    cdCanvasFillMode,
    cdCanvasFont,
    cdCanvasGetFont,
    cdCanvasNativeFont,
    cdCanvasTextAlignment,
    cdCanvasTextOrientation,
    cdCanvasMarkType,
    cdCanvasMarkSize,
    cdCanvasVectorText,
    cdCanvasMultiLineVectorText,
    cdCanvasVectorFont,
    cdCanvasVectorTextDirection,
    cdCanvasVectorTextTransform,
    cdCanvasVectorTextSize,
    cdCanvasVectorCharSize,
    cdCanvasVectorFontSize,
    cdCanvasGetVectorFontSize,
    cdCanvasGetVectorTextSize,
    cdCanvasGetVectorTextBounds,
    cdCanvasGetVectorTextBox,
    cdfCanvasVectorTextDirection,
    cdfCanvasVectorTextSize,
    cdfCanvasGetVectorTextSize,
    cdfCanvasVectorCharSize,
    cdfCanvasVectorText,
    cdfCanvasMultiLineVectorText,
    cdfCanvasGetVectorTextBounds,
    cdfCanvasGetVectorTextBox,
    cdCanvasGetFontDim,
    cdCanvasGetTextSize,
    cdCanvasGetTextBox,
    cdCanvasGetTextBounds,
    cdCanvasGetColorPlanes,
    cdCanvasPalette,
    cdCanvasGetImageRGB,
    cdCanvasPutImageRectRGB,
    cdCanvasPutImageRectRGBA,
    cdCanvasPutImageRectMap,
    cdCanvasCreateImage,
    cdKillImage,
    cdCanvasGetImage,
    cdCanvasPutImageRect,
    cdCanvasScrollArea,
    cdCreateBitmap,
    cdKillBitmap,
    cdBitmapGetData,
    cdBitmapSetRect,
    cdCanvasPutBitmap,
    cdCanvasGetBitmap,
    cdBitmapRGB2Map,
    cdEncodeColor,
    cdDecodeColor,
    cdDecodeAlpha,
    cdEncodeAlpha,
    cdRGB2Map,
    cdActivate,
    cdActiveCanvas,
    cdSimulate,
    cdFlush,
    cdClear,
    cdSaveState,
    cdRestoreState,
    cdSetAttribute,
    cdGetAttribute,
    cdGetContext,
    cdRegisterCallback,
    cdPlay,
    cdGetCanvasSize,
    cdUpdateYAxis,
    cdMM2Pixel,
    cdPixel2MM,
    cdOrigin,
    cdClip,
    cdGetClipPoly,
    cdClipArea,
    cdGetClipArea,
    cdPointInRegion,
    cdOffsetRegion,
    cdRegionBox,
    cdRegionCombineMode,
    cdPixel,
    cdMark,
    cdLine,
    cdBegin,
    cdVertex,
    cdEnd,
    cdRect,
    cdBox,
    cdArc,
    cdSector,
    cdChord,
    cdText,
    cdBackground,
    cdForeground,
    cdBackOpacity,
    cdWriteMode,
    cdLineStyle,
    cdLineStyleDashes,
    cdLineWidth,
    cdLineJoin,
    cdLineCap,
    cdInteriorStyle,
    cdHatch,
    cdStipple,
    cdGetStipple,
    cdPattern,
    cdGetPattern,
    cdFillMode,
    cdFont,
    cdGetFont,
    cdNativeFont,
    cdTextAlignment,
    cdTextOrientation,
    cdMarkType,
    cdMarkSize,
    cdMultiLineVectorText,
    cdGetVectorTextSize,
    cdGetVectorTextBounds,
    cdFontDim,
    cdImage,
    cdTextSize,
    cdTextBox,
    cdTextBounds,
    cdGetColorPlanes,
    cdPalette,
    cdGetImageRGB,
    cdPutImageRectRGB,
    cdPutImageRectRGBA,
    cdPutImageRectMap,
    cdCreateImage,
    cdGetImage,
    cdPutImageRect,
    cdScrollArea,
    cdPutBitmap,
    cdGetBitmap,
    cdVersion,
    cdVersionDate,
    cdVersionNumber,
    cdVectorText,
#    cdVectorFont,
    cdVectorTextDirection,
    cdVectorTextTransform,
    cdVectorTextSize,
    cdVectorCharSize,
    #
    # From wd.jl
    #
    wdCanvasWindow,
    wdCanvasGetWindow,
    wdCanvasViewport,
    wdCanvasGetViewport,
    wdCanvasWorld2Canvas,
    wdCanvasWorld2CanvasSize,
    wdCanvasCanvas2World,
    wdCanvasSetTransform,
    wdCanvasGetTransform,
    wdCanvasTranslate,
    wdCanvasScale,
    wdCanvasClipArea,
    wdCanvasGetClipArea,
    wdCanvasIsPointInRegion,
    wdCanvasOffsetRegion,
    wdCanvasGetRegionBox,
    wdCanvasHardcopy,
    wdCanvasPixel,
    wdCanvasMark,
    wdCanvasLine,
    wdCanvasVertex,
    wdCanvasRect,
    wdCanvasBox,
    wdCanvasArc,
    wdCanvasSector,
    wdCanvasChord,
    wdCanvasText,
    wdCanvasPutImageRect,
    wdCanvasPutImageRectRGB,
    wdCanvasPutImageRectRGBA,
    wdCanvasPutImageRectMap,
    wdCanvasPutBitmap,
    wdCanvasLineWidth,
    wdCanvasFont,
    wdCanvasGetFont,
    wdCanvasMarkSize,
    wdCanvasGetFontDim,
    wdCanvasGetTextSize,
    wdCanvasGetTextBox,
    wdCanvasGetTextBounds,
    wdCanvasStipple,
    wdCanvasPattern,
    wdCanvasVectorTextDirection,
    wdCanvasVectorTextSize,
    wdCanvasGetVectorTextSize,
    wdCanvasVectorCharSize,
    wdCanvasVectorText,
    wdCanvasMultiLineVectorText,
    wdCanvasGetVectorTextBounds,
    wdCanvasGetVectorTextBox,
    wdWindow,
    wdGetWindow,
    wdViewport,
    wdGetViewport,
    wdWorld2Canvas,
    wdWorld2CanvasSize,
    wdCanvas2World,
    wdClipArea,
    wdGetClipArea,
    wdGetClipPoly,
    wdPointInRegion,
    wdOffsetRegion,
    wdRegionBox,
    wdHardcopy,
    wdPixel,
    wdMark,
    wdLine,
    wdVertex,
    wdRect,
    wdBox,
    wdArc,
    wdSector,
    wdChord,
    wdText,
    wdPutImageRect,
    wdPutImageRectRGB,
    wdPutImageRectRGBA,
    wdPutImageRectMap,
    wdPutBitmap,
    wdLineWidth,
    wdFont,
    wdGetFont,
    wdMarkSize,
    wdFontDim,
    wdTextSize,
    wdTextBox,
    wdTextBounds,
    wdStipple,
    wdPattern,
    wdVectorTextDirection,
    wdVectorTextSize,
    wdGetVectorTextSize,
    wdVectorCharSize,
    wdVectorText,
    wdMultiLineVectorText,
    wdGetVectorTextBound,
    #
    # From cdiup.jl
    #
    cdContextIup,
    # From cdxdbuf
    cdContextDBuffer,
    #
    # From handlegraphics
    #
    guidata, setappdata, getappdata, WindowButtonDownFcn, WindowButtonMotionFcn, KeyPressFcn,
    set_current_point, get_current_point, set_current_character, get_current_character, isnada, Handles,
    gc_ref, gc_unref,
    # local
    imcdCanvasPutImage

include("libcd.jl")
include("cdiup.jl")
include("wd.jl")
include("handlegraphics.jl")

# ----------------------------------------------------------------------------------------
# Utility function to draw the image in a CD library canvas (from im_image.h)
# Works only for data_type IM_BYTE, and color spaces: IM_RGB, IM_MAP, IMGRAY and IM_BINARY.
function imcdCanvasPutImage(_canvas, _image, _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
    if (isa(_image, Ptr{imImage}))
        _image = unsafe_load(_image)
    end
    IM_RGB = 0        # Redefining here to not have to import it from libim.h (Ghrrr)
    if (_image.color_space == IM_RGB)
        data = [convert(Ptr{UInt8}, unsafe_load(_image.data,1)),
            convert(Ptr{UInt8}, unsafe_load(_image.data,2)),
            convert(Ptr{UInt8}, unsafe_load(_image.data,3))]

        if (_image.has_alpha != 0)
            alfa = convert(Ptr{UInt8}, unsafe_load(_image.data,4))
            cdCanvasPutImageRectRGBA(_canvas, _image.width, _image.height, data[1], data[2],
                                     data[3], alfa, _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
        else
            cdCanvasPutImageRectRGB(_canvas, _image.width, _image.height, data[1], data[2],
                                    data[3], _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
        end
    else
        data = [convert(Ptr{UInt8}, unsafe_load(_image.data,1))]
        cdCanvasPutImageRectMap(_canvas, _image.width, _image.height, data[1], _image.palette,
                                _x, _y, _w, _h, _xmin, _xmax, _ymin, _ymax)
    end
end


text:
# Functions for representing graphs in GraphViz's dot format
# http://www.graphviz.org/
# http://www.graphviz.org/Documentation/dotguide.pdf
# http://www.graphviz.org/pub/scm/graphviz2/doc/info/lang.html

# Write the dot representation of a graph to a file by name.
function to_dot(graph::AbstractGraph, filename::String)
    open(filename,"w") do f
        to_dot(graph, f)
    end
end

# Get the dot representation of a graph as a string.
function to_dot(graph::AbstractGraph)
    str = IOBuffer()
    to_dot(graph, str)
    takebuf_string(str)
end

# Write the dot representation of a graph to a stream.
function to_dot{G<:AbstractGraph}(graph::G, stream::IO)
    has_vertex_attrs = method_exists(attributes, (vertex_type(graph), G))
    has_edge_attrs = method_exists(attributes, (edge_type(graph), G))

    write(stream, "$(graph_type_string(graph)) graphname {\n")
    if implements_edge_list(graph)
        for edge in edges(graph)
            write(stream,"$(vertex_index(source(edge), graph)) $(edge_op(graph)) $(vertex_index(target(edge), graph))\n")
        end
    elseif implements_vertex_list(graph) && (implements_incidence_list(graph) || implements_adjacency_list(graph))
        for vertex in vertices(graph)
            if has_vertex_attrs && !isempty(attributes(vertex, graph))
                write(stream, "$(vertex_index(vertex, graph)) $(to_dot(attributes(vertex, graph)))\n")
            end
            if implements_incidence_list(graph)
                for e in out_edges(vertex, graph)
                    n = target(e, graph)
                    if is_directed(graph) || vertex_index(n, graph) > vertex_index(vertex, graph)
                        write(stream,"$(vertex_index(vertex, graph)) $(edge_op(graph)) $(vertex_index(n, graph))$(has_edge_attrs ? string(" ", to_dot(attributes(e, graph))) : "")\n")
                    end
                end
            else # implements_adjacency_list
                for n in out_neighbors(vertex, graph)
                    if is_directed(graph) || vertex_index(n, graph) > vertex_index(vertex, graph)
                        write(stream,"$(vertex_index(vertex, graph)) $(edge_op(graph)) $(vertex_index(n))\n")
                    end
                end
            end
        end
    else
        throw(ArgumentError("More graph Concepts needed: dot serialization requires iteration over edges or iteration over vertices and neighbors."))
    end
    write(stream, "}\n")
    stream
end

function to_dot(attrs::AttributeDict)
    if isempty(attrs)
        ""
    else
        string("[",join(map(to_dot,collect(attrs)),","),"]")
    end
end

to_dot(attr_tuple::(UTF8String, Any)) = "\"$(attr_tuple[1])\"=\"$(attr_tuple[2])\""

function graph_type_string(graph::AbstractGraph)
    is_directed(graph) ? "digraph" : "graph"
end

function edge_op(graph::AbstractGraph)
    is_directed(graph) ? "->" : "--"
end

function plot(g::AbstractGraph)
    stdin, proc = writesto(`neato -Tx11`)
    to_dot(g, stdin)
    close(stdin)
end


text:
"""
    maxwellian(u, ρ, U, λ)
    maxwellian(u, prim)
    maxwellian(u, v, ρ, U, V, λ)
    maxwellian(u, v, prim)
    maxwellian(u, v, w, ρ, U, V, W, λ)
    maxwellian(u, v, w, prim)

Maxwellian in discrete form

* @args: particle velocity quadrature points
* @args: density, velocity and inverse of temperature
* @return: Maxwellian distribution function

"""
maxwellian(u, ρ, U, λ) = @. ρ * sqrt(λ / π) * exp(-λ * (u - U)^2) # 1V

maxwellian(u::Union{Real,AbstractVector{T}}, prim) where {T<:Real} =
    maxwellian(u, prim[1], prim[2], prim[end]) # in case of input with length 4/5

#--- 2V ---#
maxwellian(u, v, ρ, U, V, λ) = @. ρ * (λ / π) * exp(-λ * ((u - U)^2 + (v - V)^2))

maxwellian(u, v, prim::AbstractVector{T}) where {T<:Real} =
    maxwellian(u, v, prim[1], prim[2], prim[3], prim[end]) # in case of input with length 5

#--- 3V ---#
maxwellian(u, v, w, ρ, U, V, W, λ) =
    @. ρ * sqrt((λ / π)^3) * exp(-λ * ((u - U)^2 + (v - V)^2 + (w - W)^2))

maxwellian(u, v, w, prim::AbstractVector{T}) where {T<:Real} =
    maxwellian(u, v, w, prim[1], prim[2], prim[3], prim[4], prim[5])

"""
    * Maxwellian
    maxwellian!(M, u, ρ, U, λ)
    maxwellian!(M, u, ρ, U, V, λ)
    maxwellian!(M, u, ρ, U, V, W, λ)
    maxwellian!(M, u, prim)
    maxwellian!(M, u, v, prim)
    maxwellian!(M, u, v, w, prim)

    * Rykov
    maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, u, prim, K, Kr)
    maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, u, v, prim, K, Kr)

In-place Maxwellian

* @args: particle velocity quadrature points
* @args: density, velocity and inverse of temperature
* @return: Maxwellian distribution function

"""
function maxwellian!(
    M::AbstractVector{T1},
    u::AbstractVector{T2},
    ρ,
    U,
    λ,
) where {T1<:AbstractFloat,T2<:AbstractFloat}
    @. M = ρ * sqrt(λ / π) * exp(-λ * (u - U)^2) # 1V
    return nothing
end

maxwellian!(
    M::AbstractVector{T1},
    u::AbstractVector{T2},
    prim::AbstractVector{T3},
) where {T1<:AbstractFloat,T2<:AbstractFloat,T3<:Real} =
    maxwellian!(M, u, prim[1], prim[2], prim[end])

# Rykov
function maxwellian!(
    Ht::T1,
    Bt::T1,
    Rt::T1,
    Hr::T1,
    Br::T1,
    Rr::T1,
    u::T2,
    prim::T3,
    K,
    Kr,
) where {
    T1<:AbstractVector{<:AbstractFloat},
    T2<:AbstractVector{<:AbstractFloat},
    T3<:AbstractVector{<:Real},
}
    @. Ht = prim[1] * sqrt(prim[4] / π) * exp(-prim[4] * (u - prim[2])^2)
    @. Bt = Ht * K / (2.0 * prim[4])
    @. Rt = Ht * Kr / (2.0 * prim[5])

    @. Hr = prim[1] * sqrt(prim[3] / π) * exp(-prim[3] * (u - prim[2])^2)
    @. Br = Hr * K / (2.0 * prim[3])
    @. Rr = Hr * Kr / (2.0 * prim[3])

    return nothing
end

#--- 2V ---#
function maxwellian!(
    M::T1,
    u::T2,
    v::T2,
    ρ,
    U,
    V,
    λ,
) where {T1<:AbstractArray{<:AbstractFloat,2},T2<:AbstractArray{<:AbstractFloat}}
    @. M = ρ * (λ / π) * exp(-λ * ((u - U)^2 + (v - V)^2))
    return nothing
end

maxwellian!(
    M::T1,
    u::T2,
    v::T2,
    prim::T3,
) where {
    T1<:AbstractArray{<:AbstractFloat,2},
    T2<:AbstractArray{<:AbstractFloat},
    T3<:AbstractArray{<:Real,1},
} = maxwellian!(M, u, v, prim[1], prim[2], prim[3], prim[end])

# Rykov
function maxwellian!(
    Ht::T1,
    Bt::T1,
    Rt::T1,
    Hr::T1,
    Br::T1,
    Rr::T1,
    u::T2,
    v::T2,
    prim::T3,
    K,
    Kr,
) where {
    T1<:AbstractArray{<:AbstractFloat,2},
    T2<:AbstractArray{<:AbstractFloat},
    T3<:AbstractArray{<:Real,1},
}
    @. Ht = prim[1] * (prim[5] / π) * exp(-prim[5] * ((u - prim[2])^2 + (v - prim[3])^2))
    @. Bt = Ht * K / (2.0 * prim[5])
    @. Rt = Ht * Kr / (2.0 * prim[6])

    @. Hr = prim[1] * (prim[4] / π) * exp(-prim[4] * ((u - prim[2])^2 + (v - prim[3])^2))
    @. Br = Hr * K / (2.0 * prim[4])
    @. Rr = Hr * Kr / (2.0 * prim[4])

    return nothing
end

#--- 3V ---#
function maxwellian!(
    M::T1,
    u::T2,
    v::T2,
    w::T2,
    ρ,
    U,
    V,
    W,
    λ,
) where {T1<:AbstractArray{<:AbstractFloat,3},T2<:AbstractArray{<:AbstractFloat}}
    @. M = ρ * sqrt((λ / π)^3) * exp(-λ * ((u - U)^2 + (v - V)^2 + (w - W)^2))
    return nothing
end

maxwellian!(
    M::T1,
    u::T2,
    v::T2,
    w::T2,
    prim::T3,
) where {
    T1<:AbstractArray{<:AbstractFloat,3},
    T2<:AbstractArray{<:AbstractFloat},
    T3<:AbstractArray{<:Real,1},
} = maxwellian!(M, u, v, w, prim[1], prim[2], prim[3], prim[4], prim[5])

"""
    mixture_maxwellian(u, prim)
    mixture_maxwellian(u, v, prim)
    mixture_maxwellian(u, v, w, prim)

Multi-component Maxwellian in discrete form
"""
function mixture_maxwellian(
    u::AbstractMatrix{T1},
    prim::AbstractMatrix{T2},
) where {T1<:Real,T2<:Real}

    mixM = similar(u)
    for j in axes(mixM, 2)
        mixM[:, j] .= maxwellian(u[:, j], prim[:, j])
    end

    return mixM

end

function mixture_maxwellian(
    u::X,
    v::X,
    prim::AbstractMatrix{Y},
) where {X<:AbstractArray{<:AbstractFloat,3},Y<:Real}

    mixM = similar(u)
    for k in axes(mixM, 3)
        mixM[:, :, k] .= maxwellian(u[:, :, k], v[:, :, k], prim[:, k])
    end

    return mixM

end

function mixture_maxwellian(
    u::X,
    v::X,
    w::X,
    prim::AbstractMatrix{Y},
) where {X<:AbstractArray{<:AbstractFloat,4},Y<:Real}

    mixM = similar(u)
    for l in axes(mixM, 4)
        mixM[:, :, :, l] .=
            maxwellian(u[:, :, :, l], v[:, :, :, l], w[:, :, :, l], prim[:, l])
    end

    return mixM

end

"""
    mixture_maxwellian!(M, u, prim)
    mixture_maxwellian!(M, u, v, prim)
    mixture_maxwellian!(M, u, v, w, prim)

In-place multi-component Maxwellian

"""
function mixture_maxwellian!(
    M::AbstractMatrix{T1},
    u::T2,
    prim::AbstractMatrix{T3},
) where {T1<:Real,T2<:AbstractArray{<:Real,2},T3<:Real}

    for j in axes(M, 2)
        _M = @view M[:, j]
        maxwellian!(_M, u[:, j], prim[:, j])
    end

    return nothing

end

function mixture_maxwellian!(
    M::AbstractArray{T1,3},
    u::T2,
    v::T2,
    prim::AbstractMatrix{T3},
) where {T1<:Real,T2<:AbstractArray{<:Real,3},T3<:Real}

    for k in axes(M, 3)
        _M = @view M[:, :, k]
        maxwellian!(_M, u[:, :, k], v[:, :, k], prim[:, k])
    end

    return nothing

end

function mixture_maxwellian!(
    M::AbstractArray{T1,4},
    u::T2,
    v::T2,
    w::T2,
    prim::AbstractMatrix{T3},
) where {T1<:Real,T2<:AbstractArray{<:Real,4},T3<:Real}

    for l in axes(M, 4)
        _M = @view M[:, :, :, l]
        maxwellian!(_M, u[:, :, :, l], v[:, :, :, l], w[:, :, :, l], prim[:, l])
    end

    return nothing

end


text:
module MyExample

# Write your package code here.
greet() = print("Hello, World!")

include("extrafile.jl")

export my_f

end


text:
using Test
using ThreeBodyDecay

@testset "calculation of the binary decay" begin
    (p1,p2) = four_vectors_in_binary_decay(2rand()-1, π*(2rand()-1);
        m1sq = 5^2, m2sq = 1^2, m0sq = 8^2)
    @test sum(p1+p2 .≈ [0,0,0,8]) == 4
end

@testset "calculation of the binary decay in flight" begin
    m0sq = 4^2; m1sq = 1.2; m2sq = 1.8
    p0 = [1,2,3,sqrt(m0sq+14)]
    p1, p2 = four_vectors_in_binary_decay(p0,0.3,0.5; m1sq=m1sq, m2sq=m2sq)
    @test sum(p1 + p2 .≈ p0) == 4
    @test invmasssq(p1) ≈ m1sq && invmasssq(p2) ≈ m2sq
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
@testset "Map" begin
  include("Map/FqPolyRing.jl")
  include("Map/NumberField.jl")
  include("Map/NumField.jl")
  include("Map/AbGrp.jl")
end


text:
@testset "functiontransform" begin
    rng = MersenneTwister(123546)

    @testset "Real input" begin
        t = FunctionTransform(sin)

        x = randn(rng, 4)
        x′ = map(t, x)

        @test all([t(x[n]) == sin(x[n]) for n in eachindex(x)])
        @test all([t(x[n]) ≈ x′[n] for n in eachindex(x)])
    end

    @testset "Vector input" begin
        f = x -> sin.(x)
        t = FunctionTransform(f)

        x_vecs = [randn(rng, 5) for _ in 1:6]
        x_cols = ColVecs(randn(rng, 4, 7))
        x_rows = RowVecs(randn(rng, 6, 3))

        @testset "$(typeof(x))" for x in [x_vecs, x_cols, x_rows]
            x′ = map(t, x)
            @test all([t(x[n]) ≈ f(x[n]) for n in eachindex(x)])
            @test all([t(x[n]) ≈ x′[n] for n in eachindex(x)])
        end
    end

    @test repr(FunctionTransform(sin)) == "Function Transform: $(sin)"
    f(a, x) = sin.(a .* x)
    test_ADs(x->transform(SEKernel(), FunctionTransform(y->f(x, y))), randn(rng, 3))
end