text:
import KernelDensity: UnivariateKDE
import ..UncertainValues: UncertainScalarKDE, UncertainScalarPopulation

"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedMeanWeightedResampling)

Transform index-irregularly spaced uncertain data onto a regular index-grid and
take the mean of the values in each bin. Resamples the data points in `x`
according to `resampling.weights`.

Distributions in each index bin are obtained by resampling all index values in `x`
`resampling.n` times, in proportions obeying `resampling.weights` and mapping those
index draws to the bins. Simultaneously, the values in `x` are resampled and placed
in the corresponding bins. Finally, the mean in each bin is calculated. In total,
`length(x)*resampling.n` draws are distributed among the bins to form the final mean
estimate.

Returns a vector of mean values, one for each bin.

Assumes that the points in `x` are independent.

## Example

```julia
vars = (1, 2)
npts, tstep = 100, 10
d_xind = Uniform(2.5, 15.5)
d_yind = Uniform(2.5, 15.5)
d_xval = Uniform(0.01, 0.2)
d_yval = Uniform(0.01, 0.2)

X, Y = example_uncertain_indexvalue_datasets(ar1_unidir(c_xy = 0.5), npts, vars, tstep = tstep,
d_xind = d_xind, d_yind = d_yind,
d_xval = d_xval, d_yval = d_yval);

n_draws = 10000 # draws per uncertain value
time_grid = 0:50:1000
wts = Weights(rand(length(X))) # some random weights

# Resample both X and Y so that they are both at the same time indices,
# and take the mean of each bin.
resampled_dataset = resample(X, BinnedMeanWeightedResampling(time_grid, wts, n_draws))
resampled_dataset = resample(Y, BinnedMeanWeightedResampling(time_grid, wts, n_draws))
```
"""
function resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedMeanWeightedResampling)

    # Represent entire dataset as a weighted population and sample from that
    pop_inds = UncertainValue(x.indices.indices, resampling.weights)
    pop_vals = UncertainValue(x.values.values, resampling.weights)

    # Pre-allocate an array representing each bin, and an array keeping track
    # of the values falling in that bin.
    n_bins = length(resampling.left_bin_edges) - 1
    bin_sums = fill(0.0, n_bins)
    bin_sums_n_entries = fill(0.0, n_bins)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    for (j, (pop_idx, pop_val)) in enumerate(zip(pop_inds, pop_vals))

        # Sample the j-th idx-value pair `resampling.n` times and
        # accumulate the values in the correct bin sum. Also keep
        # track of how many values there are in each bin.
        resample!(idxs, pop_idx)
        resample!(vals, pop_val)

        @inbounds for i = 1:resampling.n
            arr_idx = ceil(Int, (idxs[i] - mini) / s)

            # Because the indices of `x` are uncertain values
            # with potentially infinite support, we need to check
            # that the value falls inside the grid
            if 0 < arr_idx <= n_bins
                bin_sums[arr_idx] += vals[i]
                bin_sums_n_entries[arr_idx] += 1.0
            end
        end
    end

    # Return bin averages (entries with 0s are represented as NaNs)
    bin_avgs = bin_sums ./ bin_sums_n_entries
    bin_avgs[isapprox.(bin_avgs, 0.0)] .= NaN
    return bin_avgs
end

"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedMeanResampling)

Transform index-irregularly spaced uncertain data onto a regular index-grid and
take the mean of the values in each bin.

Distributions in each index bin are obtained by resampling all index values in `x`
`resampling.n` times, and mapping those index draws to the bins. Simultaneously, the
values in `x` are resampled and placed in the corresponding bins. Finally, the mean
in each bin is calculated. In total, `length(x)*resampling.n` draws are distributed
among the bins to form the final mean estimate.

Returns a vector of mean values, one for each bin.

Assumes that the points in `x` are independent.

## Example

```julia
vars = (1, 2)
npts, tstep = 100, 10
d_xind = Uniform(2.5, 15.5)
d_yind = Uniform(2.5, 15.5)
d_xval = Uniform(0.01, 0.2)
d_yval = Uniform(0.01, 0.2)

X, Y = example_uncertain_indexvalue_datasets(ar1_unidir(c_xy = 0.5), npts, vars, tstep = tstep,
    d_xind = d_xind, d_yind = d_yind,
    d_xval = d_xval, d_yval = d_yval);

n_draws = 10000 # draws per uncertain value
time_grid = 0:50:1000

# Resample both X and Y so that they are both at the same time indices,
# and take the mean of each bin.
resampled_dataset = resample(X, BinnedMeanResampling(time_grid, n_draws))
resampled_dataset = resample(Y, BinnedMeanResampling(time_grid, n_draws))
```
"""
function resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedMeanResampling)

    # Constrain
    # Represent entire dataset as an equally-weighted population and sample from that
    pop_inds = UncertainValue(x.indices.indices, [1 for x in x.indices])
    pop_vals = UncertainValue(x.values.values, [1 for x in x.values])

    # Pre-allocate an array representing each bin, and an array keeping track
    # of the values falling in that bin.
    n_bins = length(resampling.left_bin_edges) - 1
    bin_sums = fill(0.0, n_bins)
    bin_sums_n_entries = fill(0.0, n_bins)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    for (j, (pop_idx, pop_val)) in enumerate(zip(pop_inds, pop_vals))

        # Sample the j-th idx-value pair `resampling.n` times and
        # accumulate the values in the correct bin sum. Also keep
        # track of how many values there are in each bin.
        resample!(idxs, pop_idx)
        resample!(vals, pop_val)

        @inbounds for i = 1:resampling.n
            arr_idx = ceil(Int, (idxs[i] - mini) / s)

            # Because the indices of `x` are uncertain values
            # with potentially infinite support, we need to check
            # that the value falls inside the grid
            if 0 < arr_idx <= n_bins
                bin_sums[arr_idx] += vals[i]
                bin_sums_n_entries[arr_idx] += 1.0
            end
        end
    end

    # Return bin averages (entries with 0s are represented as NaNs)
    bin_avgs = bin_sums ./ bin_sums_n_entries
    bin_avgs[isapprox.(bin_avgs, 0.0)] .= NaN
    return bin_avgs
end

"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedResampling{UncertainScalarKDE};
        bin_repr::Symbol = :kde, nan_threshold = 0.0)

Transform index-irregularly spaced uncertain data onto a regular index-grid.
Distributions in each index bin are obtained by resampling all index values
in `x` `resampling.n` times, and mapping those index draws to the bins.
Simultaneously, the values in `x` are resampled and placed in the corresponding
bins. In total, `length(x)*resampling.n` draws are distributed among the bins to form
the final KDEs.

Returns an `UncertainIndexValueDataset`. The distribution of values in the `i`-th bin
is approximated by a kernel density estimate (KDE) over the draws falling in the
`i`-th bin.

Assumes that the points in `x` are independent.

## Example

```julia
vars = (1, 2)
npts, tstep = 100, 10
d_xind = Uniform(2.5, 15.5)
d_yind = Uniform(2.5, 15.5)
d_xval = Uniform(0.01, 0.2)
d_yval = Uniform(0.01, 0.2)

X, Y = example_uncertain_indexvalue_datasets(ar1_unidir(c_xy = 0.5), npts, vars, tstep = tstep,
    d_xind = d_xind, d_yind = d_yind,
    d_xval = d_xval, d_yval = d_yval);

n_draws = 10000 # draws per uncertain value
time_grid = 0:50:1000
resampling = BinnedResampling(time_grid, n_draws)

# Resample both X and Y so that they are both at the same time indices.
resampled_dataset = resample(X, resampling)
resampled_dataset = resample(Y, resampling)
```
"""
function resample(x::AbstractUncertainIndexValueDataset,
        resampling::BinnedResampling{UncertainScalarKDE};
        nan_threshold = 0.0)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    perminds = zeros(Int, resampling.n)
    sorted_idxs = fill(NaN, resampling.n)
    sorted_vals = fill(NaN, resampling.n)

    bin_edges = resampling.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, resampling.n*n_bins) for bv in binvecs]

    @inbounds for (j, (idx, val)) in enumerate(zip(x.indices, x.values))
        # Resample the j-th index and j-th value
        resample!(idxs, idx)
        resample!(vals, val)

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    # Estimate distributions in each bin by kernel density estimation
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarKDE}}(undef, n_bins)

    for i in 1:n_bins
        if length(binvecs[i]) > nan_threshold
            estimated_value_dists[i] = UncertainValue(binvecs[i])
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(UncertainValue.(bin_edges[1:end-1] .+ step(bin_edges)/2))
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

function resample(x::AbstractUncertainIndexValueDataset,
        resampling::BinnedResampling{UncertainScalarPopulation};
        nan_threshold = 0.0)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    perminds = zeros(Int, resampling.n)
    sorted_idxs = fill(NaN, resampling.n)
    sorted_vals = fill(NaN, resampling.n)

    bin_edges = resampling.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, resampling.n*n_bins) for bv in binvecs]

    @inbounds for (j, (idx, val)) in enumerate(zip(x.indices, x.values))
        # Resample the j-th index and j-th value
        resample!(idxs, idx)
        resample!(vals, val)

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    # Estimate distributions in each bin by kernel density estimation
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarPopulation}}(undef, n_bins)

    for i in 1:n_bins
        if length(binvecs[i]) > nan_threshold
            L = length(binvecs[i])
            estimated_value_dists[i] = UncertainValue(binvecs[i], repeat([1 / L], L))
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(UncertainValue.(bin_edges[1:end-1] .+ step(bin_edges)/2))
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

function resample(x::AbstractUncertainIndexValueDataset,
        resampling::BinnedResampling{RawValues};
        nan_threshold = 0.0)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    perminds = zeros(Int, resampling.n)
    sorted_idxs = fill(NaN, resampling.n)
    sorted_vals = fill(NaN, resampling.n)

    bin_edges = resampling.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, resampling.n*n_bins) for bv in binvecs]

    @inbounds for (j, (idx, val)) in enumerate(zip(x.indices, x.values))
        # Resample the j-th index and j-th value
        resample!(idxs, idx)
        resample!(vals, val)

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    bin_centers = bin_edges[1:end-1] .+ step(bin_edges)/2

    return bin_centers, binvecs
end

"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedWeightedResampling;
        nan_threshold = 0.0)

Transform index-irregularly spaced uncertain data onto a regular index-grid.

Distributions in each index bin are obtained by resampling all index values
in `x` `resampling.n` times, sampled according to probabilities `resampling.weights`,
and mapping those index draws to the bins. Simultaneously, the values in `x` are
resampled and placed in the corresponding bins. In total, `length(x)*resampling.n` draws
are distributed among the bins to form the final KDEs.

Returns an `UncertainIndexValueDataset`. The distribution of values in the `i`-th bin
is approximated by a kernel density estimate (KDE) over the draws falling in the
`i`-th bin.

Assumes that the points in `x` are independent.

## Example

```julia
vars = (1, 2)
npts, tstep = 100, 10
d_xind = Uniform(2.5, 15.5)
d_yind = Uniform(2.5, 15.5)
d_xval = Uniform(0.01, 0.2)
d_yval = Uniform(0.01, 0.2)

X, Y = example_uncertain_indexvalue_datasets(ar1_unidir(c_xy = 0.5), npts, vars, tstep = tstep,
    d_xind = d_xind, d_yind = d_yind,
    d_xval = d_xval, d_yval = d_yval);

left_bin_edges = 0:50:1000
n_draws = 10000
wts = Weights(rand(length(X)))
resampling = BinnedWeightedResampling(left_bin_edges, wts, 10)
resampled_dataset = resample(X, resampling)
```
"""
function resample(x::AbstractUncertainIndexValueDataset, resampling::BinnedWeightedResampling;
        nan_threshold = 0.0)

    pop_inds = UncertainValue(x.indices.indices, resampling.weights)
    pop_vals = UncertainValue(x.values.values, resampling.weights)

    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, resampling.n)
    vals = fill(NaN, resampling.n)

    perminds = zeros(Int, resampling.n)
    sorted_idxs = fill(NaN, resampling.n)
    sorted_vals = fill(NaN, resampling.n)

    bin_edges = resampling.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(resampling.left_bin_edges)
    s = step(resampling.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, resampling.n*n_bins) for bv in binvecs]

    @inbounds for (j, (idx, val)) in enumerate(zip(pop_inds, pop_vals))
        # Resample the j-th index and j-th value
        resample!(idxs, idx)
        resample!(vals, val)

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    # Estimate distributions in each bin by kernel density estimation
    estimated_value_dists = Vector{AbstractUncertainValue}(undef, n_bins)

    for i in 1:n_bins
        if length(binvecs[i]) > nan_threshold
            estimated_value_dists[i] = UncertainValue(binvecs[i])
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(UncertainValue.(bin_edges[1:end-1] .+ step(bin_edges)/2))
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end


text:
module Intern
using ..Meta
export @interned, @get!  # shouldn't need to export @get!

const new_interned = gensym("new_interned")

function replace_new(ex::Expr)
    if ex.head === :quote; return ex; end

    args = Any[replace_new(arg) for arg in ex.args]
    if (ex.head === :call) && args[1] == :new
        return Expr(:call, new_interned, args[2:end]...)
    end
    Expr(ex.head, args...)
end
replace_new(ex) = ex

macro interned(ex)
    code_interned(ex)
end
function code_interned(ex)
    @expect is_expr(ex, :type, 3)
    imm, typesig, typebody = ex.args
    typeex = (is_expr(typesig, :(<:), 2) ? typesig.args[1] : typesig)
    typename =   (is_expr(typeex, :curly) ? typeex.args[1] : typeex)::Symbol
    typeparams = (is_expr(typeex, :curly) ? typeex.args[2:end] : [])

    fields, types, sigs, defs = Symbol[], [], [], []
    needs_default_constructor = true
    for def in typebody.args
        if isa(def, Symbol)
            push!(fields, def)
            push!(types,  quot(Any))
            push!(sigs, def)
        elseif is_expr(def, :(::), 2)
            push!(fields, def.args[1])
            push!(types,  def.args[2])
            push!(sigs, def)
        elseif is_fdef(def)
            sig, body = split_fdef(def)
            if sig.args[1] === typename  # constructor
                needs_default_constructor = false
                body = replace_new(body)
                def = :($sig=$body)
            end
        end
        push!(defs, def)
    end

    objects = ObjectIdDict()
    push!(defs, :($new_interned{$(typeparams...)}($(sigs...)) =
                  @get!($(quot(objects)),($(fields...),), new{$(typeparams...)}($(fields...))) ))
    if needs_default_constructor
        push!(defs, :( $typename($(sigs...)) = $new_interned($(fields...)) ))
    end
    esc(Expr(:type, imm, typesig, Expr(:block, defs...)))
end

end

ext: jl


text:
module Startup_utils

export install_julia, read_instance_list, read_tasks_list

using DataFrames
using Master.Types

function install_julia(userhost, folder="~", sshpass=nothing, sshfile=nothing)
    if sshpass != nothing
        spawn(`sshpass -p $sshpass scp julia.tar.gz $userhost:$folder;
               sshpass -p $sshpass ssh $userhost
                   "cd $folder;
                   tar -xzf julia.tar.gz;
                   rm julia.tar.gz"`)
    elseif sshfile != nothing
        spawn(`scp -i $sshfile julia.tar.gz $userhost:$folder;
               ssh -i $sshfile $userhost
                   "cd $folder;
                   tar -xzf julia.tar.gz;
                   rm julia.tar.gz"`)
    else
        spawn(`scp julia.tar.gz $userhost:$folder;
               ssh $userhost
                   "cd $folder;
                   tar -xzf julia.tar.gz;
                   rm julia.tar.gz"`)
    end
end

function read_instance_list(
        file::AbstractString="instances.csv")::Array{Instance_info,1}
    instances = Array{Instance_info,1}()
    try
        table = readtable(file)
        run(`clear`)
        for i in 1:length(table[1]) # for each entry in the table
            println("Enter password for $(table[1][i]):")
            password = chomp(readline())
            run(`clear`)
            instance = Instance_info(table[1][i], password, table[2][i])
            push!(instances, instance)
        end
    catch e
        println(e)
        warn("The instances list failed to load.")
        instances = Array{Instance_info,1}()
    end
    return instances
end

function read_tasks_list(
        file::AbstractString="tasks.csv")::Array{Requested_task,1}
    tasks = Array{Requested_task,1}()
    try
        table = readtable(file)
        for i in 1:length(table[1]) # for each entry in the table
            task = Requested_task(table[1][i], table[2][i])
            push!(tasks, task)
        end
    catch e
        println(e)
        warn("The tasks list failed to load.")
        tasks = Array{Requested_task,1}()
    end
    return tasks
end

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
println("==============flag! fpump.jl======================================")

"""
    generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)

Generate a mip using the linear constraints of the original model
TODO: This can include quadratic constraints when the mip_solver supports them
Minimize the distance to nlp_sol and avoid using solutions inside the tabu list
"""
function generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)
    println("================flag! fpump ============= generate_mip")

    mip_optimizer = m.mip_solver.optimizer_constructor()
    mip_model = Model(m.mip_solver)
    @variable(mip_model, mx[i = 1:m.num_var],
        binary = m.var_type[i] == :Bin,
        integer = m.var_type[i] == :Int
    )

    # only add bounds for non binary variables
    for i=1:m.num_var
        if m.var_type[i] != :Bin
           @constraint(mip_model, m.l_var[i] <= mx[i] <= m.u_var[i])
        end

        if m.var_type[i] == :Bin && (m.l_var[i] > 0 || m.u_var[i] < 1)
            # must be 1
            if m.l_var[i] > 0
                @constraint(mip_model, mx[i] == 1)
            else # or 0
                @constraint(mip_model, mx[i] == 0)
            end
        end
    end

    backend = JuMP.backend(mip_model);

    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end

    @variable(mip_model, mabsx[i=1:m.num_disc_var] >= 0)
    for i=1:m.num_disc_var
        vi = m.disc2var_idx[i]
        @constraint(mip_model, mabsx[i] >= mx[vi]-nlp_sol[vi])
        @constraint(mip_model, mabsx[i] >= -mx[vi]+nlp_sol[vi])
    end

    # How long is the tabu list
    num_sols = 0
    for i=1:tabu_list.length
        if !isnan(tabu_list.sols[i][1])
            num_sols += 1
        else
            break
        end
    end

    # If there solutions in the tabu list => avoid them
    if num_sols > 0
        @variable(mip_model, z1[j=1:m.num_disc_var,k=1:num_sols], Bin)
        @variable(mip_model, z2[j=1:m.num_disc_var,k=1:num_sols], Bin)
        v = tabu_list.sols
        for k=1:num_sols, j=1:m.num_disc_var
            i = m.disc2var_idx[j]
            lbi = m.l_var[i] > typemin(Int64) ? m.l_var[i] : typemin(Int64)
            ubi = m.u_var[i] < typemax(Int64) ? m.u_var[i] : typemax(Int64)
            @constraint(mip_model, z1[j,k]+z2[j,k] <= 1)
            @constraint(mip_model, (lbi - v[k][j])*z1[j,k]+z2[j,k]+v[k][j] <= mx[i])
            @constraint(mip_model, mx[i] <= v[k][j] - z1[j,k] + (ubi-v[k][j])*z2[j,k])
        end
        for k=1:num_sols
            @constraint(mip_model, sum(z1[j,k]+z2[j,k] for j=1:m.num_disc_var) >= 1)
        end
    end

    @objective(mip_model, Min, sum(mabsx[i] for i=1:m.num_disc_var))

    # Break the mip solver if it takes too long or throw a warning when this option isn't available
    current_time = time()-start_fpump
    time_left = m.options.feasibility_pump_time_limit-current_time
    time_left < 0 && (time_left = 1.0)

    # set time limit if supported
    old_time_limit = set_time_limit!(mip_optimizer, time_left)

    status, backend = optimize_get_status_backend(mip_model)

    # reset time limit
    set_time_limit!(mip_optimizer, old_time_limit)

    obj_val = NaN
    values = fill(NaN,m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        obj_val = JuMP.objective_value(mip_model)

        # round mip values
        values = JuMP.value.(mx)
        for i=1:m.num_disc_var
            vi = m.disc2var_idx[i]
            values[vi] = round(values[vi])
        end
    end

    return status, values, obj_val
end

"""
    generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=false)

Generates the original nlp but changes the objective to minimize the distance to the mip solution
"""
function generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=false)
    nlp_optimizer = m.nl_solver.optimizer_constructor()
    nlp_model = Model(m.nl_solver)
    @variable(nlp_model, m.l_var[i] <= nx[i=1:m.num_var] <= m.u_var[i])
    if random_start
        restart_values = generate_random_restart(m)
        for i=1:m.num_var
            JuMP.set_start_value(nx[i], restart_values[i])
        end
    else
        JuMP.set_start_value.(nx[1:m.num_var],mip_sol)
    end

    # add all constraints
    backend = JuMP.backend(nlp_model);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    qlc = optimizer.quadratic_le_constraints
    qgc = optimizer.quadratic_ge_constraints
    qec = optimizer.quadratic_eq_constraints
    for constr_type in [llc, lgc, lec, qlc, qgc, qec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    for i in 1:m.num_nl_constr
        constr_expr = MOI.constraint_expr(optimizer.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, nlp_model)
        JuMP.add_NL_constraint(nlp_model, constr_expr)
    end

    @objective(nlp_model, Min, sum((nx[m.disc2var_idx[i]]-mip_sol[m.disc2var_idx[i]])^2 for i=1:m.num_disc_var))

    current_time = time()-start_fpump
    time_left = m.options.feasibility_pump_time_limit-current_time
    time_left < 0 && (time_left = 1.0)

    # set time limit if supported
    old_time_limit = set_time_limit!(nlp_optimizer, time_left)

    status, backend = optimize_get_status_backend(nlp_model)

    set_time_limit!(nlp_optimizer, old_time_limit)

    nlp_obj = NaN
    nlp_sol = fill(NaN,m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        nlp_obj = JuMP.objective_value(nlp_model)
        nlp_sol = JuMP.value.(nx)
    end

    Base.finalize(backend)

    return status, nlp_sol, nlp_obj
end

"""
    generate_real_nlp(optimizer, m, sol; random_start=false)

Generate the original nlp and get the objective for that
"""
function generate_real_nlp(optimizer, m, sol; random_start=false)
    if m.num_var == m.num_disc_var
        if optimizer.nlp_data.has_objective
            nlp_obj = MOI.eval_objective(optimizer.nlp_data.evaluator, sol)
        else
            nlp_obj = evaluate_objective(optimizer, m, sol)
        end
        status = MOI.OPTIMAL
        return status, sol, nlp_obj
    end

    rmodel = Model(m.nl_solver)

    @variable(rmodel, m.l_var[i] <= rx[i=1:m.num_var] <= m.u_var[i])
    if random_start
        restart_values = generate_random_restart(m)
        for i=1:m.num_var
            if m.var_type[i] == :Cont
                JuMP.set_start_value(rx[i], restart_values[i])
            end # discrete will be fixed anyway
        end
    end
    for i=1:m.num_disc_var
        vi = m.disc2var_idx[i]
        JuMP.fix(rx[vi], sol[vi]; force=true)
    end

    # define the objective function
    # TODO check whether it is supported
    if optimizer.nlp_data.has_objective
        obj_expr = MOI.objective_expr(optimizer.nlp_data.evaluator)
        obj_expr = expr_dereferencing(obj_expr, rmodel)
        JuMP.set_NL_objective(rmodel, optimizer.sense, obj_expr)
    else
        MOI.set(rmodel, MOI.ObjectiveFunction{typeof(optimizer.objective)}(), optimizer.objective)
        MOI.set(rmodel, MOI.ObjectiveSense(), optimizer.sense)
    end

    backend = JuMP.backend(rmodel);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    qlc = optimizer.quadratic_le_constraints
    qgc = optimizer.quadratic_ge_constraints
    qec = optimizer.quadratic_eq_constraints
    for constr_type in [llc, lgc, lec, qlc, qgc, qec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    for i in 1:m.num_nl_constr
        constr_expr = MOI.constraint_expr(optimizer.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, rmodel)
        JuMP.add_NL_constraint(rmodel, constr_expr)
    end

    status, backend = optimize_get_status_backend(rmodel)

    obj_val = NaN
    real_sol = fill(NaN,m.num_var)
    if state_is_optimal(status; allow_almost=m.options.allow_almost_solved)
        obj_val = JuMP.objective_value(rmodel)
        real_sol = JuMP.value.(rx)
    end

    Base.finalize(backend)

    return status, real_sol, obj_val
end

"""
    add!(t::TabuList, m, sol)

Add a solution to the tabu list (includes only the discrete variables)
"""
function add!(t::TabuList, m, sol)
    println("================flag! fpump ============= add!")
    t.sols[t.pointer] = [sol[m.disc2var_idx[i]] for i=1:m.num_disc_var]
    t.pointer += 1
    if t.pointer > t.length
        t.pointer = 1
    end
end

function get_fp_table(mip_obj,nlp_obj,t, fields, field_chars, catol)
    "================flag! fpump ============= get_fp_table"
    ln = ""
    i = 1
    arr = []
    for f in fields
        val = ""
        if f == :MIPobj
            if isnan(mip_obj)
                val = "-"
            else
                digits = 4
                if mip_obj < 1e-4
                    digits = convert(Int,floor(log10(1/catol)))
                end
                val = string(round(mip_obj; digits=digits))
            end
        elseif f == :NLPobj
            if isnan(nlp_obj)
                val = "-"
            else
                digits = 4
                if nlp_obj < 1e-4
                    digits = convert(Int,floor(log10(1/catol)))
                end
                val = string(round(nlp_obj; digits=digits))
            end
        elseif f == :Time
            val = string(round(t; digits=1))
        end

        if length(val) > field_chars[i]
            # too long to display shouldn't happen normally but is better than error
            # if it happens
            val = "t.l."
        end

        padding = field_chars[i]-length(val)
        ln *= repeat(" ",trunc(Int, floor(padding/2)))
        ln *= val
        ln *= repeat(" ",trunc(Int, ceil(padding/2)))
        push!(arr,val)
        i += 1
    end
    return ln, arr
end

"""
    fpump(optimizer, m)

Run the feasibility pump
"""
function fpump(optimizer, m)
    Random.seed!(1)

    if are_type_correct(m.relaxation_solution, m.var_type, m.disc2var_idx, m.options.atol)
        return Incumbent(m.relaxation_objval, m.relaxation_solution, only_almost_solved(m.status))
    end

    start_fpump = time()
    println("==============flag!=========fp==========================")
    nlp_sol = m.relaxation_solution
    nlp_obj = 1 # should be not 0 for while
    c = 1
    tabu_list = TabuList()
    mip_sols = Dict{UInt64,Bool}()
    tabu_list.length = m.options.tabu_list_length
    tabu_list.pointer = 1
    tabu_list.sols = []
    for i=1:tabu_list.length
        push!(tabu_list.sols, NaN*ones(m.num_disc_var))
    end

    last_table_arr = []
    fields = []
    field_chars = []
    ps = m.options.log_levels
    # Print table init
    if check_print(ps,[:Table])
        fields, field_chars = [:MIPobj,:NLPobj,:Time], [20,20,5]
        print_table_header(fields,field_chars)
    end

    real_status = MOI.OPTIMIZE_NOT_CALLED
    fix = false
    nlp_status = :Error
    iscorrect = false
    tl = m.options.feasibility_pump_time_limit
    # the tolerance can be changed => current atol
    catol = m.options.atol
    atol_counter = 0
    while !are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol) && time()-start_fpump < tl &&
        time()-m.start_time < m.options.time_limit

        # generate a mip or just round if no linear constraints
        if m.num_l_constr > 0
            mip_status, mip_sol, mip_obj = generate_mip(optimizer, m, nlp_sol, tabu_list, start_fpump)
        else
            # if no linear constraints just round the discrete variables
            mip_obj = NaN
            mip_sol = copy(nlp_sol)
            mip_status = MOI.OPTIMAL
            for vi=1:m.num_disc_var
                vidx = m.disc2var_idx[vi]
                mip_sol[vidx] = round(mip_sol[vidx])
            end
        end
        if mip_status != MOI.OPTIMAL
            @warn "MIP couldn't be solved to optimality. Terminated with status: "*string(mip_status)
            break
        end

        # If a cycle is detected which wasn't able to prevent by the tabu list (maybe too short)
        if haskey(mip_sols, hash(mip_sol))
            @warn "Cycle detected"
            break
        end
        add!(tabu_list, m, mip_sol)
        mip_sols[hash(mip_sol)] = true

        nlp_status, nlp_sol, nlp_obj = generate_nlp(optimizer, m, mip_sol, start_fpump)
        if !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved)
            cnlpinf = 0
            while cnlpinf < m.options.num_resolve_nlp_feasibility_pump && !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved) &&
                time()-start_fpump < tl && time()-m.start_time < m.options.time_limit
                nlp_status, nlp_sol, nlp_obj = generate_nlp(optimizer, m, mip_sol, start_fpump; random_start=true)
                cnlpinf += 1
            end
            if !state_is_optimal(nlp_status; allow_almost=m.options.allow_almost_solved)
                @warn "NLP couldn't be solved to optimality"
                if check_print(ps,[:Table])
                    print_fp_table(mip_obj, NaN, time()-start_fpump, fields, field_chars, catol)
                end
                break
            end
        end

        if check_print(ps,[:Table])
            print_fp_table(mip_obj, nlp_obj, time()-start_fpump, fields, field_chars, catol)
        end

        # if the current tolerance was nearly reached 5 times
        # => If reasonable should be an option
        if atol_counter >= m.options.feasibility_pump_tolerance_counter
            catol *= 10
            @warn "FPump tolerance changed to: ",catol
            atol_counter = 0
        end

        # if the difference is near 0 => try to improve the obj by using the original obj
        # set atol for type correct to a low value as it is checked with real_nlp anyway
        if are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol*1000) || isapprox(nlp_obj, 0.0; atol=catol)
            real_status,real_sol, real_obj = generate_real_nlp(optimizer, m, mip_sol)
            cnlpinf = 0
            while cnlpinf < m.options.num_resolve_nlp_feasibility_pump && !state_is_optimal(real_status; allow_almost=m.options.allow_almost_solved) &&
                time()-start_fpump < tl && time()-m.start_time < m.options.time_limit
                real_status,real_sol, real_obj = generate_real_nlp(optimizer, m, mip_sol; random_start=true)
                cnlpinf += 1
            end
            if state_is_optimal(real_status) || (only_almost_solved(real_status) && m.options.allow_almost_solved_integral)
                if only_almost_solved(real_status)
                    @warn "Integral feasible point only almost solved. Disable with `allow_almost_solved_integral=false`"
                end
                nlp_obj = real_obj
                nlp_sol = real_sol
                iscorrect = true
                break
            elseif are_type_correct(nlp_sol, m.var_type, m.disc2var_idx, catol)
                real_status = MOI.LOCALLY_SOLVED
                nlp_obj = evaluate_objective(optimizer, m, nlp_sol)
                iscorrect = true
                @warn "Real objective wasn't solved to optimality"
                break
            end
        end
        if !isapprox(nlp_obj, 0.0; atol=catol) && isapprox(nlp_obj, 0.0; atol=10*catol)
            atol_counter += 1
        else
            atol_counter = 0
        end
        c += 1
    end

    if check_print(ps,[:Table])
        println()
    end

    if check_print(ps,[:Info])
        println("FP: ", time()-start_fpump, " s")
        println("FP: ", c == 1 ? "$c round" : "$c rounds")
    end
    m.fpump_info = Dict{Symbol,Float64}()
    m.fpump_info[:time] = time()-start_fpump
    m.fpump_info[:rounds] = c

    if iscorrect
        check_print(ps,[:Info]) && println("FP: Obj: ", nlp_obj)
        m.fpump_info[:obj] = nlp_obj
        m.fpump_info[:gap] = abs(m.relaxation_objval-nlp_obj)/abs(nlp_obj)
        return Incumbent(nlp_obj, nlp_sol, only_almost_solved(real_status))
    end

    m.fpump_info[:obj] = NaN
    m.fpump_info[:gap] = NaN
    check_print(ps,[:Info]) && println("FP: No integral solution found")
    return nothing
end


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

@test integrate(x^5*erf(b*x), x) == :(-5//16*b^-6*erf(b*x)+1//6*x^6*erf(b*x)+1//6*pi^-1//2*b^-1*x^5*exp(-1*b^2*x^2)+5//8*x*pi^-1//2*b^-5*exp(-1*b^2*x^2)+5//12*pi^-1//2*b^-3*x^3*exp(-1*b^2*x^2))
@test integrate(x^3*erf(b*x), x) == :(-3//16*b^-4*erf(b*x)+1//4*x^4*erf(b*x)+1//4*pi^-1//2*b^-1*x^3*exp(-1*b^2*x^2)+3//8*x*pi^-1//2*b^-3*exp(-1*b^2*x^2))
@test integrate(x*erf(b*x), x) == :((1/2)*x^2*erf(b*x)+-1//4*b^-2*erf(b*x)+(1/2)*x*pi^-1//2*b^-1*exp(-1*b^2*x^2))
@test integrate(x^-3*erf(b*x), x) == :(-1*b^2*erf(b*x)+-1//2*x^-2*erf(b*x)+-1*b*pi^-1//2*x^-1*exp(-1*b^2*x^2))
@test integrate(x^-5*erf(b*x), x) == :(-1//4*x^-4*erf(b*x)+1//3*b^4*erf(b*x)+-1//6*b*pi^-1//2*x^-3*exp(-1*b^2*x^2)+1//3*pi^-1//2*b^3*x^-1*exp(-1*b^2*x^2))
@test integrate(x^-7*erf(b*x), x) == :(-4//45*b^6*erf(b*x)+-1//6*x^-6*erf(b*x)+-4//45*pi^-1//2*b^5*x^-1*exp(-1*b^2*x^2)+-1//15*b*pi^-1//2*x^-5*exp(-1*b^2*x^2)+2//45*pi^-1//2*b^3*x^-3*exp(-1*b^2*x^2))
@test integrate(x^6*erf(b*x), x) == :(1//7*x^7*erf(b*x)+6//7*pi^-1//2*b^-7*exp(-1*b^2*x^2)+1//7*pi^-1//2*b^-1*x^6*exp(-1*b^2*x^2)+3//7*pi^-1//2*b^-3*x^4*exp(-1*b^2*x^2)+6//7*pi^-1//2*b^-5*x^2*exp(-1*b^2*x^2))
@test integrate(x^4*erf(b*x), x) == :(1//5*x^5*erf(b*x)+2//5*pi^-1//2*b^-5*exp(-1*b^2*x^2)+1//5*pi^-1//2*b^-1*x^4*exp(-1*b^2*x^2)+2//5*pi^-1//2*b^-3*x^2*exp(-1*b^2*x^2))
@test integrate(x^2*erf(b*x), x) == :(1//3*x^3*erf(b*x)+1//3*pi^-1//2*b^-3*exp(-1*b^2*x^2)+1//3*pi^-1//2*b^-1*x^2*exp(-1*b^2*x^2))
@test integrate(erf(b*x), x) == :(x*erf(b*x)+pi^-1//2*b^-1*exp(-1*b^2*x^2))
@test integrate(x^-2*erf(b*x), x) == :(-1*x^-1*erf(b*x)+b*pi^-1//2*expintegralei(-1*b^2*x^2))
@test integrate(x^-4*erf(b*x), x) == :(-1//3*x^-3*erf(b*x)+-1//3*pi^-1//2*b^3*expintegralei(-1*b^2*x^2)+-1//3*b*pi^-1//2*x^-2*exp(-1*b^2*x^2))
@test integrate(x^-6*erf(b*x), x) == :(-1//5*x^-5*erf(b*x)+1//10*pi^-1//2*b^5*expintegralei(-1*b^2*x^2)+-1//10*b*pi^-1//2*x^-4*exp(-1*b^2*x^2)+1//10*pi^-1//2*b^3*x^-2*exp(-1*b^2*x^2))
@test integrate((c+d*x)^3*erf(a+b*x), x) == :(-3//16*b^-4*d^3*erf(a+b*x)+1//4*d^-1*(c+d*x)^4*erf(a+b*x)+pi^-1//2*b^-4*(b*c+-1*a*d)^3*exp(-1*(a+b*x)^2)+-3//4*d*b^-4*(b*c+-1*a*d)^2*erf(a+b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*erf(a+b*x)+pi^-1//2*b^-4*d^2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+1//4*pi^-1//2*b^-4*d^3*(a+b*x)^3*exp(-1*(a+b*x)^2)+3//8*pi^-1//2*b^-4*d^3*(a+b*x)*exp(-1*(a+b*x)^2)+pi^-1//2*b^-4*d^2*(a+b*x)^2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+3//2*d*pi^-1//2*b^-4*(b*c+-1*a*d)^2*(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate((c+d*x)^2*erf(a+b*x), x) == :(1//3*d^-1*(c+d*x)^3*erf(a+b*x)+pi^-1//2*b^-3*(b*c+-1*a*d)^2*exp(-1*(a+b*x)^2)+-1//2*d*b^-3*(b*c+-1*a*d)*erf(a+b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*erf(a+b*x)+1//3*pi^-1//2*b^-3*d^2*exp(-1*(a+b*x)^2)+1//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*exp(-1*(a+b*x)^2)+d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*exp(-1*(a+b*x)^2))
@test integrate((c+d*x)*erf(a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*erf(a+b*x)+-1//4*d*b^-2*erf(a+b*x)+pi^-1//2*b^-2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*erf(a+b*x)+(1/2)*d*pi^-1//2*b^-2*(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erf(a+b*x), x) == :(pi^-1//2*b^-1*exp(-1*(a+b*x)^2)+b^-1*(a+b*x)*erf(a+b*x))
@test integrate(x^5*erf(b*x)^2, x) == :(-5//16*b^-6*erf(b*x)^2+1//6*x^6*erf(b*x)^2+11//12*pi^-1*b^-6*exp(-2*b^2*x^2)+1//6*pi^-1*b^-2*x^4*exp(-2*b^2*x^2)+7//12*pi^-1*b^-4*x^2*exp(-2*b^2*x^2)+1//3*pi^-1//2*b^-1*x^5*erf(b*x)*exp(-1*b^2*x^2)+5//4*x*pi^-1//2*b^-5*erf(b*x)*exp(-1*b^2*x^2)+5//6*pi^-1//2*b^-3*x^3*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^3*erf(b*x)^2, x) == :(-3//16*b^-4*erf(b*x)^2+1//4*x^4*erf(b*x)^2+(1/2)*pi^-1*b^-4*exp(-2*b^2*x^2)+1//4*pi^-1*b^-2*x^2*exp(-2*b^2*x^2)+(1/2)*pi^-1//2*b^-1*x^3*erf(b*x)*exp(-1*b^2*x^2)+3//4*x*pi^-1//2*b^-3*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x*erf(b*x)^2, x) == :((1/2)*x^2*erf(b*x)^2+-1//4*b^-2*erf(b*x)^2+(1/2)*pi^-1*b^-2*exp(-2*b^2*x^2)+x*pi^-1//2*b^-1*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-3*erf(b*x)^2, x) == :(-1*b^2*erf(b*x)^2+-1//2*x^-2*erf(b*x)^2+2*pi^-1*b^2*expintegralei(-2*b^2*x^2)+-2*b*pi^-1//2*x^-1*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-5*erf(b*x)^2, x) == :(-1//4*x^-4*erf(b*x)^2+1//3*b^4*erf(b*x)^2+-4//3*pi^-1*b^4*expintegralei(-2*b^2*x^2)+-1//3*pi^-1*b^2*x^-2*exp(-2*b^2*x^2)+-1//3*b*pi^-1//2*x^-3*erf(b*x)*exp(-1*b^2*x^2)+2//3*pi^-1//2*b^3*x^-1*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-7*erf(b*x)^2, x) == :(-4//45*b^6*erf(b*x)^2+-1//6*x^-6*erf(b*x)^2+28//45*pi^-1*b^6*expintegralei(-2*b^2*x^2)+-1//15*pi^-1*b^2*x^-4*exp(-2*b^2*x^2)+2//9*pi^-1*b^4*x^-2*exp(-2*b^2*x^2)+-8//45*pi^-1//2*b^5*x^-1*erf(b*x)*exp(-1*b^2*x^2)+-2//15*b*pi^-1//2*x^-5*erf(b*x)*exp(-1*b^2*x^2)+4//45*pi^-1//2*b^3*x^-3*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^4*erf(b*x)^2, x) == :(1//5*x^5*erf(b*x)^2+-43//80*2^(1/2)*pi^-1//2*b^-5*erf(b*x*2^(1/2))+1//5*pi^-1*b^-2*x^3*exp(-2*b^2*x^2)+4//5*pi^-1//2*b^-5*erf(b*x)*exp(-1*b^2*x^2)+11//20*x*pi^-1*b^-4*exp(-2*b^2*x^2)+2//5*pi^-1//2*b^-1*x^4*erf(b*x)*exp(-1*b^2*x^2)+4//5*pi^-1//2*b^-3*x^2*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^2*erf(b*x)^2, x) == :(1//3*x^3*erf(b*x)^2+-5//12*2^(1/2)*pi^-1//2*b^-3*erf(b*x*2^(1/2))+1//3*x*pi^-1*b^-2*exp(-2*b^2*x^2)+2//3*pi^-1//2*b^-3*erf(b*x)*exp(-1*b^2*x^2)+2//3*pi^-1//2*b^-1*x^2*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(erf(b*x)^2, x) == :(x*erf(b*x)^2+-1*2^(1/2)*pi^-1//2*b^-1*erf(b*x*2^(1/2))+2*pi^-1//2*b^-1*erf(b*x)*exp(-1*b^2*x^2))
@test integrate((c+d*x)^2*erf(a+b*x)^2, x) == :(b^-3*(b*c+-1*a*d)^2*erf(a+b*x)^2*(a+b*x)+-1//2*d*b^-3*erf(a+b*x)^2*(b*c+-1*a*d)+1//3*b^-3*d^2*(a+b*x)^3*erf(a+b*x)^2+d*pi^-1*b^-3*(b*c+-1*a*d)*exp(-2*(a+b*x)^2)+d*b^-3*(a+b*x)^2*erf(a+b*x)^2*(b*c+-1*a*d)+-1*2^(1/2)*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erf(2^(1/2)*(a+b*x))+2*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erf(a+b*x)*exp(-1*(a+b*x)^2)+-5//12*2^(1/2)*pi^-1//2*b^-3*d^2*erf(2^(1/2)*(a+b*x))+1//3*pi^-1*b^-3*d^2*(a+b*x)*exp(-2*(a+b*x)^2)+2//3*pi^-1//2*b^-3*d^2*erf(a+b*x)*exp(-1*(a+b*x)^2)+2//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*erf(a+b*x)*exp(-1*(a+b*x)^2)+2*d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*erf(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erf(a+b*x)^2*(c+d*x), x) == :(-1//4*d*b^-2*erf(a+b*x)^2+b^-2*erf(a+b*x)^2*(a+b*x)*(b*c+-1*a*d)+(1/2)*d*pi^-1*b^-2*exp(-2*(a+b*x)^2)+(1/2)*d*b^-2*(a+b*x)^2*erf(a+b*x)^2+pi^-1//2*b^-2*(-2*a*d+2*b*c)*erf(a+b*x)*exp(-1*(a+b*x)^2)+-1*2^(1/2)*pi^-1//2*b^-2*(b*c+-1*a*d)*erf(2^(1/2)*(a+b*x))+d*pi^-1//2*b^-2*(a+b*x)*erf(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erf(a+b*x)^2, x) == :(b^-1*erf(a+b*x)^2*(a+b*x)+-1*2^(1/2)*pi^-1//2*b^-1*erf(2^(1/2)*(a+b*x))+2*pi^-1//2*b^-1*erf(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(x^2*erf(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*erf(d*(a+b*log(c*x^n)))+-1//3*x^3*(c*x^n)^(-3*n^-1)*erf((1/2)*b^-1*d^-1*(-3*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(9+-12*a*b*n*d^2)))
@test integrate(x*erf(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*erf(d*(a+b*log(c*x^n)))+-1//2*x^2*(c*x^n)^(-2*n^-1)*erf(b^-1*d^-1*(-1*n^-1+a*b*d^2+b^2*d^2*log(c*x^n)))*exp(b^-2*d^-2*n^-2*(1+-2*a*b*n*d^2)))
@test integrate(erf(d*(a+b*log(c*x^n))), x) == :(x*erf(d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1)*erf((1/2)*b^-1*d^-1*(-1*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(1+-4*a*b*n*d^2)))
@test integrate(x^-1*erf(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*erf(d*(a+b*log(c*x^n)))+pi^-1//2*b^-1*d^-1*n^-1*exp(-1*d^2*(a+b*log(c*x^n))^2))
@test integrate(x^-2*erf(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*erf(d*(a+b*log(c*x^n)))+x^-1*(c*x^n)^(n^-1)*erf((1/2)*b^-1*d^-1*(n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(a*b^-1*n^-1+1//4*b^-2*d^-2*n^-2))
@test integrate(x^-3*erf(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*erf(d*(a+b*log(c*x^n)))+(1/2)*x^-2*(c*x^n)^(2*n^-1)*erf(b^-1*d^-1*n^-1*(1+a*b*n*d^2+n*b^2*d^2*log(c*x^n)))*exp(b^-2*d^-2*n^-2*(1+2*a*b*n*d^2)))
@test integrate((e*x)^m*erf(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*erf(d*(a+b*log(c*x^n)))+x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*erf((1/2)*b^-1*d^-1*n^-1*(1+m+-2*a*b*n*d^2+-2*n*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(1+m)*(1+m+-4*a*b*n*d^2)))
@test integrate(erf(b*x)^2*exp(c+-1*b^2*x^2), x) == :(1//6*pi^(1/2)*b^-1*erf(b*x)^3*exp(c))
@test integrate(erf(b*x)*exp(c+-1*b^2*x^2), x) == :(1//4*pi^(1/2)*b^-1*erf(b*x)^2*exp(c))
@test integrate(erf(b*x)^-1*exp(c+-1*b^2*x^2), x) == :((1/2)*pi^(1/2)*b^-1*exp(c)*log(erf(b*x)))
@test integrate(erf(b*x)^-2*exp(c+-1*b^2*x^2), x) == :(-1//2*pi^(1/2)*b^-1*erf(b*x)^-1*exp(c))
@test integrate(erf(b*x)^-3*exp(c+-1*b^2*x^2), x) == :(-1//4*pi^(1/2)*b^-1*erf(b*x)^-2*exp(c))
@test integrate(erf(b*x)^n*exp(c+-1*b^2*x^2), x) == :((1/2)*pi^(1/2)*b^-1*(1+n)^-1*erf(b*x)^(1+n)*exp(c))
@test integrate(x^5*erf(b*x)*exp(c+d*x^2), x) == :(d^-3*erf(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^4*erf(b*x)*exp(c+d*x^2)+-1*d^-2*x^2*erf(b*x)*exp(c+d*x^2)+(1/2)*b*d^-2*(b^2+-1d)^-3//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+-1*b*d^-3*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+-3//8*b*d^-1*(b^2+-1d)^-5//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+b*pi^-1//2*d^-1*x^3*(-2d+2*b^2)^-1*exp(c+-1*x^2*(b^2+-1d))+-1*b*x*pi^-1//2*d^-2*(b^2+-1d)^-1*exp(c+-1*x^2*(b^2+-1d))+3//4*b*x*pi^-1//2*d^-1*(b^2+-1d)^-2*exp(c+-1*x^2*(b^2+-1d)))
@test integrate(x^3*erf(b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erf(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erf(b*x)*exp(c+d*x^2)+(1/2)*b*d^-2*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+-1//4*b*d^-1*(b^2+-1d)^-3//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+b*x*pi^-1//2*d^-1*(-2d+2*b^2)^-1*exp(c+-1*x^2*(b^2+-1d)))
@test integrate(x*erf(b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erf(b*x)*exp(c+d*x^2)+-1//2*b*d^-1*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c))
@test integrate(x^5*erf(b*x)*exp(c+b^2*x^2), x) == :(b^-6*erf(b*x)*exp(c+b^2*x^2)+(1/2)*b^-2*x^4*erf(b*x)*exp(c+b^2*x^2)+-1*b^-4*x^2*erf(b*x)*exp(c+b^2*x^2)+-2*x*pi^-1//2*b^-5*exp(c)+-1//5*pi^-1//2*b^-1*x^5*exp(c)+2//3*pi^-1//2*b^-3*x^3*exp(c))
@test integrate(x^3*erf(b*x)*exp(c+b^2*x^2), x) == :(-1//2*b^-4*erf(b*x)*exp(c+b^2*x^2)+x*pi^-1//2*b^-3*exp(c)+(1/2)*b^-2*x^2*erf(b*x)*exp(c+b^2*x^2)+-1//3*pi^-1//2*b^-1*x^3*exp(c))
@test integrate(x*erf(b*x)*exp(c+b^2*x^2), x) == :((1/2)*b^-2*erf(b*x)*exp(c+b^2*x^2)+-1*x*pi^-1//2*b^-1*exp(c))
@test integrate(x^5*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1*b^-6*erf(b*x)*exp(-1*b^2*x^2)+43//64*2^(1/2)*b^-6*erf(b*x*2^(1/2))+-1*b^-4*x^2*erf(b*x)*exp(-1*b^2*x^2)+-11//16*x*pi^-1//2*b^-5*exp(-2*b^2*x^2)+-1//2*b^-2*x^4*erf(b*x)*exp(-1*b^2*x^2)+-1//4*pi^-1//2*b^-3*x^3*exp(-2*b^2*x^2))
@test integrate(x^3*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1//2*b^-4*erf(b*x)*exp(-1*b^2*x^2)+5//16*2^(1/2)*b^-4*erf(b*x*2^(1/2))+-1//2*b^-2*x^2*erf(b*x)*exp(-1*b^2*x^2)+-1//4*x*pi^-1//2*b^-3*exp(-2*b^2*x^2))
@test integrate(x*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1//2*b^-2*erf(b*x)*exp(-1*b^2*x^2)+1//4*2^(1/2)*b^-2*erf(b*x*2^(1/2)))
@test integrate(x^4*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1//2*pi^-1//2*b^-5*exp(-2*b^2*x^2)+3//16*pi^(1/2)*b^-5*erf(b*x)^2+-3//4*x*b^-4*erf(b*x)*exp(-1*b^2*x^2)+-1//2*b^-2*x^3*erf(b*x)*exp(-1*b^2*x^2)+-1//4*pi^-1//2*b^-3*x^2*exp(-2*b^2*x^2))
@test integrate(x^2*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1//4*pi^-1//2*b^-3*exp(-2*b^2*x^2)+1//8*pi^(1/2)*b^-3*erf(b*x)^2+-1//2*x*b^-2*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(erf(b*x)*exp(-1*b^2*x^2), x) == :(1//4*pi^(1/2)*b^-1*erf(b*x)^2)
@test integrate(x^-2*erf(b*x)*exp(-1*b^2*x^2), x) == :(b*pi^-1//2*expintegralei(-2*b^2*x^2)+-1*x^-1*erf(b*x)*exp(-1*b^2*x^2)+-1//2*b*pi^(1/2)*erf(b*x)^2)
@test integrate(x^-4*erf(b*x)*exp(-1*b^2*x^2), x) == :(-4//3*pi^-1//2*b^3*expintegralei(-2*b^2*x^2)+-1//3*x^-3*erf(b*x)*exp(-1*b^2*x^2)+1//3*pi^(1/2)*b^3*erf(b*x)^2+-1//3*b*pi^-1//2*x^-2*exp(-2*b^2*x^2)+2//3*b^2*x^-1*erf(b*x)*exp(-1*b^2*x^2))
@test integrate(x^3*erf(a+b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erf(a+b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erf(a+b*x)*exp(c+d*x^2)+(1/2)*b*d^-2*(b^2+-1d)^-1//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1)+-1//4*b*d^-1*(b^2+-1d)^-3//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1)+b*x*pi^-1//2*d^-1*(-2d+2*b^2)^-1*exp(c+-1*a^2+-1*x^2*(b^2+-1d)+-2*a*b*x)+-1//2*a*pi^-1//2*b^2*d^-1*(b^2+-1d)^-2*exp(c+-1*a^2+-1*x^2*(b^2+-1d)+-2*a*b*x)+-1//2*a^2*b^3*d^-1*(b^2+-1d)^-5//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1))
@test integrate(x*erf(a+b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erf(a+b*x)*exp(c+d*x^2)+-1//2*b*d^-1*(b^2+-1d)^-1//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1))
@test integrate(x^-3*erf(b*x)*exp(-1*b^2*x^2)+b^2*x^-1*erf(b*x)*exp(-1*b^2*x^2), x) == :(-1*2^(1/2)*b^2*erf(b*x*2^(1/2))+-1//2*x^-2*erf(b*x)*exp(-1*b^2*x^2)+-1*b*pi^-1//2*x^-1*exp(-2*b^2*x^2))
@test integrate(x^5*erfc(b*x), x) == :(1//6*x^6*erfc(b*x)+5//16*b^-6*erf(b*x)+-5//8*x*pi^-1//2*b^-5*exp(-1*b^2*x^2)+-5//12*pi^-1//2*b^-3*x^3*exp(-1*b^2*x^2)+-1//6*pi^-1//2*b^-1*x^5*exp(-1*b^2*x^2))
@test integrate(x^3*erfc(b*x), x) == :(1//4*x^4*erfc(b*x)+3//16*b^-4*erf(b*x)+-3//8*x*pi^-1//2*b^-3*exp(-1*b^2*x^2)+-1//4*pi^-1//2*b^-1*x^3*exp(-1*b^2*x^2))
@test integrate(x*erfc(b*x), x) == :((1/2)*x^2*erfc(b*x)+1//4*b^-2*erf(b*x)+-1//2*x*pi^-1//2*b^-1*exp(-1*b^2*x^2))
@test integrate(x^-3*erfc(b*x), x) == :(b^2*erf(b*x)+-1//2*x^-2*erfc(b*x)+b*pi^-1//2*x^-1*exp(-1*b^2*x^2))
@test integrate(x^-5*erfc(b*x), x) == :(-1//3*b^4*erf(b*x)+-1//4*x^-4*erfc(b*x)+-1//3*pi^-1//2*b^3*x^-1*exp(-1*b^2*x^2)+1//6*b*pi^-1//2*x^-3*exp(-1*b^2*x^2))
@test integrate(x^-7*erfc(b*x), x) == :(-1//6*x^-6*erfc(b*x)+4//45*b^6*erf(b*x)+-2//45*pi^-1//2*b^3*x^-3*exp(-1*b^2*x^2)+1//15*b*pi^-1//2*x^-5*exp(-1*b^2*x^2)+4//45*pi^-1//2*b^5*x^-1*exp(-1*b^2*x^2))
@test integrate(x^6*erfc(b*x), x) == :(1//7*x^7*erfc(b*x)+-6//7*pi^-1//2*b^-7*exp(-1*b^2*x^2)+-6//7*pi^-1//2*b^-5*x^2*exp(-1*b^2*x^2)+-3//7*pi^-1//2*b^-3*x^4*exp(-1*b^2*x^2)+-1//7*pi^-1//2*b^-1*x^6*exp(-1*b^2*x^2))
@test integrate(x^4*erfc(b*x), x) == :(1//5*x^5*erfc(b*x)+-2//5*pi^-1//2*b^-5*exp(-1*b^2*x^2)+-2//5*pi^-1//2*b^-3*x^2*exp(-1*b^2*x^2)+-1//5*pi^-1//2*b^-1*x^4*exp(-1*b^2*x^2))
@test integrate(x^2*erfc(b*x), x) == :(1//3*x^3*erfc(b*x)+-1//3*pi^-1//2*b^-3*exp(-1*b^2*x^2)+-1//3*pi^-1//2*b^-1*x^2*exp(-1*b^2*x^2))
@test integrate(erfc(b*x), x) == :(x*erfc(b*x)+-1*pi^-1//2*b^-1*exp(-1*b^2*x^2))
@test integrate(x^-2*erfc(b*x), x) == :(-1*x^-1*erfc(b*x)+-1*b*pi^-1//2*expintegralei(-1*b^2*x^2))
@test integrate(x^-4*erfc(b*x), x) == :(-1//3*x^-3*erfc(b*x)+1//3*pi^-1//2*b^3*expintegralei(-1*b^2*x^2)+1//3*b*pi^-1//2*x^-2*exp(-1*b^2*x^2))
@test integrate(x^-6*erfc(b*x), x) == :(-1//5*x^-5*erfc(b*x)+-1//10*pi^-1//2*b^5*expintegralei(-1*b^2*x^2)+-1//10*pi^-1//2*b^3*x^-2*exp(-1*b^2*x^2)+1//10*b*pi^-1//2*x^-4*exp(-1*b^2*x^2))
@test integrate((c+d*x)^3*erfc(a+b*x), x) == :(1//4*d^-1*(c+d*x)^4*erfc(a+b*x)+3//16*b^-4*d^3*erf(a+b*x)+-1*pi^-1//2*b^-4*(b*c+-1*a*d)^3*exp(-1*(a+b*x)^2)+1//4*b^-4*d^-1*(b*c+-1*a*d)^4*erf(a+b*x)+3//4*d*b^-4*(b*c+-1*a*d)^2*erf(a+b*x)+-1*pi^-1//2*b^-4*d^2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+-3//8*pi^-1//2*b^-4*d^3*(a+b*x)*exp(-1*(a+b*x)^2)+-1//4*pi^-1//2*b^-4*d^3*(a+b*x)^3*exp(-1*(a+b*x)^2)+-1*pi^-1//2*b^-4*d^2*(a+b*x)^2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+-3//2*d*pi^-1//2*b^-4*(b*c+-1*a*d)^2*(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate((c+d*x)^2*erfc(a+b*x), x) == :(1//3*d^-1*(c+d*x)^3*erfc(a+b*x)+(1/2)*d*b^-3*(b*c+-1*a*d)*erf(a+b*x)+-1*pi^-1//2*b^-3*(b*c+-1*a*d)^2*exp(-1*(a+b*x)^2)+-1//3*pi^-1//2*b^-3*d^2*exp(-1*(a+b*x)^2)+1//3*b^-3*d^-1*(b*c+-1*a*d)^3*erf(a+b*x)+-1//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*exp(-1*(a+b*x)^2)+-1*d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*exp(-1*(a+b*x)^2))
@test integrate((c+d*x)*erfc(a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*erfc(a+b*x)+1//4*d*b^-2*erf(a+b*x)+(1/2)*b^-2*d^-1*(b*c+-1*a*d)^2*erf(a+b*x)+-1*pi^-1//2*b^-2*(b*c+-1*a*d)*exp(-1*(a+b*x)^2)+-1//2*d*pi^-1//2*b^-2*(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erfc(a+b*x), x) == :(b^-1*(a+b*x)*erfc(a+b*x)+-1*pi^-1//2*b^-1*exp(-1*(a+b*x)^2))
@test integrate(x^5*erfc(b*x)^2, x) == :(-5//16*b^-6*erfc(b*x)^2+1//6*x^6*erfc(b*x)^2+11//12*pi^-1*b^-6*exp(-2*b^2*x^2)+1//6*pi^-1*b^-2*x^4*exp(-2*b^2*x^2)+7//12*pi^-1*b^-4*x^2*exp(-2*b^2*x^2)+-5//4*x*pi^-1//2*b^-5*erfc(b*x)*exp(-1*b^2*x^2)+-5//6*pi^-1//2*b^-3*x^3*erfc(b*x)*exp(-1*b^2*x^2)+-1//3*pi^-1//2*b^-1*x^5*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^3*erfc(b*x)^2, x) == :(-3//16*b^-4*erfc(b*x)^2+1//4*x^4*erfc(b*x)^2+(1/2)*pi^-1*b^-4*exp(-2*b^2*x^2)+1//4*pi^-1*b^-2*x^2*exp(-2*b^2*x^2)+-3//4*x*pi^-1//2*b^-3*erfc(b*x)*exp(-1*b^2*x^2)+-1//2*pi^-1//2*b^-1*x^3*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x*erfc(b*x)^2, x) == :((1/2)*x^2*erfc(b*x)^2+-1//4*b^-2*erfc(b*x)^2+(1/2)*pi^-1*b^-2*exp(-2*b^2*x^2)+-1*x*pi^-1//2*b^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-3*erfc(b*x)^2, x) == :(-1*b^2*erfc(b*x)^2+-1//2*x^-2*erfc(b*x)^2+2*pi^-1*b^2*expintegralei(-2*b^2*x^2)+2*b*pi^-1//2*x^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-5*erfc(b*x)^2, x) == :(-1//4*x^-4*erfc(b*x)^2+1//3*b^4*erfc(b*x)^2+-4//3*pi^-1*b^4*expintegralei(-2*b^2*x^2)+-1//3*pi^-1*b^2*x^-2*exp(-2*b^2*x^2)+-2//3*pi^-1//2*b^3*x^-1*erfc(b*x)*exp(-1*b^2*x^2)+1//3*b*pi^-1//2*x^-3*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-7*erfc(b*x)^2, x) == :(-4//45*b^6*erfc(b*x)^2+-1//6*x^-6*erfc(b*x)^2+28//45*pi^-1*b^6*expintegralei(-2*b^2*x^2)+-1//15*pi^-1*b^2*x^-4*exp(-2*b^2*x^2)+2//9*pi^-1*b^4*x^-2*exp(-2*b^2*x^2)+-4//45*pi^-1//2*b^3*x^-3*erfc(b*x)*exp(-1*b^2*x^2)+2//15*b*pi^-1//2*x^-5*erfc(b*x)*exp(-1*b^2*x^2)+8//45*pi^-1//2*b^5*x^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^4*erfc(b*x)^2, x) == :(1//5*x^5*erfc(b*x)^2+-43//80*2^(1/2)*pi^-1//2*b^-5*erf(b*x*2^(1/2))+-4//5*pi^-1//2*b^-5*erfc(b*x)*exp(-1*b^2*x^2)+1//5*pi^-1*b^-2*x^3*exp(-2*b^2*x^2)+11//20*x*pi^-1*b^-4*exp(-2*b^2*x^2)+-4//5*pi^-1//2*b^-3*x^2*erfc(b*x)*exp(-1*b^2*x^2)+-2//5*pi^-1//2*b^-1*x^4*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^2*erfc(b*x)^2, x) == :(1//3*x^3*erfc(b*x)^2+-5//12*2^(1/2)*pi^-1//2*b^-3*erf(b*x*2^(1/2))+-2//3*pi^-1//2*b^-3*erfc(b*x)*exp(-1*b^2*x^2)+1//3*x*pi^-1*b^-2*exp(-2*b^2*x^2)+-2//3*pi^-1//2*b^-1*x^2*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(erfc(b*x)^2, x) == :(x*erfc(b*x)^2+-1*2^(1/2)*pi^-1//2*b^-1*erf(b*x*2^(1/2))+-2*pi^-1//2*b^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate((c+d*x)^2*erfc(a+b*x)^2, x) == :(b^-3*(b*c+-1*a*d)^2*erfc(a+b*x)^2*(a+b*x)+-1//2*d*b^-3*erfc(a+b*x)^2*(b*c+-1*a*d)+1//3*b^-3*d^2*(a+b*x)^3*erfc(a+b*x)^2+d*pi^-1*b^-3*(b*c+-1*a*d)*exp(-2*(a+b*x)^2)+d*b^-3*(a+b*x)^2*erfc(a+b*x)^2*(b*c+-1*a*d)+-1*2^(1/2)*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erf(2^(1/2)*(a+b*x))+-2*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erfc(a+b*x)*exp(-1*(a+b*x)^2)+-5//12*2^(1/2)*pi^-1//2*b^-3*d^2*erf(2^(1/2)*(a+b*x))+-2//3*pi^-1//2*b^-3*d^2*erfc(a+b*x)*exp(-1*(a+b*x)^2)+1//3*pi^-1*b^-3*d^2*(a+b*x)*exp(-2*(a+b*x)^2)+-2//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*erfc(a+b*x)*exp(-1*(a+b*x)^2)+-2*d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*erfc(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erfc(a+b*x)^2*(c+d*x), x) == :(-1//4*d*b^-2*erfc(a+b*x)^2+b^-2*erfc(a+b*x)^2*(a+b*x)*(b*c+-1*a*d)+(1/2)*d*pi^-1*b^-2*exp(-2*(a+b*x)^2)+(1/2)*d*b^-2*(a+b*x)^2*erfc(a+b*x)^2+-1*2^(1/2)*pi^-1//2*b^-2*(b*c+-1*a*d)*erf(2^(1/2)*(a+b*x))+-1*pi^-1//2*b^-2*(-2*a*d+2*b*c)*erfc(a+b*x)*exp(-1*(a+b*x)^2)+-1*d*pi^-1//2*b^-2*(a+b*x)*erfc(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(erfc(a+b*x)^2, x) == :(b^-1*erfc(a+b*x)^2*(a+b*x)+-1*2^(1/2)*pi^-1//2*b^-1*erf(2^(1/2)*(a+b*x))+-2*pi^-1//2*b^-1*erfc(a+b*x)*exp(-1*(a+b*x)^2))
@test integrate(x^2*erfc(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*erfc(d*(a+b*log(c*x^n)))+1//3*x^3*(c*x^n)^(-3*n^-1)*erf((1/2)*b^-1*d^-1*(-3*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(9+-12*a*b*n*d^2)))
@test integrate(x*erfc(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*erfc(d*(a+b*log(c*x^n)))+(1/2)*x^2*(c*x^n)^(-2*n^-1)*erf(b^-1*d^-1*(-1*n^-1+a*b*d^2+b^2*d^2*log(c*x^n)))*exp(b^-2*d^-2*n^-2*(1+-2*a*b*n*d^2)))
@test integrate(erfc(d*(a+b*log(c*x^n))), x) == :(x*erfc(d*(a+b*log(c*x^n)))+x*(c*x^n)^(-1*n^-1)*erf((1/2)*b^-1*d^-1*(-1*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(1+-4*a*b*n*d^2)))
@test integrate(x^-1*erfc(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*erfc(d*(a+b*log(c*x^n)))+-1*pi^-1//2*b^-1*d^-1*n^-1*exp(-1*d^2*(a+b*log(c*x^n))^2))
@test integrate(x^-2*erfc(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*erfc(d*(a+b*log(c*x^n)))+-1*x^-1*(c*x^n)^(n^-1)*erf((1/2)*b^-1*d^-1*(n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(a*b^-1*n^-1+1//4*b^-2*d^-2*n^-2))
@test integrate(x^-3*erfc(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*erfc(d*(a+b*log(c*x^n)))+-1//2*x^-2*(c*x^n)^(2*n^-1)*erf(b^-1*d^-1*n^-1*(1+a*b*n*d^2+n*b^2*d^2*log(c*x^n)))*exp(b^-2*d^-2*n^-2*(1+2*a*b*n*d^2)))
@test integrate((e*x)^m*erfc(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*erfc(d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*erf((1/2)*b^-1*d^-1*n^-1*(1+m+-2*a*b*n*d^2+-2*n*b^2*d^2*log(c*x^n)))*exp(1//4*b^-2*d^-2*n^-2*(1+m)*(1+m+-4*a*b*n*d^2)))
@test integrate(erfc(b*x)^2*exp(c+-1*b^2*x^2), x) == :(-1//6*pi^(1/2)*b^-1*erfc(b*x)^3*exp(c))
@test integrate(erfc(b*x)*exp(c+-1*b^2*x^2), x) == :(-1//4*pi^(1/2)*b^-1*erfc(b*x)^2*exp(c))
@test integrate(erfc(b*x)^-1*exp(c+-1*b^2*x^2), x) == :(-1//2*pi^(1/2)*b^-1*exp(c)*log(erfc(b*x)))
@test integrate(erfc(b*x)^-2*exp(c+-1*b^2*x^2), x) == :((1/2)*pi^(1/2)*b^-1*erfc(b*x)^-1*exp(c))
@test integrate(erfc(b*x)^-3*exp(c+-1*b^2*x^2), x) == :(1//4*pi^(1/2)*b^-1*erfc(b*x)^-2*exp(c))
@test integrate(erfc(b*x)^n*exp(c+-1*b^2*x^2), x) == :(-1//2*pi^(1/2)*b^-1*(1+n)^-1*erfc(b*x)^(1+n)*exp(c))
@test integrate(x^5*erfc(b*x)*exp(c+d*x^2), x) == :(d^-3*erfc(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^4*erfc(b*x)*exp(c+d*x^2)+-1*d^-2*x^2*erfc(b*x)*exp(c+d*x^2)+b*d^-3*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+-1//2*b*d^-2*(b^2+-1d)^-3//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+3//8*b*d^-1*(b^2+-1d)^-5//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+b*x*pi^-1//2*d^-2*(b^2+-1d)^-1*exp(c+-1*x^2*(b^2+-1d))+-1*b*pi^-1//2*d^-1*x^3*(-2d+2*b^2)^-1*exp(c+-1*x^2*(b^2+-1d))+-3//4*b*x*pi^-1//2*d^-1*(b^2+-1d)^-2*exp(c+-1*x^2*(b^2+-1d)))
@test integrate(x^3*erfc(b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erfc(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erfc(b*x)*exp(c+d*x^2)+-1//2*b*d^-2*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+1//4*b*d^-1*(b^2+-1d)^-3//2*erf(x*(b^2+-1d)^(1/2))*exp(c)+-1*b*x*pi^-1//2*d^-1*(-2d+2*b^2)^-1*exp(c+-1*x^2*(b^2+-1d)))
@test integrate(x*erfc(b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erfc(b*x)*exp(c+d*x^2)+(1/2)*b*d^-1*(b^2+-1d)^-1//2*erf(x*(b^2+-1d)^(1/2))*exp(c))
@test integrate(x^5*erfc(b*x)*exp(c+b^2*x^2), x) == :(b^-6*erfc(b*x)*exp(c+b^2*x^2)+(1/2)*b^-2*x^4*erfc(b*x)*exp(c+b^2*x^2)+-1*b^-4*x^2*erfc(b*x)*exp(c+b^2*x^2)+2*x*pi^-1//2*b^-5*exp(c)+-2//3*pi^-1//2*b^-3*x^3*exp(c)+1//5*pi^-1//2*b^-1*x^5*exp(c))
@test integrate(x^3*erfc(b*x)*exp(c+b^2*x^2), x) == :(-1//2*b^-4*erfc(b*x)*exp(c+b^2*x^2)+(1/2)*b^-2*x^2*erfc(b*x)*exp(c+b^2*x^2)+-1*x*pi^-1//2*b^-3*exp(c)+1//3*pi^-1//2*b^-1*x^3*exp(c))
@test integrate(x*erfc(b*x)*exp(c+b^2*x^2), x) == :((1/2)*b^-2*erfc(b*x)*exp(c+b^2*x^2)+x*pi^-1//2*b^-1*exp(c))
@test integrate(x^5*erfc(b*x)*exp(-1*b^2*x^2), x) == :(-1*b^-6*erfc(b*x)*exp(-1*b^2*x^2)+-43//64*2^(1/2)*b^-6*erf(b*x*2^(1/2))+-1*b^-4*x^2*erfc(b*x)*exp(-1*b^2*x^2)+-1//2*b^-2*x^4*erfc(b*x)*exp(-1*b^2*x^2)+1//4*pi^-1//2*b^-3*x^3*exp(-2*b^2*x^2)+11//16*x*pi^-1//2*b^-5*exp(-2*b^2*x^2))
@test integrate(x^3*erfc(b*x)*exp(-1*b^2*x^2), x) == :(-5//16*2^(1/2)*b^-4*erf(b*x*2^(1/2))+-1//2*b^-4*erfc(b*x)*exp(-1*b^2*x^2)+-1//2*b^-2*x^2*erfc(b*x)*exp(-1*b^2*x^2)+1//4*x*pi^-1//2*b^-3*exp(-2*b^2*x^2))
@test integrate(x*erfc(b*x)*exp(-1*b^2*x^2), x) == :(-1//2*b^-2*erfc(b*x)*exp(-1*b^2*x^2)+-1//4*2^(1/2)*b^-2*erf(b*x*2^(1/2)))
@test integrate(x^4*erfc(b*x)*exp(-1*b^2*x^2), x) == :((1/2)*pi^-1//2*b^-5*exp(-2*b^2*x^2)+-3//16*pi^(1/2)*b^-5*erfc(b*x)^2+-3//4*x*b^-4*erfc(b*x)*exp(-1*b^2*x^2)+-1//2*b^-2*x^3*erfc(b*x)*exp(-1*b^2*x^2)+1//4*pi^-1//2*b^-3*x^2*exp(-2*b^2*x^2))
@test integrate(x^2*erfc(b*x)*exp(-1*b^2*x^2), x) == :(-1//8*pi^(1/2)*b^-3*erfc(b*x)^2+1//4*pi^-1//2*b^-3*exp(-2*b^2*x^2)+-1//2*x*b^-2*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(erfc(b*x)*exp(-1*b^2*x^2), x) == :(-1//4*pi^(1/2)*b^-1*erfc(b*x)^2)
@test integrate(x^-2*erfc(b*x)*exp(-1*b^2*x^2), x) == :((1/2)*b*pi^(1/2)*erfc(b*x)^2+-1*b*pi^-1//2*expintegralei(-2*b^2*x^2)+-1*x^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^-4*erfc(b*x)*exp(-1*b^2*x^2), x) == :(-1//3*pi^(1/2)*b^3*erfc(b*x)^2+-1//3*x^-3*erfc(b*x)*exp(-1*b^2*x^2)+4//3*pi^-1//2*b^3*expintegralei(-2*b^2*x^2)+1//3*b*pi^-1//2*x^-2*exp(-2*b^2*x^2)+2//3*b^2*x^-1*erfc(b*x)*exp(-1*b^2*x^2))
@test integrate(x^3*erfc(a+b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erfc(a+b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erfc(a+b*x)*exp(c+d*x^2)+-1//2*b*d^-2*(b^2+-1d)^-1//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1)+1//4*b*d^-1*(b^2+-1d)^-3//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1)+(1/2)*a*pi^-1//2*b^2*d^-1*(b^2+-1d)^-2*exp(c+-1*a^2+-1*x^2*(b^2+-1d)+-2*a*b*x)+(1/2)*a^2*b^3*d^-1*(b^2+-1d)^-5//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1)+-1*b*x*pi^-1//2*d^-1*(-2d+2*b^2)^-1*exp(c+-1*a^2+-1*x^2*(b^2+-1d)+-2*a*b*x))
@test integrate(x*erfc(a+b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erfc(a+b*x)*exp(c+d*x^2)+(1/2)*b*d^-1*(b^2+-1d)^-1//2*erf((b^2+-1d)^-1//2*(a*b+x*(b^2+-1d)))*exp(c+d*a^2*(b^2+-1d)^-1))
@test integrate(x^-3*erfc(b*x)*exp(-1*b^2*x^2)+b^2*x^-1*erfc(b*x)*exp(-1*b^2*x^2), x) == :(2^(1/2)*b^2*erf(b*x*2^(1/2))+-1//2*x^-2*erfc(b*x)*exp(-1*b^2*x^2)+b*pi^-1//2*x^-1*exp(-2*b^2*x^2))
@test integrate(x^5*erfi(b*x), x) == :(1//6*x^6*erfi(b*x)+5//16*b^-6*erfi(b*x)+-5//8*x*pi^-1//2*b^-5*exp(b^2*x^2)+-1//6*pi^-1//2*b^-1*x^5*exp(b^2*x^2)+5//12*pi^-1//2*b^-3*x^3*exp(b^2*x^2))
@test integrate(x^3*erfi(b*x), x) == :(-3//16*b^-4*erfi(b*x)+1//4*x^4*erfi(b*x)+-1//4*pi^-1//2*b^-1*x^3*exp(b^2*x^2)+3//8*x*pi^-1//2*b^-3*exp(b^2*x^2))
@test integrate(x*erfi(b*x), x) == :((1/2)*x^2*erfi(b*x)+1//4*b^-2*erfi(b*x)+-1//2*x*pi^-1//2*b^-1*exp(b^2*x^2))
@test integrate(x^-3*erfi(b*x), x) == :(b^2*erfi(b*x)+-1//2*x^-2*erfi(b*x)+-1*b*pi^-1//2*x^-1*exp(b^2*x^2))
@test integrate(x^-5*erfi(b*x), x) == :(-1//4*x^-4*erfi(b*x)+1//3*b^4*erfi(b*x)+-1//3*pi^-1//2*b^3*x^-1*exp(b^2*x^2)+-1//6*b*pi^-1//2*x^-3*exp(b^2*x^2))
@test integrate(x^-7*erfi(b*x), x) == :(-1//6*x^-6*erfi(b*x)+4//45*b^6*erfi(b*x)+-4//45*pi^-1//2*b^5*x^-1*exp(b^2*x^2)+-2//45*pi^-1//2*b^3*x^-3*exp(b^2*x^2)+-1//15*b*pi^-1//2*x^-5*exp(b^2*x^2))
@test integrate(x^6*erfi(b*x), x) == :(1//7*x^7*erfi(b*x)+6//7*pi^-1//2*b^-7*exp(b^2*x^2)+-6//7*pi^-1//2*b^-5*x^2*exp(b^2*x^2)+-1//7*pi^-1//2*b^-1*x^6*exp(b^2*x^2)+3//7*pi^-1//2*b^-3*x^4*exp(b^2*x^2))
@test integrate(x^4*erfi(b*x), x) == :(1//5*x^5*erfi(b*x)+-2//5*pi^-1//2*b^-5*exp(b^2*x^2)+-1//5*pi^-1//2*b^-1*x^4*exp(b^2*x^2)+2//5*pi^-1//2*b^-3*x^2*exp(b^2*x^2))
@test integrate(x^2*erfi(b*x), x) == :(1//3*x^3*erfi(b*x)+1//3*pi^-1//2*b^-3*exp(b^2*x^2)+-1//3*pi^-1//2*b^-1*x^2*exp(b^2*x^2))
@test integrate(erfi(b*x), x) == :(x*erfi(b*x)+-1*pi^-1//2*b^-1*exp(b^2*x^2))
@test integrate(x^-2*erfi(b*x), x) == :(-1*x^-1*erfi(b*x)+b*pi^-1//2*expintegralei(b^2*x^2))
@test integrate(x^-4*erfi(b*x), x) == :(-1//3*x^-3*erfi(b*x)+1//3*pi^-1//2*b^3*expintegralei(b^2*x^2)+-1//3*b*pi^-1//2*x^-2*exp(b^2*x^2))
@test integrate(x^-6*erfi(b*x), x) == :(-1//5*x^-5*erfi(b*x)+1//10*pi^-1//2*b^5*expintegralei(b^2*x^2)+-1//10*b*pi^-1//2*x^-4*exp(b^2*x^2)+-1//10*pi^-1//2*b^3*x^-2*exp(b^2*x^2))
@test integrate((c+d*x)^3*erfi(a+b*x), x) == :(-3//16*b^-4*d^3*erfi(a+b*x)+1//4*d^-1*(c+d*x)^4*erfi(a+b*x)+-1*pi^-1//2*b^-4*(b*c+-1*a*d)^3*exp((a+b*x)^2)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*erfi(a+b*x)+3//4*d*b^-4*(b*c+-1*a*d)^2*erfi(a+b*x)+pi^-1//2*b^-4*d^2*(b*c+-1*a*d)*exp((a+b*x)^2)+-1//4*pi^-1//2*b^-4*d^3*(a+b*x)^3*exp((a+b*x)^2)+3//8*pi^-1//2*b^-4*d^3*(a+b*x)*exp((a+b*x)^2)+-1*pi^-1//2*b^-4*d^2*(a+b*x)^2*(b*c+-1*a*d)*exp((a+b*x)^2)+-3//2*d*pi^-1//2*b^-4*(b*c+-1*a*d)^2*(a+b*x)*exp((a+b*x)^2))
@test integrate((c+d*x)^2*erfi(a+b*x), x) == :(1//3*d^-1*(c+d*x)^3*erfi(a+b*x)+(1/2)*d*b^-3*(b*c+-1*a*d)*erfi(a+b*x)+-1*pi^-1//2*b^-3*(b*c+-1*a*d)^2*exp((a+b*x)^2)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*erfi(a+b*x)+1//3*pi^-1//2*b^-3*d^2*exp((a+b*x)^2)+-1//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*exp((a+b*x)^2)+-1*d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*exp((a+b*x)^2))
@test integrate((c+d*x)*erfi(a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*erfi(a+b*x)+1//4*d*b^-2*erfi(a+b*x)+-1*pi^-1//2*b^-2*(b*c+-1*a*d)*exp((a+b*x)^2)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*erfi(a+b*x)+-1//2*d*pi^-1//2*b^-2*(a+b*x)*exp((a+b*x)^2))
@test integrate(erfi(a+b*x), x) == :(b^-1*(a+b*x)*erfi(a+b*x)+-1*pi^-1//2*b^-1*exp((a+b*x)^2))
@test integrate(x^5*erfi(b*x)^2, x) == :(1//6*x^6*erfi(b*x)^2+5//16*b^-6*erfi(b*x)^2+11//12*pi^-1*b^-6*exp(2*b^2*x^2)+-7//12*pi^-1*b^-4*x^2*exp(2*b^2*x^2)+1//6*pi^-1*b^-2*x^4*exp(2*b^2*x^2)+-5//4*x*pi^-1//2*b^-5*erfi(b*x)*exp(b^2*x^2)+-1//3*pi^-1//2*b^-1*x^5*erfi(b*x)*exp(b^2*x^2)+5//6*pi^-1//2*b^-3*x^3*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^3*erfi(b*x)^2, x) == :(-3//16*b^-4*erfi(b*x)^2+1//4*x^4*erfi(b*x)^2+-1//2*pi^-1*b^-4*exp(2*b^2*x^2)+1//4*pi^-1*b^-2*x^2*exp(2*b^2*x^2)+-1//2*pi^-1//2*b^-1*x^3*erfi(b*x)*exp(b^2*x^2)+3//4*x*pi^-1//2*b^-3*erfi(b*x)*exp(b^2*x^2))
@test integrate(x*erfi(b*x)^2, x) == :((1/2)*x^2*erfi(b*x)^2+1//4*b^-2*erfi(b*x)^2+(1/2)*pi^-1*b^-2*exp(2*b^2*x^2)+-1*x*pi^-1//2*b^-1*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^-3*erfi(b*x)^2, x) == :(b^2*erfi(b*x)^2+-1//2*x^-2*erfi(b*x)^2+2*pi^-1*b^2*expintegralei(2*b^2*x^2)+-2*b*pi^-1//2*x^-1*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^-5*erfi(b*x)^2, x) == :(-1//4*x^-4*erfi(b*x)^2+1//3*b^4*erfi(b*x)^2+4//3*pi^-1*b^4*expintegralei(2*b^2*x^2)+-1//3*pi^-1*b^2*x^-2*exp(2*b^2*x^2)+-2//3*pi^-1//2*b^3*x^-1*erfi(b*x)*exp(b^2*x^2)+-1//3*b*pi^-1//2*x^-3*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^-7*erfi(b*x)^2, x) == :(-1//6*x^-6*erfi(b*x)^2+4//45*b^6*erfi(b*x)^2+28//45*pi^-1*b^6*expintegralei(2*b^2*x^2)+-2//9*pi^-1*b^4*x^-2*exp(2*b^2*x^2)+-1//15*pi^-1*b^2*x^-4*exp(2*b^2*x^2)+-8//45*pi^-1//2*b^5*x^-1*erfi(b*x)*exp(b^2*x^2)+-4//45*pi^-1//2*b^3*x^-3*erfi(b*x)*exp(b^2*x^2)+-2//15*b*pi^-1//2*x^-5*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^4*erfi(b*x)^2, x) == :(1//5*x^5*erfi(b*x)^2+-11//20*x*pi^-1*b^-4*exp(2*b^2*x^2)+-4//5*pi^-1//2*b^-5*erfi(b*x)*exp(b^2*x^2)+1//5*pi^-1*b^-2*x^3*exp(2*b^2*x^2)+43//80*2^(1/2)*pi^-1//2*b^-5*erfi(b*x*2^(1/2))+-2//5*pi^-1//2*b^-1*x^4*erfi(b*x)*exp(b^2*x^2)+4//5*pi^-1//2*b^-3*x^2*erfi(b*x)*exp(b^2*x^2))
@test integrate(x^2*erfi(b*x)^2, x) == :(1//3*x^3*erfi(b*x)^2+-5//12*2^(1/2)*pi^-1//2*b^-3*erfi(b*x*2^(1/2))+1//3*x*pi^-1*b^-2*exp(2*b^2*x^2)+2//3*pi^-1//2*b^-3*erfi(b*x)*exp(b^2*x^2)+-2//3*pi^-1//2*b^-1*x^2*erfi(b*x)*exp(b^2*x^2))
@test integrate(erfi(b*x)^2, x) == :(x*erfi(b*x)^2+2^(1/2)*pi^-1//2*b^-1*erfi(b*x*2^(1/2))+-2*pi^-1//2*b^-1*erfi(b*x)*exp(b^2*x^2))
@test integrate((c+d*x)^2*erfi(a+b*x)^2, x) == :(b^-3*(b*c+-1*a*d)^2*erfi(a+b*x)^2*(a+b*x)+(1/2)*d*b^-3*erfi(a+b*x)^2*(b*c+-1*a*d)+1//3*b^-3*d^2*(a+b*x)^3*erfi(a+b*x)^2+d*pi^-1*b^-3*(b*c+-1*a*d)*exp(2*(a+b*x)^2)+d*b^-3*(a+b*x)^2*erfi(a+b*x)^2*(b*c+-1*a*d)+2^(1/2)*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erfi(2^(1/2)*(a+b*x))+-2*pi^-1//2*b^-3*(b*c+-1*a*d)^2*erfi(a+b*x)*exp((a+b*x)^2)+-5//12*2^(1/2)*pi^-1//2*b^-3*d^2*erfi(2^(1/2)*(a+b*x))+1//3*pi^-1*b^-3*d^2*(a+b*x)*exp(2*(a+b*x)^2)+2//3*pi^-1//2*b^-3*d^2*erfi(a+b*x)*exp((a+b*x)^2)+-2//3*pi^-1//2*b^-3*d^2*(a+b*x)^2*erfi(a+b*x)*exp((a+b*x)^2)+-2*d*pi^-1//2*b^-3*(a+b*x)*(b*c+-1*a*d)*erfi(a+b*x)*exp((a+b*x)^2))
@test integrate(erfi(a+b*x)^2*(c+d*x), x) == :(1//4*d*b^-2*erfi(a+b*x)^2+b^-2*erfi(a+b*x)^2*(a+b*x)*(b*c+-1*a*d)+(1/2)*d*pi^-1*b^-2*exp(2*(a+b*x)^2)+(1/2)*d*b^-2*(a+b*x)^2*erfi(a+b*x)^2+2^(1/2)*pi^-1//2*b^-2*(b*c+-1*a*d)*erfi(2^(1/2)*(a+b*x))+-1*pi^-1//2*b^-2*(-2*a*d+2*b*c)*erfi(a+b*x)*exp((a+b*x)^2)+-1*d*pi^-1//2*b^-2*(a+b*x)*erfi(a+b*x)*exp((a+b*x)^2))
@test integrate(erfi(a+b*x)^2, x) == :(b^-1*erfi(a+b*x)^2*(a+b*x)+2^(1/2)*pi^-1//2*b^-1*erfi(2^(1/2)*(a+b*x))+-2*pi^-1//2*b^-1*erfi(a+b*x)*exp((a+b*x)^2))
@test integrate(x^2*erfi(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*erfi(d*(a+b*log(c*x^n)))+-1//3*x^3*(c*x^n)^(-3*n^-1)*erfi((1/2)*b^-1*d^-1*(3*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(-1//4*b^-2*d^-2*n^-2*(9+12*a*b*n*d^2)))
@test integrate(x*erfi(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*erfi(d*(a+b*log(c*x^n)))+-1//2*x^2*(c*x^n)^(-2*n^-1)*erfi(b^-1*d^-1*(n^-1+a*b*d^2+b^2*d^2*log(c*x^n)))*exp(-1*b^-2*d^-2*n^-2*(1+2*a*b*n*d^2)))
@test integrate(erfi(d*(a+b*log(c*x^n))), x) == :(x*erfi(d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1)*erfi((1/2)*b^-1*d^-1*(n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(-1//4*b^-2*d^-2*n^-2*(1+4*a*b*n*d^2)))
@test integrate(x^-1*erfi(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*erfi(d*(a+b*log(c*x^n)))+-1*pi^-1//2*b^-1*d^-1*n^-1*exp((a*d+b*d*log(c*x^n))^2))
@test integrate(x^-2*erfi(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*erfi(d*(a+b*log(c*x^n)))+x^-1*(c*x^n)^(n^-1)*erfi((1/2)*b^-1*d^-1*(-1*n^-1+2*a*b*d^2+2*b^2*d^2*log(c*x^n)))*exp(a*b^-1*n^-1+-1//4*b^-2*d^-2*n^-2))
@test integrate(x^-3*erfi(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*erfi(d*(a+b*log(c*x^n)))+(1/2)*x^-2*(c*x^n)^(2*n^-1)*erfi(b^-1*d^-1*(-1*n^-1+a*b*d^2+b^2*d^2*log(c*x^n)))*exp(-1*b^-2*d^-2*n^-2*(1+-2*a*b*n*d^2)))
@test integrate((e*x)^m*erfi(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*erfi(d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*erfi((1/2)*b^-1*d^-1*n^-1*(1+m+2*a*b*n*d^2+2*n*b^2*d^2*log(c*x^n)))*exp(-1//4*b^-2*d^-2*n^-2*(1+m)*(1+m+4*a*b*n*d^2)))
@test integrate(erfi(b*x)^2*exp(c+b^2*x^2), x) == :(1//6*pi^(1/2)*b^-1*erfi(b*x)^3*exp(c))
@test integrate(erfi(b*x)*exp(c+b^2*x^2), x) == :(1//4*pi^(1/2)*b^-1*erfi(b*x)^2*exp(c))
@test integrate(erfi(b*x)^-1*exp(c+b^2*x^2), x) == :((1/2)*pi^(1/2)*b^-1*exp(c)*log(erfi(b*x)))
@test integrate(erfi(b*x)^-2*exp(c+b^2*x^2), x) == :(-1//2*pi^(1/2)*b^-1*erfi(b*x)^-1*exp(c))
@test integrate(erfi(b*x)^-3*exp(c+b^2*x^2), x) == :(-1//4*pi^(1/2)*b^-1*erfi(b*x)^-2*exp(c))
@test integrate(erfi(b*x)^n*exp(c+b^2*x^2), x) == :((1/2)*pi^(1/2)*b^-1*(1+n)^-1*erfi(b*x)^(1+n)*exp(c))
@test integrate(x^5*erfi(b*x)*exp(c+d*x^2), x) == :(d^-3*erfi(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^4*erfi(b*x)*exp(c+d*x^2)+-1*d^-2*x^2*erfi(b*x)*exp(c+d*x^2)+-1*b*d^-3*(d+b^2)^-1//2*erfi(x*(d+b^2)^(1/2))*exp(c)+-3//8*b*d^-1*(d+b^2)^-5//2*erfi(x*(d+b^2)^(1/2))*exp(c)+-1//2*b*d^-2*(d+b^2)^-3//2*erfi(x*(d+b^2)^(1/2))*exp(c)+b*x*pi^-1//2*d^-2*(d+b^2)^-1*exp(c+x^2*(d+b^2))+-1//2*b*pi^-1//2*d^-1*x^3*(d+b^2)^-1*exp(c+x^2*(d+b^2))+3//4*b*x*pi^-1//2*d^-1*(d+b^2)^-2*exp(c+x^2*(d+b^2)))
@test integrate(x^3*erfi(b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erfi(b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erfi(b*x)*exp(c+d*x^2)+(1/2)*b*d^-2*(d+b^2)^-1//2*erfi(x*(d+b^2)^(1/2))*exp(c)+1//4*b*d^-1*(d+b^2)^-3//2*erfi(x*(d+b^2)^(1/2))*exp(c)+-1//2*b*x*pi^-1//2*d^-1*(d+b^2)^-1*exp(c+x^2*(d+b^2)))
@test integrate(x*erfi(b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erfi(b*x)*exp(c+d*x^2)+-1//2*b*d^-1*(d+b^2)^-1//2*erfi(x*(d+b^2)^(1/2))*exp(c))
@test integrate(x^5*erfi(b*x)*exp(-1*b^2*x^2), x) == :(-1*b^-6*erfi(b*x)*exp(-1*b^2*x^2)+2*x*pi^-1//2*b^-5+1//5*pi^-1//2*b^-1*x^5+2//3*pi^-1//2*b^-3*x^3+-1*b^-4*x^2*erfi(b*x)*exp(-1*b^2*x^2)+-1//2*b^-2*x^4*erfi(b*x)*exp(-1*b^2*x^2))
@test integrate(x^3*erfi(b*x)*exp(-1*b^2*x^2), x) == :(x*pi^-1//2*b^-3+-1//2*b^-4*erfi(b*x)*exp(-1*b^2*x^2)+1//3*pi^-1//2*b^-1*x^3+-1//2*b^-2*x^2*erfi(b*x)*exp(-1*b^2*x^2))
@test integrate(x*erfi(b*x)*exp(-1*b^2*x^2), x) == :(x*pi^-1//2*b^-1+-1//2*b^-2*erfi(b*x)*exp(-1*b^2*x^2))
@test integrate(x^5*erfi(b*x)*exp(c+b^2*x^2), x) == :(b^-6*erfi(b*x)*exp(c+b^2*x^2)+(1/2)*b^-2*x^4*erfi(b*x)*exp(c+b^2*x^2)+-1*b^-4*x^2*erfi(b*x)*exp(c+b^2*x^2)+-43//64*2^(1/2)*b^-6*erfi(b*x*2^(1/2))*exp(c)+-1//4*pi^-1//2*b^-3*x^3*exp(c+2*b^2*x^2)+11//16*x*pi^-1//2*b^-5*exp(c+2*b^2*x^2))
@test integrate(x^3*erfi(b*x)*exp(c+b^2*x^2), x) == :(-1//2*b^-4*erfi(b*x)*exp(c+b^2*x^2)+(1/2)*b^-2*x^2*erfi(b*x)*exp(c+b^2*x^2)+-1//4*x*pi^-1//2*b^-3*exp(c+2*b^2*x^2)+5//16*2^(1/2)*b^-4*erfi(b*x*2^(1/2))*exp(c))
@test integrate(x*erfi(b*x)*exp(c+b^2*x^2), x) == :((1/2)*b^-2*erfi(b*x)*exp(c+b^2*x^2)+-1//4*2^(1/2)*b^-2*erfi(b*x*2^(1/2))*exp(c))
@test integrate(x^4*erfi(b*x)*exp(c+b^2*x^2), x) == :((1/2)*pi^-1//2*b^-5*exp(c+2*b^2*x^2)+(1/2)*b^-2*x^3*erfi(b*x)*exp(c+b^2*x^2)+-3//4*x*b^-4*erfi(b*x)*exp(c+b^2*x^2)+-1//4*pi^-1//2*b^-3*x^2*exp(c+2*b^2*x^2)+3//16*pi^(1/2)*b^-5*erfi(b*x)^2*exp(c))
@test integrate(x^2*erfi(b*x)*exp(c+b^2*x^2), x) == :(-1//4*pi^-1//2*b^-3*exp(c+2*b^2*x^2)+(1/2)*x*b^-2*erfi(b*x)*exp(c+b^2*x^2)+-1//8*pi^(1/2)*b^-3*erfi(b*x)^2*exp(c))
@test integrate(erfi(b*x)*exp(c+b^2*x^2), x) == :(1//4*pi^(1/2)*b^-1*erfi(b*x)^2*exp(c))
@test integrate(x^-2*erfi(b*x)*exp(c+b^2*x^2), x) == :(-1*x^-1*erfi(b*x)*exp(c+b^2*x^2)+b*pi^-1//2*exp(c)*expintegralei(2*b^2*x^2)+(1/2)*b*pi^(1/2)*erfi(b*x)^2*exp(c))
@test integrate(x^-4*erfi(b*x)*exp(c+b^2*x^2), x) == :(-1//3*x^-3*erfi(b*x)*exp(c+b^2*x^2)+-2//3*b^2*x^-1*erfi(b*x)*exp(c+b^2*x^2)+-1//3*b*pi^-1//2*x^-2*exp(c+2*b^2*x^2)+1//3*pi^(1/2)*b^3*erfi(b*x)^2*exp(c)+4//3*pi^-1//2*b^3*exp(c)*expintegralei(2*b^2*x^2))
@test integrate(x^3*erfi(a+b*x)*exp(c+d*x^2), x) == :(-1//2*d^-2*erfi(a+b*x)*exp(c+d*x^2)+(1/2)*d^-1*x^2*erfi(a+b*x)*exp(c+d*x^2)+(1/2)*b*d^-2*(d+b^2)^-1//2*erfi((d+b^2)^-1//2*(a*b+x*(d+b^2)))*exp(c+d*a^2*(d+b^2)^-1)+1//4*b*d^-1*(d+b^2)^-3//2*erfi((d+b^2)^-1//2*(a*b+x*(d+b^2)))*exp(c+d*a^2*(d+b^2)^-1)+(1/2)*a*pi^-1//2*b^2*d^-1*(d+b^2)^-2*exp(c+a^2+x^2*(d+b^2)+2*a*b*x)+-1//2*b*x*pi^-1//2*d^-1*(d+b^2)^-1*exp(c+a^2+x^2*(d+b^2)+2*a*b*x)+-1//2*a^2*b^3*d^-1*(d+b^2)^-5//2*erfi((d+b^2)^-1//2*(a*b+x*(d+b^2)))*exp(c+d*a^2*(d+b^2)^-1))
@test integrate(x*erfi(a+b*x)*exp(c+d*x^2), x) == :((1/2)*d^-1*erfi(a+b*x)*exp(c+d*x^2)+-1//2*b*d^-1*(d+b^2)^-1//2*erfi((d+b^2)^-1//2*(a*b+x*(d+b^2)))*exp(c+d*a^2*(d+b^2)^-1))
@test integrate(x^-3*erfi(b*x)*exp(-1*b^2*x^2)+b^2*x^-1*erfi(b*x)*exp(-1*b^2*x^2), x) == :(-1*b*pi^-1//2*x^-1+-1//2*x^-2*erfi(b*x)*exp(-1*b^2*x^2))

max_stars_repo_path: test/8 Special functions/8.1 Error functions.jl


text:
using DependencyWalker
using Test
using ObjectFile, Pango_jll

@testset "DependencyWalker.jl" begin
    @test Library(Pango_jll.libpango_path) isa Library{<:ObjectHandle}
    @test Library("this does not exist.foo") isa Library{Missing}
end


text:
const TreesMap = BitArray{2}

function readinput(strs)
    TreesMap(transpose(hcat([[c == '#' for c ∈ l] for l ∈ strs]...)))
end

# First Half

function treesencountered₁(tmap)
    #= =========================================================================
    The trees map of the problem can be projected on cilindrical discrete plane
    where cilinder's height and perimeter are respectively mapₕ and mapₗ.
    The movement happens on a line on a discrete plane.
    The equation for a line starting at (0,0) in a cilindrical semi-infinite
    discrete plane is x = m * y mod mapₗ, where x, y are integers,
    slope m is a rational and m * y is rounded to the nearest integer.
    In our case m = mₕ / mᵥ (right slope / down slope) and the path is not
    this line, but only the dots on this line with vertical coordinate
    multiple of mᵥ (jumping).
    ========================================================================= =#

    mₕ = 3
    mᵥ = 1
    m = mₕ // mᵥ
    mapₕ, mapₗ = size(tmap)
    f(y) = tmap[y + 1, Int(m * y) % mapₗ + 1] # Arrays in Julia start from 1
    count(f, 0:mᵥ:mapₕ - 1) # Arrays in Julia start from 1
end

puzzleinput = readlines("inputs/day03.txt")
treesmap = readinput(puzzleinput)

treesencountered₁(treesmap)

# Second Half

function treesencountered₂(tmap)

    mₕ = [1, 3, 5, 7, 1]
    mᵥ = [1, 1, 1, 1, 2]
    M = collect(zip(mₕ, mᵥ))
    mapₕ, mapₗ = size(tmap)
    f(m) = y -> tmap[y + 1, Int(m[1] // m[2] * y) % mapₗ + 1]
    [count(f(m), 0:m[2]:mapₕ - 1) for m ∈ M]
end

puzzleinput = readlines("inputs/day03.txt")
treesmap = readinput(puzzleinput)

prod(treesencountered₂(treesmap))
# Multiply together the number of trees encountered on each of the listed slopes


text:
module Mean
using ValidatedNumerics
using Distributions
using Dierckx
using StatsBase

abstract BaseEstimator

function estimate(values, method :: BaseEstimator)
    error("Unsupported method")
end

immutable byKnownVariance <: BaseEstimator
    confidence_probability :: Real
    variance :: Real
end

function estimate(values, method::byKnownVariance)
    interval_confidence = getIntervalConfidenceProbability(method.confidence_probability)
    average = mean(values)
    count = length(values)
    quantile = normalQuantile(interval_confidence)
    delta = quantile * sqrt(method.variance) / sqrt(count)
    average + @interval(-delta, delta)
end

immutable byStandardError <: BaseEstimator
end

function estimate(values, method::byStandardError)
    average, variance = mean_and_var(values)
    count = length(values)
    fixed_variance = variance*count/(count - 1)
    delta = sqrt(fixed_variance/count)
    average + @interval(-delta, delta)
end

immutable byUnknownVariance <: BaseEstimator
    confidence_probability :: Real
end

function estimate(values, method::byUnknownVariance)
    average = mean(values)
    count = length(values)
    interval_confidence = getIntervalConfidenceProbability(method.confidence_probability)
    variance = mapreduce((x) -> (x - average)^2, +, values) / (count - 1)
    quantile = studentQuatile(interval_confidence, count - 1)
    delta = quantile * sqrt(variance) / sqrt(count)
    average + @interval(-delta, delta)
end

immutable byMeanAbsDeviation <: BaseEstimator
    confidence_probability :: Real
end

function estimate(values, method :: byMeanAbsDeviation)
    average = mean(values)
    count = length(values)
    interval_confidence = getIntervalConfidenceProbability(method.confidence_probability)
    mean_abs_deviation = mapreduce((x) -> abs(x - average), +, values) / count
    quantile = getMeanAbsDeviationQuantile(interval_confidence, count)
    delta = mean_abs_deviation * quantile
    average + @interval(-delta, delta)
end

immutable byInterQuartileRange <: BaseEstimator
    confidence_probability :: Real
end

function estimate(values, method :: byInterQuartileRange)
    interval_confidence = getIntervalConfidenceProbability(method.confidence_probability)
    count = length(values)
    ordered_values = sort(values)
    interquartile_range = getInterQuartileRange(ordered_values, count)
    interquertile_quitile = getInterQuartileQuintile(interval_confidence, count)
    mediam = getMediamOfSorted(ordered_values, count)
    delta = interquartile_range * interquertile_quitile
    mediam + @interval(-delta, delta)
end

function getIntervalConfidenceProbability(confidence_probability)
    (1.0 + confidence_probability) * 0.5
end

function getInterQuartileRange(values, count)
    # TODO this can be replaced with StatsBase.iqr
    quartile_index = div((count + 1), 4)
    values[3*quartile_index] - values[quartile_index]
end

function getMediamOfSorted(values, count)
    if isodd(count)
        0.5*(values[div(count, 2)] + values[div(count, 2) + 1])
    else
        values[div((count + 1), 2)]
    end
end

function getInterQuartileQuintile(value, count)
    const knots = [
    11 0.470 0.623 0.876
    15 0.400 0.514 0.678
    19 0.354 0.448 0.573
    23 0.321 0.402 0.506
    27 0.296 0.367 0.458
    31 0.276 0.341 0.422
    35 0.260 0.319 0.393
    39 0.246 0.301 0.369
    43 0.234 0.289 0.349
    47 0.224 0.273 0.332
    51 0.215 0.261 0.317
    ]
    x_knots = knots[:, 1]
    if value == 0.95
        y_knots = knots[:, 2]
    elseif value == 0.975
        y_knots = knots[:, 3]
    elseif value == 0.99
        y_knots = knots[:, 4]
    else
        error("Value=$value is not allowed")
    end
    itp = Spline1D(x_knots, y_knots, k=1, bc="extrapolate")
    evaluate(itp, count)
end

function getMeanAbsDeviationQuantile(value, count)
    if value != 0.95 && value != 0.975
        error("Value=$value is not allowed")
    end
    const knots = [
    2 12.71;
    3 3.45;
    4 2.16;
    5 1.66;
    6 1.40;
    7 1.21;
    8 1.09;
    9 1.00;
    10 0.93;
    11 0.87;
    12 0.82;
    13 0.78;
    14 0.75;
    20 0.71;
    25 0.60;
    30 0.48;
    40 0.41;
    60 0.33;
    120 0.23;
    ]
    x_knots = knots[:, 1]
    y_knots = knots[:, 2]
    itp = Spline1D(x_knots, y_knots, k=1, bc="extrapolate")
    evaluate(itp, count)
end

function studentQuatile(value, freedomDegreeNumber)
    d = TDist(freedomDegreeNumber)
    quantile(d, value)
end

function normalQuantile(value)
    d = Normal()
    quantile(d, value)
end

end


text:
module Evaluators
export Evaluator

using FactorizationMachines: FMMatrix, FMFloat
using FactorizationMachines.Common
using FactorizationMachines.Predictors: FMPredictor, predict!
using FactorizationMachines.Tasks: PredictorTask

abstract Evaluator

immutable SquaredErrorEvaluator <: Evaluator
    stat::Function

    SquaredErrorEvaluator(; stat::Function = x -> sqrt(mean(x))) = new(stat)
end
const rmse = SquaredErrorEvaluator

function evaluate{T<:PredictorTask}(evaluator::SquaredErrorEvaluator, predictor::FMPredictor{T},
        X::FMMatrix, y::Vector{FMFloat})
    predictions = zeros(X.n)
    evaluate!(evaluator, predictor, X, y, predictions)
end

function evaluate!{T<:PredictorTask}(evaluator::SquaredErrorEvaluator, predictor::FMPredictor{T},
        X::FMMatrix, y::Vector{FMFloat}, predictions::Vector{FMFloat})
    @time predict!(predictor, X, predictions)
    err = [Common.sqerr(predictions[i], y[i]) for i in 1:length(y)]
    evaluator.stat(err .* err)
end

immutable HeavisideEvaluator <: Evaluator
    stat::Function

    HeavisideEvaluator(; stat::Function = x -> mean(x)) = new(stat)
end
const heaviside = HeavisideEvaluator

function evaluate{T<:PredictorTask}(evaluator::HeavisideEvaluator, predictor::FMPredictor{T},
        X::FMMatrix, y::Vector{FMFloat})
    predictions = zeros(X.n)
    evaluate!(evaluator, predictor, X, y, predictions)
end

function evaluate!{T<:PredictorTask}(evaluator::HeavisideEvaluator, predictor::FMPredictor{T},
        X::FMMatrix, y::Vector{FMFloat}, predictions::Vector{FMFloat})
    @time predict!(predictor, X, predictions)
    err = [Common.heaviside(predictions[i], y[i]) for i in 1:length(y)]
    evaluator.stat(err)
end

end


text:
# constructors.jl

const tupleTMs = (:TaylorModel1, :RTaylorModel1)
const NumberNotSeries = TaylorSeries.NumberNotSeries

#=
Structs `TaylorModel1` and `RTaylorModel1` are essentially identical, except
the way the remainder is handled; the remainder for `TaylorModel1`
must contain 0.
=#
for TM in tupleTMs
    @eval begin
        struct $(TM){T,S} <: AbstractSeries{T}
            pol  :: Taylor1{T}    # polynomial approx (of order `ord`)
            rem  :: Interval{S}   # remainder
            x0   :: Interval{S}   # expansion point
            dom  :: Interval{S}   # interval of interest

            # Inner constructor
            function $(TM){T,S}(pol::Taylor1{T}, rem::Interval{S},
                    x0::Interval{S}, dom::Interval{S}) where {T,S}
                if $(TM) == TaylorModel1
                    @assert zero(S) ∈ rem && x0 ⊆ dom
                else
                    @assert x0 ⊆ dom
                end
                return new{T,S}(pol, rem, x0, dom)
            end
        end

        # Outer constructors
        $(TM)(pol::Taylor1{T}, rem::Interval{S},
            x0, dom::Interval{S}) where {T,S} = $(TM){T,S}(pol, rem, interval(x0), dom)

        # Constructor just chainging the remainder
        $(TM)(u::$(TM){T,S}, Δ::Interval{S}) where {T,S} = $(TM){T,S}(u.pol, Δ, u.x0, u.dom)

        # Short-cut for independent variable
        $(TM)(ord::Integer, x0, dom::Interval{T}) where {T} =
            $(TM)(x0 + Taylor1(eltype(x0), ord), zero(dom), interval(x0), dom)

        # Short-cut for a constructor expanding around midpoint by default
        $(TM)(ord::Integer, dom::Interval{T}) where {T} =
            $(TM)(ord, Interval(mid(dom)), dom)

        # Short-cut for a constant TM
        $(TM)(a::T, ord::Integer, x0::Interval{S}, dom::Interval{S}) where
            {T,S} = $(TM)(Taylor1([a], ord), zero(dom), x0, dom)

        # Functions to retrieve the order and remainder
        @inline get_order(tm::$TM) = get_order(tm.pol)
        @inline remainder(tm::$TM) = tm.rem
        @inline polynomial(tm::$TM) = tm.pol
        @inline domain(tm::$TM) = tm.dom
        @inline expansion_point(tm::$TM) = tm.x0
    end
end

@doc doc"""
    TaylorModel1{T,S}

Absolute Taylor model in 1 variable, providing a rigurous polynomial approximation
given by a Taylor polynomial `pol` (around `x0`) and an absolute remainder
`rem` for a function `f(x)` in one variable, valid in the interval `dom`.
This corresponds to definition 2.1.3 of Mioara Joldes' thesis.

Fields:
- `pol`: polynomial approximation, represented as `TaylorSeries.Taylor1`
- `rem`: the interval bound
- `x0` : expansion point
- `dom`: domain, interval over which the Taylor model is defined / valid

The approximation ``f(x) = p(x) + \Delta`` is satisfied for all
``x\in \mathcal{D}`` (``0\in \Delta``); `n` is the order (degree)
of the polynomial ``p(x)=\sum_{i=0}^n p_i (x - x_0)^i``.

""" TaylorModel1

@doc doc"""
    RTaylorModel1{T,S}

Relative Taylor model in 1 variable, providing a rigurous polynomial approximation
given by a Taylor polynomial `pol` (around `x0`) and a relative remainder
`rem` for a function `f(x)` in one variable, valid in the interval `dom`.
This corresponds to definition 2.3.2 of Mioara Joldes' thesis.

Fields:
- `pol`: polynomial approximation, represented as `TaylorSeries.Taylor1`
- `rem`: the interval bound
- `x0` : expansion point
- `dom`: domain, interval over which the Taylor model is defined / valid

The approximation ``f(x) = p(x) + \delta (x - x_0)^{n+1}`` is satisfied for all
``x\in \mathcal{D}``; `n` is the order (degree) of the polynomial
``p(x)=\sum_{i=0}^n p_i (x - x_0)^i``.

""" RTaylorModel1

# TaylorModelN's struct
"""
    TaylorModelN{N,T,S}

Taylor Models with absolute remainder for `N` independent variables.

"""
struct TaylorModelN{N,T,S} <: AbstractSeries{T}
    pol  :: TaylorN{T}        # polynomial approx (of order `ord`)
    rem  :: Interval{S}       # remainder
    x0   :: IntervalBox{N,S}  # expansion point
    dom  :: IntervalBox{N,S}  # interval of interest

    # Inner constructor
    function TaylorModelN{N,T,S}(pol::TaylorN{T}, rem::Interval{S},
            x0::IntervalBox{N,S}, dom::IntervalBox{N,S}) where {N,T<:NumberNotSeries,S<:Real}

        @assert N == get_numvars()
        @assert zero(S) ∈ rem && x0 ⊆ dom

        return new{N,T,S}(pol, rem, x0, dom)
    end
end

# Outer constructors
TaylorModelN(pol::TaylorN{T}, rem::Interval{S}, x0::IntervalBox{N,S},
    dom::IntervalBox{N,S}) where {N,T,S} = TaylorModelN{N,T,S}(pol, rem, x0, dom)

# Constructor for just changing the remainder
TaylorModelN(u::TaylorModelN{N,T,S}, Δ::Interval{S}) where {N,T,S} =
    TaylorModelN{N,T,S}(u.pol, Δ, u.x0, u.dom)

# Short-cut for independent variable
TaylorModelN(nv::Integer, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) where {N,T} =
    TaylorModelN(x0[nv] + TaylorN(Interval{T}, nv, order=ord), zero(dom[1]), x0, dom)

# Short-cut for a constant
TaylorModelN(a::Interval{T}, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) where {N,T} =
    TaylorModelN(TaylorN(a, ord), zero(dom[1]), x0, dom)
TaylorModelN(a::T, ord::Integer, x0::IntervalBox{N,T}, dom::IntervalBox{N,T}) where {N,T} =
    TaylorModelN(TaylorN(a, ord), zero(dom[1]), x0, dom)

# Functions to retrieve the order and remainder
@inline get_order(tm::TaylorModelN) = tm.pol.order
@inline remainder(tm::TaylorModelN) = tm.rem
@inline polynomial(tm::TaylorModelN) = tm.pol
@inline domain(tm::TaylorModelN) = tm.dom
@inline expansion_point(tm::TaylorModelN) = tm.x0
@inline get_numvars(::TaylorModelN{N,T,S}) where {N,T,S} = N

"""
    TMSol{N,T,V1,V2,M}

Structure containing the solution of a validated integration.

# Fields
    `time :: AbstractVector{T}`  Vector containing the expansion time of the `TaylorModel1` solutions

    `fp   :: AbstractVector{IntervalBox{N,T}}`  IntervalBox vector representing the flowpipe

    `xTMv :: AbstractMatrix{TaylorModel1{TaylorN{T},T}}`  Matrix whose entry `xTMv[i,t]` represents
    the `TaylorModel1` of the i-th dependent variable, obtained at time time[t].
"""
struct TMSol{N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
        M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
    time :: V1
    fp   :: V2
    xTM  :: M

    function TMSol(time::V1, fp::V2, xTM::M) where
            {N,T<:Real,V1<:AbstractVector{T},V2<:AbstractVector{IntervalBox{N,T}},
            M<:AbstractMatrix{TaylorModel1{TaylorN{T},T}}}
        @assert length(time) == length(fp) == size(xTM,2) && N == size(xTM,1)
        return new{N,T,V1,V2,M}(time, fp, xTM)
    end
end

@inline expansion_point(a::TMSol) = getfield(a,:time)
@inline expansion_point(a::TMSol, n::Int) = getindex(getfield(a,:time),n)
@inline flowpipe(a::TMSol) = getfield(a,:fp)
@inline flowpipe(a::TMSol, n::Int) = getindex(getfield(a,:fp),n)
@inline get_xTM(a::TMSol) = getfield(a,:xTM)
@inline get_xTM(a::TMSol, n::Int) = getindex(getfield(a,:xTM),:,n)
@inline domain(a::TMSol) = domain.(getindex(getfield(a, :xTM), 1, :)) # vector!
@inline domain(a::TMSol, n::Int) = domain(getindex(getfield(a, :xTM), 1, n))
@inline get_numvars(::TMSol{N,T,V1,V2,M}) where {N,T,V1,V2,M} = N


text:
using Dates
using Printf
using Libdl
using Sockets
function getCfg(modulename, cfgJl::String )
    println( cfgJl )
    evs =""
    cfgObj =open( cfgJl,"r")
    # skip( cfgObj,3)
    for line in Base.readlines(cfgObj)
        # println( line )
        if  strip( line ) != "" && line[1:1] != "#"
            evs = evs *line *"; "
        end
    end

    reps = Base.Meta.parse( evs )
    Base.eval(modulename,reps )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
@testset "1091.shortest-path-in-binary-matrix.jl" begin
    @test shortest_path_binary_matrix([
        0 0 0
        1 1 0
        1 1 0
    ]) == 4
    @test shortest_path_binary_matrix([
        0 1
        1 0
    ]) == 2
end