text:
using QCDNUM
using Test
using Documenter

@testset "Docs" begin
    doctest(QCDNUM)
end

@testset "Initialisation" begin

    # Standard initialisation with banner
    QCDNUM.qcinit(6, "")

    # Check next logical unit number
    lun = QCDNUM.nxtlun(1)
    @test typeof(lun) == Int32

    # Store and retrieve
    QCDNUM.qstore("write", 1, 999.9)
    out = QCDNUM.qstore("read", 1)
    @test out == 999.9

    # Set QCDNUM params
    #options = ["iter", "tlmc", "nopt", "edbg"]
    #for opt in options
    #    @test QCDNUM.setint(opt, 1) == nothing
    #end

end

@testset "Grids" begin

    # Initialisation
    QCDNUM.qcinit(-6, "")

    # x grid
    nx = QCDNUM.gxmake(Float64.([1.0e-4]), Int32.([1]), 1, 100, 3)
    @test nx == 100

    # qq grid
    nq = QCDNUM.gqmake(Float64.([2e0, 1e4]), Float64.([1e0, 1e0]), 2, 50)
    @test nq == 50

    # Indices <-> grid points
    ix = QCDNUM.ixfrmx(0.1)
    @test ix == 76

    x = QCDNUM.xfrmix(ix)
    @test x ≈ 0.1

    iq = QCDNUM.iqfrmq(1000.0)
    @test iq == 36

    q = QCDNUM.qfrmiq(iq)
    @test q == 877.3066621237417

    # Grid params
    out = QCDNUM.grpars()
    @test out[1] == 100
    @test out[2] ≈ 1e-4
    @test out[3] == 1
    @test out[4] == 50
    @test out[5] == 2.0
    @test out[6] ≈ 1e4
    @test out[7] == 3

    # Copying grids
    x_grid = QCDNUM.gxcopy(nx)
    @test x_grid[1] ≈ 1e-4

    qq_grid = QCDNUM.gqcopy(nq)
    @test qq_grid[1] ≈ 2e0

    # Weights
    for itype in [1, 2, 3]
        nw = QCDNUM.fillwt(itype)
        @test typeof(nw) == Int32
    end

end

@testset "Parameters" begin

    # Initialisation
    QCDNUM.qcinit(-6, "")

    # Order
    for iord in [1, 2, 3]
        @test QCDNUM.setord(iord) == nothing
    end

    @test QCDNUM.setalf(0.364, 2.0) == nothing

    @test QCDNUM.setabr(1.0, 0.0) == nothing

end

# Define a test PDF function and mapping
function func(ipdf, x)::Float64

    i = ipdf[]
    xb = x[]
    adbar = 0.1939875
    f = 0
    if (i == 0)
        ag = 1.7
        f = ag * xb^-0.1 * (1.0-xb)^5.0
    end
    if (i == 1)
        ad = 3.064320
        f = ad * xb^0.8 * (1.0-xb)^4.0
    end
    if (i == 2)
        au = 5.107200
        f = au * xb^0.8 * (1.0-xb)^3.0
    end
    if (i == 3)
        f = 0.0
    end
    if (i == 4)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0
    end
    if (i == 5)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
    end
    if (i == 6)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i >= 7)
        f = 0.0
    end

    return f
end

function func_sgns(ipdf, x)::Float64
    i = ipdf[]
    xb = x[]

    f = 0.0

    if (i == -1)
        iset = Int(QCDNUM.qstore("read", 1))
        iq0 = Int(QCDNUM.qstore("read", 2))
    else
        iset = Int(QCDNUM.qstore("read", 1))
        iq0 = Int(QCDNUM.qstore("read", 2))
        ix = QCDNUM.ixfrmx(xb)
        f = QCDNUM.bvalij(iset, i, ix, iq0, 1)
    end

    return f
end

def = Float64.([0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,
                0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,
                0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
                0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
                0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]);

@testset "Evolution & interpolation" begin

    # C-pointer to func
    func_c = @cfunction(func, Float64, (Ref{Int32}, Ref{Float64}))

    # Initialise
    QCDNUM.qcinit(-6, "")

    # Grids and weights
    nx = QCDNUM.gxmake(Float64.([1.0e-4]), Int32.([1]), 1, 100, 3)
    nq = QCDNUM.gqmake(Float64.([2e0, 1e4]), Float64.([1e0, 1e0]), 2, 50)
    nw = QCDNUM.fillwt(1)

    # Set parameters
    QCDNUM.setord(3)
    QCDNUM.setalf(0.364, 2.0)

    # VFNS
    iqc = QCDNUM.iqfrmq(3.0)
    iqb = QCDNUM.iqfrmq(25.0)
    QCDNUM.setcbt(0, iqc, iqb, 0)

    # Evolution
    iq0 = QCDNUM.iqfrmq(2.0)
    eps = QCDNUM.evolfg(1, func_c, def, iq0)
    @test typeof(eps) == Float64
    @test eps < 0.1

    # For SGNS
    QCDNUM.qstore("write", 1, float(1))
    QCDNUM.qstore("write", 2, float(iq0))

    # Interpolation
    x = 1e-3
    q = 1e3
    pdf = QCDNUM.allfxq(1, x, q, 0, 1)
    @test pdf[1] == 0.0
    @test size(pdf)[1] == 13

    ix = QCDNUM.ixfrmx(x)
    iq = QCDNUM.iqfrmq(q)

    for id in 0:12
        pdf_ij = QCDNUM.bvalij(1, id, ix, iq, 1)
        pdf_xq = QCDNUM.bvalxq(1, id, x, q, 1)
        @test isapprox(pdf_ij, pdf_xq, rtol=0.1)
    end

    for id in -6:6
        pdf_ij = QCDNUM.fvalij(1, id, ix, iq, 1)
        pdf_xq = QCDNUM.fvalxq(1, id, x, q, 1)
        @test isapprox(pdf_ij, pdf_xq, rtol=0.1)
    end

    # Check number of PDF tables
    ntables = QCDNUM.nptabs(1)
    @test ntables == 13

    # Check evolution type
    type = QCDNUM.ievtyp(1)
    @test type == 1

    asmz, a, b = QCDNUM.asfunc(8315.25)
    @test asmz ≈ 0.1180167650168159
    @test a == 5
    @test b == 0

end

@testset "Structure functions" begin

    # Initialise
    QCDNUM.qcinit(-6, "")

    # C-pointer to func
    func_c = @cfunction(func, Float64, (Ref{Int32}, Ref{Float64}))

    # Initialise
    QCDNUM.qcinit(-6, "")

    # Grids and weights
    nx = QCDNUM.gxmake(Float64.([1.0e-4]), Int32.([1]), 1, 100, 3)
    nq = QCDNUM.gqmake(Float64.([2e0, 1e4]), Float64.([1e0, 1e0]), 2, 50)
    nw = QCDNUM.fillwt(1)

    # ZMSTF
    ntot, nused_bf = QCDNUM.zmwords()
    nw = QCDNUM.zmfillw()
    ntot, nused_af = QCDNUM.zmwords()
    @test typeof(nw) == Int32
    @test nw == nused_af - nused_bf

    # Set parameters
    QCDNUM.setord(3)
    QCDNUM.setalf(0.364, 2.0)

    # VFNS
    iqc = QCDNUM.iqfrmq(3.0)
    iqb = QCDNUM.iqfrmq(25.0)
    QCDNUM.setcbt(0, iqc, iqb, 0)

    # Evolution
    iq0 = QCDNUM.iqfrmq(2.0)
    eps = QCDNUM.evolfg(1, func_c, def, iq0)
    @test typeof(eps) == Float64

    Au = 4.0/9.0
    Ad = 1.0/9.0
    w = [0., Ad, Au, Ad, Au, Ad, 0., Ad, Au, Ad, Au, Ad, 0.]
    x = 1e-3
    q = 1e3

    for istf in 1:4
        f = QCDNUM.zmstfun(istf, w, [x], [q], 1, 0)
        @test typeof(f[1]) == Float64
    end

end

# Define function to read from QCDNUM into spline
function func_sp(ix, iq, first)::Float64

    # deref ptrs
    ix = ix[]
    iq = iq[]

    # Read iset and ipdf
    iset = Int32(QCDNUM.dsp_uread(1))
    ipdf = Int32(QCDNUM.dsp_uread(2))

    return QCDNUM.fvalij(iset, ipdf, ix, iq, 1)
end

@testset "SPLINT" begin

        # Initialise
    QCDNUM.qcinit(-6, "")

    # C-pointer to func
    func_c = @cfunction(func, Float64, (Ref{Int32}, Ref{Float64}))

    # Initialise
    QCDNUM.qcinit(-6, "")

    # Grids and weights
    nx = QCDNUM.gxmake(Float64.([1.0e-4]), Int32.([1]), 1, 100, 3)
    nq = QCDNUM.gqmake(Float64.([2e0, 1e4]), Float64.([1e0, 1e0]), 2, 50)
    nw = QCDNUM.fillwt(1)

    # Set parameters
    QCDNUM.setord(3)
    QCDNUM.setalf(0.364, 2.0)

    # VFNS
    iqc = QCDNUM.iqfrmq(3.0)
    iqb = QCDNUM.iqfrmq(25.0)
    QCDNUM.setcbt(0, iqc, iqb, 0)

    # Evolution
    iq0 = QCDNUM.iqfrmq(2.0)
    eps = QCDNUM.evolfg(1, func_c, def, iq0)
    @test typeof(eps) == Float64

    # SPLINT
    func_sp_c = @cfunction(func_sp, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))

    iset = 1
    ipdf = 1

    # Initialise
    QCDNUM.ssp_spinit(100)

    # Store iset and ipdf
    QCDNUM.ssp_uwrite(1, Float64(iset))
    QCDNUM.ssp_uwrite(2, Float64(ipdf))

    # Make spline object
    iasp = QCDNUM.isp_s2make(5, 5)
    @test typeof(iasp) == Int32

    # Fill the spline and set no kinematic limit
    QCDNUM.ssp_s2fill(iasp, func_sp_c, 0.0)

    # Check spline type
    type = QCDNUM.isp_splinetype(iasp)
    @test type == 2

    # Check limits
    nu, u1, u2, nv, v1, v2, n = QCDNUM.ssp_splims(iasp)
    @test nu == 22
    @test u1 ≈ 1e-4
    @test u2 == 1.0
    @test nv == 12
    @test v1 == 2.0
    @test v2 ≈ 1e4
    @test n == nu * nv

    # Evalute function and integral
    x = 0.1
    q = 100.0
    f = QCDNUM.dsp_funs2(iasp, x, q, 1)
    @test isapprox(f, 0.408, rtol=0.01)

    x1 = 0.01
    x2 = 0.1
    q1 = 10.0
    q2 = 100.0
    rs = 370.0
    np = 4
    integral = QCDNUM.dsp_ints2(iasp, x1, x2, q1, q2, rs, np)
    @test isapprox(integral, 3.709, rtol = 0.01)

    # Copy nodes
    u_nodes = QCDNUM.ssp_unodes(iasp, nu, nu)
    v_nodes = QCDNUM.ssp_vnodes(iasp, nv, nv)
    @test size(u_nodes)[1] == nu
    @test size(v_nodes)[1] == nv

    @test QCDNUM.ssp_nprint(iasp) == nothing

    for i in 1:3
        @test QCDNUM.ssp_extrapu(iasp, i) == nothing
        @test QCDNUM.ssp_extrapv(iasp, i) == nothing
    end

    @test QCDNUM.ssp_erase(iasp) == nothing

    # 1D spline
    iasp = QCDNUM.isp_sxmake(5)
    @test typeof(iasp) == Int32

    iasp = QCDNUM.isp_sqmake(5)
    @test typeof(iasp) == Int32

end


text:
"""
    katz(g, alpha = 0.1)

Compute the Katz centrality for a static graph `g`.

# References:

1. L. Katz A new index derived from sociometric data analysis. Psychometrika, 18:39-43, 1953.
"""
function katz(g::AbstractStaticGraph, alpha::Real = 0.1)
    n = num_nodes(g)
    ns = nodes(g)
    v = ones(Float64, n)
    A = sparse_adjacency_matrix(g)
    spI = speye(Float64, n)
    rates = (spI - alpha * A)\v
    return [(node, rates[node_index(g, node)]) for node in ns]
end

"""
    katz(g, alpha = 0.3)
    katz(g, alpha, beta; mode = :broadcast)

Computes the katz centrality for an evolving graph `g`, where `alpha` and `beta` are scalars. `alpha` controls the influence of long walks and `beta` controls the influence of walks happened long time ago. By default, `mode = :broadcast` computes the broadcast centrality. Otherwise if `mode = :receive`, we compute the receiving centrality.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> using EvolvingGraphs.Centrality

julia> g = evolving_graph_from_arrays(["A", "B", "B", "C", "E", "A", "B", "D"], ["B", "F", "G", "E", "G", "B", "F", "F"], [1,1,1,2,2,2,2,2])
Directed EvolvingGraph 7 nodes, 8 static edges, 2 timestamps

julia> katz(g)
7-element Array{Tuple{EvolvingGraphs.Node{String},Float64},1}:
 (Node(A), 0.510104)
 (Node(B), 0.494488)
 (Node(F), 0.260257)
 (Node(G), 0.260257)
 (Node(C), 0.361757)
 (Node(E), 0.338334)
 (Node(D), 0.338334)

julia> katz(g, 0.3, 0.4, mode = :receive)
7-element Array{Tuple{EvolvingGraphs.Node{String},Float64},1}:
 (Node(A), 0.0)
 (Node(B), 0.658926)
 (Node(F), 1.0)
 (Node(G), 0.658926)
 (Node(C), 0.0)
 (Node(E), 0.204852)
 (Node(D), 0.0)
```

# References:

1. P. Grindrod, D. J. Higham, M. C. Parsons and E. Estrada Communicability across evolving networks. Physical Review E, 83 2011.

2. P. Grindrod and D. J. Higham, A matrix iteration for dynamic network summaries. SIAM Review, 55 2013.
"""
function katz(g::AbstractEvolvingGraph, alpha::Real = 0.3)
    n = num_nodes(g)
    ns = nodes(g)
    ts = unique_timestamps(g)
    ts_sorted = sort(ts, rev=true)
    v = ones(Float64, n)
    A = spzeros(Float64, n, n)
    spI = speye(Float64, n)
    for t in ts_sorted
        A =  sparse_adjacency_matrix(g,t)
        v = (spI - alpha*A)\v
        v =  v/norm(v)
    end
    return [(node, v[node_index(g, node)]) for node in ns]
end

function katz(g::AbstractEvolvingGraph, alpha::Real, beta::Real; mode::Symbol = :broadcast)
    n = num_nodes(g)
    ns = nodes(g)
    ts = unique_timestamps(g)
    S = spzeros(Float64, n, n)
    A = spzeros(Float64, n, n)
    spI = speye(Float64, n)
    v = Array{Float64}(n)
    Δt = 1.
    for t in ts
        Δt += 0.01
        A =  sparse_adjacency_matrix(g,t)
        S =  (spI + e^( -beta * Δt) * S) * (spI + alpha * A) - spI
        S =  S/norm(S,1)
    end

    v = mode == :broadcast ? S * ones(Float64,n) : S' * ones(Float64,n)

    return [(node, v[node_index(g, node)]) for node in ns]
end

ext: jl


text:
#=
Optimized search functions for UTFStr (UTF-8 encoding)

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#
function _srch_cp(::Fwd, ::CodeUnitMulti, str::Str{UTF8CSE}, cp::AbsChar, pos, len)
    @preserve str begin
        (ch = codepoint(cp)) < 0x80 &&
            return _srch_codeunit(Fwd(), pointer(str), ch%UInt8, pos, len)
        beg = pointer(str)
        pnt = beg + pos - 1
        fin = beg + len
        if ch <= 0x7ff
            b1, b2 = get_utf8_2(ch)
            while (pnt = _fwd_memchr(pnt + 1, b2, fin)) != C_NULL
                eq_bytes(pnt - 1, b1) && return Int(pnt - beg)
            end
        elseif ch <= 0xffff
            b1, b2, b3 = get_utf8_3(ch)
            pnt += 1
            while (pnt = _fwd_memchr(pnt + 1, b3, fin)) != C_NULL
                eq_bytes(pnt - 2, b1, b2) && return Int(pnt - beg) - 1
            end
        else
            b1, b2, b3, b4 = get_utf8_4(ch)
            pnt += 2
            while (pnt = _fwd_memchr(pnt + 1, b4, fin)) != C_NULL
                eq_bytes(pnt - 3, b1, b2, b3) && return Int(pnt - beg) - 2
            end
        end
        0
    end
end

function _srch_cp(::Rev, ::CodeUnitMulti, str::Str{UTF8CSE}, cp::AbsChar, pos, len)
    @preserve str begin
        (ch = codepoint(cp)) < 0x80 && return _srch_codeunit(Rev(), pointer(str), ch%UInt8, pos)
        init = beg = pointer(str)
        @inbounds pnt = beg + nextind(str, pos)
        if ch <= 0x7ff
            pnt > (beg += 1) || return 0
            b1, b2 = get_utf8_2(ch)
            while (pnt = _rev_memchr(beg, b2, pnt - beg)) != C_NULL
                eq_bytes(pnt - 1, b1) && return Int(pnt - init)
                pnt > beg || break
            end
        elseif ch <= 0xffff
            pnt > (beg += 2) || return 0
            b1, b2, b3 = get_utf8_3(ch)
            while (pnt = _rev_memchr(beg, b3, pnt - beg)) != C_NULL
                eq_bytes(pnt - 2, b1, b2) && return Int(pnt - 1 - init)
                pnt > beg || break
            end
        else
            pnt > (beg += 3) || return 0
            b1, b2, b3, b4 = get_utf8_4(ch)
            while (pnt = _rev_memchr(beg, b4, pnt - beg)) != C_NULL
                eq_bytes(pnt - 3, b1, b2, b3) && return Int(pnt - 2 - init)
                pnt > beg || break
            end
        end
    0
    end
end


text:
using MetaICVI
using Documenter

DocMeta.setdocmeta!(MetaICVI, :DocTestSetup, :(using MetaICVI); recursive=true)

makedocs(;
    modules=[MetaICVI],
    authors="Sasha Petrenko",
    repo="https://github.com/AP6YC/MetaICVI.jl/blob/{commit}{path}#{line}",
    sitename="MetaICVI.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://AP6YC.github.io/MetaICVI.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial" => [
            "Guide" => "man/guide.md",
            "Examples" => "man/examples.md",
            "Contributing" => "man/contributing.md",
            "Index" => "man/full-index.md"
        ]
    ],
)

deploydocs(;
    repo="github.com/AP6YC/MetaICVI.jl",
    devbranch = "main",
)


text:
# indicator of the L2 norm sphere with given radius

export IndSphereL2

"""
**Indicator of a Euclidean sphere**

    IndSphereL2(r=1.0)

Returns the indicator function of the set
```math
S = \\{ x : \\|x\\| = r \\},
```
where ``\\|\\cdot\\|`` is the ``L_2`` (Euclidean) norm. Parameter `r` must be positive.
"""
struct IndSphereL2{R <: Real} <: ProximableFunction
    r::R
    function IndSphereL2{R}(r::R) where {R <: Real}
        if r <= 0
            error("parameter r must be positive")
        else
            new(r)
        end
    end
end

is_set(f::IndSphereL2) = true

IndSphereL2(r::R=1.0) where {R <: Real} = IndSphereL2{R}(r)

function (f::IndSphereL2)(x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    if isapprox(norm(x), f.r, atol=eps(R), rtol=sqrt(eps(R)))
        return R(0)
    end
    return R(Inf)
end

function prox!(y::AbstractArray{T}, f::IndSphereL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    normx = norm(x)
    if normx > 0 # zero-zero?
        scal = f.r/normx
        for k in eachindex(x)
            y[k] = scal*x[k]
        end
    else
        normy = R(0)
        for k in eachindex(x)
            y[k] = randn()
            normy += y[k]*y[k]
        end
        normy = sqrt(normy)
        y .*= f.r/normy
    end
    return R(0)
end

fun_name(f::IndSphereL2) = "indicator of an L2 norm sphere"
fun_dom(f::IndSphereL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndSphereL2) = "x ↦ 0 if ||x|| = r, +∞ otherwise"
fun_params(f::IndSphereL2) = "r = $(f.r)"

function prox_naive(f::IndSphereL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    normx = norm(x)
    if normx > 0
        y = x*f.r/normx
    else
        y = randn(size(x))
        y *= f.r/norm(y)
    end
    return y, R(0)
end


text:
function main(r::Robot)
    for i in ((North, East), (North, West), (South, East), (South, West))
        go_special!(r, i)
        return_special!(r, i)
    end
    putmarker!(r)
end

function return_special!(r,side)
    while (ismarker(r))
        move!(r, inverse_side(side[1]))
        move!(r, inverse_side(side[2]))
    end
end

function go_special!(r, side::NTuple{2,HorizonSide})
    while isborder(r,side[1]) == false  &&   isborder(r,side[2]) == false
        move!(r,side[1])
        move!(r,side[2])
        putmarker!(r)
    end
end

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
module IntrinsicDecoherence
    include("Funcs.jl")
    include("Funcs_Displaced.jl")
end #module

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
module BenchFilterMapReduce

using BenchmarkTools
using Transducers

function manual_filter_map_reduce(xs)
    acc = zero(eltype(xs))
    for x in xs
        if -0.5 < x < 0.5
            acc += 2x
        end
    end
    return acc
end

suite = BenchmarkGroup()

xf = Filter(x -> -0.5 < x < 0.5) |> Map(x -> 2x)
n = 10^5

let xs = randn(10)
    @assert manual_filter_map_reduce(xs) == mapfoldl(xf, +,  xs, init=0.0)
end

suite["xf"] = @benchmarkable(
    mapfoldl($xf, +, xs, init=0.0),
    setup=(xs = randn($n)))
suite["man"] = @benchmarkable(
    manual_filter_map_reduce(xs),
    setup=(xs = randn($n)))

end  # module
BenchFilterMapReduce.suite


text:
using LinearAlgebra

abstract type BlockSolution <: AbstractBlockMatrix{Cdouble} end
struct PrimalSolution <: BlockSolution
    problem::SDPAProblem
end
getptr(X::PrimalSolution, blk) = getResultYMat(X.problem, blk)
struct VarDualSolution <: BlockSolution
    problem::SDPAProblem
end
getptr(X::VarDualSolution, blk) = getResultXMat(X.problem, blk)
nblocks(X::BlockSolution) = getBlockNumber(X.problem)
function block(X::BlockSolution, blk::Integer)
    if blk <= 0 || blk > getBlockNumber(X.problem)
        throw(BoundsError(X, blk))
    end
    n = getBlockSize(X.problem, blk)
    if getBlockType(X.problem, blk) == SDP
        unsafe_wrap(Array, getptr(X, blk), (n, n))
    else
        Diagonal(unsafe_wrap(Array, getptr(X, blk), (n,)))
    end
end
# Needed by MPB_wrapper
function Base.getindex(A::BlockSolution, i::Integer)
    block(A, i)
end


text:
using BenchmarkTools
using Contour
using StaticArrays
using Test
#using Plots

const suite = BenchmarkGroup()
suite["contour"] = BenchmarkGroup()

#circle(v, r) = sqrt(sum(v.^2)) - r

include(joinpath(@__DIR__,"../test/testdata.jl"))
v = Contour.contours(x, y, z)
@show typeof(v)

suite["contour"]["testdata"] = @benchmarkable Contour.contours($x,$y,$z)

results = run(suite)

for trial in results
    ctx = IOContext(stdout, :verbose => true, :compact => false)
    println(ctx)
    println(ctx, trial.first)
    println(ctx, trial.second)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_script/scrubber.jl
# A context used to "scrub" type specific storage, assertions, and other
# language features from a user-defined function.
#############################################################################

@context ScrubCtx

# Cassette specific functions for scrubbing udfs of objects that interfere with overloading
Cassette.overdub(ctx::ScrubCtx, ::typeof(typeassert), x::Real, type::Type) = x
function Cassette.overdub(ctx::ScrubCtx, ::typeof(zeros), t, dims...)
    if t <: AbstractFloat
        return zeros(Real, dims...)
    elseif t <: Integer
        return zeros(Real, dims...)
    end
    return zeros(t, dims...)
end
#overdub(ctx::ScrubCtx, ::typeof(zero), type::Type) = zero(type)
function Cassette.overdub(ctx::ScrubCtx, ::typeof(hcat), A...)
    vA = hcat(A...)
    sz = size(vA)
    vR = zeros(Real, sz...)
    vR[:] = vA[:]
    return vR
end
function Cassette.overdub(ctx::ScrubCtx, ::typeof(vcat), A...)
    vA = vcat(A...)
    sz = size(vA)
    vR = zeros(Real, sz...)
    vR[:] = vA[:]
    return vR
end

# scrub type assertion from definition of udf
function generated_scrubbed_method(f::Function, n::Int, inplace = false)
    ci = code_lowered(f)
    @assert(length(ci) == 1)
    function g end
    addmethod!(Tuple{typeof(g), Any, Any}, ci[1])
    addmethod!(Tuple{typeof(g), Any}, ci[1])
    addmethod!(Tuple{typeof(g),
                     SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true},
                     SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}}, ci[1])
    return g
end

"""
    scrub(f::Function, n::Int, inplace = false)

Replaces storage objects and circumvents assertions that present a UDF from
being overloaded with subtype Real objects by creating a function which
overdubs `f` in `ScrubCtx`.
"""
function scrub(f::Function, n::Int, inplace = false)
    if inplace
        if hasmethod(f, Tuple{Array{Float64,1}})
            println("scrub arc 1...")
            return false, (y,x) -> overdub(ScrubCtx(), f, y, x)
        else
            println("scrub arc 2...")
            fna = generated_scrubbed_method(f, n, true)
            return true, (y,x) -> overdub(ScrubCtx(), fna, y, x)
        end
    else
        if hasmethod(f, Tuple{Array{Float64,1}})
            return false, x -> overdub(ScrubCtx(), f, x)
        else
            fna2 = generated_scrubbed_method(f, n)
            return true, x -> overdub(ScrubCtx(), fna2, x)
        end
    end
end

"""
    scrub!(d::_NLPData)

Applies scrub to every user-defined function in the a `_NLPData` structure.
"""
function scrub!(d::_NLPData)

    # scrub multivariant
    user_ops = d.user_operators
    mvop_num = length(user_ops.multivariate_operator_evaluator)
    for i in 1:mvop_num
        evalr = user_ops.multivariate_operator_evaluator[i]
        f_old = evalr.f
        flag, fnew = scrub(evalr.f, evalr.len)
        ftemp = (out, a) -> gradient!(out, x -> f_old(x), a)
        fnew2 = (y, x) -> overdub(ScrubCtx(), ftemp, y, x)
        evalr.f = fnew
        evalr.∇f = fnew2
        user_ops.multivariate_operator_evaluator[i] = evalr
    end

    # scrub univariants
    svop_num = length(user_ops.univariate_operator_f)

    for i in 1:svop_num
        flag, fnew = scrub(user_ops.univariate_operator_f[i], 1)
        user_ops.univariate_operator_f[i] = fnew
        if flag
            f = user_ops.univariate_operator_f[i]
            fprime = x -> derivative(f, x)
            user_ops.univariate_operator_fprime[i] = fprime
            user_ops.univariate_operator_fprimeprime[i] = x -> derivative(fprime, x)
        end

        fnew = scrub(user_ops.univariate_operator_fprime[i], 1)
        user_ops.univariate_operator_fprime[i] = fnew
        fnew = scrub(user_ops.univariate_operator_fprimeprime[i], 1)
        user_ops.univariate_operator_fprimeprime[i] = fnew
    end

    # reassign evaluator
    d.user_operators = user_ops
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
module SeisAcoustic

    # the dependency of this module
    using LinearAlgebra,
          Printf,           # formated print
          DelimitedFiles,   # delimited text file
          SparseArrays,     # finite-difference sparse matrix
          DSP,              # smoothing and convolution
          FFTW,             # fourier transform
          OffsetArrays,     # User-defined index range (similar to fortran)
          Distributed       # distribute parallel computing

    # the path to the binary dependency (sparse matrix multiplication)
    const spmatveclib  = abspath(joinpath(splitdir(Base.source_path())[1],"..","deps","builds","spmatvec.so"))

    # overloading Base function
    import Base.convert, Base.show, Base.isequal, Base.zeros

    include("dataio/dataio.jl")          # read and write segy data, internally defined regular sampled data (borrow from rsf)
    include("finitediff/finitediff.jl")  # finite-difference method for acoustic wave equation
    include("imaging/imaging.jl")
    include("fwi/fwi.jl")
    include("solver/solver.jl")

end # module


text:
"""
    imshow(arg1; kw...)

Is a simple front end to the [`grdimage`](@ref)  [`grdview`](@ref) programs that accepts GMTgrid, GMTimage,
2D array of floats or strings with file names of grids or images. The normal options of the *grdimage*
and *grdview* programs also apply here but some clever guessing of suitable necessary parameters is done
if they are not provided. Contrary to other image producing modules the "show' keyword is not necessary to
display the image. Here it is set by default. If user wants to use *imshow* to create layers of a more complex
fig he can use *show=false* for the intermediate layers.

This module uses some internal logic to decide whether use `grdimge` or `grdview`. Namely, when the `view`
option is used `grdview` is choosed and a default vertical scale is assigned. However, sometimes we want
a rotated plot, optionallt tilted, but not 3D view. In that case use the option `flat=true`, which forces
the use of `grdimage`.

# Examples
```julia-repl
# Plot vertical shaded illuminated view of the Mexican hat
julia> G = gmt("grdmath -R-15/15/-15/15 -I0.3 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
julia> imshow(G, shade="+a45")

# Same as above but add automatic contours
julia> imshow(G, shade="+a45", contour=true)

# Plot a random heat map
julia> imshow(rand(128,128))

# Display a web downloaded jpeg image wrapped into a sinusoidal projection
julia> imshow("http://larryfire.files.wordpress.com/2009/07/untooned_jessicarabbit.jpg", region="d", frame="g", proj="I15", img_in="r", fmt=:jpg)
```
See also: [`grdimage`](@ref)
"""
function imshow(arg1, x=nothing, y=nothing; kw...)
	# Take a 2D array of floats and turn it into a GMTgrid or if input is a string assume it's a file name
	# In this later case try to figure if it's a grid or an image and act accordingly.
	is_image = false
	if (isa(arg1, String))		# If it's string it has to be a file name. Check extension to see if is an image
		ext = splitext(arg1)[2]
		if (ext == "" && arg1[1] != '@' && !isfile(arg1))
			G = mat2grid(arg1, x, y)
		else
			ext = lowercase(ext)
			(ext == ".jpg" || ext == ".tif" || ext == ".tiff" || ext == ".png" || ext == ".bmp" || ext == ".gif") && (is_image = true)
			G = (arg1[1] == '@') ? arg1 : gmtread(arg1)			# If it screws ...
		end
	elseif (isa(arg1, Array{UInt8}) || isa(arg1, Array{UInt16,3}))
		G = mat2img(arg1; kw...)
	else
		G = mat2grid(arg1, x, y, reg=1)			# For displaying, pixel registration is more appropriate
	end

	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true

	if (is_image)
		if (haskey(d, :D) || haskey(d, :img_in) || haskey(d, :image_in))	# OK, user set -D so don't repeat
			grdimage(G; show=see, kw...)
		else
			grdimage(G; D=1, show=see, kw...)
		end
	else
		if (isa(G, String))  grdimage(G; show=see, kw...)		# String when fname is @xxxx
		else                 imshow(G; kw...)					# Call the specialized method
		end
	end
end

function imshow(arg1::GMTgrid; kw...)
	# Here the default is to show, but if a 'show' was used let it rule
	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true
	if ((cont_opts = find_in_dict(d, [:contour])[1]) !== nothing)
		new_see = see
		see = false			# because here we know that 'see' has to wait till last command
	end
	opt_p, = parse_common_opts(d, "", [:p], true)
	til = find_in_dict(d, [:T :no_interp :tiles])[1]
	flat = (find_in_dict(d, [:flat])[1] !== nothing)		# If true, force the use of grdimage
	if (flat || (opt_p == "" && til === nothing))
		(flat && opt_p != "") && (d[:p] = opt_p[4:end])		# Restore the meanwhile deleted -p option
		R = grdimage("", arg1; show=see, d...)
	else
		zsize = ((val = find_in_dict(d, [:JZ :Jz :zscale :zsize])[1]) !== nothing) ? val : 8
		srf = ((val = find_in_dict(d, [:Q :surf :surftype])[1]) !== nothing) ? val : "i100"
		if (til !== nothing)		# This forces some others
			srf = nothing			# These are mutually exclusive
			zsize = nothing
			opt_p = " -p180/90"
		end
		R = grdview("", arg1; show=see, p=opt_p[4:end], JZ=zsize, Q=srf, T=til, d...)
	end
	if (isa(cont_opts, Bool))				# Automatic contours
		R = grdcontour!(arg1; J="", show=new_see)
	elseif (isa(cont_opts, NamedTuple))		# Expect a (cont=..., annot=..., ...)
		R = grdcontour!(arg1; J="", show=new_see, cont_opts...)
	end
	return R
end

function imshow(arg1::GMTimage; kw...)
	# Here the default is to show, but if a 'show' was used let it rule
	d = KW(kw)
	see = (!haskey(d, :show)) ? true : see = d[:show]	# No explicit 'show' keyword means show=true
	grdimage("", arg1; D=true, show=see, kw...)
end

function imshow(x, y, f::Function; kw...)
	G = mat2grid(f, x, y)
	imshow(G; kw...)
end
imshow(x, f::Function; kw...) = imshow(x, x, f::Function; kw...)
imshow(f::Function, x; kw...) = imshow(x, x, f::Function; kw...)
imshow(f::Function, x, y; kw...) = imshow(x, y, f::Function; kw...)
imshow(x, y, f::String; kw...) = imshow(f::String, x, y; kw...)
imshow(x, f::String; kw...) = imshow(f::String, x, x; kw...)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^2*sin(a+b*log(c*x^n)), x) == :(3*x^3*(9+b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x^3*(9+b^2*n^2)^-1*cos(a+b*log(c*x^n)))
@test integrate(x*sin(a+b*log(c*x^n)), x) == :(2*x^2*(4+b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x^2*(4+b^2*n^2)^-1*cos(a+b*log(c*x^n)))
@test integrate(sin(a+b*log(c*x^n)), x) == :(x*(1+b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x*(1+b^2*n^2)^-1*cos(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n)), x) == :(-1*b^-1*n^-1*cos(a+b*log(c*x^n)))
@test integrate(x^-2*sin(a+b*log(c*x^n)), x) == :(-1*x^-1*(1+b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x^-1*(1+b^2*n^2)^-1*cos(a+b*log(c*x^n)))
@test integrate(x^-3*sin(a+b*log(c*x^n)), x) == :(-2*x^-2*(4+b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x^-2*(4+b^2*n^2)^-1*cos(a+b*log(c*x^n)))
@test integrate(x^2*sin(a+b*log(c*x^n))^2, x) == :(3*x^3*(9+4*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+2*b^2*n^2*x^3*(27+12*b^2*n^2)^-1+-2*b*n*x^3*(9+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x*sin(a+b*log(c*x^n))^2, x) == :(x^2*(2+2*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+b^2*n^2*x^2*(4+4*b^2*n^2)^-1+-1*b*n*x^2*(2+2*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(sin(a+b*log(c*x^n))^2, x) == :(x*(1+4*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+2*x*b^2*n^2*(1+4*b^2*n^2)^-1+-2*b*n*x*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n))^2, x) == :((1/2)*log(x)+-1//2*b^-1*n^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-2*sin(a+b*log(c*x^n))^2, x) == :(-1*x^-1*(1+4*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+-2*b^2*n^2*x^-1*(1+4*b^2*n^2)^-1+-2*b*n*x^-1*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-3*sin(a+b*log(c*x^n))^2, x) == :(-1*x^-2*(2+2*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+-1*b^2*n^2*x^-2*(4+4*b^2*n^2)^-1+-1*b*n*x^-2*(2+2*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^2*sin(a+b*log(c*x^n))^3, x) == :(x^3*(3+3*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3+-2*b^3*n^3*x^3*(27+3*b^4*n^4+30*b^2*n^2)^-1*cos(a+b*log(c*x^n))+2*b^2*n^2*x^3*(9+b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))+-1*b*n*x^3*(3+3*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2*cos(a+b*log(c*x^n)))
@test integrate(x*sin(a+b*log(c*x^n))^3, x) == :(2*x^2*(4+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3+-6*b^3*n^3*x^2*(16+9*b^4*n^4+40*b^2*n^2)^-1*cos(a+b*log(c*x^n))+12*b^2*n^2*x^2*(16+9*b^4*n^4+40*b^2*n^2)^-1*sin(a+b*log(c*x^n))+-3*b*n*x^2*(4+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2*cos(a+b*log(c*x^n)))
@test integrate(sin(a+b*log(c*x^n))^3, x) == :(x*(1+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3+-6*x*b^3*n^3*(1+9*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))+6*x*b^2*n^2*(1+9*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))+-3*b*n*x*(1+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2*cos(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n))^3, x) == :(-1*b^-1*n^-1*cos(a+b*log(c*x^n))+1//3*b^-1*n^-1*cos(a+b*log(c*x^n))^3)
@test integrate(x^-2*sin(a+b*log(c*x^n))^3, x) == :(-1*x^-1*(1+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3+-6*b^2*n^2*x^-1*(1+9*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))+-6*b^3*n^3*x^-1*(1+9*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))+-3*b*n*x^-1*(1+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2*cos(a+b*log(c*x^n)))
@test integrate(x^-3*sin(a+b*log(c*x^n))^3, x) == :(-2*x^-2*(4+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3+-12*b^2*n^2*x^-2*(16+9*b^4*n^4+40*b^2*n^2)^-1*sin(a+b*log(c*x^n))+-6*b^3*n^3*x^-2*(16+9*b^4*n^4+40*b^2*n^2)^-1*cos(a+b*log(c*x^n))+-3*b*n*x^-2*(4+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2*cos(a+b*log(c*x^n)))
@test integrate(x^2*sin(a+b*log(c*x^n))^4, x) == :(3*x^3*(9+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^4+8*b^4*n^4*x^3*(81+64*b^4*n^4+180*b^2*n^2)^-1+36*b^2*n^2*x^3*(81+64*b^4*n^4+180*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+-24*b^3*n^3*x^3*(81+64*b^4*n^4+180*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n))+-4*b*n*x^3*(9+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n)))
@test integrate(x*sin(a+b*log(c*x^n))^4, x) == :(x^2*(2+8*b^2*n^2)^-1*sin(a+b*log(c*x^n))^4+3*b^4*n^4*x^2*(4+16*b^4*n^4+20*b^2*n^2)^-1+3*b^2*n^2*x^2*(2+8*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+-1*b*n*x^2*(1+4*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n))+-3*b^3*n^3*x^2*(2+8*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(sin(a+b*log(c*x^n))^4, x) == :(x*(1+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^4+24*x*b^4*n^4*(1+20*b^2*n^2+64*b^4*n^4)^-1+12*x*b^2*n^2*(1+20*b^2*n^2+64*b^4*n^4)^-1*sin(a+b*log(c*x^n))^2+-24*x*b^3*n^3*(1+20*b^2*n^2+64*b^4*n^4)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n))+-4*b*n*x*(1+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n))^4, x) == :(3//8*log(x)+-3//8*b^-1*n^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n))+-1//4*b^-1*n^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n)))
@test integrate(x^-2*sin(a+b*log(c*x^n))^4, x) == :(-1*x^-1*(1+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^4+-24*b^4*n^4*x^-1*(1+20*b^2*n^2+64*b^4*n^4)^-1+-12*b^2*n^2*x^-1*(1+20*b^2*n^2+64*b^4*n^4)^-1*sin(a+b*log(c*x^n))^2+-24*b^3*n^3*x^-1*(1+20*b^2*n^2+64*b^4*n^4)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n))+-4*b*n*x^-1*(1+16*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n)))
@test integrate(x^-3*sin(a+b*log(c*x^n))^4, x) == :(-1*x^-2*(2+8*b^2*n^2)^-1*sin(a+b*log(c*x^n))^4+-3*b^4*n^4*x^-2*(4+16*b^4*n^4+20*b^2*n^2)^-1+-3*b^2*n^2*x^-2*(2+8*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))^2+-1*b*n*x^-2*(1+4*b^2*n^2)^-1*sin(a+b*log(c*x^n))^3*cos(a+b*log(c*x^n))+-3*b^3*n^3*x^-2*(2+8*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(sin(log(a+b*x)), x) == :((1/2)*b^-1*(a+b*x)*sin(log(a+b*x))+-1//2*b^-1*(a+b*x)*cos(log(a+b*x)))
@test integrate(x^m*sin(a+(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n)), x) == :(-1//4*n^-1*x^(1+m)*(c*x^n)^(n^-1*(1+m))*(-1*n^-2*(1+m)^2)^-1//2*exp(a*n^-1*(-1*n^-2*(1+m)^2)^-1//2*(1+m))+(1/2)*n^-1*x^(1+m)*(c*x^n)^(-1*n^-1*(1+m))*(-1*n^-2*(1+m)^2)^-1//2*(1+m)*exp(a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)*log(x))
@test integrate(x^2*sin(a+3*(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//12*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+-1//2*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(-3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x*sin(a+2*(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//8*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(2*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+-1//2*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(-2*n^-1)*exp(a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(sin(a+(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//4*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+-1//2*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(-1*n^-1)*exp(a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-1*sin(a), x) == :(log(x)*sin(a))
@test integrate(x^-2*sin(a+(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//4*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(-1*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+(1/2)*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-3*sin(a+2*(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//8*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(-2*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+(1/2)*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(2*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^m*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^2, x) == :(x^(1+m)*(2+2m)^-1+-1*x^(1+m)*(c*x^n)^(n^-1*(1+m))*(8+8m)^-1*exp(-2*a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)+-1//4*x^(1+m)*(c*x^n)^(-1*n^-1*(1+m))*exp(2*a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)*log(x))
@test integrate(x^2*sin(a+3//2*(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :(1//6*x^3+-1//24*x^3*(c*x^n)^(3*n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))+-1//4*x^3*(c*x^n)^(-3*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x*sin(a+(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :(1//4*x^2+-1//16*x^2*(c*x^n)^(2*n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))+-1//4*x^2*(c*x^n)^(-2*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(sin(a+(1/2)*(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :((1/2)*x+-1//8*x*(c*x^n)^(n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))+-1//4*x*(c*x^n)^(-1*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-1*sin(a)^2, x) == :(sin(a)^2*log(x))
@test integrate(x^-2*sin(a+(1/2)*(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :(-1//2*x^-1+1//8*x^-1*(c*x^n)^(-1*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))+-1//4*x^-1*(c*x^n)^(n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-3*sin(a+(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :(-1//4*x^-2+1//16*x^-2*(c*x^n)^(-2*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))+-1//4*x^-2*(c*x^n)^(2*n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^m*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^3, x) == :(-4*x^(1+m)*(5+5m)^-1*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^3+8*x^(1+m)*(5+5m)^-1*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))+-4//5*n*x^(1+m)*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-2*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))+6//5*n*x^(1+m)*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-2*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^2*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n)))
@test integrate(x^2*sin(a+(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(-3//16*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(-1*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+-1//48*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(3*n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))+3//32*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+1//8*n*x^3*(-1*n^-2)^(1/2)*(c*x^n)^(-3*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x*sin(a+2//3*(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(-9//32*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(-2//3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+-1//32*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(2*n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))+9//64*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(2//3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+1//8*n*x^2*(-1*n^-2)^(1/2)*(c*x^n)^(-2*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(sin(a+1//3*(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(-9//16*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(-1//3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+-1//16*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))+9//32*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(1//3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+1//8*n*x*(-1*n^-2)^(1/2)*(c*x^n)^(-1*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-1*sin(a)^3, x) == :(sin(a)^3*log(x))
@test integrate(x^-2*sin(a+1//3*(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(-9//16*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(1//3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+-1//16*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(-1*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))+9//32*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(-1//3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+-1//8*n*x^-1*(-1*n^-2)^(1/2)*(c*x^n)^(n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-3*sin(a+2//3*(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(-9//32*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(2//3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+-1//32*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(-2*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))+9//64*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(-2//3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+-1//8*n*x^-2*(-1*n^-2)^(1/2)*(c*x^n)^(2*n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^m*sin(a+(1/2)*(-1*(1+m)^2)^(1/2)*log(c*x^2)), x) == :(-1//4*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(1/2+(1/2)*m)*exp(a*(-1*(1+m)^2)^-1//2*(1+m))+(1/2)*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(-1//2+-1//2*m)*(1+m)*exp(a*(-1*(1+m)^2)^(1/2)*(1+m)^-1)*log(x))
@test integrate(sin(a+(1/2)*im*log(c*x^2)), x) == :(-1//2*im*x*(c*x^2)^-1//2*exp(im*a)*log(x)+1//4*im*c*x^3*(c*x^2)^-1//2*exp(-1*im*a))
@test integrate(x^m*sin(a+1//4*(-1*(1+m)^2)^(1/2)*log(c*x^2))^2, x) == :(x^(1+m)*(2+2m)^-1+-1*x^(1+m)*(c*x^2)^(1/2+(1/2)*m)*(8+8m)^-1*exp(2*a*(-1*(1+m)^2)^-1//2*(1+m))+-1//4*x^(1+m)*(c*x^2)^(-1//2+-1//2*m)*exp(-2*a*(-1*(1+m)^2)^-1//2*(1+m))*log(x))
@test integrate(sin(a+1//4*im*log(c*x^2))^2, x) == :((1/2)*x+-1//4*x*(c*x^2)^-1//2*exp(2*im*a)*log(x)+-1//8*c*x^3*(c*x^2)^-1//2*exp(-2*im*a))
@test integrate(x^m*sin(a+1//6*(-1*(1+m)^2)^(1/2)*log(c*x^2))^3, x) == :(-9//32*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(1//6+1//6*m)*exp(a*(-1*(1+m)^2)^-1//2*(1+m))+1//16*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(1/2+(1/2)*m)*exp(3*a*(-1*(1+m)^2)^-1//2*(1+m))+9//16*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(-1//6+-1//6*m)*exp(a*(-1*(1+m)^2)^(1/2)*(1+m)^-1)+-1//8*x^(1+m)*(-1*(1+m)^2)^-1//2*(c*x^2)^(-1//2+-1//2*m)*(1+m)*exp(-3*a*(-1*(1+m)^2)^-1//2*(1+m))*log(x))
@test integrate(sin(a+1//6*im*log(c*x^2))^3, x) == :(-9//16*im*x*(c*x^2)^-1//6*exp(im*a)+9//32*im*x*(c*x^2)^1//6*exp(-1*im*a)+-1//16*im*c*x^3*(c*x^2)^-1//2*exp(-3*im*a)+1//8*im*x*(c*x^2)^-1//2*exp(3*im*a)*log(x))
@test integrate(x^-1*sin(a+b*log(c*x^n))^(1/2), x) == :(2*b^-1*n^-1*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sin(a+b*log(c*x^n))^3//2, x) == :(2//3*b^-1*n^-1*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2)+-2//3*b^-1*n^-1*sin(a+b*log(c*x^n))^(1/2)*cos(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n))^-1//2, x) == :(2*b^-1*n^-1*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sin(a+b*log(c*x^n))^-3//2, x) == :(-2*b^-1*n^-1*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2)+-2*b^-1*n^-1*sin(a+b*log(c*x^n))^-1//2*cos(a+b*log(c*x^n)))
@test integrate(x^-1*sin(a+b*log(c*x^n))^-5//2, x) == :(2//3*b^-1*n^-1*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2)+-2//3*b^-1*n^-1*sin(a+b*log(c*x^n))^-3//2*cos(a+b*log(c*x^n)))
@test integrate(sin(a+-2*im*log(c*x))^-3//2, x) == :((1/2)*c^-4*x^-3*sin(a+-2*im*log(c*x))^-3//2*(1+-1*c^4*x^4*exp(2*im*a))*exp(-2*im*a))
@test integrate((e*x)^m*sin(d*(a+b*log(c*x^n)))^4, x) == :(e^-1*(e*x)^(1+m)*((1+m)^2+16*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^4*(1+m)+-4*b*d*n*e^-1*(e*x)^(1+m)*((1+m)^2+16*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^3*cos(d*(a+b*log(c*x^n)))+24*b^4*d^4*e^-1*n^4*(e*x)^(1+m)*(1+m)^-1*((1+m)^2+4*b^2*d^2*n^2)^-1*((1+m)^2+16*b^2*d^2*n^2)^-1+-24*b^3*d^3*e^-1*n^3*(e*x)^(1+m)*((1+m)^2+4*b^2*d^2*n^2)^-1*((1+m)^2+16*b^2*d^2*n^2)^-1*cos(d*(a+b*log(c*x^n)))*sin(d*(a+b*log(c*x^n)))+12*b^2*d^2*e^-1*n^2*(e*x)^(1+m)*((1+m)^2+4*b^2*d^2*n^2)^-1*((1+m)^2+16*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^2*(1+m))
@test integrate((e*x)^m*sin(d*(a+b*log(c*x^n)))^3, x) == :(e^-1*(e*x)^(1+m)*((1+m)^2+9*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^3*(1+m)+-6*b^3*d^3*e^-1*n^3*(e*x)^(1+m)*((1+m)^2+b^2*d^2*n^2)^-1*((1+m)^2+9*b^2*d^2*n^2)^-1*cos(d*(a+b*log(c*x^n)))+-3*b*d*n*e^-1*(e*x)^(1+m)*((1+m)^2+9*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^2*cos(d*(a+b*log(c*x^n)))+6*b^2*d^2*e^-1*n^2*(e*x)^(1+m)*((1+m)^2+b^2*d^2*n^2)^-1*((1+m)^2+9*b^2*d^2*n^2)^-1*(1+m)*sin(d*(a+b*log(c*x^n))))
@test integrate((e*x)^m*sin(d*(a+b*log(c*x^n)))^2, x) == :(e^-1*(e*x)^(1+m)*((1+m)^2+4*b^2*d^2*n^2)^-1*sin(d*(a+b*log(c*x^n)))^2*(1+m)+2*b^2*d^2*e^-1*n^2*(e*x)^(1+m)*(1+m)^-1*((1+m)^2+4*b^2*d^2*n^2)^-1+-2*b*d*n*e^-1*(e*x)^(1+m)*((1+m)^2+4*b^2*d^2*n^2)^-1*cos(d*(a+b*log(c*x^n)))*sin(d*(a+b*log(c*x^n))))
@test integrate((e*x)^m*sin(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*((1+m)^2+b^2*d^2*n^2)^-1*(1+m)*sin(d*(a+b*log(c*x^n)))+-1*b*d*n*e^-1*(e*x)^(1+m)*((1+m)^2+b^2*d^2*n^2)^-1*cos(d*(a+b*log(c*x^n))))
@test integrate(x^2*cos(a+b*log(c*x^n)), x) == :(3*x^3*(9+b^2*n^2)^-1*cos(a+b*log(c*x^n))+b*n*x^3*(9+b^2*n^2)^-1*sin(a+b*log(c*x^n)))
@test integrate(x*cos(a+b*log(c*x^n)), x) == :(2*x^2*(4+b^2*n^2)^-1*cos(a+b*log(c*x^n))+b*n*x^2*(4+b^2*n^2)^-1*sin(a+b*log(c*x^n)))
@test integrate(cos(a+b*log(c*x^n)), x) == :(x*(1+b^2*n^2)^-1*cos(a+b*log(c*x^n))+b*n*x*(1+b^2*n^2)^-1*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n)), x) == :(b^-1*n^-1*sin(a+b*log(c*x^n)))
@test integrate(x^-2*cos(a+b*log(c*x^n)), x) == :(-1*x^-1*(1+b^2*n^2)^-1*cos(a+b*log(c*x^n))+b*n*x^-1*(1+b^2*n^2)^-1*sin(a+b*log(c*x^n)))
@test integrate(x^2*cos(a+b*log(c*x^n))^2, x) == :(3*x^3*(9+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2+2*b^2*n^2*x^3*(27+12*b^2*n^2)^-1+2*b*n*x^3*(9+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x*cos(a+b*log(c*x^n))^2, x) == :(x^2*(2+2*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2+b^2*n^2*x^2*(4+4*b^2*n^2)^-1+b*n*x^2*(2+2*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(cos(a+b*log(c*x^n))^2, x) == :(x*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2+2*x*b^2*n^2*(1+4*b^2*n^2)^-1+2*b*n*x*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^2, x) == :((1/2)*log(x)+(1/2)*b^-1*n^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-2*cos(a+b*log(c*x^n))^2, x) == :(-1*x^-1*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2+-2*b^2*n^2*x^-1*(1+4*b^2*n^2)^-1+2*b*n*x^-1*(1+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^2*cos(a+b*log(c*x^n))^3, x) == :(x^3*(3+3*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3+2*b^2*n^2*x^3*(9+b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))+2*b^3*n^3*x^3*(27+3*b^4*n^4+30*b^2*n^2)^-1*sin(a+b*log(c*x^n))+b*n*x^3*(3+3*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*sin(a+b*log(c*x^n)))
@test integrate(x*cos(a+b*log(c*x^n))^3, x) == :(2*x^2*(4+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3+6*b^3*n^3*x^2*(16+9*b^4*n^4+40*b^2*n^2)^-1*sin(a+b*log(c*x^n))+12*b^2*n^2*x^2*(16+9*b^4*n^4+40*b^2*n^2)^-1*cos(a+b*log(c*x^n))+3*b*n*x^2*(4+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*sin(a+b*log(c*x^n)))
@test integrate(cos(a+b*log(c*x^n))^3, x) == :(x*(1+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3+6*x*b^2*n^2*(1+9*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))+6*x*b^3*n^3*(1+9*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))+3*b*n*x*(1+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^3, x) == :(b^-1*n^-1*sin(a+b*log(c*x^n))+-1//3*b^-1*n^-1*sin(a+b*log(c*x^n))^3)
@test integrate(x^-2*cos(a+b*log(c*x^n))^3, x) == :(-1*x^-1*(1+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3+-6*b^2*n^2*x^-1*(1+9*b^4*n^4+10*b^2*n^2)^-1*cos(a+b*log(c*x^n))+6*b^3*n^3*x^-1*(1+9*b^4*n^4+10*b^2*n^2)^-1*sin(a+b*log(c*x^n))+3*b*n*x^-1*(1+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*sin(a+b*log(c*x^n)))
@test integrate(cos(a+b*log(c*x^n))^4, x) == :(x*(1+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))^4+24*x*b^4*n^4*(1+20*b^2*n^2+64*b^4*n^4)^-1+12*x*b^2*n^2*(1+20*b^2*n^2+64*b^4*n^4)^-1*cos(a+b*log(c*x^n))^2+4*b*n*x*(1+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3*sin(a+b*log(c*x^n))+24*x*b^3*n^3*(1+20*b^2*n^2+64*b^4*n^4)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^4, x) == :(3//8*log(x)+1//4*b^-1*n^-1*cos(a+b*log(c*x^n))^3*sin(a+b*log(c*x^n))+3//8*b^-1*n^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(cos(log(6+3x)), x) == :((1+(1/2)*x)*cos(log(6+3x))+(1+(1/2)*x)*sin(log(6+3x)))
@test integrate(x^m*cos(a+(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n)), x) == :(x^(1+m)*(c*x^n)^(n^-1*(1+m))*(4+4m)^-1*exp(a*n^-1*(-1*n^-2*(1+m)^2)^-1//2*(1+m))+(1/2)*x^(1+m)*(c*x^n)^(-1*n^-1*(1+m))*exp(a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)*log(x))
@test integrate(cos(a+(-1*n^-2)^(1/2)*log(c*x^n)), x) == :(1//4*x*(c*x^n)^(n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+(1/2)*x*(c*x^n)^(-1*n^-1)*exp(a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^m*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^2, x) == :(x^(1+m)*(2+2m)^-1+x^(1+m)*(c*x^n)^(n^-1*(1+m))*(8+8m)^-1*exp(-2*a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)+1//4*x^(1+m)*(c*x^n)^(-1*n^-1*(1+m))*exp(2*a*n*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-1)*log(x))
@test integrate(cos(a+(1/2)*(-1*n^-2)^(1/2)*log(c*x^n))^2, x) == :((1/2)*x+1//8*x*(c*x^n)^(n^-1)*exp(-2*a*n*(-1*n^-2)^(1/2))+1//4*x*(c*x^n)^(-1*n^-1)*exp(2*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^m*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^3, x) == :(-4*x^(1+m)*(5+5m)^-1*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^3+8*x^(1+m)*(5+5m)^-1*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))+4//5*n*x^(1+m)*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-2*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))+-6//5*n*x^(1+m)*(-1*n^-2*(1+m)^2)^(1/2)*(1+m)^-2*cos(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n))^2*sin(a+(1/2)*(-1*n^-2*(1+m)^2)^(1/2)*log(c*x^n)))
@test integrate(cos(a+1//3*(-1*n^-2)^(1/2)*log(c*x^n))^3, x) == :(1//16*x*(c*x^n)^(n^-1)*exp(-3*a*n*(-1*n^-2)^(1/2))+9//16*x*(c*x^n)^(-1//3*n^-1)*exp(a*n*(-1*n^-2)^(1/2))+9//32*x*(c*x^n)^(1//3*n^-1)*exp(-1*a*n*(-1*n^-2)^(1/2))+1//8*x*(c*x^n)^(-1*n^-1)*exp(3*a*n*(-1*n^-2)^(1/2))*log(x))
@test integrate(x^-1*cos(a+b*log(c*x^n))^(1/2), x) == :(2*b^-1*n^-1*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*cos(a+b*log(c*x^n))^3//2, x) == :(2//3*b^-1*n^-1*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2)+2//3*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^5//2, x) == :(6//5*b^-1*n^-1*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2)+2//5*b^-1*n^-1*cos(a+b*log(c*x^n))^3//2*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^-1//2, x) == :(2*b^-1*n^-1*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*cos(a+b*log(c*x^n))^-3//2, x) == :(-2*b^-1*n^-1*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2)+2*b^-1*n^-1*cos(a+b*log(c*x^n))^-1//2*sin(a+b*log(c*x^n)))
@test integrate(x^-1*cos(a+b*log(c*x^n))^-5//2, x) == :(2//3*b^-1*n^-1*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2)+2//3*b^-1*n^-1*cos(a+b*log(c*x^n))^-3//2*sin(a+b*log(c*x^n)))
@test integrate(cos(a+-2*im*log(c*x))^-3//2, x) == :(-1//2*c^-4*x^-3*cos(a+-2*im*log(c*x))^-3//2*(1+c^4*x^4*exp(2*im*a))*exp(-2*im*a))
@test integrate(x^m*cos(a+b*log(c*x^n))^4, x) == :(x^(1+m)*((1+m)^2+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))^4*(1+m)+4*b*n*x^(1+m)*((1+m)^2+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3*sin(a+b*log(c*x^n))+24*b^4*n^4*x^(1+m)*(1+m)^-1*((1+m)^2+4*b^2*n^2)^-1*((1+m)^2+16*b^2*n^2)^-1+12*b^2*n^2*x^(1+m)*((1+m)^2+4*b^2*n^2)^-1*((1+m)^2+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*(1+m)+24*b^3*n^3*x^(1+m)*((1+m)^2+4*b^2*n^2)^-1*((1+m)^2+16*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^m*cos(a+b*log(c*x^n))^3, x) == :(x^(1+m)*((1+m)^2+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^3*(1+m)+3*b*n*x^(1+m)*((1+m)^2+9*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*sin(a+b*log(c*x^n))+6*b^3*n^3*x^(1+m)*((1+m)^2+b^2*n^2)^-1*((1+m)^2+9*b^2*n^2)^-1*sin(a+b*log(c*x^n))+6*b^2*n^2*x^(1+m)*((1+m)^2+b^2*n^2)^-1*((1+m)^2+9*b^2*n^2)^-1*(1+m)*cos(a+b*log(c*x^n)))
@test integrate(x^m*cos(a+b*log(c*x^n))^2, x) == :(x^(1+m)*((1+m)^2+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))^2*(1+m)+2*b^2*n^2*x^(1+m)*(1+m)^-1*((1+m)^2+4*b^2*n^2)^-1+2*b*n*x^(1+m)*((1+m)^2+4*b^2*n^2)^-1*cos(a+b*log(c*x^n))*sin(a+b*log(c*x^n)))
@test integrate(x^m*cos(a+b*log(c*x^n)), x) == :(x^(1+m)*((1+m)^2+b^2*n^2)^-1*(1+m)*cos(a+b*log(c*x^n))+b*n*x^(1+m)*((1+m)^2+b^2*n^2)^-1*sin(a+b*log(c*x^n)))
@test integrate(x^-1*tan(a+b*log(c*x^n)), x) == :(-1*b^-1*n^-1*log(cos(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^2, x) == :(-1*log(x)+b^-1*n^-1*tan(a+b*log(c*x^n)))
@test integrate(x^-1*tan(a+b*log(c*x^n))^3, x) == :(b^-1*n^-1*log(cos(a+b*log(c*x^n)))+(1/2)*b^-1*n^-1*tan(a+b*log(c*x^n))^2)
@test integrate(x^-1*tan(a+b*log(c*x^n))^4, x) == :(-1*b^-1*n^-1*tan(a+b*log(c*x^n))+1//3*b^-1*n^-1*tan(a+b*log(c*x^n))^3+log(x))
@test integrate(x^-1*tan(a+b*log(c*x^n))^5, x) == :(-1*b^-1*n^-1*log(cos(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*tan(a+b*log(c*x^n))^2+1//4*b^-1*n^-1*tan(a+b*log(c*x^n))^4)
@test integrate(x^-1*tan(a+b*log(c*x^n))^5//2, x) == :(2//3*b^-1*n^-1*tan(a+b*log(c*x^n))^3//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^3//2, x) == :(2*b^-1*n^-1*tan(a+b*log(c*x^n))^(1/2)+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^(1/2), x) == :((1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^-1//2, x) == :((1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^-3//2, x) == :(-2*b^-1*n^-1*tan(a+b*log(c*x^n))^-1//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*tan(a+b*log(c*x^n))^-5//2, x) == :(-2//3*b^-1*n^-1*tan(a+b*log(c*x^n))^-3//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*tan(a+b*log(c*x^n))^(1/2)+tan(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n)), x) == :(b^-1*n^-1*log(sin(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^2, x) == :(-1*log(x)+-1*b^-1*n^-1*cot(a+b*log(c*x^n)))
@test integrate(x^-1*cot(a+b*log(c*x^n))^3, x) == :(-1*b^-1*n^-1*log(sin(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*cot(a+b*log(c*x^n))^2)
@test integrate(x^-1*cot(a+b*log(c*x^n))^4, x) == :(b^-1*n^-1*cot(a+b*log(c*x^n))+-1//3*b^-1*n^-1*cot(a+b*log(c*x^n))^3+log(x))
@test integrate(x^-1*cot(a+b*log(c*x^n))^5, x) == :(b^-1*n^-1*log(sin(a+b*log(c*x^n)))+(1/2)*b^-1*n^-1*cot(a+b*log(c*x^n))^2+-1//4*b^-1*n^-1*cot(a+b*log(c*x^n))^4)
@test integrate(x^-1*cot(a+b*log(c*x^n))^5//2, x) == :(-2//3*b^-1*n^-1*cot(a+b*log(c*x^n))^3//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^3//2, x) == :(-2*b^-1*n^-1*cot(a+b*log(c*x^n))^(1/2)+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^(1/2), x) == :((1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^-1//2, x) == :((1/2)*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^-3//2, x) == :(2*b^-1*n^-1*cot(a+b*log(c*x^n))^-1//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*cot(a+b*log(c*x^n))^-5//2, x) == :(2//3*b^-1*n^-1*cot(a+b*log(c*x^n))^-3//2+(1/2)*2^(1/2)*b^-1*n^-1*arctan(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//2*2^(1/2)*b^-1*n^-1*arctan(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2))+-1//4*2^(1/2)*b^-1*n^-1*log(1+-1*2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n)))+1//4*2^(1/2)*b^-1*n^-1*log(1+2^(1/2)*cot(a+b*log(c*x^n))^(1/2)+cot(a+b*log(c*x^n))))
@test integrate(x^-1*sec(a+b*log(c*x^n)), x) == :(b^-1*n^-1*arctanh(sin(a+b*log(c*x^n))))
@test integrate(x^-1*sec(a+b*log(c*x^n))^2, x) == :(b^-1*n^-1*tan(a+b*log(c*x^n)))
@test integrate(x^-1*sec(a+b*log(c*x^n))^3, x) == :((1/2)*b^-1*n^-1*arctanh(sin(a+b*log(c*x^n)))+(1/2)*b^-1*n^-1*sec(a+b*log(c*x^n))*tan(a+b*log(c*x^n)))
@test integrate(x^-1*sec(a+b*log(c*x^n))^4, x) == :(b^-1*n^-1*tan(a+b*log(c*x^n))+1//3*b^-1*n^-1*tan(a+b*log(c*x^n))^3)
@test integrate(-1*(1+b^2*n^2)*sec(a+b*log(c*x^n))+2*b^2*n^2*sec(a+b*log(c*x^n))^3, x) == :(-1*x*sec(a+b*log(c*x^n))+b*n*x*sec(a+b*log(c*x^n))*tan(a+b*log(c*x^n)))
@test integrate(x^m*sec(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))^3, x) == :(x^(1+m)*(2+2m)^-1*sec(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))+(1/2)*x^(1+m)*(-1*(1+m)^2)^-1//2*sec(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))*tan(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2)))))
@test integrate(x*sec(a+2*log(c*x^im))^3, x) == :(x^2*(c*x^im)^(2im)*(1+(c*x^im)^(4im)*exp(2*im*a))^-2*exp(im*a))
@test integrate(sec(a+2*log(c*x^(-1//2*im)))^3, x) == :(2*x*(c*x^(-1//2*im))^(6im)*(1+(c*x^(-1//2*im))^(4im)*exp(2*im*a))^-2*exp(3*im*a))
@test integrate(sec(a+im*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*(2+-2p)^-1*sec(a+-1*im*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+(c*x^n)^(2*n^-1*(2+-1p)^-1)*exp(2*im*a))*(2+-1p)*exp(-2*im*a))
@test integrate(sec(a+-1*im*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(2+-2p)^-1*sec(a+im*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+(c*x^n)^(-2*n^-1*(2+-1p)^-1)*exp(2*im*a))*(2+-1p))
@test integrate(x^-1*sec(a+b*log(c*x^n))^(1/2), x) == :(2*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sec(a+b*log(c*x^n))^3//2, x) == :(2*b^-1*n^-1*sec(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))+-2*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sec(a+b*log(c*x^n))^5//2, x) == :(2//3*b^-1*n^-1*sec(a+b*log(c*x^n))^3//2*sin(a+b*log(c*x^n))+2//3*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sec(a+b*log(c*x^n))^-1//2, x) == :(2*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sec(a+b*log(c*x^n))^-3//2, x) == :(2//3*b^-1*n^-1*sec(a+b*log(c*x^n))^-1//2*sin(a+b*log(c*x^n))+2//3*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*sec(a+b*log(c*x^n))^-5//2, x) == :(2//5*b^-1*n^-1*sec(a+b*log(c*x^n))^-3//2*sin(a+b*log(c*x^n))+6//5*b^-1*n^-1*cos(a+b*log(c*x^n))^(1/2)*sec(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n)), x) == :(-1*b^-1*n^-1*arctanh(cos(a+b*log(c*x^n))))
@test integrate(x^-1*csc(a+b*log(c*x^n))^2, x) == :(-1*b^-1*n^-1*cot(a+b*log(c*x^n)))
@test integrate(x^-1*csc(a+b*log(c*x^n))^3, x) == :(-1//2*b^-1*n^-1*arctanh(cos(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*cot(a+b*log(c*x^n))*csc(a+b*log(c*x^n)))
@test integrate(x^-1*csc(a+b*log(c*x^n))^4, x) == :(-1*b^-1*n^-1*cot(a+b*log(c*x^n))+-1//3*b^-1*n^-1*cot(a+b*log(c*x^n))^3)
@test integrate(-1*(1+b^2*n^2)*csc(a+b*log(c*x^n))+2*b^2*n^2*csc(a+b*log(c*x^n))^3, x) == :(-1*x*csc(a+b*log(c*x^n))+-1*b*n*x*cot(a+b*log(c*x^n))*csc(a+b*log(c*x^n)))
@test integrate(x^m*csc(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))^3, x) == :(x^(1+m)*(2+2m)^-1*csc(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))+-1//2*x^(1+m)*(-1*(1+m)^2)^-1//2*cot(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2))))*csc(a+2*log(c*x^((1/2)*(-1*(1+m)^2)^(1/2)))))
@test integrate(x*csc(a+2*log(c*x^im))^3, x) == :(-1*im*x^2*(c*x^im)^(2im)*(1+-1*(c*x^im)^(4im)*exp(2*im*a))^-2*exp(im*a))
@test integrate(csc(a+2*log(c*x^(-1//2*im)))^3, x) == :(2*im*x*(c*x^(-1//2*im))^(6im)*(1+-1*(c*x^(-1//2*im))^(4im)*exp(2*im*a))^-2*exp(3*im*a))
@test integrate(csc(a+im*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(-1*x*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*(2+-2p)^-1*csc(a+-1*im*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+-1*(c*x^n)^(2*n^-1*(2+-1p)^-1)*exp(2*im*a))*(2+-1p)*exp(-2*im*a))
@test integrate(csc(a+-1*im*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(2+-2p)^-1*csc(a+im*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+-1*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*exp(2*im*a))*(2+-1p))
@test integrate(x^-1*csc(a+b*log(c*x^n))^(1/2), x) == :(2*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n))^3//2, x) == :(-2*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*cos(a+b*log(c*x^n))+-2*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n))^5//2, x) == :(-2//3*b^-1*n^-1*csc(a+b*log(c*x^n))^3//2*cos(a+b*log(c*x^n))+2//3*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n))^-1//2, x) == :(2*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n))^-3//2, x) == :(-2//3*b^-1*n^-1*csc(a+b*log(c*x^n))^-1//2*cos(a+b*log(c*x^n))+2//3*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))
@test integrate(x^-1*csc(a+b*log(c*x^n))^-5//2, x) == :(-2//5*b^-1*n^-1*csc(a+b*log(c*x^n))^-3//2*cos(a+b*log(c*x^n))+6//5*b^-1*n^-1*csc(a+b*log(c*x^n))^(1/2)*sin(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*a+-1//4*pi+(1/2)*b*log(c*x^n),2))

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/4.7.5 x^m trig(a+b log(c x^n))^p.jl


text:
export PreconditionerDataK2
"""
Abstract type that defines preconditioners for the K2 formulation.
The available preconditioners are:
- Identity
- Jacobi
"""
abstract type PreconditionerDataK2{T <: Real, S} end

include("identity.jl")
include("jacobi.jl")


text:
module HyperlineHairColorTrigger

using ..Ahorn, Maple

@mapdef Trigger "Hyperline/HairColorTrigger" HairChangeTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16, resetOnLeave::Bool=false, hairChanges::String="")

const placements = Ahorn.PlacementDict(
    "Hair Change Trigger (Hyperline)" => Ahorn.EntityPlacement(
        HairChangeTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: lordseanington/HyperlineWiki

max_issues_repo_name: lordseanington/HyperlineWiki

max_forks_repo_name: lordseanington/HyperlineWiki


text:
"""
    jackknife_error(g::Function, x; [binsize=10])

Computes the jackknife standard deviation of `g(<x>)` by binning
x and performing leave-one-out analysis.
"""
function jackknife_error(g::Function, x::AbstractVector{T}; binsize::Int=10) where T<:Real
    rem(length(x), binsize) == 0 || throw(ArgumentError("Bin size not compatible with sample size"))

    xblock = map(mean, Base.Iterators.partition(x, binsize))
    gis = Jackknife.leaveoneout(g, xblock)
    return sqrt(Jackknife.var(g, xblock; gis=gis))
end

function jackknife_error(g::Function, x::AbstractVector{T}; binsize::Int=10) where T<:Complex
    sqrt(jackknife_error(g, real(x), binsize=binsize)^2 + jackknife_error(g, imag(x), binsize=binsize)^2)
end

function jackknife_error(g::Function, x::AbstractArray{T}; binsize::Int=10) where T<:Number
    ndimsx = ndims(x)
    mapslices(y->binning_error(y; binsize=binsize), x, ndimsx)[(Colon() for _ in 1:ndimsx-1)...,1]
end

# Data: arrays
function jackknife_error(g::Function, x::Union{AbstractArray{<:Number},AbstractVector{T} where T<:AbstractArray{<:Number}}; binsize::Int=10)
    error("Jackknife error calculation for array like data not yet supported.")
end

"""
**Jackknife** errors for functions of uncertain data, i.e. g(<x>)
"""
module Jackknife

# https://github.com/ararslan/Jackknife.jl was useful resource

"""
    leaveoneout(g::Function, x)
Estimate the `g(x)` for each subsample of `x`, systematically
omitting each index one at a time. The result is a vector with the resulting `g(x_[])`.
"""
function leaveoneout(g::Function, x::AbstractVector{T}) where T<:Real
    length(x) > 1 || throw(ArgumentError("The sample must have size > 1"))
    return [g(circshift(x, -i)[2:end]) for i in 0:length(x)-1]
end
# function leaveoneout(g::Function, x::AbstractVector{T}) where T<:Real
#     length(x) > 1 || throw(ArgumentError("The sample must have size > 1"))
#     inds = eachindex(x)
#     return [g(mean(x[filter(j -> j != i, inds)])) for i in inds]
# end

"""
    var(g::Function, x)
Compute the jackknife estimate of the variance of `g(<x>)`, where `g` is given as a
function that computes a point estimate when passed a real-valued vector `x`.
"""
function var(g::Function, x::AbstractVector{T}; gis::AbstractVector{T}=leaveoneout(g, x)) where T<:Real
    n = length(x)
    return Base.var(gis) * (n - 1)^2 / n # Eq. (3.35) in QMC Methods book
end

"""
    bias(g::Function, x)
Compute the jackknife estimate of the bias of `g`, which is given as a
function that computes a point estimate when passed a real-valued vector `x`.
"""
function bias(g::Function, x::AbstractVector{T}; gis::AbstractVector{T}=leaveoneout(g, x)) where T<:Real
    return (length(x) - 1) * (mean(gis) - g(x)) # Basically Eq. (3.33)
end

"""
    estimate(g::Function, x)
Compute the bias-corrected jackknife estimate of the parameter `g(<x>)`.
"""
function estimate(g::Function, x::AbstractVector{T}; gis::AbstractVector{T}=leaveoneout(g, x)) where T<:Real
    n = length(x)
    return n * g(x) - (n - 1) * mean(gis) # Eq. (3.34) in QMC Methods book
end

end # module

export Jackknife


text:
module AttemptAtQNM

include("SpectralSolver.jl")
include("NewtonSolver.jl")
include("SchwarszchildModes.jl")
include("Interface.jl")

using .Interface

# Write your package code here.

struct Potato
    Root::Float64
end

print(GetModes(2,2,2,2))

export Potato
export GetModes

end

max_stars_repo_name: Potatoasad/AttemptAtQNM

max_issues_repo_name: Potatoasad/AttemptAtQNM

max_forks_repo_name: Potatoasad/AttemptAtQNM


text:
module TiledViews
using IndexFunArrays # needed for the default window function.
using NDTools # for linear_index

# using NDTools
export TiledView, get_num_tiles, TiledWindowView, tile_centers, get_window, tiled_processing
export get_num_tiles, eachtile, eachtilenumber, eachtilerelpos

tuple_len(::NTuple{N, Any}) where {N} = Val{N}()

 # T refers to the result type. N to the dimensions of the final array, and M to the dimensions of the raw array
struct TiledView{T, N, M, AA<:AbstractArray{T, M}} <: AbstractArray{T, N}
    # stores the data.
    parent::AA
    # output size of the array
    tile_size::NTuple{M, Int}
    tile_period::NTuple{M, Int}
    tile_offset::NTuple{M, Int}  # the distance from the wrapping arrray to the start of stored data
    pad_value::T

    # Constructor function
    function TiledView{T, N, M}(data::AA; tile_size::NTuple{M,Int}, tile_period::NTuple{M,Int}, tile_offset::NTuple{M,Int}, pad_value=nothing) where {T,M,N,AA}
        if isnothing(pad_value)
        if T <: NTuple
            pad_value = T(Base.Iterators.repeated(0))
        # elseif T <: AbstractArray
        else
            pad_value = convert(T,0)  # this may crash, but then the user should specify a valid pad_value
        end
        end
        return new{T, N, M, AA}(data, tile_size, tile_period, tile_offset, pad_value)
    end
end

function center(data)
    return size(data) .÷2 .+1
end

"""
    TiledView(data::F, tile_size::NTuple{N,Int}, tile_overlap::NTuple{N,Int}, tile_center::NTuple{M,Int} = (mod.(tile_size,2) .+1); pad_value::T, keep_center=true)

Creates an 2N dimensional view of the data by tiling the N-dimensional data as
specified by tile_size, tile_overlap and optionally tile_center.

`data`. the inputdata to decompose into a TiledView. No copies are made for the TiledView and
the raw data can be accessed via myview.parent.

`tile_size`. A Tuple describing the size of each tile. This size will form the first N dimensions of the
result of size(myview). The second N dimensions refer to N-dimensional tile numbering.

`tile_overlap`. Tuple specifying the overlap between successive tiles in voxels. This implicitely
defines the pitch between tiles as (tile_size .- tile_overlap).

`pad_value`. Specifies the answer that is returned when get_index is applied to a position outside the source array.

`keep_center`.  This boolean specifies whether the center of the parant `data` will be aligned with the center of the central tile. If `false`, the first tile starts at offset zero.
`tile_center`.  Only used if `keep_center` is true. It defines the center position in the central tile. The default is `tile_size .÷ 2 .+1`.

# Examples
```jldoctest
julia> a = TiledView(reshape(1:49,(7,7)), (4, 4),(1, 1));
julia> a.parent
7×7 reshape(::UnitRange{Int64}, 7, 7) with eltype Int64:
 1   8  15  22  29  36  43
 2   9  16  23  30  37  44
 3  10  17  24  31  38  45
 4  11  18  25  32  39  46
 5  12  19  26  33  40  47
 6  13  20  27  34  41  48
 7  14  21  28  35  42  49
 julia> size(a)
(4, 4, 3, 3)
```
"""
function TiledView(data::AbstractArray{T,M}, tile_size::NTuple{M,Int}, tile_overlap::NTuple{M,Int}=tile_size .* 0,
                   tile_center::NTuple{M,Int} = (tile_size .÷ 2 .+1); pad_value=nothing, keep_center=true) where {T, M}
    # Note that N refers to the original number of dimensions
    tile_period = tile_size .- tile_overlap
    if keep_center
        data_center = center(data)
        tile_offset = mod.((data_center .- tile_center), tile_period)
    else
        tile_offset = tile_period .* 0
    end
    N = 2*M
    return TiledView{T,N,M}(data; tile_size=tile_size, tile_period=tile_period, tile_offset=tile_offset, pad_value=pad_value)
end

function get_num_tiles(data::TiledView)
    num_tiles = ((size(data.parent) .+ data.tile_offset .- 1) .÷ data.tile_period) .+ 1
    return num_tiles
end

# define AbstractArray function to allow to treat the generator as an array
# See https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array
function Base.size(A::TiledView)
    return (A.tile_size...,(get_num_tiles(A))...)
end

function zeros_like(A::TiledView, ::Type{T}=eltype(A.parent)) where {T}
    TiledView{T,ndims(A),length(A.tile_size)}(zeros(T, size(A.parent)); tile_size=A.tile_size, tile_period=A.tile_period, tile_offset=A.tile_offset, pad_value=A.pad_value)
end

function ones_like(A::TiledView, ::Type{T}=eltype(A.parent)) where {T}
    TiledView{T,ndims(A),length(A.tile_size)}(ones(T, size(A.parent)); tile_size=A.tile_size, tile_period=A.tile_period, tile_offset=A.tile_offset, pad_value=A.pad_value)
end

# Note that the similar function below will most of the times expand the overall required datasize
function Base.similar(A::TiledView, ::Type{T}=eltype(A.parent), dims::Dims=size(A)) where {T}
    # The first N coordinates are position within a tile, the second N coordinates are tile number
    #= N = length(A.tile_size)
    new_tile_sz = dims[1:N]
    new_num_tiles = dims[N + 1:end]
    new_tile_period = A.tile_period .+ new_tile_sz .- A.tile_size
    new_core_sz = new_num_tiles .* new_tile_period .- A.tile_offset  # keep the overlap the same as before
    TiledView{T,ndims(A),length(A.tile_size)}(similar(A.parent, T, new_core_sz); tile_size=new_tile_sz, tile_period=new_tile_period, tile_offset=A.tile_offset, pad_value=A.pad_value)
    =#
    similar(A.parent, T, dims) # this returns an ordinary array, but this seems the only way to handle cases like: my_tiled_array[:,:,2,3] which expect a 2D array
end
# Array{eltype(A)}(undef, size...)

# %24 = Base.getproperty(A, :parent)::AbstractMatrix{Float64}

# calculate the entry according to the index
# Base.getindex(A::IndexFunArray{T,N}, I::Vararg{B, N}) where {T,N, B} = return A.generator(I)

function pos_from_tile(A::TiledView{T,N}, TilePos::NTuple{M,Int}, TileNum::NTuple{M,Int}) where {T,N,M}
    Tuple(TilePos[n] - A.tile_offset[n] + (TileNum[n]-1) * A.tile_period[n]  for n in 1:M)
end

# calculate the entry according to the index
Base.@propagate_inbounds function Base.getindex(A::TiledView{T,N,M,AA}, I::Vararg{Int, N})::T where {T,N,M,AA}
    @boundscheck checkbounds(A, I...)
    @inbounds pos = (I[n] - A.tile_offset[n] + (I[n+M].-1) * A.tile_period[n]  for n in 1:M)
    if Base.checkbounds(Bool, A.parent, pos...)
        return Base.getindex(A.parent, pos...)::T
    else
        return A.pad_value :: T;
    end
end

Base.setindex!(A::TiledView{T,N,M,AA}, v, I::Vararg{Int,N}) where {T,N,M,AA} = begin
    @boundscheck checkbounds(A, I...)
    @inbounds pos = (I[n] - A.tile_offset[n] + (I[n+M].-1) * A.tile_period[n]  for n in 1:M)
    # pos = TilePos .- A.tile_offset .+ (TileNum.-1) .* A.tile_period
    if Base.checkbounds(Bool, A.parent, pos...)
        return setindex!(A.parent, v, pos... )
    else
        return convert(T,0)
    end
end

## Some functions for generating useful tilings

"""
    get_window(A::TiledView; window_function=window_hanning, get_norm=false, verbose=false, offset = CtrFT);

Calculates a window matching to the `TiledView`.

`window_function`. The window function as defined in IndexFunArrays to apply to the TiledView.
The result is currently not any longer a view as it is unclear how to wrap the multiplication into a view.
For this reason the TiledView without the window applied is also returned and can be used for assignments.
By default a von Hann window (window_hanning) is used. For even sizes the window is centered at the integer coordinate right of the middle position (`CtrFT`).

`get_norm`. An optional Boolean argument allowing to also obtain the normalization map for the boarder pixels, which
not necessarily undergo all required window operations. In a future version it may be possible
to automatically lay out the windowing such that this effect can be avoided.

`verbose`. If true, diagnostic information on the window layout is printed.

`offset`. defines where the center of the window is placed. See `IndexFunArrays.jl` for details.

# Returns
`matching_window`. a window that can be applied to the view via multiplication myview.*matching_window
This is intentionally not provided as a product to separate the features conceptually
when it comes to write access.

`normalized`. only returned for get_norm=true. Contains an array with the normalization information by
mapping the window back to the original data. This is useful for incomplete coverage of the tiles
as well as using windows which do not sum up to one in the tiling process.

# Examples
```jldoctest
julia> data = ones(10,10).+0.0;
julia> myview = TiledView(data, (5, 5), (2,2));
julia> win = get_window(myview, verbose=true);
Tiles with pitch (3, 3) overlap by (2, 2) pixels.
Window starts at (0.5, 0.5) and ends at (2.5, 2.5).

julia> win
5×5 IndexFunArrays.IndexFunArray{Float64, 2, IndexFunArrays.var"#329#331"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.0214466  0.125     0.146447  0.125     0.0214466
 0.125      0.728553  0.853553  0.728553  0.125
 0.146447   0.853553  1.0       0.853553  0.146447
 0.125      0.728553  0.853553  0.728553  0.125
 0.0214466  0.125     0.146447  0.125     0.0214466

 see TiledWindowView() for more examples.
"""
function get_window(A::TiledView; window_function=window_hanning, get_norm=false, verbose=false, offset=CtrFT)
    tile_size = A.tile_size
    tile_pitch = A.tile_period
    tile_overlap = tile_size .- tile_pitch

    winend = (tile_size ./ 2.0)
    winstart = (winend .- tile_overlap)
    if verbose
        print("Tiles with pitch $tile_pitch overlap by $tile_overlap pixels.\n")
        print("Window starts at $winstart and ends at $winend.\n")
    end
    if get_norm == false
        return window_function(tile_size;
            scale=ScaUnit, offset=offset,
            border_in=winstart, border_out= winend)
    else
        my_view = ones_like(A)
        normalization = A.parent
        normalization .= 0
        my_view .+= A .*window_function(tile_size;scale=ScaUnit, offset=offset, border_in=winstart, border_out= winend)
        return (window_function(tile_size;
            scale=ScaUnit, offset=offset,
            border_in=winstart, border_out= winend), normalization)
    end
end

"""
function TiledWindowView(data::AbstractArray{T,M}, tile_size::NTuple{M,Int};
    rel_overlap::NTuple{M,Float64}=tile_size .*0 .+ 1.0,
    window_function=window_hanning, get_norm=false, verbose=false, offset=CtrFT) where {T, M}

Creates an 2N dimensional view of the data by tiling the N-dimensional data as
specified by tile_size, tile_overlap and optionally tile_center.
Additionally a window is applied to this view. If the window_type as defined in
IndexFunArrays sums up to one,
which is the case for window_linear and window_hanning, a linear decomposition of the
data is obtained apart from possible border effects.
`data`. the inputdata to decompose into a TiledView. No copies are made for the TiledView and
the raw data can be accessed via myview.parent.

`tile_size`. A Tuple describing the size of each tile. This size will form the first N dimensions of the
result of size(myview). The second N dimensions refer to N-dimensional tile numbering.

`rel_overlap`. Tuple specifying the relative overlap between successive tiles. The absolute overlap is then calculated as `round.(Int,tile_size./2.0 .* rel_overlap)`.

`window_function`. The window function as defined in IndexFunArrays to apply to the TiledView.
The result is currently not any longer a view as it is unclear how to wrap the multiplication into a view.
For this reason the TiledView without the window applied is also returned and can be used for assignments.
By default a von Hann window (window_hanning) is used.

`get_norm`. An optional Boolean argument allowing to also obtain the normalization map for the boarder pixels, which
not necessarily undergo all required window operations. In a future version it may be possible
to automatically lay out the windowing such that this effect can be avoided.

`verbose`. If true, diagnostic information on the window layout is printed.

`offset`. defines where the center of the window is placed. See `IndexFunArrays.jl` for details.

# Returns
myview, matching_window = TiledWindowView ...
a Tuple of two or three (get_norm=true) with

`myview`. the TiledView of the data without the window which can also be written to.

`matching_window`. a window that can be applied to the view via multiplication myview.*matching_window
This is intentionally not provided as a product to separate the features conceptually
when it comes to write access.

`normalized`. only returned for get_norm=true. Contains an array with the normalization information by
mapping the window back to the original data. This is useful for incomplete coverage of the tiles
as well as using windows which do not sum up to one in the tiling process.

Note that it may be dangerous to directly access the view via a simple .+= operation as
it is not entirely clear, whether it is always garanteed that there could not be any running
conditions with read-write operations, since some points in the referenced array are accessed
multiple times. To avoid such an effect, you can, for example, only acess every second tile along each dimension
in one call.
# Examples
```jldoctest
julia> data = ones(10,10).+0.0;
julia> myview, matching_window = TiledWindowView(data, (5, 5);verbose=true);
Tiles with pitch (3, 3) overlap by (2, 2) pixels.
Window starts at (0.5, 0.5) and ends at (2.5, 2.5).
julia> size(myview)
(5, 5, 4, 4)
julia> matching_window
5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#199#200"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.0214466  0.125     0.146447  0.125     0.0214466
 0.125      0.728553  0.853553  0.728553  0.125
 0.146447   0.853553  1.0       0.853553  0.146447
 0.125      0.728553  0.853553  0.728553  0.125
 0.0214466  0.125     0.146447  0.125     0.0214466
julia> windowed = collect(myview .* matching_window);
julia> myview[:,:,:,:].=0  # cleares the original array
julia> myview.parent
10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
julia> myview .+= windowed  # writes the windowed data back into the array
julia> data # lets see if the weigths correctly sum to one?
10×10 Matrix{Float64}:
 0.728553  0.853553  0.853553  0.853553  0.853553  0.853553  0.853553  0.853553  0.853553  0.853553
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
 0.853553  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
```
# This result may also be used for subsequent normalization but can also be directly obtained by
julia> myview, matching_window, normalized = TiledWindowView(rand(10,10).+0, (5, 5);get_norm=true);
"""
function TiledWindowView(data::AbstractArray{T,M}, tile_size::NTuple{M,Int};
                         rel_overlap::NTuple{M,Float64}=tile_size .*0 .+ 1.0,
                         window_function=window_hanning, get_norm=false, verbose=false, keep_center=true, offset=CtrFT) where {T, M}
    tile_overlap = round.(Int,tile_size./2.0 .* rel_overlap)
    changeable = TiledView(data,tile_size, tile_overlap, keep_center=keep_center);
    win = get_window(changeable, window_function= window_function, get_norm=get_norm, verbose=verbose, offset=offset)
    if get_norm
        return (changeable, win...)
    else
        return (changeable, win)
    end
end

"""
    tile_centers(A, scale=nothing)
    returns the relative center coordinates of integer tile centers with respect to the integer center `1 .+ size(A) .÷ 2 `
    The tuple `scale` is used to multiply the relative position with a physical pixelsize.
    See also: `eachtilerelpos` for a corresponding iterator
"""
function tile_centers(A, scale=nothing)
    return collect(eachtilerelpos(A, scale))
end

"""
    eachtilerelpos(A, scale=nothing)
    returns a generator that iterates through the relative distance of each tile center `1 .+ size(A).÷2` to the center of the the untiled parent array `1 .+ size(parent).÷2`
    The tuple `scale` is used to multiply the relative position with a physical pixelsize.
"""
function eachtilerelpos(A, scale=nothing)
    nd = ndims(A)/2
    ctr_array = (size(A.parent) .÷ 2) .+ 1 # center of the parent array
    num_tiles = size(A)[end-nd+1:end]
    # ctr_tile = (num_tiles.÷2 .+1)
    tile_ctr = (size(A)[1:nd] .÷ 2) .+ 1
    if isnothing(scale)
        (pos_from_tile(A, tile_ctr, Tuple(idx)) .- ctr_array for idx in CartesianIndices(num_tiles))  # Only the "[" generate a 2D array
        # [(Tuple(idx).-1).* A.tile_period .+ ctr for idx in CartesianIndices(num_tiles)]
    else
        (scale .* (pos_from_tile(A, tile_ctr, Tuple(idx)) .- ctr_array) for idx in CartesianIndices(num_tiles))
    end
end

"""
    eachtile(tiled_view::TiledView)
    returns an iterator which iterates through all tiles. Depending on your application you may also want to use
    `tiled_processing` for a convenient way to apply a function to each tile and join all tiles back together.
    If you need simultaneous access to the tiles and tile numbers, you can also use `eachtilenumber`.
"""
function eachtile(tiled_view::TiledView)
    nd = ndims(tiled_view)÷2
    nz = ((size(tiled_view)[1:nd])..., prod(size(tiled_view)[nd+1:end]))
    reshaped = reshape(tiled_view, nz)
    return eachslice(reshaped, dims=nd+1)
end

"""
    eachtilenumber(tiled_view::TiledView)
    returns an iterator iterating though all the tile numbers. If you need access to the tiles themselves, use
    `eachtile`
"""
function eachtilenumber(tiled_view::TiledView)
    return (Tuple(tn) for tn in CartesianIndices(get_num_tiles(tiled_view)))
end

"""
    tiled_processing(tiled_view::TiledView, fct; verbose=true, dtype=eltype(tiled_view.parent), window_function=window_hanning)
    processes a raw dataset using tiled views by submitting each tile to the function `fct` and merging the results via the `window_function`.
"""
function tiled_processing(tiled_view::TiledView, fct; verbose=true, dtype=eltype(tiled_view.parent), window_function=window_hanning)
    @time res = zeros_like(tiled_view, dtype)
    res.parent .= zero(dtype)
    @time win = get_window(tiled_view, window_function=window_function)
    ttn = get_num_tiles(tiled_view)
    for (src, dest, tn) in zip(eachtile(tiled_view), eachtile(res), eachtilenumber(res))
        if verbose
            perc = round(100 * (linear_index(tn, ttn)-1) ./ prod(ttn))
            print("processing tile $(tn) out of $(ttn), $(perc)%\n")
        end
        size(src)
        res_tile = fct(collect(src))
        dest .+= win .* res_tile
    end
    return res
end

"""
    tiled_processing(data, fct; verbose=true, dtype=eltype(data), window_function=window_hanning)
    processes a raw dataset using tiled views by submitting each tile to the function `fct` and merging the results via the `window_function`.
"""
function tiled_processing(data, fct, tile_size, tile_overlap; verbose=true, dtype=eltype(data), keep_center=false, window_function=window_hanning)
    tiles = TiledView(data, tile_size, tile_overlap, keep_center=keep_center);
    return tiled_processing(tiles,fct; verbose=verbose, dtype=dtype, window_function=window_function)
end

end # module


text:
function test_ag(amat::NamedArray; m=Symbol[], c=Symbol[], debug=false)

  function show_differences(res)
    if length(res) > 0
      debug && println()
      debug && display(res)
      debug && println()
    end
  end

  a = amat'
  vars = names(a, 1)
  zna = NamedArray(zeros(Int, size(a)), (vars, vars), ("Rows", "Cols"));

  s = StructuralCausalModels.update_s(a, c, debug)

  ar = copy(a);
  count = 0
  while true
    count += 1
    debug && println("\ncount = $count\n")
    debug && println()

    at = copy(ar)

    debug && println("\nupdate_21\n")
    ar21 = StructuralCausalModels.update_21(ar, m, debug);
    res = findall(ar21 .!== ar)
    debug && show_differences(res)

    ar = copy(ar21);
    debug && println("\nupdate_22\n")
    ar22 = StructuralCausalModels.update_22(ar, s, debug)
    res = findall(ar22 .!== zna)
    debug && show_differences(res)

    ar .+= ar22

    debug && println("\nupdate_23\n")
    ar23 = StructuralCausalModels.update_23(ar, m, debug)
    res = findall(ar23 .!== zna)
    debug && show_differences(res)

    ar .+= ar23

    debug && println("\nupdate_24\n")
    ar24 = StructuralCausalModels.update_24(ar, m, debug)
    res = findall(ar24 .!== zna)
    debug && show_differences(res)

    ar .+= ar24

    debug && println("\nupdate_25\n")
    ar25 = StructuralCausalModels.update_25(ar, m, debug)
    res = findall(ar25 .!== zna)
    debug && show_differences(res)

    ar .+= ar25

    debug && println("\nupdate_26\n")
    ar26 = StructuralCausalModels.update_26(ar, m, debug)
    res = findall(ar26 .!== zna)
    debug && show_differences(res)

    ar .+= ar26

    debug && println("\nupdate_27\n")
    ar27 = StructuralCausalModels.update_27(ar, s, debug)
    res = findall(ar27 .!== zna)
    debug && show_differences(res)

    ar .+= ar27

    debug && println("\nupdate_28\n")
    ar28 = StructuralCausalModels.update_28(ar, s, debug)
    res = findall(ar28 .!== zna)
    debug && show_differences(res)

    ar .+= ar28

    debug && println("\nupdate_29\n")
    ar29 = StructuralCausalModels.update_29(ar, m, debug)
    res = findall(ar29 .!== zna)
    debug && show_differences(res)

    ar .+= ar29

    debug && println("\nupdate_30\n")
    ar30 = StructuralCausalModels.update_30(ar, m, debug)
    res = findall(ar30 .!== zna)
    debug && show_differences(res)

    ar .+= ar30

    ar == at && break

  end

  debug && println("\nEnd of 1st while loop ($count iterations).\n")

  ar_old = copy(ar)
  debug && println("\nupdate_a\n")
  ar = StructuralCausalModels.update_a(ar, debug)
  res = findall(ar .!== ar_old)
  debug && show_differences(res)

  ar_old = copy(ar)
  debug && println("\nupdate_b\n")
  ar = StructuralCausalModels.update_b(ar, s, debug)
  res = findall(ar .!== ar_old)
  debug && show_differences(res)

  an = copy(ar)
  count = 0
  while true
    count += 1
    debug && println("\ncount = $count\n")
    debug && println()

    at = copy(ar)

    debug && println("\nupdate_27n\n")
    a27n = StructuralCausalModels.update_27n(ar, an)
    res = findall(a27n .!== zna)
    debug && show_differences(res)

    an .+= a27n
    ar .+= a27n

    debug && println("\nupdate_22n\n")
    a22n = StructuralCausalModels.update_22n(ar, an)
    res = findall(a22n .!== zna)
    debug && show_differences(res)

    an .+= a22n
    ar .+= a22n

    at == ar && break

  end

  debug && println("\nEnd of 2nd while loop ($count iterations).\n")

  debug && println("\nupdate_i\n")
  ai = StructuralCausalModels.update_i(ar)
  res = findall(ar .!== ai)
  debug && show_differences(res)

  ar = copy(ai)

  select = setdiff(vars, vcat(c, m))
  debug && println(select)
  fr = ar[select, select]

  #println("/ntest_ag() completed.")

  fr
end

function test_ag(d::DAG; m=Symbol[], c=Symbol[], debug=false)
  test_ag(d.e; m=m, c=c, debug=debug)
end


text:
# Automatically generated using Clang.jl

@cenum cutensorOperator_t::UInt32 begin
    CUTENSOR_OP_IDENTITY = 1
    CUTENSOR_OP_SQRT = 2
    CUTENSOR_OP_RELU = 8
    CUTENSOR_OP_CONJ = 9
    CUTENSOR_OP_RCP = 10
    CUTENSOR_OP_ADD = 3
    CUTENSOR_OP_MUL = 5
    CUTENSOR_OP_MAX = 6
    CUTENSOR_OP_MIN = 7
    CUTENSOR_OP_UNKNOWN = 126
end

@cenum cutensorStatus_t::UInt32 begin
    CUTENSOR_STATUS_SUCCESS = 0
    CUTENSOR_STATUS_NOT_INITIALIZED = 1
    CUTENSOR_STATUS_ALLOC_FAILED = 3
    CUTENSOR_STATUS_INVALID_VALUE = 7
    CUTENSOR_STATUS_ARCH_MISMATCH = 8
    CUTENSOR_STATUS_MAPPING_ERROR = 11
    CUTENSOR_STATUS_EXECUTION_FAILED = 13
    CUTENSOR_STATUS_INTERNAL_ERROR = 14
    CUTENSOR_STATUS_NOT_SUPPORTED = 15
    CUTENSOR_STATUS_LICENSE_ERROR = 16
    CUTENSOR_STATUS_CUBLAS_ERROR = 17
    CUTENSOR_STATUS_CUDA_ERROR = 18
    CUTENSOR_STATUS_INSUFFICIENT_WORKSPACE = 19
    CUTENSOR_STATUS_INSUFFICIENT_DRIVER = 20
end

@cenum cutensorAlgo_t::Int32 begin
    CUTENSOR_ALGO_TGETT = -7
    CUTENSOR_ALGO_GETT = -6
    CUTENSOR_ALGO_LOG_TENSOR_OP = -5
    CUTENSOR_ALGO_LOG = -4
    CUTENSOR_ALGO_TTGT_TENSOR_OP = -3
    CUTENSOR_ALGO_TTGT = -2
    CUTENSOR_ALGO_DEFAULT = -1
end

@cenum cutensorWorksizePreference_t::UInt32 begin
    CUTENSOR_WORKSPACE_MIN = 1
    CUTENSOR_WORKSPACE_RECOMMENDED = 2
    CUTENSOR_WORKSPACE_MAX = 3
end

const cutensorHandle_t = Ptr{Cvoid}
const cutensorTensorDescriptor_t = Ptr{Cvoid}
const CUTENSOR_MAJOR = 0
const CUTENSOR_MINOR = 2
const CUTENSOR_PATCH = 2


text:
@testset "InverseSolve" begin
    @testset "testModels.jl" begin
        include("testModels.jl")
    end
    @testset "testMisfit.jl" begin
        include("testMisfit.jl")
    end
    @testset "testRegularizers.jl" begin
       include("testRegularizers.jl")
    end
    @testset "testMisfitParams.jl" begin
       include("testMisfitParams.jl")
    end
    @testset "testLeastSquares.jl" begin
       include("testLeastSquares.jl")
    end
    @testset "testGetHessian.jl" begin
       include("testGetHessian.jl")
    end
    @testset "testHessMatVec.jl" begin
        include("testHessMatVec.jl")
    end
end