text:
module VerificationData

export create_dicts

using PyCall
using Serialization

const sdf = PyNULL()
const sh = PyNULL()

function __init__()
    copy!(sdf, pyimport("sdf"))
    copy!(sh, pyimport("sdf_helper"))
end

function convert_it(read_only_numpy_arr::PyObject)
    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pyarray.jl#L14-L24

    # instead of PyBUF_ND_STRIDED =  Cint(PyBUF_WRITABLE | PyBUF_FORMAT | PyBUF_ND | PyBUF_STRIDES)
    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pybuffer.jl#L113
    pybuf = PyBuffer(read_only_numpy_arr, PyCall.PyBUF_FORMAT | PyCall.PyBUF_ND | PyCall.PyBUF_STRIDES)

    T, native_byteorder = PyCall.array_format(pybuf)
    sz = size(pybuf)
    strd = PyCall.strides(pybuf)
    length(strd) == 0 && (sz = ())
    N = length(sz)
    isreadonly = pybuf.buf.readonly==1
    info = PyCall.PyArray_Info{T,N}(native_byteorder, sz, strd, pybuf.buf.buf, isreadonly, pybuf)

    # See https://github.com/JuliaPy/PyCall.jl/blob/master/src/pyarray.jl#L123-L126
    PyCall.PyArray{T,N}(read_only_numpy_arr, info)
end

function filereader(path; convert=false)
    sdf.read(path, convert=convert)
end

convert_variable(py_data::Tuple) = Array.(convert_it.(py_data))
convert_variable(py_data::PyObject) = Array(convert_it(py_data))
convert_variable(py_data) = py_data

function get_grid(py_obj, data_size)
    grid = convert_it.(py_obj.grid.data)
    if data_size[1] == length.(grid)[1]
        return Array.(grid)
    else
        grid_mid = Array.(convert_it.(py_obj.grid_mid.data))
        return grid_mid
    end
end

function create_dicts(fn)
    fr = filereader(fn)
    header = fr.Header
    variables = filter!(k->!startswith(k, "__"), String.(keys(fr)))
    filter!(k->k≠"Header", variables)
    data = Dict{String, Any}()
    grids = Dict{String, NTuple{3,Vector}}()
    units = Dict{String, Union{String,Tuple{String},NTuple{3,String}}}()
    for v in variables
        sdf_var = getproperty(fr, v)
        id = sdf_var.id
        if hasproperty(sdf_var, :data)
            data[id] = convert_variable(sdf_var.data)
            if hasproperty(sdf_var, :grid)
                grids[id] = get_grid(sdf_var, size(data[id]))
            end
        end
        if hasproperty(sdf_var, :units)
            units[id] = sdf_var.units
        end
    end

    return header, data, grids, units
end

function create_dicts()
    serialize("0002.jl", create_dicts("0002.sdf"))
end

end  # module VerificationData


text:
# This is an example of how to implement attention in Mill
# Since attention reduce set of vectors to a single vector,
# it is an aggregation and as such, it belongs to aggregation layer
using Mill
using Flux
using Transformers
using Transformers.Basic: MultiheadAttention

""" BagAttention uses `f` to convert the input and `a` to create the attention"""
struct BagMultiheadAttention{F,A}
	mh::F
	agg::A
end

# remove agg_a from trainable parameters to keep the one in handling missing intact
Flux.@functor BagMultiheadAttention

function (m::BagMultiheadAttention)(ds, bags)
	xs = map(b -> m.mh(ds.data[:,b],ds.data[:,b],ds.data[:,b]), filter(!isempty, bags.bags))
	ArrayNode(m.agg(reduce(hcat, xs), bags))
end

function BagMultiheadAttention(input_dim::Int, nheads::Int, hidden_dim::Int, odim::Int)
		BagMultiheadAttention(
			MultiheadAttention(nheads, input_dim, hidden_dim, odim; future=false, pdrop = 0),
			SegmentedMeanMax(odim))
end

#Let's create an absolutely dummy dataset
ds = BagNode(ArrayNode(randn(2,5)), [1:2,2:5,0:-1])

#try it with reflectinmodel
model = reflectinmodel(ds,
	d -> Dense(d, 4, selu),
	d -> BagMultiheadAttention(d, 3, 3, 5),
	)

model(ds)
gradient(() -> sum(sin.(model(ds).data)), Flux.params(model))

ext: jl


text:
# ----------------------------------------------------------------------------------- #
# Copyright [c] 2016 Varnerlab
# School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14853

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

function Balances(t,x,dxdt_vector,data_dictionary)

#-------------------------------------------
# COMPLEMENT SPECIES
#-------------------------------------------
Initiator2          = x[1];   #   Initiator complex 2, lectin pathway
C4                  = x[2];   #   C4 protein
C2                  = x[3];   #   C2 protein
C4a                 = x[4];   #   C4a protein
C4b                 = x[5];   #   C4b protein
C2a                 = x[6];   #   C2a protein
C2b                 = x[7];   #   C2b protein
C3                  = x[8];   #   C3 protein
C3b                 = x[9];   #   C3b protein
C4bC2a              = x[10];  #   Classical/Lectin C3 convertase
C3Convertase2       = x[11];  #   Alternate C3 convertase
C4bC2aC3b           = x[12];  #   Classical C5 convertase
C5Convertase2       = x[13];  #   Alternate C5 convertase
C5                  = x[14];  #   C5 protein
C5a                 = x[15];  #   C5a protein
C5b                 = x[16];  #   C5b protein
C4BP                = x[17];  #   C4 binding protein
FactorH             = x[18];  #   Factor H
C3a                 = x[19];  #   C3a protein

if (Initiator2<1e-3)
  Initiator2 = 0.0
end

if (C5<0)
  C5 = 0.0
end

if (C3<0)
  C3 = 0.0
end

#-------------------------------------------
# parameter_array
#-------------------------------------------
parameter_array = data_dictionary["PARAMETER_ARRAY"]
k_Initiator2_C4             = 1.0*parameter_array[1];
Km_Initiator2_C4            = 1.0*parameter_array[2];
k_Initiator2_C2             = 1.0*parameter_array[3];
Km_Initiator2_C2            = 1.0*parameter_array[4];

n_In2_C4                    = parameter_array[5];
n_In2_C2                    = parameter_array[6];
alpha_In2_C4                = parameter_array[7];
alpha_In2_C2                = parameter_array[8];

n_C5a                       = parameter_array[9];
n_C3a                       = parameter_array[10];
k_C4bC2a                    = parameter_array[11];
k_C3b_basal                 = parameter_array[12];
k_C3Convertase2             = parameter_array[13];
k_C5Convertase2             = parameter_array[14];
k_C5_C4bC2a_bind            = parameter_array[15];
k_C5_conversion             = parameter_array[16];
Km_C5_conversion            = parameter_array[17];
k_C5_conversion_alternate   = parameter_array[18];
Km_C5_conversion_alternate  = parameter_array[19];
k_inhibit_C4bC2a            = parameter_array[20];
k_inhibit_C3Convertase2     = parameter_array[21];
k_cat_C4bC2a                = parameter_array[22];
Km_C4bC2a                   = parameter_array[23];
k_cat_C3convertase2         = parameter_array[24];
Km_C3Convertase2            = parameter_array[25];
k_degradationC3a            = parameter_array[26];
kC5inhibit                  = parameter_array[27];
k_degradationC5a            = parameter_array[28];

#-------------------------------------------
# REACTIONS
#-------------------------------------------

# Formation Reactions
rV = zeros(14)
rV[1]  = k_Initiator2_C4*Initiator2*C4/(Km_Initiator2_C4+C4);                                             #   Lectin Pathway initiator generating C4a and C4b from C4
rV[2]  = k_Initiator2_C2*Initiator2*C2/(Km_Initiator2_C2+C2);                                             #   Lectin Pathway generating C2a and C2b from C2
rV[3]  = k_C4bC2a*C4b*C2a;                                                                                #   Formation of classical/lectin C3 convertase
rV[4]  = k_C3b_basal*C3;                                                                                #   Formation of basal C3b and C3a from alternate pathway
rV[5]  = k_C3Convertase2*C3b;                                                                             #   Formation of alternate C3 convertase
rV[6]  = k_C5Convertase2*C3Convertase2*C3b;                                                               #   Formation of alternate C5 convertase
rV[7]  = k_C5_C4bC2a_bind*C4bC2a*C3b;                                                                     #   Formation of classical C5 convertase
rV[8] = k_C5_conversion*C4bC2aC3b*C5^n_C5a/(Km_C5_conversion^n_C5a+C5^n_C5a);                             #   Formation of C5a and C5b through classical C5 convertase
rV[9] = k_C5_conversion_alternate*C5Convertase2*C5/(Km_C5_conversion_alternate+C5);                       #   Formation of C5a and C5b through alternate C5 convertase

# Inhibition Reactions
rV[10] = k_inhibit_C4bC2a*C4b*C4BP;                                                                       #   Inhibition of C4b through C4BP
rV[11] = k_inhibit_C3Convertase2*C3Convertase2*FactorH;                                                   #   Inhibition of C3b through Factor H
rV[12] = kC5inhibit*C4bC2aC3b*C4BP;                                                                       #   Inhibition of Classical C5 Convertase by C4BP

# Formation of C3a and C3b from C3 Convertase
rV[13] = k_cat_C4bC2a*C4bC2a*C3^n_C3a/(Km_C4bC2a^n_C3a+C3^n_C3a);                                         #   Formation of C3a and C3b from Lectin C3 Convertase
rV[14] = k_cat_C3convertase2*C3Convertase2*C3/(Km_C3Convertase2+C3);                                      #   Formation of C3a and C3b from Alternate C3 Convertase

#Modified Initiator Steps -- Includes control functions
rV_modified = zeros(2)

# Setup logistics function -
control_vector = zeros(2)
# if (Initiator2 > 1e-5)
#
#   L = 1.0
#
#   alpha_C2 = alpha_In2_C2
#   alpha_C4 = alpha_In2_C4
#   k_C2 = n_In2_C2
#   k_C4 = n_In2_C4
#
#   control_vector[1] = L/(1+exp(-k_C4*(Initiator2 - alpha_C4)))
#   control_vector[2] = L/(1+exp(-k_C2*(Initiator2 - alpha_C2)))
#
# end

control_vector[1] = (Initiator2^(n_In2_C4))/(Initiator2^(n_In2_C4)+alpha_In2_C4^(n_In2_C4))
control_vector[2] = (Initiator2^(n_In2_C2))/(Initiator2^(n_In2_C2)+alpha_In2_C2^(n_In2_C2))

# @show control_vector,Initiator2
tau_C5 = 1.0
if (Initiator2>2e-5)
  tau_C5 = Initiator2;
end

rV_modified[1] = rV[1]*control_vector[1]                                                               #   Modified reaction rate for C4 conversion
rV_modified[2] = rV[2]*control_vector[2]                                                               #   Modified reaction rate for C2 conversion

# if (Initiator2 > 1e-5)
#   tau_C5 = control_vector[2]
# end

# rV_modified[1] = rV[1]*control_vector[1]
# rV_modified[1] = rV[2]*control_vector[2]

#-------------------------------------------
# MATERIAL BALANCES
#-------------------------------------------

dxdt_vector[1]   = 0;                                                             # Lectin Pathway Initiator
dxdt_vector[2]   = -rV_modified[1];                                               # C4
dxdt_vector[3]   = -rV_modified[2];                                               # C2
dxdt_vector[4]   = rV_modified[1];                                                # C4a
dxdt_vector[5]   = rV_modified[1]-rV[3];                                          # C4b
dxdt_vector[6]   = rV_modified[2]-rV[3];                                          # C2a
dxdt_vector[7]   = rV_modified[2];                                                # C2b
dxdt_vector[8]   = -rV[4]-rV[13]-rV[14];                                          # C3
dxdt_vector[9]   = rV[4]-rV[5]-rV[6]-rV[7]+rV[13]+rV[14];                         # C3b
dxdt_vector[10]  = rV[3]-rV[7]-rV[10];                                            # C4bC2a
dxdt_vector[11]  = rV[5]-rV[6]-rV[11];                                            # C3Convertase2
dxdt_vector[12]  = rV[7]-rV[12];                                                  # C4bC2aC3b
dxdt_vector[13]  = rV[6];                                                         # C5Convertase2
dxdt_vector[14]  = -1*(rV[8]+rV[9]);                                              # C5
dxdt_vector[15]  = tau_C5*(rV[8]+rV[9]-k_degradationC5a*C5a);                     # C5a
dxdt_vector[16]  = tau_C5*(rV[8]+rV[9]);                                          # C5b
dxdt_vector[17]  = -rV[10]-rV[12];                                                # C4BP
dxdt_vector[18]  = -rV[11];                                                       # FactorH
dxdt_vector[19]  = rV[4]+rV[13]+rV[14]-k_degradationC3a*C3a;                      # C3a

end


text:
using MyLib
using Base.Test

try
    ptms = new_tms()
    println("Start")
    times = libc_times(ptms)
    tms_utime = get_tms_utime(ptms)
    tms_stime = get_tms_stime(ptms)
    tms_cutime = get_tms_cutime(ptms)
    tms_cstime = get_tms_cstime(ptms)
    println("times: ", times)
    println("tms_utime: ", tms_utime)
    println("tms_stime: ", tms_stime)
    println("tms_cutime: ", tms_cutime)
    println("tms_cstime: ", tms_cstime)
    println("End")
    times = libc_times(ptms)
    tms_utime = get_tms_utime(ptms)
    tms_stime = get_tms_stime(ptms)
    tms_cutime = get_tms_cutime(ptms)
    tms_cstime = get_tms_cstime(ptms)
    println("times: ", times)
    println("tms_utime: ", tms_utime)
    println("tms_stime: ", tms_stime)
    println("tms_cutime: ", tms_cutime)
    println("tms_cstime: ", tms_cstime)
    delete_tms(ptms)
catch
    @test false
end

@test true


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpdateFirewallRuleProperties <: SwaggerModel
    startIpAddress::Any # spec type: Union{ Nothing, String } # spec name: startIpAddress
    endIpAddress::Any # spec type: Union{ Nothing, String } # spec name: endIpAddress

    function UpdateFirewallRuleProperties(;startIpAddress=nothing, endIpAddress=nothing)
        o = new()
        validate_property(UpdateFirewallRuleProperties, Symbol("startIpAddress"), startIpAddress)
        setfield!(o, Symbol("startIpAddress"), startIpAddress)
        validate_property(UpdateFirewallRuleProperties, Symbol("endIpAddress"), endIpAddress)
        setfield!(o, Symbol("endIpAddress"), endIpAddress)
        o
    end
end # type UpdateFirewallRuleProperties

const _property_map_UpdateFirewallRuleProperties = Dict{Symbol,Symbol}(Symbol("startIpAddress")=>Symbol("startIpAddress"), Symbol("endIpAddress")=>Symbol("endIpAddress"))
const _property_types_UpdateFirewallRuleProperties = Dict{Symbol,String}(Symbol("startIpAddress")=>"String", Symbol("endIpAddress")=>"String")
Base.propertynames(::Type{ UpdateFirewallRuleProperties }) = collect(keys(_property_map_UpdateFirewallRuleProperties))
Swagger.property_type(::Type{ UpdateFirewallRuleProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpdateFirewallRuleProperties[name]))}
Swagger.field_name(::Type{ UpdateFirewallRuleProperties }, property_name::Symbol) =  _property_map_UpdateFirewallRuleProperties[property_name]

function check_required(o::UpdateFirewallRuleProperties)
    true
end

function validate_property(::Type{ UpdateFirewallRuleProperties }, name::Symbol, val)
end


text:
module MeasureInterface

using DensityInterface
using FillArrays
using Static
using Random
using MappedArrays
import Random

abstract type AbstractMeasure end

include("primitive.jl")
include("basemeasure.jl")
include("dirac.jl")
include("domains.jl")
include("kernel.jl")
include("lebesgue.jl")
include("product.jl")
include("power.jl")
include("restrict.jl")
include("weighted.jl")

end


text:
include("lib_for_practic_10.jl")

#r = Robot()
#r = BorderRobot{CoordRobot}(r)
 function border_movements!(r,side::HorizonSide, side_part :: HorizonSide)

    flag = false
    count = 0

    while isborder(r,side_part) == false
        while isborder(r,side) ==  false
            move!(r,side)
            if (flag == false && isborder(r, side_part))
                flag = true
                count += 1
            elseif (flag == true && !isborder(r, side_part))
                flag = false
            end
        end

        move!(r, side_part)
        side = inverse(side)
    end

    return count

end

function border_count(r)
    x = movements!(r,West)
    y = movements!(r,Sud)

    ans = border_movements!(r, Ost, Nord)

    movements!(r,West)
    movements!(r,Sud)
    movements!(r,Nord, y)
    movements!(r,Ost, x)

    println(ans)
end

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# World containing a collection of shapes

"""
    World

Alias of `Vector{Shape}`, to store a list of [`Shape`](@ref).
"""
const World = Vector{Shape}

"""
    ray_intersection(ray::Ray, world::World)

Intersect a [`Ray`](@ref) with each [`Shape`](@ref) in [`World`](@ref) and return the nearest hit point.
"""
function ray_intersection(ray::Ray, world::World)
    hit = nothing
    for shape ∈ world
        last_hit = ray_intersection(ray, shape)
        last_hit === nothing && continue
        if hit === nothing
            hit = last_hit
        else
            last_hit.t/norm(last_hit.ray.dir) >= hit.t/norm(hit.ray.dir) && continue
            hit = last_hit
        end
    end
    hit
end

"""
    is_point_visible(world::World, point::Point, observer_pos::Point)

Tells if a particular [`Point`](@ref) in a [`World`](@ref) filled with [`Shape`](@ref) is visible from the observer position.
"""
function is_point_visible(world::World, point::Point, observer_pos::Point)
    direction = point - observer_pos

    ray = Ray(observer_pos, direction, 1f-2, 1f0, 0)

    return !any(shape -> quick_ray_intersection(ray, shape), world)
end


text:
function new(args::Dict)
  pkg_name = args["<pkg_name>"]
  license = args["<license>"]

  if license == "mit" ; license = "MIT" ; end

  new(pkg_name, license)
end

function new(pkg_name, license="null")

  initial_dir = pwd()

  cleaned_pkg_name = get_class_name(pkg_name)
  cleaned_pkg_name = replace(cleaned_pkg_name, ".jl", "")

  PkgDev.generate(cleaned_pkg_name, license, path=initial_dir)

  package_dir = "$initial_dir/$cleaned_pkg_name"

  cd(package_dir)
  init()
  cd(initial_dir)

  mv(package_dir, "$package_dir.jl")

end


text:
@testset "Neighborhoods" begin
  @testset "NormBall" begin
    # Euclidean metric
    ball = NormBall(T(1/2))
    @test evaluate(metric(ball), T[0], T[0]) ≤ radius(ball)
    @test evaluate(metric(ball), T[0], T[1]) > radius(ball)

    ball = NormBall(T(1))
    @test evaluate(metric(ball), T[0,0], T[0,0]) ≤ radius(ball)
    @test evaluate(metric(ball), T[0,0], T[1,0]) ≤ radius(ball)
    @test evaluate(metric(ball), T[0,0], T[0,1]) ≤ radius(ball)

    # Chebyshev metric
    ball = NormBall(T(1/2), Chebyshev())
    @test evaluate(metric(ball), T[0], T[0]) ≤ radius(ball)
    @test evaluate(metric(ball), T[0], T[1]) > radius(ball)

    for r in [1.,2.,3.,4.,5.]
      ball = NormBall(r, Chebyshev())
      for i in 0.0:1.0:r, j in 0.0:1.0:r
        @test evaluate(metric(ball), T[0,0], T[i,j]) ≤ radius(ball)
      end
    end

    ball = NormBall(T(1))
    if T == Float32
      @test sprint(show, ball) == "NormBall{Float32}(1.0, Euclidean(0.0))"
    else
      @test sprint(show, ball) == "NormBall{Float64}(1.0, Euclidean(0.0))"
    end
  end

  @testset "Ellipsoid" begin
    # 2D ellipsoid rotated 45 degrees in GSLIB convention
    ellipsoid = Ellipsoid(T[2,1], T[45], convention=GSLIB)

    # tests along main semiaxes, slightly below threshold
    @test evaluate(metric(ellipsoid), T[0,0], T(1.9) * T[√2/2, √2/2]) ≤ T(1)
    @test evaluate(metric(ellipsoid), T[0,0], T(0.9) * T[√2/2,-√2/2]) ≤ T(1)

    # tests along main semiaxes, slightly above threshold
    @test evaluate(metric(ellipsoid), T[0,0], T(2.1) * T[√2/2, √2/2]) > T(1)
    @test evaluate(metric(ellipsoid), T[0,0], T(1.1) * T[√2/2,-√2/2]) > T(1)

    # 3D ellipsoid rotated (45, -45, 0) in GSLIB convention
    ellipsoid = Ellipsoid(T[3,2,1], T[45,-45,0], convention=GSLIB)

    # tests along main semiaxes, slightly below threshold
    @test evaluate(metric(ellipsoid), T[0,0,0], T(2.9) * T[0.5,0.5,-√2/2]) ≤ T(1)
    @test evaluate(metric(ellipsoid), T[0,0,0], T(1.9) * T[√2/2,-√2/2,0.0]) ≤ T(1)
    @test evaluate(metric(ellipsoid), T[0,0,0], T(0.9) * T[0.5,0.5,√2/2]) ≤ T(1)

    # Tests along main semiaxes, slightly above threshold
    @test evaluate(metric(ellipsoid), T[0,0,0], T(3.1) * T[0.5,0.5,-√2/2]) > T(1)
    @test evaluate(metric(ellipsoid), T[0,0,0], T(2.1) * T[√2/2,-√2/2,0.0]) > T(1)
    @test evaluate(metric(ellipsoid), T[0,0,0], T(1.1) * T[0.5,0.5,√2/2]) > T(1)

    # 2D simple test of default convention
    d₁ = metric(Ellipsoid([1.,1.], [0.]))
    d₂ = metric(Ellipsoid([1.,2.], [0.]))
    @test evaluate(d₁, [1.,0.], [0.,0.]) == evaluate(d₁, [0.,1.], [0.,0.])
    @test evaluate(d₂, [1.,0.], [0.,0.]) != evaluate(d₂, [0.,1.], [0.,0.])

    # 3D simple test of default convention
    d₃ = metric(Ellipsoid([1.,.5,.5], [π/4,0.,0.]))
    @test evaluate(d₃, [1.,1.,0.], [0.,0.,0.]) ≈ √2
    @test evaluate(d₃, [-1.,1.,0.], [0.,0.,0.]) ≈ √8

    # test of intrinsic conventions
    gslib = metric(Ellipsoid([50.,25.,5.], [30.,-30.,30.], convention=GSLIB))
    tait  = metric(Ellipsoid([25.,50.,5.], [-π/6,-π/6,π/6], convention=TaitBryanIntr))
    euler = metric(Ellipsoid([50.,25.,5.], deg2rad.([-78,-41,-50]), convention=EulerIntr))
    lpf   = metric(Ellipsoid([50.,25.,5.], [78.,41.,50.], convention=Leapfrog))
    dm    = metric(Ellipsoid([50.,25.,5.], [78.,41.,50.], convention=Datamine))

    @test evaluate(gslib, [1.,0.,0.], [0.,0.,0.]) ≈ 0.1325707358356285
    @test evaluate(gslib, [0.,1.,0.], [0.,0.,0.]) ≈ 0.039051248379533283
    @test evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) ≈ 0.15132745950421558

    @test evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(tait, [0.,0.,1.], [0.,0.,0.])
    @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(dm, [0.,0.,1.], [0.,0.,0.])
    @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(lpf, [0.,0.,1.], [0.,0.,0.])
    @test evaluate(euler, [0.,0.,1.], [0.,0.,0.]) - evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) < 10^-3

    # test for https://github.com/JuliaEarth/GeoStats.jl/issues/197
    gslib = metric(Ellipsoid([1.0, 0.5, 0.3], [100., -10., -20.], convention = GSLIB))

    @test evaluate(gslib, [1.,0.,0.], [0.,0.,0.]) ≈ 1.233956165693094
    @test evaluate(gslib, [0.,1.,0.], [0.,0.,0.]) ≈ 2.14219475359467
    @test evaluate(gslib, [0.,0.,1.], [0.,0.,0.]) ≈ 3.1621930570302914

    # test of extrinsic conventions
    xtait  = metric(Ellipsoid([50.,25.,5.], [π,0,π/2], convention=TaitBryanExtr))
    xeuler = metric(Ellipsoid([50.,25.,5.], [-π/2,-π/2,-π/2], convention=EulerExtr))

    @test evaluate(xtait, [1.,0.,0.], [0.,0.,0.]) ≈ 0.20
    @test evaluate(xtait, [0.,1.,0.], [0.,0.,0.]) ≈ 0.04
    @test evaluate(xtait, [0.,0.,1.], [0.,0.,0.]) ≈ 0.02

    @test evaluate(xtait, [1.,0.,0.], [0.,0.,0.]) ≈ evaluate(xeuler, [1.,0.,0.], [0.,0.,0.])
    @test evaluate(xtait, [0.,1.,0.], [0.,0.,0.]) ≈ evaluate(xeuler, [0.,1.,0.], [0.,0.,0.])
    @test evaluate(xtait, [0.,0.,1.], [0.,0.,0.]) ≈ evaluate(xeuler, [0.,0.,1.], [0.,0.,0.])

    ellipsoid = Ellipsoid(T.((2,1)), T.((45,)), convention=GSLIB)
    if T == Float32
      @test sprint(show, ellipsoid) == "Ellipsoid{Float32}((2.0f0, 1.0f0), (45.0f0,), GSLIB)"
    else
      @test sprint(show, ellipsoid) == "Ellipsoid{Float64}((2.0, 1.0), (45.0,), GSLIB)"
    end
  end
end


text:
using FSimBase
using DifferentialEquations
using ComponentArrays
using Test
using LinearAlgebra
using DataFrames

function main()
    state0 = [1.0, 2.0]
    p = 1
    tf = 1.0
    Δt = 0.01
    @Loggable function dynamics!(dx, x, p, t)
        @log t
        @log x
        dx .= -p.*x
    end
    simulator = Simulator(
                          state0, dynamics!, p;
                          Problem=:ODE,
                          solver=Tsit5(),
                          tf=tf,
                         )
    # solve approach (automatically reinitialised)
    @time _df = solve(simulator; savestep=Δt)
    # interactive simulation
    ## step!
    reinit!(simulator)
    step!(simulator, Δt)
    @test simulator.integrator.t ≈ Δt
    ## step_until! (callback-like)
    ts_weird = 0:Δt:tf+Δt
    df_ = DataFrame()
    reinit!(simulator)
    @time for t in ts_weird
        flag = step_until!(simulator, t)  # flag == false if step is inappropriate
        if simulator.integrator.u[1] < 5e-1
            break
        else
            push!(simulator, df_, flag)  # push data only when flag == true
        end
    end
    println(df_[end-5:end, :])
    ## step_until!
    df = DataFrame()
    reinit!(simulator)
    @time for t in ts_weird
        step_until!(simulator, t)
        push!(simulator, df)  # flag=true is default
        # safer way:
        # flag = step_until!(simulator, t)
        # push!(simulator, df, flag)
        # or, equivalently,
        # push!(simulator, df, step_until!(simulator, t))  # compact form
    end
    println(df[end-5:end, :])
    @test norm(_df.sol[end].x - df.sol[end].x) < 1e-6
    @test simulator.integrator.t ≈ tf
end

@testset "minimal" begin
    main()
end


text:
@testset "3.1.2 (d x)^m (a+b log(c x^n))^p" begin
    (a, b, c, d, d1, d2, m, m1, m2, n, p, q, q1, q2, x, ) = @variables a b c d d1 d2 m m1 m2 n p q q1 q2 x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(d*x)^m*(a+b*log(c*x^n))^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*log(c*x)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*log(c*x), x, 1, -x^4/16 + (x^4*log(c*x))/4]
    @test_int [x^2*log(c*x), x, 1, -x^3/9 + (x^3*log(c*x))/3]
    @test_int [x^1*log(c*x), x, 1, -x^2/4 + (x^2*log(c*x))/2]
    @test_int [x^0*log(c*x), x, 1, -x + x*log(c*x)]
    @test_int [log(c*x)/x^1, x, 1, log(c*x)^2/2]
    @test_int [log(c*x)/x^2, x, 1, -x^(-1) - log(c*x)/x]
    @test_int [log(c*x)/x^3, x, 1, -1/(4*x^2) - log(c*x)/(2*x^2)]

    @test_int [x^3*log(c*x)^2, x, 2, x^4/32 - (x^4*log(c*x))/8 + (x^4*log(c*x)^2)/4]
    @test_int [x^2*log(c*x)^2, x, 2, (2*x^3)/27 - (2*x^3*log(c*x))/9 + (x^3*log(c*x)^2)/3]
    @test_int [x^1*log(c*x)^2, x, 2, x^2/4 - (x^2*log(c*x))/2 + (x^2*log(c*x)^2)/2]
    @test_int [x^0*log(c*x)^2, x, 2, 2*x - 2*x*log(c*x) + x*log(c*x)^2]
    @test_int [log(c*x)^2/x^1, x, 2, log(c*x)^3/3]
    @test_int [log(c*x)^2/x^2, x, 2, -2/x - (2*log(c*x))/x - log(c*x)^2/x]
    @test_int [log(c*x)^2/x^3, x, 2, -1/(4*x^2) - log(c*x)/(2*x^2) - log(c*x)^2/(2*x^2)]

    @test_int [x^3*log(c*x)^3, x, 3, (-3*x^4)/128 + (3*x^4*log(c*x))/32 - (3*x^4*log(c*x)^2)/16 + (x^4*log(c*x)^3)/4]
    @test_int [x^2*log(c*x)^3, x, 3, (-2*x^3)/27 + (2*x^3*log(c*x))/9 - (x^3*log(c*x)^2)/3 + (x^3*log(c*x)^3)/3]
    @test_int [x^1*log(c*x)^3, x, 3, (-3*x^2)/8 + (3*x^2*log(c*x))/4 - (3*x^2*log(c*x)^2)/4 + (x^2*log(c*x)^3)/2]
    @test_int [x^0*log(c*x)^3, x, 3, -6*x + 6*x*log(c*x) - 3*x*log(c*x)^2 + x*log(c*x)^3]
    @test_int [log(c*x)^3/x^1, x, 2, log(c*x)^4/4]
    @test_int [log(c*x)^3/x^2, x, 3, -6/x - (6*log(c*x))/x - (3*log(c*x)^2)/x - log(c*x)^3/x]
    @test_int [log(c*x)^3/x^3, x, 3, -3/(8*x^2) - (3*log(c*x))/(4*x^2) - (3*log(c*x)^2)/(4*x^2) - log(c*x)^3/(2*x^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3/log(c*x), x, 2, ExpIntegralEi(4*log(c*x))/c^4]
    @test_int [x^2/log(c*x), x, 2, ExpIntegralEi(3*log(c*x))/c^3]
    @test_int [x^1/log(c*x), x, 2, ExpIntegralEi(2*log(c*x))/c^2]
    @test_int [x^0/log(c*x), x, 1, LogIntegral(c*x)/c]
    @test_int [1/(x^1*log(c*x)), x, 2, log(log(c*x))]
    @test_int [1/(x^2*log(c*x)), x, 2, c*ExpIntegralEi(-log(c*x))]
    @test_int [1/(x^3*log(c*x)), x, 2, c^2*ExpIntegralEi(-2*log(c*x))]

    @test_int [x^3/log(c*x)^2, x, 3, (4*ExpIntegralEi(4*log(c*x)))/c^4 - x^4/log(c*x)]
    @test_int [x^2/log(c*x)^2, x, 3, (3*ExpIntegralEi(3*log(c*x)))/c^3 - x^3/log(c*x)]
    @test_int [x^1/log(c*x)^2, x, 3, (2*ExpIntegralEi(2*log(c*x)))/c^2 - x^2/log(c*x)]
    @test_int [x^0/log(c*x)^2, x, 2, -(x/log(c*x)) + LogIntegral(c*x)/c]
    @test_int [1/(x^1*log(c*x)^2), x, 2, -log(c*x)^(-1)]
    @test_int [1/(x^2*log(c*x)^2), x, 3, -(c*ExpIntegralEi(-log(c*x))) - 1/(x*log(c*x))]
    @test_int [1/(x^3*log(c*x)^2), x, 3, -2*c^2*ExpIntegralEi(-2*log(c*x)) - 1/(x^2*log(c*x))]

    @test_int [x^3/log(c*x)^3, x, 4, (8*ExpIntegralEi(4*log(c*x)))/c^4 - x^4/(2*log(c*x)^2) - (2*x^4)/log(c*x)]
    @test_int [x^2/log(c*x)^3, x, 4, (9*ExpIntegralEi(3*log(c*x)))/(2*c^3) - x^3/(2*log(c*x)^2) - (3*x^3)/(2*log(c*x))]
    @test_int [x^1/log(c*x)^3, x, 4, (2*ExpIntegralEi(2*log(c*x)))/c^2 - x^2/(2*log(c*x)^2) - x^2/log(c*x)]
    @test_int [x^0/log(c*x)^3, x, 3, -(x/(2*log(c*x)^2)) - x/(2*log(c*x)) + LogIntegral(c*x)/(2*c)]
    @test_int [1/(x^1*log(c*x)^3), x, 2, -1/(2*log(c*x)^2)]
    @test_int [1/(x^2*log(c*x)^3), x, 4, (1/2)*c*ExpIntegralEi(-log(c*x)) - 1/(2*x*log(c*x)^2) + 1/(2*x*log(c*x))]
    @test_int [1/(x^3*log(c*x)^3), x, 4, 2*c^2*ExpIntegralEi(-2*log(c*x)) - 1/(2*x^2*log(c*x)^2) + 1/(x^2*log(c*x))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*log(c*x^n))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*(a + b*log(c*x^n)), x, 1, -(b*n*x^4)/16 + (x^4*(a + b*log(c*x^n)))/4]
    @test_int [x^2*(a + b*log(c*x^n)), x, 1, -(b*n*x^3)/9 + (x^3*(a + b*log(c*x^n)))/3]
    @test_int [x^1*(a + b*log(c*x^n)), x, 1, -(b*n*x^2)/4 + (x^2*(a + b*log(c*x^n)))/2]
    @test_int [x^0*(a + b*log(c*x^n)), x, 2, a*x - b*n*x + b*x*log(c*x^n)]
    @test_int [(a + b*log(c*x^n))/x^1, x, 1, (a + b*log(c*x^n))^2/(2*b*n)]
    @test_int [(a + b*log(c*x^n))/x^2, x, 1, -((b*n)/x) - (a + b*log(c*x^n))/x]
    @test_int [(a + b*log(c*x^n))/x^3, x, 1, -(b*n)/(4*x^2) - (a + b*log(c*x^n))/(2*x^2)]

    @test_int [x^3*(a + b*log(c*x^n))^2, x, 2, (b^2*n^2*x^4)/32 - (b*n*x^4*(a + b*log(c*x^n)))/8 + (x^4*(a + b*log(c*x^n))^2)/4]
    @test_int [x^2*(a + b*log(c*x^n))^2, x, 2, (2*b^2*n^2*x^3)/27 - (2*b*n*x^3*(a + b*log(c*x^n)))/9 + (x^3*(a + b*log(c*x^n))^2)/3]
    @test_int [x^1*(a + b*log(c*x^n))^2, x, 2, (b^2*n^2*x^2)/4 - (b*n*x^2*(a + b*log(c*x^n)))/2 + (x^2*(a + b*log(c*x^n))^2)/2]
    @test_int [x^0*(a + b*log(c*x^n))^2, x, 3, -2*a*b*n*x + 2*b^2*n^2*x - 2*b^2*n*x*log(c*x^n) + x*(a + b*log(c*x^n))^2]
    @test_int [(a + b*log(c*x^n))^2/x^1, x, 2, (a + b*log(c*x^n))^3/(3*b*n)]
    @test_int [(a + b*log(c*x^n))^2/x^2, x, 2, (-2*b^2*n^2)/x - (2*b*n*(a + b*log(c*x^n)))/x - (a + b*log(c*x^n))^2/x]
    @test_int [(a + b*log(c*x^n))^2/x^3, x, 2, -(b^2*n^2)/(4*x^2) - (b*n*(a + b*log(c*x^n)))/(2*x^2) - (a + b*log(c*x^n))^2/(2*x^2)]

    @test_int [x^3*(a + b*log(c*x^n))^3, x, 3, (-3*b^3*n^3*x^4)/128 + (3*b^2*n^2*x^4*(a + b*log(c*x^n)))/32 - (3*b*n*x^4*(a + b*log(c*x^n))^2)/16 + (x^4*(a + b*log(c*x^n))^3)/4]
    @test_int [x^2*(a + b*log(c*x^n))^3, x, 3, (-2*b^3*n^3*x^3)/27 + (2*b^2*n^2*x^3*(a + b*log(c*x^n)))/9 - (b*n*x^3*(a + b*log(c*x^n))^2)/3 + (x^3*(a + b*log(c*x^n))^3)/3]
    @test_int [x^1*(a + b*log(c*x^n))^3, x, 3, (-3*b^3*n^3*x^2)/8 + (3*b^2*n^2*x^2*(a + b*log(c*x^n)))/4 - (3*b*n*x^2*(a + b*log(c*x^n))^2)/4 + (x^2*(a + b*log(c*x^n))^3)/2]
    @test_int [x^0*(a + b*log(c*x^n))^3, x, 4, 6*a*b^2*n^2*x - 6*b^3*n^3*x + 6*b^3*n^2*x*log(c*x^n) - 3*b*n*x*(a + b*log(c*x^n))^2 + x*(a + b*log(c*x^n))^3]
    @test_int [(a + b*log(c*x^n))^3/x^1, x, 2, (a + b*log(c*x^n))^4/(4*b*n)]
    @test_int [(a + b*log(c*x^n))^3/x^2, x, 3, (-6*b^3*n^3)/x - (6*b^2*n^2*(a + b*log(c*x^n)))/x - (3*b*n*(a + b*log(c*x^n))^2)/x - (a + b*log(c*x^n))^3/x]
    @test_int [(a + b*log(c*x^n))^3/x^3, x, 3, (-3*b^3*n^3)/(8*x^2) - (3*b^2*n^2*(a + b*log(c*x^n)))/(4*x^2) - (3*b*n*(a + b*log(c*x^n))^2)/(4*x^2) - (a + b*log(c*x^n))^3/(2*x^2)]
    @test_int [(a + b*log(c*x^n))^3/x^4, x, 3, (-2*b^3*n^3)/(27*x^3) - (2*b^2*n^2*(a + b*log(c*x^n)))/(9*x^3) - (b*n*(a + b*log(c*x^n))^2)/(3*x^3) - (a + b*log(c*x^n))^3/(3*x^3)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

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

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

    @test_int [x^3/(a + b*log(c*x^n))^3, x, 4, (8*x^4*ExpIntegralEi((4*(a + b*log(c*x^n)))/(b*n)))/(b^3*ℯ^((4*a)/(b*n))*n^3*(c*x^n)^(4/n)) - x^4/(2*b*n*(a + b*log(c*x^n))^2) - (2*x^4)/(b^2*n^2*(a + b*log(c*x^n)))]
    @test_int [x^2/(a + b*log(c*x^n))^3, x, 4, (9*x^3*ExpIntegralEi((3*(a + b*log(c*x^n)))/(b*n)))/(2*b^3*ℯ^((3*a)/(b*n))*n^3*(c*x^n)^(3/n)) - x^3/(2*b*n*(a + b*log(c*x^n))^2) - (3*x^3)/(2*b^2*n^2*(a + b*log(c*x^n)))]
    @test_int [x^1/(a + b*log(c*x^n))^3, x, 4, (2*x^2*ExpIntegralEi((2*(a + b*log(c*x^n)))/(b*n)))/(b^3*ℯ^((2*a)/(b*n))*n^3*(c*x^n)^(2/n)) - x^2/(2*b*n*(a + b*log(c*x^n))^2) - x^2/(b^2*n^2*(a + b*log(c*x^n)))]
    @test_int [x^0/(a + b*log(c*x^n))^3, x, 4, (x*ExpIntegralEi((a + b*log(c*x^n))/(b*n)))/(2*b^3*ℯ^(a/(b*n))*n^3*(c*x^n)^n^(-1)) - x/(2*b*n*(a + b*log(c*x^n))^2) - x/(2*b^2*n^2*(a + b*log(c*x^n)))]
    @test_int [1/(x^1*(a + b*log(c*x^n))^3), x, 2, -1/(2*b*n*(a + b*log(c*x^n))^2)]
    @test_int [1/(x^2*(a + b*log(c*x^n))^3), x, 4, (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-((a + b*log(c*x^n))/(b*n))))/(2*b^3*n^3*x) - 1/(2*b*n*x*(a + b*log(c*x^n))^2) + 1/(2*b^2*n^2*x*(a + b*log(c*x^n)))]
    @test_int [1/(x^3*(a + b*log(c*x^n))^3), x, 4, (2*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi((-2*(a + b*log(c*x^n)))/(b*n)))/(b^3*n^3*x^2) - 1/(2*b*n*x^2*(a + b*log(c*x^n))^2) + 1/(b^2*n^2*x^2*(a + b*log(c*x^n)))]
    @test_int [1/(x^4*(a + b*log(c*x^n))^3), x, 4, (9*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)*ExpIntegralEi((-3*(a + b*log(c*x^n)))/(b*n)))/(2*b^3*n^3*x^3) - 1/(2*b*n*x^3*(a + b*log(c*x^n))^2) + 3/(2*b^2*n^2*x^3*(a + b*log(c*x^n)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*x)^(m/2)*(a+b*log(c*x^n))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(d*x)^(5/2)*(a + b*log(c*x^n)), x, 1, (-4*b*n*(d*x)^(7/2))/(49*d) + (2*(d*x)^(7/2)*(a + b*log(c*x^n)))/(7*d)]
    @test_int [(d*x)^(3/2)*(a + b*log(c*x^n)), x, 1, (-4*b*n*(d*x)^(5/2))/(25*d) + (2*(d*x)^(5/2)*(a + b*log(c*x^n)))/(5*d)]
    @test_int [(d*x)^(1/2)*(a + b*log(c*x^n)), x, 1, (-4*b*n*(d*x)^(3/2))/(9*d) + (2*(d*x)^(3/2)*(a + b*log(c*x^n)))/(3*d)]
    @test_int [(a + b*log(c*x^n))/(d*x)^(1/2), x, 1, (-4*b*n*sqrt(d*x))/d + (2*sqrt(d*x)*(a + b*log(c*x^n)))/d]
    @test_int [(a + b*log(c*x^n))/(d*x)^(3/2), x, 1, (-4*b*n)/(d*sqrt(d*x)) - (2*(a + b*log(c*x^n)))/(d*sqrt(d*x))]
    @test_int [(a + b*log(c*x^n))/(d*x)^(5/2), x, 1, (-4*b*n)/(9*d*(d*x)^(3/2)) - (2*(a + b*log(c*x^n)))/(3*d*(d*x)^(3/2))]

    @test_int [(d*x)^(5/2)*(a + b*log(c*x^n))^2, x, 2, (16*b^2*n^2*(d*x)^(7/2))/(343*d) - (8*b*n*(d*x)^(7/2)*(a + b*log(c*x^n)))/(49*d) + (2*(d*x)^(7/2)*(a + b*log(c*x^n))^2)/(7*d)]
    @test_int [(d*x)^(3/2)*(a + b*log(c*x^n))^2, x, 2, (16*b^2*n^2*(d*x)^(5/2))/(125*d) - (8*b*n*(d*x)^(5/2)*(a + b*log(c*x^n)))/(25*d) + (2*(d*x)^(5/2)*(a + b*log(c*x^n))^2)/(5*d)]
    @test_int [(d*x)^(1/2)*(a + b*log(c*x^n))^2, x, 2, (16*b^2*n^2*(d*x)^(3/2))/(27*d) - (8*b*n*(d*x)^(3/2)*(a + b*log(c*x^n)))/(9*d) + (2*(d*x)^(3/2)*(a + b*log(c*x^n))^2)/(3*d)]
    @test_int [(a + b*log(c*x^n))^2/(d*x)^(1/2), x, 2, (16*b^2*n^2*sqrt(d*x))/d - (8*b*n*sqrt(d*x)*(a + b*log(c*x^n)))/d + (2*sqrt(d*x)*(a + b*log(c*x^n))^2)/d]
    @test_int [(a + b*log(c*x^n))^2/(d*x)^(3/2), x, 2, (-16*b^2*n^2)/(d*sqrt(d*x)) - (8*b*n*(a + b*log(c*x^n)))/(d*sqrt(d*x)) - (2*(a + b*log(c*x^n))^2)/(d*sqrt(d*x))]
    @test_int [(a + b*log(c*x^n))^2/(d*x)^(5/2), x, 2, (-16*b^2*n^2)/(27*d*(d*x)^(3/2)) - (8*b*n*(a + b*log(c*x^n)))/(9*d*(d*x)^(3/2)) - (2*(a + b*log(c*x^n))^2)/(3*d*(d*x)^(3/2))]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [(d*x)^(5/2)/(a + b*log(c*x^n)), x, 2, ((d*x)^(7/2)*ExpIntegralEi((7*(a + b*log(c*x^n)))/(2*b*n)))/(b*d*ℯ^((7*a)/(2*b*n))*n*(c*x^n)^(7/(2*n)))]
    @test_int [(d*x)^(3/2)/(a + b*log(c*x^n)), x, 2, ((d*x)^(5/2)*ExpIntegralEi((5*(a + b*log(c*x^n)))/(2*b*n)))/(b*d*ℯ^((5*a)/(2*b*n))*n*(c*x^n)^(5/(2*n)))]
    @test_int [(d*x)^(1/2)/(a + b*log(c*x^n)), x, 2, ((d*x)^(3/2)*ExpIntegralEi((3*(a + b*log(c*x^n)))/(2*b*n)))/(b*d*ℯ^((3*a)/(2*b*n))*n*(c*x^n)^(3/(2*n)))]
    @test_int [1/((d*x)^(1/2)*(a + b*log(c*x^n))), x, 2, (sqrt(d*x)*ExpIntegralEi((a + b*log(c*x^n))/(2*b*n)))/(b*d*ℯ^(a/(2*b*n))*n*(c*x^n)^(1/(2*n)))]
    @test_int [1/((d*x)^(3/2)*(a + b*log(c*x^n))), x, 2, (ℯ^(a/(2*b*n))*(c*x^n)^(1/(2*n))*ExpIntegralEi(-(a + b*log(c*x^n))/(2*b*n)))/(b*d*n*sqrt(d*x))]
    @test_int [1/((d*x)^(5/2)*(a + b*log(c*x^n))), x, 2, (ℯ^((3*a)/(2*b*n))*(c*x^n)^(3/(2*n))*ExpIntegralEi((-3*(a + b*log(c*x^n)))/(2*b*n)))/(b*d*n*(d*x)^(3/2))]

    @test_int [(d*x)^(5/2)/(a + b*log(c*x^n))^2, x, 3, (7*(d*x)^(7/2)*ExpIntegralEi((7*(a + b*log(c*x^n)))/(2*b*n)))/(2*b^2*d*ℯ^((7*a)/(2*b*n))*n^2*(c*x^n)^(7/(2*n))) - (d*x)^(7/2)/(b*d*n*(a + b*log(c*x^n)))]
    @test_int [(d*x)^(3/2)/(a + b*log(c*x^n))^2, x, 3, (5*(d*x)^(5/2)*ExpIntegralEi((5*(a + b*log(c*x^n)))/(2*b*n)))/(2*b^2*d*ℯ^((5*a)/(2*b*n))*n^2*(c*x^n)^(5/(2*n))) - (d*x)^(5/2)/(b*d*n*(a + b*log(c*x^n)))]
    @test_int [(d*x)^(1/2)/(a + b*log(c*x^n))^2, x, 3, (3*(d*x)^(3/2)*ExpIntegralEi((3*(a + b*log(c*x^n)))/(2*b*n)))/(2*b^2*d*ℯ^((3*a)/(2*b*n))*n^2*(c*x^n)^(3/(2*n))) - (d*x)^(3/2)/(b*d*n*(a + b*log(c*x^n)))]
    @test_int [1/((d*x)^(1/2)*(a + b*log(c*x^n))^2), x, 3, (sqrt(d*x)*ExpIntegralEi((a + b*log(c*x^n))/(2*b*n)))/(2*b^2*d*ℯ^(a/(2*b*n))*n^2*(c*x^n)^(1/(2*n))) - sqrt(d*x)/(b*d*n*(a + b*log(c*x^n)))]
    @test_int [1/((d*x)^(3/2)*(a + b*log(c*x^n))^2), x, 3, -(ℯ^(a/(2*b*n))*(c*x^n)^(1/(2*n))*ExpIntegralEi(-(a + b*log(c*x^n))/(2*b*n)))/(2*b^2*d*n^2*sqrt(d*x)) - 1/(b*d*n*sqrt(d*x)*(a + b*log(c*x^n)))]
    @test_int [1/((d*x)^(5/2)*(a + b*log(c*x^n))^2), x, 3, (-3*ℯ^((3*a)/(2*b*n))*(c*x^n)^(3/(2*n))*ExpIntegralEi((-3*(a + b*log(c*x^n)))/(2*b*n)))/(2*b^2*d*n^2*(d*x)^(3/2)) - 1/(b*d*n*(d*x)^(3/2)*(a + b*log(c*x^n)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*log(c*x^n))^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(a + b*log(c*x^n))^(1/2), x, 4, ((-(1/2))*sqrt(b)*sqrt(n)*sqrt(pi)*x*erfi(sqrt(a + b*log(c*x^n))/(sqrt(b)*sqrt(n))))/(ℯ^(a/(b*n))*(c*x^n)^n^(-1)) + x*sqrt(a + b*log(c*x^n))]

    @test_int [x^3*sqrt(log(a*x^n)), x, 4, ((-(1/16))*sqrt(n)*sqrt(pi)*x^4*erfi((2*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(4/n) + (1/4)*x^4*sqrt(log(a*x^n))]
    @test_int [x^2*sqrt(log(a*x^n)), x, 4, ((-(1/6))*sqrt(n)*sqrt(pi/3)*x^3*erfi((sqrt(3)*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(3/n) + (1/3)*x^3*sqrt(log(a*x^n))]
    @test_int [x^1*sqrt(log(a*x^n)), x, 4, ((-(1/4))*sqrt(n)*sqrt(pi/2)*x^2*erfi((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(2/n) + (1/2)*x^2*sqrt(log(a*x^n))]
    @test_int [x^0*sqrt(log(a*x^n)), x, 4, ((-(1/2))*sqrt(n)*sqrt(pi)*x*erfi(sqrt(log(a*x^n))/sqrt(n)))/(a*x^n)^n^(-1) + x*sqrt(log(a*x^n))]
    @test_int [sqrt(log(a*x^n))/x^1, x, 2, (2*log(a*x^n)^(3/2))/(3*n)]
    @test_int [sqrt(log(a*x^n))/x^2, x, 4, (sqrt(n)*sqrt(pi)*(a*x^n)^(1/n)*erf(sqrt(log(a*x^n))/sqrt(n)))/(2*x) - sqrt(log(a*x^n))/x]
    @test_int [sqrt(log(a*x^n))/x^3, x, 4, (sqrt(n)*sqrt(pi/2)*(a*x^n)^(2/n)*erf((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(4*x^2) - sqrt(log(a*x^n))/(2*x^2)]

    @test_int [x^3*log(a*x^n)^(3/2), x, 5, ((3/128)*n^(3/2)*sqrt(pi)*x^4*erfi((2*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(4/n) - (3/32)*n*x^4*sqrt(log(a*x^n)) + (1/4)*x^4*log(a*x^n)^(3/2)]
    @test_int [x^2*log(a*x^n)^(3/2), x, 5, ((1/12)*n^(3/2)*sqrt(pi/3)*x^3*erfi((sqrt(3)*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(3/n) - (1/6)*n*x^3*sqrt(log(a*x^n)) + (1/3)*x^3*log(a*x^n)^(3/2)]
    @test_int [x^1*log(a*x^n)^(3/2), x, 5, ((3/16)*n^(3/2)*sqrt(pi/2)*x^2*erfi((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(a*x^n)^(2/n) - (3/8)*n*x^2*sqrt(log(a*x^n)) + (1/2)*x^2*log(a*x^n)^(3/2)]
    @test_int [x^0*log(a*x^n)^(3/2), x, 5, ((3/4)*n^(3/2)*sqrt(pi)*x*erfi(sqrt(log(a*x^n))/sqrt(n)))/(a*x^n)^n^(-1) - (3/2)*n*x*sqrt(log(a*x^n)) + x*log(a*x^n)^(3/2)]
    @test_int [log(a*x^n)^(3/2)/x^1, x, 2, (2*log(a*x^n)^(5/2))/(5*n)]
    @test_int [log(a*x^n)^(3/2)/x^2, x, 5, (3*n^(3/2)*sqrt(pi)*(a*x^n)^(1/n)*erf(sqrt(log(a*x^n))/sqrt(n)))/(4*x) - (3*n*sqrt(log(a*x^n)))/(2*x) - log(a*x^n)^(3/2)/x]
    @test_int [log(a*x^n)^(3/2)/x^3, x, 5, (3*n^(3/2)*sqrt(pi/2)*(a*x^n)^(2/n)*erf((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(16*x^2) - (3*n*sqrt(log(a*x^n)))/(8*x^2) - log(a*x^n)^(3/2)/(2*x^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3/sqrt(log(a*x^n)), x, 3, (sqrt(pi)*x^4*erfi((2*sqrt(log(a*x^n)))/sqrt(n)))/(2*sqrt(n)*(a*x^n)^(4/n))]
    @test_int [x^2/sqrt(log(a*x^n)), x, 3, (sqrt(pi/3)*x^3*erfi((sqrt(3)*sqrt(log(a*x^n)))/sqrt(n)))/(sqrt(n)*(a*x^n)^(3/n))]
    @test_int [x^1/sqrt(log(a*x^n)), x, 3, (sqrt(pi/2)*x^2*erfi((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(sqrt(n)*(a*x^n)^(2/n))]
    @test_int [x^0/sqrt(log(a*x^n)), x, 3, (sqrt(pi)*x*erfi(sqrt(log(a*x^n))/sqrt(n)))/(sqrt(n)*(a*x^n)^n^(-1))]
    @test_int [1/(x^1*sqrt(log(a*x^n))), x, 2, (2*sqrt(log(a*x^n)))/n]
    @test_int [1/(x^2*sqrt(log(a*x^n))), x, 3, (sqrt(pi)*(a*x^n)^n^(-1)*erf(sqrt(log(a*x^n))/sqrt(n)))/(sqrt(n)*x)]
    @test_int [1/(x^3*sqrt(log(a*x^n))), x, 3, (sqrt(pi/2)*(a*x^n)^(2/n)*erf((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(sqrt(n)*x^2)]

    @test_int [x^3/log(a*x^n)^(3/2), x, 4, (4*sqrt(pi)*x^4*erfi((2*sqrt(log(a*x^n)))/sqrt(n)))/(n^(3/2)*(a*x^n)^(4/n)) - (2*x^4)/(n*sqrt(log(a*x^n)))]
    @test_int [x^2/log(a*x^n)^(3/2), x, 4, (2*sqrt(3*pi)*x^3*erfi((sqrt(3)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^(3/n)*n^(3/2)) - (2*x^3)/(n*sqrt(log(a*x^n)))]
    @test_int [x^1/log(a*x^n)^(3/2), x, 4, (2*sqrt(2*pi)*x^2*erfi((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^(2/n)*n^(3/2)) - (2*x^2)/(n*sqrt(log(a*x^n)))]
    @test_int [x^0/log(a*x^n)^(3/2), x, 4, (2*sqrt(pi)*x*erfi(sqrt(log(a*x^n))/sqrt(n)))/(n^(3/2)*(a*x^n)^n^(-1)) - (2*x)/(n*sqrt(log(a*x^n)))]
    @test_int [1/(x^1*log(a*x^n)^(3/2)), x, 2, -2/(n*sqrt(log(a*x^n)))]
    @test_int [1/(x^2*log(a*x^n)^(3/2)), x, 4, -((2*sqrt(pi)*(a*x^n)^(1/n)*erf(sqrt(log(a*x^n))/sqrt(n)))/(n^(3/2)*x)) - 2/(n*x*sqrt(log(a*x^n)))]
    @test_int [1/(x^3*log(a*x^n)^(3/2)), x, 4, -((2*sqrt(2*pi)*(a*x^n)^(2/n)*erf((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(n^(3/2)*x^2)) - 2/(n*x^2*sqrt(log(a*x^n)))]

    @test_int [x^3/log(a*x^n)^(5/2), x, 5, (32*sqrt(pi)*x^4*erfi((2*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^(4/n)*(3*n^(5/2))) - (2*x^4)/(3*n*log(a*x^n)^(3/2)) - (16*x^4)/(3*n^2*sqrt(log(a*x^n)))]
    @test_int [x^2/log(a*x^n)^(5/2), x, 5, (4*sqrt(3*pi)*x^3*erfi((sqrt(3)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^(3/n)*n^(5/2)) - (2*x^3)/(3*n*log(a*x^n)^(3/2)) - (4*x^3)/(n^2*sqrt(log(a*x^n)))]
    @test_int [x^1/log(a*x^n)^(5/2), x, 5, (8*sqrt(2*pi)*x^2*erfi((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^(2/n)*(3*n^(5/2))) - (2*x^2)/(3*n*log(a*x^n)^(3/2)) - (8*x^2)/(3*n^2*sqrt(log(a*x^n)))]
    @test_int [x^0/log(a*x^n)^(5/2), x, 5, (4*sqrt(pi)*x*erfi(sqrt(log(a*x^n))/sqrt(n)))/((a*x^n)^n^(-1)*(3*n^(5/2))) - (2*x)/(3*n*log(a*x^n)^(3/2)) - (4*x)/(3*n^2*sqrt(log(a*x^n)))]
    @test_int [1/(x^1*log(a*x^n)^(5/2)), x, 2, -(2/(3*n*log(a*x^n)^(3/2)))]
    @test_int [1/(x^2*log(a*x^n)^(5/2)), x, 5, (4*sqrt(pi)*(a*x^n)^(1/n)*erf(sqrt(log(a*x^n))/sqrt(n)))/(3*n^(5/2)*x) - 2/(3*n*x*log(a*x^n)^(3/2)) + 4/(3*n^2*x*sqrt(log(a*x^n)))]
    @test_int [1/(x^3*log(a*x^n)^(5/2)), x, 5, (8*sqrt(2*pi)*(a*x^n)^(2/n)*erf((sqrt(2)*sqrt(log(a*x^n)))/sqrt(n)))/(3*n^(5/2)*x^2) - 2/(3*n*x^2*log(a*x^n)^(3/2)) + 8/(3*n^2*x^2*sqrt(log(a*x^n)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*x)^m*(a+b*log(c*x^n))^p*with*m*symbolic=#

    @test_int [(d*x)^m*(a + a*(m + 1)/n*log(c*x^n)), x, 1, (a*(d*x)^(1 + m)*log(c*x^n))/(d*n)]

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

    @test_int [(d*x)^(n - 1)*log(c*x^n)^3, x, 3, -((6*(d*x)^n)/(d*n)) + (6*(d*x)^n*log(c*x^n))/(d*n) - (3*(d*x)^n*log(c*x^n)^2)/(d*n) + ((d*x)^n*log(c*x^n)^3)/(d*n)]
    @test_int [(d*x)^(n - 1)*log(c*x^n)^2, x, 2, (2*(d*x)^n)/(d*n) - (2*(d*x)^n*log(c*x^n))/(d*n) + ((d*x)^n*log(c*x^n)^2)/(d*n)]
    @test_int [(d*x)^(n - 1)*log(c*x^n)^1, x, 1, -((d*x)^n/(d*n)) + ((d*x)^n*log(c*x^n))/(d*n)]
    @test_int [(d*x)^(n - 1)/log(c*x^n)^1, x, 3, (x^(1 - n)*(d*x)^(-1 + n)*LogIntegral(c*x^n))/(c*n)]
    @test_int [(d*x)^(n - 1)/log(c*x^n)^2, x, 4, -((d*x)^n/(d*n*log(c*x^n))) + (x^(1 - n)*(d*x)^(-1 + n)*LogIntegral(c*x^n))/(c*n)]
    @test_int [(d*x)^(n - 1)/log(c*x^n)^3, x, 5, -((d*x)^n/(2*d*n*log(c*x^n)^2)) - (d*x)^n/(2*d*n*log(c*x^n)) + (x^(1 - n)*(d*x)^(-1 + n)*LogIntegral(c*x^n))/(2*c*n)]

    @test_int [x^m*log(a*x^n)^(3/2), x, 5, (3*n^(3/2)*sqrt(pi)*x^(1 + m)*erfi((sqrt(1 + m)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^((1 + m)/n)*(4*(1 + m)^(5/2))) - (3*n*x^(1 + m)*sqrt(log(a*x^n)))/(2*(1 + m)^2) + (x^(1 + m)*log(a*x^n)^(3/2))/(1 + m)]
    @test_int [x^m*log(a*x^n)^(1/2), x, 4, -((sqrt(n)*sqrt(pi)*x^(1 + m)*erfi((sqrt(1 + m)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^((1 + m)/n)*(2*(1 + m)^(3/2)))) + (x^(1 + m)*sqrt(log(a*x^n)))/(1 + m)]
    @test_int [x^m/log(a*x^n)^(1/2), x, 3, (sqrt(pi)*x^(1 + m)*erfi((sqrt(1 + m)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^((1 + m)/n)*(sqrt(1 + m)*sqrt(n)))]
    @test_int [x^m/log(a*x^n)^(3/2), x, 4, (2*sqrt(1 + m)*sqrt(pi)*x^(1 + m)*erfi((sqrt(1 + m)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^((1 + m)/n)*n^(3/2)) - (2*x^(1 + m))/(n*sqrt(log(a*x^n)))]
    @test_int [x^m/log(a*x^n)^(5/2), x, 5, (4*(1 + m)^(3/2)*sqrt(pi)*x^(1 + m)*erfi((sqrt(1 + m)*sqrt(log(a*x^n)))/sqrt(n)))/((a*x^n)^((1 + m)/n)*(3*n^(5/2))) - (2*x^(1 + m))/(3*n*log(a*x^n)^(3/2)) - (4*(1 + m)*x^(1 + m))/(3*n^2*sqrt(log(a*x^n)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*x)^m*(a+b*log(a*x^n))^p*with*p*symbolic=#

    @test_int [(d*x)^m*(a + b*log(c*x^n))^p, x, 2, ((d*x)^(1 + m)*Gamma(1 + p, -(((1 + m)*(a + b*log(c*x^n)))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^((a*(1 + m))/(b*n))*(c*x^n)^((1 + m)/n)*(-(((1 + m)*(a + b*log(c*x^n)))/(b*n)))^p*(d*(1 + m)))]

    @test_int [x^2*(a + b*log(c*x^n))^p, x, 2, (3^(-1 - p)*x^3*Gamma(1 + p, -((3*(a + b*log(c*x^n)))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)*(-((a + b*log(c*x^n))/(b*n)))^p)]
    @test_int [x^1*(a + b*log(c*x^n))^p, x, 2, (2^(-1 - p)*x^2*Gamma(1 + p, -((2*(a + b*log(c*x^n)))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*(-((a + b*log(c*x^n))/(b*n)))^p)]
    @test_int [x^0*(a + b*log(c*x^n))^p, x, 2, (x*Gamma(1 + p, -((a + b*log(c*x^n))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^(a/(b*n))*(c*x^n)^n^(-1)*(-((a + b*log(c*x^n))/(b*n)))^p)]
    @test_int [(a + b*log(c*x^n))^p/x^1, x, 2, (a + b*log(c*x^n))^(1 + p)/(b*n*(1 + p))]
    @test_int [(a + b*log(c*x^n))^p/x^2, x, 2, -((ℯ^(a/(b*n))*(c*x^n)^(1/n)*Gamma(1 + p, (a + b*log(c*x^n))/(b*n))*(a + b*log(c*x^n))^p)/(((a + b*log(c*x^n))/(b*n))^p*x))]
    @test_int [(a + b*log(c*x^n))^p/x^3, x, 2, -((2^(-1 - p)*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*Gamma(1 + p, (2*(a + b*log(c*x^n)))/(b*n))*(a + b*log(c*x^n))^p)/(((a + b*log(c*x^n))/(b*n))^p*x^2))]
    @test_int [(a + b*log(c*x^n))^p/x^4, x, 2, -((3^(-1 - p)*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)*Gamma(1 + p, (3*(a + b*log(c*x^n)))/(b*n))*(a + b*log(c*x^n))^p)/(((a + b*log(c*x^n))/(b*n))^p*x^3))]

    @test_int [(d*x)^m*(a + b*log(c*x))^p, x, 2, ((c*x)^(-1 - m)*(d*x)^(1 + m)*Gamma(1 + p, -(((1 + m)*(a + b*log(c*x)))/b))*(a + b*log(c*x))^p)/(ℯ^((a*(1 + m))/b)*(-(((1 + m)*(a + b*log(c*x)))/b))^p*(d*(1 + m)))]

    @test_int [x^2*(a + b*log(c*x))^p, x, 2, (3^(-1 - p)*Gamma(1 + p, -((3*(a + b*log(c*x)))/b))*(a + b*log(c*x))^p)/(ℯ^((3*a)/b)*(-((a + b*log(c*x))/b))^p*c^3)]
    @test_int [x^1*(a + b*log(c*x))^p, x, 2, (2^(-1 - p)*Gamma(1 + p, -((2*(a + b*log(c*x)))/b))*(a + b*log(c*x))^p)/(ℯ^((2*a)/b)*(-((a + b*log(c*x))/b))^p*c^2)]
    @test_int [x^0*(a + b*log(c*x))^p, x, 2, (Gamma(1 + p, -((a + b*log(c*x))/b))*(a + b*log(c*x))^p)/(ℯ^(a/b)*(-((a + b*log(c*x))/b))^p*c)]
    @test_int [(a + b*log(c*x))^p/x^1, x, 2, (a + b*log(c*x))^(1 + p)/(b*(1 + p))]
    @test_int [(a + b*log(c*x))^p/x^2, x, 2, ((-c)*ℯ^(a/b)*Gamma(1 + p, (a + b*log(c*x))/b)*(a + b*log(c*x))^p)/((a + b*log(c*x))/b)^p]
    @test_int [(a + b*log(c*x))^p/x^3, x, 2, ((-2^(-1 - p))*c^2*ℯ^((2*a)/b)*Gamma(1 + p, (2*(a + b*log(c*x)))/b)*(a + b*log(c*x))^p)/((a + b*log(c*x))/b)^p]
    @test_int [(a + b*log(c*x))^p/x^4, x, 2, ((-3^(-1 - p))*c^3*ℯ^((3*a)/b)*Gamma(1 + p, (3*(a + b*log(c*x)))/b)*(a + b*log(c*x))^p)/((a + b*log(c*x))/b)^p]

    @test_int [(d*x)^m*(a + b*log(c*sqrt(x)))^p, x, 2, ((d*x)^(1 + m)*Gamma(1 + p, -((2*(1 + m)*(a + b*log(c*sqrt(x))))/b))*(a + b*log(c*sqrt(x)))^p)/(2^p*ℯ^((2*a*(1 + m))/b)*(c*sqrt(x))^(2*(1 + m))*(-(((1 + m)*(a + b*log(c*sqrt(x))))/b))^p*(d*(1 + m)))]

    @test_int [x^2*(a + b*log(c*sqrt(x)))^p, x, 2, (3^(-1 - p)*Gamma(1 + p, -((6*(a + b*log(c*sqrt(x))))/b))*(a + b*log(c*sqrt(x)))^p)/(2^p*ℯ^((6*a)/b)*(-((a + b*log(c*sqrt(x)))/b))^p*c^6)]
    @test_int [x^1*(a + b*log(c*sqrt(x)))^p, x, 2, (2^(-1 - 2*p)*Gamma(1 + p, -((4*(a + b*log(c*sqrt(x))))/b))*(a + b*log(c*sqrt(x)))^p)/(ℯ^((4*a)/b)*(-((a + b*log(c*sqrt(x)))/b))^p*c^4)]
    @test_int [x^0*(a + b*log(c*sqrt(x)))^p, x, 2, (Gamma(1 + p, -((2*(a + b*log(c*sqrt(x))))/b))*(a + b*log(c*sqrt(x)))^p)/(2^p*ℯ^((2*a)/b)*(-((a + b*log(c*sqrt(x)))/b))^p*c^2)]
    @test_int [(a + b*log(c*sqrt(x)))^p/x^1, x, 2, (2*(a + b*log(c*sqrt(x)))^(1 + p))/(b*(1 + p))]
    @test_int [(a + b*log(c*sqrt(x)))^p/x^2, x, 2, ((-2^(-p))*c^2*ℯ^((2*a)/b)*Gamma(1 + p, (2*(a + b*log(c*sqrt(x))))/b)*(a + b*log(c*sqrt(x)))^p)/((a + b*log(c*sqrt(x)))/b)^p]
    @test_int [(a + b*log(c*sqrt(x)))^p/x^3, x, 2, ((-2^(-1 - 2*p))*c^4*ℯ^((4*a)/b)*Gamma(1 + p, (4*(a + b*log(c*sqrt(x))))/b)*(a + b*log(c*sqrt(x)))^p)/((a + b*log(c*sqrt(x)))/b)^p]
    @test_int [(a + b*log(c*sqrt(x)))^p/x^4, x, 2, ((-2^(-p))*3^(-1 - p)*c^6*ℯ^((6*a)/b)*Gamma(1 + p, (6*(a + b*log(c*sqrt(x))))/b)*(a + b*log(c*sqrt(x)))^p)/((a + b*log(c*sqrt(x)))/b)^p]

    @test_int [x^(n - 1)*(a + b*log(c*x^n))^p, x, 2, (Gamma(1 + p, -((a + b*log(c*x^n))/b))*(a + b*log(c*x^n))^p)/(ℯ^(a/b)*(-((a + b*log(c*x^n))/b))^p*(c*n))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*x^q)^m*(a+b*log(c*x^n))^p=#

    @test_int [(d*x^q)^m*(a + b*log(c*x^n))^p, x, 3, (x*(d*x^q)^m*Gamma(1 + p, -(((1 + m*q)*(a + b*log(c*x^n)))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^((a + a*m*q)/(b*n))*(c*x^n)^((1 + m*q)/n)*(-(((1 + m*q)*(a + b*log(c*x^n)))/(b*n)))^p*(1 + m*q))]

    @test_int [(d1*x^q1)^m1*(d2*x^q2)^m2*(a + b*log(c*x^n))^p, x, 4, (x*(d1*x^q1)^m1*(d2*x^q2)^m2*Gamma(1 + p, -(((1 + m1*q1 + m2*q2)*(a + b*log(c*x^n)))/(b*n)))*(a + b*log(c*x^n))^p)/(ℯ^((a*(1 + m1*q1 + m2*q2))/(b*n))*(c*x^n)^((1 + m1*q1 + m2*q2)/n)*(-(((1 + m1*q1 + m2*q2)*(a + b*log(c*x^n)))/(b*n)))^p*(1 + m1*q1 + m2*q2))]
    end

max_stars_repo_path: test/3 Logarithms/3.1.2 (d x)^m (a+b log(c x^n))^p.jl


text:
"""
Sub-module of the module `DynamicalSystemsBase`, which contains pre-defined
famous systems.
"""
module Systems
using DynamicalSystemsBase

using StaticArrays
const twopi = 2π

using DynamicalSystemsBase: DDS
include("discrete_famous_systems.jl")

using DynamicalSystemsBase: CDS
include("continuous_famous_systems.jl")

end# Systems module


text:
using BioSimulator
using Test, Statistics

for model in readdir("test-models")
  include(joinpath("test-models", model))
end

include(joinpath("interface", "well-mixed", "network.jl"))

include(joinpath("data-structures", "dep_graph.jl"))
include(joinpath("data-structures", "priority_queue.jl"))

include(joinpath("model", "reaction_system.jl"))

include(joinpath("simulators", "mean_convergence.jl"))

#### spatial simulations
include(joinpath("lattice.jl"))
include(joinpath("ips-interface.jl"))
include(joinpath("execute.jl"))

#### simulation output
include(joinpath("output", "SamplePath.jl"))


text:
#=
# 5. Generalized edge conditions

In the previous part, we demonstrated the means of annihilating the (nearly) singular behavior at edges on a discretized surface. In some cases, our desire is not to annihilate this behavior, but simply to keep it within some bounds. In the analytical treatment of potential flow problems, this objective is served by placing an inequality constraint on the edge suction parameter. That parameter is proportional to the coefficient on the bound vortex sheet strength's singularity, so in this discrete setting, in which we have extracted the singular part of $\mathfrak{f}$ in the form of $\mathfrak{f}_0$, we expect the suction parameter to be related to the value of $\tilde{\mathfrak{f}}$ at the edge. In fact, by simple comparison, it can be shown that

$\mathfrak{e}_{k}^T\tilde{\mathfrak{f}} = -\frac{2\pi c}{\Gamma_0} \sigma_{k}$

for a flat plate of length $c$, where $\sigma_{k}$ is the suction parameter at the edge corresponding to point $k$.

Let $\sigma_{k}^{\mathrm{min}}$ and $\sigma_{k}^{\mathrm{max}}$ denote the minimum and maximum tolerable values of $\sigma_{k}$ at edge $k$. We then seek to confine the suction parameter to the range $\sigma_{k}^{\mathrm{min}} \leq \sigma_{k} \leq \sigma_{k}^{\mathrm{max}}$. This generalized edge constraint is placed on the suction parameter of the intermediate sheet $\tilde{\mathfrak{f}}^*$. To avoid confusion, we will redefine the bounds based on this smooth part of the vortex sheet rather than $\sigma_{k}$ itself; for this, we define $f^{\mathrm{min}}_{k} = -\frac{2\pi c}{\Gamma_0} \sigma_{k}^{\mathrm{max}}$ and $f^{\mathrm{max}}_{k} = -\frac{2\pi c}{\Gamma_0} \sigma_{k}^{\mathrm{min}}$ if $\Gamma_0$ is positive or $f^{\mathrm{max}}_{k} = -\frac{2\pi c}{\Gamma_0} \sigma_{k}^{\mathrm{max}}$ and $f^{\mathrm{min}}_{k} = -\frac{2\pi c}{\Gamma_0} \sigma_{k}^{\mathrm{min}}$ if $\Gamma_0$ is negative. Thus, we inspect whether the value $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}}^*$ lies in the range

$f^{\mathrm{min}}_{k} \leq \mathfrak{e}_{k}^T\tilde{\mathfrak{f}}^* \leq f^{\mathrm{max}}_{k}.$

If $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}}^*$ lies within this range, then no new vortex is created near the edge (or equivalently, a new vortex of zero strength is created); if $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}}^* > f^{\mathrm{max}}_{k}$, then we create a new vortex so that $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}} = f^{\mathrm{max}}_{k}$; and if $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}}^* < f^{\mathrm{min}}_{k}$, then we do the same, but now so that $\mathfrak{e}_{k}^T\tilde{\mathfrak{f}} = f^{\mathrm{min}}_{k}$. Note that the Kutta condition simply corresponds to setting $f^{\mathrm{min}}_{k} = f^{\mathrm{max}}_{k} = 0$.

=#

# Let's illustrate this again with the example of a flat plate. To clearly show the effect of varying the suction parameter, we will run a simulation of the vortex shedding for a couple of time steps for three different suction parameter ranges. First we create the grid and the flate plate.
#md # ```@setup 5.-Generalized-edge-conditions
#md # using GridPotentialFlow
#md # using Plots
#md # ```
#!md using GridPotentialFlow
#!md using Plots
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);
c = Lx/2 # chord length
α = -π/3 # angle of attack
plate = Plate(c,4*cellsize(g));
Tr = RigidTransform((0.0,0.0),α)
Tr(plate)
Δs = dlengthmid(plate);

# We create three instances of `PotentialFlowBody`, each one with a different suction parameter range. We only vary the suction parameter at the leading edge and keep it at zero (Kutta condition) at the trailing edge.
σLE_list = [0.0,0.05,0.1];
pfb_list = [PotentialFlowBody(plate,edges=[1,length(plate)],σ=[SuctionParameter(σLE),SuctionParameter(0.0)]) for σLE in σLE_list];

# The initial point vortices are the same in each case.
Δt = 2e-2
vLE = Vortex(plate.x[1]+3Δt*plate.len*cos(plate.α+π/2),plate.y[1]+3Δt*plate.len*sin(plate.α+π/2),0.0);
vTE = Vortex(plate.x[end]+3Δt*cos(α+π/2),plate.y[end]+3Δt*sin(α+π/2),0.0);

# At each time step, we will need to insert new point vortices near the regularized edges to simulate the vortex shedding. A good strategy is to place them one-third of the way from the edge to the last released vortex from that edge. This will happen at each time step, so we create a function for it.
function createsheddedvortices(plate,oldvortices)

    vLE = Vortex(2/3*plate.x[1]+1/3*oldvortices[end-1].x,2/3*plate.y[1]+1/3*oldvortices[end-1].y,0.0)
    vTE = Vortex(2/3*plate.x[end]+1/3*oldvortices[end].x,2/3*plate.y[end]+1/3*oldvortices[end].y,0.0)

    return vLE, vTE
end

# We create a vortex models with the initial point vortices.
model_list = [VortexModel(g,bodies=[pfb],vortices=[vLE,vTE],U∞=(1.0,0.0)) for pfb in pfb_list];

# Now we can advance the solution in time.
T = 0:Δt:0.2
for t in T
    for vm in model_list
        X = getvortexpositions(vm)
        Ẋ = vortexvelocities!(vm)
        X .= X .+ Ẋ*Δt
        setvortexpositions!(vm, X)
        vLEnew, vTEnew = createsheddedvortices(plate,vm.vortices[end-1:end])
        pushvortices!(vm,vLEnew,vTEnew)
    end
end

# By plotting the positions of the point vortices emanating from the leading edge, we see that as $\sigma_{k}^{\mathrm{max}}$ increases, the stream of point vortices is swept back from the edge.
colors = [:red,:blue,:green];
plot(plate,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2,xlabel="x",ylabel="y")
for i in 1:length(model_list)
    plot!(model_list[i].vortices.x[4:2:end],model_list[i].vortices.y[4:2:end],color=colors[i],marker=:circle,markersize=2)
    plot!(model_list[i].vortices.x[3:2:end],model_list[i].vortices.y[3:2:end],color=colors[i],marker=:circle,markersize=2)
    scatter!(model_list[i].vortices.x[1:2],model_list[i].vortices.y[1:2],color=colors[i],marker=:circle,markersize=2,label="σLE=$(σLE_list[i])")
end
plot!()

#
plot(xlabel="body point index",ylabel="f̃")
for i in 1:length(model_list)
    sol = solve(model_list[i]);
    plot!(sol.f./model_list[i].system.f₀,linecolor=colors[i],label="σLE=$(σLE_list[i])")
end
plot!()


text:
# Code in this file inspired by NetworkX.

"""
    cycle_basis(g, root=nothing)

Return a list of cycles which form a basis for cycles of graph `g`, optionally starting at node `root`.

A basis for cycles of a network is a minimal collection of
cycles such that any cycle in the network can be written
as a sum of cycles in the basis.  Here summation of cycles
is defined as "exclusive or" of the edges. Cycle bases are
useful, e.g. when deriving equations for electric circuits
using Kirchhoff's Laws.

Example:
```jldoctest
julia> elist = [(1,2),(2,3),(2,4),(3,4),(4,1),(1,5)]
julia> g = SimpleGraph(SimpleEdge.(elist))
julia> cycle_basis(g)
2-element Array{Array{Int64,1},1}:
 [2, 3, 4]
 [2, 1, 3]
```

### References
* Paton, K. An algorithm for finding a fundamental set of cycles of a graph. Comm. ACM 12, 9 (Sept 1969), 514-518. [https://dl.acm.org/citation.cfm?id=363232]
"""
function cycle_basis(g::AbstractSimpleGraph, root=nothing)
    gnodes = Set(vertices(g))
    cycles = Vector{Vector{eltype(g)}}()
    while !isempty(gnodes)
        if root == nothing
            root = pop!(gnodes)
        end
        stack = [root]
        pred = Dict(root => root)
        keys_pred = Set(root)
        used = Dict(root => [])
        keys_used = Set(root)
        while !isempty(stack)
            z = pop!(stack)
            zused = used[z]
            for nbr in neighbors(g,z)
                if !in(nbr, keys_used)
                    pred[nbr] = z
                    push!(keys_pred, nbr)
                    push!(stack,nbr)
                    used[nbr] = [z]
                    push!(keys_used, nbr)
                elseif nbr == z
                    push!(cycles, [z])
                elseif !in(nbr, zused)
                    pn = used[nbr]
                    cycle = [nbr,z]
                    p = pred[z]
                    while !in(p, pn)
                        push!(cycle, p)
                        p = pred[p]
                    end
                    push!(cycle,p)
                    push!(cycles,cycle)
                    push!(used[nbr], z)
                end
            end
        end
        setdiff!(gnodes,keys_pred)
        root = nothing
    end
    return cycles
end


text:
using DataStructures

mutable struct NonTerminalInstance
    node::BNFNode
    start_index::Int
    end_index::Union{Nothing, Int}
    contents::Vector

    function NonTerminalInstance(node, start=0)
        new(node, start, nothing, [])
    end
end

# Reduce: When we reach the end of a NonTerminalInstance we set it's
# end index, pop it from the Parser's stack, and add it to the
# contents of the new top of stack.

struct Parser
    input::String
    input_index::Int
    stack::Stack{NonTerminalInstance}

    function Parser(input::String; start=0)
        new(input, start, Stack())
    end
end

function start(parser::Parser, target::BNFNode)
    push(stack, NonTerminalInstance(target, parser.input_index))
    parser
end

function finish(parser::Parser, at::Int)
    nt = pop(parser.stack)
    nt.end_index = at
    push!(first(parser.stack)contents, nt)
    # *** How do we validate whether nt.node is a v alid consitituent
    # *** of first(parser.stack).node?
    parser
end

function parse(parser::Parser)
end


text:
"""

    iterateCSteps(setting, subsetindices, h)

Perform a concentration step for a given subset of a regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `subsetindices::Array{Int, 1}`: Indicies of observations in the initial subset.
- `h::Int`: A constant at least half of the number of observations.

# Notes
    This function is a helper for the lts function. A concentration step starts with a
    initial subset. The size of the subset is enlarged to h, a constant at least half of the
    observations. Please refer to the citation given below.

# References
Rousseeuw, Peter J., and Katrien Van Driessen. "An algorithm for positive-breakdown
regression based on concentration steps." Data Analysis.
Springer, Berlin, Heidelberg, 2000. 335-346.
"""
function iterateCSteps(setting::RegressionSetting, subsetindices::Array{Int,1}, h::Int)
    X, y = @extractRegressionSetting setting
    return iterateCSteps(X, y, subsetindices, h)
end

function iterateCSteps(X::Array{Float64,2}, y::Array{Float64,1}, subsetindices::Array{Int,1}, h::Int)
    starterset = copy(subsetindices)
    oldobjective = Inf
    objective = Inf
    iter = 0
    maxiter = 10000
    n, p = size(X)
    while iter < maxiter
        try
            Xsub = X[subsetindices, :]
            ysub = y[subsetindices]
            olsreg = ols(Xsub, ysub)
            betas = coef(olsreg)
            res = [y[i] - sum(X[i,:] .* betas) for i in 1:n]
            sortedresindices = sortperm(abs.(res))
            subsetindices = sortedresindices[1:h]
            objective = sum(sort(res.^2.0)[1:h])
            if oldobjective == objective
                break
            end
            oldobjective = objective
            iter += 1
        catch er
            @warn er
            return (objective, subsetindices)
        end
    end
    if iter >= maxiter
        @warn "in c-step stage of LTS, a h-subset is not converged for starting indices " starterset
    end
    return (objective, subsetindices)
end

function iterateCSteps(setting::RegressionSetting, initialBetas::Array{Float64,1}, h::Int)
    X = designMatrix(setting)
    y = responseVector(setting)
    return iterateCSteps(X, y, initialBetas, h)
end

function iterateCSteps(X::Array{Float64,2}, y::Array{Float64,1}, initialBetas::Array{Float64,1}, h::Int)
    n, p = size(X)
    res = [y[i] - sum(X[i,:] .* initialBetas) for i in 1:n]
    sortedresindices = sortperm(abs.(res))
    subsetindices = sortedresindices[1:p]
    return iterateCSteps(X, y, subsetindices, h)
end

"""

    lts(setting; iters, crit)

Perform the Fast-LTS (Least Trimmed Squares) algorithm for a given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `iters::Int`: Number of iterations.
- `crit::Float64`: Critical value.

# Description
The algorithm searches for estimations of regression parameters which minimize the sum of first h
ordered squared residuals where h is Int(floor((n + p + 1.0) / 2.0)). Specifically, our implementation,
uses the algorithm Fast-LTS in which concentration steps are used for enlarging a basic
subset to subset of clean observation of size h.

# Output
- `["betas"]`: Estimated regression coefficients
- `["S"]`: Standard error of regression
- `["hsubset"]`: Best subset of clean observation of size h.
- `["outliers"]`: Array of indices of outliers
- `["scaled.residuals"]`: Array of scaled residuals
- `["objective"]`: LTS objective value.

# Examples
```julia-repl
julia> reg = createRegressionSetting(@formula(calls ~ year), phones);
julia> lts(reg)
Dict{Any,Any} with 6 entries:
  "betas"            => [-56.5219, 1.16488]
  "S"                => 1.10918
  "hsubset"          => [11, 10, 5, 6, 23, 12, 13, 9, 24, 7, 3, 4, 8]
  "outliers"         => [14, 15, 16, 17, 18, 19, 20, 21]
  "scaled.residuals" => [2.41447, 1.63472, 0.584504, 0.61617, 0.197052, -0.222066, -0.551027, -0.970146, -0.397538, -0.185558  …  …
  "objective"        => 3.43133
```

# References
Rousseeuw, Peter J., and Katrien Van Driessen. "An algorithm for positive-breakdown
regression based on concentration steps." Data Analysis.
Springer, Berlin, Heidelberg, 2000. 335-346.
"""
function lts(setting::RegressionSetting; iters=nothing, crit=2.5)
    X = designMatrix(setting)
    y = responseVector(setting)
    return lts(X, y, iters=iters, crit=crit)
end

function lts(X::Array{Float64,2}, y::Array{Float64,1}; iters=nothing, crit=2.5)
    n, p = size(X)
    h = Int(floor((n + p + 1.0) / 2.0))
    if iters === nothing
        iters = minimum([500 * p, 3000])
    end
    allindices = collect(1:n)
    bestobjective = Inf
    besthsubset = []
    for iter in 1:iters
        subsetindices = sample(allindices, p, replace=false)
        objective, hsubsetindices = iterateCSteps(X, y, subsetindices, h)
        if objective < bestobjective
            bestobjective = objective
            besthsubset = hsubsetindices
        end
    end
    ltsreg = ols(X[besthsubset, :], y[besthsubset])
    ltsbetas = coef(ltsreg)
    ltsres = [y[i] - sum(X[i,:] .* ltsbetas) for i in 1:n]
    ltsS = sqrt(sum((ltsres.^2.0)[1:h]) / (h - p))
    ltsresmean = mean(ltsres[besthsubset])
    ltsScaledRes = (ltsres .- ltsresmean) / ltsS
    outlierindices = filter(i -> abs(ltsScaledRes[i]) > crit, 1:n)
    result = Dict()
    result["objective"] = bestobjective
    result["hsubset"] = besthsubset
    result["betas"] = ltsbetas
    result["S"] = ltsS
    result["outliers"] = outlierindices
    result["scaled.residuals"] = ltsScaledRes
    return result
end

max_stars_repo_name: jbytecode/LinRegOutliers

max_issues_repo_name: jbytecode/LinRegOutliers

max_forks_repo_name: jbytecode/LinRegOutliers


text:
# Deligne tensor product of different sectors: ⊠
#==============================================================================#
const SectorTuple = Tuple{Vararg{Sector}}

struct ProductSector{T<:SectorTuple} <: Sector
    sectors::T
end
_sectors(::Type{Tuple{}}) = ()
_sectors(::Type{T}) where {T<:SectorTuple} =
    (Base.tuple_type_head(T),  _sectors(Base.tuple_type_tail(T))...)

Base.IteratorSize(::Type{SectorValues{ProductSector{T}}}) where {T<:SectorTuple} =
    Base.IteratorSize(Base.Iterators.product(map(values, _sectors(T))...))
Base.size(::SectorValues{ProductSector{T}}) where {T<:SectorTuple} =
    map(s->length(values(s)), _sectors(T))
Base.length(P::SectorValues{<:ProductSector}) = *(size(P)...)

function Base.iterate(::SectorValues{ProductSector{T}}, args...) where {T<:SectorTuple}
    next = iterate(product(values.(_sectors(T))...), args...)
    next === nothing && return nothing
    val, state = next
    return ProductSector{T}(val), state
end
function Base.getindex(P::SectorValues{ProductSector{T}}, i::Int) where {T<:SectorTuple}
    Base.IteratorSize(P) isa IsInfinite &&
        throw(ArgumentError("cannot index into infinite product sector"))
    ProductSector{T}(getindex.(values.(_sectors(T)), Tuple(CartesianIndices(size(P))[i])))
end
function findindex(P::SectorValues{ProductSector{T}}, c::ProductSector{T}) where
                                                                    {T<:SectorTuple}
    Base.IteratorSize(P) isa IsInfinite &&
        throw(ArgumentError("cannot index into infinite product sector"))
    LinearIndices(size(P))[CartesianIndex(findindex.(values.(_sectors(T)), c.sectors))]
end

ProductSector{T}(args...) where {T<:SectorTuple} = ProductSector{T}(args)
Base.convert(::Type{ProductSector{T}}, t::Tuple) where {T<:SectorTuple} =
    ProductSector{T}(convert(T, t))

Base.one(::Type{ProductSector{T}}) where {I<:Sector, T<:Tuple{I}} = ProductSector((one(I),))
Base.one(::Type{ProductSector{T}}) where {I<:Sector, T<:Tuple{I, Vararg{Sector}}} =
    one(I) ⊠ one(ProductSector{Base.tuple_type_tail(T)})

Base.conj(p::ProductSector) = ProductSector(map(conj, p.sectors))
function ⊗(p1::P, p2::P) where {P<:ProductSector}
    if FusionStyle(P) isa UniqueFusion
        (P(first(product(map(⊗, p1.sectors, p2.sectors)...))),)
    else
        return SectorSet{P}(product(map(⊗, p1.sectors, p2.sectors)...))
    end
end

Nsymbol(a::P, b::P, c::P) where {P<:ProductSector} =
    prod(map(Nsymbol, a.sectors, b.sectors, c.sectors))

_firstsector(x::ProductSector) = x.sectors[1]
_tailsector(x::ProductSector) = ProductSector(tail(x.sectors))

function Fsymbol(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProductSector}
    heads = map(_firstsector, (a, b, c, d, e, f))
    tails = map(_tailsector, (a, b, c, d, e, f))
    f1 = Fsymbol(heads...)
    f2 = Fsymbol(tails...)
    if f1 isa Number || f2 isa Number
        f1 * f2
    else
        _kron(f1, f2)
    end
end
Fsymbol(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProductSector{<:Tuple{Sector}}} =
    Fsymbol(map(_firstsector, (a, b, c, d, e, f))...)

function Rsymbol(a::P, b::P, c::P) where {P<:ProductSector}
    heads = map(_firstsector, (a, b, c))
    tails = map(_tailsector, (a, b, c))
    r1 = Rsymbol(heads...)
    r2 = Rsymbol(tails...)
    if r1 isa Number || r2 isa Number
        r1 * r2
    else
        _kron(r1, r2)
    end
end
Rsymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} =
    Rsymbol(map(_firstsector, (a, b, c))...)

function Bsymbol(a::P, b::P, c::P) where {P<:ProductSector}
    heads = map(_firstsector, (a, b, c))
    tails = map(_tailsector, (a, b, c))
    b1 = Bsymbol(heads...)
    b2 = Bsymbol(tails...)
    if b1 isa Number || b2 isa Number
        b1 * b2
    else
        _kron(b1, b2)
    end
end
Bsymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} =
    Bsymbol(map(_firstsector, (a, b, c))...)

function Asymbol(a::P, b::P, c::P) where {P<:ProductSector}
    heads = map(_firstsector, (a, b, c))
    tails = map(_tailsector, (a, b, c))
    a1 = Asymbol(heads...)
    a2 = Asymbol(tails...)
    if a1 isa Number || a2 isa Number
        a1 * a2
    else
        _kron(a1, a2)
    end
end
Asymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} =
    Asymbol(map(_firstsector, (a, b, c))...)

frobeniusschur(p::ProductSector) = prod(map(frobeniusschur, p.sectors))

function fusiontensor(a::P, b::P, c::P) where {P<:ProductSector}
    return _kron(fusiontensor(map(_firstsector, (a,b,c))...),
                    fusiontensor(map(_tailsector, (a,b,c))...))
end

fusiontensor(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} =
    fusiontensor(map(_firstsector, (a,b,c))...)

FusionStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} =
    Base.:&(map(FusionStyle, _sectors(T))...)
BraidingStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} =
    Base.:&(map(BraidingStyle, _sectors(T))...)
Base.isreal(::Type{<:ProductSector{T}}) where {T<:SectorTuple} = _isreal(T)
_isreal(::Type{Tuple{}}) = true
_isreal(T::Type{<:SectorTuple}) =
    isreal(Base.tuple_type_head(T)) && _isreal(Base.tuple_type_tail(T))

fermionparity(P::ProductSector) = mapreduce(fermionparity, xor, P.sectors)

dim(p::ProductSector) = *(dim.(p.sectors)...)

Base.isequal(p1::ProductSector, p2::ProductSector) = isequal(p1.sectors, p2.sectors)
Base.hash(p::ProductSector, h::UInt) = hash(p.sectors, h)
Base.isless(p1::ProductSector{T}, p2::ProductSector{T}) where {T} =
    isless(reverse(p1.sectors), reverse(p2.sectors))

# Default construction from tensor product of sectors
#-----------------------------------------------------
⊠(s1, s2, s3, s4...) = ⊠(⊠(s1, s2), s3, s4...)
const deligneproduct = ⊠

"""
    ⊠(s₁::Sector, s₂::Sector)
    deligneproduct(s₁::Sector, s₂::Sector)

Given two sectors `s₁` and `s₂`, which label an isomorphism class of simple objects in a
fusion category ``C₁`` and ``C₂``, `s1 ⊠ s2` (obtained as `\boxtimes+TAB`) labels the
isomorphism class of simple objects in the Deligne tensor product category ``C₁ ⊠ C₂``.

The Deligne tensor product also works in the type domain and for spaces and tensors. For
group representations, we have `Irrep[G₁] ⊠ Irrep[G₂] == Irrep[G₁ × G₂]`.
"""
⊠(s1::Sector, s2::Sector) = ProductSector((s1, s2))
⊠(s1::Trivial, s2::Trivial) = s1
⊠(s1::Sector, s2::Trivial) = s1
⊠(s1::Trivial, s2::Sector) = s2
⊠(p1::ProductSector, s2::Trivial) = p1
⊠(p1::ProductSector, s2::Sector) = ProductSector(tuple(p1.sectors..., s2))
⊠(s1::Trivial, p2::ProductSector) = p2
⊠(s1::Sector, p2::ProductSector) = ProductSector(tuple(s1, p2.sectors...))
⊠(p1::ProductSector, p2::ProductSector) =
    ProductSector(tuple(p1.sectors..., p2.sectors...))

# grow types from the left using Base.tuple_type_cons
⊠(I1::Type{Trivial}, I2::Type{Trivial}) = Trivial
⊠(I1::Type{Trivial}, I2::Type{<:Sector}) = I2
⊠(I1::Type{<:Sector}, I2::Type{<:Trivial}) = I1
⊠(I1::Type{<:Sector}, I2::Type{<:Sector}) = ProductSector{Tuple{I1, I2}}

⊠(I1::Type{<:ProductSector}, I2::Type{Trivial}) = I1
⊠(I1::Type{<:ProductSector}, I2::Type{<:Sector}) = I1 ⊠ ProductSector{Tuple{I2}}

⊠(::Type{Trivial}, P::Type{ProductSector{T}}) where {T<:SectorTuple} = P
⊠(I::Type{<:Sector}, ::Type{ProductSector{T}}) where {T<:SectorTuple} =
    ProductSector{Base.tuple_type_cons(I, T)}

⊠(::Type{ProductSector{Tuple{I}}},
    ::Type{ProductSector{T}}) where {I<:Sector, T<:SectorTuple} =
    ProductSector{Base.tuple_type_cons(I, T)}

⊠(::Type{ProductSector{T1}},
        I2::Type{ProductSector{T2}}) where {T1<:SectorTuple, T2<:SectorTuple} =
    Base.tuple_type_head(T1) ⊠ (ProductSector{Base.tuple_type_tail(T1)} ⊠ I2)

function Base.show(io::IO, P::ProductSector)
    sectors = P.sectors
    compact = get(io, :typeinfo, nothing) === typeof(P)
    sep = compact ? ", " : " ⊠ "
    print(io, "(")
    for i = 1:length(sectors)
        i == 1 || print(io, sep)
        io2 = compact ? IOContext(io, :typeinfo => typeof(sectors[i])) : io
        print(io2, sectors[i])
    end
    print(io, ")")
end

function type_repr(P::Type{<:ProductSector})
    sectors = P.parameters[1].parameters
    if length(sectors) == 1
        s = "ProductSector{Tuple{" * type_repr(sectors[1]) * "}}"
    else
        s = "("
        for i in 1:length(sectors)
            if i != 1
                s *= " ⊠ "
            end
            s *= type_repr(sectors[i])
        end
        s *= ")"
    end
    return s
end

# TODO: Do we want custom printing for product of Irreps
function type_repr(::Type{ProductSector{T}}) where {T<:Tuple{Vararg{AbstractIrrep}}}
    sectors = T.parameters
    s = "Irrep["
    for i in 1:length(sectors)
        if i != 1
            s *= " × "
        end
        s *= type_repr(supertype(sectors[i]).parameters[1])
    end
    s *= "]"
    return s
end

function Base.getindex(::IrrepTable, ::Type{ProductGroup{Gs}}) where {Gs<:GroupTuple}
    G1 = tuple_type_head(Gs)
    Grem = tuple_type_tail(Gs)
    ProductSector{Tuple{Irrep[G1]}} ⊠ Irrep[ProductGroup{tuple_type_tail(Gs)}]
end
function Base.getindex(::IrrepTable, ::Type{ProductGroup{Tuple{G}}}) where {G<:Group}
    ProductSector{Tuple{Irrep[G]}}
end


text:
using CGP
using Logging
using ArgParse

CGP.Config.init("cfg/test.yaml")

function xorfit(c::Chromosome, nbits::Int64, nsamples::Int64)
    accuracy = 0
    for d in 1:nsamples
        ins = bitrand(nbits)
        outputs = process(c, Float64.(ins))
        if outputs[1] == xor(ins...)
            accuracy += 1
        end
    end
    accuracy /= nsamples
    accuracy
end

function get_args()
    argtable = ArgParseSettings()
    @add_arg_table(
        argtable,
        "--seed", arg_type = Int, default = 0,
        "--log", arg_type = String, default = "xor.log",
        "--nbits", arg_type = Int, default = 2,
        "--nsamples", arg_type = Int, default = 10
    )
    parse_args(argtable)
end

function run_xor(args::Dict)
    srand(args["seed"])
    Logging.configure(filename="xor.log", level=INFO)
    nin = args["nbits"]
    nout = 1
    fitness = x->xorfit(x, args["nbits"], args["nsamples"])
    maxfit, best = oneplus(nin, nout, fitness; id="xor")
end

if ~isinteractive()
    run_xor(get_args())
end


text:
################################################################################################
# OpenGL interop support for OpenCV
#
# http://docs.opencv.org/ref/master/df/d60/namespacecv_1_1ogl.html
# "opengl.hpp"
################################################################################################

# namespace cv::ogl

# Default constructor
oglArrays() = @cxx cv::ogl::Arrays()

# Binds all vertex arrays
oglBind(oglArr) = @cxx oglArr->bind()

# Release all inner buffers
oglRelease(oglArr) = @cxx oglArr->release()

# Resets vertex colors
oglResetColor(oglArr) = @cxx oglArr->resetColorArray()

# Resets vertex normals
oglResetNormal(oglArr) = @cxx oglArr->resetNormalArray()

# Resets vertex texture coordinates
oglResetTexCoord(oglArr) = @cxx oglArr->resetTexCoordArray()

# Resets vertex coordinates
oglResetVertex(oglArr) = @cxx oglArr->resetVertexArray()

# Sets auto release mode all inner buffers
oglSetAutoRelease(oglArr, flag::Bool) = @cxx oglArr->setAutoRelease(flag)

# Sets an array of vertex colors
oglSetColor(oglArr, color) = @cxx oglArr->setColorArray(color)
# color::InputArray, can be both host and device memory

# Sets an array of vertex normals
oglSetNormal(oglArr, normal) = @cxx oglArr->setNormalArray(normal)
# normal array with vertex normals, can be both host and device memory

# Sets an array of vertex texture coordinates
oglSetTexCoord(oglArr, texCoord) = @cxx oglArr->setTexCoordArray(texCoord)
# texCoord	array with vertex texture coordinates, can be both host and device memory

# Sets an array of vertex coordinates
oglSetVertexArray(oglArr, vertex) = @cxx oglArr->setVertexArray(vertex)
# vertex array with vertex coordinates, can be both host and device memory.

# Returns the vertex count
oglSize(oglArr) = @cxx oglArr->size()

# OpenGL Buffers
# Buffer Objects are OpenGL objects that store an array of unformatted memory allocated by the OpenGL context.
# These can be used to store vertex data, pixel data retrieved from images or the framebuffer, and a variety of other things.
# ogl::Buffer has interface similar with Mat interface and represents 2D array memory.
# ogl::Buffer supports memory transfers between host and device and also can be mapped to CUDA memory.

# basic constructor
# Creates empty ogl::Buffer object, creates ogl::Buffer object from existed buffer ( abufId parameter),
# allocates memory for ogl::Buffer object or copies from host/device memory
oglBuffer() = @cxx cv::ogl::Buffer()

# OpenGL buffer constructors
# Parameters
  # arows	Number of rows in a 2D array
  # acols	Number of columns in a 2D array
  # atype	Array type ( CV_8UC1, ..., CV_64FC4 )
  # abufI:: Buffer object name
     #  READ_ONLY, WRITE_ONLY, READ_WRITE
  # autoRelease	Auto release mode (if true, release will be called in object's destructor)
oglBuffer(arows::Int, acols::Int, atype::Int, abufId=READ_WRITE, autoRelease = false) = @cxx cv::ogl::Buffer::Buffer(arows, acols, atype, abufId, autoRelease)

# using size::cvSize (2D array size)
oglBuffer(size, atype::Int, abufId=READ_WRITE, autoRelease = false) = @cxx cv::ogl::Buffer::Buffer(size, atype, abufId, autoRelease)

# using defined buffer usage, TARGET= ARRAY_BUFFER
# ARRAY_BUFFER
# The buffer will be used as a source for vertex data.
# ELEMENT_ARRAY_BUFFER
# The buffer will be used for indices (in glDrawElements, for example)
# PIXEL_PACK_BUFFER
# The buffer will be used for reading from OpenGL textures.
# PIXEL_UNPACK_BUFFER
# The buffer will be used for writing to OpenGL textures.

oglBuffer(arows::Int, acols::Int, atype::Int, target = ARRAY_BUFFER, autoRelease = false) = @cxx cv::ogl::Buffer::Buffer(arows, acols, atype, target, autoRelease)

# takes input array
oglBuffer(arr, atype::Int, target = ARRAY_BUFFER, autoRelease = false) = @cxx cv::ogl::Buffer::Buffer(arr, target, autoRelease)

# TO DO: Texture2D
oglTexture2D() = @cxx cv::ogl::Texture2D()


text:
function DerivativeMatrixSingle(OrdPoly)
(wmat,xw)=GaussLobattoQuad(OrdPoly);
DS=zeros(OrdPoly+1,OrdPoly+1);
for i=1:OrdPoly+1
  for j=1:OrdPoly+1
    DS[i,j]=DLagrange(xw[i],xw,j);
  end
end
w = vec(wmat)
DW=-inv(diagm(w))*DS'*diagm(w);
DV=2*DS;
DV[1,1]=2*DS[1,1]+1/w[1];
DV[OrdPoly+1,OrdPoly+1]=2*DS[OrdPoly+1,OrdPoly+1]-1/w[OrdPoly+1];
DVT=DV';
B=zeros(OrdPoly+1,OrdPoly+1);
B[1,1]=-1/w[1];
B[OrdPoly+1,OrdPoly+1]=1/w[OrdPoly+1];
return (DW,DS,DV,DVT,B)
end


text:
using DataFrames
using CSV
df = DataFrame(CSV.File("../../datasets/abalone.csv", header=0))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
struct Item
    id::String
    label::Union{Nothing,String}
end

struct Menu
    header::String
    items::Array{Union{Nothing,Item},1}
end

struct MyStruct
    menu::Menu
end


text:
@testset "Metric" begin
    @testset "MeterProvider" begin
        global_meter_provider!(global_meter_provider())  # smoke test

        m = Meter("test")
        c = Counter{UInt}("cc", m)
        @test provider(c) == global_meter_provider()
        @test resource(c) == resource(global_meter_provider())
    end

    @testset "instruments" begin
        m1 = Meter("test")

        counter = Counter{Int}("counter", m1)
        counter()
        counter(1)
        counter(3; a = 1, m = "b", c = 3.0)
        @test_throws ArgumentError counter(-1)

        hist = Histogram{Float64}("hist", m1)
        hist(1)
        hist(3; a = 1, b = "b", c = 3.0)

        up_down_counter = UpDownCounter{Int}("up_down_counter", m1)
        up_down_counter(1)
        up_down_counter(3; a = 1, b = "b", c = 3.0)

        obs_counter = ObservableCounter{Int}("obs_counter", m1) do
            rand(1:10)
        end
        obs_counter()

        broken_obs_counter = ObservableCounter{Int}("obs_counter", m1) do
            -1
        end

        @test_throws ArgumentError broken_obs_counter()

        obs_gauge = ObservableGauge{Int}("obs_gauge", m1) do
            rand(1:10)
        end
        obs_gauge()

        obs_up_down_counter = ObservableUpDownCounter{Int}("obs_up_down_counter", m1) do
            rand(1:10)
        end
        obs_up_down_counter()
    end
end