text:
"""Functions for computing CIB models."""

abstract type AbstractCIBModel{T<:Real} <: AbstractForegroundModel end

"""
    CIB_Planck2013{T}(; kwargs...)

Define CIB model parameters. Defaults are from Viero et al. 2013. All numbers
not typed are converted to type T. This model has the following parameters and default values:

* `nside::Int64 = 4096`
* `hod::String = "shang"`
* `Inu_norm = 0.3180384`
* `min_redshift = 0.0`
* `max_redshift = 5.0`
* `min_mass = 1e12`
* `box_size = 40000`
* `shang_zplat = 2.0`
* `shang_Td = 20.7`
* `shang_betan = 1.6`
* `shang_eta = 2.4`
* `shang_alpha = 0.2`
* `shang_Mpeak = 10^12.3`
* `shang_sigmaM = 0.3`
* `shang_Msmin = 1e11`
* `shang_Mmin = 1e10`
* `shang_I0 = 46`
* `jiang_gamma_1 = 0.13`
* `jiang_alpha_1 = -0.83`
* `jiang_gamma_2 = 1.33`
* `jiang_alpha_2 = -0.02`
* `jiang_beta_2 = 5.67`
* `jiang_zeta = 1.19`

"""
@with_kw struct CIB_Planck2013{T<:Real} <: AbstractCIBModel{T} @deftype T
    nside::Int64    = 4096
    hod::String     = "shang"
    Inu_norm     = 0.3180384
    min_redshift = 0.0
    max_redshift = 5.0
    min_mass     = 1e12
    box_size     = 40000

    # shang HOD
    shang_zplat  = 2.0
    shang_Td     = 20.7
    shang_beta   = 1.6
    shang_eta    = 2.4
    shang_alpha  = 0.2
    shang_Mpeak  = 10^12.3
    shang_sigmaM = 0.3
    shang_Msmin  = 1e11
    shang_Mmin   = 1e10
    shang_I0     = 46

    # jiang
    jiang_gamma_1    = 0.13
    jiang_alpha_1    = -0.83
    jiang_gamma_2    = 1.33
    jiang_alpha_2    = -0.02
    jiang_beta_2     = 5.67
    jiang_zeta       = 1.19
end

function jiang_shmf(m, M_halo, model)
    dndm = (((model.jiang_gamma_1*((m/M_halo)^model.jiang_alpha_1))+
             (model.jiang_gamma_2*((m/M_halo)^model.jiang_alpha_2)))*
             (exp(-(model.jiang_beta_2)*((m/M_halo)^model.jiang_zeta))))
    return dndm
end

"""
Build a linear interpolation function which maps log(M_h) to N_sat.
"""
function build_shang_interpolator(
    min_log_M::T, max_log_M::T, model::AbstractCIBModel;
    n_bin=1000) where T

    x_m = LinRange(min_log_M, max_log_M, 1000)
    N_sat_i = zero(x_m)

    function integrand_m(lm, lM_halo)
        # lm is natural log of m
        m = exp(lm)
        M_halo = exp(lM_halo)
        dns_dm = jiang_shmf(m, M_halo, model)
        return dns_dm
    end

    for i in 1:size(x_m,1)
        N_sat_i[i], err = quadgk( t->integrand_m(t, x_m[i]),
                log(model.shang_Msmin), x_m[i], rtol=1e-8)
        N_sat_i[i] = convert(T, max(0.0, N_sat_i[i]))
    end

    return LinearInterpolation(x_m, N_sat_i)
end

function sigma_cen(m::T, model::AbstractCIBModel) where T
    return (exp( -(log10(m) - log10(model.shang_Mpeak))^2 /
        (T(2)*model.shang_sigmaM) ) * m) / sqrt(T(2π) * model.shang_sigmaM)
end

function nu2theta(nu::T, z::T, model::AbstractCIBModel) where T
    phys_h = T(6.62606957e-27)   # erg.s
    phys_k = T(1.3806488e-16)    # erg/K
    Td = model.shang_Td * (one(T)+z)^model.shang_alpha
    xnu = phys_h * nu / phys_k / Td
    return xnu^(T(4) + model.shang_beta) / expm1(xnu) / nu / model.shang_I0
end

"""
<L_sat> interpolation values
"""
function integrand_L(lm, lM_halo, model::AbstractCIBModel)
    m = exp(lm)
    return sigma_cen(m, model) * jiang_shmf(m, exp(lM_halo), model)
end

"""
Build a linear interpolator that takes in ln(M_halo) and returns sigma.
"""
function build_sigma_sat_ln_interpolator(
        max_ln_M::T, model::AbstractCIBModel; n_bin=1000) where T
    x = LinRange(log(model.shang_Mmin), max_ln_M, n_bin)
    L_mean = zero(x)

    for i in 1:n_bin
        L_mean[i], err = quadgk( t->integrand_L(t, x[i], model),
                log(model.shang_Mmin), x[i], rtol=1e-6)
    end
    return LinearInterpolation(T.(x), T.(L_mean), extrapolation_bc=zero(T))
end

function build_muofn_interpolator(model;
        min_mu::T=-6.0f0, max_mu::T=0.0f0, nbin=1000) where T
    mu = T(10) .^ LinRange(min_mu, max_mu, nbin)
    n  = zero(mu)
    for i in 1:nbin
        n[i], err = quadgk( lmu->jiang_shmf(exp(lmu), one(T), model),
                log(mu[i]), 0.0f0, rtol=1.0f-6)
    end
    return LinearInterpolation(T.(reverse(n)), T.(reverse(mu)))
end

"""
Compute redshift evolution factor for LF.
"""
function shang_z_evo(z::T, model::AbstractCIBModel) where T
    return (one(T) + min(z, model.shang_zplat))^model.shang_eta
end

"""
Construct the necessary interpolator set.
"""
function get_interpolators(model::AbstractCIBModel, cosmo::Cosmology.FlatLCDM{T},
    min_halo_mass::T, max_halo_mass::T) where T
    return (
        r2z = build_r2z_interpolator(
            model.min_redshift, model.max_redshift, cosmo),
        hod_shang = build_shang_interpolator(
            log(min_halo_mass), log(max_halo_mass), model),
        c_lnm2r = build_c_lnm2r_interpolator(),
        sigma_sat = build_sigma_sat_ln_interpolator(
            log(max_halo_mass), model),
        muofn = build_muofn_interpolator(model)
    )
end

# Fill up arrays with information related to CIB central sources.
function process_centrals!(
    model::AbstractCIBModel{T}, cosmo::Cosmology.FlatLCDM{T}, Healpix_res::Resolution;
    interp, hp_ind_cen, dist_cen, redshift_cen, theta_cen, phi_cen,
    lum_cen, n_sat_bar, n_sat_bar_result,
    halo_pos, halo_mass) where T

    N_halos = size(halo_mass, 1)

    Threads.@threads for i = 1:N_halos
        # location information for centrals
        hp_ind_cen[i] = Healpix.vec2pixRing(Healpix_res,
            halo_pos[1,i], halo_pos[2,i], halo_pos[3,i])
        theta_cen[i], phi_cen[i] = Healpix.vec2ang(halo_pos[1,i], halo_pos[2,i], halo_pos[3,i])
        dist_cen[i] = sqrt(halo_pos[1,i]^2 + halo_pos[2,i]^2 + halo_pos[3,i]^2)
        redshift_cen[i] = interp.r2z(dist_cen[i])

        # compute HOD
        n_sat_bar[i] = interp.hod_shang(log(halo_mass[i]))
        n_sat_bar_result[i] = rand(Distributions.Poisson(Float64.(n_sat_bar[i])))

        # get central luminosity
        lum_cen[i] = sigma_cen(halo_mass[i], model) * shang_z_evo(
            redshift_cen[i], model)
    end
end

# Fill up arrays with information related to CIB satellites.
function process_sats!(
        model::AbstractCIBModel{T}, cosmo::Cosmology.FlatLCDM{T},
        Healpix_res::Resolution;
        interp, hp_ind_sat, dist_sat, redshift_sat, theta_sat, phi_sat,
        lum_sat, cumsat,
        halo_mass, halo_pos, redshift_cen, n_sat_bar, n_sat_bar_result) where T

    N_halos = size(halo_mass, 1)
    Threads.@threads for i_halo = 1:N_halos
        r_cen = m2r(halo_mass[i_halo], cosmo)
        c_cen = mz2c(halo_mass[i_halo], redshift_cen[i_halo], cosmo)
        for j in 1:n_sat_bar_result[i_halo]
            # i is central halo index, j is index of satellite within each halo
            i_sat = cumsat[i_halo]+j # index of satellite in satellite arrays

            log_msat_inner = max(log(rand(T)), T(-7.0))
            r_sat = r_cen * interp.c_lnm2r(
                c_cen, log_msat_inner) * T(200.0^(-1.0/3.0))
            m_sat = (interp.muofn(rand(T) * n_sat_bar[i_halo])
                * halo_mass[i_halo])

            phi = random_phi(T)
            theta = random_theta(T)
            x_sat = halo_pos[1,i_halo] + r_sat * sin(theta) * cos(phi)
            y_sat = halo_pos[2,i_halo] + r_sat * sin(theta) * sin(phi)
            z_sat = halo_pos[3,i_halo] + r_sat * cos(theta)
            dist_sat[i_sat] = sqrt(x_sat^2 + y_sat^2 + z_sat^2)
            redshift_sat[i_sat] = interp.r2z(dist_sat[i_sat])
            theta_sat[i_sat], phi_sat[i_sat] = Healpix.vec2ang(x_sat, y_sat, z_sat)

            # lum_sat[i_sat] = interp.sigma_sat(log(m_sat)) * shang_z_evo(
            #     redshift_sat[i], model)
            lum_sat[i_sat] = sigma_cen(m_sat, model) * shang_z_evo(
                redshift_sat[i_sat], model)
            hp_ind_sat[i_sat] = Healpix.vec2pixRing(
                Healpix_res, x_sat, y_sat, z_sat)

        end
    end
end

"""
    generate_sources(model, cosmo, halo_pos_inp, halo_mass_inp; verbose=true)

Produce a source catalog from a model and halo catalog. This converts the
halo arrays into the type specified by `model`.

# Arguments:
- `model::AbstractCIBModel{T}`: source model parameters
- `cosmo::Cosmology.FlatLCDM{T}`: background cosmology
- `Healpix_res::Resolution`: Healpix map resolution
- `halo_pos_inp::AbstractArray{TH,2}`: halo positions with dims (3, nhalos)
- `halo_mass_inp::AbstractArray{TH,1}`: halo masses

# Keywords
- `verbose::Bool=true`: print out progress details
"""
function generate_sources(
        model::AbstractCIBModel{T}, cosmo::Cosmology.FlatLCDM{T},
        halo_pos_inp::AbstractArray{TH,2}, halo_mass_inp::AbstractArray{TH,1};
        verbose=true) where {T, TH}

    # make sure halo inputs are the CIB type
    halo_pos = convert(Array{T,2}, halo_pos_inp)
    halo_mass = convert(Array{T,1}, halo_mass_inp)

    # set up basics
    N_halos = size(halo_mass, 1)
    interp = get_interpolators( model, cosmo, minimum(halo_mass), maximum(halo_mass))
    res = Resolution(model.nside)

    verbose && println("Allocating for $(N_halos) centrals.")
    hp_ind_cen = Array{Int64}(undef, N_halos)  # healpix index of halo
    lum_cen = Array{T}(undef, N_halos)  # Lum of central w/o ν-dependence
    redshift_cen = Array{T}(undef, N_halos)
    theta_cen = Array{T}(undef, N_halos)
    phi_cen = Array{T}(undef, N_halos)
    dist_cen = Array{T}(undef, N_halos)
    n_sat_bar = Array{T}(undef, N_halos)
    n_sat_bar_result = Array{Int32}(undef, N_halos)

    # STEP 1: compute central properties -----------------------------------
    verbose && println("Processing centrals on $(Threads.nthreads()) threads.")
    process_centrals!(model, cosmo, res,
        interp=interp, hp_ind_cen=hp_ind_cen, dist_cen=dist_cen,
        redshift_cen=redshift_cen, theta_cen=theta_cen, phi_cen=phi_cen,
        lum_cen=lum_cen, n_sat_bar=n_sat_bar,
        n_sat_bar_result=n_sat_bar_result,
        halo_pos=halo_pos, halo_mass=halo_mass)

    # STEP 2: Generate satellite arrays -----------------------------
    cumsat = generate_subhalo_offsets(n_sat_bar_result)
    total_n_sat = cumsat[end]
    hp_ind_sat = Array{Int64}(undef, total_n_sat)  # healpix index of halo
    lum_sat = Array{T}(undef, total_n_sat)  # Lum of central w/o ν-dependence
    redshift_sat = Array{T}(undef, total_n_sat)
    theta_sat = Array{T}(undef, total_n_sat)
    phi_sat = Array{T}(undef, total_n_sat)
    dist_sat = Array{T}(undef, total_n_sat)

    # STEP 3: compute satellite properties -----------------------------------
    verbose && println("Processing $(total_n_sat) satellites.")
    process_sats!(model, cosmo, res,
        interp=interp, hp_ind_sat=hp_ind_sat, dist_sat=dist_sat,
        redshift_sat=redshift_sat, theta_sat=theta_sat, phi_sat=phi_sat,
        lum_sat=lum_sat, cumsat=cumsat,
        halo_mass=halo_mass, halo_pos=halo_pos, redshift_cen=redshift_cen,
        n_sat_bar=n_sat_bar, n_sat_bar_result=n_sat_bar_result)

    return (
        hp_ind_cen=hp_ind_cen, lum_cen=lum_cen,
        redshift_cen=redshift_cen, theta_cen=theta_cen, phi_cen=phi_cen, dist_cen=dist_cen,
        hp_ind_sat=hp_ind_sat, lum_sat=lum_sat,
        redshift_sat=redshift_sat, theta_sat=theta_sat, phi_sat=phi_sat, dist_sat=dist_sat,
        N_cen=N_halos, N_sat=total_n_sat
    )
end

"""
    paint!(result_map, nu_obs, model, sources, fluxes_cen, fluxes_sat)

Paint a source catalog onto a map, recording the fluxes in
`fluxes_cen` and `fluxes_sat`.

# Arguments:
- `result_map::HealpixMap{T_map, RingOrder}`: Healpix map to paint
- `nu_obs`: frequency in Hz
- `model::AbstractCIBModel{T}`: source model parameters
- `sources`: NamedTuple containing source information from generate_sources
- `fluxes_cen::AbstractArray`: buffer for writing fluxes of centrals
- `fluxes_sat::AbstractArray`: buffer for writing fluxes of satellites
"""
function paint!(result_map::HealpixMap{T_map, RingOrder},
        nu_obs, model::AbstractCIBModel{T}, sources,
        fluxes_cen::AbstractArray, fluxes_sat::AbstractArray) where {T_map, T}

    pixel_array = result_map.pixels
    fill!(pixel_array, zero(T))  # prepare the frequency map

    # process centrals for this frequency
    Threads.@threads for i in 1:sources.N_cen
        nu = (one(T) + sources.redshift_cen[i]) * nu_obs
        fluxes_cen[i] = l2f(
            sources.lum_cen[i] * nu2theta(
                nu, sources.redshift_cen[i], model),
            sources.dist_cen[i], sources.redshift_cen[i])
        pixel_array[sources.hp_ind_cen[i]] += fluxes_cen[i]
    end

    # process satellites for this frequency
    Threads.@threads for i in 1:sources.N_sat
        nu = (one(T) + sources.redshift_sat[i]) * nu_obs
        fluxes_sat[i] = l2f(
            sources.lum_sat[i] * nu2theta(
                nu, sources.redshift_sat[i], model),
            sources.dist_sat[i], sources.redshift_sat[i])
        pixel_array[sources.hp_ind_sat[i]] += fluxes_sat[i]
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(result_map.resolution.nside)
    pixel_array .*= per_pixel_steradian
end

"""
Paint a source catalog onto a map.

This function creates the arrays for you.
"""
function paint!(result_map::HealpixMap{T,RingOrder},
        nu_obs::T, model::AbstractCIBModel, sources) where T

    fluxes_cen = Array{T, 1}(undef, sources.N_cen)
    fluxes_sat = Array{T, 1}(undef, sources.N_sat)

    paint!(result_map, nu_obs, model, sources,
        fluxes_cen, fluxes_sat)

    return fluxes_cen, fluxes_sat
end

max_stars_repo_name: xzackli/XGPaint

max_issues_repo_name: xzackli/XGPaint

max_forks_repo_name: xzackli/XGPaint


text:
module Interpreter

using ..DataFlow
using ..DataFlow: Call, Constant, constant, value, spliceinputs, Frame, Line,
  Lambda, Split, inputs
using MacroTools: @q

export Context, mux, stack, interpret,
  interpv, ituple, ilambda, iconst, iline, ilinev, iargs,
  @icatch, @ithrow

mux(f) = f
mux(m, f) = (xs...) -> m(f, xs...)
mux(ms...) = foldr(mux, ms)

struct Context{T}
  interp::T
  cache::IdDict{Any, Any}
  stack::Vector{Any}
  data::Dict{Symbol,Any}
end

Context(interp; kws...) = Context(interp, IdDict{Any, Any}(), [], Dict{Symbol,Any}(kws))

Base.getindex(ctx::Context, k::Symbol) = ctx.data[k]
Base.setindex!(ctx::Context, v, k::Symbol) = ctx.data[k] = v

function stack(c::Context)
  stk = []
  isempty(c.stack) && return stk
  frame = nothing
  for i = 1:length(c.stack)
    isa(c.stack[i], Frame) || continue
    i > 1 && isa(c.stack[i-1], Line) && pushfirst!(stk, (frame, c.stack[i-1]))
    frame = c.stack[i].f
  end
  isa(c.stack[end], Line) && pushfirst!(stk, (frame, c.stack[end]))
  return stk
end

function interpv(ctx::Context, graph::IVertex)
  haskey(ctx.cache, graph) && return ctx.cache[graph]
  ctx.cache[graph] = ctx.interp(ctx, value(graph), inputs(graph)...)
end

function interpret(ctx::Context, graph::IVertex, args::IVertex...)
  graph = spliceinputs(graph, args...)
  interpv(ctx, graph)
end

interpret(ctx::Context, graph::IVertex, args...) =
  interpret(ctx, graph, map(constant, args)...)

# The `ifoo` convention denotes a piece of interpreter middleware

iconst(f, ctx::Context, c::Constant) = c.value

function iline(f, ctx::Context, l::Union{Line,Frame}, v)
  push!(ctx.stack, l)
  val = interpv(ctx, v)
  pop!(ctx.stack)
  return val
end

ilinev(f, ctx::Context, l::Union{Line,Frame}, v) = vertex(l, iline(f, ctx, l, v))

ilambda(f, ctx::Context, λ::Lambda, vars...) =
  (xs...) -> interpret(ctx, λ.body, vars..., xs...)

iargs(cb, ctx::Context, f, xs...) = cb(ctx, f, interpv.(ctx, xs)...)

function ituple(f, ctx::Context, s::Split, xs)
  isa(xs, Vertex) && value(xs) == tuple ? inputs(xs)[s.n] :
  isa(xs, Tuple) ? xs[s.n] :
    f(ctx, s, xs)
end

ituple(f, ctx::Context, ::Call, ::typeof(tuple), xs...) = tuple(xs...)

for m in :[iconst, iline, ilinev, ilambda, ituple].args
  @eval $m(f, args...) = f(args...)
end

interpeval = mux(iline, ilambda, iconst, iargs,
                 ituple, (ctx, f, xs...) -> f(xs...))

interpret(graph::IVertex, args...) =
  interpret(Context(interpeval), graph, args...)

# Error Handling

import Juno: errmsg, errtrace

framename(f::Function) = typeof(f).name.mt.name
framename(f::Base.Nothing) = Symbol("<none>")
framename(x) = Symbol(string(typeof(x)))

totrace(stack) = [StackTraces.StackFrame(framename(f), Symbol(line.file), line.line)
                  for (f, line) in stack]

Base.stacktrace(c::Context) = totrace(stack(c))

struct Exception{T}
  err::T
  trace::StackTraces.StackTrace
end

errmsg(e::Exception) = errmsg(e.err)
errtrace(e::Exception, bt) = errtrace(e.err, [e.trace..., bt...])

function Base.showerror(io::IO, e::Exception)
  showerror(io, e.err)
  for t in e.trace
    println(io)
    show(io, t)
  end
end

struct InterpError
  e
end

function trimtrace(trace)
  trace = reverse(trace)
  comp = reverse(StackTraces.stacktrace())
  n = 1
  while n < length(trace) && n < length(comp) && trace[n] == comp[n]
    n += 1
  end
  reverse(trace[n:end])
end

macro icatch(ctx, ex)
  :(try
      $(esc(ex))
    catch e
      isa(e, InterpError) && rethrow()
      ctx = $(esc(ctx))
      trace = vcat(trimtrace(catch_stacktrace()), stacktrace(ctx))
      throw(InterpError(Exception(e, trace)))
    end)
end

macro ithrow(ex)
  @q try
    $(esc(ex))
  catch e
    isa(e, InterpError) || rethrow()
    throw(e.e)
  end
end

end

ext: jl


text:
#=
Simulation:
- Julia version:
- Author: inaba
- Date: 2021-11-20
=#
module Simulation
using Dates: format, now
using LightGraphs
using Random
using Statistics: mean

mutable struct Agent
    id::Int
    is_cooperator::Bool
    is_punisher::Bool
    next_is_cooperator::Bool
    next_is_punisher::Bool
    payoff::Float64

    Agent(id::Int) = new(id, rand([true, false]), rand([true, false]), false, false, 0.0)
end

struct Model
    agents::Vector{Agent}
    graph::SimpleGraph{Int64}
    hop::Int  # neighbor parameter
    neighbours::Vector{Vector{Int64}}  # neighbors of every nodes

    g::Float64  # game interaction locality parameter
    p::Float64  # punishment locality parameter
    a::Float64  # adaptation (learning) locality parameter

    b::Float64  # benefit of a token
    c::Float64  # game cost
    f::Float64  # fine
    s::Float64  # sanction for punishment

    μ::Float64  # mutation rate

    Model(agents::Vector{Agent}, graph::SimpleGraph; g::Float64, p::Float64, a::Float64) = new(
        agents, graph, 2,
        [[_n for _n in neighborhood(graph, i, 2) if _n != i] for i in vertices(graph)],
        g, p, a,
        2.0, 1.0, 6.0, 3.0, 0.01
    )
end

cooperator_rate(model::Model)::Float64 = mean([agent.is_cooperator for agent in model.agents])
punisher_rate(model::Model)::Float64 = mean([agent.is_punisher for agent in model.agents])

function calc_payoffs!(model::Model)
    n_size = length(model.agents)
    global_cooperator_rate = cooperator_rate(model)
    global_punisher_rate = punisher_rate(model)

    for agent in model.agents
        local_cooperator_rate = mean([model.agents[neighbor].is_cooperator for neighbor in model.neighbours[agent.id]])
        local_punisher_rate = mean([model.agents[neighbor].is_punisher for neighbor in model.neighbours[agent.id]])

        _global_cooperator_rate = (global_cooperator_rate * n_size - Int(agent.is_cooperator)) / (n_size - 1)
        _global_punisher_rate = (global_punisher_rate * n_size - Int(agent.is_punisher)) / (n_size - 1)

        # ゲームの中で協力者が占める割合
        _cooperator_rate = local_cooperator_rate * (1. - model.g) + _global_cooperator_rate * model.g

        # benefit (ゲームの中で協力者が占める割合 (自分を含む))
        _b = (_cooperator_rate * model.hop * 2 + Int(agent.is_cooperator)) / (model.hop * 2 + 1)
        # game cost
        _c = agent.is_cooperator ? 1 : 0
        # sanction (罰を課すコスト)
        _s = agent.is_punisher ? (1 - local_cooperator_rate) * (1 - model.p) + (1 - _global_cooperator_rate) * model.p : 0.0
        # Fine (罰金)
        _f = agent.is_cooperator ? 0.0 : local_punisher_rate * (1 - model.p) + _global_punisher_rate * model.p

        agent.payoff = _b * model.b - _c * model.c - _s * model.s - _f * model.f
    end
end

function set_next_strategies!(model::Model)
    for agent in model.agents
        if rand() < model.μ
            agent.next_is_cooperator = rand([true, false])
            agent.next_is_punisher = rand([true, false])
        else
            y = if rand() < 0.5 # a
                rand([_a.id for _a in model.agents if _a.id != agent.id])
            else
                rand(model.neighbours[agent.id])
            end

            if model.agents[y].payoff > agent.payoff
                agent.next_is_cooperator = model.agents[y].is_cooperator
                agent.next_is_punisher = model.agents[y].is_punisher
            else
                agent.next_is_cooperator = agent.is_cooperator
                agent.next_is_punisher = agent.is_punisher
            end
        end
    end
end

function update_strategies!(model::Model)
    for agent in model.agents
        agent.is_cooperator = agent.next_is_cooperator
        agent.is_punisher = agent.next_is_punisher
    end
end

function run(;
    n_size::Int = 100,
    generation::Int = 100,
    trial::Int = 100,
    gs::Vector = [0.0, 0.25, 0.5, 0.75, 1.0],
    ps::Vector = [0.0, 0.25, 0.5, 0.75, 1.0],
    as::Vector = [0.0, 0.25, 0.5, 0.75, 1.0],
    file_name = "out/jl$(format(now(), "yyyymmdd_HHMMSS")).csv")

    for g = gs, p = ps, a = as
        cooperator_rate_vec = []
        punisher_rate_vec = []
        for t in 1:trial
            agents = [Agent(id) for id in 1:n_size]
            graph = cycle_graph(n_size)
            model = Model(agents, graph; g, p, a)

            for step = 1:generation
                calc_payoffs!(model)
                set_next_strategies!(model)
                update_strategies!(model)
            end

            push!(cooperator_rate_vec, cooperator_rate(model))
            push!(punisher_rate_vec, punisher_rate(model))
        end

        csv = "$(g),$(p),$(a),$(mean(cooperator_rate_vec)),$(mean(punisher_rate_vec))"
        println(csv)
        open(file_name, "a") do out
            println(out, csv)
        end
    end
end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Simulation
    Simulation.run()
    # Simulation.run(gs = [0.0], ps = [1.0], as = [0.0])
end


text:
## subs
"""
Substitute values into a symbolic expression.

Examples
```
@syms x y
ex = x^2 + y^2
subs(ex, x, 1) # 1 + y^2
subs(ex, (x, 1)) # 1 + y^2
subs(ex, (x, 1), (y,x)) # 1 + x^2, values are substituted left to right.
subs(ex, x=>1)  # alternate to subs(x, (x,1))
subs(ex, x=>1, y=>1) # ditto
```
"""
function subs{T<:SymbolicType, S<:SymbolicType}(ex::T, var::S, val)
    s = Basic()
    ccall((:basic_subs2, libsymengine), Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Basic}, Ptr{Basic}), &s, &ex, &var, &val)
    return s
end
function subs{T<:SymbolicType}(ex::T, d::CMapBasicBasic)
    s = Basic()
    ccall((:basic_subs, libsymengine), Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Void}), &s, &ex, d.ptr)
    return s
end

subs{T<:SymbolicType}(ex::T, d::Dict) = subs(ex, CMapBasicBasic(d))
subs{T <: SymbolicType, S<:SymbolicType}(ex::T, y::Tuple{S, Any}) = subs(ex, y[1], y[2])
subs{T <: SymbolicType, S<:SymbolicType}(ex::T, y::Tuple{S, Any}, args...) = subs(subs(ex, y), args...)
subs{T <: SymbolicType}(ex::T, d::Pair...) = subs(ex, [(p.first, p.second) for p in d]...)

## Allow an expression to be called, as with ex(2). When there is more than one symbol, one can rely on order of `free_symbols` or
## be explicit by passing in pairs : `ex(x=>1, y=>2)` or a dict `ex(Dict(x=>1, y=>2))`.
## This uses `eval` to avoid having to work around different styles in v0.5 and v0.4
call_v0_4 = quote
  function Base.call{T <: Basic}(ex::T, args...)
      xs = free_symbols(ex)
      subs(ex, collect(zip(xs, args))...)
  end
  Base.call(ex::Basic, x::Dict) = subs(ex, x)
  Base.call(ex::Basic, x::Pair...) = subs(ex, x...)
end

call_v0_5 = quote
  function (ex::Basic)(args...)
      xs = free_symbols(ex)
      subs(ex, collect(zip(xs, args))...)
  end
  (ex::Basic)(x::Dict) = subs(ex, x)
  (ex::Basic)(x::Pair...) = subs(ex, x...)
end
VERSION < v"0.5.0" ? eval(call_v0_4) : eval(call_v0_5)

## Lambdify

## Map symengine classes to function names
fn_map = Dict(
              :Add => :+,
              :Sub => :-,
              :Mul => :*,
              :Div => :/,
              :Pow => :^,
              :re  => :real,
              :im  => :imag,
              :Abs => :abs # not really needed as default in map_fn covers this
              )

map_fn(key, fn_map) = haskey(fn_map, key) ? fn_map[key] : Symbol(lowercase(string(key)))

function walk_expression(ex)
    fn = get_symengine_class(ex)

    if fn == :Symbol
        return Symbol(toString(ex))
    elseif (fn in number_types) || (fn == :Constant)
        return N(ex)
    end

    as = get_args(ex)

    Expr(:call, map_fn(fn, fn_map), [walk_expression(a) for a in as]...)
end

## evaluate symbolless expression or return afunction
function lambdify(ex::Basic)
    vars = free_symbols(ex)
    if length(vars) == 0
        _lambdify(ex)
    else
        _lambdify(ex, vars)
    end
end
lambdify(ex::BasicType) = lambdify(Basic(ex))

## return a number
function _lambdify(ex)
    body = walk_expression(ex)
    eval(body)
end

## return a function
function _lambdify(ex::Basic, vars)
    body = walk_expression(ex)

    try
        fn = eval(Expr(:function,
                  Expr(:call, gensym(), map(Symbol,vars)...),
                       body))
        (args...) -> invokelatest(fn, args...) # https://github.com/JuliaLang/julia/pull/19784
    catch err
        throw(ArgumentError("Expression does not lambdify"))
    end
end


text:
"""
    DelayFilm

The delay film stores a 1D-histogram of observed intensity as a function of depth in the
scene. This can then be post-processed into e.g. reflectivity as a function of time delay.

"""
struct DelayFilm <: Film
    wavelength::Float64
    area::Float64
    resX::Int64
    resY::Int64
    nbins::Int64
    tmin::Float64
    tmax::Float64
    histogram::Array{Float64,1}
    counts::Array{Int64,1}
end

function DelayFilm(lmd::Real, window::Array{Float64,1}, resX::Integer, resY::Integer, nbins::Integer, tmin::Real, tmax::Real)
    area = (window[2] - window[1]) * (window[4] - window[3])
    DelayFilm(lmd, area, resX, resY, nbins, tmin, tmax, zeros(Float64, nbins), zeros(Int64,nbins+1))
end

"""
    find_bin(F::DelayFilm, d::Real)

For a given distance `d`, find the corresponding bin in the `DelayFilm` histogram.

"""
function find_bin(F::DelayFilm, d::Real)
    nd = (d - F.tmin) / (F.tmax - F.tmin)
    return ceil(Int64, nd * F.nbins)
end

"""
    add_sample!(F::DelayFilm, sample::Sample, L::Spectrum, isect::Union{Intersection,Nothing})

Add a ray sample to the `DelayFilm`. The ray intersection `isect` is used to find the depth
in the scene where the reflection happened. The spectral intensity `L` will be added to the
bin corresponding to that distance. The camera sample `sample` is not used, but it's part
of the `Film` API.

"""
function add_sample!(F::DelayFilm, sample::Sample, L::Spectrum, isect::Union{Intersection,Nothing})
    F.counts[end] += 1
    if isect == nothing
        return nothing
    end
    if isect.tmin > F.tmax || isect.tmin < F.tmin
        return nothing
    end
    idx = find_bin(F, isect.tmin)
    F.counts[idx] += 1

    if isblack(L)
        return nothing
    end
    I = interpolate(L, F.wavelength)
    F.histogram[idx] += I
    return nothing
end

"""
    reset!(F::DelayFilm)

Reset all the pixel values of the `DelayFilm` to zero. This is mainly useful if you want to
use the same film in multiple computations, writing out the data and resetting in between.

"""
function reset!(F::DelayFilm)
    F.counts .= 0
    F.histogram .= 0.0
    nothing
end

"""
    write_txt(F::DelayFilm, fname::String)

Write the `DelayFilm` histogram into a text file. The file starts with seven lines of
metadata, followed by the values in the histogram vector.

# Metadata header example

```
# Distance and delay
# lambda = 532.0
# tmin = 0.0
# tmax = 4.0
# nbins = 1000
# step = 0.004
# counts = 2359296
```

"""
function write_txt(F::DelayFilm, fname::String)
    f = open(fname, "w")
    write(f, "# Distance and delay\n")
    write(f, "# lambda = $(F.wavelength)\n")
    write(f, "# tmin = $(F.tmin)\n")
    write(f, "# tmax = $(F.tmax)\n")
    write(f, "# nbins = $(F.nbins)\n")
    write(f, "# step = $((F.tmax - F.tmin) / F.nbins)\n")
    write(f, "# counts = $(F.counts[end])\n")
    sum_counts = F.counts[end]
    for i = 1:F.nbins
        area_in_bin = F.area * F.counts[i] / sum_counts
        value = F.counts[i] > 0 ? area_in_bin * F.histogram[i] / F.counts[i] : 0.0
        write(f, string(value))
        write(f, "\n")
    end
    close(f)
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using Documenter
using KernelFunctions

makedocs(
    sitename = "KernelFunctions",
    format = Documenter.HTML(),
    modules = [KernelFunctions],
    pages = ["Home"=>"index.md",
             "User Guide" => "userguide.md",
             "Examples"=>"example.md",
             "Kernel Functions"=>"kernels.md",
             "Transform"=>"transform.md",
             "Metrics"=>"metrics.md",
             "Theory"=>"theory.md",
             "API"=>"api.md"]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/theogf/KernelFunctions.jl.git",
    target = "build"
)


text:
using Plots
using MeshCat
using GeometryTypes
using CoordinateTransformations
using FileIO
using MeshIO
using Random

Random.seed!(123)

# Solver options
N = 101 # 201
integration = :rk4
opts = SolverOptions()
opts.verbose = false
opts.square_root = true
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-4
opts.outer_loop_update_type = :feedback

dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

# Obstacle Avoidance
model,obj_uncon = TrajectoryOptimization.Dynamics.quadrotor
r_quad = 3.0
n = model.n
m = model.m
obj_con = TrajectoryOptimization.Dynamics.quadrotor_3obs[2]
spheres = TrajectoryOptimization.Dynamics.quadrotor_3obs[3]
n_spheres = length(spheres)

solver_uncon = Solver(model,obj_uncon,integration=integration,N=N,opts=opts)
solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

U_hover = 0.5*9.81/4.0*ones(solver_uncon.model.m, solver_uncon.N-1)
X_hover = rollout(solver_uncon,U_hover)
@time results_uncon, stats_uncon = solve(solver_uncon,U_hover)
# @time results_uncon_dircol, stats_uncon_dircol = TrajectoryOptimization.solve_dircol(solver_uncon, X_hover, U_hover, options=dircol_options)

@time results_con, stats_con = solve(solver_con,U_hover)
# @time results_dircol, stats_dircol = TrajectoryOptimization.solve_dircol(solver_con, X_hover, U_hover, options=dircol_options)

# Trajectory plots
t_array = range(0,stop=solver_con.obj.tf,length=solver_con.N)
plot(t_array[1:end-1],to_array(results_uncon.U)',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="")
plot(t_array,to_array(results_con.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)

# plot(t_array[1:end-1],Array(results_uncon_dircol.U)[:,1:end-1]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="dircol")
# plot(t_array,Array(results_uncon_dircol.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="dircol")

t_array = range(0,stop=solver_con.obj.tf,length=solver_con.N)
plot(t_array[1:end-1],Array(results_dircol.U)[:,1:end-1]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="control",labels="dircol")
plot(t_array,Array(results_dircol.X)[1:3,:]',title="Quadrotor Obstacle Avoidance",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)
@assert max_violation(results_con) <= opts.constraint_tolerance

# Constraint convergence plot
plot(log.(stats_con["c_max"]),title="Quadrotor Obstacle Avoidance",xlabel="iteration",ylabel="log(max constraint violation)",label="")

@show stats_con["iterations"]
@show stats_con["outer loop iterations"]
@show stats_con["c_max"][end]
@show stats_con["cost"][end]


text:
# Tests using TaylorModel1 and RTaylorModel1

using TaylorModels
using LinearAlgebra: norm
using Test

const _num_tests = 1000
const α_mid = TaylorModels.α_mid

setformat(:full)

function check_containment(ftest, tma::T) where {T<:Union{TaylorModel1, RTaylorModel1}}
    x0 = expansion_point(tma)
    xfp = diam(tma.dom)*(rand()-0.5) + mid(x0)
    xbf = big(xfp)
    range = tma((xfp .. xfp)-x0)
    bb = ftest(xbf) ∈ range
    bb || @show(ftest, xfp, xbf, ftest(xbf), range)
    return bb
end

@testset "Test `bound_taylor1`" begin
    x0 = Interval(0.0)
    ii0 = Interval(-0.5, 0.5)

    tpol = exp( Taylor1(2) )
    @test TaylorModels.bound_taylor1( tpol, ii0) == tpol(ii0.lo) .. tpol(ii0.hi)
    @test TaylorModels.bound_taylor1( exp( Taylor1(Interval{Float64}, 2) ),
        ii0) == tpol(ii0.lo) .. tpol(ii0.hi)

    # An uncomfortable example from Makino
    t = Taylor1(5)
    f(x) = 1 - x^4 + x^5
    @test interval(1-4^4/5^5,1) ⊆ TaylorModels.bound_taylor1(f(t), 0..1)
    tm = TaylorModel1(5, x0, ii0)
    @test interval(1-4^4/5^5,1) ⊆ TaylorModels.bound_taylor1(f(tm))
    @test interval(1-4^4/5^5,1) ⊆ TaylorModels.bound_taylor1(f(tm), 0..1)
end

@testset "Tests for TaylorModel1 " begin
    x0 = Interval(0.0)
    ii0 = Interval(-0.5, 0.5)
    x1 = Interval(1.0)
    ii1 = Interval(0.5, 1.5)

    @testset "TaylorModel1 constructors" begin
        tv = TaylorModel1{Interval{Float64},Float64}(Taylor1(Interval{Float64},5), x0, x0, ii0)
        @test tv == TaylorModel1(Taylor1(Interval{Float64},5), x0, x0, ii0)
        @test tv == TaylorModel1(5, x0, ii0)
        @test tv == TaylorModel1(5, ii0)
        @test tv == TaylorModel1(5, 0.0, ii0)
        @test TaylorModel1(x1, 5, x0, ii0) == TaylorModel1(Taylor1(x1, 5), x0, x0, ii0)
        @test TaylorModel1(5, 0.7, ii1) == TaylorModel1(5, interval(0.7), ii1)

        @test TaylorModel1(tv, ii0) == TaylorModel1(Taylor1(Interval{Float64}, 5), ii0, x0, ii0)
        @test TaylorModel1(5, x0, ii0) == TaylorModel1(tv, x0)
        @test TaylorModel1(5, ii0) == TaylorModel1(tv, x0)

        @test isa(tv, AbstractSeries)
        @test TaylorModel1{Interval{Float64},Float64} <: AbstractSeries{Interval{Float64}}

        # Test errors in construction
        @test_throws AssertionError TaylorModel1(Taylor1(Interval{Float64},5), x1, x0, ii0)
        @test_throws AssertionError TaylorModel1(5, x1, ii0)

        # Tests for get_order, remainder, polynomial and domain
        @test get_order(tv) == 5
        @test remainder(tv) == interval(0.0)
        @test polynomial(tv) == Taylor1(Interval{Float64},5)
        @test domain(tv) == ii0
        @test expansion_point(tv) == x0
        @test constant_term(tv) == interval(0.0)
        @test linear_polynomial(tv) == Taylor1(Interval{Float64},5)

        # Tests related to fixorder
        a = TaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        b = TaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
        aa, bb = TaylorModels.fixorder(a, b)
        @test get_order(aa) == get_order(bb) == 1
        @test isa(aa, TaylorModel1) == isa(bb, TaylorModel1) == true
        @test aa == a
        @test bb == TaylorModel1(Taylor1([1.0, 1, 0]), -1 .. 2, 0..0, -1 .. 1)
        # a and b remain the same
        @test a == TaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        @test b == TaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
    end

    @testset "Arithmetic operations" begin
        Δ = interval(-0.25, 0.25)
        a = TaylorModel1(x1+Taylor1(5), Δ, x1, ii1)
        tv = TaylorModel1(5, x1, ii1)
        @test zero(a) == TaylorModel1(zero(a.pol), 0..0, x1, ii1)
        @test one(a) == TaylorModel1(one(a.pol), 0..0, x1, ii1)
        @test a+x1 == TaylorModel1(2*x1+Taylor1(5), Δ, x1, ii1)
        @test a+a == TaylorModel1(2*(x1+Taylor1(5)), 2*Δ, x1, ii1)
        @test a-x1 == TaylorModel1(zero(x1)+Taylor1(5), Δ, x1, ii1)
        @test a-a == TaylorModel1(zero(a.pol), 2*Δ, x1, ii1)
        b = a * tv
        @test b == TaylorModel1(a.pol*tv.pol, a.rem*tv.pol(ii1-x1), x1, ii1)
        @test remainder(b/tv) ⊆ Interval(-0.78125, 0.84375)
        b = a * a.pol[0]
        @test b == a
        @test constant_term(a) == x1
        @test linear_polynomial(a) == Taylor1(5)

        a = TaylorModel1(x0, 5, x0, ii0)
        @test a^0 == TaylorModel1(x0^0, 5, x0, ii0)
        @test a^1 == TaylorModel1(x0^1, 5, x0, ii0)
        @test a^2 == TaylorModel1(x0^2, 5, x0, ii0)
        @test a^3 == TaylorModel1(x0^3, 5, x0, ii0)
        a = TaylorModel1(x1, 5, x1, ii1)
        @test a^0 == TaylorModel1(x1^0, 5, x1, ii1)
        @test a^1 == TaylorModel1(x1^1, 5, x1, ii1)
        @test a^2 == TaylorModel1(x1^2, 5, x1, ii1)
        @test a^3 == TaylorModel1(x1^3, 5, x1, ii1)

        # Tests involving TM1s with different orders
        a = TaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        b = TaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
        aa, bb = TaylorModels.fixorder(a, b)
        @test get_order(aa) == get_order(bb)
        @test get_order(bb) == 1
        @test a + b == aa + bb
        @test a - b == aa - bb
        res1 = a * b
        res2 = aa * bb
        @test res1 == TaylorModel1(Taylor1([1.0, 2]), -2 .. 9 , 0..0, -1 .. 1)
        @test res2 == TaylorModel1(Taylor1([1.0, 2]), -3 .. 9 , 0..0, -1 .. 1)
        res1 = a / b
        res2 = aa / bb
        @test res1 == TaylorModel1(Taylor1([1.0, 0]), entireinterval() , 0..0, -1 .. 1)
        @test res2 == res1
        # a and b remain the same
        @test a == TaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        @test b == TaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)

        @test_throws AssertionError a+TaylorModel1(a.pol, a.rem, 1..1, -1..1)
        @test_throws AssertionError a+TaylorModel1(a.pol, a.rem, 0..0, -2..2)
    end

    @testset "TM1's with TaylorN coefficients" begin
        # Tests for TM1's with TaylorN coefficients
        orderT = 4
        orderQ = 5
        ξ = set_variables("ξ", order = 2 * orderQ, numvars=1)
        q0 = [0.5]
        δq0 = IntervalBox(-0.1 .. 0.1, Val(1))
        qaux = normalize_taylor(q0[1] + TaylorN(1, order=orderQ),
                                δq0, true)
        symIbox = IntervalBox(-1 .. 1, Val(1))
        t = Taylor1([qaux, 1], orderT)
        dom = 0 .. 1
        x00 = mid(dom)
        tm = TaylorModel1(deepcopy(t), 0 .. 0, x00, dom)

        f(x) = x + x^2
        g(x) = x

        fT = f(tm)
        gT = g(tm)

        for ind = 1:_num_tests
            fgTM1 = fT * gT
            xξ = rand(domain(fgTM1))
            q0ξ = (q0 .+ rand(δq0))[1]
            t = Taylor1(orderT) + q0ξ
            fgT1 = f(t) * g(t)
            @test fgT1(xξ - q0ξ) ∈ fgTM1(xξ-fgTM1.x0)(symIbox)
        end

    end

    @testset "RPAs, functions and remainders" begin
        @test rpa(x->5+zero(x), TaylorModel1(4, x0, ii0)) ==
            TaylorModel1(interval(5.0), 4, x0, ii0)
        @test rpa(x->5+one(x), TaylorModel1(4, x1, ii1)) ==
            TaylorModel1(5+x1, 4, x1, ii1)
        @test rpa(x->5*x, TaylorModel1(4, x1, ii1)) ==
            TaylorModel1(5.0*(x1+Taylor1(4)), x0, x1, ii1)
        @test rpa(x->5*x^0, TaylorModel1(4, x0, ii0)) == 5*TaylorModel1(4, x0, ii0)^0
        @test rpa(x->5*x^0, TaylorModel1(4, x0, ii0)) ==
            TaylorModel1( interval(5.0)*Taylor1(4)^0, x0, x0, ii0)
        @test rpa(x->5*x^1, TaylorModel1(4, x0, ii0)) == 5*TaylorModel1(4, x0, ii0)^1
        @test rpa(x->5*x^1, TaylorModel1(4, x0, ii0)) ==
            TaylorModel1( interval(5.0)*Taylor1(4)^1, x0, x0, ii0)
        @test rpa(x->5*x^2, TaylorModel1(4, x0, ii0)) == 5*TaylorModel1(4, x0, ii0)^2
        @test rpa(x->5*x^2, TaylorModel1(4, x0, ii0)) ==
            TaylorModel1( interval(5.0)*Taylor1(4)^2, x0, x0, ii0)
        @test rpa(x->5*x^4, TaylorModel1(4, x0, ii0)) ==
            TaylorModel1( interval(5.0)*Taylor1(4)^4, x0, x0, ii0)
        @test rpa(x->5*x^4, TaylorModel1(3, x0, ii0)) ==
            TaylorModel1( Taylor1(x0, 3), 5*ii0^4, x0, ii0)

        # Testing remainders of an RPA
        order = 2
        ii = ii0
        xx = x0
        ftest = x -> exp(x)-1
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        # test for TM with scalar coefficients
        @test fp_rpa(tmc) == tmc

        order = 2
        ii = ii1
        xx = x1
        ftest = x -> exp(x)
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 3
        ii = ii0
        xx = x0
        ftest = x -> sin(x)
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 2
        ii = ii1
        xx = x1
        ftest = x -> sqrt(x)
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 5
        ii = ii1
        xx = x1
        ftest = x -> inv(x)
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.hi)-tmc.pol(ii.hi-ξ0),
                        ftest(ii.lo)-tmc.pol(ii.lo-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        # Example of Makino's thesis (page 98 and fig 4.2)
        order = 8
        ii = interval(-0.5, 1.0)
        xx = interval(mid(ii))
        ftest = x -> x*(x-1.1)*(x+2)*(x+2.2)*(x+2.5)*(x+3)*sin(1.7*x+0.5)
        tm = TaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test remainder(tmb) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))
    end

    @testset "Composition of functions and their inverses" begin
        tv = TaylorModel1(2, x0, ii0)

        tma = exp(tv)
        tmb = log(tma)
        @test tmb == log(exp(tv))
        @test tmb.pol == tv.pol

        tma = sin(tv)
        tmb = asin(tma)
        @test tmb == asin(sin(tv))
        @test tmb.pol == tv.pol

        tma = asin(tv)
        tmb = sin(tma)
        @test tmb == sin(asin(tv))
        @test tmb.pol == tv.pol

        tma = acos(tv)
        tmb = cos(tma)
        @test tmb == cos(acos(tv))
        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-16

        tma = tan(tv)
        tmb = atan(tma)
        @test tmb == atan(tan(tv))
        @test tmb.pol == tv.pol

        tma = atan(tv)
        tmb = tan(tma)
        @test tmb == tan(atan(tv))
        @test tmb.pol == tv.pol

        ####
        tv = TaylorModel1(2, x1, ii1)

        tma = log(tv)
        tmb = exp(tma)
        @test tmb == exp(log(tv))
        @test tmb.pol == tv.pol

        tma = cos(tv)
        tmb = acos(tma)
        @test tmb == acos(cos(tv))
        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-15
    end

    @testset "Tests for integrate" begin
        order = 4
        tm = TaylorModel1(order, x0, ii0)

        integ_res = integrate(exp(tm), 1..1)
        exact_res = exp(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem ⊆ integ_res.rem
        for ind = 1:_num_tests
            @test check_containment(exp, integ_res)
        end

        integ_res = integrate(cos(tm))
        exact_res = sin(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem ⊆ integ_res.rem
        for ind = 1:_num_tests
            @test check_containment(sin, integ_res)
        end

        integ_res = integrate(-sin(tm), 1..1)
        exact_res = cos(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem ⊆ integ_res.rem
        for ind = 1:_num_tests
            @test check_containment(cos, integ_res)
        end

        integ_res = integrate(1/(1+tm^2))
        exact_res = atan(tm)
        @test exact_res.pol == integ_res.pol
        # @test exact_res.rem ⊆ integ_res.rem
        for ind = 1:_num_tests
            @test check_containment(atan, integ_res)
        end
    end

    @testset "Display" begin
        tm = TaylorModel1(2, x1, ii1)
        use_show_default(true)
        @test string(exp(tm)) == "TaylorModel1{Interval{Float64},Float64}" *
            "(Taylor1{Interval{Float64}}(Interval{Float64}" *
            "[Interval(2.718281828459045, 2.7182818284590455), Interval(2.718281828459045, 2.7182818284590455), " *
            "Interval(1.3591409142295225, 1.3591409142295228)], 2), Interval(-0.05020487208677604, 0.06448109909211741), " *
            "Interval(1.0, 1.0), Interval(0.5, 1.5))"
        use_show_default(false)
        @test string(tm^3) == " Interval(1.0, 1.0) + Interval(3.0, 3.0) t + " *
            "Interval(3.0, 3.0) t² + Interval(-0.125, 0.125)"
        @test string(exp(tm)) == " Interval(2.718281828459045, 2.7182818284590455) + " *
            "Interval(2.718281828459045, 2.7182818284590455) t + " *
            "Interval(1.3591409142295225, 1.3591409142295228) t² + " *
            "Interval(-0.05020487208677604, 0.06448109909211741)"
    end

    @testset "Tests for bounders" begin
        @testset "Tests for linear dominated bounder" begin
            order = 3

            f(x) = 1 + x^5 - x^4
            D = 0.9375 .. 1
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_interval = f(D)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            @test diam(bound_ldb) <= diam(bound_interval)
            @test bound_ldb ⊆ bound_naive_tm

            D = 0.75 .. 0.8125
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_interval = f(D)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            @test diam(bound_ldb) <= diam(bound_interval)
            @test bound_ldb ⊆ bound_naive_tm

            f(x) = x^2 * sin(x)
            D = -1.875 .. -1.25
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_interval = f(D)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            @test diam(bound_ldb) <= diam(bound_interval)
            @test bound_ldb ⊆ bound_naive_tm

            D = 1.25 .. 1.875
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_interval = f(D)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            @test diam(bound_ldb) <= diam(bound_interval)
            @test bound_ldb ⊆ bound_naive_tm
        end

        @testset "Tests for quadratic fast bounder" begin
            order = 3

            f(x) = 1 + x^5 - x^4
            D = 0.75 .. 0.8125
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            bound_qfb = quadratic_fast_bounder(fT)
            @test bound_qfb ⊆ bound_naive_tm
            # @test diam(bound_qfb) <= diam(bound_ldb)

            f(x) = x^2 * sin(x)
            D = -2.5 .. -1.875
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            bound_qfb = quadratic_fast_bounder(fT)
            @test bound_qfb ⊆ bound_naive_tm
            @test diam(bound_qfb) <= diam(bound_ldb)

            f(x) = x^3 * cos(x) + x
            D = 3.75 .. 4.375
            x0 = mid(D)
            tm = TaylorModel1(order, x0, D)
            fT = f(tm)
            bound_naive_tm = fT(fT.dom - fT.x0)
            bound_ldb = linear_dominated_bounder(fT)
            bound_qfb = quadratic_fast_bounder(fT)
            @test bound_qfb ⊆ bound_naive_tm
            # @test diam(bound_qfb) <= diam(bound_ldb)
        end
    end
end

@testset "Tests for RTaylorModel1 " begin
    x0 = Interval(0.0)
    ii0 = Interval(-0.5, 0.5)
    x1 = Interval(1.0)
    ii1 = Interval(0.5, 1.5)

    @testset "RTaylorModel1 constructors" begin
        tv = RTaylorModel1{Interval{Float64},Float64}(Taylor1(Interval{Float64},5), x0, x0, ii0)
        @test tv == RTaylorModel1(Taylor1(Interval{Float64},5), x0, x0, ii0)
        @test tv == RTaylorModel1(5, x0, ii0)
        @test tv == RTaylorModel1(5, ii0)
        @test tv == RTaylorModel1(5, 0.0, ii0)
        @test RTaylorModel1(x1, 5, x0, ii0) == RTaylorModel1(Taylor1(x1, 5), x0, x0, ii0)
        @test RTaylorModel1(5, 0.7, ii1) == RTaylorModel1(5, interval(0.7), ii1)

        @test RTaylorModel1(tv, ii0) == RTaylorModel1(Taylor1(Interval{Float64}, 5), ii0, x0, ii0)
        @test RTaylorModel1(5, x0, ii0) == RTaylorModel1(tv, x0)
        @test RTaylorModel1(5, ii0) == RTaylorModel1(tv, x0)

        @test isa(tv, AbstractSeries)
        @test RTaylorModel1{Interval{Float64},Float64} <: AbstractSeries{Interval{Float64}}

        # Zero may not be contained in the remainder of a RTaylorModel1
        @test 0 ∉ remainder(RTaylorModel1(Taylor1(Interval{Float64},5), x1, x0, ii0))

        # Test errors in construction
        @test_throws AssertionError RTaylorModel1(5, x1, ii0)

        # Tests for get_order and remainder
        @test get_order(tv) == 5
        @test remainder(tv) == interval(0.0)
        @test polynomial(tv) == Taylor1(Interval{Float64},5)
        @test domain(tv) == ii0
        @test expansion_point(tv) == x0
        @test constant_term(tv) == interval(0.0)
        @test linear_polynomial(tv) == Taylor1(Interval{Float64},5)

        # Tests related to fixorder
        a = RTaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        b = RTaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
        aa, bb = TaylorModels.fixorder(a, b)
        @test get_order(aa) == get_order(bb) == 1
        @test isa(aa, RTaylorModel1) == isa(bb, RTaylorModel1) == true
        @test aa == a
        @test bb == RTaylorModel1(Taylor1([1.0, 1]), -1 .. 2, 0..0, -1 .. 1)
        # a and b remain the same
        @test a == RTaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        @test b == RTaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
    end

    @testset "Arithmetic operations" begin
        Δ = interval(-0.25, 0.25)
        a = RTaylorModel1(x1+Taylor1(5), Δ, x1, ii1)
        tv = RTaylorModel1(5, x1, ii1)

        @test zero(a) == RTaylorModel1(zero(a.pol), 0..0, x1, ii1)
        @test one(a) == RTaylorModel1(one(a.pol), 0..0, x1, ii1)
        @test a+x1 == RTaylorModel1(2*x1+Taylor1(5), Δ, x1, ii1)
        @test a+a == RTaylorModel1(2*(x1+Taylor1(5)), 2*Δ, x1, ii1)
        @test a-x1 == RTaylorModel1(zero(x1)+Taylor1(5), Δ, x1, ii1)
        @test a-a == RTaylorModel1(zero(a.pol), 2*Δ, x1, ii1)
        b = a * tv
        @test b == RTaylorModel1(a.pol*tv.pol, a.rem*tv.pol(ii1-x1), x1, ii1)
        @test remainder(b/tv) ⊆ Interval(-2.75, 4.75)
        b = a * a.pol[0]
        @test b == a
        @test constant_term(a) == x1
        @test linear_polynomial(a) == Taylor1(5)

        a = RTaylorModel1(x0, 5, x0, ii0)
        @test a^0 == RTaylorModel1(x0^0, 5, x0, ii0)
        @test a^1 == RTaylorModel1(x0^1, 5, x0, ii0)
        @test a^2 == RTaylorModel1(x0^2, 5, x0, ii0)
        @test a^3 == RTaylorModel1(x0^3, 5, x0, ii0)
        a = RTaylorModel1(x1, 5, x1, ii1)
        @test a^0 == RTaylorModel1(x1^0, 5, x1, ii1)
        @test a^1 == RTaylorModel1(x1^1, 5, x1, ii1)
        @test a^2 == RTaylorModel1(x1^2, 5, x1, ii1)
        @test a^3 == RTaylorModel1(x1^3, 5, x1, ii1)

        # Tests involving TM1s with different orders
        a = RTaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        b = RTaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)
        aa, bb = TaylorModels.fixorder(a, b)
        @test get_order(aa) == get_order(bb)
        @test get_order(bb) == 1
        @test a + b == aa + bb
        @test a - b == aa - bb
        res1 = a * b
        res2 = aa * bb
        @test res1 == RTaylorModel1(Taylor1([1.0, 2]), -1 .. 9 , 0..0, -1 .. 1)
        @test res2 == RTaylorModel1(Taylor1([1.0, 2]), -2 .. 9 , 0..0, -1 .. 1)
        res1 = a / b
        res2 = aa / bb
        @test res1 == RTaylorModel1(Taylor1([1.0, 0]), entireinterval() , 0..0, -1 .. 1)
        @test res2 == res1
        # a and b remain the same
        @test a == RTaylorModel1(Taylor1([1.0, 1]), 0..1, 0..0, -1 .. 1)
        @test b == RTaylorModel1(Taylor1([1.0, 1, 0, 1]), 0..1, 0..0, -1 .. 1)

        @test_throws AssertionError a+RTaylorModel1(a.pol, a.rem, 1..1, -1..1)
        @test_throws AssertionError a+RTaylorModel1(a.pol, a.rem, 0..0, -2..2)
    end

    @testset "RPAs, functions and remainders" begin
        @test rpa(x->5+zero(x), RTaylorModel1(4, x0, ii0)) ==
            RTaylorModel1(interval(5.0), 4, x0, ii0)
        @test rpa(x->5+one(x), RTaylorModel1(4, x1, ii1)) ==
            RTaylorModel1(5+x1, 4, x1, ii1)
        @test rpa(x->5*x, RTaylorModel1(4, x1, ii1)) ==
            RTaylorModel1(5.0*(x1+Taylor1(4)), x0, x1, ii1)
        @test rpa(x->5*x^0, RTaylorModel1(4, x0, ii0)) == 5*RTaylorModel1(4, x0, ii0)^0
        @test rpa(x->5*x^0, RTaylorModel1(4, x0, ii0)) ==
            RTaylorModel1( interval(5.0)*Taylor1(4)^0, x0, x0, ii0)
        @test rpa(x->5*x^1, RTaylorModel1(4, x0, ii0)) == 5*RTaylorModel1(4, x0, ii0)^1
        @test rpa(x->5*x^1, RTaylorModel1(4, x0, ii0)) ==
            RTaylorModel1( interval(5.0)*Taylor1(4)^1, x0, x0, ii0)
        @test rpa(x->5*x^2, RTaylorModel1(4, x0, ii0)) == 5*RTaylorModel1(4, x0, ii0)^2
        @test rpa(x->5*x^2, RTaylorModel1(4, x0, ii0)) ==
            RTaylorModel1( interval(5.0)*Taylor1(4)^2, x0, x0, ii0)
        @test rpa(x->5*x^4, RTaylorModel1(4, x0, ii0)) ==
            RTaylorModel1( interval(5.0)*Taylor1(4)^4, x0, x0, ii0)
        @test rpa(x->5*x^4, RTaylorModel1(3, x0, ii0)) ==
            RTaylorModel1( Taylor1(x0, 3), interval(5), x0, ii0)

        # Testing remainders and inclusion of RPAs
        order = 2
        ii = ii0
        xx = x0
        ftest = x -> exp(x)-1
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        # fT, Δ, ξ0, δ = fp_rpa(tma)
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)*(ii-ξ0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        # test for TM with scalar coefficients
        @test fp_rpa(tmc) == tmc

        order = 2
        ii = ii1
        xx = x1
        ftest = x -> exp(x)
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        # fT, Δ, ξ0, δ = fp_rpa(tma)
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)*(ii-ξ0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 3
        ii = ii0
        xx = x0
        ftest = x -> sin(x)
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        # fT, Δ, ξ0, δ = fp_rpa(tma)
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)*(ii-ξ0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 2
        ii = ii1
        xx = x1
        ftest = x -> sqrt(x)
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        # fT, Δ, ξ0, δ = fp_rpa(tma)
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.lo)-tmc.pol(ii.lo-ξ0),
                        ftest(ii.hi)-tmc.pol(ii.hi-ξ0)) ⊆ remainder(tma)*(ii-ξ0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        order = 5
        ii = ii1
        xx = x1
        ftest = x -> inv(x)
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test tma == tmb
        # fT, Δ, ξ0, δ = fp_rpa(tma)
        ξ0 = mid(xx, α_mid)
        tmc = fp_rpa(tma)
        @test interval(ftest(ii.hi)-tmc.pol(ii.hi-ξ0),
                        ftest(ii.lo)-tmc.pol(ii.lo-ξ0)) ⊆ remainder(tma)*(ii-ξ0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))

        # Example of Makino's thesis (page 98 and fig 4.2)
        order = 8
        ii = interval(-0.5, 1.0)
        xx = interval(mid(ii))
        ftest = x -> x*(x-1.1)*(x+2)*(x+2.2)*(x+2.5)*(x+3)*sin(1.7*x+0.5)
        tm = RTaylorModel1(order, xx, ii)
        tma = rpa(ftest, tm)
        tmb = ftest(tm)
        @test remainder(tmb) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, tma)
        end
        @test_throws AssertionError tmb(ii.hi+1.0)
        @test_throws AssertionError tmb(ii+Interval(1))
    end

    @testset "Composition of functions and their inverses" begin
        tv = RTaylorModel1(2, x0, ii0)

        tma = exp(tv)
        tmb = log(tma)
        @test tmb == log(exp(tv))
        @test tmb.pol == tv.pol

        tma = sin(tv)
        tmb = asin(tma)
        @test tmb == asin(sin(tv))
        @test tmb.pol == tv.pol

        tma = asin(tv)
        tmb = sin(tma)
        @test tmb == sin(asin(tv))
        @test tmb.pol == tv.pol

        tma = acos(tv)
        tmb = cos(tma)
        @test tmb == cos(acos(tv))
        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-16

        tma = tan(tv)
        tmb = atan(tma)
        @test tmb == atan(tan(tv))
        @test tmb.pol == tv.pol

        tma = atan(tv)
        tmb = tan(tma)
        @test tmb == tan(atan(tv))
        @test tmb.pol == tv.pol

        ####
        tv = RTaylorModel1(2, x1, ii1)

        tma = log(tv)
        tmb = exp(tma)
        @test tmb == exp(log(tv))
        @test tmb.pol == tv.pol

        tma = cos(tv)
        tmb = acos(tma)
        @test tmb == acos(cos(tv))
        @test sup(norm(tmb.pol - tv.pol, Inf)) < 1.0e-15
    end

    @testset "Tests for integrate" begin
        order = 4
        tm = RTaylorModel1(order, x0, ii0)

        integ_res = integrate(exp(tm), 1..1)
        exact_res = exp(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem*(ii0-x0)^(order+1) ⊆ integ_res.rem*(ii0-x0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(exp, integ_res)
        end

        integ_res = integrate(cos(tm))
        exact_res = sin(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem*(ii0-x0)^(order+1) ⊆ integ_res.rem*(ii0-x0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(sin, integ_res)
        end

        integ_res = integrate(-sin(tm), 1..1)
        exact_res = cos(tm)
        @test exact_res.pol == integ_res.pol
        @test exact_res.rem*(ii0-x0)^(order+1) ⊆ integ_res.rem*(ii0-x0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(cos, integ_res)
        end

        integ_res = integrate(1/(1+tm^2))
        exact_res = atan(tm)
        @test exact_res.pol == integ_res.pol
        # @test exact_res.rem*(ii0-x0)^(order+1) ⊆ integ_res.rem*(ii0-x0)^(order+1)
        for ind = 1:_num_tests
            @test check_containment(atan, integ_res)
        end
    end

    @testset "Display" begin
        tm = RTaylorModel1(3, x1, ii1)
        use_show_default(true)
        @test string(exp(tm)) == "RTaylorModel1{Interval{Float64},Float64}" *
            "(Taylor1{Interval{Float64}}(Interval{Float64}" *
            "[Interval(2.718281828459045, 2.7182818284590455), Interval(2.718281828459045, 2.7182818284590455), " *
            "Interval(1.3591409142295225, 1.3591409142295228), Interval(0.45304697140984085, 0.45304697140984096)], 3), " *
            "Interval(0.10281598943126369, 0.1256036426541982), Interval(1.0, 1.0), Interval(0.5, 1.5))"
        use_show_default(false)
        @test string(tm^3) == " Interval(1.0, 1.0) + Interval(3.0, 3.0) t + " *
            "Interval(3.0, 3.0) t² + Interval(1.0, 1.0) t³ + Interval(0.0, 0.0) t⁴"
        @test string(exp(tm)) == " Interval(2.718281828459045, 2.7182818284590455) + " *
            "Interval(2.718281828459045, 2.7182818284590455) t + Interval(1.3591409142295225, 1.3591409142295228) t² + " *
            "Interval(0.45304697140984085, 0.45304697140984096) t³ + Interval(0.10281598943126369, 0.1256036426541982) t⁴"
    end
end


text:
const TREE_DATA = let
  t = [Symbol("t$i") for i in 1:11]
  [
  (n = 3,
  id = "{1,2}",orbit_size=3,
  tree = [t[1] t[4] t[5];
         t[4] t[2] t[5];
         t[5] t[5] t[3]]),
            # 4
  #1,2,10
  (n = 4,
   id = "{1,2,3}",orbit_size=4,
   tree = [t[1] t[5] t[5] t[6];
     t[5] t[2] t[5] t[6];
     t[5] t[5] t[3] t[6];
     t[6] t[6] t[6] t[4]]),

  (n = 4,
   id = "{1,2}",orbit_size=6,
   tree = [t[1] t[5] t[6] t[6];
           t[5] t[2] t[6] t[6];
           t[6] t[6] t[3] t[6];
           t[6] t[6] t[6] t[4]]),
  (n = 4,
   id = "{1,2},{3,4}",orbit_size=3,
   tree = [t[1] t[5] t[7] t[7];
         t[5] t[2] t[7] t[7];
         t[7] t[7] t[3] t[6];
         t[7] t[7] t[6] t[4]]),

  (n = 4,
    id = "{1,2},{1,2,3}",orbit_size=12,
    tree = [t[1] t[5] t[6] t[7];
         t[5] t[2] t[6] t[7];
         t[6] t[6] t[3] t[7];
         t[7] t[7] t[7] t[4]]),
  #1,3,25
  (n=5,
   id="{1,2,3,4}",orbit_size=5,
 tree=[t[1] t[6] t[6] t[6] t[7];
   t[6] t[2] t[6] t[6] t[7];
   t[6] t[6] t[3] t[6] t[7];
   t[6] t[6] t[6] t[4] t[7];
   t[7] t[7] t[7] t[7] t[5]]),

  (n=5,
   id="{1,2}",orbit_size=10,
 tree=[t[1] t[6] t[7] t[7] t[7];
   t[6] t[2] t[7] t[7] t[7];
   t[7] t[7] t[3] t[7] t[7];
   t[7] t[7] t[7] t[4] t[7];
   t[7] t[7] t[7] t[7] t[5]]),

  (n=5,
   id="{1,2,3}",orbit_size=10,
 tree=[t[1] t[6] t[6] t[7] t[7];
   t[6] t[2] t[6] t[7] t[7];
   t[6] t[6] t[3] t[7] t[7];
   t[7] t[7] t[7] t[4] t[7];
   t[7] t[7] t[7] t[7] t[5]]),

  # 2,5,105
  (n=5,
   id="{1,2},{3,4,5}",orbit_size=10,
 tree=[t[1] t[6] t[8] t[8] t[8];
   t[6] t[2] t[8] t[8] t[8];
   t[8] t[8] t[3] t[7] t[7];
   t[8] t[8] t[7] t[4] t[7];
   t[8] t[8] t[7] t[7] t[5]]),

  (n=5,
   id="{1,2},{3,4}",orbit_size=15,
 tree=[t[1] t[6] t[8] t[8] t[8];
   t[6] t[2] t[8] t[8] t[8];
   t[8] t[8] t[3] t[7] t[8];
   t[8] t[8] t[7] t[4] t[8];
   t[8] t[8] t[8] t[8] t[5]]),

  (n=5,
   id="{1,2,3},{1,2,3,4}",orbit_size=20,
 tree=[t[1] t[6] t[6] t[7] t[8];
   t[6] t[2] t[6] t[7] t[8];
   t[6] t[6] t[3] t[7] t[8];
   t[7] t[7] t[7] t[4] t[8];
   t[8] t[8] t[8] t[8] t[5]]),

  (n=5,
   id="{1,2},{1,2,3}",orbit_size=30,
 tree=[t[1] t[6] t[7] t[8] t[8];
   t[6] t[2] t[7] t[8] t[8];
   t[7] t[7] t[3] t[8] t[8];
   t[8] t[8] t[8] t[4] t[8];
   t[8] t[8] t[8] t[8] t[5]]),

  (n=5,
   id="{1,2},{1,2,3,4}",orbit_size=30,
 tree=[t[1] t[6] t[7] t[7] t[8];
   t[6] t[2] t[7] t[7] t[8];
   t[7] t[7] t[3] t[7] t[8];
   t[7] t[7] t[7] t[4] t[8];
   t[8] t[8] t[8] t[8] t[5]]),

  # 3,3,105
  (n=5,
   id="{1,2},{3,4},{1,2,3,4}",orbit_size=15,
 tree=[t[1] t[6] t[8] t[8] t[9];
   t[6] t[2] t[8] t[8] t[9];
   t[8] t[8] t[3] t[7] t[9];
   t[8] t[8] t[7] t[4] t[9];
   t[9] t[9] t[9] t[9] t[5]]),

  (n=5,
   id="{1,2},{3,4},{1,2,5}",orbit_size=30,
 tree=[t[1] t[6] t[9] t[9] t[8];
   t[6] t[2] t[9] t[9] t[8];
   t[9] t[9] t[3] t[7] t[9];
   t[9] t[9] t[7] t[4] t[9];
   t[8] t[8] t[9] t[9] t[5]]),

  (n=5,
   id="{1,2},{1,2,3},{1,2,3,4}",orbit_size=60,
 tree=[t[1] t[6] t[7] t[8] t[9];
   t[6] t[2] t[7] t[8] t[9];
   t[7] t[7] t[3] t[8] t[9];
   t[8] t[8] t[8] t[4] t[9];
   t[9] t[9] t[9] t[9] t[5]]),
  #         1.726
  #
  #           6
  #
  # 1,4,56
  (n=6,
   id="{1,2,3,4,5}",orbit_size=6,
 tree=[t[1] t[7] t[7] t[7] t[7] t[8];
   t[7] t[2] t[7] t[7] t[7] t[8];
   t[7] t[7] t[3] t[7] t[7] t[8];
   t[7] t[7] t[7] t[4] t[7] t[8];
   t[7] t[7] t[7] t[7] t[5] t[8];
   t[8] t[8] t[8] t[8] t[8] t[6]]),

  (n=6,
   id="{1,2}",orbit_size=15,
 tree=[t[1] t[7] t[8] t[8] t[8] t[8];
   t[7] t[2] t[8] t[8] t[8] t[8];
   t[8] t[8] t[3] t[8] t[8] t[8];
   t[8] t[8] t[8] t[4] t[8] t[8];
   t[8] t[8] t[8] t[8] t[5] t[8];
   t[8] t[8] t[8] t[8] t[8] t[6]]),

  (n=6,
   id="{1,2,3,4}",orbit_size=15,
 tree=[t[1] t[7] t[7] t[7] t[8] t[8];
   t[7] t[2] t[7] t[7] t[8] t[8];
   t[7] t[7] t[3] t[7] t[8] t[8];
   t[7] t[7] t[7] t[4] t[8] t[8];
   t[8] t[8] t[8] t[8] t[5] t[8];
   t[8] t[8] t[8] t[8] t[8] t[6]]),

  (n=6,
   id="{1,2,3}",orbit_size=20,
 tree=[t[1] t[7] t[7] t[8] t[8] t[8];
   t[7] t[2] t[7] t[8] t[8] t[8];
   t[7] t[7] t[3] t[8] t[8] t[8];
   t[8] t[8] t[8] t[4] t[8] t[8];
   t[8] t[8] t[8] t[8] t[5] t[8];
   t[8] t[8] t[8] t[8] t[8] t[6]]),

  # 2,10,490
  (n=6,
   id="{1,2,3},{4,5,6}",orbit_size=10,
 tree=[t[1] t[7] t[7] t[9] t[9] t[9];
   t[7] t[2] t[7] t[9] t[9] t[9];
   t[7] t[7] t[3] t[9] t[9] t[9];
   t[9] t[9] t[9] t[4] t[8] t[8];
   t[9] t[9] t[9] t[8] t[5] t[8];
   t[9] t[9] t[9] t[8] t[8] t[6]]),

  (n=6,
   id="{1,2},{3,4,5,6}",orbit_size=15,
 tree=[t[1] t[7] t[9] t[9] t[9] t[9];
   t[7] t[2] t[9] t[9] t[9] t[9];
   t[9] t[9] t[3] t[8] t[8] t[8];
   t[9] t[9] t[8] t[4] t[8] t[8];
   t[9] t[9] t[8] t[8] t[5] t[8];
   t[9] t[9] t[8] t[8] t[8] t[6]]),

  (n=6,
   id="{1,2,3,4},{1,2,3,4,5}",orbit_size=30,
 tree=[t[1] t[7] t[7] t[7] t[8] t[9];
   t[7] t[2] t[7] t[7] t[8] t[9];
   t[7] t[7] t[3] t[7] t[8] t[9];
   t[7] t[7] t[7] t[4] t[8] t[9];
   t[8] t[8] t[8] t[8] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4}",orbit_size=45,
 tree=[t[1] t[7] t[9] t[9] t[9] t[9];
   t[7] t[2] t[9] t[9] t[9] t[9];
   t[9] t[9] t[3] t[8] t[9] t[9];
   t[9] t[9] t[8] t[4] t[9] t[9];
   t[9] t[9] t[9] t[9] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{1,2,3}",orbit_size=60,
 tree=[t[1] t[7] t[8] t[9] t[9] t[9];
   t[7] t[2] t[8] t[9] t[9] t[9];
   t[8] t[8] t[3] t[9] t[9] t[9];
   t[9] t[9] t[9] t[4] t[9] t[9];
   t[9] t[9] t[9] t[9] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4,5}",orbit_size=60,
 tree=[t[1] t[7] t[9] t[9] t[9] t[9];
   t[7] t[2] t[9] t[9] t[9] t[9];
   t[9] t[9] t[3] t[8] t[8] t[9];
   t[9] t[9] t[8] t[4] t[8] t[9];
   t[9] t[9] t[8] t[8] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{1,2,3,4,5}",orbit_size=60,
 tree=[t[1] t[7] t[8] t[8] t[8] t[9];
   t[7] t[2] t[8] t[8] t[8] t[9];
   t[8] t[8] t[3] t[8] t[8] t[9];
   t[8] t[8] t[8] t[4] t[8] t[9];
   t[8] t[8] t[8] t[8] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2,3},{1,2,3,4}",orbit_size=60,
 tree=[t[1] t[7] t[7] t[8] t[9] t[9];
   t[7] t[2] t[7] t[8] t[9] t[9];
   t[7] t[7] t[3] t[8] t[9] t[9];
   t[8] t[8] t[8] t[4] t[9] t[9];
   t[9] t[9] t[9] t[9] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2,3},{1,2,3,4,5}",orbit_size=60,
 tree=[t[1] t[7] t[7] t[8] t[8] t[9];
   t[7] t[2] t[7] t[8] t[8] t[9];
   t[7] t[7] t[3] t[8] t[8] t[9];
   t[8] t[8] t[8] t[4] t[8] t[9];
   t[8] t[8] t[8] t[8] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{1,2,3,4}",orbit_size=90,
 tree=[t[1] t[7] t[8] t[8] t[9] t[9];
   t[7] t[2] t[8] t[8] t[9] t[9];
   t[8] t[8] t[3] t[8] t[9] t[9];
   t[8] t[8] t[8] t[4] t[9] t[9];
   t[9] t[9] t[9] t[9] t[5] t[9];
   t[9] t[9] t[9] t[9] t[9] t[6]]),

  # 3,12,1260
  (n=6,
   id="{1,2},{3,4},{5,6}",orbit_size=15,
 tree=[t[1] t[7] t[10] t[10] t[10] t[10];
   t[7] t[2] t[10] t[10] t[10] t[10];
   t[10] t[10] t[3] t[8] t[10] t[10];
   t[10] t[10] t[8] t[4] t[10] t[10];
   t[10] t[10] t[10] t[10] t[5] t[9] ;
   t[10] t[10] t[10] t[10] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,3,4}",orbit_size=45,
 tree=[t[1] t[7] t[9] t[9] t[10] t[10];
   t[7] t[2] t[9] t[9] t[10] t[10];
   t[9] t[9] t[3] t[8] t[10] t[10];
   t[9] t[9] t[8] t[4] t[10] t[10];
   t[10] t[10] t[10] t[10] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{1,2,3},{4,5,6}",orbit_size=60,
 tree=[t[1] t[7] t[8] t[10] t[10] t[10];
   t[7] t[2] t[8] t[10] t[10] t[10];
   t[8] t[8] t[3] t[10] t[10] t[10];
   t[10] t[10] t[10] t[4] t[9] t[9] ;
   t[10] t[10] t[10] t[9] t[5] t[9] ;
   t[10] t[10] t[10] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4,5},{3,4,5,6}",orbit_size=60,
 tree=[t[1] t[7] t[10] t[10] t[10] t[10];
   t[7] t[2] t[10] t[10] t[10] t[10];
   t[10] t[10] t[3] t[8] t[8] t[9] ;
   t[10] t[10] t[8] t[4] t[8] t[9] ;
   t[10] t[10] t[8] t[8] t[5] t[9] ;
   t[10] t[10] t[9] t[9] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4,5},{1,2,3,4,5}",orbit_size=60,
 tree=[t[1] t[7] t[9] t[9] t[9] t[10];
   t[7] t[2] t[9] t[9] t[9] t[10];
   t[9] t[9] t[3] t[8] t[8] t[10];
   t[9] t[9] t[8] t[4] t[8] t[10];
   t[9] t[9] t[8] t[8] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,5,6}",orbit_size=90,
 tree=[t[1] t[7] t[10] t[10] t[9] t[9] ;
   t[7] t[2] t[10] t[10] t[9] t[9] ;
   t[10] t[10] t[3] t[8] t[10] t[10];
   t[10] t[10] t[8] t[4] t[10] t[10];
   t[9] t[9] t[10] t[10] t[5] t[9] ;
   t[9] t[9] t[10] t[10] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,3,4,5}",orbit_size=90,
 tree=[t[1] t[7] t[9] t[9] t[9] t[10];
   t[7] t[2] t[9] t[9] t[9] t[10];
   t[9] t[9] t[3] t[8] t[9] t[10];
   t[9] t[9] t[8] t[4] t[9] t[10];
   t[9] t[9] t[9] t[9] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2,3},{1,2,3,4},{1,2,3,4,5}",orbit_size=120,
 tree=[t[1] t[7] t[7] t[8] t[9] t[10];
   t[7] t[2] t[7] t[8] t[9] t[10];
   t[7] t[7] t[3] t[8] t[9] t[10];
   t[8] t[8] t[8] t[4] t[9] t[10];
   t[9] t[9] t[9] t[9] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,5}",orbit_size=180,
 tree=[t[1] t[7] t[10] t[10] t[9] t[10];
   t[7] t[2] t[10] t[10] t[9] t[10];
   t[10] t[10] t[3] t[8] t[10] t[10];
   t[10] t[10] t[8] t[4] t[10] t[10];
   t[9] t[9] t[10] t[10] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{1,2,3},{1,2,3,4}",orbit_size=180,
 tree=[t[1] t[7] t[8] t[9] t[10] t[10];
   t[7] t[2] t[8] t[9] t[10] t[10];
   t[8] t[8] t[3] t[9] t[10] t[10];
   t[9] t[9] t[9] t[4] t[10] t[10];
   t[10] t[10] t[10] t[10] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{1,2,3},{1,2,3,4,5}",orbit_size=180,
 tree=[t[1] t[7] t[8] t[9] t[9] t[10];
   t[7] t[2] t[8] t[9] t[9] t[10];
   t[8] t[8] t[3] t[9] t[9] t[10];
   t[9] t[9] t[9] t[4] t[9] t[10];
   t[9] t[9] t[9] t[9] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  (n=6,
   id="{1,2},{1,2,3,4},{1,2,3,4,5}",orbit_size=180,
 tree=[t[1] t[7] t[8] t[8] t[9] t[10];
   t[7] t[2] t[8] t[8] t[9] t[10];
   t[8] t[8] t[3] t[8] t[9] t[10];
   t[8] t[8] t[8] t[4] t[9] t[10];
   t[9] t[9] t[9] t[9] t[5] t[10];
   t[10] t[10] t[10] t[10] t[10] t[6]]),

  #4,6,945
  (n=6,
   id="{1,2},{3,4},{5,6},{1,2,3,4}",orbit_size=45,
 tree=[t[1] t[7] t[10] t[10] t[11] t[11];
   t[7] t[2] t[10] t[10] t[11] t[11];
   t[10] t[10] t[3] t[8] t[11] t[11];
   t[10] t[10] t[8] t[4] t[11] t[11];
   t[11] t[11] t[11] t[11] t[5] t[9] ;
   t[11] t[11] t[11] t[11] t[9] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,5},{3,4,6}",orbit_size=90,
 tree=[t[1] t[7] t[11] t[11] t[9] t[11];
   t[7] t[2] t[11] t[11] t[9] t[11];
   t[11] t[11] t[3] t[8] t[11] t[10];
   t[11] t[11] t[8] t[4] t[11] t[10];
   t[9] t[9] t[11] t[11] t[5] t[11];
   t[11] t[11] t[10] t[10] t[11] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,3,4},{1,2,3,4,5}",orbit_size=90,
 tree=[t[1] t[7] t[9] t[9] t[10] t[11];
   t[7] t[2] t[9] t[9] t[10] t[11];
   t[9] t[9] t[3] t[8] t[10] t[11];
   t[9] t[9] t[8] t[4] t[10] t[11];
   t[10] t[10] t[10] t[10] t[5] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,5},{1,2,5,6}",orbit_size=180,
 tree=[t[1] t[7] t[11] t[11] t[9] t[10];
   t[7] t[2] t[11] t[11] t[9] t[10];
   t[11] t[11] t[3] t[8] t[11] t[11];
   t[11] t[11] t[8] t[4] t[11] t[11];
   t[9] t[9] t[11] t[11] t[5] t[10];
   t[10] t[10] t[11] t[11] t[10] t[6]]),

  (n=6,
   id="{1,2},{3,4},{1,2,5},{1,2,3,4,5}",orbit_size=180,
 tree=[t[1] t[7] t[10] t[10] t[9] t[11];
   t[7] t[2] t[10] t[10] t[9] t[11];
   t[10] t[10] t[3] t[8] t[10] t[11];
   t[10] t[10] t[8] t[4] t[10] t[11];
   t[9] t[9] t[10] t[10] t[5] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6]]),

  (n=6,
   id="{1,2},{1,2,3},{1,2,3,4},{1,2,3,4,5}",orbit_size=360,
 tree=[t[1] t[7] t[8] t[9] t[10] t[11];
   t[7] t[2] t[8] t[9] t[10] t[11];
   t[8] t[8] t[3] t[9] t[10] t[11];
   t[9] t[9] t[9] t[4] t[10] t[11];
   t[10] t[10] t[10] t[10] t[5] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6]]),

    #        215.768
    #
    #           7
    #
    # 1,5,119
    (n=7,
     id="{1,2,3,4,5,6}",orbit_size=7,
   tree=[t[1] t[8] t[8] t[8] t[8] t[8] t[9];
     t[8] t[2] t[8] t[8] t[8] t[8] t[9];
     t[8] t[8] t[3] t[8] t[8] t[8] t[9];
     t[8] t[8] t[8] t[4] t[8] t[8] t[9];
     t[8] t[8] t[8] t[8] t[5] t[8] t[9];
     t[8] t[8] t[8] t[8] t[8] t[6] t[9];
     t[9] t[9] t[9] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2}",orbit_size=21,
   tree=[t[1] t[8] t[9] t[9] t[9] t[9] t[9];
     t[8] t[2] t[9] t[9] t[9] t[9] t[9];
     t[9] t[9] t[3] t[9] t[9] t[9] t[9];
     t[9] t[9] t[9] t[4] t[9] t[9] t[9];
     t[9] t[9] t[9] t[9] t[5] t[9] t[9];
     t[9] t[9] t[9] t[9] t[9] t[6] t[9];
     t[9] t[9] t[9] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2,3,4,5}",orbit_size=21,
   tree=[t[1] t[8] t[8] t[8] t[8] t[9] t[9];
     t[8] t[2] t[8] t[8] t[8] t[9] t[9];
     t[8] t[8] t[3] t[8] t[8] t[9] t[9];
     t[8] t[8] t[8] t[4] t[8] t[9] t[9];
     t[8] t[8] t[8] t[8] t[5] t[9] t[9];
     t[9] t[9] t[9] t[9] t[9] t[6] t[9];
     t[9] t[9] t[9] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2,3}",orbit_size=35,
   tree=[t[1] t[8] t[8] t[9] t[9] t[9] t[9];
     t[8] t[2] t[8] t[9] t[9] t[9] t[9];
     t[8] t[8] t[3] t[9] t[9] t[9] t[9];
     t[9] t[9] t[9] t[4] t[9] t[9] t[9];
     t[9] t[9] t[9] t[9] t[5] t[9] t[9];
     t[9] t[9] t[9] t[9] t[9] t[6] t[9];
     t[9] t[9] t[9] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2,3,4}",orbit_size=35,
   tree=[t[1] t[8] t[8] t[8] t[9] t[9] t[9];
     t[8] t[2] t[8] t[8] t[9] t[9] t[9];
     t[8] t[8] t[3] t[8] t[9] t[9] t[9];
     t[8] t[8] t[8] t[4] t[9] t[9] t[9];
     t[9] t[9] t[9] t[9] t[5] t[9] t[9];
     t[9] t[9] t[9] t[9] t[9] t[6] t[9];
     t[9] t[9] t[9] t[9] t[9] t[9] t[7]]),

    # 2,16,1918
    (n=7,
     id="{1,2},{3,4,5,6,7}",orbit_size=21,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[10];
     t[8] t[2] t[10] t[10] t[10] t[10] t[10];
     t[10] t[10] t[3] t[9] t[9] t[9] t[9] ;
     t[10] t[10] t[9] t[4] t[9] t[9] t[9] ;
     t[10] t[10] t[9] t[9] t[5] t[9] t[9] ;
     t[10] t[10] t[9] t[9] t[9] t[6] t[9] ;
     t[10] t[10] t[9] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2,3},{4,5,6,7}",orbit_size=35,
   tree=[t[1] t[8] t[8] t[10] t[10] t[10] t[10];
     t[8] t[2] t[8] t[10] t[10] t[10] t[10];
     t[8] t[8] t[3] t[10] t[10] t[10] t[10];
     t[10] t[10] t[10] t[4] t[9] t[9] t[9] ;
     t[10] t[10] t[10] t[9] t[5] t[9] t[9] ;
     t[10] t[10] t[10] t[9] t[9] t[6] t[9] ;
     t[10] t[10] t[10] t[9] t[9] t[9] t[7]]),

    (n=7,
     id="{1,2,3,4,5},{1,2,3,4,5,6}",orbit_size=42,
   tree=[t[1] t[8] t[8] t[8] t[8] t[9] t[10];
     t[8] t[2] t[8] t[8] t[8] t[9] t[10];
     t[8] t[8] t[3] t[8] t[8] t[9] t[10];
     t[8] t[8] t[8] t[4] t[8] t[9] t[10];
     t[8] t[8] t[8] t[8] t[5] t[9] t[10];
     t[9] t[9] t[9] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3},{4,5,6}",orbit_size=70,
   tree=[t[1] t[8] t[8] t[10] t[10] t[10] t[10];
     t[8] t[2] t[8] t[10] t[10] t[10] t[10];
     t[8] t[8] t[3] t[10] t[10] t[10] t[10];
     t[10] t[10] t[10] t[4] t[9] t[9] t[10];
     t[10] t[10] t[10] t[9] t[5] t[9] t[10];
     t[10] t[10] t[10] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4}",orbit_size=105,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[10];
     t[8] t[2] t[10] t[10] t[10] t[10] t[10];
     t[10] t[10] t[3] t[9] t[10] t[10] t[10];
     t[10] t[10] t[9] t[4] t[10] t[10] t[10];
     t[10] t[10] t[10] t[10] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{1,2,3}",orbit_size=105,
   tree=[t[1] t[8] t[9] t[10] t[10] t[10] t[10];
     t[8] t[2] t[9] t[10] t[10] t[10] t[10];
     t[9] t[9] t[3] t[10] t[10] t[10] t[10];
     t[10] t[10] t[10] t[4] t[10] t[10] t[10];
     t[10] t[10] t[10] t[10] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5,6}",orbit_size=105,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[10];
     t[8] t[2] t[10] t[10] t[10] t[10] t[10];
     t[10] t[10] t[3] t[9] t[9] t[9] t[10];
     t[10] t[10] t[9] t[4] t[9] t[9] t[10];
     t[10] t[10] t[9] t[9] t[5] t[9] t[10];
     t[10] t[10] t[9] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{1,2,3,4,5,6}",orbit_size=105,
   tree=[t[1] t[8] t[9] t[9] t[9] t[9] t[10];
     t[8] t[2] t[9] t[9] t[9] t[9] t[10];
     t[9] t[9] t[3] t[9] t[9] t[9] t[10];
     t[9] t[9] t[9] t[4] t[9] t[9] t[10];
     t[9] t[9] t[9] t[9] t[5] t[9] t[10];
     t[9] t[9] t[9] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3,4},{1,2,3,4,5}",orbit_size=105,
   tree=[t[1] t[8] t[8] t[8] t[9] t[10] t[10];
     t[8] t[2] t[8] t[8] t[9] t[10] t[10];
     t[8] t[8] t[3] t[8] t[9] t[10] t[10];
     t[8] t[8] t[8] t[4] t[9] t[10] t[10];
     t[9] t[9] t[9] t[9] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3,4},{1,2,3,4,5,6}",orbit_size=105,
   tree=[t[1] t[8] t[8] t[8] t[9] t[9] t[10];
     t[8] t[2] t[8] t[8] t[9] t[9] t[10];
     t[8] t[8] t[3] t[8] t[9] t[9] t[10];
     t[8] t[8] t[8] t[4] t[9] t[9] t[10];
     t[9] t[9] t[9] t[9] t[5] t[9] t[10];
     t[9] t[9] t[9] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3},{1,2,3,4}",orbit_size=140,
   tree=[t[1] t[8] t[8] t[9] t[10] t[10] t[10];
     t[8] t[2] t[8] t[9] t[10] t[10] t[10];
     t[8] t[8] t[3] t[9] t[10] t[10] t[10];
     t[9] t[9] t[9] t[4] t[10] t[10] t[10];
     t[10] t[10] t[10] t[10] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3},{1,2,3,4,5,6}",orbit_size=140,
   tree=[t[1] t[8] t[8] t[9] t[9] t[9] t[10];
     t[8] t[2] t[8] t[9] t[9] t[9] t[10];
     t[8] t[8] t[3] t[9] t[9] t[9] t[10];
     t[9] t[9] t[9] t[4] t[9] t[9] t[10];
     t[9] t[9] t[9] t[9] t[5] t[9] t[10];
     t[9] t[9] t[9] t[9] t[9] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5}",orbit_size=210,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[10];
     t[8] t[2] t[10] t[10] t[10] t[10] t[10];
     t[10] t[10] t[3] t[9] t[9] t[10] t[10];
     t[10] t[10] t[9] t[4] t[9] t[10] t[10];
     t[10] t[10] t[9] t[9] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{1,2,3,4}",orbit_size=210,
   tree=[t[1] t[8] t[9] t[9] t[10] t[10] t[10];
     t[8] t[2] t[9] t[9] t[10] t[10] t[10];
     t[9] t[9] t[3] t[9] t[10] t[10] t[10];
     t[9] t[9] t[9] t[4] t[10] t[10] t[10];
     t[10] t[10] t[10] t[10] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{1,2,3,4,5}",orbit_size=210,
   tree=[t[1] t[8] t[9] t[9] t[9] t[10] t[10];
     t[8] t[2] t[9] t[9] t[9] t[10] t[10];
     t[9] t[9] t[3] t[9] t[9] t[10] t[10];
     t[9] t[9] t[9] t[4] t[9] t[10] t[10];
     t[9] t[9] t[9] t[9] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3},{1,2,3,4,5}",orbit_size=210,
   tree=[t[1] t[8] t[8] t[9] t[9] t[10] t[10];
     t[8] t[2] t[8] t[9] t[9] t[10] t[10];
     t[8] t[8] t[3] t[9] t[9] t[10] t[10];
     t[9] t[9] t[9] t[4] t[9] t[10] t[10];
     t[9] t[9] t[9] t[9] t[5] t[10] t[10];
     t[10] t[10] t[10] t[10] t[10] t[6] t[10];
     t[10] t[10] t[10] t[10] t[10] t[10] t[7]]),

    # 3,29,9450
    (n=7,
     id="{1,2,3},{4,5,6},{1,2,3,4,5,6}",orbit_size=70,
   tree=[t[1] t[8] t[8] t[10] t[10] t[10] t[11];
     t[8] t[2] t[8] t[10] t[10] t[10] t[11];
     t[8] t[8] t[3] t[10] t[10] t[10] t[11];
     t[10] t[10] t[10] t[4] t[9] t[9] t[11];
     t[10] t[10] t[10] t[9] t[5] t[9] t[11];
     t[10] t[10] t[10] t[9] t[9] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4},{5,6}",orbit_size=105,
   tree=[t[1] t[8] t[11] t[11] t[11] t[11] t[11];
     t[8] t[2] t[11] t[11] t[11] t[11] t[11];
     t[11] t[11] t[3] t[9] t[11] t[11] t[11];
     t[11] t[11] t[9] t[4] t[11] t[11] t[11];
     t[11] t[11] t[11] t[11] t[5] t[10] t[11];
     t[11] t[11] t[11] t[11] t[10] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4},{5,6,7}",orbit_size=105,
   tree=[t[1] t[8] t[11] t[11] t[11] t[11] t[11];
     t[8] t[2] t[11] t[11] t[11] t[11] t[11];
     t[11] t[11] t[3] t[9] t[11] t[11] t[11];
     t[11] t[11] t[9] t[4] t[11] t[11] t[11];
     t[11] t[11] t[11] t[11] t[5] t[10] t[10];
     t[11] t[11] t[11] t[11] t[10] t[6] t[10];
     t[11] t[11] t[11] t[11] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4},{1,2,3,4}",orbit_size=105,
   tree=[t[1] t[8] t[10] t[10] t[11] t[11] t[11];
     t[8] t[2] t[10] t[10] t[11] t[11] t[11];
     t[10] t[10] t[3] t[9] t[11] t[11] t[11];
     t[10] t[10] t[9] t[4] t[11] t[11] t[11];
     t[11] t[11] t[11] t[11] t[5] t[11] t[11];
     t[11] t[11] t[11] t[11] t[11] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{1,2,3},{4,5,6,7}",orbit_size=105,
   tree=[t[1] t[8] t[9] t[11] t[11] t[11] t[11];
     t[8] t[2] t[9] t[11] t[11] t[11] t[11];
     t[9] t[9] t[3] t[11] t[11] t[11] t[11];
     t[11] t[11] t[11] t[4] t[10] t[10] t[10];
     t[11] t[11] t[11] t[10] t[5] t[10] t[10];
     t[11] t[11] t[11] t[10] t[10] t[6] t[10];
     t[11] t[11] t[11] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5,6},{3,4,5,6,7}",orbit_size=105,
   tree=[t[1] t[8] t[11] t[11] t[11] t[11] t[11];
     t[8] t[2] t[11] t[11] t[11] t[11] t[11];
     t[11] t[11] t[3] t[9] t[9] t[9] t[10];
     t[11] t[11] t[9] t[4] t[9] t[9] t[10];
     t[11] t[11] t[9] t[9] t[5] t[9] t[10];
     t[11] t[11] t[9] t[9] t[9] t[6] t[10];
     t[11] t[11] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5,6},{1,2,3,4,5,6}",orbit_size=105,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[11];
     t[8] t[2] t[10] t[10] t[10] t[10] t[11];
     t[10] t[10] t[3] t[9] t[9] t[9] t[11];
     t[10] t[10] t[9] t[4] t[9] t[9] t[11];
     t[10] t[10] t[9] t[9] t[5] t[9] t[11];
     t[10] t[10] t[9] t[9] t[9] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2,3},{4,5,6},{1,2,3,7}",orbit_size=140,
   tree=[t[1] t[8] t[8] t[11] t[11] t[11] t[10];
     t[8] t[2] t[8] t[11] t[11] t[11] t[10];
     t[8] t[8] t[3] t[11] t[11] t[11] t[10];
     t[11] t[11] t[11] t[4] t[9] t[9] t[11];
     t[11] t[11] t[11] t[9] t[5] t[9] t[11];
     t[11] t[11] t[11] t[9] t[9] t[6] t[11];
     t[10] t[10] t[10] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4},{1,2,5,6,7}",orbit_size=210,
   tree=[t[1] t[8] t[11] t[11] t[10] t[10] t[10];
     t[8] t[2] t[11] t[11] t[10] t[10] t[10];
     t[11] t[11] t[3] t[9] t[11] t[11] t[11];
     t[11] t[11] t[9] t[4] t[11] t[11] t[11];
     t[10] t[10] t[11] t[11] t[5] t[10] t[10];
     t[10] t[10] t[11] t[11] t[10] t[6] t[10];
     t[10] t[10] t[11] t[11] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5},{1,2,6,7}",orbit_size=210,
   tree=[t[1] t[8] t[11] t[11] t[11] t[10] t[10];
     t[8] t[2] t[11] t[11] t[11] t[10] t[10];
     t[11] t[11] t[3] t[9] t[9] t[11] t[11];
     t[11] t[11] t[9] t[4] t[9] t[11] t[11];
     t[11] t[11] t[9] t[9] t[5] t[11] t[11];
     t[10] t[10] t[11] t[11] t[11] t[6] t[10];
     t[10] t[10] t[11] t[11] t[11] t[10] t[7]]),

    (n=7,
     id="{1,2},{3,4,5},{1,2,3,4,5}",orbit_size=210,
   tree=[t[1] t[8] t[10] t[10] t[10] t[11] t[11];
     t[8] t[2] t[10] t[10] t[10] t[11] t[11];
     t[10] t[10] t[3] t[9] t[9] t[11] t[11];
     t[10] t[10] t[9] t[4] t[9] t[11] t[11];
     t[10] t[10] t[9] t[9] t[5] t[11] t[11];
     t[11] t[11] t[11] t[11] t[11] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4,5},{3,4,5,6,7}",orbit_size=210,
   tree=[t[1] t[8] t[11] t[11] t[11] t[11] t[11];
     t[8] t[2] t[11] t[11] t[11] t[11] t[11];
     t[11] t[11] t[3] t[9] t[9] t[10] t[10];
     t[11] t[11] t[9] t[4] t[9] t[10] t[10];
     t[11] t[11] t[9] t[9] t[5] t[10] t[10];
     t[11] t[11] t[10] t[10] t[10] t[6] t[10];
     t[11] t[11] t[10] t[10] t[10] t[10] t[7]]),

    (n=7,
     id="{1,2,3,4},{1,2,3,4,5},{1,2,3,4,5,6}",orbit_size=210,
   tree=[t[1] t[8] t[8] t[8] t[9] t[10] t[11];
     t[8] t[2] t[8] t[8] t[9] t[10] t[11];
     t[8] t[8] t[3] t[8] t[9] t[10] t[11];
     t[8] t[8] t[8] t[4] t[9] t[10] t[11];
     t[9] t[9] t[9] t[9] t[5] t[10] t[11];
     t[10] t[10] t[10] t[10] t[10] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4},{1,2,3,4,5}",orbit_size=315,
   tree=[t[1] t[8] t[10] t[10] t[10] t[11] t[11];
     t[8] t[2] t[10] t[10] t[10] t[11] t[11];
     t[10] t[10] t[3] t[9] t[10] t[11] t[11];
     t[10] t[10] t[9] t[4] t[10] t[11] t[11];
     t[10] t[10] t[10] t[10] t[5] t[11] t[11];
     t[11] t[11] t[11] t[11] t[11] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

    (n=7,
     id="{1,2},{3,4},{1,2,3,4,5,6}",orbit_size=315,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[11];
     t[8] t[2] t[10] t[10] t[10] t[10] t[11];
     t[10] t[10] t[3] t[9] t[10] t[10] t[11];
     t[10] t[10] t[9] t[4] t[10] t[10] t[11];
     t[10] t[10] t[10] t[10] t[5] t[10] t[11];
     t[10] t[10] t[10] t[10] t[10] t[6] t[11];
     t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3},{4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[9] t[11] t[11] t[11] t[11];
   t[8] t[2] t[9] t[11] t[11] t[11] t[11];
   t[9] t[9] t[3] t[11] t[11] t[11] t[11];
   t[11] t[11] t[11] t[4] t[10] t[10] t[11];
   t[11] t[11] t[11] t[10] t[5] t[10] t[11];
   t[11] t[11] t[11] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3},{1,2,3,4}",orbit_size=420,
   tree=[t[1] t[8] t[9] t[10] t[11] t[11] t[11];
   t[8] t[2] t[9] t[10] t[11] t[11] t[11];
   t[9] t[9] t[3] t[10] t[11] t[11] t[11];
   t[10] t[10] t[10] t[4] t[11] t[11] t[11];
   t[11] t[11] t[11] t[11] t[5] t[11] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3},{1,2,3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[9] t[10] t[10] t[10] t[11];
   t[8] t[2] t[9] t[10] t[10] t[10] t[11];
   t[9] t[9] t[3] t[10] t[10] t[10] t[11];
   t[10] t[10] t[10] t[4] t[10] t[10] t[11];
   t[10] t[10] t[10] t[10] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{3,4,5},{3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[11] t[11] t[11] t[11] t[11];
   t[8] t[2] t[11] t[11] t[11] t[11] t[11];
   t[11] t[11] t[3] t[9] t[9] t[10] t[11];
   t[11] t[11] t[9] t[4] t[9] t[10] t[11];
   t[11] t[11] t[9] t[9] t[5] t[10] t[11];
   t[11] t[11] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{3,4,5},{1,2,3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[10] t[10] t[10] t[10] t[11];
   t[8] t[2] t[10] t[10] t[10] t[10] t[11];
   t[10] t[10] t[3] t[9] t[9] t[10] t[11];
   t[10] t[10] t[9] t[4] t[9] t[10] t[11];
   t[10] t[10] t[9] t[9] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3,4,5},{1,2,3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[9] t[9] t[9] t[10] t[11];
   t[8] t[2] t[9] t[9] t[9] t[10] t[11];
   t[9] t[9] t[3] t[9] t[9] t[10] t[11];
   t[9] t[9] t[9] t[4] t[9] t[10] t[11];
   t[9] t[9] t[9] t[9] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2,3},{1,2,3,4},{1,2,3,4,5}",orbit_size=420,
   tree=[t[1] t[8] t[8] t[9] t[10] t[11] t[11];
   t[8] t[2] t[8] t[9] t[10] t[11] t[11];
   t[8] t[8] t[3] t[9] t[10] t[11] t[11];
   t[9] t[9] t[9] t[4] t[10] t[11] t[11];
   t[10] t[10] t[10] t[10] t[5] t[11] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2,3},{1,2,3,4},{1,2,3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[8] t[9] t[10] t[10] t[11];
   t[8] t[2] t[8] t[9] t[10] t[10] t[11];
   t[8] t[8] t[3] t[9] t[10] t[10] t[11];
   t[9] t[9] t[9] t[4] t[10] t[10] t[11];
   t[10] t[10] t[10] t[10] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2,3},{1,2,3,4,5},{1,2,3,4,5,6}",orbit_size=420,
   tree=[t[1] t[8] t[8] t[9] t[9] t[10] t[11];
   t[8] t[2] t[8] t[9] t[9] t[10] t[11];
   t[8] t[8] t[3] t[9] t[9] t[10] t[11];
   t[9] t[9] t[9] t[4] t[9] t[10] t[11];
   t[9] t[9] t[9] t[9] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{3,4},{1,2,5}",orbit_size=630,
   tree=[t[1] t[8] t[11] t[11] t[10] t[11] t[11];
   t[8] t[2] t[11] t[11] t[10] t[11] t[11];
   t[11] t[11] t[3] t[9] t[11] t[11] t[11];
   t[11] t[11] t[9] t[4] t[11] t[11] t[11];
   t[10] t[10] t[11] t[11] t[5] t[11] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{3,4},{1,2,5,6}",orbit_size=630,
   tree=[t[1] t[8] t[11] t[11] t[10] t[10] t[11];
   t[8] t[2] t[11] t[11] t[10] t[10] t[11];
   t[11] t[11] t[3] t[9] t[11] t[11] t[11];
   t[11] t[11] t[9] t[4] t[11] t[11] t[11];
   t[10] t[10] t[11] t[11] t[5] t[10] t[11];
   t[10] t[10] t[11] t[11] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3},{1,2,3,4,5}",orbit_size=630,
   tree=[t[1] t[8] t[9] t[10] t[10] t[11] t[11];
   t[8] t[2] t[9] t[10] t[10] t[11] t[11];
   t[9] t[9] t[3] t[10] t[10] t[11] t[11];
   t[10] t[10] t[10] t[4] t[10] t[11] t[11];
   t[10] t[10] t[10] t[10] t[5] t[11] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3,4},{1,2,3,4,5}",orbit_size=630,
   tree=[t[1] t[8] t[9] t[9] t[10] t[11] t[11];
   t[8] t[2] t[9] t[9] t[10] t[11] t[11];
   t[9] t[9] t[3] t[9] t[10] t[11] t[11];
   t[9] t[9] t[9] t[4] t[10] t[11] t[11];
   t[10] t[10] t[10] t[10] t[5] t[11] t[11];
   t[11] t[11] t[11] t[11] t[11] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]]),

  (n=7,
   id="{1,2},{1,2,3,4},{1,2,3,4,5,6}",orbit_size=630,
   tree=[t[1] t[8] t[9] t[9] t[10] t[10] t[11];
   t[8] t[2] t[9] t[9] t[10] t[10] t[11];
   t[9] t[9] t[3] t[9] t[10] t[10] t[11];
   t[9] t[9] t[9] t[4] t[10] t[10] t[11];
   t[10] t[10] t[10] t[10] t[5] t[10] t[11];
   t[10] t[10] t[10] t[10] t[10] t[6] t[11];
   t[11] t[11] t[11] t[11] t[11] t[11] t[7]])
  ]
end


text:
# ---
# title: 415. Add Strings
# id: problem415
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/add-strings/description/>
# hidden: true
# ---
#
# Given two non-negative integers `num1` and `num2` represented as string,
# return the sum of `num1` and `num2`.
#
# **Note:**
#
#   1. The length of both `num1` and `num2` is < 5100.
#   2. Both `num1` and `num2` contains only digits `0-9`.
#   3. Both `num1` and `num2` does not contain any leading zero.
#   4. You **must not use any built-in BigInteger library** or **convert the inputs to integer** directly.
#
#
## @lc code=start
using LeetCode

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


text:
using DFTShims
using Base.Test
using Unitful
using AxisArrays

const DD = DFTShims.Dispatch.Dimensions
const DH = DFTShims.Dispatch.Hartree
const UH = DFTShims.UnitfulHartree

@testset "Dispatch" begin
    include("dispatch.jl")
end
@testset "Traits" begin
    include("Traits.jl")
end
@testset "Initialization" begin
    include("Initialization.jl")
end
@testset "Constant arrays" begin
    include("ConstantArrays.jl")
end
nothing

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims