text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: personalize_runtime

using AWS.Compat
using AWS.UUIDs
"""
    GetPersonalizedRanking()

Re-ranks a list of recommended items for the given user. The first item in the list is deemed the most likely item to be of interest to the user.  The solution backing the campaign must have been created using a recipe of type PERSONALIZED_RANKING.

# Required Parameters
- `campaignArn`: The Amazon Resource Name (ARN) of the campaign to use for generating the personalized ranking.
- `inputList`: A list of items (by itemId) to rank. If an item was not included in the training dataset, the item is appended to the end of the reranked list. The maximum is 500.
- `userId`: The user for which you want the campaign to provide a personalized ranking.

# Optional Parameters
- `context`: The contextual metadata to use when getting recommendations. Contextual metadata includes any interaction information that might be relevant when getting a user's recommendations, such as the user's current location or device type.
- `filterArn`: The Amazon Resource Name (ARN) of a filter you created to include or exclude items from recommendations for a given user.
"""

get_personalized_ranking(campaignArn, inputList, userId; aws_config::AWSConfig=global_aws_config()) = personalize_runtime("POST", "/personalize-ranking", Dict{String, Any}("campaignArn"=>campaignArn, "inputList"=>inputList, "userId"=>userId); aws_config=aws_config)
get_personalized_ranking(campaignArn, inputList, userId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = personalize_runtime("POST", "/personalize-ranking", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("campaignArn"=>campaignArn, "inputList"=>inputList, "userId"=>userId), args)); aws_config=aws_config)

"""
    GetRecommendations()

Returns a list of recommended items. The required input depends on the recipe type used to create the solution backing the campaign, as follows:   RELATED_ITEMS - itemId required, userId not used   USER_PERSONALIZATION - itemId optional, userId required    Campaigns that are backed by a solution created using a recipe of type PERSONALIZED_RANKING use the API.

# Required Parameters
- `campaignArn`: The Amazon Resource Name (ARN) of the campaign to use for getting recommendations.

# Optional Parameters
- `context`: The contextual metadata to use when getting recommendations. Contextual metadata includes any interaction information that might be relevant when getting a user's recommendations, such as the user's current location or device type.
- `filterArn`: The ARN of the filter to apply to the returned recommendations. For more information, see Using Filters with Amazon Personalize. When using this parameter, be sure the filter resource is ACTIVE.
- `itemId`: The item ID to provide recommendations for. Required for RELATED_ITEMS recipe type.
- `numResults`: The number of results to return. The default is 25. The maximum is 500.
- `userId`: The user ID to provide recommendations for. Required for USER_PERSONALIZATION recipe type.
"""

get_recommendations(campaignArn; aws_config::AWSConfig=global_aws_config()) = personalize_runtime("POST", "/recommendations", Dict{String, Any}("campaignArn"=>campaignArn); aws_config=aws_config)
get_recommendations(campaignArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = personalize_runtime("POST", "/recommendations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("campaignArn"=>campaignArn), args)); aws_config=aws_config)


text:
using ArgParse
using ReadCoverage

println("####################\n   ReadCoverage.jl\n####################")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "relcov"
            help = "Calculates relative gene body coverage from a BAM file.\n(Almost same as genebody_coverage.py in RSeQC.)"
            action = :command
        "abcov"
            help = "Cacluates absolute gene body coverage from a BAM file."
            action = :command
        "coverage"
            help = "Cacluates absolute gene body coverage from a BAM file."
            action = :command
    end

    @add_arg_table! s["relcov"] begin
        "path_bam"
            help = "Path to a BAM file."
            arg_type = String
            required = true
        "path_bed12"
            help = "Path to a reference gene model file with BED12 format."
            arg_type = String
            required = true
        "output_prefix"
            help = "Prefix for output files. If not specified or set to \"\" (defalut), no output files are saved."
            arg_type = String
            required = true
        "--max_depth"
            help = "Maximum depth for a position. Depth more than `max_depth` is cut to `max_depth`. If this keyword is set to 0 (defalut), no cut is occurred."
            arg_type = Int
            default = 0
    end

    @add_arg_table! s["abcov"] begin
        "path_bam"
            help = "Path to a BAM file."
            arg_type = String
            required = true
        "path_bed12"
            help = "Path to a reference gene model file with BED12 format."
            arg_type = String
            required = true
        "output_prefix"
            help = "refix for output files. If not specified or set to \"\" (defalut), no output files are saved."
            arg_type = String
            required = true
        "--bin_size"
            help = "Bin size for transcripts (Defalut: 100)."
            arg_type = Int
            default = 100
    end
    @add_arg_table! s["coverage"] begin
        "path_bam"
            help = "Path to a BAM file."
            arg_type = String
            required = true
        "chrom"
            help = "Chromosome name for a genomic interval."
            arg_type = String
            required = true
        "leftpos"
            help = "Left position for a genomic interval."
            arg_type = Int64
            required = true
        "rightpos"
            help = "Right position for a genomic interval."
            arg_type = Int64
            required = true
        "output_prefix"
            help = "Prefix for output files. If not specified or set to \"\" (defalut), no output files are saved."
            arg_type = String
            required = true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    if parsed_args["%COMMAND%"] == "relcov"
        relative_genebodycoverage(
            parsed_args["relcov"]["path_bam"],
            parsed_args["relcov"]["path_bed12"],
            output_prefix = parsed_args["relcov"]["output_prefix"],
            max_depth=parsed_args["relcov"]["max_depth"]
        );
    elseif parsed_args["%COMMAND%"] == "abcov"
        absolute_genebodycoverage(
            parsed_args["abcov"]["path_bam"],
            parsed_args["abcov"]["path_bed12"],
            output_prefix = parsed_args["abcov"]["output_prefix"],
            bin_size=parsed_args["abcov"]["bin_size"]
        );
    elseif parsed_args["%COMMAND%"] == "coverage"
        readcoverage_bam(
            parsed_args["coverage"]["path_bam"],
            parsed_args["coverage"]["chrom"],
            parsed_args["coverage"]["leftpos"],
            parsed_args["coverage"]["rightpos"],
            output_prefix = parsed_args["coverage"]["output_prefix"]
        );
    end

end

main()

ext: jl


text:
##################### Binary Individual Adaptation ####################

#################### Types and Constructors ####################

mutable struct BIATune <: SamplerTune
  logf::Union{Function, Missing}
  A::Vector{Float64}
  D::Vector{Float64}

  epsilon::Float64
  decay::Float64
  target::Float64

  iter::Int

  BIATune() = new()

  function BIATune(x::Vector, logf::Union{Function, Missing};
                   A::Vector{Float64} = fill(1, size(x)) ./ length(x),
                   D::Vector{Float64} = fill(1, size(x)) ./ length(x),
                   epsilon::Real = 0.01 / length(x), decay::Real = 0.55,
                   target::Real = 0.45)
    new(logf, A, D, epsilon, decay, target, 0)
  end
end

BIATune(x::Vector; args...) =
   BIATune(x, missing; args...)

const BIAVariate = SamplerVariate{BIATune}

function validate(v::BIAVariate)
  n = length(v)
  epsilon = v.tune.epsilon
  0.0 < epsilon < 0.5 ||
    throw(ArgumentError("epsilon is not in (0, 0.5)"))
  length(v.tune.A) == n && all(epsilon .< v.tune.A .< 1 - epsilon) ||
    throw(ArgumentError("A must contain $n probabilities in ($epsilon, $(1 - epsilon))"))
  length(v.tune.D) == n && all(epsilon .< v.tune.D .< 1 - epsilon) ||
    throw(ArgumentError("D must contain $n probabilities in ($epsilon, $(1 - epsilon))"))
  0.5 < v.tune.decay <= 1.0 ||
    throw(ArgumentError("decay is not in (0.5, 1]"))
  0.0 < v.tune.target < 1.0 ||
    throw(ArgumentError("target is not in (0, 1)"))
  validatebinary(v)
end

#################### Sampler Constructor ####################
"""
    BIA(params::ElementOrVector{Symbol}; args...)
Construct a `Sampler` object for BIA sampling. Parameters are assumed to have
   binary numerical values (0 or 1).

Returns a `Sampler{BIATune}` type object.

* `params`: stochastic node(s) to be updated with the sampler.

* `args...`: additional keyword arguments to be passed to the `BIAVariate` constructor.
"""
function BIA(params::ElementOrVector{Symbol}; args...)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block)
    v = SamplerVariate(block; args...)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, BIATune())
end

#################### Sampling Functions ####################

sample!(v::BIAVariate) = sample!(v, v.tune.logf)
"""
    sample!(v::BIAVariate, logf::Function)
Draw one sample from a target distribution using the BIA sampler. Parameters
 are assumed to have binary numerical values (0 or 1).

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::BIAVariate, logf::Function)
  tune = v.tune
  tune.iter += 1

  p = length(v)

  ## proposal
  x = v[:]
  added = zeros(Int8, p)
  deleted = zeros(Int8, p)
  q_ratio = 1.0
  for j in 1:p
    if v.value[j] == 0.0
      if rand() < tune.A[j]
        x[j] = 1.0
        added[j] = 1
        q_ratio *= tune.D[j] / tune.A[j]
      end
    else
      if rand() < tune.D[j]
        x[j] = 0.0
        deleted[j] = 1
        q_ratio *= tune.A[j] / tune.D[j]
      end
    end
  end

  ## M-H acceptance probability
  alpha = min(1.0, exp(logf(x) - logf(v.value)) * q_ratio)

  ## Adapt tuning parameters
  for j in 1:p
    ## new A[j]
    C = log((tune.A[j] - tune.epsilon) / (1.0 - tune.A[j] - tune.epsilon)) +
        tune.iter^(-tune.decay) * added[j] * (alpha - tune.target)
    tune.A[j] = (exp(C) * (1.0 - tune.epsilon) + tune.epsilon) / (1.0 + exp(C))

    ## new D[j]
    C = log((tune.D[j] - tune.epsilon) / (1.0 - tune.D[j] - tune.epsilon)) +
        tune.iter^(-tune.decay) * deleted[j] * (alpha - tune.target)
    tune.D[j] = (exp(C) * (1.0 - tune.epsilon) + tune.epsilon) / (1.0 + exp(C))
  end

  ## Accept or reject proposal
  if rand() < alpha
    v[:] = x
  end

  v
end


text:
using Flux

W = rand(2, 5)
b = rand(2)

predict(x) = W*x .+ b

function loss(x, y)
  y1 = predict(x)
  sum((y .- y1).^2)
end

x, y = rand(5), rand(2)
@show(loss(x, y))

gs = gradient(() -> loss(x, y), params(W, b))
W1 = gs[W]
W .-= 0.1 .* W1
@show(loss(x, y))


text:
module SIMeasurements

export Dimensions, Unit, Quantity, unit, conversion

using FixedSizeArrays: FixedVectorNoTuple
import Base: *, ^, /, //, +, -

include("units.jl")
include("quantity.jl")
include("quantity_ops.jl")
end # module


text:
module MDL

using LinearAlgebra
import Statistics: mean
import StatsBase: fit, Histogram
import ..LMCLUS: Manifold, indim, outdim, projection, threshold, points, hasfullbasis, project, LMCLUSResult

# Various types for MDL calculation
abstract type MethodType end

struct Raw <: MethodType end
struct ZeroDim <: MethodType end

struct Uniform <: MethodType end
struct Gaussian <: MethodType end
struct Empirical <: MethodType end
struct OptimalQuant <: MethodType end    # ICPR-2016 Eq. 6
struct SizeIndependent <: MethodType end # ICPR-2016 Eq. 8

"""Optimal number of bins given constant C over intervals"""
function optbins(intervals::Vector{T}, C::T) where T<:AbstractFloat
    ns = ceil.(intervals * exp( (C - sum(intervals))/length(intervals) ), digits=0)
    nsi = zeros(UInt, size(ns))
    for i in 1:length(nsi)
        nsi[i] = if ns[i] >= typemax(UInt)
            typemax(UInt)
        else
            round(UInt, ns[i])
        end
    end
    return nsi
end

"""Variance of a uniform distribution over interval"""
univar(interval::Vector{T}) where T<:AbstractFloat = (1//12)*(interval).^2

"""Quantizaton error"""
quanterror(intervals, N_k_opt) = sum(univar(intervals./N_k_opt))

"""Optimal quantization of the interval"""
function optquant(intervals::Vector{T}, ɛ::Real; tot::Int=1000, tol::Real=1e-6) where {T<:AbstractFloat}
    intervals[isnan.(intervals)]  .= eps(T) # remove nans
    intervals[intervals .< eps()] .= eps(T) # remove 0s
    intervals[isinf.(intervals)]  .= 1.0   # remove inf

    # Setup C bounds
    K = length(intervals)
    C = zero(T)
    Cmin = zero(T)
    Cmax = broadcast(+, log.(typemax(UInt)./intervals).*K, sum(log.(intervals))) |> minimum

    i = 1
    N_k_opt = ones(UInt, K)
    err_opt = Inf
    while i < tot
        C = Cmax - (Cmax - Cmin)/2
        N_k_opt = optbins(intervals, C)
        err_opt = quanterror(intervals, N_k_opt)
        err_diff = err_opt - ɛ^2
        if abs(err_diff) < tol
            break # error is small enough to stop search
        elseif err_diff > 0.
            # error is too big
            Cmin = C
        elseif err_diff < 0.
            # error is too small
            Cmax = C
        end
        i += 1
    end
    return N_k_opt, err_opt, C, i
end

"""Calculate LM cluster bounding box using PCA

    Ref: Jonathon Shlens, A Tutorial on Principal Component Analysis, 2005
"""
function boundingbox(X, m)
    n, l = size(X)

    # PCA data
    F = svdfact(X./sqrt(n))

    # Take an orthogonal compliment subspace (OCS) basis
    r = if m < min(n,l)
        (m+1):min(n,l)
    else
        1:min(n,l)
    end
    BC = F[:U][:,r]

    # Project data to OCS
    Y = BC'*X

    # Calculate data spread in OCS
    Ymin = vec(minimum(Y, dims=2))
    Ymax = vec(maximum(Y, dims=2))
    range = abs.(Ymax - Ymin)

    return range, Y, Ymin, Ymax
end

function boundingbox(Y::AbstractMatrix{T}) where T<:AbstractFloat
    Ymin = vec(minimum(Y, dims=2))
    Ymax = vec(maximum(Y, dims=2))
    Yrange = abs.(Ymax - Ymin)
    return Yrange, Ymin, Ymax
end

#=
    Minimum Description Length Encoding
=#

"""Raw representation of the data"""
modeldl(::Type{Raw}, M::Manifold, X::AbstractMatrix, P::Int) = P*length(X)
datadl(::Type{Raw}, aggs...) = 0

function raw(M::Manifold, Pm::Int)
    return Pm*size(M)*indim(M)
end

"""Zero dimensional model description length: L(ZD)

    From: O. Georgieva, K. Tschumitschew, and F. Klawonn,
    “Cluster validity measures based on the minimum description length principle”
"""
function modeldl(::Type{ZeroDim}, M::Manifold, X::AbstractMatrix, P::Int)
    return P*size(X,1)
end

"""Zero dimensional data description length: L(X|ZD)

    Ref: O. Georgieva, K. Tschumitschew, and F. Klawonn, “Cluster validity
    measures based on the minimum description length principle”
"""
function datadl(::Type{ZeroDim}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    N = size(X,1)  # space dimension
    μ = mean(C)    # manifold translation
    n = size(C)  # size of cluster

    E = 0.0
    for i in 1:n
        for j in 1:N
            E += log2(abs(X[j,i] - μ[j]) + 1)
        end
    end
    return round(Int, E)
end

#=
    Various models of MDL calculation (model MDL)
=#

#=
""" General model description length (v1)

    Only encodes linear manifold basis vectors (no orthogonal complement space)
"""
function modeldl{MT<:MethodType}(::Type{MT}, C::Manifold, X::Matrix, P::Int)
    N = size(X,1)  # space dimension
    M = outdim(C)   # manifold dimension
    return P*(N + M*(N - (M+1)>>1))
end
=#

""" General model description length (v2)

    Encodes linear manifold & orthogonal complement space basis vectors
"""
function modeldl(::Type{MT}, C::Manifold, X::AbstractMatrix, P::Int) where MT <: MethodType
    N = size(X,1)  # space dimension
    return if outdim(C) != 0
        (P*N*(N+1))>>1
    else # for 0D manifold no basis encoding required
        P*N
    end
end

"""Size independent data description length: L(X|SI)
"""
function datadl(::Type{SizeIndependent}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    M = outdim(C)  # manifold dimension
    N = size(X,1)  # space dimension

    # Project points to the orthogonal compliment space (OCS) of the manifold
    Y = project(C, X)

    # Get OCS bounding box
    rngs, _ = boundingbox(Y)

    # Estimate number of bins required for each component of OSC
    bins, _ = optquant(rngs, ɛ, tot=tot, tol=tol)

    return convert(Int, P*(N + 2*sum(bins)))
end

#=
    Various models of MDL calculation (data MDL)
=#

""" Data MDL: Uniform encoding

    Points in the orthogonal complement space of the linear manifold cluster
    are uniformly distributed around the cluster manifold.
"""
function datadl(::Type{Uniform}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    M = outdim(C)  # manifold dimension
    n = size(X,2)  # size of cluster

    # Encoding size of point projected on manifold
    PR = P*M

    # Entropy of points in the orthogonal complement space
    E = -float(n)*log(threshold(C)[1])

    # Number of bits of two parts for every point
    return round(Int, (PR + E)*n)
end

""" Data MDL: Gaussian encoding

    Points in the orthogonal complement space of the linear manifold cluster
    are normally distributed around the cluster manifold.
"""
function datadl(::Type{Gaussian}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    N, n = size(X)  # space dimension & size
    M = outdim(C)   # manifold dimension
    μ = mean(C)     # manifold translation
    B = projection(C) # manifold basis

    # Encoding size of point projected on manifold
    PR = P*M

    # Change points basis to of orthogonal complement subspace
    OP = (Matrix{T}(I, N, N) - B*B')*(X.-μ)

    # multivariate normal distribution entropy
    Σ = OP*OP'
    Σd = det(Σ)
    E = (M/2)*(1+log(2π)) + (Σd > eps() ? log(Σd) : 0.0)/2

    # Number of bits of two parts for every point
    return round(Int, (PR + E)*n)
end

""" Data MDL: Empirical entropy encoding

    Distribution of points in the orthogonal complement space of the linear manifold cluster,
    constructed by optimal quantizing with a fixed length bin, is used to
    calculate entropy, thus number of bits.
"""
function datadl(::Type{Empirical}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    N, n = size(X) # space dimension & size
    M = outdim(C)  # manifold dimension

    # Encoding size of point projected on manifold
    PR = P*M

    # Project points to the orthogonal compliment space (OCS) of the manifold
    Y = project(C, X)

    # Get OCS bounding box
    rngs, Ymin, Ymax = boundingbox(Y)

    # Estimate number of bins required for each competent of OSC
    if ɛ < 1. # using optimal quantization
        bins, _ = optquant(rngs, ɛ, tot=tot, tol=tol)
    else      # constant number of bins per component
        bins = fill(round(Int,ɛ), N-M)
    end

    # Cumulative entropy of data spread in OSC (calculated from quantization)
    E = 0.0
    for i in 1:N-M
        Yb = range(Ymin[i], stop=Ymax[i], length=bins[i]+1)
        H = fit(Histogram, vec(Y[i,:]), Yb, closed=:left)
        E += -sum(map(x-> x > 0. ? (x/(n-1))*log2(x/(n-1)) : 0., H.weights))
    end

    # Number of bits of two parts for every point
    return round(Int, (PR + E)*n)
end

""" Data MDL: Optimal quantization encoding

    Distribution of points in the orthogonal complement space of the linear manifold cluster,
    constructed by optimal quantizing with a fixed length bin, is used to
    calculate entropy, thus number of bits.
"""
function datadl(::Type{OptimalQuant}, C::Manifold{T}, X::AbstractMatrix{T},
                P::Int, ɛ::Real, tot::Int, tol::Real) where T<:AbstractFloat
    M = outdim(C)  # manifold dimension
    n = size(X,2)  # size of cluster

    # Encoding size of a point projected on manifold
    PR = P*M

    # Project points to the orthogonal compliment space (OCS) of the manifold
    Y = project(C, X)

    # Get OCS bounding box
    rngs, _ = boundingbox(Y)

    # Estimate number of bins required for each component of OSC
    bins, sqerr, C = optquant(rngs, ɛ, tot=tot, tol=tol)

    # Number of bins required to encode data spread in OSC given optimal quantization
    # C = Σ_k log(N_k)
    # where N_k number of bins in quantized interval of component K of OSC
    E = C/log(2)

    # Number of bits of two parts for every point
    return round(Int, (PR + E)*n)
end

#=
    Main MDL call
=#

"Calculate MDL for the manifold"
function calculate(::Type{MT}, M::Manifold{T}, X::AbstractMatrix{T}, Pm::Int, Pd::Int;
                   ɛ::Real = 1e-3, tot::Int=1000, tol::Real=1e-18) where {MT<:MethodType, T<:AbstractFloat}
    return modeldl(MT, M, X, Pm) + datadl(MT, M, X, Pd, ɛ, tot, tol)
end

"Calculate MDL for the clustering"
function calculate(::Type{MT}, Ms::Vector{Manifold}, X::AbstractMatrix{T}, Pm::Int, Pd::Int;
                   ɛ::Real=1e-3, tot::Int=1000, tol::Real=1e-18) where {MT<:MethodType, T<:AbstractFloat}
    return sum(calculate(MT, m, X[:,points(m)], Pm, Pd, ɛ=ɛ, tot=tot, tol=tol) for m in Ms)
end

function calculate!(::Type{MT}, C::LMCLUSResult, X::AbstractMatrix{T}, Pm::Int, Pd::Int;
                   ɛ::Real=1e-3, d=-1) where {MT<:MethodType, T<:AbstractFloat}
    mdl = 0
    for (i, m) in enumerate(C.manifolds)
        pts = points(m)
        MX = view(X, :, pts)
        μ = mean(MX, dims=2)
        Svd = svd(MX .- μ, full=true)
        U = Svd.U
        mtmp = Manifold((d >= 0 ? d : m.d), vec(μ), U, pts)
        mdl += calculate(MT, mtmp, MX, Pm, Pd, ɛ=ɛ)
    end
    mdl
end

function compression(::Type{MT}, C::LMCLUSResult, X::AbstractMatrix{T}, Pm::Int,
                     Pd::Int; kwargs...) where {MT<:MethodType, T<:AbstractFloat}
    raw = Pm * length(X)
    mdl = calculate!(MT, C, X, Pm, Pd; kwargs...)
    raw/mdl
end

end


text:
using DataFrames
using PowerFlowData
using Tables
using Test

@testset "PowerFlowData.jl" begin

    @testset "Tables interface" begin
        struct TestRecs <: PowerFlowData.Records
            x1::Vector{Int}
            yy::Vector{Float64}
        end
        recs = TestRecs([1, 2, 3], [3.14, 1.72, 0.05])

        @test Tables.istable(TestRecs)

        @test Tables.columnaccess(TestRecs)
        @test Tables.columns(recs) === recs
        @test Tables.columnnames(recs) == (:x1, :yy)

        @test recs.x1 == [1, 2, 3]
        @test Tables.getcolumn(recs, :x1) == [1, 2, 3]
        @test Tables.getcolumn(recs, 1) == [1, 2, 3]

        df = DataFrame(recs)
        @test df isa DataFrame
        @test size(df) == (3, 2)

        for T in (Buses, Loads, Generators, Branches)
            @test T <: PowerFlowData.Records
            @test Tables.istable(T)
        end

        caseid = CaseID(0, 100.0)
        @test caseid[1] == caseid[:ic]
        @test caseid[2] == caseid[:sbase]
        @test NamedTuple(caseid) == (ic=0, sbase=100.0)
    end

    @testset "v30 file" begin
        net1 = parse_network("testfiles/synthetic_data_v30.raw")
        @test net1 isa Network

        caseid = net1.caseid
        @test caseid.ic == 0
        @test caseid.sbase == 100.0

        # Test first and last columns are parsed as expected
        buses = net1.buses
        @test buses.i == [111, 112, 113]
        @test buses.owner == [1, 2, 2]

        # Test string column as expected
        @test buses.name[2] == "D2JK  "

        loads = net1.loads
        @test loads.i == [111, 113]
        @test loads.owner == [1, 2]

        gens = net1.generators
        @test gens.i == [111, -112, 113]
        @test gens.fi == [1.0, 1.0, 1.0]
        @test gens.id[1] == "ST"

        branches = net1.branches
        @test branches.i == [111, 111, 112]
        @test branches.j == [112, -113, 113]  # negative numbers should be allowed
        @test branches.fi == [1.0, 1.0, 1.0]
        @test branches.ckt[1] == "3 "
    end

    @testset "v29 file" begin
        net2 = parse_network("testfiles/synthetic_data_v29.raw")

        caseid = net2.caseid
        @test caseid.ic == 0
        @test caseid.sbase == 100.0

        # Test first and last columns are parsed as expected
        buses = net2.buses
        @test buses.i == [1, 222222]
        @test buses.owner == [1, 7]

        # Test string column as expected
        @test buses.name[2] == "PRPR C D    "

        loads = net2.loads
        @test loads.i == [1, 222222]
        @test loads.owner == [1, 7]

        gens = net2.generators
        @test gens.i == [104]    # first col
        @test gens.fi == [1.0]   # last col
        @test gens.id[1] == "1 " # string col

        branches = net2.branches
        @test branches.i == [1, 2, 222222]
        @test branches.j == [-543210, 9, 333333]  # negative numbers should be allowed
        @test branches.fi == [1.0, 1.0, 1.0]
        @test branches.ckt[2] == "6 "

    end
end


text:
using st16SseqJuliaTools

if abspath(PROGRAM_FILE) == @__FILE__
    mergejc(ARGS[1], ARGS[2])
end


text:
module Bezier

export bezier

# Linear interpolation between two points
function bezier(t::Real, P1::Array{T}, P2::Array{T})::T where T <: AbstractFloat
    if 0 ≤ t ≤ 1
        return P1 + t * (P2 - P1)
    else
        # t must be bounded by [0, 1]
        # TODO is this the right way to notify users?
        throw(DomainError(t, "t must be bounded by [0, 1]"))
    end
end

# Higher-order interpolation
function bezier(t::Real, points::Vector{Vector{T}})::T where T <: AbstractFloat
    dimension = size(points, 1)
    if dimension ≤ 1
        return points
    else
        interpolated_points = Array{Vector{T}}(dimension - 1)
        for i in 2:dimension
            interpolated_points[i - 1] = bezier(t, points[i - 1], points[i])
        end
        return bezier(t, interpolated_points)
    end
end

# Each row of the matrix is a point
function bezier(t::Real, points::Matrix{T})::T where T <: AbstractFloat
    rows, cols = size(points)
    if rows ≤ 1
        return points
    else
        interpolated_points = Array{T}(rows - 1, cols)
        for r in 2:rows
            interpolated_point = bezier(t, points[r - 1, :], points[r, :])
            for c in 1:cols
                interpolated_points[r - 1, c] = interpolated_point[c]
            end
        end
        return bezier(t, interpolated_points)
    end
end

end # module


text:
#=
    Computes and array of magnitudes for 2D or 3D vector fields. Each component of
    the vectors field is provided in a separate array:
        magnitudemap( U, V ) or magnitudemap( U, V, W )
=#
function magnitudemap( arrays::Array{T,N}... ; typ=Float32 ) where {T,N}
    magnitudes = zeros( typ, size( arrays[1] ) )
    magnitudemap!( magnitudes, arrays... )
    return magnitudes
end

function magnitudemap!( magnitudes::Array{U,N}, arrays::Array{T,N}... ) where {T,U,N}
    numcomponents = length( arrays );
    @inbounds for e in 1:length( arrays[1] )
        sum2 = 0.0;
        @simd for c in 1:numcomponents
            sum2 += arrays[c][e]*arrays[c][e]
        end
        magnitudes[e] = convert( U, sqrt( sum2 ) )
    end
end

#=
    Computes a boolean mask, where vectors whose magnitude is > n * std are set to
    true. Each component of the vectors field is provided in a separate array.
=#
function stdmask( n, arrays::Array{T,N}... ) where {T,N}
    mask = zeros( Bool, size( arrays[1] ) )
    stdmask!( n, mask, arrays... );
    return mask
end

function stdmask!( n, mask::Array{Bool,N}, arrays::Array{T,N}... ) where {T,N}
    magnitudes = magnitudemap( arrays... );
    mean = Statistics.mean( magnitudes );
    std  = Statistics.std( magnitudes, mean=mean ) ;
    @simd for e in 1:length( magnitudes )
        @inbounds mask[ e ] = abs( magnitudes[e] - mean ) > n*std
    end
end

#=
    Replacement functions accept a coordinate and a radius argument ( a part from
    the vector field components ). zeroRepalce does not use the coordinate or
    radius, but mean and median require them.
=#
function zeroReplace( coords::NTuple{N,Int64}, radius::Int64,
                      arrays::Array{T,N}... ) where {T,N}
    return zeros( T, length(arrays) )
end

function medianReplace( coords::NTuple{N,Int64}, radius::Int64,
                        arrays::Array{T,N}... ) where {T,N}
    mincoords = max.( 1, coords .- radius );
    maxcoords = min.( size(arrays[1]), coords .+ radius );
    medians   = zeros( T, length( arrays ) );

    for c in 1:length( arrays )
        medians[ c ] = Statistics.median( arrays[c][ UnitRange.( mincoords, maxcoords )... ] );
    end
    return medians
end

function meanReplace( coords::NTuple{N,Int64}, radius::Int64,
                      arrays::Array{T,N}... ) where {T,N}
    mincoords = max.( 1, coords .- radius );
    maxcoords = min.( size(arrays[1]), coords .+ radius );
    means     = zeros( T, length( arrays ) );

    for c in 1:length( arrays )
        means[ c ] = Statistics.mean( arrays[c][ UnitRange.( mincoords, maxcoords )... ] )
    end
    return means
end

"""
    Replaces vectors whose magnitude is > n * std. The replacement function can be
    specified, as well as the radius of the neighbouring vectors that are used to
    compute median and mean vectors.

    Implemented replacement functions are.
        PIV3D.zeroReplace ----> replaces outlier vectors by a vector of zeros
        PIV3D.medianReplace --> replaces outlider vector by the median vector
        PIV3D.meanReplace ----> replaces outlier vectors by the mean vector
"""
function stdFilter( n, arrays::Array{T,N}...;
                    replace::Function=zeroReplace, radius=0 ) where {T,N}

    filtered = [ copy( arrays[idx] ) for idx in 1:length(arrays) ];

    if n > 0
        mask   = stdmask( n, arrays... );
        coords = CartesianIndices( mask );

        for e in 1:length(mask)
            if mask[e]
                replacement = replace( Tuple(coords[e]), radius, arrays... )
                for c in 1:length(arrays)
                    filtered[c][e] = replacement[c]
                end
            end
        end
    end

    return filtered
end

"""
    Replaces vectors whose magnitude is > n * std. The replacement function can be
    specified, as well as the radius of the neighbouring vectors that are used to
    compute median and mean vectors.

    Implemented replacement functions are.
        PIV3D.zeroReplace ----> replaces outlier vectors by a vector of zeros
        PIV3D.medianReplace --> replaces outlider vector by the median vector
        PIV3D.meanReplace ----> replaces outlier vectors by the mean vector
"""
function magnitudeFilter( threshold, u::Array{T,N}, v::Array{T,N} ) where {T<:Real,N}
    return magnitudeFilter( >, threshold, u, v )
end

function magnitudeFilter( f, threshold, arrays::Array{T,N}...;
                          replace::Function=zeroReplace, radius=0 ) where {T,N}

    magnitudes = magnitudemap( arrays... );
    coords     = CartesianIndices( magnitudes );
    filtered   = [ copy( arrays[idx] ) for idx in 1:length(arrays) ];

    for e in 1:length(magnitudes)
        if f( magnitudes[e], threshold )
            replacement = replace( Tuple(coords[e]), radius, arrays... )
            for c in 1:length(arrays)
                filtered[c][e] = replacement[c]
            end
        end
    end
    return filtered
end

""" Spatial averaging for 2D vector fields"""
function spaceAveraging( avg_radius::Int64, u::Array{T,2}, v::Array{T,2};
                         th=0.1 ) where {T<:Real}
    u_avg = zeros( T, size(u) );
    v_avg = zeros( T, size(v) );
    h, w  = size( u )

    for col in 1:w
        cmin = max( 1, col-avg_radius );
        cmax = min( w, col+avg_radius );

        for row in 1:h

            mag = sqrt( u[ row, col ]^2 + v[ row, col ]^2 );
            if mag < th
                continue
            end

            rmin = max( 1, row-avg_radius );
            rmax = min( h, row+avg_radius );
            N = length(cmin:cmax)*length(rmin:rmax);

            mean_u = 0.0;
            mean_v = 0.0;
            @inbounds for x in cmin:cmax
                @simd for y in rmin:rmax
                    mean_u += u[y,x]
                    mean_v += v[y,x]
                end
            end

            u_avg[ row, col ] = mean_u/N;
            v_avg[ row, col ] = mean_v/N;
        end
    end
    return u_avg, v_avg
end

""" Spatial averaging for 3D vector fields"""
function spaceAveraging( avg_radius::Int64, U::Array{T,3}, V::Array{T,3}, W::Array{T,3} ) where {T<:Real}

    u_avg   = zeros( T, size(U) );
    v_avg   = zeros( T, size(V) );
    w_avg   = zeros( T, size(W) );
    h, w, d = size( U )

    for zet in 1:d
        zmin = max( 1, zet-avg_radius );
        zmax = min( d, zet+avg_radius );

        for col in 1:w
            cmin = max( 1, col-avg_radius );
            cmax = min( w, col+avg_radius );

            for row in 1:h

                mag = sqrt( U[row,col,zet]^2 + V[row,col,zet]^2 + W[row,col,zet]^2 )
                if mag == 0
                    continue
                end

                rmin = max( 1, row-avg_radius );
                rmax = min( h, row+avg_radius );
                N = length(rmin:rmax)*length(cmin:cmax)*length(zmin:zmax);

                mean_u = 0.0;
                mean_v = 0.0;
                mean_w = 0.0;

                # Looping over neighbours
                for z in zmin:zmax
                    for x in cmin:cmax
                        for y in rmin:rmax
                            mean_u += U[y,x,z]
                            mean_v += V[y,x,z]
                            mean_w += W[y,x,z]
                        end
                    end
                end

                u_avg[ row, col, zet ] = mean_u/N;
                v_avg[ row, col, zet ] = mean_v/N;
                w_avg[ row, col, zet ] = mean_w/N;
            end
        end
    end

    return u_avg, v_avg, w_avg
end

""" space-time averaging """ # fancier circular mask in notes_16Oct20
function spaceTimeAveraging( avg_rs::Int64, avg_rt::Int64, u::Array{T,3}, v::Array{T,3}; th=0.0 )  where {T<:Real}

    u_avg = zeros( T, size(u) );
    v_avg = zeros( T, size(v) );

    h, w, time = size( u )

    for t in 1:time
        tmin = max(   1 , t - avg_rt );
        tmax = min( time, t + avg_rt );

        for c in 1:w
            cmin = max( 1, c - avg_rs );
            cmax = min( w, c + avg_rs );

            for r in 1:h

                magnitude = sqrt( u[r,c,t]^2 + v[r,c,t]^2 )
                if  magnitude < th
                    continue
                end

                rmin = max( 1, r - avg_rs );
                rmax = min( h, r + avg_rs );
                N = length(tmin:tmax) + length(rmin:rmax)*length(cmin:cmax) - 1

                # central vector is counted twice, so we account for it by substracting it once
                mean_u = -u[r,c,t];
                mean_v = -v[r,c,t];

                # spatial averaging
                for x in cmin:cmax
                    for y in rmin:rmax
                        mean_u += u[y,x,t]
                        mean_v += v[y,x,t]
                    end
                end

                # temporal averaging
                for tp in tmin:tmax
                    mean_u += u[r,c,tp]
                    mean_v += v[r,c,tp]
                end

                u_avg[r,c,t] = mean_u/N;
                v_avg[r,c,t] = mean_v/N;
            end
        end
    end

    return u_avg, v_avg
end

""" space-time averaging """ # fancier circular mask in notes_16Oct20
function spaceTimeAveraging( avg_rs::Int64, avg_rt::Int64, u::Array{T,4}, v::Array{T,4}, w::Array{T,4}; th=0.0 )  where {T<:Real}

    u_avg = zeros( T, size(u) );
    v_avg = zeros( T, size(v) );
    w_avg = zeros( T, size(w) );

    h, w, d, time = size( u )

    for t in 1:time
        tmin = max(   1 , t - avg_rt );
        tmax = min( time, t + avg_rt );

        for s in 1:d
            smin = max( 1, s - avg_rs );
            smax = min( d, s + avg_rs );

            for c in 1:w
                cmin = max( 1, c - avg_rs );
                cmax = min( w, c + avg_rs );

                for r in 1:h

                    magnitude = sqrt( u[r,c,s,t]^2 + v[r,c,s,t]^2 + w[r,c,s,t]^2 )
                    if  magnitude < th
                        continue
                    end

                    rmin = max( 1, r - avg_rs );
                    rmax = min( h, r + avg_rs );
                    N = length(tmin:tmax) + length(smin:smax)*length(cmin:cmax)*length(rmin:rmax)

                    # central vector is counted twice, so we account for it by substracting it once
                    mean_u = -u[ r, c, s, t ];
                    mean_v = -v[ r, c, s, t ];
                    mean_w = -w[ r, c, s, t ];

                    # spatial averaging, t is constant
                    for z in smin:smax
                        for x in cmin:cmax
                            for y in rmin:rmax
                                mean_u += u[y,x,z,t]
                                mean_v += v[y,x,z,t]
                                mean_w += w[y,x,z,t]
                            end
                        end
                    end

                    # temporal averaging, (r,c,s) are constant
                    for tp in tmin:tmax
                        mean_u += u[r,c,s,tp]
                        mean_v += v[r,c,s,tp]
                        mean_w += w[r,c,s,tp]
                    end

                    u_avg[ r,c,s,t ] = mean_u/N;
                    v_avg[ r,c,s,t ] = mean_v/N;
                    w_avg[ r,c,s,t ] = mean_w/N;
                end
            end
        end
    end

    return u_avg, v_avg, w_avg
end

""" Spatial averaging + similarity thresholding combo """
function similarityAveraging3D( avg_radius::Int64, U::Array{T,3}, V::Array{T,3}, W::Array{T,3};
                                st=0.0 ) where {T<:Real}

    u_avg   = zeros( T, size(U) );
    v_avg   = zeros( T, size(V) );
    w_avg   = zeros( T, size(W) );
    h, w, d = size( U )

    for zet in 1:d
        zmin = max( 1, zet-avg_radius );
        zmax = min( d, zet+avg_radius );

        for col in 1:w
            cmin = max( 1, col-avg_radius );
            cmax = min( w, col+avg_radius );

            for row in 1:h

                u1 = U[row,col,zet]
                v1 = V[row,col,zet]
                w1 = W[row,col,zet]
                mag = sqrt( u1*u1 + v1*v1 + w1*w1 )
                if mag == 0
                    continue
                end

                nu1 = U[row,col,zet]/mag
                nv1 = V[row,col,zet]/mag
                nw1 = W[row,col,zet]/mag

                rmin = max( 1, row-avg_radius );
                rmax = min( h, row+avg_radius );

                len  = length(rmin:rmax)*length(cmin:cmax)*length(zmin:zmax);

                # variables to hold the averaged vector component
                mean_u = 0.0;
                mean_v = 0.0;
                mean_w = 0.0;

                # counter of similar neighbouring vectors
                n  = 0;

                for z in zmin:zmax
                    for x in cmin:cmax
                        for y in rmin:rmax

                            mag2 = sqrt( U[y,x,z]*U[y,x,z] + V[y,x,z]*V[y,x,z] + W[y,x,z]*W[y,x,z] )
                            nu2  = U[y,x,z]/mag2
                            nv2  = V[y,x,z]/mag2
                            nw2  = W[y,x,z]/mag2
                            dot  = nu1*nu2 + nv1*nv2 + nw1*nw2

                            if dot > st
                                mean_u += U[y,x,z]
                                mean_v += V[y,x,z]
                                mean_w += W[y,x,z]
                                n += 1;
                            end
                        end
                    end
                end

                sim  = (n/len)^2;
                mmag = sqrt( mean_u*mean_u + mean_v*mean_v + mean_w*mean_w );

                u_avg[ row, col, zet ] = mean_u/mmag * sim;
                v_avg[ row, col, zet ] = mean_v/mmag * sim;
                w_avg[ row, col, zet ] = mean_w/mmag * sim;
            end
        end
    end

    return u_avg, v_avg, w_avg
end

""" Spatial averaging """
function similarityMap( avg_radius::Int64, U::Array{T,3}, V::Array{T,3}, W::Array{T,3};
                        st=0.0 ) where {T<:Real}

    similty = zeros( Float32, size(U) );
    h, w, d = size( U )

    for zet in 1:d
        zmin = max( 1, zet-avg_radius );
        zmax = min( d, zet+avg_radius );

        for col in 1:w
            cmin = max( 1, col-avg_radius );
            cmax = min( w, col+avg_radius );

            for row in 1:h

                mag = sqrt( U[row,col,zet]^2 + V[row,col,zet]^2 + W[row,col,zet]^2 )
                if mag == 0
                    continue
                end

                nu1 = U[row,col,zet]/mag
                nv1 = V[row,col,zet]/mag
                nw1 = W[row,col,zet]/mag

                rmin = max( 1, row-avg_radius );
                rmax = min( h, row+avg_radius );

                len  = length(rmin:rmax)*length(cmin:cmax)*length(zmin:zmax);

				len  = avg_radius*avg_radius*avg_radius;

                n = 0;
                for z in zmin:zmax
                    for x in cmin:cmax
                        for y in rmin:rmax
                            mag2 = sqrt( U[y,x,z]^2 + V[y,x,z]^2 + W[y,x,z]^2 )
                            nu2  = U[y,x,z]/mag2
                            nv2  = V[y,x,z]/mag2
                            nw2  = W[y,x,z]/mag2
                            dot  = nu1*nu2 + nv1*nv2 + nw1*nw2
                            if dot > st
                                n += 1;
                            end
                        end
                    end
                end

                similty[row,col,zet] = (n/len);
            end
        end
    end

    return similty
end

function velocityMap( U::Array{T,3}, V::Array{T,3}, W::Array{T,3} ) where {T<:Real}

    mags = zeros( Float32, size( U ) )
    for li in 1:length( U )
        mags[li] = sqrt( U[li]*U[li] + V[li]*V[li] + W[li]*W[li] )
    end
    return mags
end

# 2D
function sink( radius::Int64, pos::II )
    return [ [ (pos[1] - y )/sqrt( (pos[1]-y)^2 + (pos[2]-x)^2 ),
               (pos[2] - x )/sqrt( (pos[1]-y)^2 + (pos[2]-x)^2 ) ] for y in pos[1]-radius:pos[1]+radius,
                                                                       x in pos[2]-radius:pos[2]+radius ]
end

# 3D
function sink( radius::Int64, pos::III )
    return [ [ (pos[1] - y )/sqrt( (pos[1]-y)^2 + (pos[2]-x)^2 + (pos[3]-z)^2 ),
               (pos[2] - x )/sqrt( (pos[1]-y)^2 + (pos[2]-x)^2 + (pos[3]-z)^2 ),
               (pos[3] - z )/sqrt( (pos[1]-y)^2 + (pos[2]-x)^2 + (pos[3]-z)^2 ) ] for y in pos[1]-radius:pos[1]+radius,
                                                                                x in pos[2]-radius:pos[2]+radius,
                                                                                z in pos[3]-radius:pos[3]+radius ]
end

function crossCorrelateSink!( u, v, w, sink, corr )

    hi, wi, di = size(sink) # interrogation area
    hs, ws, ds = size(u)    # search area

    @inbounds begin
    for z in di:-1:1
        for c in wi:-1:1
            for r in hi:-1:1
                px = sink[ r, c, z ];
                iy = hi-r+1;
                ix = wi-c+1;
                iz = di-z+1;
                corr[ iy:(iy-1+hs),
                      ix:(ix-1+ws),
                      iz:(iz-1+ds) ] .+= px[1] .* u .+ px[2] .* v .+ px[3] .* w
            end
        end
    end
    end # @inbounds
end

function divergenceMap( r::Int64, U::Array{T,3}, V::Array{T,3}, W::Array{T,3}  ) where {T<:Real}

    interSink = sink( r, (0,0,0) );

    for idx in 1:length(interSink)
        isnan(interSink[idx][1]) && ( interSink[idx][1] = 0.0 );
        isnan(interSink[idx][2]) && ( interSink[idx][2] = 0.0 );
        isnan(interSink[idx][3]) && ( interSink[idx][3] = 0.0 );
    end

    # normalizing vector fields
    NU = zeros( T, size( U ) );
    NV = zeros( T, size( V ) );
    NW = zeros( T, size( W ) );

    for e in 1:length( U )
        mag   = sqrt( U[e]*U[e] + V[e]*V[e] + W[e]*W[e] )
        NU[e] = ( mag == 0 ) ? 0.0 : U[e]/mag
        NV[e] = ( mag == 0 ) ? 0.0 : V[e]/mag
        NW[e] = ( mag == 0 ) ? 0.0 : W[e]/mag
    end

    cmatrix = zeros( Float64, size(NU) .+ 2*r )

    crossCorrelateSink!( NU, NV, NW, interSink, cmatrix )

    return cmatrix
end

""" PIV Trajectories, work started by Michelle Gottlieb  """

function PIVtrajectories( U::Array{P,4}, V::Array{P,4}, W::Array{P,4},
                          T0, T1, numpoints; facF=1,
				          subr=( -1: -1, -1: -1, -1: -1 ) ) where {P<:Real}
    return PIVtrajectories( U, V, W, T0, T1, numpoints, 1:size(U,1), 1:size(U,2), 1:size(U,3), subr=subr, facF=facF )
end

function PIVtrajectories( U::Array{P,4}, V::Array{P,4}, W::Array{P,4},
                          T0, T1, numpoints, vrange, lrange, drange;
                          subr=( -1: -1, -1: -1, -1: -1 ), facF=1 ) where {P<:Real}

    # Vertical, lateral and depth axes.
    VV = collect( vrange );
    LL = collect( lrange );
    DD = collect( drange );

    dims = ( length(VV), length(LL), length(DD) );

    # I am still not sure of what these scaling factors do
    endV = VV[ end ]
    endH = LL[ end ]
    endD = DD[ end ]
    facV = 1 # endV/dims[1]
    facL = 1 # endH/dims[2]
    facD = 1 # endD/dims[3]

    numT = T1 - T0;
    TrajectoriesV = zeros( Float64, numT, numpoints )
    TrajectoriesL = zeros( Float64, numT, numpoints )
    TrajectoriesD = zeros( Float64, numT, numpoints )

    for pidx in 1:numpoints
        # Each iteration in this loop corresponds to a new particle
        # A particle is nothing but a pair of coordinates
        # Random coordinates of the particle inside VV, HH and DD
        c1 = rand( subr[1].start == -1 ? (2:dims[1]-1) : subr[1] )
        c2 = rand( subr[2].start == -1 ? (2:dims[2]-1) : subr[2] )
        c3 = rand( subr[3].start == -1 ? (2:dims[3]-1) : subr[3] )

        # Trajectory starts at the ( VV, HH ) position
        TrajectoriesV[ 1, pidx ] = VV[ c1 ] + 0.0
        TrajectoriesL[ 1, pidx ] = LL[ c2 ] + 0.0
        TrajectoriesD[ 1, pidx ] = DD[ c3 ] + 0.0

        # Translation converted into indices on VV and HH
        Ut = facF * U[ c1, c2, c3, T0 ]
        Vt = facF * V[ c1, c2, c3, T0 ]
        Wt = facF * W[ c1, c2, c3, T0 ]

        # Time-iteration
        for t in 2:numT

            TrajectoriesV[t,pidx] = TrajectoriesV[t-1, pidx] + Ut
            TrajectoriesL[t,pidx] = TrajectoriesL[t-1, pidx] + Vt
            TrajectoriesD[t,pidx] = TrajectoriesD[t-1, pidx] + Wt

            c1 = round( Int64, TrajectoriesV[t,pidx] / facV )
            c2 = round( Int64, TrajectoriesL[t,pidx] / facL )
            c3 = round( Int64, TrajectoriesD[t,pidx] / facD )

            # If new position is out of the coordinate space, stop
            if c1 > dims[1] || c1 < 1 ||
               c2 > dims[2] || c2 < 1 ||
               c3 > dims[3] || c3 < 1
                TrajectoriesV[ t:end, pidx ] .= TrajectoriesV[ t-1, pidx ]
                TrajectoriesL[ t:end, pidx ] .= TrajectoriesL[ t-1, pidx ]
                TrajectoriesD[ t:end, pidx ] .= TrajectoriesD[ t-1, pidx ]
                break
            end

            # Using nextV and nextH to sample the translation from U and V
            Ut = facF * U[ c1, c2, c3, T0+t-1 ]
            Vt = facF * V[ c1, c2, c3, T0+t-1 ]
            Wt = facF * W[ c1, c2, c3, T0+t-1 ]

        end # Time-iteration

    end

    return TrajectoriesV, TrajectoriesL, TrajectoriesD
end


text:
module test_millboard_regex

using Millboard
using Test

Millboard.set_table_mode(:markdown)

str = "The quick brown fox jumps over the lazy dog α,β,γ"
regexz = [r"b[\w]*n", r"[\w]{4,}"]
functions = [match]
board = [func(regex, str) for regex=regexz, func=functions]
@test sprint(show, table(board, colnames=functions, rownames=regexz)) == """
|             | match               |
|-------------|---------------------|
| r"b[\\w]*n"  | RegexMatch("brown") |
| r"[\\w]{4,}" | RegexMatch("quick") |"""

end # module test_millboard_regex


text:
# compare FastqRead from gz files and not gz files
function test_fastq_gz_reading()
    data_dir = joinpath(dirname(@__FILE__),"..","testdata","fq")
    # open the gz file
    gzstream = fastq_open(joinpath(data_dir, "R1.fq.gz"))
    # open the fq file, whose content is same as the gz file
    fqstream = fastq_open(joinpath(data_dir, "R1.fq"))

    while (fq = fastq_read(fqstream)) != false
        gzfq = fastq_read(gzstream)
        if fq.sequence != gzfq.sequence || fq.name != gzfq.name || fq.quality.qual != gzfq.quality.qual
            return false
        end
    end
    return true
end

function test_fastq_writing()
    data_dir = joinpath(dirname(@__FILE__),"..","testdata","fq")

    # open the fq file, whose content is same as the gz file
    instream = fastq_open(joinpath(data_dir, "R1.fq"))
    outstream = fastq_open(joinpath(data_dir, "R1.out.fq"), "w")

    while (fq = fastq_read(instream)) != false
        fastq_write(outstream, fq)
    end
    close(outstream)

    instream1 = fastq_open(joinpath(data_dir, "R1.fq"))
    instream2 = fastq_open(joinpath(data_dir, "R1.out.fq"))

    while (fq1 = fastq_read(instream1)) != false
        fq2 = fastq_read(instream2)
        if fq1.sequence != fq2.sequence || fq1.name != fq2.name || fq1.quality.qual != fq2.quality.qual
            return false
        end
    end

    rm(joinpath(data_dir, "R1.out.fq"))

    return true
end

@test test_fastq_gz_reading()
@test test_fastq_writing()

include("vcf.jl")


text:
# Julia1_3.2.jl  Looping Control Structures
#
# while loops keeps evaluating body of loop
#   until condition is false
# as with other blocks value defined in block
#    is not available in outer blocks
#    but is available to a nested inner blocks
i=0
while i<26
    global i
    print('A'+i)
    i+=1
end

# for loop
# i steps through <range> 0 to 25 by 1
# loop index not available outside loop block
for i=0:25
    print('a'+i)
end

# i steps through range by 2s*
for i=0:2:25
    print('a'+i)
end

# for can loop over any container
for i in [2,0,19]
    print('a'+i)
end

# The word <break> will exit a <while> or <for> loop
# This is useful for iterations that can terminate
#   early either because it is diverging or
#   has converged below an error bound.
i=0
while i<26
    global i
    print('A'+i)
    i+=1
    if i>5; break; end
end

# <continue> key word skips to next loop start
# in a for loop
for i=0:25
    2<i<23 && continue
    print('a'+i)
end

###

# nested <for> loops for matrices and in general
M=Array{Int,2}(undef,3,3)
(d1,d2)=size(M)
for i=1:d1, j=1:d2
    global M
    M[i,j]=i*j
end
println(M)

# exception handling
# try <do something>
# catch <if error, handle it>
# end
function fsqrt(x)
    try; sqrt(x)
    catch; sqrt(Complex(x,0))
    end
end
fsqrt(-4)

# <throw> your own error
# there is a list of 35 Built-in exceptions
#  but you can add your own if you don't find
#  it on the list
f(x)= x≥0 ? exp(-x) : throw(DomainError(x, "arg must be non-negative"))
#
f(-1)
f(1)


text:
@testset "proteins.jl" begin
    @testset "fasta_records" begin
        fasta_file = IOBuffer(
        """
        >seqA some description
        QIKDLLVSSSTDLDTTLKMK
        ILELPFASGDLSM
        >seqB
        VLMALGMTDLFIPSANLTG*
        """)
        records = fasta_records(fasta_file)
        @test length(records) == 2
        rec = records[1]
        @test FASTA.identifier(rec) == "seqA"
        @test FASTA.description(rec) == "some description"
        @test FASTA.sequence(rec) == AminoAcidSequence("QIKDLLVSSSTDLDTTLKMKILELPFASGDLSM")
        rec = records[2]
        @test FASTA.identifier(rec) == "seqB"
        @test FASTA.sequence(rec) == AminoAcidSequence("VLMALGMTDLFIPSANLTG*")
    end
    @testset "signature" begin
        record = FASTA.Record(
        """
        >seqA some description
        QIKDLLVSSSTDLDTTLKMK
        ILELPFASGDLSM
        """)
        sig = signature(record, 5, 2)
        @test signature(sig) == UInt32[0x05094006, 0x0afea36a]
    end
    @testset "fasta_signatures" begin
        fasta_file = IOBuffer(
        """
        >seqA some description
        QIKDLLVSSSTDLDTTLKMK
        ILELPFASGDLSM
        >seqB
        VLMALGMTDLFIPSANLTG*
        """)
        sigs = fasta_signatures(fasta_file, 5, 2)
        sig = sigs[1]
        identifier(sig) == "seqA"
        signature(sig) == UInt32[0x05094006, 0x0afea36a]
        sig = sigs[2]
        identifier(sig) == "seqB"
        signature(sig) == UInt32[0x05094006, 0x0afea36a]
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)

function spmatmul_witheps_and_flops{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}, eps;
                         sortindices::Symbol = :sortcols)
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    flops = 0

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Array(Ti, nB+1)
    rowvalC = Array(Ti, nnzC)
    nzvalC = Array(Tv, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end

                    flops += 2
                end
            end
            idx = colptrC[i]
            for vp in colptrC[i]:(ip - 1)
                col = rowvalC[vp]
                if col == i || abs(x[col]) > eps
                  nzvalC[idx] = x[col]
                  rowvalC[idx] = col
                  idx += 1
                end
            end
            ip = idx
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = Base.SparseMatrix.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C, flops
end

function spmatmul_witheps{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}, eps;
                         sortindices::Symbol = :sortcols)
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Array(Ti, nB+1)
    rowvalC = Array(Ti, nnzC)
    nzvalC = Array(Tv, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
            idx = colptrC[i]
            for vp in colptrC[i]:(ip - 1)
                col = rowvalC[vp]
                if col == i || abs(x[col]) > eps
                  nzvalC[idx] = x[col]
                  rowvalC[idx] = col
                  idx += 1
                end
            end
            ip = idx
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = Base.SparseMatrix.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C
end

@doc """ Using Gershgorin disc to estimate the bounds of eigen value. """
function gershgorin(A :: SparseMatrixCSC)
  hsize = size(A, 1)
  eMin = 10000
  eMax = -10000

  for i=1:hsize
    sumM = 0
    sumP = 0
    for j=A.colptr[i]:A.colptr[i+1]-1
      hx = abs(A.nzval[j])
      sumM += hx
      if A.rowval[j] == i
        sumP = A.nzval[j]
        sumM -= hx
      end
    end
    eMax = max(eMax, sumP + sumM)
    eMin = min(eMin, sumP - sumM)
  end

  return eMax, eMin
end

function count_CoSP2_flop(X)
  m = size(X, 1)
  occ = m/2
  idemTol = 1e-14
  eps = 1e-5

  eMax, eMin = gershgorin(X')
  X = SparseMatrixCSC{Float64, Int32}((eMax*speye(m) - X)/(eMax - eMin))
  assert(issym(X))
  idempErr = 0
  idempErr1 = 0
  idempErr2 = 0
  BreakLoop = 0
  iter = 0
  trXOLD = 0

  flops = 0
  spgemm_flops = 0

  while BreakLoop == 0
    trX = trace(X)
    flops += m

    X2, f = spmatmul_witheps_and_flops(X, X', eps) # CoSP2 does approximate spgemm
    flops += f
    spgemm_flops += f

    trX2 = trace(X2)
    flops += m
    trXOLD = trX

    limDiff = abs(trX2 - occ) - abs(2*trX - trX2 - occ)
    if limDiff > idemTol

      X = 2*X - X2
      flops += nnz(X) + nnz(X2)

      trX = 2*trX - trX2
    elseif limDiff < -idemTol
      X = X2 # output X is dead here, so we can free
      trX = trX2
    else
      trX = trXOLD
      BreakLoop = 1
    end

    idempErr2 = idempErr1
    idempErr1 = idempErr
    idempErr = abs(trX - trXOLD)
    #println(idempErr)

    iter += 1

    if iter >= 25 && idempErr >= idempErr2
      BreakLoop = 1
    end
  end

  X = 2*X
  flops += nnz(X)

  println("X sum = $(sum(X)), max = $(maximum(X))")
  println("Number of iterations = $iter")
  #println("$flops $spgemm_flops")
  return flops
end

function CoSP2_ref(X, flops)
  set_matrix_property(Dict(
      :X => SA_SYMM_VALUED,
    )
  )

  m = size(X, 1)
  occ = m/2
  idemTol = 1e-14
  eps = 1e-5

  eMax, eMin = gershgorin(X')
  X = SparseMatrixCSC{Float64, Int32}((eMax*speye(m) - X)/(eMax - eMin))
  assert(issym(X))
  idempErr = 0
  idempErr1 = 0
  idempErr2 = 0
  BreakLoop = 0
  iter = 0
  trXOLD = 0

  total_time = -time()
  spgemm_time = 0
  spadd_time = 0

  while BreakLoop == 0
    trX = trace(X)

    spgemm_time -= time()
    #X2 = X*X'
    X2 = spmatmul_witheps(X, X', eps) # CoSP2 does approximate spgemm
    spgemm_time += time()

    trX2 = trace(X2)
    trXOLD = trX

    limDiff = abs(trX2 - occ) - abs(2*trX - trX2 - occ)
    if limDiff > idemTol

      spadd_time -= time()
      X = 2*X - X2
      spadd_time += time()

      trX = 2*trX - trX2
    elseif limDiff < -idemTol
      X = X2 # output X is dead here, so we can free
      trX = trX2
    else
      trX = trXOLD
      BreakLoop = 1
    end

    idempErr2 = idempErr1
    idempErr1 = idempErr
    idempErr = abs(trX - trXOLD)
    #println(idempErr)

    iter += 1

    if iter >= 25 && idempErr >= idempErr2
      BreakLoop = 1
    end
  end

  X = 2*X

  total_time += time()

  println("X sum = $(sum(X)), max = $(maximum(X))")
  println("Number of iterations = $iter")
  println("Total time $total_time ($(flops/total_time/1e9) gflops), SpGEMM time $spgemm_time, SpADD time $spadd_time")
end

function CoSP2_call_replacement(X, flops)
  m = size(X, 1)
  occ = m/2
  idemTol = 1e-14
  eps = 1e-5

  eMax, eMin = gershgorin(X')
  X = SparseMatrixCSC{Float64, Int32}((eMax*speye(m) - X)/(eMax - eMin))
  assert(issym(X))
  idempErr = 0
  idempErr1 = 0
  idempErr2 = 0
  BreakLoop = 0
  iter = 0
  trXOLD = 0

  total_time = -time()
  spgemm_time = 0
  spadd_time = 0

  while BreakLoop == 0
    #trX = trace(X)
    trX = Sparso.trace(X)

    spgemm_time -= time()
    #X2 = X*X'
    #X2 = spmatmul_witheps(X, X', eps) # CoSP2 does approximate spgemm
    X2 = Sparso.SpSquareWithEps(X, eps)
    spgemm_time += time()

    # trX2 = trace(X2)
    trX2 = Sparso.trace(X2)
    trXOLD = trX

    limDiff = abs(trX2 - occ) - abs(2*trX - trX2 - occ)
    if limDiff > idemTol

      spadd_time -= time()
      #X = 2*X - X2 # input X is dead here, so we can free
      X = Sparso.SpAdd(2, X, -1, X2)
      spadd_time += time()

      trX = 2*trX - trX2
    elseif limDiff < -idemTol
      X = X2 # output X is dead here, so we can free
      trX = trX2
    else
      trX = trXOLD
      BreakLoop = 1
    end

    idempErr2 = idempErr1
    idempErr1 = idempErr
    idempErr = abs(trX - trXOLD)
    #println(idempErr)

    iter += 1

    if iter >= 25 && idempErr >= idempErr2
      BreakLoop = 1
    end
  end

  X = 2*X

  total_time += time()

  println("X sum = $(sum(X)), max = $(maximum(X))")
  println("Number of iterations = $iter")
  println("Total time $total_time ($(flops/total_time/1e9) gflops), SpGEMM time $spgemm_time, SpADD time $spadd_time")
end

function CoSP2_call_replacement_and_context_opt(X, flops)
  m = size(X, 1)
  occ = m/2
  idemTol = 1e-14
  eps = 1e-5

  eMax, eMin = gershgorin(X')
  X = SparseMatrixCSC{Float64, Int32}((eMax*speye(m) - X)/(eMax - eMin))
  assert(issym(X))
  idempErr = 0
  idempErr1 = 0
  idempErr2 = 0
  BreakLoop = 0
  iter = 0
  trXOLD = 0

  total_time = -time()
  spgemm_time = 0
  spadd_time = 0

  mknobX = (Sparso.new_matrix_knob)(:X, false, false, true, true, false, false) # X is symmetric

  fknob_spgemm = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(mknobX, fknob_spgemm)

  while BreakLoop == 0
    #trX = trace(X)
    trX = Sparso.trace(X)

    spgemm_time -= time()
    #X2 = X*X'
    #X2 = spmatmul_witheps(X, X', eps) # CoSP2 does approximate spgemm
    X2 = Sparso.SpSquareWithEps(X, eps, fknob_spgemm)
    spgemm_time += time()

    # trX2 = trace(X2)
    trX2 = Sparso.trace(X2)
    trXOLD = trX

    limDiff = abs(trX2 - occ) - abs(2*trX - trX2 - occ)
    if limDiff > idemTol

      spadd_time -= time()
      #X = 2*X - X2 # input X is dead here, so we can free
      X = Sparso.SpAdd(2, X, -1, X2)
      spadd_time += time()

      trX = 2*trX - trX2
    elseif limDiff < -idemTol
      X = X2 # output X is dead here, so we can free
      trX = trX2
    else
      trX = trXOLD
      BreakLoop = 1
    end

    idempErr2 = idempErr1
    idempErr1 = idempErr
    idempErr = abs(trX - trXOLD)
    #println(idempErr)

    iter += 1

    if iter >= 25 && idempErr >= idempErr2
      BreakLoop = 1
    end
  end

  X = 2*X

  total_time += time()

  println("X sum = $(sum(X)), max = $(maximum(X))")
  println("Number of iterations = $iter")
  println("Total time $total_time ($(flops/total_time/1e9) gflops), SpGEMM time $spgemm_time, SpADD time $spadd_time")
end

file_name = "hmatrix.512.mtx"
if length(ARGS) >= 1
  file_name = ARGS[1]
end

X = matrix_market_read(file_name)
assert(issym(X))

flops = count_CoSP2_flop(X)

# currently, something wrong with spmatmul_witheps used in CoSP2_ref.
# So, ignore the results of CoSP2_ref
println("\nOriginal:")
CoSP2_ref(X, flops)
CoSP2_ref(X, flops)
println("End original.")

# Expected results for hmatrix.512.mtx: X sum = 6106.4049873666145, max = 1.280873477930676
# Number of iterations = 25
# Expected results for hmatrix.1024.mtx: X sum = 12212.785128790038, max = 1.2808837088549991
# Number of iterations = 25
#println("\nCoSP2_call_replacement:")
#CoSP2_call_replacement(X, flops)
#CoSP2_call_replacement(X, flops)
#println("End CoSP2_call_replacement.")

#println("\nCoSP2_call_replacement_and_context_opt:")
#CoSP2_call_replacement_and_context_opt(X, flops)
#CoSP2_call_replacement_and_context_opt(X, flops)
#println("End CoSP2_call_replacement_and_context_opt.")

println("\nAccelerated:")
@acc CoSP2_ref(X, flops)
@acc CoSP2_ref(X, flops)
println("End accelerated.")

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
function orthogonal_rows(axis::AbstractVector)
    if norm(axis) > 0
        axis = normalize(axis)
    end
    A = svd(skew(axis)).Vt
    inds = SA[1; 2; 3]
    V1 = A[1,inds]'
    V2 = A[2,inds]'
    V3 = axis' # instead of A[3,:] for correct sign: abs(axis) = abs(A[3,:])

    return V1, V2, V3
end

function orthogonal_columns(axis::AbstractVector)
    V1, V2, V3 = orthogonal_rows(axis)
    return V1', V2', V3'
end


text:
@testset "Rosenbrock" begin
    using LinearAlgebra: Symmetric

    nls = Rosenbrock()

    @test parameter_size(nls) == 2
    @test residue_size(nls) == 2

    θ=Float64[1;0]

    (r, J) = eval_r_J(nls,θ)

    @test r ≈ Float64[0; -10]
    @test J ≈ Float64[-1 0; -20 10]
    @test eval_nls_fobj(r) ≈ 50

    grad = eval_nls_∇fobj(r,J)
    @test grad ≈ [200.0; -100.0]

    H = eval_nls_∇∇fobj(J)
    @test H ≈ Float64[401.0 -200.0; -200.0 100.0]
end


text:
export Target, targettrain!, difftargettrain!;

mutable struct Target{F, S, L}
	f::F
	dual_f::S
	loss::L
	σ::Real
	in::Array
	out::Union{Array, TrackedArray}
	debuglog::Dict{String, Array}
end

treelike(Target);

function Target(f, dual_f, loss; σ::Real = 0.2)::Target
	return Target(f, dual_f, loss, σ, Array{Float32, 0}(), TrackedArray(Array{Float32, 0}()), Dict{String, Array}());
end

function (a::Target)(x)
	a.in = data(x);
	a.out = a.f(a.in);
	return a.out;
end

function Base.show(io::IO, a::Target)
	print(io, "Target(");
	print(io, a.f);
	print(io, ", ");
	print(io, a.dual_f);
	print(io, ")");
end

# targetprop

function targetprop!(a::Target, targetTuple; debug::Array = [])
	(target, realgrad) = targetTuple;
	function debuglog(name, value)
		if name in debug
			if !haskey(a.debuglog, name)
				a.debuglog[name] = Array{typeof(data(value)), 1}();
			end
			push!(a.debuglog[name], data(value));
		end
	end

	if "difference" in debug
		inc = similar(a.in, (size(a.in, 1), size(a.in, 2) - 1));
		for i in 1:size(inc,2)
			inc[:,i] = a.in[:, i + 1] - a.in[:, i];
		end
		ind = mapreduce(i->inc[:, i], +, 1:size(inc,2))
		debuglog("difference", norm(a.f(ind)));
	end
	if "average" in debug
		avg = mean(map(i->a.in[:, i], 1:size(a.in, 2)));
		debuglog("average", mean(map(i->norm(data(a.f(a.in[:, i] - avg))), 1:size(a.in, 2))));
	end

	retgrad = [];
	if "angle" in debug
		paramin = param(a.in);
		fcopy = deepcopy(a.f);
		back!(fcopy(paramin), realgrad);
		retgrad = paramin.grad;
	end
	if isa(a.out, TrackedArray)
		l1 = a.loss(target, a.out); # TODO: Regularisation
		debuglog("Layer-local loss function", l1);
		back!(l1);
	end
	ϵ = a.σ * randn(size(a.in));
	l2 = a.loss(a.dual_f(data(a.f(a.in .+ ϵ))), a.in .+ ϵ);
	#l2 = a.loss(a.dual_f(a.f(a.in .+ ϵ)), a.in .+ ϵ); # Non-standard approach
	debuglog("Dual layer-local loss function", l2);
	if "Reverse auto-encoder loss" in debug
		l2i = a.loss(data(a.f(data(a.dual_f(data(a.out))))), data(a.out));
		debuglog("Reverse auto-encoder loss", l2i);
	end
	back!(l2);
	if isa(a.out,TrackedArray)
		veclength(a) = sqrt(sum(a.^2));
		vecangle(a, b) = begin
			acosd(min(dot(a, b) / (veclength(a) * veclength(b)), 1.0))
		end
		if "angle" in debug
			debuglog("angle", vecangle(vcat(map(i->vec(i.grad), params(a.f))...), vcat(map(i->vec(i.grad), params(fcopy))...)));
		end
		debuglog("jacobian", jacobian(a.f, a.in[:, rand(1:size(a.in, 2))]));
	end
	return (data(a.dual_f(data(target))), retgrad);
end

function targetprop!(a::Chain, target; debug::Array = [])
	foldl((m, x) -> targetprop!(x, m; debug = debug), target, reverse(a.layers))
	return target;
end

function targettrain!(model, modelloss, data, opt; η::Real = 0.001, cb = () -> (), debug::Array = [])
	cb = Optimise.runall(cb);
	opt = Optimise.runall(opt);
	@progress for d in data
		y_hat = Flux.data(model(d[1]));
		grad = param(y_hat);
		back!(modelloss(grad, d[2]));
		target = @fix y_hat - η * length(d[2]) * grad.grad;
		Optimise.@interrupts targetprop!(model, (target, grad.grad); debug = debug);
		opt();
		cb() == :stop && break;
	end
end

# difftargetprop

function difftargetprop!(a::Target, packedTarget; debug::Array = [])
	(target, realgrad, last) = packedTarget
	if !last
		target += data(a.out);
	end
	(nextTarget, nextGrad) = targetprop!(a, (target, realgrad); debug = debug);
	nextTarget -= data(a.dual_f(data(a.out)));
	return (nextTarget, nextGrad, false);
end

function difftargetprop!(a::Chain, target; debug::Array = [])
	foldl((m, x) -> difftargetprop!(x, m; debug = debug), target, reverse(a.layers))
	return target;
end

function difftargettrain!(model, modelloss, data, opt; η::Real = 0.001, cb = () -> (), debug::Array = [])
	cb = Optimise.runall(cb);
	opt = Optimise.runall(opt);
	@progress for d in data
		y_hat = Flux.data(model(d[1]));
		grad = param(y_hat);
		back!(modelloss(grad, d[2]));
		target = @fix y_hat - η * length(d[2]) * grad.grad;
		Optimise.@interrupts difftargetprop!(model, (target, grad.grad, true); debug = debug);
		opt();
		cb() == :stop && break;
	end
end


text:
using Test

using BioGenerics
using FormatSpecimens
using GenomicFeatures
using XAM

import BioAlignments: Alignment, AlignmentAnchor, OP_START, OP_MATCH, OP_DELETE
import BGZFStreams: BGZFStream
import BioGenerics.Exceptions: MissingFieldException
import BioSequences: @dna_str, @aa_str

# Generate a random range within `range`.
function randrange(range)
    x = rand(range)
    y = rand(range)
    if x < y
        return x:y
    else
        return y:x
    end
end

include("test_sam.jl")
include("test_bam.jl")
include("test_crosscheck.jl")


text:
int_rules_4_3_3_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.3.1*(a+b*tan)^m*(c+d*tan)^n*(A+B*tan) =#
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~a * ~c) / ~f) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m - 1) * (~c + ~d * ~x) ^ (~n - 1) * (~A + ~B * ~x), ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 + (~b) ^ 2, 0))
    @apply_utils Antiderivative(((~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => ((~B * ~d) / ~b) * Antiderivative(tan(~e + ~f * ~x), ~x) + (1 / ~b) * Antiderivative(Simp(~A * ~b * ~c + (~A * ~b * ~d + ~B * (~b * ~c - ~a * ~d)) * tan(~e + ~f * ~x), ~x) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (-((~A * ~b - ~a * ~B)) * (~a * ~c + ~b * ~d) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m) / (2 * (~a) ^ 2 * ~f * ~m) + (1 / (2 * ~a * ~b)) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * Simp(~A * ~b * ~c + ~a * ~B * ~c + ~a * ~A * ~d + ~b * ~B * ~d + 2 * ~a * ~B * ~d * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~m, -1) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~b * ~c - ~a * ~d) * (~A * ~b - ~a * ~B) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * Simp(((~a * ~A * ~c + ~b * ~B * ~c + ~A * ~b * ~d) - ~a * ~B * ~d) - (((~A * ~b * ~c - ~a * ~B * ~c) - ~a * ~A * ~d) - ~b * ~B * ~d) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~m, -1) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~B * ~d * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1)) + Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * Simp((~A * ~c - ~B * ~d) + (~B * ~c + ~A * ~d) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && Not(LeQ(~m, -1)))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (-((~a) ^ 2) * (~B * ~c - ~A * ~d) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~b * ~c + ~a * ~d) * (~n + 1)) - (~a / (~d * (~b * ~c + ~a * ~d) * (~n + 1))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp((~A * ~b * ~d * ((~m - ~n) - 2) - ~B * (~b * ~c * (~m - 1) + ~a * ~d * (~n + 1))) + (~a * ~A * ~d * (~m + ~n) - ~B * (~a * ~c * (~m - 1) + ~b * ~d * (~n + 1))) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (GtQ(~m, 1) && LtQ(~n, -1))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~b * ~B * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n)) + (1 / (~d * (~m + ~n))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp((~a * ~A * ~d * (~m + ~n) + ~B * (~a * ~c * (~m - 1) - ~b * ~d * (~n + 1))) - (~B * (~b * ~c - ~a * ~d) * (~m - 1) - ~d * (~A * ~b + ~a * ~B) * (~m + ~n)) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (GtQ(~m, 1) && Not(LtQ(~n, -1)))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (-((~A * ~b - ~a * ~B)) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n) / (2 * ~a * ~f * ~m) + (1 / (2 * (~a) ^ 2 * ~m)) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n - 1) * Simp((~A * (~a * ~c * ~m + ~b * ~d * ~n) - ~B * (~b * ~c * ~m + ~a * ~d * ~n)) - ~d * (~b * ~B * (~m - ~n) - ~a * ~A * (~m + ~n)) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~m, 0) && GtQ(~n, 0))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~a * ~A + ~b * ~B) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (2 * ~f * ~m * (~b * ~c - ~a * ~d)) + (1 / (2 * ~a * ~m * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(~A * (~b * ~c * ~m - ~a * ~d * (2 * ~m + ~n + 1)) + ~B * (~a * ~c * ~m - ~b * ~d * (~n + 1)) + ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 1) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~m, 0) && Not(GtQ(~n, 0)))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~B * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n) / (~f * (~m + ~n)) + (1 / (~a * (~m + ~n))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n - 1) * Simp((~a * ~A * ~c * (~m + ~n) - ~B * (~b * ~c * ~m + ~a * ~d * ~n)) + (~a * ~A * ~d * (~m + ~n) - ~B * (~b * ~d * ~m - ~a * ~c * ~n)) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~n, 0)))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A * ~d - ~B * ~c) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~a * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp((~A * (~b * ~d * ~m - ~a * ~c * (~n + 1)) - ~B * (~b * ~c * ~m + ~a * ~d * (~n + 1))) - ~a * (~B * ~c - ~A * ~d) * (~m + ~n + 1) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~b * ~B) / ~f) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m - 1) * (~c + ~d * ~x) ^ ~n, ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && EqQ(~A * ~b + ~a * ~B, 0)))
    @apply_utils Antiderivative(((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b + ~a * ~B) / (~b * ~c + ~a * ~d)) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m, ~x) - ((~B * ~c - ~A * ~d) / (~b * ~c + ~a * ~d)) * Antiderivative(((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~a - ~b * tan(~e + ~f * ~x))) / (~c + ~d * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ(~A * ~b + ~a * ~B, 0)))
    #= Antiderivative(((~a)+(~b')*tan((~e')+(~f')*(~x)))^(~m)*((~c')+(~d')*tan((~e')+(~f')*(~x)))^(~n)*((~A')+ (~B')*tan((~e')+(~f')*(~x))),~x) :=  (A*b-a*B)/b*Antiderivative((a+b*tan(e+f*x))^m*(c+d*tan(e+f*x))^n,x) + B/b*Antiderivative((a+b*tan(e+f*x))^(m+1)*(c+d*tan(e+f*x))^n,x) <--  FreeQ([a,b,c,d,e,f,A,B,m],x) && NeQ(b*c-a*d,0) && EqQ(a^2+b^2,0) &&  NeQ(c^2+d^2,0) =#
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b + ~a * ~B) / ~b) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n, ~x) - (~B / ~b) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (~a - ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ(~A * ~b + ~a * ~B, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~A + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A) ^ 2 / ~f) * Subst(Antiderivative(((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n) / (~A - ~B * ~x), ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (Not(IntegersQ(2 * ~m, 2 * ~n)) && EqQ((~A) ^ 2 + (~B) ^ 2, 0))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A + im * ~B) / 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 - im * tan(~e + ~f * ~x)), ~x) + ((~A - im * ~B) / 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 + im * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && (Not(IntegersQ(2 * ~m, 2 * ~n)) && NeQ((~A) ^ 2 + (~B) ^ 2, 0))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ 2 * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (-((~B * ~c - ~A * ~d)) * (~b * ~c - ~a * ~d) ^ 2 * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~d) ^ 2 * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) + (1 / (~d * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(~B * (~b * ~c - ~a * ~d) ^ 2 + ~A * ~d * (((~a) ^ 2 * ~c - (~b) ^ 2 * ~c) + 2 * ~a * ~b * ~d) + ~d * (~B * (((~a) ^ 2 * ~c - (~b) ^ 2 * ~c) + 2 * ~a * ~b * ~d) + ~A * ((2 * ~a * ~b * ~c - (~a) ^ 2 * ~d) + (~b) ^ 2 * ~d)) * tan(~e + ~f * ~x) + (~b) ^ 2 * ~B * ((~c) ^ 2 + (~d) ^ 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && LtQ(~n, -1))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~b * ~c - ~a * ~d) * (~B * ~c - ~A * ~d) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~d * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 2) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(((~a * ~A * ~d * (~b * ~d * (~m - 1) - ~a * ~c * (~n + 1)) + (~b * ~B * ~c - (~A * ~b + ~a * ~B) * ~d) * (~b * ~c * (~m - 1) + ~a * ~d * (~n + 1))) - ~d * ((~a * ~A - ~b * ~B) * (~b * ~c - ~a * ~d) + (~A * ~b + ~a * ~B) * (~a * ~c + ~b * ~d)) * (~n + 1) * tan(~e + ~f * ~x)) - ~b * (~d * ((~A * ~b * ~c + ~a * ~B * ~c) - ~a * ~A * ~d) * (~m + ~n) - ~b * ~B * ((~c) ^ 2 * (~m - 1) - (~d) ^ 2 * (~n + 1))) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 1) && (LtQ(~n, -1) && (IntegerQ(~m) || IntegersQ(2 * ~m, 2 * ~n)))))))
    @apply_utils Antiderivative(((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ 2 * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)), ~x) => ((~b) ^ 2 * ~B * tan(~e + ~f * ~x)) / (~d * ~f) + (1 / ~d) * Antiderivative((((~a) ^ 2 * ~A * ~d - (~b) ^ 2 * ~B * ~c) + (2 * ~a * ~A * ~b + ~B * ((~a) ^ 2 - (~b) ^ 2)) * ~d * tan(~e + ~f * ~x) + (~A * (~b) ^ 2 * ~d - ~b * ~B * (~b * ~c - 2 * ~a * ~d)) * tan(~e + ~f * ~x) ^ 2) / (~c + ~d * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~b * ~B * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n)) + (1 / (~d * (~m + ~n))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 2) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp((((~a) ^ 2 * ~A * ~d * (~m + ~n) - ~b * ~B * (~b * ~c * (~m - 1) + ~a * ~d * (~n + 1))) + ~d * (~m + ~n) * (2 * ~a * ~A * ~b + ~B * ((~a) ^ 2 - (~b) ^ 2)) * tan(~e + ~f * ~x)) - (~b * ~B * (~b * ~c - ~a * ~d) * (~m - 1) - ~b * (~A * ~b + ~a * ~B) * ~d * (~m + ~n)) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 1) && ((IntegerQ(~m) || IntegersQ(2 * ~m, 2 * ~n)) && Not(IGtQ(~n, 1) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0))))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A * ~b - ~a * ~B) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n) / (~f * (~m + 1) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / (~b * (~m + 1) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n - 1) * Simp(((~b * ~B * (~b * ~c * (~m + 1) + ~a * ~d * ~n) + ~A * ~b * (~a * ~c * (~m + 1) - ~b * ~d * ~n)) - ~b * (~A * (~b * ~c - ~a * ~d) - ~B * (~a * ~c + ~b * ~d)) * (~m + 1) * tan(~e + ~f * ~x)) - ~b * ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 1) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (LtQ(~m, -1) && (LtQ(0, ~n, 1) && (IntegerQ(~m) || IntegersQ(2 * ~m, 2 * ~n)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~b * (~A * ~b - ~a * ~B) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~b * ~B * (~b * ~c * (~m + 1) + ~a * ~d * (~n + 1)) + ~A * (~a * (~b * ~c - ~a * ~d) * (~m + 1) - (~b) ^ 2 * ~d * (~m + ~n + 2))) - (~A * ~b - ~a * ~B) * (~b * ~c - ~a * ~d) * (~m + 1) * tan(~e + ~f * ~x)) - ~b * ~d * (~A * ~b - ~a * ~B) * (~m + ~n + 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (LtQ(~m, -1) && ((IntegerQ(~m) || IntegersQ(2 * ~m, 2 * ~n)) && Not(ILtQ(~n, -1) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0))))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => (~B * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n) / (~f * (~m + ~n)) + (1 / (~m + ~n)) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n - 1) * Simp((~a * ~A * ~c * (~m + ~n) - ~B * (~b * ~c * ~m + ~a * ~d * ~n)) + ((~A * ~b * ~c + ~a * ~B * ~c + ~a * ~A * ~d) - ~b * ~B * ~d) * (~m + ~n) * tan(~e + ~f * ~x) + (~A * ~b * ~d * (~m + ~n) + ~B * (~a * ~d * ~m + ~b * ~c * ~n)) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (LtQ(0, ~m, 1) && LtQ(0, ~n, 1)))))
    @apply_utils Antiderivative((~(A') + ~(B') * tan(~(e') + ~(f') * ~x)) / ((~a + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x))), ~x) => ((~B * (~b * ~c + ~a * ~d) + ~A * (~a * ~c - ~b * ~d)) * ~x) / (((~a) ^ 2 + (~b) ^ 2) * ((~c) ^ 2 + (~d) ^ 2)) + ((~b * (~A * ~b - ~a * ~B)) / ((~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~b - ~a * tan(~e + ~f * ~x)) / (~a + ~b * tan(~e + ~f * ~x)), ~x) + ((~d * (~B * ~c - ~A * ~d)) / ((~b * ~c - ~a * ~d) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~d - ~c * tan(~e + ~f * ~x)) / (~c + ~d * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative((sqrt(~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(Simp((~A * (~a * ~c + ~b * ~d) + ~B * (~b * ~c - ~a * ~d)) - (~A * (~b * ~c - ~a * ~d) - ~B * (~a * ~c + ~b * ~d)) * tan(~e + ~f * ~x), ~x) / sqrt(~c + ~d * tan(~e + ~f * ~x)), ~x) - (((~b * ~c - ~a * ~d) * (~B * ~a - ~A * ~b)) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((1 + tan(~e + ~f * ~x) ^ 2) / ((~a + ~b * tan(~e + ~f * ~x)) * sqrt(~c + ~d * tan(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative(((~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp((~a * ~A + ~b * ~B) - (~A * ~b - ~a * ~B) * tan(~e + ~f * ~x), ~x), ~x) + ((~b * (~A * ~b - ~a * ~B)) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 + tan(~e + ~f * ~x) ^ 2)) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative((sqrt(~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x))) / sqrt(~(c') + ~(d') * tan(~(e') + ~(f') * ~x)), ~x) => Antiderivative(Simp((~a * ~A - ~b * ~B) + (~A * ~b + ~a * ~B) * tan(~e + ~f * ~x), ~x) / (sqrt(~a + ~b * tan(~e + ~f * ~x)) * sqrt(~c + ~d * tan(~e + ~f * ~x))), ~x) + ~b * ~B * Antiderivative((1 + tan(~e + ~f * ~x) ^ 2) / (sqrt(~a + ~b * tan(~e + ~f * ~x)) * sqrt(~c + ~d * tan(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    #= Antiderivative(((~A')+(~B')*tan((~e')+(~f')*(~x)))/(sqrt((~a')+(~b')*tan((~e')+(~f')*(~x)))*sqrt( (~c')+(~d')*tan((~e')+(~f')*(~x)))),~x) :=  A^2/f*Subst(Antiderivative(1/((A-B*x)*sqrt(a+b*x)*sqrt(c+d*x)),x),x,tan(e+f*x)) <--  FreeQ([a,b,c,d,e,f,A,B],x) && NeQ(b*c-a*d,0) && NeQ(a^2+b^2,0) &&  NeQ(c^2+d^2,0) && EqQ(A^2+B^2,0) =#
    #= Antiderivative(((~A')+(~B')*tan((~e')+(~f')*(~x)))/(sqrt((~a')+(~b')*tan((~e')+(~f')*(~x)))*sqrt( (~c')+(~d')*tan((~e')+(~f')*(~x)))),~x) :=  (A+im*B)/2*Antiderivative((1-im*tan(e+f*x))/(sqrt(a+b*tan(e+f*x))*sqrt(c+d*tan(e+f* x))),x) + (A-im*B)/2*Antiderivative((1+im*tan(e+f*x))/(sqrt(a+b*tan(e+f*x))*sqrt(c+d*tan(e+ f*x))),x) <-- FreeQ([a,b,c,d,e,f,A,B],x) && NeQ(b*c-a*d,0) &&  NeQ(a^2+b^2,0) && NeQ(c^2+d^2,0) && NeQ(A^2+B^2,0) =#
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A) ^ 2 / ~f) * Subst(Antiderivative(((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n) / (~A - ~B * ~x), ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && EqQ((~A) ^ 2 + (~B) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x)), ~x) => ((~A + im * ~B) / 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 - im * tan(~e + ~f * ~x)), ~x) + ((~A - im * ~B) / 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 + im * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~A) ^ 2 + (~B) ^ 2, 0)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.3.1 (a+b tan)^m (c+d tan)^n (A+B tan).jl


text:
using NISTCrossSectionData
const P = NISTCrossSectionData
using Test
using PeriodicTable
using Unitful: cm, g, mg, MeV, eV

@testset "implementation details" begin
    @test P.linterpol(1.7, 1=>10, 2=>20) ≈ 17
    @test P.linterpol(1., 1=>10, 2=>20) ≈ 10
    @test P.linterpol(2., 1=>10, 2=>20) ≈ 20
end

@testset "ESTAR" begin
    unit = g/cm^2
    O = elements[8]
    E = 1MeV
    @test lookup(O, Electron(E), CSDA()) == 4.961E-01unit
    @test lookup(ESTAR, O, Electron(E), CSDA())                  == 4.961E-01unit

    E_max = 1e3MeV
    E_min = 1e-2*MeV
    @test lookup(ESTAR, O, Electron(E_max), CSDA()) == 9.454E+01unit
    @test lookup(ESTAR, O, Electron(E_min), CSDA()) == 2.950E-04unit

    Cf = elements[98]
    @test lookup(ESTAR, Cf, Electron(E), CSDA()) == 8.026E-01unit

    pt = Electron(1MeV)
    mat = :Si
    for (proc, value) in [
        (Collision(), 1.507MeV*cm^2/g)
        (Radiative(), 2.360E-02*MeV*cm^2/g)
        (Total(), 1.531E+00*MeV*cm^2/g)
        (CSDA(), 5.386E-01*g*cm^-2)
        (RadiationYield(), 8.275E-03)
        (DensityEffect(), 3.424E-01)
            ]

        @test lookup(mat, pt, proc) ≈ value
        @test lookup(ESTAR, mat, pt, proc) ≈ value
    end
end

@testset "XAAMDI" begin
    unit = cm^2/g
    O = elements[8]
    E = 1MeV
    @test mass_coeff(O, E, EnergyLoss(),datasource=XAAMDI     ) == 0.02794unit
    @test mass_coeff(O, E, EnergyLoss()) == 0.02794unit

    @test mass_coeff(:H2O, 1MeV, EnergyLoss()) == 3.103E-02*cm^2/g
end
@testset "XCOM" begin
    unit = cm^2/g
    V = elements[23]
    E = 2.044MeV
    for (PT, value) in [
        (CoherentScattering                       , 8.268E-05*cm^2/g)
        (IncoherentScattering                     , 3.934E-02*cm^2/g)
        (PhotoelectricAbsorption                  , 5.986E-05*cm^2/g)
        (PairProductionNuclearField               , 1.248E-03*cm^2/g)
        (PairProductionElectronField              , 0.000E+00*cm^2/g)
        (TotalAttenuation   , 4.073E-02*cm^2/g)
        (TotalAttenuationWithoutCoherentScattering, 4.065E-02*cm^2/g)]
        @test mass_coeff(V,E,PT(),datasource=XCOM) == value
        @test mass_coeff(V,E,PT()) == value
    end
    Fm = elements[100]
    @test mass_coeff(Fm, 1MeV, TotalAttenuation(), datasource=XCOM) == 0.08941*cm^2/g

    @test mass_coeff(:H2O, 1MeV, TotalAttenuation()) == 0.07072cm^2/g
end

@testset "mass_coeff" begin
    unit = cm^2/g
    O = elements[8]
    @test mass_coeff(O, 1MeV, EnergyLoss()) == 0.02794unit
    @test mass_coeff(O, 1MeV, EnergyLoss()) == 0.02794unit

    U = elements[92]
    # extremal energies
    E_min = 0.001MeV
    E_max = 20MeV
    @test mass_coeff(U, E_min, TotalAttenuation(), datasource=XAAMDI) == 6626.0unit
    @test mass_coeff(U, E_min, EnergyLoss()) == 6612.0unit

    @test mass_coeff(U, E_max, TotalAttenuation(), datasource=XAAMDI) == 6.512e-2unit
    @test mass_coeff(U, E_max, EnergyLoss()) == 3.662e-2unit

    # absorption edge
    E_K = 0.115606MeV
    @test_throws ArgumentError mass_coeff(U, E_K, TotalAttenuation(), datasource=XAAMDI)
    E_K₋ = E_K - 1eV
    E_K₊ = E_K + 1eV
    @test mass_coeff(U, E_K₋, TotalAttenuation(), datasource=XAAMDI) ≈ 1.378unit rtol=1e-2
    @test mass_coeff(U, E_K₋, EnergyLoss(), datasource=XAAMDI)       ≈ 1.027unit rtol=1e-2
    @test mass_coeff(U, E_K₊, TotalAttenuation(), datasource=XAAMDI) ≈ 4.893unit rtol=1e-2
    @test mass_coeff(U, E_K₊, EnergyLoss(), datasource=XAAMDI)       ≈ 1.382unit rtol=1e-2
end

@testset "api" begin
    oxygen = elements[8];
    @inferred mass_coeff(      oxygen,  1MeV, EnergyLoss())
    @inferred mean_free_path(  oxygen,  1MeV, TotalAttenuation())
    @inferred mean_free_path(  oxygen,  1MeV, TotalAttenuation(), density=1mg*cm^-3)
    @inferred transmission(  oxygen,  1MeV, TotalAttenuation(), length=1cm)
    @inferred transmission(  :O,  1MeV, TotalAttenuation(), length=1cm)
    @inferred attenuation_coeff(oxygen, 10MeV, EnergyLoss())
    @inferred cross_section(  oxygen,  1MeV, TotalAttenuation())
    @inferred lookup( oxygen, 1MeV, CoherentScattering())

    proc = CoherentScattering()
    pt = 1MeV
    @test lookup(oxygen, pt, proc) == lookup(:O, pt, proc)
    @test lookup(oxygen, pt, proc) == lookup("oxygen", pt, proc)
    @test lookup(oxygen, pt, proc) == lookup(8, pt, proc)
    @test lookup(oxygen, pt, TotalAttenuation()) ==
        lookup(8, pt, TotalAttenuation())
end