text:
# Calibration, smoothing, then jump to the OrdinaryDiffEq._postamble!
function OrdinaryDiffEq.postamble!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})
    OrdinaryDiffEq._postamble!(integ)
    # For some unknown reason, the following is necessary
    copyat_or_push!(integ.sol.k, integ.saveiter_dense, integ.k)

    # Add the final timepoint to the solution
    pn_solution_endpoint_match_cur_integrator!(integ)

    if isstatic(integ.cache.diffusionmodel) # Calibrate
        # @warn "sol.log_likelihood is not correct for static diffusion models!"
        integ.sol.log_likelihood = NaN
        final_diff = integ.sol.diffusions[end]
        @simd ivdep for s in integ.sol.x_filt
            # s.Σ .*= final_diff
            copy!(s.Σ, apply_diffusion(s.Σ, final_diff))
        end

        if isempty(size(final_diff))
            integ.sol.diffusions .= final_diff
        else
            [(d .= final_diff) for d in integ.sol.diffusions]
        end
    end

    if integ.alg.smooth
        smooth_all!(integ)
        for (pu, x) in zip(integ.sol.pu, integ.sol.x_smooth)
            mul!(pu, integ.cache.SolProj, x)
        end
        integ.sol.interp = set_smooth(integ.sol.interp)
        [(su[:] .= pu) for (su, pu) in zip(integ.sol.u, integ.sol.pu.μ)]
    end
    @assert (length(integ.sol.u) == length(integ.sol.pu))

    return nothing
end

function pn_solution_endpoint_match_cur_integrator!(integ)
    # Inspired from OrdinaryDiffEq.solution_match_cur_integrator!
    if integ.opts.save_end
        if integ.alg.smooth
            OrdinaryDiffEq.copyat_or_push!(
                integ.sol.x_filt,
                integ.saveiter_dense,
                integ.cache.x,
            )
        end

        OrdinaryDiffEq.copyat_or_push!(
            integ.sol.pu,
            integ.saveiter,
            mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x),
        )
    end
end

function DiffEqBase.savevalues!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    force_save=false,
    reduce_size=true,
)

    # Do whatever OrdinaryDiffEq would do
    out = OrdinaryDiffEq._savevalues!(integ, force_save, reduce_size)

    # Save our custom stuff that we need for the posterior
    # TODO If we don't want dense output, we might not want to save these!
    # It's not completely clear how to specify that though; They are also needed for sampling.
    if integ.alg.smooth
        OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, integ.saveiter, integ.cache.x)
    end
    OrdinaryDiffEq.copyat_or_push!(
        integ.sol.diffusions,
        integ.saveiter,
        integ.cache.global_diffusion,
    )
    if integ.opts.save_everystep
        OrdinaryDiffEq.copyat_or_push!(
            integ.sol.pu,
            integ.saveiter,
            mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x),
        )
    end

    return out
end


text:
"""Returns seed from parameters"""
getseed(params::Parameters) = params.random_seed == 0 ? time_ns() : params.random_seed

"""
Sample uniformly k integers from the integer range 1:n, making sure that
the same integer is not sampled twice. the function returns an integer vector
containing the sampled integers.
"""
function randpermset(n::Int, k::Int, rng::AbstractRNG)
    if n < k
        error("Sample size cannot be grater then range")
    end
    sample = Set{UInt}()
    sample_size = 0
    while sample_size < k
        selected = rand(rng, 1:n, k-sample_size)
        union!(sample, selected)
        sample_size = length(sample)
    end
    return sample
end

"""
Sample randomly lm_dim+1 points from the dataset, making sure
that the same point is not sampled twice.
The function will return a index vector, specifying the index of the sampled points.
"""
function sample_points(X::AbstractMatrix{T}, n::Int;
                       rng::AbstractRNG=GLOBAL_RNG,
                       resample::Int=3) where T <: Real
    if n <= 0
        error("Sample size must be positive")
    end
    data_rows, data_cols = size(X)
    idxs = Set{UInt}()
    hashes = Set{UInt}()
    hashes_size = 0
    resample_count = 0

    # get sample indexes
    while length(idxs) != n && resample_count < resample
        J = randpermset(data_cols, n, rng)
        for idx in setdiff(J, idxs)
            h = hash(view(X, :, idx))
            push!(hashes, h)
            hashes_size += 1
            # Find duplicate row in sample and delete it
            if length(hashes) != hashes_size
                hashes_size -= 1
            else
                push!(idxs, idx)
            end
        end
        resample_count += 1
        @logmsg(DEBUG_SAMPLING,
            "Sampling iteration $resample_count",
            sample_indexes = idxs,
            sampled =J
        )
    end

    # If at this point we do not have proper sample
    # then our dataset doesn't have enough unique rows
    if length(idxs) < n
        @logmsg(DEBUG_SAMPLING,
            "Dataset doesn't have enough unique points, decrease number of sampled points",
            sample_size = n,
            sampled = length(idxs)
        )
        return UInt[]
    end
    return collect(idxs)
end

"""Reservoir sampling"""
function reservoir_sample_points(X::AbstractMatrix{T}, k::Int, r::MersenneTwister) where T <: Real
    N, n = size(X)
    if n < k
        @warn "Not enough samples to construct manifold" maxlog=5
        return Int[]
    end

    idxs = collect(1:k)
    for i in (k+1):n
        j = trunc(Int, rand(r)*i)+1
        if j <= k
            idxs[j] = i
        end
    end

    for i in 1:(k-1)
        for j in (i+1):k
            if all([X[k,idxs[i]] == X[k,idxs[j]] for k in 1:N])
                @warn "Sample is not unique: X[:,$(idxs[i])] == X[:,$(idxs[j])]"
                return Int[]
            end
        end
    end

    return idxs
end

"""Modified Gram-Schmidt orthogonalization algorithm"""
function orthogonalize(vecs::AbstractMatrix{T}) where T<:Real
    m, n = size(vecs)
    basis = zeros(T, m, n)
    for j = 1:n
        v_j = vecs[:,j]
        for i = 1:(j-1)
            q_i = basis[:,i]
            r_ij = dot(q_i, v_j)
            v_j -= q_i*r_ij
        end
        r_jj = norm(v_j)
        basis[:,j] = r_jj != 0.0 ? v_j/r_jj : v_j
    end
    basis
end

# Forming basis from sample. The idea is to pick a point (origin) from the sampled points
# and generate the basis vectors by subtracting all other points from the origin,
# creating a basis matrix with one less vector than the number of sampled points.
# Then perform orthogonalization through Gram-Schmidt process.
# Note: Resulting basis is transposed.
function form_basis(X::AbstractMatrix{T}, sample_indices::Vector{S}) where {T<:AbstractFloat, S<:Integer}
    origin = X[:,sample_indices[1]]
    basis = X[:,sample_indices[2:end]] .- origin
    vec(origin), orthogonalize(basis)
end

function filter_separated(selected_points, X, O, B, S; ocss=false)
    cluster_points = Int[]
    removed_points  = Int[]
    θ = threshold(S)
    ii = ocss ? 2 : 1

    for i in eachindex(selected_points)
        idx = selected_points[i]
        # point i has distances less than the threshold value
        d = distance_to_manifold(X[:, idx], O, B)
        if d[ii] < θ
            push!(cluster_points, idx)
        else
            push!(removed_points, idx)
        end
    end

    return cluster_points, removed_points
end

function adjustbasis!(M::Manifold, X::AbstractMatrix{T};
                      adjust_dim::Bool=false, adjust_dim_ratio::Real=0.99) where {T<:AbstractFloat}
    pts = points(M)
    MX = view(X, :, pts)
    d,n = size(MX)
    μ = mean(MX, dims=2)
    Svd = svd(MX .- μ)
    v = Svd.S
    U = Svd.U
    map!(x->abs2(x)/length(pts),v,v)
    if adjust_dim
        thres = sum(v)*convert(T,adjust_dim_ratio)
        k = findfirst(v->v > thres, cumsum(v))
        M.d = k === nothing ? d : k
    end
    M.μ = vec(μ)
    M.basis =  U
    return M
end

function origstats(H::Vector{T}) where T<:Real
    N = length(H)

    # calculate threshold
    S = zeros(N,6)

    # recursive defintions
    S[1,1] = H[1]
    S[N-1,2] = H[N]
    S[1,3] = 1.0
    S[N-1,4] = (H[N] == 0 ? 0 : N-1)
    S[1,5] = 0
    S[N-1,6] = 0
    i = 2
    j = N-2
    while i <= N-1
        S[i,1] = S[i-1,1] + H[i]
        if S[i,1] != 0
            S[i,3] = ((S[i-1,3] * S[i-1,1]) + ((i-1) * H[i])) / S[i,1]
            S[i,5] = (S[i-1,1] *
                       (S[i-1,5] + (S[i-1,3]-S[i,3]) * (S[i-1,3]-S[i,3])) +
                        H[i] * ((i-1) - S[i,3]) * ((i-1) - S[i,3]) ) / S[i,1]
        end

        S[j,2] = S[j+1,2] + H[j+1]
        if S[j+1,2] != 0
            S[j,4] = ((S[j+1,4] * S[j+1,2]) + (j * H[j+1])) / S[j,2]
            S[j,6] = (S[j+1,2] *
                       (S[j+1,6] + (S[j+1,4]-S[j,4]) * (S[j+1,4]-S[j,4])) +
                        H[j+1] * (j - S[j,4]) * (j - S[j,4]) ) / S[j,2]
        end

        i += 1
        j -= 1
    end

    return S
end

function refstats(H::Vector{T}) where T<:Real
    x = 1:length(H)
	A=cumsum(H)
	B=cumsum(H.*x)
	C=cumsum(H.*x.^2)

	p=A./A[end]
	q=(A[end] .- A)./A[end]

	u=B./A
    v=(B[end] .- B)./(A[end] .- A)
    v[isnan.(v)] .= 0.0

	s2=C./A - u.^2
    t2=(C[end] .- C)./(A[end].-A) - v.^2
    t2[isnan.(t2)] .= 0.0

	return hcat(p, q, u, v, s2, t2)
end

function histstats(H::Vector{T}) where T<:Real
    N = length(H)
    S = zeros(N,6)

    S[1, 1] = H[1]
    S[1, 3] = 1.0
    S[1, 2] = 1.0 - S[1, 1]
    S[N-1, 2] = H[N]
    S[N-1, 4] = N
    j = N
    for i in 2:N
        P = H[i]
        S[i, 1] = S[i-1, 1] + P
        A = S[i-1, 1]
        if (A + H[i]) != 0
            S[i, 3] = (S[i-1, 3]*A + H[i]*i)/(A + H[i])
            S[i, 5] = ((S[i-1, 5] + S[i-1, 3]*S[i-1, 3])*A + H[i]*i*i)/(A + H[i]) - S[i, 3]*S[i, 3]
        end

        S[j-1, 2] = S[j, 2] + H[j]
        ΔA = S[j-1, 2]
        if ΔA != 0
            S[j-1, 4] = (S[j, 4]*S[j, 2] + H[j]*j)/ΔA
            S[j-1, 6] = ((S[j, 6] + S[j, 4]*S[j, 4])*S[j, 2] + H[j]*j*j)/ΔA - S[j-1, 4]*S[j-1, 4]
        end
        j -= 1
    end

	return S
end

function find_global_min(J::Vector{T}, tol::T) where T<:Real
    N = length(J)

    # Mark minima
    M = zeros(Bool,N)
    if N > 1
        prev = J[2] - J[1]
        curr = 0.0
        for i=2:(N-1)
            curr = J[i+1] - J[i]
            M[i] = prev<=0 && curr>=0
            prev=curr
        end
    end
    @logmsg DEBUG_SEPARATION "Local minima" M=M'

    # Special case: flat minimum so separate at the first bin
    all(M[2:end-1]) && length(M)>2 && return (Inf, 1)

    # Find global minima of criterion funtion if exists
    lmin = findfirst(isequal(true),M)
    lmin === nothing && throw(LMCLUSException("No minimum found. Presume unimode histogram."))

    depth = 0.0
    global_min = 0
    if lmin == N
        throw(LMCLUSException("No minimum found, unimode histogram"))
    else
        while lmin < N
            # Detect flat
            rmin = lmin
            while rmin<N && M[rmin]
                rmin += 1
            end
            loc_min=( lmin + rmin - 1 ) >> 1

            # Monotonically ascend to the left
            lheight = loc_min
            while lheight > 1 && J[lheight-1] >= J[lheight]
                lheight -= 1
            end

            # Monotonically ascend to the right
            rheight = loc_min
            while rheight < N && J[rheight] <= J[rheight+1]
                rheight += 1
            end

            # Compute depth
            local_depth = 0
            local_depth = (J[lheight] < J[rheight] ? J[lheight] : J[rheight]) - J[loc_min]

            if local_depth > depth
                depth = local_depth
                global_min = loc_min
            end

            lmin = rmin
            while lmin<N && !M[lmin]
                lmin += 1
            end
        end
    end

    if depth < tol
        throw(LMCLUSException("Separation depth ($depth) beyond the tolerance level ($tol). Presume unimode histogram."))
    end

    depth, global_min
end

"Projection of the data to the manifold or its orthoganal comppliment"
function project(MC::Manifold{T}, X::AbstractMatrix{T}; ocs = true) where T<:Real
    M = outdim(MC)
    N = min(size(X)...)

    @assert hasfullbasis(MC) "$MC must have full basis"
    @assert M < N "Manifold dimension must be less then of the full space"

    # change basis for data
    r = ocs ? ((M+1):N) : (1:M)
    BC = MC.basis[:,r]
    Y = transpose(BC) * (X.-mean(MC))

    return Y
end

ext: jl


text:
using Plots

""" Helper Function for 2D contour plot for linear functions f(x) = c'x"""
function plot_lin_contour!(c, grid_res; kwargs...)
    f(x) = dot(c,x)
    contour!(grid_res, grid_res, (x,y) -> f([x,y]); kwargs...)
end

""" Helper Function for Plotting an arrow in 2D"""
function plot_vec!(x, label="", x_start= [0.0,0.0]; kwargs...)
    x_end = x
    x_pos = [x_start[1], x_end[1]]
    y_pos = [x_start[2], x_end[2]]
    plot!(x_pos, y_pos, arrow=(1,:black); kwargs...)
end


text:
"""
    noise_psd(x::Vector{<:Real}, smoothing, samplerate, nfft, hop=nfft÷2)
    noise_psd(X::Matrix{<:Complex}, smoothing, framerate)

Estimate the noise power-spectral-density of a noisy time-domain signal `x`,
or an STFT `X` using minimum statistics. This is similar to the approach in [1], but assumes a constant noise PSD. To reduce variance the signal will be
smoothed with a gaussian window `smoothing` seconds long. Genrally longer windows will provide a less-biased estimate, but it must be short enough that there will be some regions of the signal where the the smoothing window has none of the target signal in it.

`samplerate` is used for time-domain signals, vs `framerate` for STFTs. You can
convert between them with `framerate = samplerate / hop`. You can set the
samplerate/framerate to 1 to give smoothing in samples of frames, respectively.

The result will be of size `nfft÷2+1` for a time signal `x`, or the same height
as an STFT `X`.

[1] Martin, R., Spectral subtraction based on minimum statistics, 1994
"""
function noise_psd(x::AbstractVector{<:Real},
                   smoothing, samplerate,
                   nfft, hop=nfft÷2)
    noise_psd(stft2(x, nfft, hop), smoothing, samplerate/hop)
end

function noise_psd(X::AbstractMatrix{<:Complex}, smoothing, framerate)
    L = round(Int, smoothing * framerate)
    if iseven(L)
        L += 1
    end
    win = gaussian(L, 0.14)
    win ./= sum(win)
    map(1:size(X, 1)) do band
        # we skip the last STFT frame because it is often lower-power due to
        # being only partially-filled. Note that it would be more correct to
        # actually take the STFT hop and window sizes into account, but in
        # in practice we're often using 50% overlap and this effect shouldn't be
        # that important.
        sm = conv(abs2.(X[band, 1:end-1]), win)
        minimum(sm[L:end-L+1])
    end
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
import Pkg
Pkg.add(path="https://github.com/RoyCCWang/RKHSRegularization")
Pkg.add("Genie")

using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Json

import RKHSRegularization

#Genie.config.run_as_server = true

# test case to see if server is up.
route("/hello") do
  "Hello World"
end

function errorcheck1Dregression(x::Vector{T}, y::Vector{T}, xq::Vector{T}, θ_len::T, σ²::T)::Bool where T
    # check length.
    if length(x) != length(y)
        return false
    end

    # check finite input.
    if !(all(isfinite.(x)) & all(isfinite.(y)) & all(isfinite.(xq)))
        return false
    end

    # check if positive numbers for the hyperparameters.
    if !(0 < θ_len < Inf && 0 < σ² < Inf)
        return false
    end

    return true
end

# unrecognized inputs.
function errorcheck1Dregression(x, y, xq, θ_len, σ²)::Bool
    return false
end

function geterrormessage1Dregression()
    return "Please check your inputs. The input format is:
    input_array = Vector{Vector{Float64}}(undef, 4)
    input_array[1] = x
    input_array[2] = y
    input_array[3] = xq
    input_array[4] = [θ_len; σ²]"
end

function run1Dregression!(yq::Vector{T}, x::Vector{T}, y::Vector{T}, xq::Vector{T}, θ_len::T, σ²::T)::Nothing where T

    θ = RKHSRegularization.Spline34KernelType(θ_len)

    N = length(x)
    Nq = length(xq)

    # do GP regression.
    X = collect( [x[n]] for n = 1:N )

    η = RKHSRegularization.RKHSProblemType( zeros(Float64,length(X)),
                     X,
                     θ,
                     σ²)
    RKHSRegularization.fitRKHS!(η, y)

    # query.
    resize!(yq, Nq)
    Xq = collect( [xq[n]] for n = 1:Nq )
    RKHSRegularization.query!(yq, Xq, η)

    return nothing
end

route("/jsonpayload", method = POST) do
  #@show jsonpayload()
  #@show rawpayload()

  incoming_array = jsonpayload()["input_array"]
  #println("typeof(incoming_array) = ", incoming_array)
  x = convert(Vector{Float64}, incoming_array[1])
  y = convert(Vector{Float64}, incoming_array[2])
  xq = convert(Vector{Float64}, incoming_array[3])

  hyperparams = convert(Vector{Float64}, incoming_array[4])
  θ_len = hyperparams[1]
  σ² = hyperparams[2]

  # error check.
  flag = errorcheck1Dregression(x, y, xq, θ_len, σ²)

  yq = Vector{Float64}(undef, 0)
  if flag

    run1Dregression!(yq, x, y, xq, θ_len, σ²)

    # assemble output.
    println("Done yq")
  else
    println("Bad input")
    out = geterrormessage1Dregression()
    #push!(yq, Inf)
  end

  json(yq)
end

up(async = false) # block the execution of this script.
#up()
#Genie.startup()

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
#=

File defines the abstract types used within the package.

=#

abstract type AbstractProblem{T} end

"""
    abstract type AbstractPollDirectionGenerator end

Parent type for any type used for implementing a direction generator.

Generally any direction generators should return vectors with unit length of
one

"""
abstract type AbstractPoll end

abstract type AbstractSearch end

"""
    abstract type AbstractMesh end

Parent type of any struct implementing the construction of a mesh. To maintain
compatibility with other aspects of the package, the naming convention for
variables within structs must be followed. These respect the notation used
within Audet & Dennis 2006.
"""
abstract type AbstractMesh end

abstract type AbstractConstraint end

abstract type AbstractCache end


text:
using FermiDiracIntegrals
using Test

@testset "FermiDiracIntegrals.jl" begin
    for x in -10:10
        @test isapprox(F(0,x),log(1+exp(x)))
    end

    for x in -100:0.1:100
        @test isapprox(F(0.5,x),F(Val(1/2),x),rtol=0.03)
    end

    ## https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1950.0183
    @test isapprox(F(1,-4),0.0182324,atol=1e-7)
    @test isapprox(F(2,-4),0.0182739,atol=1e-7)
    @test isapprox(F(3,-4),0.0182947,atol=1e-7)
    @test isapprox(F(4,-4),0.0183052,atol=1e-7)

    @test isapprox(F(1,-1.3),0.2559184,atol=1e-7)
    @test isapprox(F(2,-1.3),0.2639215,atol=1e-7)
    @test isapprox(F(3,-1.3),0.2681202,atol=1e-7)
    @test isapprox(F(4,-1.3),0.2702891,atol=1e-7)

    @test isapprox(F(1,0),0.8224670,atol=1e-7)
    @test isapprox(F(2,0),0.9015427,atol=1e-7)
    @test isapprox(F(3,0),0.9470328,atol=1e-7)
    @test isapprox(F(4,0),0.9721198,atol=1e-7)
end


text:
@testset "893.groups-of-special-equivalent-strings.jl" begin
    @test num_special_equiv_groups(["abc","acb","bac","bca","cab","cba"]) == 3
    @test num_special_equiv_groups(["abcd","cdab","cbad","xyzz","zzxy","zzyx"]) == 3
end


text:
# Copyright 2017 Raytheon BBN Technologies
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the license at
#
#		http://www.apache.org/licenses/LICENSE-2.0

using LsqFit

## TODO: Devise a way to get specify a fit as a string...
## TODO: Fix types...

"""
FitResult type that stores the output of a fit.
"""
struct FitResult
    fit_params::Dict{String,Float64}
    sq_error::Float64
    Nσ::Float64
    errors::Dict{String,Float64}
    fit_curve::Function
    model_str::String
end

"""
`generic_fit(xpts, ypts, model, initial_guess, fit_params, model_string; yvars=[])`
Helper function to return a fit to a model.
# Arguments
* xpts: Independent fit variable.
* ypts: Dependent fit variable.
* model: Fit function(x, p) of data x and parameter vector p.
* initial_guess: Initial guess of fit parameters.
* fit_params::function(Array{Float64} => Dict{String, Float64}): Result dict of fit parameters.
* model_string::String: String indicating fit function.
* yvars: Variance of y points. Defaults to empty vector.
# Returns
* A FitResult object for the specified fit.
"""
function generic_fit(xpts, ypts, model, initial_guess, fit_params, model_string::String; yvars=[])
    @assert length(xpts) == length(ypts) "X and Y data length must match."
    #try
        if isempty(yvars)
          result = curve_fit(model, xpts, ypts, initial_guess)
          errors = margin_error(result)
          sq_error = sum(((model(xpts, result.param) - ypts) .^ 2))
        else
          @assert length(ypts) == length(yvars) "Y data and Y variance lengths must match."
          result = curve_fit(model, xpts, ypts, 1 ./ sqrt(yvars), initial_guess)
          errors = margin_error(result)
          sq_error = sum(((model(xpts, result.param) - ypts) .^ 2) ./ yvars)
        end
    #catch
    #    result =
    # Compute badness of fit:
    # Under the null hypothesis (that the model is valid and that the observations
    # do indeed have Gaussian statistics), the mean squared error is χ² distributed
    # with `dof` degrees of freedom. We can quantify badness-of-fit in terms of how
    # far the observed MSE is from the expected value, in units of σ = 2dof (the expected
    # standard deviation for the χ² distribution)
    dof = length(xpts)-length(initial_guess)
    Nσ = sq_error/sqrt(2dof) - dof/sqrt(2dof)
    return FitResult( fit_params(result.param),
              sq_error,
              Nσ,
              fit_params(errors),
              xpts->model(xpts,result.param),
              model_string)
end

""" `fit_t1(xpts, ypts, yvars=[])`
Fit to exponential T decay y = a exp(-t/T) + b.
"""
function fit_t1(xpts, ypts, yvars=[])
    T1_fit_dict(p) = Dict("a" => p[1], "T" => p[2], "b" => p[3])
    model(t, p) = p[1]*exp.(-t ./ p[2]) .+ p[3]
    p_guess = [ypts[1]-ypts[end], xpts[end]/3., ypts[end]]
    return generic_fit(xpts, ypts, model, p_guess, T1_fit_dict,
              "a * exp(-t/T) + b", yvars=yvars)
end

""" `fit_double_exp(xpts, ypts, yvars=[])`
Fit to double exponential decay y = a exp(-t/Ta) + b exp(-t/Tb) + c.
"""
function fit_double_exp(xpts, ypts, yvars=[])
    T1_fit_dict(p) = Dict("a" => p[1], "Ta" => p[2], "b" => p[3], "Tb" => p[4], "c" => p[5])
    model(t, p) = p[1]*exp.(-t ./ p[2]) .+ p[3]*exp.(-t ./ p[4]) .+ p[5]
    p_guess = [0.5*(ypts[1]-ypts[end]), xpts[end]/3., 0.5*(ypts[1]-ypts[end]), xpts[end]/3., ypts[end]]
    return generic_fit(xpts, ypts, model, p_guess, T1_fit_dict,
              "a * exp(-t/Ta) + b * exp(-t/Tb) + c", yvars=yvars)
end

""" `fit_line(xpts, ypts, yvars=[])`
Fit to linear model y = ax + b.
"""
function fit_line(xpts, ypts, yvars=[])
    model(t, p) = p[1]*t .+ p[2]
    p_guess = [-ypts[findmin(abs.(xpts))[2]] / xpts[findmin(abs.(ypts))[2]], ypts[findmin(abs.(xpts))[2]]]
    return generic_fit(xpts, ypts, model, p_guess,
    x->Dict("a" => x[1], "b"=>x[2]), "ax + b", yvars=yvars)
end

""" `fit_ramsey(xpts, ypts, yvars=[])`
Fit data to a Ramsey decay of the form a*exp(-t ./ T).*cos(2πf .* t + ϕ) + b
with optional weights for each point.
"""
function fit_ramsey(xpts, ypts, yvars=[])

    fit_dict1(p) = Dict("a"=>p[1], "T"=>p[2], "f"=>p[3], "ϕ"=>p[4], "b"=>p[5])
    model(t, p) = p[1]*exp.(-t ./ p[2]).*cos.(2π*p[3] .* t .+ p[4]) .+ p[5]
    # Use KT estimation to get a guess for the fit
    freqs,Ts,amps = KT_estimation(ypts .- mean(ypts), xpts[2]-xpts[1], 1)
    p_guess = [abs(amps[1]), Ts[1], freqs[1], angle(amps[1]), mean(ypts)];

    return generic_fit(xpts, ypts, model, p_guess,
              fit_dict1, "a*exp(-t/T).*cos(2πf .* t + ϕ) + b)", yvars=yvars)
end

""" `fit_twofreq_ramsey(xpts, ypts, yvars=[])`

Fit the Ramsey data to a one-frequency model a*exp(-t ./ T).*cos(2πf .* t + ϕ) + b
or a two-frequency model
a₁*exp(-t ./ T₁).*cos(2πf₁ .* t + ϕ₁) + a₂*exp(-t ./ T₂).*cos(2πf₂ .* t + ϕ₂) + b
optionally taking variances into account, and report which should be
used based on Akaike's Information Criterion (AIC). This assumes
Gaussian statistics for each observation.
"""
function fit_twofreq_ramsey(xpts, ypts, yvars=[])
    fit_dict2(p) = Dict("a₁"=>p[1], "T₁"=>p[2], "f₁"=>p[3], "ϕ₁"=>p[4],
    "a₂"=>p[5], "T₂"=>p[6], "f₂"=>p[7], "ϕ₂"=>p[8],
    "b"=>p[9])
    model2(t, p) = ( p[1] * exp.(-t ./ p[2]) .* cos.(2π * p[3] .*t .+ p[4]) .+
    p[5] * exp.(-t ./ p[6]) .* cos.(2π * p[7] .*t .+ p[8]) .+ p[9] )

    #Use KT estimation to get a guess for the fit
    freqs,Ts,amps = KT_estimation(ypts.-mean(ypts), xpts[2].-xpts[1], 2)
    phases = angle.(amps)
    amps = abs.(amps)
    p_guess = [1, Ts[1], freqs[1], phases[1], 1, Ts[2], freqs[2], phases[2], mean(ypts)]

    fitresult2 = generic_fit(xpts, ypts, model2, p_guess, fit_dict2,
    "a₁*exp(-t ./ T₁).*cos(2πf₁ .* t .+ ϕ₁) .+ a₂*exp(-t ./ T₂).*cos(2πf₂ .* t .+ ϕ₂) .+ b",
    yvars=yvars)

    fitresult1 = fit_ramsey(xpts, ypts, yvars)

    # AIC computation
    k1 = 5
    k2 = 9
    corr(k,n) = (k+1)*(k+1)/(n-k-2)
    aicc(e,k,n) = 2 * k + e + corr(k,n)

    aic = aicc(fitresult2.sq_error,k2,length(xpts)) - aicc(fitresult1.sq_error,k1,length(xpts))

    println(aic)
    return fitresult2 #(aic > 0) ? fitresult2 : fitresult1
end

"""`fit_sin(xpts, ypts, yvars)`
Fit to a sine a*sin(2πf t) + b
"""
function fit_sin(xpts, ypts, yvars=[])
    model(t, p) = p[1]*sin.(2π*p[2] .* t) .+ p[3]
    # Use KT estimation to get a guess for the fit
    freqs,Ts,amps = KT_estimation(ypts, xpts[2]-xpts[1], 2)
    idx = argmax(abs.(amps))
    p_guess = [abs(amps[idx]), Ts[idx], freqs[idx], mean(ypts)];
    return generic_fit(xpts, ypts, model, p_guess,
              x->Dict("a"=>x[1], "f"=>x[2], "b"=>x[3]),
              "a*sin(2πf t) + b", yvars=yvars)
end

## TODO: Change this to use generic_fit
"""`fit_photon_ramsey(xpts, ypts, params)`
Fit function in McClure et al., Phys. Rev. App. 2016. Params are:
1. cavity decay rate κ (MHz)
2. detuning Δ (MHz)
3. dispersive shift 2χ (MHz)
4. Ramsey decay time T₂* (μs)
5. exp(-t_meas/T₁) (μs)
6. initial qubit state (0/1)
"""
function fit_photon_ramsey(xpts, ypts, params)
    params[1:3]*=2*pi #convert to angular frequencies
    model_0(t, p) = (-imag.(exp.(-(1/params[4]+params[2]*1im).*t + (p[1]-p[2]*params[3]*(1-exp.(-((params[1] + params[3]*1im).*t)))/(params[1]+params[3]*1im))*1im)))
    function model(t, p)
        if params[6] == 1
            return params[5]*model_0(t, p) + (1-params[5])*model_0(t, [p[1]+pi; p[2:end]])
        else
            return model_0(t, p)
        end
    end
    p_guess = [0., 1.]
    result = curve_fit(model, xpts, ypts, p_guess)
    errors = margin_error(result)
    xfine = range(xpts[1], stop=xpts[end],length=1001)
    fit_curve = (xfine, model(xfine, result.param))
    return (result.param[2], errors[2], fit_curve)
end

"""
analyzeRB(ypts, seqlengths; purity=False)
analyzeRB Analyzes a randomized benchmarking experiment
seqlengths, example: [4,8,16,32,64,128,256,384,512,640,768,896,1024,1152,1280,1408,1536,1664]
if purity = true, fit for incoherent errors. See J.J. Wallman et al., New J. Phys. 17, 113020 (2015)
Assuming data of the form <Z>s1	, <Z>s2, ...., <Z>sN, <X>s1, ..., <X>sN, ..., <Y>s1, ..., <Y>sN
with seqlengths = [s1, s2, ..., sN]
"""
function analyzeRB(ypts, seqlengths; purity=false)
    # figure out how many sequences of each length we have
    num_repeats = length(ypts) ÷ length(seqlengths) ÷ (purity ? 3 : 1)
    xpts = repeat(seqlengths, inner=num_repeats)
    if purity
        # compute length of Bloch vector
        data = vec(sum(reshape(ypts, (length(ypts) ÷ 3),3).^2,  2))
    else
        # otherwise convert <Z> to prob of 0
        data = .5 * (1 + vec(ypts))
    end
    model(n, p) = p[1] * (1-p[2]).^n + p[3]
    fit = curve_fit(model, xpts-purity, data, [0.5, .01, 0.5]) #fit to ...^(n-1) for purity
    xfine = range(seqlengths[1],stop=seqlengths[end],length=1001)
    fit_curve = (xfine, model(xfine, fit.param))
    errors = margin_error(fit)
    if purity
        println("ϵ_inc = $(0.5*(1-sqrt(1-fit.param[2]))*100)")
    else
        println("ϵ = $(fit.param[2]/2*100)")
    end
    return (xpts, data, fit.param, fit_curve, errors)
end

"""`fit_RB(xpts, ypts, yvars=[])`
Fit to RB decay a*(1-b)^x + c
"""
function fit_RB(xpts, ypts, yvars=[])
    RB_fit_dict(p) = Dict("a" => p[1], "b" => p[2], "c" => p[3])
    model(n, p) = p[1] * (1-p[2]).^n + p[3]
    p_guess = [0.5, .01, 0.5]
    return generic_fit(xpts, ypts, model, p_guess, RB_fit_dict, "a * exp(-n/b) + c", yvars=yvars)
end

function calculate_residueCR(params, t, rvec)
    rx        = params[1]; # Rabi frequency along x
    ry        = params[2]; # Rabi frequency along y
    delta     = params[3]; # detuning
    x0        = params[4]; # initial x
    y0        = params[5]; # initial y
    z0        = params[6]; # initial z
    omega = sqrt(rx^2+ry^2+delta^2);
    x = (rx*delta*(-z0+x0*rx+y0*ry)+(rx*(z0*delta-y0*ry)+x0*(delta^2+ry^2)).*cos.(t*omega)+omega*(y0*delta+z0*ry).*sin.(t*omega))/omega^2;
    y = (ry*delta*(-z0+x0*rx+y0*ry)+(ry*(z0*delta-x0*rx)+y0*(delta^2+rx^2)).*cos.(t*omega)-omega*(x0*delta+z0*rx).*sin.(t*omega))/omega^2;
    z = (delta*(z0*delta-x0*rx-y0*ry)+(x0*delta*rx+y0*delta*ry+z0*(rx^2+ry^2)).*cos.(t*omega)+omega*(y0*rx-x0*ry).*sin.(omega*t))/omega^2;
    rmodel = [x,y,z];
    diffvec = rvec-rmodel;
    return sum(hcat(diffvec...).^2), rmodel
end

function fit_CR_hamiltonian(xpoints, datavec, params_estimate = [5,5,0,0,0,1])
    result = optimize(x-> calculate_residueCR(x, xpoints, datavec)[1], params_estimate)
    return Optim.minimizer(result)
end


text:
# using Base.Test
using IndexedArrays

const pdb_list = ["2HHB", "1IRD", "1DN3"]

print("""

Tests for IndexedArray
=======================
""")

print("""

Creation
--------
""")
const pdbs = IndexedArray(pdb_list)
@test pdbs.items == pdb_list
@test_throws IndexedArrayError IndexedArray([1,2,3,3])

print("""

Mapping
-------
""")

print("""
Test selectvalue
""")
for i in 1:3
  @test pdbs[i] == pdb_list[i]
end
@test_throws BoundsError pdbs[4]

@test pdbs[ [2,3] ]  == ["1IRD", "1DN3"]
@test pdbs[ Bool[false, true, true] ] == ["1IRD", "1DN3"]
@test pdbs[ [1,2,3] .> 1 ] == ["1IRD", "1DN3"]

@test pdbs[:] == pdb_list

print("""
Test selectindex
""")
i = 1
for pdb in pdb_list
  @test findfirst(pdbs, pdb) == i
  i += 1
end

@test_throws KeyError findfirst(pdbs, "2TRX")

@test [ findfirst(pdbs, i) for i in ["2HHB", "1IRD"] ] == Int[1,2]

@test [ findfirst(pdbs, i) for i in pdbs ] == Int[1,2,3]

print("""

Comparisons and getindex
------------------------
""")

@test IndexedArray(pdbs[[2,3]]) == IndexedArray(["1IRD", "1DN3"])

@test IndexedArray(pdbs[[false, true, true]]) == IndexedArray(["1IRD", "1DN3"])
@test IndexedArray(pdbs[[true, true, false]]) != IndexedArray(["1IRD", "1DN3"])

print("""

Iterations
----------
""")

@test collect(pdbs) == ["2HHB", "1IRD", "1DN3"]

i = 1
for element in pdbs
  @test element == pdb_list[i]
  i += 1
end

print("""

swap! and copy
--------------
""")

const copyofpdbs = copy(pdbs)

swap!(copyofpdbs, 2, 3)
@test findfirst(copyofpdbs, "1IRD") == 3
@test findfirst(copyofpdbs, "1DN3") == 2

@test copyofpdbs[2] == pdbs[3]
@test copyofpdbs[3] == pdbs[2]
@test findfirst(copyofpdbs, "1IRD") == findfirst(pdbs, "1DN3")
@test findfirst(copyofpdbs, "1DN3") == findfirst(pdbs, "1IRD")

print("""

push!
-----
""")
push!(pdbs, "2TRX")
@test_throws IndexedArrayError push!(pdbs, "2TRX")
@test pdbs[4] == "2TRX"
@test findfirst(pdbs, "2TRX") == 4


text:
function OrekitFixedStepHandlerMultiplexer()
    return OrekitFixedStepHandlerMultiplexer(())
end

function add(obj::OrekitFixedStepHandlerMultiplexer, arg0::OrekitFixedStepHandler)
    return jcall(obj, "add", void, (OrekitFixedStepHandler,), arg0)
end

function handle_step(obj::OrekitFixedStepHandlerMultiplexer, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "handleStep", void, (SpacecraftState, jboolean), arg0, arg1)
end

function init(obj::OrekitFixedStepHandlerMultiplexer, arg0::SpacecraftState, arg1::AbsoluteDate, arg2::jdouble)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate, jdouble), arg0, arg1, arg2)
end


text:
@testset "Check Constructor" begin

    @testset "No SE" begin
        for x in [1, .1, "1", π]
            @test FormattedNumber(x).val == x
            @test FormattedNumber(x).star == 0
        end
        @test FormattedNumber(1).format     == "{:,n}"
        @test FormattedNumber(.1).format    == "{:.3f}"
        @test FormattedNumber(1e-3).format  == "{:.3f}"
        @test FormattedNumber(1e-4).format  == "{:.3e}"
        @test FormattedNumber(1e3).format   == "{:.3f}"
        @test FormattedNumber(1e6).format   == "{:.3e}"

        @test_throws ErrorException FormattedNumber(Complex(1, 1))
    end

    @testset "SE" begin
        for x in [1, "1"]
            @test_throws AssertionError FormattedNumber(x, .2)
        end

        @test FormattedNumber(1.0, 0.1).format    == "{:.3f}"
        @test FormattedNumber(1.0, 0.1).format_se == "{:.3f}"

        @test FormattedNumber(1.0e4, 0.1).format    == "{:.3f}"
        @test FormattedNumber(1.0e4, 0.1).format_se == "{:.3f}"

        @test FormattedNumber(1.0e5, 0.1).format    == "{:.3e}"
        @test FormattedNumber(1.0e5, 0.1).format_se == "{:.3f}"

        @test FormattedNumber(1.0e-4, 0.1).format    == "{:.3e}"
        @test FormattedNumber(1.0e-4, 0.1).format_se == "{:.3f}"

        @test FormattedNumber(1.0, 1e4).format_se  == "{:.3f}"
        @test FormattedNumber(1.0, 1e5).format_se  == "{:.3e}"

        @test FormattedNumber(1.0, 1e-3).format_se  == "{:.3f}"
        @test FormattedNumber(1.0, 1e-4).format_se  == "{:.3e}"

        x = FormattedNumber(1.0, .1)
        @test x.val  == 1.0
        @test x.se   == 0.1
        @test x.star == 0
        star!(x, 1)
        @test x.star == 1

        # Check automatic conversion to FNumSE
        @testset "Promotion"  begin
            x = FormattedNumber(1.0)
            y = FormattedNumber(1.0, 0.1)

            x1, y1 = promote(x, y)
            @test y1 == y
            @test x1.val == x.val
            @test isnan(x1.se)
        end
    end
end

@testset "Showing Formatted Numbers" begin

    @test sprint(show, FormattedNumber(.1)) == "0.100"
    @test sprint(show, FormattedNumber(1)) == "1"
    @test sprint(show, FormattedNumber("test")) == "test"
    @test sprint(show, FormattedNumber(true)) == "true"

    @test sprint(show, FormattedNumber(1.0, 0.1))   == "1.000 (0.100)"
    @test sprint(show, FormattedNumber(1.0, 1e-3))  == "1.000 (0.001)"
    @test sprint(show, FormattedNumber(1.0, 1e-4))  == "1.000 (1.000e-04)"
    @test sprint(show, FormattedNumber(1.0, 1e4))   == "1.000 (10000.000)"
    @test sprint(show, FormattedNumber(1.0, 1e5))   == "1.000 (1.000e+05)"

    @test sprint(show, FormattedNumber(1.0, NaN))   == "1.000"
end


text:
function get_ecef(obj::GLONASSAnalyticalPropagator)
    return jcall(obj, "getECEF", Frame, ())
end

function get_eci(obj::GLONASSAnalyticalPropagator)
    return jcall(obj, "getECI", Frame, ())
end

function get_frame(obj::GLONASSAnalyticalPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_glonass_orbital_elements(obj::GLONASSAnalyticalPropagator)
    return jcall(obj, "getGLONASSOrbitalElements", GLONASSOrbitalElements, ())
end

function get_mu(::Type{GLONASSAnalyticalPropagator})
    return jcall(GLONASSAnalyticalPropagator, "getMU", jdouble, ())
end

function propagate_in_ecef(obj::GLONASSAnalyticalPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagateInEcef", PVCoordinates, (AbsoluteDate,), arg0)
end

function reset_initial_state(obj::GLONASSAnalyticalPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end


text:
type RacosOptimization
    best_solution
    algorithm

    function RacosOptimization()
        return new(Nullable(), Nullable())
    end
end

# General optimization function, it will choose optimization algorithm according to parameter.get_sequential()
# If user hasn't define uncertain_bits in parameter, set_ub() will set uncertain_bits automatically according to dim
# in objective
function opt!(ro::RacosOptimization, objective, parameter)
    ro_clear!(ro)
    uncertain_bits = set_ub(objective)
    if parameter.sequential == true
        ro.algorithm = SRacos()
        ro.best_solution = sracos_opt!(ro.algorithm, objective, parameter, ub=uncertain_bits)
    else
        ro.algorithm = Racos()
        ro.best_solution = racos_opt!(ro.algorithm, objective, parameter, ub=uncertain_bits)
    end
    return ro.best_solution
end

function ro_clear!(ro::RacosOptimization)
    ro.best_solution = Nullable()
    ro.algorithm = Nullable()
end

# Set uncertain_bits
function set_ub(objective)
    dim = objective.dim
    dim_size = dim.size
    discrete = is_discrete(dim)
    if discrete==false
        if dim_size <= 100
            ub = 1
        elseif dim_size <= 1000
            ub = 2
        else
            ub = 3
        end
    else
        if dim_size <= 10
            ub = 1
        elseif dim_size <= 50
            ub = 2
        elseif dim_size <= 100
            ub = 3
        elseif dim_size <= 1000
            ub = 4
        else
            ub = 5
        end
    end
    return ub
end


text:
#!/bin/bash
#=
exec julia --color=yes make.jl
=#

push!(LOAD_PATH, "../src/")
using Documenter, StochasticArithmetic

makedocs(
    modules = [StochasticArithmetic],
    format = :html,
    checkdocs = :exports,
    sitename = "StochasticArithmetic.jl",
    # doctest = :fix,
    pages = Any["Home" => "index.md",
                "Error-Free Transformations" => "EFT.md",
                "Rounding modes" => "rounding.md",
                "Stochastic FP types" => "sfloat.md"]
)

deploydocs(
    repo = "github.com/ffevotte/StochasticArithmetic.jl.git",
)

# Local Variables:
# mode: julia
# End:


text:
module SmolyakInterpolation

export
    # chebyshev
    Chebyshev, set_length, set_range, gridpoint_set, evaluate, evaluate!,
    # basis
    HomogeneousBasis, degrees_of_freedom, basis_matrix_and_coordinates, interpolate,
    interpolated_basis!, interpolated_basis

using ArgCheck: @argcheck
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using UnPack: @unpack
using StaticArrays: SVector

include("utilities.jl")
include("chebyshev.jl")
include("basis.jl")

end # module


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

@test integrate(x^2*(a*x^2+b*x^3+c*x^4), x) == :(1//5*a*x^5+1//6*b*x^6+1//7*c*x^7)
@test integrate(x*(a*x^2+b*x^3+c*x^4), x) == :(1//4*a*x^4+1//5*b*x^5+1//6*c*x^6)
@test integrate(a*x^2+b*x^3+c*x^4, x) == :(1//3*a*x^3+1//4*b*x^4+1//5*c*x^5)
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4), x) == :((1/2)*a*x^2+1//3*b*x^3+1//4*c*x^4)
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4), x) == :(a*x+(1/2)*b*x^2+1//3*c*x^3)
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^2, x) == :(1//7*a^2*x^7+1//9*x^9*(b^2+2*a*c)+1//11*c^2*x^11+1//4*a*b*x^8+1//5*b*c*x^10)
@test integrate(x*(a*x^2+b*x^3+c*x^4)^2, x) == :(1//6*a^2*x^6+1//8*x^8*(b^2+2*a*c)+1//10*c^2*x^10+2//7*a*b*x^7+2//9*b*c*x^9)
@test integrate((a*x^2+b*x^3+c*x^4)^2, x) == :(1//5*a^2*x^5+1//7*x^7*(b^2+2*a*c)+1//9*c^2*x^9+1//3*a*b*x^6+1//4*b*c*x^8)
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^2, x) == :(1//4*a^2*x^4+1//6*x^6*(b^2+2*a*c)+1//8*c^2*x^8+2//5*a*b*x^5+2//7*b*c*x^7)
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^2, x) == :(1//3*a^2*x^3+1//5*x^5*(b^2+2*a*c)+1//7*c^2*x^7+(1/2)*a*b*x^4+1//3*b*c*x^6)
@test integrate(x^5*(a*x^2+b*x^3+c*x^4)^-1, x) == :((1/2)*c^-1*x^2+(1/2)*c^-3*(b^2+-1*a*c)*log(a+b*x+c*x^2)+-1*b*x*c^-2+b*c^-3*(b^2+-4*a*c)^-1//2*(b^2+-3*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^4*(a*x^2+b*x^3+c*x^4)^-1, x) == :(x*c^-1+-1//2*b*c^-2*log(a+b*x+c*x^2)+-1*c^-2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^3*(a*x^2+b*x^3+c*x^4)^-1, x) == :((1/2)*c^-1*log(a+b*x+c*x^2)+b*c^-1*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^-1, x) == :(-2*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^-1, x) == :(a^-1*log(x)+-1//2*a^-1*log(a+b*x+c*x^2)+b*a^-1*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate((a*x^2+b*x^3+c*x^4)^-1, x) == :(-1*a^-1*x^-1+(1/2)*b*a^-2*log(a+b*x+c*x^2)+-1*b*a^-2*log(x)+-1*a^-2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^-1, x) == :(-1//2*a^-1*x^-2+b*a^-2*x^-1+a^-3*(b^2+-1*a*c)*log(x)+-1//2*a^-3*(b^2+-1*a*c)*log(a+b*x+c*x^2)+b*a^-3*(b^2+-4*a*c)^-1//2*(b^2+-3*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^-1, x) == :(-1//3*a^-1*x^-3+(1/2)*b*a^-2*x^-2+-1*a^-3*x^-1*(b^2+-1*a*c)+(1/2)*b*a^-4*(b^2+-2*a*c)*log(a+b*x+c*x^2)+-1*b*a^-4*(b^2+-2*a*c)*log(x)+-1*a^-4*(b^2+-4*a*c)^-1//2*(b^4+2*a^2*c^2+-4*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^8*(a*x^2+b*x^3+c*x^4)^-2, x) == :(-1*b*c^-3*log(a+b*x+c*x^2)+x*c^-2*(b^2+-4*a*c)^-1*(2*b^2+-6*a*c)+x^3*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(2a+b*x)+-1*b*c^-1*x^2*(b^2+-4*a*c)^-1+-1*c^-3*(b^2+-4*a*c)^-3//2*(2*b^4+12*a^2*c^2+-12*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^7*(a*x^2+b*x^3+c*x^4)^-2, x) == :((1/2)*c^-2*log(a+b*x+c*x^2)+x^2*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(2a+b*x)+-1*b*x*c^-1*(b^2+-4*a*c)^-1+b*c^-2*(b^2+-4*a*c)^-3//2*(b^2+-6*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^6*(a*x^2+b*x^3+c*x^4)^-2, x) == :(4*a*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x))+x*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(2a+b*x))
@test integrate(x^5*(a*x^2+b*x^3+c*x^4)^-2, x) == :((b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(2a+b*x)+-2*b*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^4*(a*x^2+b*x^3+c*x^4)^-2, x) == :(-1*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b+2*c*x)+4*c*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^3*(a*x^2+b*x^3+c*x^4)^-2, x) == :(a^-2*log(x)+-1//2*a^-2*log(a+b*x+c*x^2)+a^-1*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b^2+-2*a*c+b*c*x)+b*a^-2*(b^2+-4*a*c)^-3//2*(b^2+-6*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x)))
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^-2, x) == :(b*a^-3*log(a+b*x+c*x^2)+-2*b*a^-3*log(x)+a^-2*x^-1*(b^2+-4*a*c)^-1*(-2*b^2+6*a*c)+-1*a^-3*(b^2+-4*a*c)^-3//2*(2*b^4+12*a^2*c^2+-12*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x))+a^-1*x^-1*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b^2+-2*a*c+b*c*x))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^-2, x) == :(a^-4*(3*b^2+-2*a*c)*log(x)+-1//2*a^-4*(3*b^2+-2*a*c)*log(a+b*x+c*x^2)+-1//2*a^-2*x^-2*(b^2+-4*a*c)^-1*(3*b^2+-8*a*c)+b*a^-4*(b^2+-4*a*c)^-3//2*(3*b^4+30*a^2*c^2+-20*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x))+b*a^-3*x^-1*(b^2+-4*a*c)^-1*(3*b^2+-11*a*c)+a^-1*x^-2*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b^2+-2*a*c+b*c*x))
@test integrate((a*x^2+b*x^3+c*x^4)^-2, x) == :(b*a^-5*(2*b^2+-3*a*c)*log(a+b*x+c*x^2)+-1*a^-5*(b^2+-4*a*c)^-3//2*(4*b^6+-20*a^3*c^3+-30*a*c*b^4+60*a^2*b^2*c^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x))+-1*a^-4*x^-1*(b^2+-4*a*c)^-1*(4*b^4+10*a^2*c^2+-18*a*c*b^2)+-2*b*a^-5*(2*b^2+-3*a*c)*log(x)+-1//3*a^-2*x^-3*(b^2+-4*a*c)^-1*(4*b^2+-10*a*c)+b*a^-3*x^-2*(b^2+-4*a*c)^-1*(2*b^2+-7*a*c)+a^-1*x^-3*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b^2+-2*a*c+b*c*x))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^-2, x) == :(a^-6*(5*b^4+3*a^2*c^2+-12*a*c*b^2)*log(x)+-1//2*a^-6*(5*b^4+3*a^2*c^2+-12*a*c*b^2)*log(a+b*x+c*x^2)+-1//2*a^-4*x^-2*(b^2+-4*a*c)^-1*(5*b^4+12*a^2*c^2+-22*a*c*b^2)+-1//4*a^-2*x^-4*(b^2+-4*a*c)^-1*(5*b^2+-12*a*c)+b*a^-6*(b^2+-4*a*c)^-3//2*(5*b^6+-70*a^3*c^3+-42*a*c*b^4+105*a^2*b^2*c^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x))+b*a^-5*x^-1*(b^2+-4*a*c)^-1*(5*b^4+29*a^2*c^2+-27*a*c*b^2)+a^-1*x^-4*(b^2+-4*a*c)^-1*(a+b*x+c*x^2)^-1*(b^2+-2*a*c+b*c*x)+1//3*b*a^-3*x^-3*(b^2+-4*a*c)^-1*(5*b^2+-17*a*c))
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1//240*x*c^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(7*b^2+-16*a*c)+-1//1920*c^-4*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(105*b^4+256*a^2*c^2+-460*a*c*b^2)+1//40*c^-1*x^2*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+8*c*x)+1//960*b*c^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(35*b^2+-116*a*c)+1//256*b*x*c^-9//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-4*a*c)*(7*b^2+-12*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1//96*c^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(5*b^2+-12*a*c)+1//24*x*c^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+6*c*x)+1//192*b*c^-3*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(15*b^2+-52*a*c)+-1//128*x*c^-7//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-4*a*c)*(5*b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate((a*x^2+b*x^3+c*x^4)^(1/2), x) == :(1//3*c^-1*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(a+b*x+c*x^2)+-1//8*b*c^-2*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+2*c*x)+1//16*b*c^-5//2*x^-1*(a+b*x+c*x^2)^-1//2*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(1//4*c^-1*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+2*c*x)+-1//8*x*c^-3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+-1*x*a^(1/2)*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x))+(1/2)*b*x*c^-1//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-3*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)+x*c^(1/2)*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x))+-1//2*b*x*a^-1//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x)))
@test integrate(x^-4*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1//2*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)+1//8*a^-3//2*(b^2+-4*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+-1//4*b*a^-1*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2))
@test integrate(x^-5*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1//3*x^-4*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//12*b*a^-1*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//16*b*a^-5//2*(b^2+-4*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+1//24*a^-2*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(3*b^2+-8*a*c))
@test integrate(x^-6*(a*x^2+b*x^3+c*x^4)^(1/2), x) == :(-1//4*x^-5*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//24*b*a^-1*x^-4*(a*x^2+b*x^3+c*x^4)^(1/2)+1//96*a^-2*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(5*b^2+-12*a*c)+1//128*a^-7//2*(b^2+-4*a*c)*(5*b^2+-4*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+-1//192*b*a^-3*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(15*b^2+-52*a*c))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//286720*c^-5*(a*x^2+b*x^3+c*x^4)^(1/2)*(1155*b^6+-6720*a^3*c^3+-8988*a*c*b^4+18896*a^2*b^2*c^2)+-1//4480*c^-2*x^3*(a*x^2+b*x^3+c*x^4)^(1/2)*(b*(11*b^2+68*a*c)+10*c*x*(11*b^2+-28*a*c))+1//112*x*c^-1*(a*x^2+b*x^3+c*x^4)^3//2*(3b+14*c*x)+1//35840*c^-3*x^2*(a*x^2+b*x^3+c*x^4)^(1/2)*(99*b^4+560*a^2*c^2+-568*a*c*b^2)+-1//71680*b*x*c^-4*(a*x^2+b*x^3+c*x^4)^(1/2)*(231*b^4+2416*a^2*c^2+-1560*a*c*b^2)+-1//573440*b*c^-6*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(3465*b^6+-58816*a^3*c^3+-30660*a*c*b^4+81648*a^2*b^2*c^2)+3//32768*x*c^-13//2*(b^2+-4*a*c)^2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(33*b^4+16*a^2*c^2+-72*a*c*b^2)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate((a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//7*x*(a*x^2+b*x^3+c*x^4)^3//2+-1//17920*b*c^-4*(a*x^2+b*x^3+c*x^4)^(1/2)*(105*b^4+1168*a^2*c^2+-728*a*c*b^2)+1//280*c^-1*x^3*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+24*a*c+10*b*c*x)+1//35840*c^-5*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(315*b^6+-2048*a^3*c^3+-2520*a*c*b^4+5488*a^2*b^2*c^2)+-1//2240*b*c^-2*x^2*(a*x^2+b*x^3+c*x^4)^(1/2)*(9*b^2+-44*a*c)+1//4480*x*c^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(3*b^2+-4*a*c)*(7*b^2+-32*a*c)+-3//2048*b*x*c^-11//2*(b^2+-4*a*c)^2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(3*b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//3840*c^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(35*b^4+240*a^2*c^2+-216*a*c*b^2)+-1//960*x*c^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(b*(7*b^2+12*a*c)+6*c*x*(7*b^2+-20*a*c))+1//60*c^-1*x^-1*(a*x^2+b*x^3+c*x^4)^3//2*(3b+10*c*x)+-1//7680*b*c^-4*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(105*b^4+1296*a^2*c^2+-760*a*c*b^2)+1//1024*x*c^-9//2*(b^2+-4*a*c)^2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(7*b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//5*c^-1*x^-5*(a*x^2+b*x^3+c*x^4)^5//2+-1//16*b*c^-2*x^-3*(a*x^2+b*x^3+c*x^4)^3//2*(b+2*c*x)+3//128*b*c^-3*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+2*c*x)*(b^2+-4*a*c)+-3//256*b*x*c^-7//2*(b^2+-4*a*c)^2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-3*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//8*c^-1*x^-3*(a*x^2+b*x^3+c*x^4)^3//2*(b+2*c*x)+-1//64*c^-2*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+2*c*x)*(3*b^2+-12*a*c)+3//128*x*c^-5//2*(b^2+-4*a*c)^2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-4*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(1//3*x^-3*(a*x^2+b*x^3+c*x^4)^3//2+1//8*c^-1*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+8*a*c+2*b*c*x)+-1*x*a^3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x))+-1//16*b*x*c^-3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-12*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-5*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(-1*x^-4*(a*x^2+b*x^3+c*x^4)^3//2+1//4*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(9b+6*c*x)+-3//2*b*x*a^(1/2)*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x))+1//8*x*c^-1//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(3*b^2+12*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-6*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(-1//2*x^-5*(a*x^2+b*x^3+c*x^4)^3//2+-1//4*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(3b+-6*c*x)+-1//8*x*a^-1//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(3*b^2+12*a*c)*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x))+3//2*b*x*c^(1/2)*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^-7*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(-1//3*x^-6*(a*x^2+b*x^3+c*x^4)^3//2+-1//4*b*a^-1*x^-5*(a*x^2+b*x^3+c*x^4)^3//2+1//8*a^-1*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+-8*a*c+2*b*c*x)+x*c^3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x))+1//16*b*x*a^-3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-12*a*c)*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x)))
@test integrate(x^-8*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(-1//4*x^-7*(a*x^2+b*x^3+c*x^4)^3//2+-3//128*a^-5//2*(b^2+-4*a*c)^2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+-1//8*x^-4*(a*x^2+b*x^3+c*x^4)^(1/2)*(b+6*c*x)+-1//32*a^-1*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+-12*a*c)+1//64*b*a^-2*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(3*b^2+-20*a*c))
@test integrate(x^-9*(a*x^2+b*x^3+c*x^4)^3//2, x) == :(-1//5*x^-8*(a*x^2+b*x^3+c*x^4)^3//2+-1//40*x^-5*(a*x^2+b*x^3+c*x^4)^(1/2)*(3b+12*c*x)+-1//80*a^-1*x^-4*(a*x^2+b*x^3+c*x^4)^(1/2)*(b^2+-8*a*c)+-1//640*a^-3*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2)*(15*b^4+128*a^2*c^2+-100*a*c*b^2)+3//256*b*a^-7//2*(b^2+-4*a*c)^2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+1//320*b*a^-2*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)*(5*b^2+-28*a*c))
@test integrate(x^3*(a*x^2+b*x^3+c*x^4)^-1//2, x) == :((1/2)*c^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+-3//4*b*c^-2*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+1//8*x*c^-5//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(3*b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^-1//2, x) == :(c^-1*x^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//2*b*x*c^-3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^-1//2, x) == :(x*c^-1//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate((a*x^2+b*x^3+c*x^4)^-1//2, x) == :(-1*a^-1//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x)))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^-1//2, x) == :((1/2)*b*a^-3//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+-1*a^-1*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2))
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^-1//2, x) == :(-1//2*a^-1*x^-3*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//8*a^-5//2*(3*b^2+-4*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+3//4*b*a^-2*x^-2*(a*x^2+b*x^3+c*x^4)^(1/2))
@test integrate(x^7*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :((1/2)*c^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(5*b^2+-12*a*c)+2*x^4*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x)+-2*b*x*c^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+-1//4*b*c^-3*x^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(15*b^2+-52*a*c)+1//8*x*c^-7//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*(15*b^2+-12*a*c)*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^6*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(-2*b*c^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)+2*x^3*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x)+c^-2*x^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(3*b^2+-8*a*c)+-3//2*b*x*c^-5//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x)))
@test integrate(x^5*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(2*x^2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x)+x*c^-3//2*(a+b*x+c*x^2)^(1/2)*(a*x^2+b*x^3+c*x^4)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x+c*x^2)^-1//2*(b+2*c*x))+-2*b*c^-1*x^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2))
@test integrate(x^4*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(2*x*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))
@test integrate(x^3*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(-2*x*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(b+2*c*x))
@test integrate(x^2*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(-1*a^-3//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+2*x*a^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(b^2+-2*a*c+b*c*x))
@test integrate(x*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(3//2*b*a^-5//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+a^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2*b^2+-4*a*c+2*b*c*x)+-1*a^-2*x^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(3*b^2+-8*a*c))
@test integrate((a*x^2+b*x^3+c*x^4)^-3//2, x) == :(-1//8*a^-7//2*(15*b^2+-12*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+a^-1*x^-1*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2*b^2+-4*a*c+2*b*c*x)+-1//2*a^-2*x^-3*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(5*b^2+-12*a*c)+1//4*b*a^-3*x^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(15*b^2+-52*a*c))
@test integrate(x^-1*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(5//16*b*a^-9//2*(7*b^2+-12*a*c)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+a^-1*x^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2*b^2+-4*a*c+2*b*c*x)+-1//3*a^-2*x^-4*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(7*b^2+-16*a*c)+-1//24*a^-4*x^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(105*b^4+256*a^2*c^2+-460*a*c*b^2)+1//12*b*a^-3*x^-3*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(35*b^2+-116*a*c))
@test integrate(x^-2*(a*x^2+b*x^3+c*x^4)^-3//2, x) == :(-1//128*a^-11//2*(315*b^4+240*a^2*c^2+-840*a*c*b^2)*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x))+a^-1*x^-3*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^-1//2*(2*b^2+-4*a*c+2*b*c*x)+-1//4*a^-2*x^-5*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(9*b^2+-20*a*c)+-1//32*a^-4*x^-3*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(105*b^4+240*a^2*c^2+-448*a*c*b^2)+1//8*b*a^-3*x^-4*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(21*b^2+-68*a*c)+1//64*b*a^-5*x^-2*(b^2+-4*a*c)^-1*(a*x^2+b*x^3+c*x^4)^(1/2)*(315*b^4+1808*a^2*c^2+-1680*a*c*b^2))
@test integrate(x^m*(a*x+b*x^3+c*x^5), x) == :(a*x^(2+m)*(2+m)^-1+b*x^(4+m)*(4+m)^-1+c*x^(6+m)*(6+m)^-1)
@test integrate(x^2*(a*x+b*x^3+c*x^5), x) == :(1//4*a*x^4+1//6*b*x^6+1//8*c*x^8)
@test integrate(x*(a*x+b*x^3+c*x^5), x) == :(1//3*a*x^3+1//5*b*x^5+1//7*c*x^7)
@test integrate(a*x+b*x^3+c*x^5, x) == :((1/2)*a*x^2+1//4*b*x^4+1//6*c*x^6)
@test integrate(x^-1*(a*x+b*x^3+c*x^5), x) == :(a*x+1//3*b*x^3+1//5*c*x^5)
@test integrate(x^-2*(a*x+b*x^3+c*x^5), x) == :(a*log(x)+(1/2)*b*x^2+1//4*c*x^4)
@test integrate(x^-3*(a*x+b*x^3+c*x^5), x) == :(b*x+-1*a*x^-1+1//3*c*x^3)
@test integrate(x^m*(a*x+b*x^3+c*x^5)^2, x) == :(a^2*x^(3+m)*(3+m)^-1+c^2*x^(11+m)*(11+m)^-1+x^(7+m)*(7+m)^-1*(b^2+2*a*c)+2*a*b*x^(5+m)*(5+m)^-1+2*b*c*x^(9+m)*(9+m)^-1)
@test integrate(x^2*(a*x+b*x^3+c*x^5)^2, x) == :(1//5*a^2*x^5+1//9*x^9*(b^2+2*a*c)+1//13*c^2*x^13+2//7*a*b*x^7+2//11*b*c*x^11)
@test integrate(x*(a*x+b*x^3+c*x^5)^2, x) == :(1//4*a^2*x^4+1//8*x^8*(b^2+2*a*c)+1//12*c^2*x^12+1//3*a*b*x^6+1//5*b*c*x^10)
@test integrate((a*x+b*x^3+c*x^5)^2, x) == :(1//3*a^2*x^3+1//7*x^7*(b^2+2*a*c)+1//11*c^2*x^11+2//5*a*b*x^5+2//9*b*c*x^9)
@test integrate(x^-1*(a*x+b*x^3+c*x^5)^2, x) == :((1/2)*a^2*x^2+1//6*x^6*(b^2+2*a*c)+1//10*c^2*x^10+(1/2)*a*b*x^4+1//4*b*c*x^8)
@test integrate(x^-2*(a*x+b*x^3+c*x^5)^2, x) == :(x*a^2+1//5*x^5*(b^2+2*a*c)+1//9*c^2*x^9+2//3*a*b*x^3+2//7*b*c*x^7)
@test integrate(x^8*(a*x+b*x^3+c*x^5)^-1, x) == :(1//4*c^-1*x^4+-1//2*b*c^-2*x^2+1//4*c^-3*(b^2+-1*a*c)*log(a+b*x^2+c*x^4)+(1/2)*b*c^-3*(b^2+-4*a*c)^-1//2*(b^2+-3*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^7*(a*x+b*x^3+c*x^5)^-1, x) == :(1//3*c^-1*x^3+-1*b*x*c^-2+(1/2)*2^(1/2)*c^-5//2*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-1*a*c+b*(b^2+-4*a*c)^-1//2*(b^2+-3*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+(1/2)*2^(1/2)*c^-5//2*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-1*a*c+-1*b*(b^2+-4*a*c)^-1//2*(b^2+-3*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^6*(a*x+b*x^3+c*x^5)^-1, x) == :((1/2)*c^-1*x^2+-1//4*b*c^-2*log(a+b*x^2+c*x^4)+-1//2*c^-2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^5*(a*x+b*x^3+c*x^5)^-1, x) == :(x*c^-1+-1//2*2^(1/2)*c^-3//2*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b+(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+-1//2*2^(1/2)*c^-3//2*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b+-1*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^4*(a*x+b*x^3+c*x^5)^-1, x) == :(1//4*c^-1*log(a+b*x^2+c*x^4)+(1/2)*b*c^-1*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^3*(a*x+b*x^3+c*x^5)^-1, x) == :((1/2)*2^(1/2)*c^-1//2*(b+(b^2+-4*a*c)^(1/2))^(1/2)*(b^2+-4*a*c)^-1//2*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+-1//2*2^(1/2)*c^-1//2*(b+-1*(b^2+-4*a*c)^(1/2))^(1/2)*(b^2+-4*a*c)^-1//2*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^2*(a*x+b*x^3+c*x^5)^-1, x) == :(-1*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x*(a*x+b*x^3+c*x^5)^-1, x) == :(2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1//2*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2)+-1*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1//2*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate((a*x+b*x^3+c*x^5)^-1, x) == :(a^-1*log(x)+-1//4*a^-1*log(a+b*x^2+c*x^4)+(1/2)*b*a^-1*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^-1*(a*x+b*x^3+c*x^5)^-1, x) == :(-1*a^-1*x^-1+-1//2*2^(1/2)*a^-1*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(1+-1*b*(b^2+-4*a*c)^-1//2)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+-1//2*2^(1/2)*a^-1*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(1+b*(b^2+-4*a*c)^-1//2)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^-2*(a*x+b*x^3+c*x^5)^-1, x) == :(-1//2*a^-1*x^-2+-1*b*a^-2*log(x)+1//4*b*a^-2*log(a+b*x^2+c*x^4)+-1//2*a^-2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^11*(a*x+b*x^3+c*x^5)^-2, x) == :(-1//2*b*c^-3*log(a+b*x^2+c*x^4)+c^-2*x^2*(b^2+-4*a*c)^-1*(b^2+-3*a*c)+x^6*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+-1*c^-3*(b^2+-4*a*c)^-3//2*(b^4+6*a^2*c^2+-6*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2))+-1//2*b*c^-1*x^4*(b^2+-4*a*c)^-1)
@test integrate(x^10*(a*x+b*x^3+c*x^5)^-2, x) == :(x^5*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+(1/2)*x*c^-2*(b^2+-4*a*c)^-1*(3*b^2+-10*a*c)+-1//2*b*c^-1*x^3*(b^2+-4*a*c)^-1+-1//4*2^(1/2)*c^-5//2*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1*(3*b^3+(b^2+-4*a*c)^-1//2*(3*b^4+20*a^2*c^2+-19*a*c*b^2)+-13*a*b*c)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+-1//4*2^(1/2)*c^-5//2*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1*(3*b^3+-1*(b^2+-4*a*c)^-1//2*(3*b^4+20*a^2*c^2+-19*a*c*b^2)+-13*a*b*c)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^9*(a*x+b*x^3+c*x^5)^-2, x) == :(1//4*c^-2*log(a+b*x^2+c*x^4)+x^4*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+-1//2*b*c^-1*x^2*(b^2+-4*a*c)^-1+(1/2)*b*c^-2*(b^2+-4*a*c)^-3//2*(b^2+-6*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^8*(a*x+b*x^3+c*x^5)^-2, x) == :(x^3*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+-1//2*b*x*c^-1*(b^2+-4*a*c)^-1+1//4*2^(1/2)*c^-3//2*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1*(b^2+-6*a*c+b*(b^2+-4*a*c)^-1//2*(b^2+-8*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+1//4*2^(1/2)*c^-3//2*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1*(b^2+-6*a*c+-1*b*(b^2+-4*a*c)^-1//2*(b^2+-8*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^7*(a*x+b*x^3+c*x^5)^-2, x) == :(2*a*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2))+x^2*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2))
@test integrate(x^6*(a*x+b*x^3+c*x^5)^-2, x) == :(x*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+1//4*2^(1/2)*c^-1//2*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(b^2+b*(b^2+-4*a*c)^(1/2)+4*a*c)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+1//4*2^(1/2)*c^-1//2*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-1*(b+-1*(b^2+-4*a*c)^-1//2*(b^2+4*a*c))*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^5*(a*x+b*x^3+c*x^5)^-2, x) == :((2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(2a+b*x^2)+-1*b*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^4*(a*x+b*x^3+c*x^5)^-2, x) == :(-1*x*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b+2*c*x^2)+(1/2)*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(-1*(b^2+-4*a*c)^(1/2)+2b)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2)+-1//2*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*((b^2+-4*a*c)^(1/2)+2b)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^3*(a*x+b*x^3+c*x^5)^-2, x) == :(-1*(2*b^2+-8*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b+2*c*x^2)+2*c*(b^2+-4*a*c)^-3//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate(x^2*(a*x+b*x^3+c*x^5)^-2, x) == :((1/2)*x*a^-1*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2)+-1//4*2^(1/2)*a^-1*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(b^2+-1*b*(b^2+-4*a*c)^(1/2)+-12*a*c)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+1//4*2^(1/2)*a^-1*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(b^2+b*(b^2+-4*a*c)^(1/2)+-12*a*c)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x*(a*x+b*x^3+c*x^5)^-2, x) == :(a^-2*log(x)+-1//4*a^-2*log(a+b*x^2+c*x^4)+(1/2)*a^-1*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2)+(1/2)*b*a^-2*(b^2+-4*a*c)^-3//2*(b^2+-6*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2)))
@test integrate((a*x+b*x^3+c*x^5)^-2, x) == :(-1//2*a^-2*x^-1*(b^2+-4*a*c)^-1*(3*b^2+-10*a*c)+(1/2)*a^-1*x^-1*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2)+-1//4*2^(1/2)*a^-2*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(3*b^3+(b^2+-4*a*c)^(1/2)*(3*b^2+-10*a*c)+-16*a*b*c)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2)+1//4*2^(1/2)*a^-2*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(3*b^3+-1*(b^2+-4*a*c)^(1/2)*(3*b^2+-10*a*c)+-16*a*b*c)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^-1*(a*x+b*x^3+c*x^5)^-2, x) == :((1/2)*b*a^-3*log(a+b*x^2+c*x^4)+-2*b*a^-3*log(x)+-1*a^-3*(b^2+-4*a*c)^-3//2*(b^4+6*a^2*c^2+-6*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2))+-1*a^-2*x^-2*(b^2+-4*a*c)^-1*(b^2+-3*a*c)+(1/2)*a^-1*x^-2*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2))
@test integrate(x^-2*(a*x+b*x^3+c*x^5)^-2, x) == :(-1//6*a^-2*x^-3*(b^2+-4*a*c)^-1*(5*b^2+-14*a*c)+(1/2)*b*a^-3*x^-1*(b^2+-4*a*c)^-1*(5*b^2+-19*a*c)+(1/2)*a^-1*x^-3*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2)+-1//4*2^(1/2)*a^-3*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(5*b^4+28*a^2*c^2+-1*b*(b^2+-4*a*c)^(1/2)*(5*b^2+-19*a*c)+-29*a*c*b^2)*arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2)+1//4*2^(1/2)*a^-3*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2*(b^2+-4*a*c)^-3//2*(5*b^4+28*a^2*c^2+b*(b^2+-4*a*c)^(1/2)*(5*b^2+-19*a*c)+-29*a*c*b^2)*arctan(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2))
@test integrate(x^-3*(a*x+b*x^3+c*x^5)^-2, x) == :(a^-4*(3*b^2+-2*a*c)*log(x)+-1//4*a^-4*(3*b^2+-2*a*c)*log(a+b*x^2+c*x^4)+-1//4*a^-2*x^-4*(b^2+-4*a*c)^-1*(3*b^2+-8*a*c)+(1/2)*b*a^-4*(b^2+-4*a*c)^-3//2*(3*b^4+30*a^2*c^2+-20*a*c*b^2)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*x^2))+(1/2)*b*a^-3*x^-2*(b^2+-4*a*c)^-1*(3*b^2+-11*a*c)+(1/2)*a^-1*x^-4*(b^2+-4*a*c)^-1*(a+b*x^2+c*x^4)^-1*(b^2+-2*a*c+b*c*x^2))
@test integrate(x^3//2*(a*x+b*x^3+c*x^5)^(1/2), x) == :(1//15*c^-1*x^(1/2)*(a*x+b*x^3+c*x^5)^(1/2)*(b+3*c*x^2)+-1//15*c^-3//2*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(a*x+b*x^3+c*x^5)^-1//2*(2*b^2+-6*a*c)*(a+b*x^2+c*x^4)+-1//30*a^1//4*c^-7//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(2*b^2+-6*a*c+b*a^(1/2)*c^(1/2))*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+2//15*a^1//4*c^-7//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(b^2+-3*a*c)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^(1/2)*(a*x+b*x^3+c*x^5)^(1/2), x) == :(1//8*c^-1*x^-1//2*(a*x+b*x^3+c*x^5)^(1/2)*(b+2*c*x^2)+-1//16*c^-3//2*x^(1/2)*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-4*a*c)*arctanh((1/2)*c^-1//2*(a+b*x^2+c*x^4)^-1//2*(b+2*c*x^2)))
@test integrate(x^-1//2*(a*x+b*x^3+c*x^5)^(1/2), x) == :(1//3*x^(1/2)*(a*x+b*x^3+c*x^5)^(1/2)+1//3*b*c^-1//2*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a+b*x^2+c*x^4)+-1//3*b*a^1//4*c^-3//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+1//6*a^1//4*c^-3//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(b+2*a^(1/2)*c^(1/2))*(a^(1/2)+c^(1/2)*x^2)*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^-3//2*(a*x+b*x^3+c*x^5)^(1/2), x) == :((1/2)*x^-1//2*(a*x+b*x^3+c*x^5)^(1/2)+-1//2*a^(1/2)*x^(1/2)*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*arctanh((1/2)*a^-1//2*(a+b*x^2+c*x^4)^-1//2*(2a+b*x^2))+1//4*b*c^-1//2*x^(1/2)*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x^2+c*x^4)^-1//2*(b+2*c*x^2)))
@test integrate(x^3//2*(a*x+b*x^3+c*x^5)^3//2, x) == :(-1//640*c^-2*x^3//2*(a*x+b*x^3+c*x^5)^(1/2)*(b*(5*b^2+-4*a*c)+4*c*x^2*(5*b^2+-16*a*c))+1//80*c^-1*x^(1/2)*(a*x+b*x^3+c*x^5)^3//2*(3b+8*c*x^2)+1//1280*c^-3*x^-1//2*(a*x+b*x^3+c*x^5)^(1/2)*(15*b^4+128*a^2*c^2+-100*a*c*b^2)+-3//512*b*c^-7//2*x^(1/2)*(b^2+-4*a*c)^2*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x^2+c*x^4)^-1//2*(b+2*c*x^2)))
@test integrate(x^(1/2)*(a*x+b*x^3+c*x^5)^3//2, x) == :(-1//315*c^-2*x^(1/2)*(a*x+b*x^3+c*x^5)^(1/2)*(b*(4*b^2+-9*a*c)+6*c*x^2*(2*b^2+-7*a*c))+1//63*c^-1*x^-1//2*(a*x+b*x^3+c*x^5)^3//2*(3b+7*c*x^2)+1//315*c^-5//2*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a+b*x^2+c*x^4)*(8*b^4+84*a^2*c^2+-57*a*c*b^2)+-1//315*a^1//4*c^-11//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(8*b^4+84*a^2*c^2+-57*a*c*b^2)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+1//630*a^1//4*c^-11//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(8*b^4+84*a^2*c^2+-57*a*c*b^2+4*b*a^(1/2)*c^(1/2)*(b^2+-6*a*c))*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^-1//2*(a*x+b*x^3+c*x^5)^3//2, x) == :(1//16*c^-1*x^-3//2*(a*x+b*x^3+c*x^5)^3//2*(b+2*c*x^2)+-1//128*c^-2*x^-1//2*(a*x+b*x^3+c*x^5)^(1/2)*(b+2*c*x^2)*(3*b^2+-12*a*c)+3//256*c^-5//2*x^(1/2)*(b^2+-4*a*c)^2*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x^2+c*x^4)^-1//2*(b+2*c*x^2)))
@test integrate(x^-3//2*(a*x+b*x^3+c*x^5)^3//2, x) == :(1//7*x^-1//2*(a*x+b*x^3+c*x^5)^3//2+1//35*c^-1*x^(1/2)*(a*x+b*x^3+c*x^5)^(1/2)*(b^2+10*a*c+3*b*c*x^2)+-2//35*b*c^-3//2*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-8*a*c)*(a+b*x^2+c*x^4)+-1//70*a^1//4*c^-7//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(2*b*(b^2+-8*a*c)+a^(1/2)*c^(1/2)*(b^2+-20*a*c))*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+2//35*b*a^1//4*c^-7//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(b^2+-8*a*c)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^3//2*(a*x+b*x^3+c*x^5)^-1//2, x) == :((1/2)*c^-1//2*x^(1/2)*(a+b*x^2+c*x^4)^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*arctanh((1/2)*c^-1//2*(a+b*x^2+c*x^4)^-1//2*(b+2*c*x^2)))
@test integrate(x^(1/2)*(a*x+b*x^3+c*x^5)^-1//2, x) == :((1/2)*a^-1//4*c^-1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^-1//2*(a*x+b*x^3+c*x^5)^-1//2, x) == :(-1//2*a^-1//2*arctanh((1/2)*a^-1//2*x^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(2a+b*x^2)))
@test integrate(x^-3//2*(a*x+b*x^3+c*x^5)^-1//2, x) == :(-1*a^-1*x^-3//2*(a*x+b*x^3+c*x^5)^(1/2)+a^-1*c^(1/2)*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a+b*x^2+c*x^4)+(1/2)*a^-3//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+-1*a^-3//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^3//2*(a*x+b*x^3+c*x^5)^-3//2, x) == :(a^-1*x^3//2*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-2*a*c+b*c*x^2)+-1*b*a^-1*c^(1/2)*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a+b*x^2+c*x^4)+-1//2*a^-3//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(b+-2*a^(1/2)*c^(1/2))^-1*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+b*a^-3//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^(1/2)*(a*x+b*x^3+c*x^5)^-3//2, x) == :(-1//2*a^-3//2*arctanh((1/2)*a^-1//2*x^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(2a+b*x^2))+a^-1*x^(1/2)*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-2*a*c+b*c*x^2))
@test integrate(x^-1//2*(a*x+b*x^3+c*x^5)^-3//2, x) == :(a^-1*x^-1//2*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-2*a*c+b*c*x^2)+-1*a^-2*x^-3//2*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^(1/2)*(2*b^2+-6*a*c)+2*a^-2*c^(1/2)*x^3//2*(a^(1/2)+c^(1/2)*x^2)^-1*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-3*a*c)*(a+b*x^2+c*x^4)+(1/2)*a^-7//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(2*b^2+-6*a*c+b*a^(1/2)*c^(1/2))*Elliptic.F(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2)+-2*a^-7//4*c^1//4*x^(1/2)*((a^(1/2)+c^(1/2)*x^2)^-2*(a+b*x^2+c*x^4))^(1/2)*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(a^(1/2)+c^(1/2)*x^2)*(b^2+-3*a*c)*Elliptic.E(2*arctan(x*a^-1//4*c^1//4),1/2+-1//4*b*a^-1//2*c^-1//2))
@test integrate(x^-3//2*(a*x+b*x^3+c*x^5)^-3//2, x) == :(3//4*b*a^-5//2*arctanh((1/2)*a^-1//2*x^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(2a+b*x^2))+a^-1*x^-3//2*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^-1//2*(b^2+-2*a*c+b*c*x^2)+-1//2*a^-2*x^-5//2*(b^2+-4*a*c)^-1*(a*x+b*x^3+c*x^5)^(1/2)*(3*b^2+-8*a*c))
@test integrate(x^(-3//2+3//2*n)*(a*x^(-1+n)+b*x^n+c*x^(1+n))^-3//2, x) == :(-2*x^(-1//2+(1/2)*n)*(b^2+-4*a*c)^-1*(a*x^(-1+n)+b*x^n+c*x^(1+n))^-1//2*(b+2*c*x))
@test integrate((x^6+-3*x^4+3*x^2)^-1//2, x) == :(-1//6*3^(1/2)*arctanh(1//6*x*3^(1/2)*(x^6+-3*x^4+3*x^2)^-1//2*(6+-3*x^2)))
@test integrate((x^2*(3+x^4+-3*x^2))^-1//2, x) == :(-1//6*3^(1/2)*arctanh(1//6*x*3^(1/2)*(x^6+-3*x^4+3*x^2)^-1//2*(6+-3*x^2)))
@test integrate((1+-1*(1+-1*x^2)^3)^-1//2, x) == :(-1//6*3^(1/2)*arctanh(1//6*x*3^(1/2)*(x^6+-3*x^4+3*x^2)^-1//2*(6+-3*x^2)))
@test integrate((x^6+-3*x^4+3*x^2)^(1/2), x) == :(-1//8*x^-1*(x^6+-3*x^4+3*x^2)^(1/2)*(3+-2*x^2)+-3//16*x^-1*(3+x^4+-3*x^2)^-1//2*(x^6+-3*x^4+3*x^2)^(1/2)*arcsinh(1//3*3^(1/2)*(3+-2*x^2)))
@test integrate((x^2*(3+x^4+-3*x^2))^(1/2), x) == :(-1//8*x^-1*(x^6+-3*x^4+3*x^2)^(1/2)*(3+-2*x^2)+-3//16*x^-1*(3+x^4+-3*x^2)^-1//2*(x^6+-3*x^4+3*x^2)^(1/2)*arcsinh(1//3*3^(1/2)*(3+-2*x^2)))
@test integrate((1+-1*(1+-1*x^2)^3)^(1/2), x) == :(-1//8*x^-1*(x^6+-3*x^4+3*x^2)^(1/2)*(3+-2*x^2)+-3//16*x^-1*(3+x^4+-3*x^2)^-1//2*(x^6+-3*x^4+3*x^2)^(1/2)*arcsinh(1//3*3^(1/2)*(3+-2*x^2)))
@test integrate(x^-1*(a+b*x+c*x^2)^-1//2, x) == :(-1*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*x+c*x^2)^-1//2*(2a+b*x)))
@test integrate((x^2*(a+b*x+c*x^2))^-1//2, x) == :(-1*a^-1//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^3+c*x^4)^-1//2*(2a+b*x)))
@test integrate(x^-1//2*(x*(a+b*x+c*x^2))^-1//2, x) == :(-1*a^-1//2*arctanh((1/2)*a^-1//2*x^(1/2)*(a*x+b*x^2+c*x^3)^-1//2*(2a+b*x)))
@test integrate(x^(1/2)*(x^3*(a+b*x+c*x^2))^-1//2, x) == :(-1*a^-1//2*arctanh((1/2)*a^-1//2*x^3//2*(a*x^3+b*x^4+c*x^5)^-1//2*(2a+b*x)))
@test integrate(x^-1*(a+b*x^2+c*x^4)^-1//2, x) == :(-1//2*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*x^2+c*x^4)^-1//2*(2a+b*x^2)))
@test integrate((x^2*(a+b*x^2+c*x^4))^-1//2, x) == :(-1//2*a^-1//2*arctanh((1/2)*x*a^-1//2*(a*x^2+b*x^4+c*x^6)^-1//2*(2a+b*x^2)))
@test integrate(x^-1//2*(x*(a+b*x^2+c*x^4))^-1//2, x) == :(-1//2*a^-1//2*arctanh((1/2)*a^-1//2*x^(1/2)*(a*x+b*x^3+c*x^5)^-1//2*(2a+b*x^2)))
@test integrate(x^(1/2)*(x^3*(a+b*x^2+c*x^4))^-1//2, x) == :(-1//2*a^-1//2*arctanh((1/2)*a^-1//2*x^3//2*(a*x^3+b*x^5+c*x^7)^-1//2*(2a+b*x^2)))
@test integrate(x^-1*(3+x^4+-3*x^2)^-1//2, x) == :(-1//6*3^(1/2)*arctanh((1/2)*3^(1/2)*(3+x^4+-3*x^2)^-1//2*(2+-1*x^2)))
@test integrate((x^2*(3+x^4+-3*x^2))^-1//2, x) == :(-1//6*3^(1/2)*arctanh(1//6*x*3^(1/2)*(x^6+-3*x^4+3*x^2)^-1//2*(6+-3*x^2)))
@test integrate(x^-1//2*(x*(3+x^2+-3x))^-1//2, x) == :(-1//3*3^(1/2)*arctanh((1/2)*3^(1/2)*x^(1/2)*(x^3+-3*x^2+3x)^-1//2*(2+-1x)))
@test integrate(x^(-1+(1/2)*q)*(a*x^q+b*x^n+c*x^(-1q+2n))^-1//2, x) == :(-1*a^-1//2*(n+-1q)^-1*arctanh((1/2)*a^-1//2*x^((1/2)*q)*(a*x^q+b*x^n+c*x^(-1q+2n))^-1//2*(2a+b*x^(n+-1q))))

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.4 Improper/1.2.4.2 (d x)^m (a x^q+b x^n+c x^(2 n-q))^p.jl


text:
# https://github.com/stan-dev/blob/master/basic_estimators/normal_mixture.stan

# normal mixture, unknown proportion and means, known variance
# p(y|mu,theta) = theta * Normal(y|mu[1],1) + (1-theta) * Normal(y|mu[2],1);

const simplenormalmixturemodel = "
data {
  int<lower=0>  N;
  real y[N];
}
parameters {
  real<lower=0,upper=1> theta;
  real mu[2];
}
transformed parameters {
  real log_theta;
  real log_one_minus_theta;

  log_theta <- log(theta);
  log_one_minus_theta <- log(1.0 - theta);
}
model {
  theta ~ uniform(0,1); // equivalently, ~ beta(1,1);
  for (k in 1:2)
    mu[k] ~ normal(0,10);
  for (n in 1:N)
    increment_log_prob(
      log_sum_exp(log_theta + normal_log(y[n],mu[1],1.0),
                  log_one_minus_theta + normal_log(y[n],mu[2],1.0)));
}
"


text:
# TODO
# - find need allocations in attractive case
# - decimate allocations here
# - I guess maybe it's additional Diagonal's etc?

################################################################################
### Type and general utility
################################################################################

struct BlockDiagonal{T, N, AT <: AbstractMatrix{T}} <: AbstractMatrix{T}
    blocks::NTuple{N, AT}

    function BlockDiagonal(blocks::NTuple{N, AT}) where {T, N, AT <: AbstractMatrix{T}}
        n, m = size(first(blocks))
        n == m || throw(ErrorException("Each block must be square."))
        all(b -> size(b) == (n, n), blocks) || throw(ErrorException(
            "All Blocks must have the same size."
        ))
        # One matrix with zeros is faster than multiple 3x3 matrices on my
        # machine (avx2). This will most likely depend on AVX version...
        # if n < 4
        #     full = zeros(T, n*N, n*N)
        #     for i in 1:N
        #         @views copyto!(full[(i-1)*n+1 : i*n, (i-1)*n+1 : i*n], blocks[i])
        #     end
        #     # new{T, 1, AT}((full,))
        #     # Is this bad? - Yes, because the rest of the code doesn't know about Matrix returns
        #     return full
        # else
            return new{T, N, AT}(blocks)
        # end
    end
end

BlockDiagonal(blocks::AT...) where {T, AT <: AbstractMatrix{T}} = BlockDiagonal(blocks)
function BlockDiagonal{T, N, AT}(::UndefInitializer, n, m) where {T, N, AT}
    n == m || throw(ErrorException("Each block must be square."))
    (n, m) = divrem(n, N)
    m == 0 || throw(ErrorException(
        "The size of the Matrix is incompatible with the blocking ($n % $N != 0)"
    ))
    BlockDiagonal([AT(undef, n, n) for _ in 1:N]...)
end
function BlockDiagonal{T, N, AT}(::UniformScaling, n, m) where {T, N, AT}
    n == m || throw(ErrorException("Each block must be square."))
    (n, m) = divrem(n, N)
    m == 0 || throw(ErrorException(
        "The size of the Matrix is incompatible with the blocking ($n % $N != 0)"
    ))
    BlockDiagonal([AT(I, n, n) for _ in 1:N]...)
end

function Base.Matrix(B::BlockDiagonal{T, N}) where {T, N}
    n = size(B.blocks[1], 1)
    M = N * n
    output = zeros(T, M, M)
    for i in 1:N
        @views copyto!(output[(i-1)*n+1 : i*n, (i-1)*n+1 : i*n], B.blocks[i])
    end
    output
end

function Base.checkbounds(B::BlockDiagonal{T, N}, i, j) where {T, N}
    @inbounds n = size(B.blocks[1], 1)
    (1 ≤ i ≤ N*n) && (1 ≤ j ≤ N*n)
end

function Base.getindex(B::BlockDiagonal{T, N}, i, j) where {T, N}
    @boundscheck checkbounds(B, i, j)
    @inbounds n = size(B.blocks[1], 1)

    b, k = divrem(i-1, n)
    _b, l = divrem(j-1, n)

    if b != _b
        return zero(T)
    else
        @inbounds return B.blocks[b+1][k+1, l+1]
    end
end

function Base.copyto!(B1::BlockDiagonal{T, N}, B2::BlockDiagonal{T, N}) where {T, N}
    @inbounds for i in 1:N
        copyto!(B1.blocks[i], B2.blocks[i])
    end
    return B1
end
function Base.copyto!(B::BlockDiagonal{T, N}, ::UniformScaling) where {T, N}
    @inbounds for i in 1:N
        copyto!(B.blocks[i], I)
    end
    return B
end
function Base.copyto!(B::BlockDiagonal{T, N}, D::Diagonal{T}) where {T, N}
    @inbounds n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        b = B.blocks[i]
        offset = (i-1)*n
        @turbo for j in 1:n, k in 1:n
            b[j, k] = ifelse(j == k, D.diag[offset+j], zero(T))
        end
    end
    return B
end

# for Base.show
function Base.size(B::BlockDiagonal{T, N}) where {T, N}
    @inbounds n = size(B.blocks[1], 1)
    (n*N, n*N)
end

# Needed for tests
Base.copy(B::BlockDiagonal) = deepcopy(B)

# Needed for stack building (shouldn't be performance critical)
function Base.:(*)(x::T, B::BlockDiagonal{T, N, AT}) where {T<:Number, N, AT<:AbstractMatrix{T}}
    BlockDiagonal(map(block -> x * block, B.blocks)...)
end
Base.:(*)(B::BlockDiagonal{T, N, AT}, x::T) where {T<:Number, N, AT<:AbstractMatrix{T}} = x * B
function Base.:(*)(B1::BlockDiagonal{T, N, AT}, B2::BlockDiagonal{T, N, AT}) where {
        T<:Number, N, AT<:AbstractMatrix{T}
    }
    BlockDiagonal(map(*, B1.blocks, B2.blocks)...)
end

function Base.exp(B::BlockDiagonal)
    BlockDiagonal(map(block -> exp(block), B.blocks)...)
end

function Base.log(B::BlockDiagonal)
    BlockDiagonal(map(block -> log(block), B.blocks)...)
end

function LinearAlgebra.det(B::BlockDiagonal{T}) where {T}
    output = T(1)
    for b in B.blocks
        output *= det(b)
    end
    output
end

# I thought this would be needed for greens(k, l), but it's not?
function LinearAlgebra.transpose!(A::BlockDiagonal{T, N}, B::BlockDiagonal{T, N}) where {T, N}
    @inbounds for i in 1:N
        transpose!(A.blocks[i], B.blocks[i])
    end
    A
end
# This however is used for greens(k, l)
function LinearAlgebra.rmul!(B::BlockDiagonal{T, N}, f::Number) where {T, N}
    @inbounds for i in 1:N
        rmul!(B.blocks[i], f)
    end
    B
end

################################################################################
### AVX-powered matrix multiplications
################################################################################

function vmul!(C::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, B::BlockDiagonal{T, N}) where {T, N}
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A.blocks[i], B.blocks[i])
    end
end
function vmul!(C::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, B::Diagonal{T}) where {T<:Real, N}
    # Assuming correct size
    @inbounds n = size(A.blocks[1], 1)
    @inbounds for i in 1:N
        c = C.blocks[i]
        a = A.blocks[i]
        @turbo for k in 1:n, l in 1:n
            c[k,l] = a[k,l] * B.diag[(i-1)*n + l]
        end
    end
end
function vmul!(C::BlockDiagonal{T, N}, A::Diagonal{T}, B::BlockDiagonal{T, N}) where {T<:Real, N}
    # Assuming correct size
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        c = C.blocks[i]
        b = B.blocks[i]
        @turbo for k in 1:n, l in 1:n
            c[k,l] = A.diag[(i-1)*n + k] * b[k,l]
        end
    end
end
function vmul!(C::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, X::Adjoint{T}) where {T<:Real, N}
    B = X.parent
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        a = A.blocks[i]
        b = B.blocks[i]
        c = C.blocks[i]
        @turbo for k in 1:n, l in 1:n
            Ckl = zero(eltype(c))
            for m in 1:n
                Ckl += a[k,m] * b[l, m]
            end
            c[k,l] = Ckl
        end
    end
end
function vmul!(C::BlockDiagonal{T, N}, X::Adjoint{T}, B::BlockDiagonal{T, N}) where {T<:Real, N}
    A = X.parent
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        a = A.blocks[i]
        b = B.blocks[i]
        c = C.blocks[i]
        @turbo for k in 1:n, l in 1:n
            Ckl = zero(eltype(c))
            for m in 1:n
                Ckl += a[m,k] * b[m,l]
            end
            c[k,l] = Ckl
        end
    end
end
function vmul!(C::BD, X1::Adjoint{T, BD}, X2::Adjoint{T, BD}) where {T <: Real, N, BD <: BlockDiagonal{T, N}}
    A = X1.parent
    B = X2.parent
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        a = A.blocks[i]
        b = B.blocks[i]
        c = C.blocks[i]
        @turbo for k in 1:n, l in 1:n
            Ckl = zero(eltype(c))
            for m in 1:n
                Ckl += a[m,k] * b[l, m]
            end
            c[k,l] = Ckl
        end
    end
end
function rvmul!(A::BlockDiagonal{T, N}, B::Diagonal) where {T, N}
    # Assuming correct size
    @inbounds n = size(A.blocks[1], 1)
    @inbounds for i in 1:N
        @views rvmul!(A.blocks[i], Diagonal(B.diag[(i-1)*n+1 : i*n]))
    end
end
function lvmul!(A::Diagonal, B::BlockDiagonal{T, N}) where {T, N}
    # Assuming correct size
    @inbounds n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        @views lvmul!(Diagonal(A.diag[(i-1)*n+1 : i*n]), B.blocks[i])
    end
end

# used in greens(k, l)
function rvadd!(A::BlockDiagonal{T, N}, B::BlockDiagonal{T, N}) where {T <: Real, N}
    @inbounds for i in 1:N
        a = A.blocks[i]
        b = B.blocks[i]
        @turbo for j in axes(a, 1), k in axes(a, 2)
            a[j, k] = a[j, k] + b[j, k]
        end
    end
end
# used in equal time greens
function rvadd!(B::BlockDiagonal{T, N}, D::Diagonal{T}) where {T<:Real, N}
    # Assuming correct size
    @inbounds n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        a = B.blocks[i]
        offset = (i-1) * n
        @turbo for j in 1:n
            a[j, j] = a[j, j] + D.diag[j+offset]
        end
    end
end
# used in CombinedGreensIterator
function vsub!(O::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, ::UniformScaling) where {T<:Real, N}
    @inbounds n = size(O.blocks[1], 1)
    T1 = one(T)
    @inbounds for i in 1:N
        a = A.blocks[i]
        o = O.blocks[i]
        @turbo for j in 1:n, k in 1:n
            o[j, k] = a[j, k]
        end
        @turbo for j in 1:n
            o[j, j] -= T1
        end
    end
end

function rdivp!(A::BD, T::BD, O::BD, pivot) where {ET<:Real, N, BD <: BlockDiagonal{ET, N}}
    # assume Diagonal is ±1!
    @inbounds begin
        n = size(A.blocks[1], 1)

        # Apply pivot
        for b in 1:N
            o = O.blocks[b]
            a = A.blocks[b]
            t = T.blocks[b]
            offset = (b-1)*n

            for j in 1:n
                p = pivot[j + offset]
                @turbo for i in 1:n
                    o[i, j] = a[i, p]
                end
            end

            # do the rdiv
            # @turbo will segfault on `k in 1:0`, so pull out first loop
            invt11 = 1.0 / t[1, 1]
            @turbo for i in 1:n
                a[i, 1] = o[i, 1] * invt11
            end
            for j in 2:n
                invtjj = 1.0 / t[j, j]
                @turbo for i in 1:n
                    x = o[i, j]
                    for k in 1:j-1
                        x -= a[i, k] * t[k, j]
                    end
                    a[i, j] = x * invtjj
                end
            end
        end
    end
    A
end

################################################################################
### Fallbacks
################################################################################

function vmul!(C::BlockDiagonal{T, N}, A::Diagonal{T}, B::BlockDiagonal{T, N}) where {T, N}
    # Assuming correct size
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        @views vmul!(C.blocks[i], Diagonal(A.diag[(i-1)*n+1 : i*n]), B.blocks[i])
    end
end
function vmul!(C::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, B::Diagonal) where {T, N}
    # Assuming correct size
    @inbounds n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        @views vmul!(C.blocks[i], A.blocks[i], Diagonal(B.diag[(i-1)*n+1 : i*n]))
    end
end
function vmul!(C::BlockDiagonal{T, N}, A::BlockDiagonal{T, N}, X::Adjoint{T}) where {T, N}
    B = X.parent
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A.blocks[i], adjoint(B.blocks[i]))
    end
end
function vmul!(C::BlockDiagonal{T, N}, X::Adjoint{T}, B::BlockDiagonal{T, N}) where {T, N}
    A = X.parent
    @inbounds for i in 1:N
        vmul!(C.blocks[i], adjoint(A.blocks[i]), B.blocks[i])
    end
end
function rvadd!(B::BlockDiagonal{T, N}, D::Diagonal{T}) where {T, N}
    # Assuming correct size
    @inbounds n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        @views rvadd!(B.blocks[i], Diagonal(D.diag[(i-1)*n+1 : i*n]))
    end
end
function rvadd!(A::BlockDiagonal{T, N}, B::BlockDiagonal{T, N}) where {T, N}
    @inbounds for i in 1:N
        a = A.blocks[i]
        b = B.blocks[i]
        rvadd!(a, b)
    end
end
function rdivp!(A::BD, T::BD, O::BD, pivot) where {ET, N, BD <: BlockDiagonal{ET, N}}
    # Assuming correct size
    @inbounds n = size(T.blocks[1], 1)
    @inbounds for i in 1:N
        @views rdivp!(A.blocks[i], T.blocks[i], O.blocks[i], pivot[(i-1)*n+1 : i*n])
    end
    A
end

################################################################################
### UDT
################################################################################

function udt_AVX_pivot!(
        U::BlockDiagonal{Float64, N, AT},
        D::AbstractArray{Float64, 1},
        input::BlockDiagonal{Float64, N, AT},
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{Float64, 1} = Vector{Float64}(undef, length(D)),
        apply_pivot::Val = Val(true)
    ) where {N, AT <: AbstractMatrix{Float64}}
    # Assuming correct size
    @inbounds n = size(input.blocks[1], 1)
    @inbounds for i in 1:N
        @views udt_AVX_pivot!(
            U.blocks[i], D[(i-1)*n+1 : i*n], input.blocks[i],
            pivot[(i-1)*n+1 : i*n], temp[(i-1)*n+1 : i*n], apply_pivot
        )
        # pivot[(i-1)*n+1 : i*n] .+= (i-1)*n
    end
end
function udt_AVX_pivot!(
        U::BlockDiagonal{ComplexF64, N, AT},
        D::AbstractArray{Float64, 1},
        input::BlockDiagonal{ComplexF64, N, AT},
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{ComplexF64, 1} = Vector{ComplexF64}(undef, length(D)),
        apply_pivot::Val = Val(true)
    ) where {N, AT <: AbstractMatrix{ComplexF64}}
    # Assuming correct size
    @inbounds n = size(input.blocks[1], 1)
    @inbounds for i in 1:N
        @views udt_AVX_pivot!(
            U.blocks[i], D[(i-1)*n+1 : i*n], input.blocks[i],
            pivot[(i-1)*n+1 : i*n], temp[(i-1)*n+1 : i*n], apply_pivot
        )
        # pivot[(i-1)*n+1 : i*n] .+= (i-1)*n
    end
end


text:
@testset "planes" begin
    p = Point([1, 1])
    t1 = Polygon([[[0, 0, 0], [2, 0, 0], [1, 2, 2], [0, 0, 0]]])

    @test planepoint(p, t1) == 1.
end


text:
using DeformableRegistration: Distance, ImageProcessing, Transformation, Interpolation, regOptions
 using Base.Test
 using Logging

include("helpers/checkDerivative.jl")

@testset "Distance: distance derivatives - SSD" begin
    # create reference image
    refImg = createImage(100*rand(50,60))

    # measure distance and check derivative (nonparametric)
    centeredGrid = getCellCenteredGrid(refImg)
     centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
     D,dD,d2D = ssdDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true)
     Dfunc(x) = ssdDistance(refImg,refImg,x)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
     @test checkErrorDecay(errquad)

    options = regOptions()
     options.matrixFree = true;
     D,dD,d2D = ssdDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
     Dfunc(x) = ssdDistance(refImg,refImg,x)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data)
     @test checkErrorDecay(errquad)

    # measure distance and check derivative (parametic)
    centeredGrid = getCellCenteredGrid(refImg)
     options = DeformableRegistration.regOptions()
     options.parametricOnly=true
     evaluationPoint = [1,0,0,0,1,0]+0.1*rand(6)
     D,dD,d2D = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
     Dfunc(p) = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint)
     @test checkErrorDecay(errquad)

    options.matrixFree = true;
     D,dD,d2D = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
     Dfunc(p) = ssdDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
     errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint)
     @test checkErrorDecay(errquad)
end

@testset "Distance: distance derivatives - NGF - parametric" begin
    refImg = createImage(100*rand(50,60))

    # measure distance and check derivative (nonparametric)
    centeredGrid = getCellCenteredGrid(refImg)
    centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
    D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true)
    Dfunc(x) = ngfDistance(refImg,refImg,x)[1]
    errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=false)
    @test checkErrorDecay(errquad)

    evaluationPoint = [1,0,0,0,1,0.0]+0.1*rand(6)
    centeredGrid = getCellCenteredGrid(refImg)
    options = DeformableRegistration.regOptions()
    options.edgeParameterR = 0.1
    options.edgeParameterT = 0.1
    options.matrixFree = false
    options.parametricOnly = true
    D,dD,d2D = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
    Dfunc(p) = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
    errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint,doPlot=true)
    @test checkErrorDecay(errquad)

    evaluationPoint = [1,0,0,0,1,0.0]+0.1*rand(6)
    centeredGrid = getCellCenteredGrid(refImg)
    options = DeformableRegistration.regOptions()
    options.edgeParameterR = 0.1
    options.edgeParameterT = 0.1
    options.matrixFree = true
    options.parametricOnly = true
    D,dD,d2D = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,evaluationPoint).data,doDerivative=true,doHessian=true,options=options)
    Dfunc(p) = ngfDistance(refImg,refImg,transformGridAffine(centeredGrid,p).data,options=options)[1]
    errlin,errquad = checkDerivative(Dfunc,dD',evaluationPoint,doPlot=true)
    @test checkErrorDecay(errquad)
end

@testset "Distance: distance derivatives - NGF - nonparametric" begin
    refImg = createImage(100*rand(50,60))

    ##
    centeredGrid = getCellCenteredGrid(refImg)
    centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
    options = DeformableRegistration.regOptions()
    options.edgeParameterR = 0.1
    options.edgeParameterT = 0.1
    options.matrixFree = false
    options.parametricOnly = false
    D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
    Dfunc(grid) = ngfDistance(refImg,refImg,grid,options=options)[1]
    errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=true)
    @test checkErrorDecay(errquad)

    centeredGrid = getCellCenteredGrid(refImg)
    centeredGrid.data[:] = centeredGrid.data[:] + 0.3 *rand(size(centeredGrid.data))
    options = DeformableRegistration.regOptions()
    options.edgeParameterR = 0.1
    options.edgeParameterT = 0.1
    options.matrixFree = true
    options.parametricOnly = false
    D,dD,d2D = ngfDistance(refImg,refImg,centeredGrid.data,doDerivative=true,doHessian=true,options=options)
    Dfunc(grid) = ngfDistance(refImg,refImg,grid,options=options)[1]
    errlin,errquad = checkDerivative(Dfunc,dD',centeredGrid.data,doPlot=true)
    @test checkErrorDecay(errquad)

end

@testset "Distance: ngf epsilon estimation" begin
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    refImg = loadImage(testimage)
    for r = refImg.data
         r = r + 0.01*rand()
    end
    Logging.warn("Skipping ")
    @test_skip estimateNGFEpsilon(refImg,cutoffPercent=80)[1] >= 0 # seems to run forever
end


text:
module BiodiversityObservationNetworks
    using SimpleSDMLayers
    using Distributions
    using NeutralLandscapes
    using Random
    using HaltonSequences
    using StatsBase
    using Base: @kwdef

    include("types.jl")
    export SpatialSampler

    include("sampler.jl")
    export rand, rand!

    include("balancedacceptance.jl")
    export BalancedAcceptance

    include("_helpers.jl")
    export makesdm, makeoccurrence, makebon, makeenv

end

# makebon(makesdm())


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MiniZinc"

version = v"2.6.2"

sources = [
    GitSource(
        "https://github.com/MiniZinc/libminizinc.git",
        "a56602765b4294b796c063664733b28f5a663af7",
    ),
    DirectorySource("./bundled"),
]

script = raw"""
cd $WORKSPACE/srcdir/libminizinc

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fixes.patch

# Patch for MinGW toolchain
find .. -type f -exec sed -i 's/Windows.h/windows.h/g' {} +

mkdir -p build
cd build

cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ..

if [[ "${target}" == *-linux-* ]]; then
        make -j ${nproc}
else
    if [[ "${target}" == *-mingw* ]]; then
        cmake --build . --config Release
    else
        cmake --build . --config Release --parallel
    fi
fi
make install
"""

products = [
    ExecutableProduct("minizinc", :minizinc),
]

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# TODO(odow): fix build issues on Windows
platforms = filter(!Sys.iswindows, platforms)

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version = v"6",
    julia_compat = "1.6",
)


text:
function save_noise!(W::NoiseApproximation)

end

function interpolate!(W::NoiseApproximation,t)
  W.Z != nothing ? z = W.source2.sol(t) : z = nothing
  W.source1.sol(t),z
end

function interpolate!(out1,out2,W::NoiseApproximation,t)
  W.source1.sol(out1,t), W.source2.sol(out2,t)
end

function calculate_step!(W::NoiseApproximation,dt,u,p)
  W.dt = dt
  t = W.curt+dt

  add_tstop!(W.source1,t)
  step!(W.source1)
  if W.Z != nothing
    add_tstop!(W.source2,t)
    step!(W.source2)
  end

  if isinplace(W)
    W.dW .= W.source1.u .- W.curW
    if W.Z != nothing
      W.dZ .= W.source2.u .- W.curZ
    end
  else
    W.dW = W.source1.u - W.curW
    if W.Z != nothing
      W.dZ = W.source2.u - W.curZ
    end
  end
  return nothing
end

function accept_step!(W::NoiseApproximation,dt,u,p,setup_next=true)
  if isinplace(W)
    W.curW .+= W.dW
  else
    W.curW += W.dW
  end
  W.curt += W.dt
  if W.Z != nothing
    if isinplace(W)
      W.curZ .+= W.dZ
    else
      W.curZ += W.dZ
    end
  end

  W.dt = dt #dtpropose
  if setup_next
    calculate_step!(W,dt,u,p)
  end
  return nothing
end

function reject_step!(W::NoiseApproximation,dtnew,u,p)
  W.dt = dtnew
  if isinplace(W)
    W.source1(W.dW,W.curt+dtnew)
    W.dW .-= W.curW
    if W.Z != nothing
      W.source2(W.dZ,W.curt+dtnew)
      W.dZ .-= W.curZ
    end
  else
    W.dW = W.source1(W.curt+dtnew) - W.curW
    if W.Z != nothing
      W.dZ = W.source2(W.curt+dtnew) - W.curZ
    end
  end
  return nothing
end

function setup_next_step!(W::NoiseApproximation,u,p)
  calculate_step!(W,W.dt,u,p)
  return nothing
end