text:
### abstract type hierarchy

abstract type Basis end
abstract type Basislike <: Basis end

# the basis in which derivatives are sparse for a given field
struct DerivBasis <: Basislike end
const Ð! = DerivBasis
const Ð  = DerivBasis

# the basis in which lensing is a pixel remapping for a given field
struct LenseBasis <: Basislike end
const Ł! = LenseBasis
const Ł  = LenseBasis

# the nearest harmonic basis, e.g. for anything QU its QUFourier and
# for anything EB its EBFourier
abstract type HarmonicBasis <: Basislike end

## fields
abstract type Field{B<:Basis,T} <: AbstractVector{T} end
abstract type ImplicitField{B<:Basis,T} <: Field{B,T} end

## linear operators
abstract type ImplicitOp{T} <: AbstractMatrix{T} end
const DiagOp{F<:Field,T} = Diagonal{T,F}
# includes Diagonal and Adjoint wrappers, which can't be made <:ImplicitOp directly
const FieldOp{T} = Union{ImplicitOp{T},Adjoint{T,<:ImplicitOp{T}},Diagonal{T,<:Field{<:Any,T}}}

# scalars which are allowed in our expressions must be real because we
# implicitly assume our maps are real, and addition/multiplication by a complex
# number, even of the fourier transform, would break this.
const Scalar = Real

# other useful unions
const FieldOrOp   = Union{Field,FieldOp}
const FieldOpScal = Union{Field,FieldOp,Scalar}

## basis types
abstract type BasisProd{Bs} <: Basis end

struct Basis2Prod{B₁,B₂}    <: BasisProd{Tuple{B₁,B₂}} end
struct Basis3Prod{B₁,B₂,B₃} <: BasisProd{Tuple{B₁,B₂,B₃}} end

abstract type S0Basis <: Basis end
abstract type PolBasis <: Basis end

struct Map     <: S0Basis end
struct Fourier <: S0Basis end
struct 𝐈       <: PolBasis end
struct 𝐐𝐔      <: PolBasis end
struct 𝐄𝐁      <: PolBasis end

const QUMap      = Basis2Prod{    𝐐𝐔, Map     }
const EBMap      = Basis2Prod{    𝐄𝐁, Map     }
const QUFourier  = Basis2Prod{    𝐐𝐔, Fourier }
const EBFourier  = Basis2Prod{    𝐄𝐁, Fourier }
const IQUMap     = Basis3Prod{ 𝐈, 𝐐𝐔, Map     }
const IEBMap     = Basis3Prod{ 𝐈, 𝐄𝐁, Map     }
const IQUFourier = Basis3Prod{ 𝐈, 𝐐𝐔, Fourier }
const IEBFourier = Basis3Prod{ 𝐈, 𝐄𝐁, Fourier }

# handy for picking out anything Map/Fourier
const SpatialBasis{B,I,P} = Union{B, Basis2Prod{P,B}, Basis3Prod{I,P,B}}

## applying and converting bases

LenseBasis(                            ::S0Basis)  = Map
LenseBasis(::Basis2Prod{   <:PolBasis, <:S0Basis}) = QUMap
LenseBasis(::Basis3Prod{𝐈, <:PolBasis, <:S0Basis}) = IQUMap
DerivBasis(                            ::S0Basis)  = Fourier
DerivBasis(::Basis2Prod{   <:PolBasis, <:S0Basis}) = QUFourier
DerivBasis(::Basis3Prod{𝐈, <:PolBasis, <:S0Basis}) = IQUFourier
HarmonicBasis(                         ::S0Basis)  = Fourier
HarmonicBasis(::Basis2Prod{𝐐𝐔,         <:S0Basis}) = QUFourier
HarmonicBasis(::Basis2Prod{𝐄𝐁,         <:S0Basis}) = EBFourier

# B(::Basis) converts the basis type, e.g. Map(QUFourier()) = QUMap
(::Type{B})(::B′) where {B<:S0Basis, B′<:S0Basis} = B
(::Type{B})(::Basis2Prod{  Pol,B′}) where {Pol, B<:S0Basis, B′<:S0Basis} = Basis2Prod{  Pol,B}
(::Type{B})(::Basis3Prod{𝐈,Pol,B′}) where {Pol, B<:S0Basis, B′<:S0Basis} = Basis3Prod{𝐈,Pol,B}

# B(f::Field) where B is a basis converts f to that basis. This is the
# fallback if the field is already in the right basis.
(::Type{B})(f::Field{B}) where {B<:Basis} = f
# This is the fallback if no conversion is explicilty defined, which
# tries to first convert the basis, e.g. Map(f::BaseQUFourier) becomes
# Map(QUFourier())(f) which calls QUMap(f). This is also used for
# Basislike objects like LenseBasis(...)
(::Type{B})(f::Field{B′}) where {B′<:Basis,B<:Basis} = B(B′())(f)

# In-place version of the above
(::Type{B})(f′::Field, f::Field{B′}) where {B′<:Basis,B<:Basis} = B(B′())(f′,f)
# And the case where its already in the right basis (but note, we
# never actually set f′ in this case, which is more efficient, but
# necessitates some care when using this construct)
(::Type{B})(f′::Field{B}, f::Field{B}) where {B<:Basis} = f

# Basis conversion automatically maps over arrays
(::Type{B})(a::AbstractArray{<:Field}...) where {B<:Basis} = B.(a...)

# The abstract `Basis` type means "any basis", hence this conversion rule:
Basis(f::Field) = f

# used in make_field_aliases below
basis_aliases = OrderedDict(
    "Map"        => Map,
    "Fourier"    => Fourier,
    "QUMap"      => QUMap,
    "QUFourier"  => QUFourier,
    "EBMap"      => EBMap,
    "EBFourier"  => EBFourier,
    "IQUMap"     => IQUMap,
    "IQUFourier" => IQUFourier,
    "IEBMap"     => IEBMap,
    "IEBFourier" => IEBFourier,
    "S0"         => S0Basis,
    "QU"         => Basis2Prod{   𝐐𝐔        , <:S0Basis},
    "EB"         => Basis2Prod{   𝐄𝐁        , <:S0Basis},
    "S2Map"      => Basis2Prod{   <:PolBasis, Map},
    "S2Fourier"  => Basis2Prod{   <:PolBasis, Fourier},
    "S2"         => Basis2Prod{   <:PolBasis, <:S0Basis},
    "IQU"        => Basis3Prod{𝐈, 𝐐𝐔        , <:S0Basis},
    "IEB"        => Basis3Prod{𝐈, 𝐄𝐁        , <:S0Basis},
    "S02Map"     => Basis3Prod{𝐈, <:PolBasis, Map},
    "S02Fourier" => Basis3Prod{𝐈, <:PolBasis, Fourier},
    "S02"        => Basis3Prod{𝐈, <:PolBasis, <:S0Basis},
    "Field"      => Any
)

# Enumerates all the fields types like FlatMap, FlatFourier, etc...
# for a field_name like "Flat" and a bound on the M type parameter in
# BaseField{B,M,T,A}. Note: the seemingly-redundant <:AbstractArray{T}
# in the argument (which is enforced in BaseField anyway) is there to
# help prevent method ambiguities
function make_field_aliases(field_root, M_bound; export_names=true, extra_aliases=Dict())
    for (basis_alias, B) in merge(basis_aliases, extra_aliases)
        F = Symbol(field_root, basis_alias)
        if isconcretetype(B)
            @eval const $F{       M<:$M_bound, T, A<:AbstractArray{T}} = BaseField{$B, M, T, A}
        else
            @eval const $F{B<:$B, M<:$M_bound, T, A<:AbstractArray{T}} = BaseField{B,  M, T, A}
        end
        if export_names
            @eval export $F
        end
    end
end

# for printing
for (alias,B) in basis_aliases
    if isconcretetype(B)
        @eval typealias(::Type{$B}) = $alias
    end
end

## generic promotion rules
basis_promotion_rules = Dict(
    # spin-0
    (Map,        Fourier)     => Map,
    # spin-2
    (QUMap,      QUFourier)   => QUMap,
    (EBMap,      EBFourier)   => EBFourier,
    (QUMap,      EBMap)       => QUMap,
    (QUFourier,  EBFourier)   => QUFourier,
    (QUMap,      EBFourier)   => QUMap,
    (QUFourier,  EBMap)       => QUFourier,
    # spin-(0,2)
    (IQUMap,     IQUFourier)  => IQUMap,
    (IEBMap,     IEBFourier)  => IEBFourier,
    (IQUMap,     IEBMap)      => IQUMap,
    (IQUFourier, IEBFourier)  => IQUFourier,
    (IQUMap,     IEBFourier)  => IQUMap,
    (IQUFourier, IEBMap)      => IQUFourier
)
for ((B₁,B₂),B) in basis_promotion_rules
    @eval promote_basis_generic_rule(::$B₁, ::$B₂) = $B()
end
promote_basis_generic_rule(::B, ::B) where {B<:Basis} = B()
promote_basis_generic_rule(::Any, ::Any) = Unknown()
promote_basis_generic(x, y) = select_known_rule(promote_basis_generic_rule, x, y)
unknown_rule_error(::typeof(promote_basis_generic_rule), ::B₁, ::B₂) where {B₁, B₂} =
    error("Can't promote fields in $(typealias(B₁)) and $(typealias(B₂)) bases.")

## stricter rules used only in broadcasting
promote_basis_strict_rule(::B,   ::B )        where {B <:Basis}                                                    = B()
promote_basis_strict_rule(::B,   ::Basislike) where {B <:Basis}                                                    = B()
promote_basis_strict_rule(::B₂,  ::B₀)        where {B₀<:Union{Map,Fourier}, B₂ <:Basis2Prod{  <:Union{𝐐𝐔,𝐄𝐁},B₀}} = B₂()
promote_basis_strict_rule(::B₀₂, ::B₀)        where {B₀<:Union{Map,Fourier}, B₀₂<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁},B₀}} = B₀₂()
promote_basis_strict_rule(::Any, ::Any)                                                                            = Unknown()
promote_basis_strict(x, y) = select_known_rule(promote_basis_strict_rule, x, y)
unknown_rule_error(::typeof(promote_basis_strict_rule), ::B₁, ::B₂) where {B₁, B₂} =
    error("Can't broadcast fields in $(typealias(B₁)) and $(typealias(B₂)) bases.")

# convenience "getter" functions for the B basis type parameter
basis(f::F) where {F<:Field} = basis(F)
basis(::Type{<:Field{B}}) where {B<:Basis} = B
basis(::Type{<:Field}) = Basis
basis(::AbstractVector) = Basis

### printing
typealias(::Type{B}) where {B<:Basis} = string(B)
Base.show_datatype(io::IO, t::Type{<:Union{Field,FieldOp}}) = print(io, typealias(t))
Base.isempty(::ImplicitOp) = true
Base.isempty(::ImplicitField) = true
Base.isempty(::Diagonal{<:Any,<:ImplicitField}) = true
function Base.summary(io::IO, x::FieldOp)
    try
        print(io, join(size(x), "×"), " ")
    catch err
        @assert err isa MethodError
        print(io, "⍰×⍰ ")
    end
    Base.showarg(io, x, true)
end
function Base.summary(io::IO, x::ImplicitField)
    try
        print(io, size(x,1), "-element ")
    catch err
        @assert err isa MethodError
        print(io, "⍰-element ")
    end
    Base.showarg(io, x, true)
end

### logdet

"""
    logdet(L::FieldOp, θ)

If `L` depends on `θ`, evaluates `logdet(L(θ))` offset by its fiducial value at
`L()`. Otherwise, returns 0.
"""
logdet(L::FieldOp, θ) = depends_on(L,θ) ? logdet(L()\L(θ)) : 0
logdet(L::Int, θ) = 0 # the default returns Float64 which unwantedly poisons the backprop to Float64
logdet(L, θ) = logdet(L)

### Simulation

@doc doc"""
    simulate(Σ; rng=global_rng_for(Σ), seed=nothing)

Draw a simulation from the covariance matrix `Σ`, i.e. draw a random vector
$\xi$ such that the covariance $\langle \xi \xi^\dagger \rangle = \Sigma$.

The random number generator `rng` will be used and advanced in the proccess, and
is by default the appropriate one depending on if `Σ` is backed by `Array` or
`CuArray`.

The `seed` argument can also be used to seed the `rng`.
"""
function simulate(Σ; rng=global_rng_for(Σ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    simulate(rng, Σ; kwargs...)
end
function white_noise(ξ; rng=global_rng_for(ξ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    white_noise(ξ, rng; kwargs...)
end
function fixed_white_noise(ξ; rng=global_rng_for(Σ), seed=nothing, kwargs...)
    isnothing(seed) || seed!(rng, seed)
    fixed_white_noise(ξ, rng; kwargs...)
end

global_rng_for(x::T) where {T<:AbstractArray} = global_rng_for(T)
global_rng_for(::Type{<:Array}) = Random.default_rng()

### spin adjoints
# represents a field which is adjoint over just the "spin" indices.
# multiplying such a field by a non-adjointed one should be the inner
# product over just the spin indices, hence return a spin-0 field,
# rather than a scalar. note: these are really only lightly used in
# one place in LenseFlow, so they have almost no real functionality
# for now.
struct SpinAdjoint{F<:Field}
    f :: F
end
spin_adjoint(f::Field) = SpinAdjoint(f)

### misc

function ≈(x::Field{<:Any,T}, y::Field{<:Any,S}; atol=0, rtol=Base.rtoldefault(T,S,atol)) where {T,S}
    _norm(x) = norm(unbatch(norm(x)))
    _norm(x - y) <= max(atol, rtol*max(_norm(x), _norm(y)))
end

# allows L(θ) or L(ϕ) to work and be a no-op for things which don't
# depend on parameters or a field
(L::Union{FieldOp,UniformScaling})(::Union{Field,NamedTuple}) = L

# allow using `I` as a lensing operator to represent no lensing
alloc_cache(x, ::Any) = x
cache(x, ::Any) = x
cache!(x, ::Any) = x

# caching for adjoints
cache(L::Adjoint, f) = cache(L',f)'
cache!(L::Adjoint, f) = cache!(L',f)'

# todo: fix this
*(::UniformScaling{Bool}, L::FieldOp) = L

# we use Field cat'ing mainly for plotting, e.g. plot([f f; f f]) plots a 2×2
# matrix of maps. the following definitions make it so that Fields aren't
# splatted into a giant matrix when doing [f f; f f] (which they would othewise
# be since they're Arrays)
hvcat(rows::Tuple{Vararg{Int}}, values::Field...) = hvcat(rows, ([x] for x in values)...)
hcat(values::Field...) = hcat(([x] for x in values)...)

### printing
print_array(io::IO, f::Field) = !isempty(f) && print_array(io, f[:])
show_vector(io::IO, f::Field) = !isempty(f) && show_vector(io, f[:])
Base.has_offset_axes(::Field) = false # needed for Diagonal(::Field) if the Field is implicitly-sized

# addition/subtraction works between any fields and scalars, promotion is done
# automatically if fields are in different bases
for op in (:+,:-), (T1,T2,promote) in ((:Field,:Scalar,false),(:Scalar,:Field,false),(:Field,:Field,true))
    @eval ($op)(a::$T1, b::$T2) = broadcast($op, ($promote ? promote(a,b) : (a,b))...)
end

# multiplication/division is not strictly defined for abstract vectors, but
# make it work anyway if the two fields are exactly the same type, in which case
# its clear we wanted broadcasted multiplication/division.
for op in (:*, :/)
    @eval ($op)(a::Field{B}, b::Field{B}) where {B} = broadcast($op, a, b)
end

# package code should use batch_index(f, 1), but for interactive work its
# convenient to be able to do f[!,1]
getindex(x::Union{Real,Field,FieldOp}, ::typeof(!), I) = batch_index(x, I)

# adapt_structure(to, x::Union{ImplicitOp,ImplicitField}) = x
# adapt_structure(to, arr::Array{<:Field}) = map(f -> adapt(to,f), arr)

one(f::Field) = fill!(similar(f), one(eltype(f)))
norm(f::Field) = sqrt(dot(f,f))
# sum_kbn(f::Field) = sum_kbn(f[:])

@init @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    # never try to auto-convert Fields or FieldOps to Python arrays
    PyCall.PyObject(x::Union{FieldOp,Field}) = PyCall.pyjlwrap_new(x)
end


text:
cₚ = 1004.703 # Specific Heat at constant pressure
kₜ = 0.0246295028571 #Thermal conductivity
k_cyl = kₜ * 4
kᵦ = 1.38064852e-23 # Boltzmann constant (m² kg/(s² K))

density = 0.000210322
R = kᵦ * 6.0221409e23 # kg⋅m²/(s²*K*mol)
mol_mass = 28.96 # g/mol
μ = 28.96 / 6.0221409e23 # mean molecular mass (g)
R₀ = R / (mol_mass / 1000)

# Heat diffusion constant in fluid
k₁ = kₜ / (density * cₚ)

# Heat diffusion constant in cylinder
k₂ = k_cyl / (density * cₚ)

kᵨ = 1e-3

ν = 1.716e-5 # 0.0005081150545826

kᵥ = ν# / density

γ = 1 + R₀ / (cₚ - R₀) # Mayer's formula + definition of adiabatic constant
println("γ: $γ")
# T = (γ - 1) * (mol_mass / 1000) / R₀ * e
#e2t = (μ/1000) / kᵦ
#e2t = (γ - 1) * mol_mass / R₀
e2t = 1/(cₚ * density)
@show e2t
e2p = e2t * R₀
#e2p = (γ - 1)
t2e = 1/e2t

ext: jl

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
@testset "8.4 Trig integral functions" begin
    (a, b, c, d, e, m, n, x, ) = @variables a b c d e m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Trig*Integral*Functions=#

    #= ::Section::Closed:: =#
    #=Sine*integral*function=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*SinIntegral(b*x)^n=#

    @test_int [x^m*SinIntegral(b*x), x, 5, (x^m*Gamma(1 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b*(1 + m))) + (x^m*Gamma(1 + m, im*b*x))/((im*b*x)^m*(2*b*(1 + m))) + (x^(1 + m)*SinIntegral(b*x))/(1 + m)]

    @test_int [x^3*SinIntegral(b*x), x, 6, -((3*x*cos(b*x))/(2*b^3)) + (x^3*cos(b*x))/(4*b) + (3*sin(b*x))/(2*b^4) - (3*x^2*sin(b*x))/(4*b^2) + (1/4)*x^4*SinIntegral(b*x)]
    @test_int [x^2*SinIntegral(b*x), x, 5, -((2*cos(b*x))/(3*b^3)) + (x^2*cos(b*x))/(3*b) - (2*x*sin(b*x))/(3*b^2) + (1/3)*x^3*SinIntegral(b*x)]
    @test_int [x^1*SinIntegral(b*x), x, 4, (x*cos(b*x))/(2*b) - sin(b*x)/(2*b^2) + (1/2)*x^2*SinIntegral(b*x)]
    @test_int [x^0*SinIntegral(b*x), x, 1, cos(b*x)/b + x*SinIntegral(b*x)]
    @test_int [SinIntegral(b*x)/x^1, x, 1, (1/2)*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-im)*b*x) + (1/2)*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], im*b*x)]
    @test_int [SinIntegral(b*x)/x^2, x, 4, b*CosIntegral(b*x) - sin(b*x)/x - SinIntegral(b*x)/x]
    @test_int [SinIntegral(b*x)/x^3, x, 5, -((b*cos(b*x))/(4*x)) - sin(b*x)/(4*x^2) - (1/4)*b^2*SinIntegral(b*x) - SinIntegral(b*x)/(2*x^2)]

    @test_int [x^m*SinIntegral(b*x)^2, x, 0, CannotIntegrate(x^m*SinIntegral(b*x)^2, x)]

    @test_int [x^3*SinIntegral(b*x)^2, x, 19, x^2/(2*b^2) + (3*CosIntegral(2*b*x))/(2*b^4) - (3*log(x))/(2*b^4) - (x*cos(b*x)*sin(b*x))/b^3 + (2*sin(b*x)^2)/b^4 - (x^2*sin(b*x)^2)/(4*b^2) - (3*x*cos(b*x)*SinIntegral(b*x))/b^3 + (x^3*cos(b*x)*SinIntegral(b*x))/(2*b) + (3*sin(b*x)*SinIntegral(b*x))/b^4 - (3*x^2*sin(b*x)*SinIntegral(b*x))/(2*b^2) + (1/4)*x^4*SinIntegral(b*x)^2]
    @test_int [x^2*SinIntegral(b*x)^2, x, 15, (5*x)/(6*b^2) - (5*cos(b*x)*sin(b*x))/(6*b^3) - (x*sin(b*x)^2)/(3*b^2) - (4*cos(b*x)*SinIntegral(b*x))/(3*b^3) + (2*x^2*cos(b*x)*SinIntegral(b*x))/(3*b) - (4*x*sin(b*x)*SinIntegral(b*x))/(3*b^2) + (1/3)*x^3*SinIntegral(b*x)^2 + (2*SinIntegral(2*b*x))/(3*b^3)]
    @test_int [x^1*SinIntegral(b*x)^2, x, 10, -(CosIntegral(2*b*x)/(2*b^2)) + log(x)/(2*b^2) - sin(b*x)^2/(2*b^2) + (x*cos(b*x)*SinIntegral(b*x))/b - (sin(b*x)*SinIntegral(b*x))/b^2 + (1/2)*x^2*SinIntegral(b*x)^2]
    @test_int [x^0*SinIntegral(b*x)^2, x, 6, (2*cos(b*x)*SinIntegral(b*x))/b + x*SinIntegral(b*x)^2 - SinIntegral(2*b*x)/b]
    @test_int [SinIntegral(b*x)^2/x^1, x, 0, CannotIntegrate(SinIntegral(b*x)^2/x, x)]
    @test_int [SinIntegral(b*x)^2/x^2, x, 0, CannotIntegrate(SinIntegral(b*x)^2/x^2, x)]
    @test_int [SinIntegral(b*x)^2/x^3, x, 0, CannotIntegrate(SinIntegral(b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*SinIntegral(a+b*x)^n=#

    @test_int [x^m*SinIntegral(a + b*x), x, 1, -((b*CannotIntegrate((x^(1 + m)*sin(a + b*x))/(a + b*x), x))/(1 + m)) + (x^(1 + m)*SinIntegral(a + b*x))/(1 + m)]

    @test_int [x^3*SinIntegral(a + b*x), x, 14, (a*cos(a + b*x))/(2*b^4) - (a^3*cos(a + b*x))/(4*b^4) - (3*x*cos(a + b*x))/(2*b^3) + (a^2*x*cos(a + b*x))/(4*b^3) - (a*x^2*cos(a + b*x))/(4*b^2) + (x^3*cos(a + b*x))/(4*b) + (3*sin(a + b*x))/(2*b^4) - (a^2*sin(a + b*x))/(4*b^4) + (a*x*sin(a + b*x))/(2*b^3) - (3*x^2*sin(a + b*x))/(4*b^2) - (a^4*SinIntegral(a + b*x))/(4*b^4) + (1/4)*x^4*SinIntegral(a + b*x)]
    @test_int [x^2*SinIntegral(a + b*x), x, 10, -((2*cos(a + b*x))/(3*b^3)) + (a^2*cos(a + b*x))/(3*b^3) - (a*x*cos(a + b*x))/(3*b^2) + (x^2*cos(a + b*x))/(3*b) + (a*sin(a + b*x))/(3*b^3) - (2*x*sin(a + b*x))/(3*b^2) + (a^3*SinIntegral(a + b*x))/(3*b^3) + (1/3)*x^3*SinIntegral(a + b*x)]
    @test_int [x^1*SinIntegral(a + b*x), x, 7, -((a*cos(a + b*x))/(2*b^2)) + (x*cos(a + b*x))/(2*b) - sin(a + b*x)/(2*b^2) - (a^2*SinIntegral(a + b*x))/(2*b^2) + (1/2)*x^2*SinIntegral(a + b*x)]
    @test_int [x^0*SinIntegral(a + b*x), x, 1, cos(a + b*x)/b + ((a + b*x)*SinIntegral(a + b*x))/b]
    @test_int [SinIntegral(a + b*x)/x^1, x, 0, CannotIntegrate(SinIntegral(a + b*x)/x, x)]
    @test_int [SinIntegral(a + b*x)/x^2, x, 7, (b*CosIntegral(b*x)*sin(a))/a + (b*cos(a)*SinIntegral(b*x))/a - (b*SinIntegral(a + b*x))/a - SinIntegral(a + b*x)/x]
    @test_int [SinIntegral(a + b*x)/x^3, x, 11, (b^2*cos(a)*CosIntegral(b*x))/(2*a) - (b^2*CosIntegral(b*x)*sin(a))/(2*a^2) - (b*sin(a + b*x))/(2*a*x) - (b^2*cos(a)*SinIntegral(b*x))/(2*a^2) - (b^2*sin(a)*SinIntegral(b*x))/(2*a) + (b^2*SinIntegral(a + b*x))/(2*a^2) - SinIntegral(a + b*x)/(2*x^2)]

    @test_int [x^m*SinIntegral(a + b*x)^2, x, 0, CannotIntegrate(x^m*SinIntegral(a + b*x)^2, x)]

    @test_int [x^2*SinIntegral(a + b*x)^2, x, 39, (2*x)/(3*b^2) - (a*cos(2*a + 2*b*x))/(3*b^3) + (x*cos(2*a + 2*b*x))/(6*b^2) + (a*CosIntegral(2*a + 2*b*x))/b^3 - (a*log(a + b*x))/b^3 - (2*cos(a + b*x)*sin(a + b*x))/(3*b^3) - sin(2*a + 2*b*x)/(12*b^3) - (4*cos(a + b*x)*SinIntegral(a + b*x))/(3*b^3) + (2*a^2*cos(a + b*x)*SinIntegral(a + b*x))/(3*b^3) - (2*a*x*cos(a + b*x)*SinIntegral(a + b*x))/(3*b^2) + (2*x^2*cos(a + b*x)*SinIntegral(a + b*x))/(3*b) + (2*a*sin(a + b*x)*SinIntegral(a + b*x))/(3*b^3) - (4*x*sin(a + b*x)*SinIntegral(a + b*x))/(3*b^2) + (a^2*(a + b*x)*SinIntegral(a + b*x)^2)/(3*b^3) - (a*x*(a + b*x)*SinIntegral(a + b*x)^2)/(3*b^2) + (x^2*(a + b*x)*SinIntegral(a + b*x)^2)/(3*b) + (2*SinIntegral(2*a + 2*b*x))/(3*b^3) - (a^2*SinIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*SinIntegral(a + b*x)^2, x, 17, cos(2*a + 2*b*x)/(4*b^2) - CosIntegral(2*a + 2*b*x)/(2*b^2) + log(a + b*x)/(2*b^2) - (a*cos(a + b*x)*SinIntegral(a + b*x))/b^2 + (x*cos(a + b*x)*SinIntegral(a + b*x))/b - (sin(a + b*x)*SinIntegral(a + b*x))/b^2 - (a*(a + b*x)*SinIntegral(a + b*x)^2)/(2*b^2) + (x*(a + b*x)*SinIntegral(a + b*x)^2)/(2*b) + (a*SinIntegral(2*a + 2*b*x))/b^2]
    @test_int [x^0*SinIntegral(a + b*x)^2, x, 5, (2*cos(a + b*x)*SinIntegral(a + b*x))/b + ((a + b*x)*SinIntegral(a + b*x)^2)/b - SinIntegral(2*a + 2*b*x)/b]
    @test_int [SinIntegral(a + b*x)^2/x^1, x, 0, CannotIntegrate(SinIntegral(a + b*x)^2/x, x)]
    @test_int [SinIntegral(a + b*x)^2/x^2, x, 0, CannotIntegrate(SinIntegral(a + b*x)^2/x^2, x)]
    @test_int [SinIntegral(a + b*x)^2/x^3, x, 0, CannotIntegrate(SinIntegral(a + b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*SinIntegral(d*(a+b*log(c*x^n)))=#

    @test_int [x^2*SinIntegral(d*(a + b*log(c*x^n))), x, 7, ((-(1/6))*im*x^3*ExpIntegralEi(((3 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) + ((1/6)*im*x^3*ExpIntegralEi(((3 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) + (1/3)*x^3*SinIntegral(d*(a + b*log(c*x^n)))]
    @test_int [x^1*SinIntegral(d*(a + b*log(c*x^n))), x, 7, ((-(1/4))*im*x^2*ExpIntegralEi(((2 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) + ((1/4)*im*x^2*ExpIntegralEi(((2 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) + (1/2)*x^2*SinIntegral(d*(a + b*log(c*x^n)))]
    @test_int [x^0*SinIntegral(d*(a + b*log(c*x^n))), x, 7, ((-(1/2))*im*x*ExpIntegralEi(((1 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^(a/(b*n))*(c*x^n)^n^(-1)) + ((1/2)*im*x*ExpIntegralEi(((1 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^(a/(b*n))*(c*x^n)^n^(-1)) + x*SinIntegral(d*(a + b*log(c*x^n)))]
    @test_int [SinIntegral(d*(a + b*log(c*x^n)))/x^1, x, 3, cos(d*(a + b*log(c*x^n)))/(b*d*n) + ((a + b*log(c*x^n))*SinIntegral(d*(a + b*log(c*x^n))))/(b*n)]
    @test_int [SinIntegral(d*(a + b*log(c*x^n)))/x^2, x, 7, -((im*ℯ^(a/(b*n))*(c*x^n)^(1/n)*ExpIntegralEi(-(((1 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x)) + (im*ℯ^(a/(b*n))*(c*x^n)^(1/n)*ExpIntegralEi(-(((1 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x) - SinIntegral(d*(a + b*log(c*x^n)))/x]
    @test_int [SinIntegral(d*(a + b*log(c*x^n)))/x^3, x, 7, -((im*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2)) + (im*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2) - SinIntegral(d*(a + b*log(c*x^n)))/(2*x^2)]

    @test_int [(e*x)^m*SinIntegral(d*(a + b*log(c*x^n))), x, 7, -((im*x*(e*x)^m*ExpIntegralEi(((1 + m - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((a*(1 + m))/(b*n))*(c*x^n)^((1 + m)/n)*(2*(1 + m)))) + (im*x*(e*x)^m*ExpIntegralEi(((1 + m + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(ℯ^((a*(1 + m))/(b*n))*(c*x^n)^((1 + m)/n)*(2*(1 + m))) + ((e*x)^(1 + m)*SinIntegral(d*(a + b*log(c*x^n))))/(e*(1 + m))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(b*x)*SinIntegral(b*x)^n=#

    @test_int [sin(b*x)*SinIntegral(b*x)/x^3, x, 14, b^2*CosIntegral(2*b*x) - (b*cos(b*x)*sin(b*x))/(2*x) - sin(b*x)^2/(4*x^2) - (b*sin(2*b*x))/(4*x) - (b*cos(b*x)*SinIntegral(b*x))/(2*x) - (sin(b*x)*SinIntegral(b*x))/(2*x^2) - (1/4)*b^2*SinIntegral(b*x)^2]
    @test_int [sin(b*x)*SinIntegral(b*x)/x^2, x, 5, b*CannotIntegrate((cos(b*x)*SinIntegral(b*x))/x, x) - sin(b*x)^2/x - (sin(b*x)*SinIntegral(b*x))/x + b*SinIntegral(2*b*x)]
    @test_int [sin(b*x)*SinIntegral(b*x)/x, x, 1, (1/2)*SinIntegral(b*x)^2]
    @test_int [sin(b*x)*SinIntegral(b*x), x, 5, -((cos(b*x)*SinIntegral(b*x))/b) + SinIntegral(2*b*x)/(2*b)]
    @test_int [x*sin(b*x)*SinIntegral(b*x), x, 9, CosIntegral(2*b*x)/(2*b^2) - log(x)/(2*b^2) + sin(b*x)^2/(2*b^2) - (x*cos(b*x)*SinIntegral(b*x))/b + (sin(b*x)*SinIntegral(b*x))/b^2]
    @test_int [x^2*sin(b*x)*SinIntegral(b*x), x, 14, -((5*x)/(4*b^2)) + (5*cos(b*x)*sin(b*x))/(4*b^3) + (x*sin(b*x)^2)/(2*b^2) + (2*cos(b*x)*SinIntegral(b*x))/b^3 - (x^2*cos(b*x)*SinIntegral(b*x))/b + (2*x*sin(b*x)*SinIntegral(b*x))/b^2 - SinIntegral(2*b*x)/b^3]
    @test_int [x^3*sin(b*x)*SinIntegral(b*x), x, 18, -(x^2/b^2) - (3*CosIntegral(2*b*x))/b^4 + (3*log(x))/b^4 + (2*x*cos(b*x)*sin(b*x))/b^3 - (4*sin(b*x)^2)/b^4 + (x^2*sin(b*x)^2)/(2*b^2) + (6*x*cos(b*x)*SinIntegral(b*x))/b^3 - (x^3*cos(b*x)*SinIntegral(b*x))/b - (6*sin(b*x)*SinIntegral(b*x))/b^4 + (3*x^2*sin(b*x)*SinIntegral(b*x))/b^2]

    @test_int [cos(b*x)*SinIntegral(b*x)/x^3, x, 12, -((b*cos(2*b*x))/(4*x)) - (1/2)*b^2*CannotIntegrate((cos(b*x)*SinIntegral(b*x))/x, x) + (b*sin(b*x)^2)/(2*x) - sin(2*b*x)/(8*x^2) - (cos(b*x)*SinIntegral(b*x))/(2*x^2) + (b*sin(b*x)*SinIntegral(b*x))/(2*x) - b^2*SinIntegral(2*b*x)]
    @test_int [cos(b*x)*SinIntegral(b*x)/x^2, x, 7, b*CosIntegral(2*b*x) - sin(2*b*x)/(2*x) - (cos(b*x)*SinIntegral(b*x))/x - (1/2)*b*SinIntegral(b*x)^2]
    @test_int [cos(b*x)*SinIntegral(b*x)/x, x, 0, CannotIntegrate((cos(b*x)*SinIntegral(b*x))/x, x)]
    @test_int [cos(b*x)*SinIntegral(b*x), x, 5, CosIntegral(2*b*x)/(2*b) - log(x)/(2*b) + (sin(b*x)*SinIntegral(b*x))/b]
    @test_int [x*cos(b*x)*SinIntegral(b*x), x, 9, -(x/(2*b)) + (cos(b*x)*sin(b*x))/(2*b^2) + (cos(b*x)*SinIntegral(b*x))/b^2 + (x*sin(b*x)*SinIntegral(b*x))/b - SinIntegral(2*b*x)/(2*b^2)]
    @test_int [x^2*cos(b*x)*SinIntegral(b*x), x, 13, -(x^2/(4*b)) - CosIntegral(2*b*x)/b^3 + log(x)/b^3 + (x*cos(b*x)*sin(b*x))/(2*b^2) - (5*sin(b*x)^2)/(4*b^3) + (2*x*cos(b*x)*SinIntegral(b*x))/b^2 - (2*sin(b*x)*SinIntegral(b*x))/b^3 + (x^2*sin(b*x)*SinIntegral(b*x))/b]
    @test_int [x^3*cos(b*x)*SinIntegral(b*x), x, 20, (4*x)/b^3 - x^3/(6*b) - (4*cos(b*x)*sin(b*x))/b^4 + (x^2*cos(b*x)*sin(b*x))/(2*b^2) - (2*x*sin(b*x)^2)/b^3 - (6*cos(b*x)*SinIntegral(b*x))/b^4 + (3*x^2*cos(b*x)*SinIntegral(b*x))/b^2 - (6*x*sin(b*x)*SinIntegral(b*x))/b^3 + (x^3*sin(b*x)*SinIntegral(b*x))/b + (3*SinIntegral(2*b*x))/b^4]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(b*x)*SinIntegral(d*x)^n=#

    @test_int [sin(5*x)*SinIntegral(2*x), x, 6, (-(1/5))*cos(5*x)*SinIntegral(2*x) - (1/10)*SinIntegral(3*x) + (1/10)*SinIntegral(7*x)]

    @test_int [cos(5*x)*SinIntegral(2*x), x, 6, (-(1/10))*CosIntegral(3*x) + (1/10)*CosIntegral(7*x) + (1/5)*sin(5*x)*SinIntegral(2*x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(a+b*x)*SinIntegral(a+b*x)^n=#

    #= [x^3*sin(a + b*x)*SinIntegral(a + b*x), x, 32, (7*a*x)/(4*b^3) - x^2/b^2 + (3*cos(a + b*x)^2)/b^4 - (a^2*cos(a + b*x)^2)/(2*b^4) - (3*CosIntegral(2*a + 2*b*x))/b^4 + (3*a^2*CosIntegral(2*a + 2*b*x))/(2*b^4) + (3*log(a + b*x))/b^4 - (3*a^2*log(a + b*x))/(2*b^4) - (7*a*cos(a + b*x)*sin(a + b*x))/(4*b^4) + (2*x*cos(a + b*x)*sin(a + b*x))/b^3 - sin(a + b*x)^2/b^4 - (a*x*sin(a + b*x)^2)/(2*b^3) + (x^2*sin(a + b*x)^2)/(2*b^2) + (6*x*cos(a + b*x)*SinIntegral(a + b*x))/b^3 - (x^3*cos(a + b*x)*SinIntegral(a + b*x))/b - (6*sin(a + b*x)*SinIntegral(a + b*x))/b^4 + (3*x^2*sin(a + b*x)*SinIntegral(a + b*x))/b^2 + (3*a*SinIntegral(2*a + 2*b*x))/b^4 - (a^3*SinIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*sin(a + b*x)*SinIntegral(a + b*x), x, 21, -(x/b^2) + (a*cos(2*a + 2*b*x))/(4*b^3) - (x*cos(2*a + 2*b*x))/(4*b^2) - (a*CosIntegral(2*a + 2*b*x))/b^3 + (a*log(a + b*x))/b^3 + (cos(a + b*x)*sin(a + b*x))/b^3 + sin(2*a + 2*b*x)/(8*b^3) + (2*cos(a + b*x)*SinIntegral(a + b*x))/b^3 - (x^2*cos(a + b*x)*SinIntegral(a + b*x))/b + (2*x*sin(a + b*x)*SinIntegral(a + b*x))/b^2 - SinIntegral(2*a + 2*b*x)/b^3 + (a^2*SinIntegral(2*a + 2*b*x))/(2*b^3)]
    @test_int [x^1*sin(a + b*x)*SinIntegral(a + b*x), x, 11, -(cos(2*a + 2*b*x)/(4*b^2)) + CosIntegral(2*a + 2*b*x)/(2*b^2) - log(a + b*x)/(2*b^2) - (x*cos(a + b*x)*SinIntegral(a + b*x))/b + (sin(a + b*x)*SinIntegral(a + b*x))/b^2 - (a*SinIntegral(2*a + 2*b*x))/(2*b^2)]
    @test_int [x^0*sin(a + b*x)*SinIntegral(a + b*x), x, 4, -((cos(a + b*x)*SinIntegral(a + b*x))/b) + SinIntegral(2*a + 2*b*x)/(2*b)]
    @test_int [sin(a + b*x)*SinIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((sin(a + b*x)*SinIntegral(a + b*x))/x, x)]

    #= [x^3*cos(a + b*x)*SinIntegral(a + b*x), x, 32, (4*x)/b^3 - (a^2*x)/(2*b^3) + (a*x^2)/(4*b^2) - x^3/(6*b) - (3*a*cos(a + b*x)^2)/(2*b^4) + (3*a*CosIntegral(2*a + 2*b*x))/b^4 - (a^3*CosIntegral(2*a + 2*b*x))/(2*b^4) - (3*a*log(a + b*x))/b^4 + (a^3*log(a + b*x))/(2*b^4) - (4*cos(a + b*x)*sin(a + b*x))/b^4 + (a^2*cos(a + b*x)*sin(a + b*x))/(2*b^4) - (a*x*cos(a + b*x)*sin(a + b*x))/(2*b^3) + (x^2*cos(a + b*x)*sin(a + b*x))/(2*b^2) + (a*sin(a + b*x)^2)/(4*b^4) - (2*x*sin(a + b*x)^2)/b^3 - (6*cos(a + b*x)*SinIntegral(a + b*x))/b^4 + (3*x^2*cos(a + b*x)*SinIntegral(a + b*x))/b^2 - (6*x*sin(a + b*x)*SinIntegral(a + b*x))/b^3 + (x^3*sin(a + b*x)*SinIntegral(a + b*x))/b + (3*SinIntegral(2*a + 2*b*x))/b^4 - (3*a^2*SinIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*cos(a + b*x)*SinIntegral(a + b*x), x, 21, (a*x)/(2*b^2) - x^2/(4*b) + cos(2*a + 2*b*x)/(2*b^3) - CosIntegral(2*a + 2*b*x)/b^3 + (a^2*CosIntegral(2*a + 2*b*x))/(2*b^3) + log(a + b*x)/b^3 - (a^2*log(a + b*x))/(2*b^3) - (a*cos(a + b*x)*sin(a + b*x))/(2*b^3) + (x*cos(a + b*x)*sin(a + b*x))/(2*b^2) - sin(a + b*x)^2/(4*b^3) + (2*x*cos(a + b*x)*SinIntegral(a + b*x))/b^2 - (2*sin(a + b*x)*SinIntegral(a + b*x))/b^3 + (x^2*sin(a + b*x)*SinIntegral(a + b*x))/b + (a*SinIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*cos(a + b*x)*SinIntegral(a + b*x), x, 12, -(x/(2*b)) - (a*CosIntegral(2*a + 2*b*x))/(2*b^2) + (a*log(a + b*x))/(2*b^2) + (cos(a + b*x)*sin(a + b*x))/(2*b^2) + (cos(a + b*x)*SinIntegral(a + b*x))/b^2 + (x*sin(a + b*x)*SinIntegral(a + b*x))/b - SinIntegral(2*a + 2*b*x)/(2*b^2)]
    @test_int [x^0*cos(a + b*x)*SinIntegral(a + b*x), x, 4, CosIntegral(2*a + 2*b*x)/(2*b) - log(a + b*x)/(2*b) + (sin(a + b*x)*SinIntegral(a + b*x))/b]
    @test_int [cos(a + b*x)*SinIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((cos(a + b*x)*SinIntegral(a + b*x))/x, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(a+b*x)*SinIntegral(c+d*x)^n=#

    #= [x^2*sin(a + b*x)*SinIntegral(c + d*x), x, 46, -((c*cos(a - c + (b - d)*x))/(2*b*(b - d)*d)) + (x*cos(a - c + (b - d)*x))/(2*b*(b - d)) + (c*cos(a + c + (b + d)*x))/(2*b*d*(b + d)) - (x*cos(a + c + (b + d)*x))/(2*b*(b + d)) + (CosIntegral(((b - d)*(c + d*x))/d)*(2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d)))/(2*b^3*d^2) - (CosIntegral(((b + d)*(c + d*x))/d)*(2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d)))/(2*b^3*d^2) - sin(a - c + (b - d)*x)/(2*b*(b - d)^2) - sin(a - c + (b - d)*x)/(b^2*(b - d)) + sin(a + c + (b + d)*x)/(2*b*(b + d)^2) + sin(a + c + (b + d)*x)/(b^2*(b + d)) + (((2 - b^2*x^2)*cos(a + b*x) + 2*b*x*sin(a + b*x))*SinIntegral(c + d*x))/b^3 - (((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d))*SinIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) + (((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d))*SinIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*sin(a + b*x)*SinIntegral(c + d*x), x, 24, cos(a - c + (b - d)*x)/(2*b*(b - d)) - cos(a + c + (b + d)*x)/(2*b*(b + d)) - (cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b*d) - (c*CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b*d) + (c*cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) - (x*cos(a + b*x)*SinIntegral(c + d*x))/b + (sin(a + b*x)*SinIntegral(c + d*x))/b^2 - (c*cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) - (sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2)]
    @test_int [x^0*sin(a + b*x)*SinIntegral(c + d*x), x, 9, -((CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b)) + (CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b) - (cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (cos(a + b*x)*SinIntegral(c + d*x))/b + (cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [sin(a + b*x)*SinIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((sin(a + b*x)*SinIntegral(c + d*x))/x, x)]

    #= [x^2*cos(a + b*x)*SinIntegral(c + d*x), x, 46, -(cos(a - c + (b - d)*x)/(2*b*(b - d)^2)) - cos(a - c + (b - d)*x)/(b^2*(b - d)) + cos(a + c + (b + d)*x)/(2*b*(b + d)^2) + cos(a + c + (b + d)*x)/(b^2*(b + d)) - (CosIntegral(((b - d)*(c + d*x))/d)*((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d)))/(2*b^3*d^2) + (CosIntegral(((b + d)*(c + d*x))/d)*((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d)))/(2*b^3*d^2) + (c*sin(a - c + (b - d)*x))/(2*b*(b - d)*d) - (x*sin(a - c + (b - d)*x))/(2*b*(b - d)) - (c*sin(a + c + (b + d)*x))/(2*b*d*(b + d)) + (x*sin(a + c + (b + d)*x))/(2*b*(b + d)) + ((2*b*x*cos(a + b*x) - (2 - b^2*x^2)*sin(a + b*x))*SinIntegral(c + d*x))/b^3 - ((2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d))*SinIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) + ((2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d))*SinIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*cos(a + b*x)*SinIntegral(c + d*x), x, 24, (c*cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) - (c*cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) + (CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b^2) - (CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b^2) - sin(a - c + (b - d)*x)/(2*b*(b - d)) + sin(a + c + (b + d)*x)/(2*b*(b + d)) + (cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) - (c*sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (cos(a + b*x)*SinIntegral(c + d*x))/b^2 + (x*sin(a + b*x)*SinIntegral(c + d*x))/b - (cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d)]
    @test_int [x^0*cos(a + b*x)*SinIntegral(c + d*x), x, 9, -((cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b)) + (cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b) + (sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b) + (sin(a + b*x)*SinIntegral(c + d*x))/b - (sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [cos(a + b*x)*SinIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((cos(a + b*x)*SinIntegral(c + d*x))/x, x)]

    #= ::Section::Closed:: =#
    #=Cosine*integral*function=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*CosIntegral(b*x)^n=#

    @test_int [x^m*CosIntegral(b*x), x, 5, (x^(1 + m)*CosIntegral(b*x))/(1 + m) + (im*x^m*Gamma(1 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b*(1 + m))) - (im*x^m*Gamma(1 + m, im*b*x))/((im*b*x)^m*(2*b*(1 + m)))]

    @test_int [x^3*CosIntegral(b*x), x, 6, (3*cos(b*x))/(2*b^4) - (3*x^2*cos(b*x))/(4*b^2) + (1/4)*x^4*CosIntegral(b*x) + (3*x*sin(b*x))/(2*b^3) - (x^3*sin(b*x))/(4*b)]
    @test_int [x^2*CosIntegral(b*x), x, 5, -((2*x*cos(b*x))/(3*b^2)) + (1/3)*x^3*CosIntegral(b*x) + (2*sin(b*x))/(3*b^3) - (x^2*sin(b*x))/(3*b)]
    @test_int [x^1*CosIntegral(b*x), x, 4, -(cos(b*x)/(2*b^2)) + (1/2)*x^2*CosIntegral(b*x) - (x*sin(b*x))/(2*b)]
    @test_int [x^0*CosIntegral(b*x), x, 1, x*CosIntegral(b*x) - sin(b*x)/b]
    @test_int [CosIntegral(b*x)/x^1, x, 1, (-(1/2))*im*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-im)*b*x) + (1/2)*im*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], im*b*x) + EulerGamma*log(x) + (1/2)*log(b*x)^2]
    @test_int [CosIntegral(b*x)/x^2, x, 4, -(cos(b*x)/x) - CosIntegral(b*x)/x - b*SinIntegral(b*x)]
    @test_int [CosIntegral(b*x)/x^3, x, 5, -(cos(b*x)/(4*x^2)) - (1/4)*b^2*CosIntegral(b*x) - CosIntegral(b*x)/(2*x^2) + (b*sin(b*x))/(4*x)]

    @test_int [x^m*CosIntegral(b*x)^2, x, 0, CannotIntegrate(x^m*CosIntegral(b*x)^2, x)]

    @test_int [x^3*CosIntegral(b*x)^2, x, 19, x^2/(4*b^2) + (3*cos(b*x)^2)/(8*b^4) + (3*cos(b*x)*CosIntegral(b*x))/b^4 - (3*x^2*cos(b*x)*CosIntegral(b*x))/(2*b^2) + (1/4)*x^4*CosIntegral(b*x)^2 - (3*CosIntegral(2*b*x))/(2*b^4) - (3*log(x))/(2*b^4) + (x*cos(b*x)*sin(b*x))/b^3 + (3*x*CosIntegral(b*x)*sin(b*x))/b^3 - (x^3*CosIntegral(b*x)*sin(b*x))/(2*b) - (13*sin(b*x)^2)/(8*b^4) + (x^2*sin(b*x)^2)/(4*b^2)]
    @test_int [x^2*CosIntegral(b*x)^2, x, 15, x/(2*b^2) - (4*x*cos(b*x)*CosIntegral(b*x))/(3*b^2) + (1/3)*x^3*CosIntegral(b*x)^2 + (5*cos(b*x)*sin(b*x))/(6*b^3) + (4*CosIntegral(b*x)*sin(b*x))/(3*b^3) - (2*x^2*CosIntegral(b*x)*sin(b*x))/(3*b) + (x*sin(b*x)^2)/(3*b^2) - (2*SinIntegral(2*b*x))/(3*b^3)]
    @test_int [x^1*CosIntegral(b*x)^2, x, 10, -((cos(b*x)*CosIntegral(b*x))/b^2) + (1/2)*x^2*CosIntegral(b*x)^2 + CosIntegral(2*b*x)/(2*b^2) + log(x)/(2*b^2) - (x*CosIntegral(b*x)*sin(b*x))/b + sin(b*x)^2/(2*b^2)]
    @test_int [x^0*CosIntegral(b*x)^2, x, 6, x*CosIntegral(b*x)^2 - (2*CosIntegral(b*x)*sin(b*x))/b + SinIntegral(2*b*x)/b]
    @test_int [CosIntegral(b*x)^2/x^1, x, 0, CannotIntegrate(CosIntegral(b*x)^2/x, x)]
    @test_int [CosIntegral(b*x)^2/x^2, x, 0, CannotIntegrate(CosIntegral(b*x)^2/x^2, x)]
    @test_int [CosIntegral(b*x)^2/x^3, x, 0, CannotIntegrate(CosIntegral(b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*CosIntegral(a+b*x)^n=#

    @test_int [x^m*CosIntegral(a + b*x), x, 1, (x^(1 + m)*CosIntegral(a + b*x))/(1 + m) - (b*CannotIntegrate((x^(1 + m)*cos(a + b*x))/(a + b*x), x))/(1 + m)]

    @test_int [x^3*CosIntegral(a + b*x), x, 14, (3*cos(a + b*x))/(2*b^4) - (a^2*cos(a + b*x))/(4*b^4) + (a*x*cos(a + b*x))/(2*b^3) - (3*x^2*cos(a + b*x))/(4*b^2) - (a^4*CosIntegral(a + b*x))/(4*b^4) + (1/4)*x^4*CosIntegral(a + b*x) - (a*sin(a + b*x))/(2*b^4) + (a^3*sin(a + b*x))/(4*b^4) + (3*x*sin(a + b*x))/(2*b^3) - (a^2*x*sin(a + b*x))/(4*b^3) + (a*x^2*sin(a + b*x))/(4*b^2) - (x^3*sin(a + b*x))/(4*b)]
    @test_int [x^2*CosIntegral(a + b*x), x, 10, (a*cos(a + b*x))/(3*b^3) - (2*x*cos(a + b*x))/(3*b^2) + (a^3*CosIntegral(a + b*x))/(3*b^3) + (1/3)*x^3*CosIntegral(a + b*x) + (2*sin(a + b*x))/(3*b^3) - (a^2*sin(a + b*x))/(3*b^3) + (a*x*sin(a + b*x))/(3*b^2) - (x^2*sin(a + b*x))/(3*b)]
    @test_int [x^1*CosIntegral(a + b*x), x, 7, -(cos(a + b*x)/(2*b^2)) - (a^2*CosIntegral(a + b*x))/(2*b^2) + (1/2)*x^2*CosIntegral(a + b*x) + (a*sin(a + b*x))/(2*b^2) - (x*sin(a + b*x))/(2*b)]
    @test_int [x^0*CosIntegral(a + b*x), x, 1, ((a + b*x)*CosIntegral(a + b*x))/b - sin(a + b*x)/b]
    @test_int [CosIntegral(a + b*x)/x^1, x, 0, CannotIntegrate(CosIntegral(a + b*x)/x, x)]
    @test_int [CosIntegral(a + b*x)/x^2, x, 7, (b*cos(a)*CosIntegral(b*x))/a - (b*CosIntegral(a + b*x))/a - CosIntegral(a + b*x)/x - (b*sin(a)*SinIntegral(b*x))/a]
    @test_int [CosIntegral(a + b*x)/x^3, x, 11, -((b*cos(a + b*x))/(2*a*x)) - (b^2*cos(a)*CosIntegral(b*x))/(2*a^2) + (b^2*CosIntegral(a + b*x))/(2*a^2) - CosIntegral(a + b*x)/(2*x^2) - (b^2*CosIntegral(b*x)*sin(a))/(2*a) - (b^2*cos(a)*SinIntegral(b*x))/(2*a) + (b^2*sin(a)*SinIntegral(b*x))/(2*a^2)]

    @test_int [x^m*CosIntegral(a + b*x)^2, x, 0, CannotIntegrate(x^m*CosIntegral(a + b*x)^2, x)]

    @test_int [x^2*CosIntegral(a + b*x)^2, x, 39, (2*x)/(3*b^2) + (a*cos(2*a + 2*b*x))/(3*b^3) - (x*cos(2*a + 2*b*x))/(6*b^2) + (2*a*cos(a + b*x)*CosIntegral(a + b*x))/(3*b^3) - (4*x*cos(a + b*x)*CosIntegral(a + b*x))/(3*b^2) + (a^2*(a + b*x)*CosIntegral(a + b*x)^2)/(3*b^3) - (a*x*(a + b*x)*CosIntegral(a + b*x)^2)/(3*b^2) + (x^2*(a + b*x)*CosIntegral(a + b*x)^2)/(3*b) - (a*CosIntegral(2*a + 2*b*x))/b^3 - (a*log(a + b*x))/b^3 + (2*cos(a + b*x)*sin(a + b*x))/(3*b^3) + (4*CosIntegral(a + b*x)*sin(a + b*x))/(3*b^3) - (2*a^2*CosIntegral(a + b*x)*sin(a + b*x))/(3*b^3) + (2*a*x*CosIntegral(a + b*x)*sin(a + b*x))/(3*b^2) - (2*x^2*CosIntegral(a + b*x)*sin(a + b*x))/(3*b) + sin(2*a + 2*b*x)/(12*b^3) - (2*SinIntegral(2*a + 2*b*x))/(3*b^3) + (a^2*SinIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*CosIntegral(a + b*x)^2, x, 17, -(cos(2*a + 2*b*x)/(4*b^2)) - (cos(a + b*x)*CosIntegral(a + b*x))/b^2 - (a*(a + b*x)*CosIntegral(a + b*x)^2)/(2*b^2) + (x*(a + b*x)*CosIntegral(a + b*x)^2)/(2*b) + CosIntegral(2*a + 2*b*x)/(2*b^2) + log(a + b*x)/(2*b^2) + (a*CosIntegral(a + b*x)*sin(a + b*x))/b^2 - (x*CosIntegral(a + b*x)*sin(a + b*x))/b - (a*SinIntegral(2*a + 2*b*x))/b^2]
    @test_int [x^0*CosIntegral(a + b*x)^2, x, 5, ((a + b*x)*CosIntegral(a + b*x)^2)/b - (2*CosIntegral(a + b*x)*sin(a + b*x))/b + SinIntegral(2*a + 2*b*x)/b]
    @test_int [CosIntegral(a + b*x)^2/x^1, x, 0, CannotIntegrate(CosIntegral(a + b*x)^2/x, x)]
    @test_int [CosIntegral(a + b*x)^2/x^2, x, 0, CannotIntegrate(CosIntegral(a + b*x)^2/x^2, x)]
    @test_int [CosIntegral(a + b*x)^2/x^3, x, 0, CannotIntegrate(CosIntegral(a + b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*CosIntegral(d*(a+b*log(c*x^n)))=#

    @test_int [x^2*CosIntegral(d*(a + b*log(c*x^n))), x, 7, (x^3*CosIntegral(d*(a + b*log(c*x^n))))/3 - (x^3*ExpIntegralEi(((3 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) - (x^3*ExpIntegralEi(((3 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n))]
    @test_int [x^1*CosIntegral(d*(a + b*log(c*x^n))), x, 7, (x^2*CosIntegral(d*(a + b*log(c*x^n))))/2 - (x^2*ExpIntegralEi(((2 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) - (x^2*ExpIntegralEi(((2 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n))]
    @test_int [x^0*CosIntegral(d*(a + b*log(c*x^n))), x, 7, x*CosIntegral(d*(a + b*log(c*x^n))) - (x*ExpIntegralEi(((1 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1)) - (x*ExpIntegralEi(((1 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1))]
    @test_int [CosIntegral(d*(a + b*log(c*x^n)))/x^1, x, 3, (CosIntegral(d*(a + b*log(c*x^n)))*(a + b*log(c*x^n)))/(b*n) - sin(d*(a + b*log(c*x^n)))/(b*d*n)]
    @test_int [CosIntegral(d*(a + b*log(c*x^n)))/x^2, x, 7, -(CosIntegral(d*(a + b*log(c*x^n)))/x) + (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x) + (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x)]
    @test_int [CosIntegral(d*(a + b*log(c*x^n)))/x^3, x, 7, -CosIntegral(d*(a + b*log(c*x^n)))/(2*x^2) + (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 - im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2) + (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 + im*b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2)]

    @test_int [(e*x)^m*CosIntegral(d*(a + b*log(c*x^n))), x, 7, ((e*x)^(1 + m)*CosIntegral(d*(a + b*log(c*x^n))))/(e*(1 + m)) - (x*(e*x)^m*ExpIntegralEi(((1 + m - im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n)) - (x*(e*x)^m*ExpIntegralEi(((1 + m + im*b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(b*x)*CosIntegral(b*x)^n=#

    @test_int [sin(b*x)*CosIntegral(b*x)/x^3, x, 12, -((b*cos(b*x)^2)/(2*x)) - (b*cos(2*b*x))/(4*x) - (b*cos(b*x)*CosIntegral(b*x))/(2*x) - (1/2)*b^2*CannotIntegrate((CosIntegral(b*x)*sin(b*x))/x, x) - (CosIntegral(b*x)*sin(b*x))/(2*x^2) - sin(2*b*x)/(8*x^2) - b^2*SinIntegral(2*b*x)]
    @test_int [sin(b*x)*CosIntegral(b*x)/x^2, x, 7, (1/2)*b*CosIntegral(b*x)^2 + b*CosIntegral(2*b*x) - (CosIntegral(b*x)*sin(b*x))/x - sin(2*b*x)/(2*x)]
    @test_int [sin(b*x)*CosIntegral(b*x)/x, x, 0, CannotIntegrate((CosIntegral(b*x)*sin(b*x))/x, x)]
    @test_int [sin(b*x)*CosIntegral(b*x), x, 5, -((cos(b*x)*CosIntegral(b*x))/b) + CosIntegral(2*b*x)/(2*b) + log(x)/(2*b)]
    @test_int [x*sin(b*x)*CosIntegral(b*x), x, 9, x/(2*b) - (x*cos(b*x)*CosIntegral(b*x))/b + (cos(b*x)*sin(b*x))/(2*b^2) + (CosIntegral(b*x)*sin(b*x))/b^2 - SinIntegral(2*b*x)/(2*b^2)]
    @test_int [x^2*sin(b*x)*CosIntegral(b*x), x, 13, x^2/(4*b) + cos(b*x)^2/(4*b^3) + (2*cos(b*x)*CosIntegral(b*x))/b^3 - (x^2*cos(b*x)*CosIntegral(b*x))/b - CosIntegral(2*b*x)/b^3 - log(x)/b^3 + (x*cos(b*x)*sin(b*x))/(2*b^2) + (2*x*CosIntegral(b*x)*sin(b*x))/b^2 - sin(b*x)^2/b^3]
    @test_int [x^3*sin(b*x)*CosIntegral(b*x), x, 20, -((5*x)/(2*b^3)) + x^3/(6*b) + (x*cos(b*x)^2)/(2*b^3) + (6*x*cos(b*x)*CosIntegral(b*x))/b^3 - (x^3*cos(b*x)*CosIntegral(b*x))/b - (4*cos(b*x)*sin(b*x))/b^4 + (x^2*cos(b*x)*sin(b*x))/(2*b^2) - (6*CosIntegral(b*x)*sin(b*x))/b^4 + (3*x^2*CosIntegral(b*x)*sin(b*x))/b^2 - (3*x*sin(b*x)^2)/(2*b^3) + (3*SinIntegral(2*b*x))/b^4]

    @test_int [cos(b*x)*CosIntegral(b*x)/x^3, x, 14, -(cos(b*x)^2/(4*x^2)) - (cos(b*x)*CosIntegral(b*x))/(2*x^2) - (1/4)*b^2*CosIntegral(b*x)^2 - b^2*CosIntegral(2*b*x) + (b*cos(b*x)*sin(b*x))/(2*x) + (b*CosIntegral(b*x)*sin(b*x))/(2*x) + (b*sin(2*b*x))/(4*x)]
    @test_int [cos(b*x)*CosIntegral(b*x)/x^2, x, 5, -(cos(b*x)^2/x) - (cos(b*x)*CosIntegral(b*x))/x - b*CannotIntegrate((CosIntegral(b*x)*sin(b*x))/x, x) - b*SinIntegral(2*b*x)]
    @test_int [cos(b*x)*CosIntegral(b*x)/x, x, 1, (1/2)*CosIntegral(b*x)^2]
    @test_int [cos(b*x)*CosIntegral(b*x), x, 5, (CosIntegral(b*x)*sin(b*x))/b - SinIntegral(2*b*x)/(2*b)]
    @test_int [x*cos(b*x)*CosIntegral(b*x), x, 9, (cos(b*x)*CosIntegral(b*x))/b^2 - CosIntegral(2*b*x)/(2*b^2) - log(x)/(2*b^2) + (x*CosIntegral(b*x)*sin(b*x))/b - sin(b*x)^2/(2*b^2)]
    @test_int [x^2*cos(b*x)*CosIntegral(b*x), x, 14, -((3*x)/(4*b^2)) + (2*x*cos(b*x)*CosIntegral(b*x))/b^2 - (5*cos(b*x)*sin(b*x))/(4*b^3) - (2*CosIntegral(b*x)*sin(b*x))/b^3 + (x^2*CosIntegral(b*x)*sin(b*x))/b - (x*sin(b*x)^2)/(2*b^2) + SinIntegral(2*b*x)/b^3]
    @test_int [x^3*cos(b*x)*CosIntegral(b*x), x, 18, -(x^2/(2*b^2)) - (3*cos(b*x)^2)/(4*b^4) - (6*cos(b*x)*CosIntegral(b*x))/b^4 + (3*x^2*cos(b*x)*CosIntegral(b*x))/b^2 + (3*CosIntegral(2*b*x))/b^4 + (3*log(x))/b^4 - (2*x*cos(b*x)*sin(b*x))/b^3 - (6*x*CosIntegral(b*x)*sin(b*x))/b^3 + (x^3*CosIntegral(b*x)*sin(b*x))/b + (13*sin(b*x)^2)/(4*b^4) - (x^2*sin(b*x)^2)/(2*b^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(b*x)*CosIntegral(d*x)^n=#

    @test_int [sin(5*x)*CosIntegral(2*x), x, 6, (-(1/5))*cos(5*x)*CosIntegral(2*x) + (1/10)*CosIntegral(3*x) + (1/10)*CosIntegral(7*x)]

    @test_int [cos(5*x)*CosIntegral(2*x), x, 6, (1/5)*CosIntegral(2*x)*sin(5*x) - (1/10)*SinIntegral(3*x) - (1/10)*SinIntegral(7*x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(a+b*x)*CosIntegral(a+b*x)^n=#

    #= [x^3*sin(a + b*x)*CosIntegral(a + b*x), x, 32, -((5*x)/(2*b^3)) + (a^2*x)/(2*b^3) - (a*x^2)/(4*b^2) + x^3/(6*b) - (7*a*cos(a + b*x)^2)/(4*b^4) + (x*cos(a + b*x)^2)/(2*b^3) + (6*x*cos(a + b*x)*CosIntegral(a + b*x))/b^3 - (x^3*cos(a + b*x)*CosIntegral(a + b*x))/b + (3*a*CosIntegral(2*a + 2*b*x))/b^4 - (a^3*CosIntegral(2*a + 2*b*x))/(2*b^4) + (3*a*log(a + b*x))/b^4 - (a^3*log(a + b*x))/(2*b^4) - (4*cos(a + b*x)*sin(a + b*x))/b^4 + (a^2*cos(a + b*x)*sin(a + b*x))/(2*b^4) - (a*x*cos(a + b*x)*sin(a + b*x))/(2*b^3) + (x^2*cos(a + b*x)*sin(a + b*x))/(2*b^2) - (6*CosIntegral(a + b*x)*sin(a + b*x))/b^4 + (3*x^2*CosIntegral(a + b*x)*sin(a + b*x))/b^2 - (3*x*sin(a + b*x)^2)/(2*b^3) + (3*SinIntegral(2*a + 2*b*x))/b^4 - (3*a^2*SinIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*sin(a + b*x)*CosIntegral(a + b*x), x, 21, -((a*x)/(2*b^2)) + x^2/(4*b) + cos(a + b*x)^2/(4*b^3) + cos(2*a + 2*b*x)/(2*b^3) + (2*cos(a + b*x)*CosIntegral(a + b*x))/b^3 - (x^2*cos(a + b*x)*CosIntegral(a + b*x))/b - CosIntegral(2*a + 2*b*x)/b^3 + (a^2*CosIntegral(2*a + 2*b*x))/(2*b^3) - log(a + b*x)/b^3 + (a^2*log(a + b*x))/(2*b^3) - (a*cos(a + b*x)*sin(a + b*x))/(2*b^3) + (x*cos(a + b*x)*sin(a + b*x))/(2*b^2) + (2*x*CosIntegral(a + b*x)*sin(a + b*x))/b^2 + (a*SinIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*sin(a + b*x)*CosIntegral(a + b*x), x, 12, x/(2*b) - (x*cos(a + b*x)*CosIntegral(a + b*x))/b - (a*CosIntegral(2*a + 2*b*x))/(2*b^2) - (a*log(a + b*x))/(2*b^2) + (cos(a + b*x)*sin(a + b*x))/(2*b^2) + (CosIntegral(a + b*x)*sin(a + b*x))/b^2 - SinIntegral(2*a + 2*b*x)/(2*b^2)]
    @test_int [x^0*sin(a + b*x)*CosIntegral(a + b*x), x, 4, -((cos(a + b*x)*CosIntegral(a + b*x))/b) + CosIntegral(2*a + 2*b*x)/(2*b) + log(a + b*x)/(2*b)]
    @test_int [sin(a + b*x)*CosIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((sin(a + b*x)*CosIntegral(a + b*x))/x, x)]

    #= [x^3*cos(a + b*x)*CosIntegral(a + b*x), x, 32, (5*a*x)/(4*b^3) - x^2/(2*b^2) - (15*cos(a + b*x)^2)/(4*b^4) + (a^2*cos(a + b*x)^2)/(2*b^4) - (6*cos(a + b*x)*CosIntegral(a + b*x))/b^4 + (3*x^2*cos(a + b*x)*CosIntegral(a + b*x))/b^2 + (3*CosIntegral(2*a + 2*b*x))/b^4 - (3*a^2*CosIntegral(2*a + 2*b*x))/(2*b^4) + (3*log(a + b*x))/b^4 - (3*a^2*log(a + b*x))/(2*b^4) + (7*a*cos(a + b*x)*sin(a + b*x))/(4*b^4) - (2*x*cos(a + b*x)*sin(a + b*x))/b^3 - (6*x*CosIntegral(a + b*x)*sin(a + b*x))/b^3 + (x^3*CosIntegral(a + b*x)*sin(a + b*x))/b + sin(a + b*x)^2/(4*b^4) + (a*x*sin(a + b*x)^2)/(2*b^3) - (x^2*sin(a + b*x)^2)/(2*b^2) - (3*a*SinIntegral(2*a + 2*b*x))/b^4 + (a^3*SinIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*cos(a + b*x)*CosIntegral(a + b*x), x, 21, -(x/b^2) - (a*cos(2*a + 2*b*x))/(4*b^3) + (x*cos(2*a + 2*b*x))/(4*b^2) + (2*x*cos(a + b*x)*CosIntegral(a + b*x))/b^2 + (a*CosIntegral(2*a + 2*b*x))/b^3 + (a*log(a + b*x))/b^3 - (cos(a + b*x)*sin(a + b*x))/b^3 - (2*CosIntegral(a + b*x)*sin(a + b*x))/b^3 + (x^2*CosIntegral(a + b*x)*sin(a + b*x))/b - sin(2*a + 2*b*x)/(8*b^3) + SinIntegral(2*a + 2*b*x)/b^3 - (a^2*SinIntegral(2*a + 2*b*x))/(2*b^3)]
    @test_int [x^1*cos(a + b*x)*CosIntegral(a + b*x), x, 11, cos(2*a + 2*b*x)/(4*b^2) + (cos(a + b*x)*CosIntegral(a + b*x))/b^2 - CosIntegral(2*a + 2*b*x)/(2*b^2) - log(a + b*x)/(2*b^2) + (x*CosIntegral(a + b*x)*sin(a + b*x))/b + (a*SinIntegral(2*a + 2*b*x))/(2*b^2)]
    @test_int [x^0*cos(a + b*x)*CosIntegral(a + b*x), x, 4, (CosIntegral(a + b*x)*sin(a + b*x))/b - SinIntegral(2*a + 2*b*x)/(2*b)]
    @test_int [cos(a + b*x)*CosIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((cos(a + b*x)*CosIntegral(a + b*x))/x, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Trig(a+b*x)*CosIntegral(c+d*x)^n=#

    #= [x^2*sin(a + b*x)*CosIntegral(c + d*x), x, 46, cos(a - c + (b - d)*x)/(2*b*(b - d)^2) + cos(a - c + (b - d)*x)/(b^2*(b - d)) + cos(a + c + (b + d)*x)/(2*b*(b + d)^2) + cos(a + c + (b + d)*x)/(b^2*(b + d)) + (CosIntegral(((b - d)*(c + d*x))/d)*((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d)))/(2*b^3*d^2) + (CosIntegral(((b + d)*(c + d*x))/d)*((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d)))/(2*b^3*d^2) + (CosIntegral(c + d*x)*((2 - b^2*x^2)*cos(a + b*x) + 2*b*x*sin(a + b*x)))/b^3 - (c*sin(a - c + (b - d)*x))/(2*b*(b - d)*d) + (x*sin(a - c + (b - d)*x))/(2*b*(b - d)) - (c*sin(a + c + (b + d)*x))/(2*b*d*(b + d)) + (x*sin(a + c + (b + d)*x))/(2*b*(b + d)) + ((2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d))*SinIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) + ((2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d))*SinIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*sin(a + b*x)*CosIntegral(c + d*x), x, 24, -((c*cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d)) - (x*cos(a + b*x)*CosIntegral(c + d*x))/b - (c*cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) - (CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b^2) - (CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b^2) + (CosIntegral(c + d*x)*sin(a + b*x))/b^2 + sin(a - c + (b - d)*x)/(2*b*(b - d)) + sin(a + c + (b + d)*x)/(2*b*(b + d)) - (cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (c*sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) - (cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d)]
    @test_int [x^0*sin(a + b*x)*CosIntegral(c + d*x), x, 9, (cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (cos(a + b*x)*CosIntegral(c + d*x))/b + (cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b) - (sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [sin(a + b*x)*CosIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((CosIntegral(c + d*x)*sin(a + b*x))/x, x)]

    #= [x^2*cos(a + b*x)*CosIntegral(c + d*x), x, 46, -((c*cos(a - c + (b - d)*x))/(2*b*(b - d)*d)) + (x*cos(a - c + (b - d)*x))/(2*b*(b - d)) - (c*cos(a + c + (b + d)*x))/(2*b*d*(b + d)) + (x*cos(a + c + (b + d)*x))/(2*b*(b + d)) + (CosIntegral(((b - d)*(c + d*x))/d)*(2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d)))/(2*b^3*d^2) + (CosIntegral(((b + d)*(c + d*x))/d)*(2*b*c*d*cos(a - (b*c)/d) - (b^2*c^2 - 2*d^2)*sin(a - (b*c)/d)))/(2*b^3*d^2) + (CosIntegral(c + d*x)*(2*b*x*cos(a + b*x) - (2 - b^2*x^2)*sin(a + b*x)))/b^3 - sin(a - c + (b - d)*x)/(2*b*(b - d)^2) - sin(a - c + (b - d)*x)/(b^2*(b - d)) - sin(a + c + (b + d)*x)/(2*b*(b + d)^2) - sin(a + c + (b + d)*x)/(b^2*(b + d)) - (((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d))*SinIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) - (((b^2*c^2 - 2*d^2)*cos(a - (b*c)/d) + 2*b*c*d*sin(a - (b*c)/d))*SinIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*cos(a + b*x)*CosIntegral(c + d*x), x, 24, cos(a - c + (b - d)*x)/(2*b*(b - d)) + cos(a + c + (b + d)*x)/(2*b*(b + d)) - (cos(a - (b*c)/d)*CosIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (cos(a + b*x)*CosIntegral(c + d*x))/b^2 - (cos(a - (b*c)/d)*CosIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b*d) + (c*CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b*d) + (x*CosIntegral(c + d*x)*sin(a + b*x))/b + (c*cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (sin(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (c*cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) + (sin(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2)]
    @test_int [x^0*cos(a + b*x)*CosIntegral(c + d*x), x, 9, -((CosIntegral((c*(b - d))/d + (b - d)*x)*sin(a - (b*c)/d))/(2*b)) - (CosIntegral((c*(b + d))/d + (b + d)*x)*sin(a - (b*c)/d))/(2*b) + (CosIntegral(c + d*x)*sin(a + b*x))/b - (cos(a - (b*c)/d)*SinIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (cos(a - (b*c)/d)*SinIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [cos(a + b*x)*CosIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((cos(a + b*x)*CosIntegral(c + d*x))/x, x)]
    end

max_stars_repo_path: test/8 Special functions/8.4 Trig integral functions.jl


text:
module MyModelMenagerie

using ModelingToolkit, OrdinaryDiffEq

include("models/CircadianOscillator.jl")
include("models/MassSpringOscillator.jl")
include("models/MichaelisMenten.jl")
include("models/STG_Liu.jl")
include("models/NFKB.jl")

export CircadianOscillator, CalciumNeuron, MassSpringOscillator, NFKBModel, NFKB_output_map

end


text:
module Gui

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf
using JSON

const DISPLAY_RATIO = 16.0 / 9.0
const GUI_WIDTH = 1024 * 2 + 512 * 2
const GUI_HEIGHT = UInt32(Float64(GUI_WIDTH) / DISPLAY_RATIO)

using ..Model

include("init.jl")
include("gui_data.jl")
include("graphs/graphs.jl")

using .Graphs
using ..Comm

include("main_panel.jl")
include("popup_window.jl")

function run(gui_data::GuiData, app_data::Model.AppData, sock::Comm.SocClient)
    while !GLFW.WindowShouldClose(gui_data.window)

        GLFW.PollEvents()

        begin_render()

        draw_main_panel(gui_data, app_data, sock)

        # Graph rendering -----------------------------------------------------------
        vert_pos = 25

        # Renders boths stimulus and cell output spikes.
        Graphs.draw(Graphs.spikes_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        # vert_pos += Graphs.GRAPH_WINDOW_HEIGHT + 20
        # Graphs.draw(Graphs.synapse_input_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT + 20
        Graphs.draw(Graphs.synapse_surge_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT
        Graphs.draw(Graphs.synapse_psp_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT
        Graphs.draw(Graphs.synapse_weights_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT
        Graphs.draw(Graphs.soma_psp_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT
        Graphs.draw(Graphs.soma_apFast_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

        vert_pos += Graphs.GRAPH_WINDOW_HEIGHT
        Graphs.draw(Graphs.soma_apSlow_graph, gui_data::GuiData, app_data.model, app_data.samples, vert_pos)

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

        end_render(gui_data)

        # App communications
        Comm.read_channel(sock)
    end
end

function begin_render()
    # start the Dear ImGui frame
    ImGui_ImplOpenGL3_NewFrame()
    ImGui_ImplGlfw_NewFrame()
    CImGui.NewFrame()
end

function end_render(data::GuiData)
        # rendering
    CImGui.Render()
    GLFW.MakeContextCurrent(data.window)

    glViewport(0, 0, data.display_w, data.display_h)
    glClearColor(data.clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)

    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

    GLFW.MakeContextCurrent(data.window)
    GLFW.SwapBuffers(data.window)
end

function shutdown(data::GuiData)
    println("Shutting down...")

    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(data.ctx)

    GLFW.DestroyWindow(data.window)
end

end # Module ---------------------------------------------


text:
export phasespace, plot_linear_regions
export plot_dataset, plot_trajectory
using PyPlot

"""
    phasespace(ds, limits, density::Int, t::Int; kwargs...)
    phasespace(ds, ics, t::Int; kwargs...)
Plot a phasespace of a 2D discrete dynamical system, by evolving it for `t` time.

Initial conditions can be given in the form of `limits, density` with
`limits = [(x0, xf), (y0, yf)]` or directly
as a vector of initial conditions using the second method.

## Keyword Arguments:
* `ax=PyPlot.gca()` : the axes to plot on.
* `kwargs...` : Keyword arguments passed into `ax[:plot]()`.
"""
function phasespace(ds::DiscreteDynamicalSystem{IIP, S, 2}, limits,
    density, t;
    ax = PyPlot.gca(), kwargs...) where {IIP, S}

    data = Dataset{2, eltype(S)}()
    for x in range(limits[1][1], stop=limits[1][2], length=density)
        for y in range(limits[2][1], stop=limits[2][2], length=density)
            append!(data, trajectory(ds, t, SVector{2}(x,y)))
        end
    end
    ax[:plot](data[:,1], data[:,2];
    marker = "s", ms = 0.5, color="black", kwargs..., lw=0)
    return nothing
end
#
function phasespace(ds::DiscreteDynamicalSystem{IIP, S, 2},
    ics::Vector{<:AbstractVector}, t;
    ax = PyPlot.gca(), kwargs...) where {IIP, S}

    data = Dataset{2, T}()
    for ic in ics
        append!(data, trajectory(ds, t, SVector{2}(ic)))
    end
    ax[:plot](data[:,1], data[:,2];
    marker = "s", ms = 0.5, color="black", kwargs..., lw=0)
    return nothing
end

function _plot_lrs(x, y, lrs, tangents)
  for i ∈ 1:length(lrs)-1
    PyPlot.plot(x[lrs[i]:lrs[i+1]], y[lrs[i]:lrs[i+1]])
  end
end

"""
    plot_linear_regions(x, y; dxi = 1, tol = 0.2)
Visualize the outcome of calling `linear_regions` by plotting each
linear segment with a different color (on the current axes).
"""
function plot_linear_regions(x, y; dxi = 1, tol = 0.2)
    PyPlot.plot(x, y; ms= 5, marker="o", color = "black", lw=0)
    _plot_lrs(x, y, linear_regions(x, y;  dxi = dxi, tol = tol)...)
end

"""
    plot_dataset(data::Dataset, tvec)
Plot each column of a dataset as a timeseries vs `tvec`.
"""
function plot_dataset(data::Dataset, tvec)
    for (i, col) in enumerate(columns(data))
        PyPlot.plot(tvec, col, label = "var. $i")
    end
    PyPlot.legend()
end

"""
    plot_trajectory(ds::DynamicalSystem, T [, dt])
Plot each column of the `trajectory` as a timeseries vs `tvec`.
"""
function plot_trajectory(ds::DDS, T, dt = 1)
    data = trajectory(ds, T, dt=dt)
    tvec = inittime(ds):dt:(T+inittime(ds))
    plot_dataset(data, tvec)
end
function plot_trajectory(ds::CDS, T, dt = 0.01)
    data = trajectory(ds, T, dt=dt)
    tvec = inittime(ds):dt:(T+inittime(ds))
    plot_dataset(data, tvec)
end


text:
module ROSCar1D
export main

using RobotOS
using PyCall
roslib = pyimport("roslib")
roslib.load_manifest("amrl_msgs")
@rosimport amrl_msgs.msg: AckermannCurvatureDriveMsg
@rosimport nav_msgs.msg: Odometry
rostypegen(@__MODULE__)
#import .amrl_msgs.msg: AckermannCurvatureDriveMsg
#import .nav_msgs.msg: Odometry
using .amrl_msgs.msg: AckermannCurvatureDriveMsg
using .nav_msgs.msg: Odometry
using StaticArrays
using Unrolled
using Random

using OneVision
using OneVision.Car1DExample

global x_cur = CarX(0.0, 0.0) # CarX needs to be immutalble - will find a container to wrap this around.

function odom_callback(msg::Odometry)
    x_cur.pos = sqrt(square(msg.position.position.x) + square(msg.position.position.y))
    x_cur.velocity = 2.0 * atan(msg.pose.pose.orientation.z, msg.pose.pose.orientation.w)
end

function ros_publish_actuation(pub::Publisher, act::CarU, state::CarX)
    # interperate
    next_vel = state.velocity + act.acc*0.05
    msg = AckermannCurvatureDriveMsg()
    msg.velocity = next_vel
    publish(pub, msg)
    publish(pub, msg)
    publish(pub, msg)
end

"The simple wall dynamics model used by OneVision."
struct WallObsModel <: ObsDynamics end

OneVision.obs_forward(dy::WallObsModel, x::CarX, z::CarZ, t::𝕋) = z

function make_ros_controllers(times, freq::ℝ; noise=0.0, plot_result=false, log_prediction=false, N=1)
    delta_t = 1 / freq
    times::Vector{𝕋} = collect(times)
    idx_to_time(xs) = (xs .- 1) .* delta_t
    t0, t_end = times[1], times[end]
    rng = MersenneTwister(1234)

    agent_info(id) = begin
        acc_noise = randn(rng, ℝ, 1 + t_end - t0) * noise
        sys_dy = car_system(delta_t, t -> CarX(0.0, acc_noise[1 + t - t0]))
        obs_dy =
            if id == 1; WallObsDynamics(wall_position=30.0, detector_range=8.0)
            else WallObsModel() end

        sys_dy, obs_dy
    end

    world_dynamics = WorldDynamics([agent_info(i) for i in 1:N])
    init_states = fill((CarX(0.0, 0.0), CarZ(0.0, 0.0), CarU(0.0)), N)
    delays = DelayModel(obs=2, act=3, com=4)  # DelayModel(obs=1, act=2, com=5)
    world_model = WorldDynamics(
        fill((car_system(delta_t), WallObsModel()), N))

    comps = ["pos", "velocity", "acceleration", "obstacle"]
    function record_f(xs, zs, us)
        [(x -> x.pos).(xs) (x -> x.velocity).(xs) (u -> u.acc).(us) (z -> z.distance).(zs)]
    end

    H = 20
    πc = LeaderFollowerControl(warm_up_time=delays.act)

    u_weights = fill(CarU(1.0), N)
    x_weights = fill(CarX(1.0, 1.0), N)

    controller_framework = OvCF{N,CarX{ℝ},CarZ{ℝ},CarU{ℝ},H}(
            πc,
            world_model, delays, x_weights, u_weights,
            FuncT(Tuple{ℕ,𝕋,CarX{ℝ},CarZ{ℝ}}, Bool) do (id, t, _, _)
                false
            end
        )
    return world_dynamics, delays, controller_framework, init_states, times
end

function loop(pub::Publisher)
    freq = 20.0
    wd, delay_model, cf, init_status, times = make_ros_controllers(1:20, freq)
    (controllers, msg_qs) = make_controllers(cf, init_status, times[1])
    println(controllers)
    make_agent(id::ℕ)::AgentState = begin
        (x₀, z₀, u₀) = init_status[id]
        AgentState(
            state_queue=constant_queue(x₀, delay_model.obs),
            obs_queue=constant_queue(z₀, delay_model.obs),
            act_queue=constant_queue(u₀, delay_model.act),
            msg_queue=msg_qs[id],
            controller=controllers[id],
        )
    end
    N = 1
    agents = map(make_agent, SVector{N}(1:N))
    init_states = SVector{N}(init_status)
    #a::SVector{N,Tuple{CarX,CarZ,CarU}} = init_states
    loop_rate = Rate(freq)
    tmp = first(agents[1].msg_queue)
    transmition = MMatrix{N,N,typeof(tmp)}(undef)
    #print(a)
    Each = MVector{N}
    xs, zs, us = @unzip(MVector(init_states), Each{Tuple{CarX{Float64},CarZ{Float64},CarU{Float64}}})

    for t in 1:20
        if is_shutdown()
            error("FATAL: ROS node has been shutdown")
            exit(0)
        end
        rossleep(loop_rate)
        for i in 1:length(agents)  # unroll away dynamic dispatch at compile time!
            a = agents[i]
            x = xs[i]
            u = us[i]
            z = zs[i]
            ms = first(a.msg_queue)
            u, ms′ = control!(a.controller, x, z, ms)
            #transmition[:,i] = ms′
            ros_publish_actuation(pub, u, x)
        end

    end
    """
    loop_rate = Rate(20.0)
    for t in 1:20
        if is_shutdown()
            error("FATAL: ROS node has been shutdown")
        end
        rossleep(loop_rate) # handle callbacks here - busy loop. Potentially very inefficient
        # at this point, messages should be handled, - run the simulation framework
        print(xs[0])
        xs[0] = x_cur # assign the current observation

        # observe, control, and send messages
        transmition = MMatrix{N,N,Msg}(undef)  # indexed by (receiver, sender)
        for i in 1:length(agents)  # unroll away dynamic dispatch at compile time!
            a = agents[i]
            x = pushpop!(a.state_queue, xs[i])
            z = pushpop!(a.obs_queue, zs[i])
            ms = first(a.msg_queue)
            u, ms′ = control!(a.controller, x, z, ms)
            us[i] = pushpop!(a.act_queue, u)
            transmition[:,i] = ms′
        end
        # receive messagees and update world states
        for i in 1:length(agents)
            a = agents[i]
            pushpop!(a.msg_queue, convert(Vector{Msg}, transmition[i,:]))
            ros_publish_actuation(pub, act, cur_state) # needs to be obtained from ROS
            zs[i] = zs[i] # do not update the wall pos for now #obs_forward(world_dynamics.obs_dynamics[i], xs[i], zs[i], t)
        end
    end
    """
    print("finished loop")
end

function main()
    # ROS related initialization
    init_node("rosjl_example")
    pub = Publisher{AckermannCurvatureDriveMsg}("ackermann_curvature_drive", queue_size=1)
    sub = Subscriber{Odometry}("odom", odom_callback,  queue_size=10)
    # simulation initialization
    loop(pub)
end
end


text:
using Dates
using DataFrames

include("mysql_db_utils.jl")
include("sqlite_db_utils.jl")

# This function takes a single quote and replaces it with
# two single quotes. This is what MySQL requires
db_clean_string(str) = replace(str, "'" => "''")

"""
    assemble_vals(data_values)

Given a dictionary containg (:column=>value) return a single string properly
formatted for a MySQL insert. E.g MySQL requires CHAR or
other non-numeric values be passed with single quotes around them.

"""
function assemble_cols_and_vals(data_values::Dict{Symbol, T}) where T
    vals_single_quotes = Vector{Any}(undef, length(data_values))        # put values in Array to be joined
    table_cols_backticks = Vector{String}(undef, length(data_values))

    for (i, (key,val)) in enumerate(data_values)
        table_cols_backticks[i] = string("`", key, "`")
        if typeof(val) <: Number && !ismissing(val)
            vals_single_quotes[i] = val
        elseif ismissing(val)
            vals_single_quotes[i] = "NULL"
        elseif val == nothing
            vals_single_quotes[i] = "NULL"
        elseif isa(val, Date)
            vals_single_quotes[i] = string("'", val, "'")
        else
            vals_single_quotes[i] = string("'", db_clean_string(val), "'")
        end
    end

    cols_string = join(table_cols_backticks, ", ")
    value_string = string("(", join(vals_single_quotes, ", "), ")")

    return cols_string, value_string
end

"""
    assemble_vals(data_values, column_names)

Given a Dict of values and the column names, return a single string properly
formatted for a MySQL INSERT. E.g MySQL requires CHAR or
other non-numeric values be passed with single quotes around them.
"""
function assemble_vals(data_values::Dict{Symbol, T}, column_names::Array{Symbol, 1}) where T
    vals_single_quotes = Vector{Any}()        # put values in Array to be joined

    for k in column_names
        if typeof(data_values[k]) <: Number
            push!(vals_single_quotes, data_values[k])
        elseif data_values[k] == nothing
            push!(vals_single_quotes, "NULL")
        elseif isa(data_values[k], Date)
            push!(vals_single_quotes, string("'", data_values[k], "'"))
        else
            push!(vals_single_quotes, string("'", db_clean_string(data_values[k]), "'"))
        end
    end
    value_string = string("(", join(vals_single_quotes, ", "), ")")
    return value_string
end

"""
    assemble_vals(data_values)

Given a dictionary containg (:column=>value), return a single string properly
formatted for a MySQL SELECT. E.g MySQL requires CHAR or
other non-numeric values be passed with single quotes around them.

"""
function assemble_cols_and_vals_string(data_values::Dict{Symbol, T}, op = "AND") where T
    val_single_quotes::Any = nothing        # put values in Array to be joined
    col_backticks::String = ""
    select_string_array = Array{String}(undef, length(data_values))

    for (i, (key,val)) in enumerate(data_values)
        col_backticks = string("`", key, "`")
        if typeof(val) <: Number && !ismissing(val)
            val_single_quotes=val
        elseif val == nothing || ismissing(val)
            val_single_quotes = "NULL"
        elseif isa(val, Date)
            val_single_quotes=string("'", val, "'")
        else
            val_single_quotes=string("'", db_clean_string(val), "'")
        end

        select_string_array[i] = string(col_backticks, "=" , val_single_quotes)
    end
    select_string = join(select_string_array, string(" ", op, " "))
    return select_string
end

function assemble_cols_and_vals_select(data_values::Dict{Symbol, T}, op = "AND") where T
    assemble_cols_and_vals_string(data_values, "AND")
end

"""
    assemble_cols(data_values::DataFrame)
Given a DataFrame, returns a column name string formatted for an insert/load statement
"""
function assemble_cols(data_values::DataFrame)
    col_str = ""
    for col in getfield(data_values, :colindex).names
        col_str *= string(col) * ","
    end
    return col_str[1:end-1]
end

#*****************
# colname_dict
#*****************
"""
    colname_dict_(con)
Return a dictionary maping tables and their columns for a given
MySQL-connection/SQLite-database
"""
function colname_dict(con)

    tables_query = select_all_tables(con)
    colname_dict = Dict{Symbol, Array{Symbol, 1}}()

    for table in tables_query
        cols_query = select_columns(con, table)
        cols = [symbol(c) for c in cols_query]
        colname_dict[symbol(table)] = cols
    end

    return colname_dict
end

#*****************
# select
#*****************
"""
    select_(con, colnames, tablename, data_values)

Perform: SELECT colnames tablename WHERE keys(data_values)=values(data_values)
"""
function db_select(con, colnames, tablename, data_values::Dict{Symbol, T}) where T
    select_cols_backticks = [string("`", x, "`") for x in colnames]
    select_cols_string = join(select_cols_backticks, ", ")
    select_string = assemble_cols_and_vals_select(data_values)
    sel = db_query(con, "SELECT $select_cols_string FROM `$tablename` WHERE $select_string;")
    sel
end

"""
    col_match(con, tablename, data_values)
Checks if each column in the dataframe has a matching column in the table.
"""
function col_match(con, tablename::String, data_values::DataFrame)
    cols = String.(getfield(data_values, :colindex).names)
    col_match(con, tablename, cols)
end

"""
    col_match(con, tablename, col_names)
Checks if each column in the csv/data frame has a matching column in the table.
"""
function col_match(con, tablename::String, col_names::Vector{String})
    all_match = true

    table_cols = select_columns(con, tablename)

    for col in col_names
        this_match = false
        for tc in table_cols
            if tc == string(col)
                this_match = true
                break
            end
        end
        if this_match == false
            all_match = false
            break
        end
    end

    return all_match
end


text:
using StructArrays: StructArray
using ComponentArrays: ComponentArrays
using LinearAlgebra: diagm, mul!
using NQCDynamics: RingPolymers

function set_quantum_derivative! end

function calculate_density_matrix_propagator!(sim::Simulation, v)
    V = sim.method.density_propagator

    V .= diagm(sim.calculator.eigen.values)
    for I in eachindex(v)
        @. V -= im * v[I] * sim.calculator.nonadiabatic_coupling[I]
    end
    return V
end

function calculate_density_matrix_propagator!(sim::RingPolymerSimulation, v)
    V = sim.method.density_propagator
    centroid_v = RingPolymers.get_centroid(v)

    V .= diagm(sim.calculator.centroid_eigen.values)
    for I in eachindex(centroid_v)
        @. V -= im * centroid_v[I] * sim.calculator.centroid_nonadiabatic_coupling[I]
    end
    return V
end

function commutator!(out, A, B, tmp)
    mul!(out, A, B)
    mul!(tmp, B, A)
    out .-= tmp
    return nothing
end

get_quantum_subsystem(u::ComponentArrays.ComponentVector{T}) where {T} =
    StructArray{Complex{T}}((u.σreal, u.σimag))


text:
module CSL

end


text:
using Flux
using Zygote
using MAT
#using Plots
using Statistics
using Distributions
using LinearAlgebra
using Random
using BSON
using JLD
#using PyPlot
using CUDA
#using Base.Iterators: partition

#gpu_id = 1
#CUDA.device!(gpu_id)

#****************** load data ******************
mf = matopen("../../Data/Data_processed_36Vol_NoCrop.mat")
Xall = read(mf, "Xtrain");
Yall = read(mf, "Ytrain");
close(mf)

for idx_subj = 1:size(Xall,5)
    Xall[:,:,:,:,idx_subj] = Xall[:,:,:,:,idx_subj]./maximum(Xall[:,:,:,:,idx_subj]);
    Yall[:,:,:,:,idx_subj] = Yall[:,:,:,:,idx_subj]./maximum(Yall[:,:,:,:,idx_subj]);
end

Xall = Xall[:,:,37:136,:,:];
Yall = Yall[:,:,37:136,:,:];

Sall = size(Xall);
Xall = reshape(Xall,Sall[1],Sall[2],1,Sall[3]*Sall[4],Sall[5]);

Sall = size(Yall);
Yall = reshape(Yall,Sall[1],Sall[2],1,Sall[3]*Sall[4],Sall[5]);
NVol = Sall[4] # number of volumes, 36

mf = matopen("../../Data/T1wBrain_processed_NoCrop.mat")
T1wBrain_resize_ds = read(mf, "T1wBrain_resize_ds");
close(mf)

for idx_subj = 1:size(T1wBrain_resize_ds,4)
    T1wBrain_resize_ds[:,:,:,idx_subj] = T1wBrain_resize_ds[:,:,:,idx_subj]./maximum(T1wBrain_resize_ds[:,:,:,idx_subj]);
end
T1wBrain_resize_ds = T1wBrain_resize_ds[:,:,37:136,:];

Xall_new = zeros(Float32,size(Xall,1),size(Xall,2),2,size(Xall,4),size(Xall,5));
Xall_new[:,:,1,:,:] = Xall;
Xall_new[:,:,2,:,:] = reshape(repeat(T1wBrain_resize_ds,1,1,NVol,1),size(Xall,1),size(Xall,2),1,size(Xall,4),size(Xall,5));

Xall = Xall_new;
Xall_new=nothing;

Sall = size(Xall);
Xtrain = reshape(Xall[:,:,:,:,[2,3,4]],Sall[1],Sall[2],Sall[3],Sall[4]*3);
Ytrain = reshape(Yall[:,:,:,:,[2,3,4]],Sall[1],Sall[2],1,Sall[4]*3);
idx_shuffle = shuffle(collect(1:size(Xtrain,4)));
Xtrain = Xtrain[:,:,:,idx_shuffle];
Ytrain = Ytrain[:,:,:,idx_shuffle];

# for noT1w only
#=
Xtrain_new = zeros(size(Ytrain));
Xtrain_new[:,:,1,:] = Xtrain[:,:,1,:];
Xtrain = Xtrain_new;
Xtrain_new = nothing;
=#

#Ytrain[:,:,1,:] = (Ytrain[:,:,1,:]-Xtrain[:,:,1,:])./(2.0).+(0.5)

Nim = size(Xtrain,1);

XValid = reshape(Xall[:,:,:,:,5],Sall[1],Sall[2],Sall[3],Sall[4]);
YValid = reshape(Yall[:,:,:,:,5],Sall[1],Sall[2],1,Sall[4]);

Xall = nothing;
Yall = nothing;

# for noT1w only
#=
XValid_new = zeros(size(YValid));
XValid_new[:,:,1,:] = XValid[:,:,1,:];
XValid = XValid_new;
XValid_new = nothing;
=#

# put data on GPU
Xtrain = gpu(Xtrain);
Ytrain = gpu(Ytrain);
XValid = gpu(XValid);
YValid = gpu(YValid);

# Unet define
#****************** neural network ******************
channels_in=size(Xtrain,3);
channels_out=1;
kernelsize = (3, 3);

BatchNormWrap(out_ch) = BatchNorm(out_ch)

UNetConvBlock(in_chs, out_chs, kernel = kernelsize) =
    Chain(Conv(kernel, in_chs=>out_chs,pad = (1, 1);init=Flux.kaiming_normal),
	BatchNormWrap(out_chs),
	x->relu.(x))

ConvDown(in_chs,out_chs,kernel = (2,2)) = Conv(kernel,in_chs=>out_chs,stride=(2,2);init=Flux.kaiming_normal)

UNetUpBlock(in_chs, out_chs; kernel = (2,2)) = ConvTranspose(kernel, in_chs=>out_chs,stride=(2, 2);init=Flux.kaiming_normal)

HP_all = load("HyperParsN50ForDeEncoder.jld")

bs_list = HP_all["bs_list"]
li_list = HP_all["li_list"]
ldf_list = HP_all["ldf_list"]
lds_list = HP_all["lds_list"]
NpoolList = HP_all["NpoolList"]
NchList = HP_all["NchList"]

Nsample = 50;

loss_valid = zeros(Float32,Nsample)
loss_train = zeros(Float32,Nsample)

for idx_randsearch = 1:Nsample

channels_out=1;

Nch_init = NchList[idx_randsearch];
batch_size = bs_list[idx_randsearch];
Nepochs = 30;
LR_initial = li_list[idx_randsearch];
LR_decay_factor = ldf_list[idx_randsearch];
LR_decay_step = lds_list[idx_randsearch];
Npool = NpoolList[idx_randsearch];

println("RandomSearch No.",idx_randsearch,":")
println("  Nch_init        = ",Nch_init)
println("  Npool           = ",Npool)
println("  batch_size      = ",batch_size)
println("  LR_initial      = ",LR_initial)
println("  LR_decay_factor = ",LR_decay_factor)
println("  LR_decay_step   = ",LR_decay_step)

if Npool == 2
r = SkipConnection(Chain(
       UNetConvBlock(channels_in, Nch_init),
       UNetConvBlock(Nch_init, Nch_init), #op
       SkipConnection(Chain(ConvDown(Nch_init,Nch_init),
                            UNetConvBlock(Nch_init, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#x1
                            SkipConnection(Chain(ConvDown(Nch_init*2,Nch_init*2),
                                                 UNetConvBlock(Nch_init*2, Nch_init*4),
                                                 UNetConvBlock(Nch_init*4, Nch_init*4),#x2
                                                 UNetUpBlock(Nch_init*4,Nch_init*2)),
                                         (mx,x)->cat(mx,x,dims=3)),
                            UNetConvBlock(Nch_init*4, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#up_x3
                            UNetUpBlock(Nch_init*2,Nch_init)),
                      (mx,x)->cat(mx,x,dims=3)),
       UNetConvBlock(Nch_init*2, Nch_init),
       UNetConvBlock(Nch_init, Nch_init),#up_x5
       Conv((1, 1), Nch_init=>channels_out;init=Flux.kaiming_normal)),
       (mx,x)->reshape(mx[:,:,1,:]+x[:,:,1,:],Nim,Nim,1,size(x,4)))
elseif Npool == 3

r = SkipConnection(Chain(
       UNetConvBlock(channels_in, Nch_init),
       UNetConvBlock(Nch_init, Nch_init), #op
       SkipConnection(Chain(ConvDown(Nch_init,Nch_init),
                            UNetConvBlock(Nch_init, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#x1
                            SkipConnection(Chain(ConvDown(Nch_init*2,Nch_init*2),
                                                 UNetConvBlock(Nch_init*2, Nch_init*4),
                                                 UNetConvBlock(Nch_init*4, Nch_init*4),#x2
                                                 SkipConnection(Chain(ConvDown(Nch_init*4,Nch_init*4),
                                                                      UNetConvBlock(Nch_init*4, Nch_init*8),
                                                                      UNetConvBlock(Nch_init*8, Nch_init*8),#x3
                                                                      UNetUpBlock(Nch_init*8,Nch_init*4)),
                                                                (mx,x)->cat(mx,x,dims=3)),
                                                 UNetConvBlock(Nch_init*8, Nch_init*4),
                                                 UNetConvBlock(Nch_init*4, Nch_init*4),#up_x2
                                                 UNetUpBlock(Nch_init*4,Nch_init*2)),
                                           (mx,x)->cat(mx,x,dims=3)),
                            UNetConvBlock(Nch_init*4, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#up_x3
                            UNetUpBlock(Nch_init*2,Nch_init)),
                      (mx,x)->cat(mx,x,dims=3)),
       UNetConvBlock(Nch_init*2, Nch_init),
       UNetConvBlock(Nch_init, Nch_init),#up_x5
       Conv((1, 1), Nch_init=>channels_out;init=Flux.kaiming_normal)),
       (mx,x)->reshape(mx[:,:,1,:]+x[:,:,1,:],Nim,Nim,1,size(x,4)))

else
r = SkipConnection(Chain(
       UNetConvBlock(channels_in, Nch_init),
       UNetConvBlock(Nch_init, Nch_init), #op
       SkipConnection(Chain(ConvDown(Nch_init,Nch_init),
                            UNetConvBlock(Nch_init, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#x1
                            SkipConnection(Chain(ConvDown(Nch_init*2,Nch_init*2),
                                                 UNetConvBlock(Nch_init*2, Nch_init*4),
                                                 UNetConvBlock(Nch_init*4, Nch_init*4),#x2
                                                 SkipConnection(Chain(ConvDown(Nch_init*4,Nch_init*4),
                                                                      UNetConvBlock(Nch_init*4, Nch_init*8),
                                                                      UNetConvBlock(Nch_init*8, Nch_init*8),#x3
                                                                      SkipConnection(Chain(ConvDown(Nch_init*8,Nch_init*8),
                                                                                UNetConvBlock(Nch_init*8, Nch_init*16),
                                                                                UNetConvBlock(Nch_init*16, Nch_init*16),#x3
                                                                                UNetUpBlock(Nch_init*16,Nch_init*8)),
                                                                                (mx,x)->cat(mx,x,dims=3)),
                                                                      UNetConvBlock(Nch_init*16, Nch_init*8),
                                                                      UNetConvBlock(Nch_init*8, Nch_init*8),#up_x2
                                                                      UNetUpBlock(Nch_init*8,Nch_init*4)),
                                                                (mx,x)->cat(mx,x,dims=3)),
                                                 UNetConvBlock(Nch_init*8, Nch_init*4),
                                                 UNetConvBlock(Nch_init*4, Nch_init*4),#up_x2
                                                 UNetUpBlock(Nch_init*4,Nch_init*2)),
                                           (mx,x)->cat(mx,x,dims=3)),
                            UNetConvBlock(Nch_init*4, Nch_init*2),
                            UNetConvBlock(Nch_init*2, Nch_init*2),#up_x3
                            UNetUpBlock(Nch_init*2,Nch_init)),
                      (mx,x)->cat(mx,x,dims=3)),
       UNetConvBlock(Nch_init*2, Nch_init),
       UNetConvBlock(Nch_init, Nch_init),#up_x5
       Conv((1, 1), Nch_init=>channels_out;init=Flux.kaiming_normal)),
       (mx,x)->reshape(mx[:,:,1,:]+x[:,:,1,:],Nim,Nim,1,size(x,4)))

end

r = gpu(r)
#****************** Data loader and Train/Loss function ******************
train_loader = Flux.Data.DataLoader(Xtrain,Ytrain, batchsize=batch_size,shuffle=true)
loss(x,y) = Flux.mse(r(x),y)
validationloss() = mean(Flux.mse(r(XValid[:,:,:,(idx-1)*100+1:idx*100]),YValid[:,:,:,(idx-1)*100+1:idx*100]) for idx=1:round(Int,size(XValid,4)/100))
trainingloss() = mean(Flux.mse(r(Xtrain[:,:,:,(idx-1)*100+1:idx*100]),Ytrain[:,:,:,(idx-1)*100+1:idx*100]) for idx=1:round(Int,size(Xtrain,4)/100))
#trainingloss() = Flux.mse(r(Xtrain),Ytrain)
#validationloss() = Flux.mse(r(XValid),YValid)

#****************** Training ******************
loss_train = zeros(Float32,Nepochs);
loss_valid = zeros(Float32,Nepochs);
#NepochStop = 0;
for epoch in 1:Nepochs
  learn_rate = LR_initial*exp(-LR_decay_factor*floor((epoch-1)./LR_decay_step))
  opt = ADAM(learn_rate)
  @time Flux.train!(loss, Flux.params(r), train_loader, opt)
  #BSON.@save string("Results/model_resnet_ptx_processed_7dir_gpu_epoch",string(epoch),".bson") model
  #CUDA.unsafe_free!(u)
  #u = gpu(uc)
  loss_valid[epoch] = validationloss()

  println("Training epoch No.",epoch,", loss_valid = ", loss_valid[epoch])

end

end

#JLD.save(string("Results/Loss_origresnet_processed_36Vol_T1W_L1_NestedCV_124_5_N1to20.jld"),"loss_train",loss_train,"loss_valid",loss_valid)


text:
@inline xorshift_rotl(x::UInt64, k::Int) = (x >>> (0x3f & -k)) | (x << (0x3f & k))
@inline xorshift_rotl(x::UInt32, k::Int) = (x >>> (0x1f & -k)) | (x << (0x1f & k))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

touch(joinpath(@__DIR__, "deps.jl"))

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
include("hooks.jl")
include("stop_conditions.jl")
include("run.jl")
include("experiment.jl")


text:
using SamplingHS3
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@testset "Adaptive MCMC" begin include("test_adaptive_mcmc.jl") end
@testset "HR Tests" begin include("test_hr.jl") end
@testset "MH Tests" begin include("test_mh.jl") end
@testset "NS Integration" begin include("test_ns_integration_01.jl") end
@testset "SamplingHS3 Utils 01" begin include("test_utils_01.jl") end

#@testset "HR Tests" begin include("test_hr.jl") end


text:
#InspectDR: Draw legends with Cairo layer.
#-------------------------------------------------------------------------------

#==
===============================================================================#

#NOTE: linelength & textoffset are absolute values
function legend_renderitem(ctx::CairoContext, wfrm::DWaveform, x::Float64, y::Float64,
	linelength::Float64, textoffset::Float64)
	linestart = x+textoffset
	lineend = linestart+linelength

Cairo.save(ctx)
	setlinestyle(ctx, LineStyle(wfrm.line))
	drawline(ctx, Point2D(linestart, y), Point2D(lineend, y))
	drawglyph_safe(ctx, wfrm, Point2D(linestart+linelength/2, y))
Cairo.restore(ctx)
	x = lineend+textoffset #Compute new x
	render(ctx, wfrm.id, Point2D(x, y), align=ALIGN_VCENTER|ALIGN_LEFT)
end

function render_legend(ctx::CairoContext, rstrip::RStrip2D, display_data::Vector,
		lyt::PlotLayout, bbleg::BoundingBox)
	Cairo.save(ctx)
	drawrectangle(ctx, bbleg, lyt.frame_legend)
	Cairo.restore(ctx)

	Cairo.save(ctx)
	setclip(ctx, bbleg)
	setfont(ctx, lyt.font_legend)
	(w, h) = textextents_wh(ctx, "M")
	ypitch = h*(1+lyt.valloc_legenditemsp)

	#Compute absolute values:
	textoffset = lyt.hoffset_legendtext*w

	y = bbleg.ymin + h/2
	for d in display_data
		if d.strip != rstrip.istrip; continue; end
		if "" == d.id; continue; end
		legend_renderitem(ctx, d, bbleg.xmin, y, lyt.halloc_legendlineseg, textoffset)
		y += ypitch
	end
	Cairo.restore(ctx)
end

function render_colorscale(ctx::CairoContext, rstrip::RStrip2D, strip::GraphStrip, lyt::PlotLayout)
	colorscale = strip.colorscale

	#Compute bounding box of colorscale:
	xleft = rstrip.bb.xmax + lyt.halloc_right
	bb = BoundingBox(xleft, xleft+lyt.halloc_colorscale, rstrip.bb.ymin, rstrip.bb.ymax)

	#Draw color scale itself:
	N = length(colorscale.e)
	Δ = bb.ymin - bb.ymax
	ystart = bb.ymax #start at bottom (ymax)
	Cairo.save(ctx)
	for i in 1:N
		yend = bb.ymax + Δ * (i/N)
		cbb = BoundingBox(bb.xmin, bb.xmax, ystart, yend)
		drawrectangle(ctx, cbb, colorscale.e[i])
		ystart = yend
	end
	Cairo.restore(ctx)

	#Draw frame:
	Cairo.save(ctx)
	zscale = strip.zscale; zext = strip.zext
	zlines = gridlines(zscale, zext.min, zext.max, true, true, zscale.tgtmajor, zscale.tgtminor)
	render_zticks(ctx, rstrip, lyt, bb, zext, zlines, zscale)
	drawrectangle(ctx, bb, lyt.frame_colorscale)
	Cairo.restore(ctx)
	return
end

function render_legends(ctx::CairoContext, rplot::RPlot2D, plot::Plot2D)
	lyt = plot.layout.values #WANTCONST

	if lyt.enable_colorscale
		for rstrip in rplot.strips
			render_colorscale(ctx, rstrip, plot.strips[rstrip.istrip], lyt)
		end
	end
	if lyt.enable_legend
		legend_left = rplot.bb.xmax - lyt.halloc_legend
		for rstrip in rplot.strips
			bbleg = BoundingBox(legend_left, rplot.bb.xmax, rstrip.bb.ymin, rstrip.bb.ymax)
			render_legend(ctx, rstrip, plot.display_data, lyt, bbleg)
		end
	end
end

#Last line