text:
using DelimitedFiles
using Statistics

open("data_3.dat") do f

  # line_number
  line = 1

  data = zeros(1000, 12)
  # read till end of file
  while ! eof(f)

     # read a new / next line for every itera
     s = readline(f)
     for ii in 1:12
         data[line, ii] = parse(Int, s[ii])
     end
     line += 1
  end
  gamma = []
  epsilon = []
  for ii in 1:12
      if mean(data[:,ii]) < 0.5
          append!(gamma, 0)
          append!(epsilon, 1)
      else
          append!(gamma, 1)
          append!(epsilon, 0)
      end
   end

   binary_g = 0
   for g in gamma
      binary_g = binary_g * 2 + g
   end

   binary_e = 0
   for e in epsilon
       binary_e = binary_e * 2 + e
   end
   println(binary_e * binary_g)
end


text:
module SweepContractor
    import Base: show, print, summary
    export Tensor, TensorNetwork, LabelledTensorNetwork
    # Define the types
    include("tensor_network.jl")

    import Base: By, Lt, Ordering, Forward, ForwardOrdering, ReverseOrdering, lt
    include("sorted_data_structs/data_structures.jl")
    include("sorted_data_structs/avl_tree.jl")
    import LinearAlgebra
    import Base: issorted, sortperm, sort!, isless, intersect
    # Utilities to make sure networks are suitably formatted for sweep contraction
    include("tn_utilities.jl")

    export sweep_contract!, sweep_contract
    export sweep_contract_dangling, sweep_contract_dangling!
    # The sweep contraction algorithm
    include("sweep_contract.jl")
    include("sweep_contract_dangling.jl")
end

ext: jl


text:
# ps3000_save_streaming_data.jl

export ps3000_save_streaming_data

function ps3000_save_streaming_data(handle :: Int16,
                                    callback :: Ptr{Void},
                                    data_buffer_size :: Int16 = 0x4000)
  buffer = Array(Int16,data_buffer_size)
  undefined_return_value = ccall((:ps3000_save_streaming_data, ps3000driver), Int16,
    (Int16, Ptr{Void}, Ptr{Int16}, Int16),
    handle, callback, buffer, data_buffer_size)
  return nothing
end


text:
### Image Processing

import Base: hist, scale

# Export methods

export  loadImage,
        saveImage,
        isImageIOAvailable,
        colorspace,
        gray2rgb,
        rgb2gray,
        rgb2hsv,
        rgb2ycbcr,
        ycbcr2rgb,
        hsv2rgb,
        regions,
        SAT,
        bilateral,
        maxfilt,
        meanshift,
        medfilt,
        minfilt,
        sobel,
        histequal,
        resize,
        rotate,
        skew,
        transform,
        transformCoordinates,
        translate,
        dilate,
        erode,
        dilate3d,
        erode3d,
        gaussiankernel

# Export constants

export  AF_GRAY,
        AF_RGB,
        AF_HSV,
        AF_YCbCr,
        AF_YCC_601,
        AF_YCC_709,
        AF_YCC_2020,
        AF_CONNECTIVITY_4,
        AF_CONNECTIVITY_8,
        AF_PAD_ZERO,
        AF_PAD_SYM,
        AF_INTERP_NEAREST,
        AF_INTERP_LINEAR,
        AF_INTERP_BILINEAR,
        AF_INTERP_CUBIC,
        AF_INTERP_LOWER

# Colorspaces and Constants

AF_GRAY = 0
AF_RGB = 1
AF_HSV = 2
AF_YCbCr = 3

AF_YCC_601 = 601
AF_YCC_709 = 709
AF_YCC_2020 = 2020

AF_CONNECTIVITY_4 = 4
AF_CONNECTIVITY_8 = 8

AF_PAD_ZERO = 0
AF_PAD_SYM = 1

AF_INTERP_NEAREST = 0
AF_INTERP_LINEAR = 1
AF_INTERP_BILINEAR = 2
AF_INTERP_CUBIC = 3
AF_INTERP_LOWER = 4

# Read and write images

function loadImage(path::AbstractString; color = false)
    out = new_ptr()
    af_load_image(out, path, color)
    AFArray{backend_eltype(out[])}(out[])
end

function saveImage(path::AbstractString, a::AFArray)
    af_save_image(path, a)
end

function isImageIOAvailable()
    out = Base.Ref{Bool}(0)
    af_is_image_io_available(out)
    out[]
end

# Colorspace conversions

function colorspace(img::AFArray, to::Int, from::Int)
    out = new_ptr()
    af_colorspace(out, img, to, from)
    AFArray{backend_eltype(out[])}(out[])
end

function gray2rgb(img::AFArray; rFactor::Cdouble = 1., gFactor::Cdouble = 1., bFactor::Cdouble = 1.)
    out = new_ptr()
    af_gray2rgb(out, img, rFactor, gFactor, bFactor)
    AFArray{backend_eltype(out[])}(out[])
end

function hsv2rgb(img::AFArray)
    out = new_ptr()
    af_hsv2rgb(out, img)
    AFArray{backend_eltype(out[])}(out[])
end

function rgb2gray(img::AFArray; rPercent::Cfloat = 1., gPercent::Cfloat = 1., bFactor::Cfloat = 1.)
    out = new_ptr()
    af_rgb2gray(out, img)
    AFArray{backend_eltype(out[])}(out[])
end

function rgb2hsv(img::AFArray)
    out = new_ptr()
    af_rgb2hsv(out, img)
    AFArray{backend_eltype(out[])}(out[])
end

function rgb2ycbcr(img::AFArray; std::Int = AF_YCC_601)
    out = new_ptr()
    af_rgb2ycbcr(out, img, std)
    AFArray{backend_eltype(out[])}(out[])
end

function ycbcr2rgb(img::AFArray; std::Int = AF_YCC_601)
    out = new_ptr()
    af_ycbcr2rgb(out, img, std)
    AFArray{backend_eltype(out[])}(out[])
end

# Regions

function regions(a::AFArray; conn = AF_CONNECTIVITY_4, typ = Float32)
    out = new_ptr()
    af_regions(out, a, conn, typ)
    AFArray{typ}(out[])
end

# Filters

function SAT(a::AFArray)
    out = new_ptr()
    af_sat(out, a)
    AFArray{backend_eltype(out[])}(out[])
end

function bilateral(a::AFArray, spatial_sigma::Real, chromatic_sigma::Real; color = false)
    out = new_ptr()
    af_bilateral(out, a, spatial_sigma, chromatic_sigma, color)
    AFArray{backend_eltype(out[])}(out[])
end

function maxfilt(a::AFArray; wind_length = 3, wind_width = 3, edge_pad = AF_PAD_ZERO)
    out = new_ptr()
    af_maxfilt(out, a, wind_length, wind_width, edge_pad)
    AFArray{backend_eltype(out[])}(out[])
end

function meanshift(a::AFArray, spatial_sigma::Real, chromatic_sigma::Real, iter::Int; color = false)
    out = new_ptr()
    af_meanshift(out, a, spatial_sigma, chromatic_sigma, Cuint(iter), color)
    AFArray{backend_eltype(out[])}(out[])
end

function medfilt(a::AFArray; wind_length = 3, wind_width = 3, edge_pad = AF_PAD_ZERO)
    out = new_ptr()
    af_medfilt(out, a, wind_length, wind_width, edge_pad)
    AFArray{backend_eltype(out[])}(out[])
end

function minfilt(a::AFArray; wind_length = 3, wind_width = 3, edge_pad = AF_PAD_ZERO)
    out = new_ptr()
    af_minfilt(out, a, wind_length, wind_width, edge_pad)
    AFArray{backend_eltype(out[])}(out[])
end

function sobel(a::AFArray; ker_size = 3)
    dx, dy = new_ptr(), new_ptr()
    af_sobel_operator(dx, dy, a, Cuint(ker_size))
    AFArray{backend_eltype(dx[])}(dx[]), AFArray{backend_eltype(dy[])}(dy[])
end

# Histograms

function histequal(a::AFArray, hist::AFArray)
    out = new_ptr()
    af_hist_equal(out, a, hist)
    AFArray{backend_eltype(out[])}(out[])
end

function hist(a::AFArray, nbins::Int, minval::Cdouble, maxval::Cdouble)
    out = new_ptr()
    af_histogram(out, a, Cuint(nbins), minval, maxval)
    AFArray{backend_eltype(out[])}(out[])
end

# Image Transformations

function resize(a::AFArray, dim1::Int, dim2::Int;
                method = AF_INTERP_NEAREST)
    out = new_ptr()
    af_resize(out, a, dim1, dim2, method)
    AFArray{backend_eltype(out[])}(out[])
end

function rotate(a::AFArray, theta::Real; crop::Bool = true,
                method = AF_INTERP_NEAREST)
    out = new_ptr()
    af_rotate(out, a, theta, crop, method)
    AFArray{backend_eltype(out[])}(out[])
end

function scale(a::AFArray, scale0::Real, scale1::Real;
                dim1 = 0, dim2 = 0, method = AF_INTERP_NEAREST)
    out = new_ptr()
    af_scale(out, a, scale0, scale1, dim1, dim2, method)
    AFArray{backend_eltype(out[])}(out[])
end

function skew(a::AFArray, skew1::Real, skew2::Real;
                dim1 = 0, dim2 = 0, method = AF_INTERP_NEAREST,
                inverse = true)
    out = new_ptr()
    af_skew(out, a, skew1, skew2, dim1, dim2, method, inverse)
    AFArray{backend_eltype(out[])}(out[])
end

function transform(a::AFArray, transform::AFArray;
                    dim1 = 0, dim2 = 0, method = AF_INTERP_NEAREST,
                    inverse = true)
    out = new_ptr()
    af_transform(out, a, transform, dim1, dim2, method, inverse)
    AFArray{backend_eltype(out[])}(out[])
end

function transformCoordinates(a::AFArray, d1::Real, d2::Real)
    out = new_ptr()
    af_transform_coordinates(out, a, d1, d2)
    AFArray{backend_eltype(out[])}(out[])
end

function translate(a::AFArray, trans1::Real, trans2::Real;
                    dim1 = 0, dim2 =0, method = AF_INTERP_NEAREST)
    out = new_ptr()
    af_translate(out, a, trans1, trans2, dim1, dim2, method)
    AFArray{backend_eltype(out[])}(out[])
end

# Morphological Operations

for (op, fn) in ((:dilate, :af_dilate), (:dilate3d, :af_dilate3d),
                (:erode, :af_erode), (:erode3d, :af_erode3d))

    @eval function ($op)(a::AFArray, mask::AFArray)
        out = new_ptr()
        eval($fn)(out, a, mask)
        AFArray{backend_eltype(out[])}(out[])
    end

end

# Gaussian Kernel

function gaussiankernel(rows::Int, cols::Int; sigma_r = 0., sigma_c = 0.)
    out = new_ptr()
    af_gaussian_kernel(out, rows, cols, sigma_r, sigma_c)
    AFArray{backend_eltype(out[])}(out[])
end


text:
"""
This struct contains data for the REC market product.
    name: REC.
    certificates: Expected output (number of certificates) for REC market bids
    rec_bid: REC market bidding price
"""
mutable struct REC <: Product
    name::Symbol
    certificates::Float64
    rec_bid::Float64
end

# Certificates and REC bids only returned when product is of type REC
get_rec_certificates(prod::Product) = nothing
get_rec_certificates(prod::REC) = prod.certificates
get_rec_bid(prod::Product) = nothing
get_rec_bid(prod::REC) = prod.rec_bid

# REC bids and certificates only set when product is of type REC
function set_rec_certificates!(prod::T, certificates) where T <: Product
    return
end

function set_rec_certificates!(prod::REC, certificates)
    prod.certificates = certificates
    return
end

function set_rec_bid!(prod::T, rec_bid) where T <: Product
    return
end

function set_rec_bid!(prod::REC, rec_bid)
    prod.rec_bid = rec_bid
    return
end


text:
import Observables
import AbstractPlotting

# Map from S² into its upper hemisphere
s2tos2map(b::S²) = begin
    p = Geographic(b)
    r = sqrt((1 - sin(p.θ)) / 2)
    Geographic(p.r, r * cos(p.ϕ), r * sin(p.ϕ))
end

scene = AbstractPlotting.Scene()
number = rand(5:10)
points = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
top = U1(rand() * 2pi - pi)
bottom = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
color = AbstractPlotting.RGBAf0(rand(4)...)
transparency = rand(1:2) == 1 ? true : false
s2tos3map = rand() > 0.5 ? σmap : τmap
scale = 1.0 + rand()
whirl = Whirl(scene,
              points,
              s2tos3map,
              s2tos2map,
              top = top,
              bottom = bottom,
              s3rotation = s3rotation,
              config = config,
              segments = segments,
              color = color,
              transparency = transparency,
              scale = scale)

## update points

points2 = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
update(whirl, points2)

@test isapprox(whirl.points, points2)

# update maps f: S² → S³ and g: S² → S²

s2tos3map = rand() > 0.5 ? σmap : τmap
s2tos2map′(x) = Geographic(0.5 * x.r, 0.5 * x.ϕ, 0.5 * x.θ)
update(whirl, s2tos3map, s2tos2map′)
p = Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2)

@test isapprox(whirl.s2tos3map(p), s2tos3map(p))
@test isapprox(whirl.s2tos2map(p), s2tos2map′(p))

## update top and bottom

top2 = U1(rand() * 2pi - pi)
bottom2 = U1(rand() * 2pi - pi)
update(whirl, top2, bottom2)

@test isapprox(whirl.top, top2)
@test isapprox(whirl.bottom, bottom2)

## update S³ rotation

s3rotation2 = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
update(whirl, s3rotation2)

@test isapprox(whirl.s3rotation, s3rotation2)

## update configuration

config2 = config * Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
value1 = getsurface(whirl.observable, whirl.segments, length(whirl.points))
update(whirl, config2)
value2 = getsurface(whirl.observable, whirl.segments, length(whirl.points))

@test isapprox(whirl.config, config2)
@test isapprox(value1, value2) == false

color1 = AbstractPlotting.RGBAf0(rand(4)...)
update(whirl, color1)
color2 = Observables.to_value(whirl.color)[1] # Select element 1

@test isapprox(color1, color2)

## bulk update for reducing extra computation

points3 = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
s3rotation3 = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
top3 = U1(rand() * 2pi - pi)
bottom3 = U1(rand() * 2pi - pi)
update(whirl, points3, s3rotation3, top3, bottom3)

@test isapprox(whirl.points, points3)
@test isapprox(whirl.s3rotation, s3rotation3)
@test isapprox(whirl.top, top3)
@test isapprox(whirl.bottom, bottom3)


text:
# Simple Parsers
# ==============

"""
    EP.SatisfiesPredicate(predicate_func[, errormessage])

Construct an ExprParser which checks whether the given `predicate_func` returns true.
If so, the to-be-parsed value is returned as such, otherwise an [`ParseError`](@ref) is thrown as usual.
If `errormessage` is given, it will be appended to the default error message.

# Example
```jldoctest
julia> using ExprParsers;

julia> parser = EP.SatisfiesPredicate(isodd);

julia> parse_expr(parser, 3)
3
julia> parse_expr(parser, 4)
ERROR: ParseError: Predicate `isodd` returned false on expr `4`.
julia> is44(x) = x==44;

julia> parser2 = EP.SatisfiesPredicate(is44, "It should be 44.");

julia> parse_expr(parser2, 44)
44
julia> parse_expr(parser2, 4)
ERROR: ParseError: Predicate `is44` returned false on expr `4`. It should be 44.
```
"""
struct SatisfiesPredicate{Func} <: ExprParser
  predicate::Func
  errormessage
end
SatisfiesPredicate(func) = SatisfiesPredicate(func, "")
function parse_expr(parser::SatisfiesPredicate, expr)
  parser.predicate(expr) ? expr : throw(ParseError(
    "Predicate `$(parser.predicate)` returned false on expr `$expr`. $(parser.errormessage)"
  ))
end

"""
    EP.Isa(T::Type)

Constructs an ExprParser which checks whether a value is of the given type.

# Example
```jldoctest
julia> using ExprParsers

julia> parser = EP.Isa(Symbol);

julia> parse_expr(parser, :thisisasymbol)
:thisisasymbol
julia> parse_expr(parser, 42)
ERROR: ParseError: Expected type `Symbol`, got `42` of type `Int64`.
```

`EP.Isa(Symbol)` is so common that there is a special constant for it [`anysymbol`](@ref).
"""
struct Isa{T} <: ExprParser
  # this way we cannot forget to instantiate the Type
  Isa(::Base.Type{T}) where T = new{T}()
end
parse_expr(parser::Isa{T}, expr::T) where T = expr
parse_expr(parser::Isa{T}, other::S) where {T, S} = throw(ParseError("Expected type `$T`, got `$other` of type `$S`."))

"""
    EP.anything = Isa(Any)

Special constant ExprParser which matches literally anything.

# Examples
```jldoctest
julia> using ExprParsers

julia> parse_expr(EP.anything, 42)
42
julia> parse_expr(EP.anything, :whatever)
:whatever
```
"""
const anything = Isa(Any)

"""
    EP.anysymbol = Isa(Symbol)

Special constant ExprParser which matches Symbols, and only Symbols.

# Examples
```jldoctest
julia> using ExprParsers

julia> parse_expr(EP.anysymbol, :asymbol)
:asymbol
julia> parse_expr(EP.anysymbol, 42)
ERROR: ParseError: Expected type `Symbol`, got `42` of type `Int64`.
```
"""
const anysymbol = Isa(Symbol)

"""
    EP.AnyOf(parser1, parser2, parser3, ...; errormessage = "")

Constructs an ExprParser from multiple given parsers. When given a value it first tries to match
`parser1`, and if that fails with a ParseError, then `parser2`, and so forth.
The result from the first parser which matches will be returned.
If no parser matches, a dedicated ParseError is raised.

If `errormessage` is given, it will be appended to the default error message in case of ParseError.

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.AnyOf(EP.anysymbol, EP.Isa(String), EP.SatisfiesPredicate(isodd),
                         errormessage="My error message.");

julia> parse_expr(parser, :hi)
:hi
julia> parse_expr(parser, 3)
3
julia> parse_expr(parser, "something")
"something"
julia> parse_expr(parser, 4)
ERROR: ParseError: AnyOf could not parse expr `4` with any of the parsers `(ExprParsers.Isa{Symbol}(), ExprParsers.Isa{String}(), ExprParsers.SatisfiesPredicate{typeof(isodd)}(isodd, ""))`. My error message.
```
"""
struct AnyOf{T} <: ExprParser
  several::T
  errormessage
  AnyOf(several...; errormessage = "") = new{typeof(several)}(several, errormessage)
end

function parse_expr(parser::AnyOf, expr)
  for subparser in parser.several
    try
      return parse_expr(subparser, expr)
    catch e
      # only continue if ParseError
      e isa ParseError || rethrow()
    end
  end
  throw(ParseError("AnyOf could not parse expr `$(expr)` with any of the parsers `$(parser.several)`. $(parser.errormessage)"))
end

"""
    EP.AllOf(parser1, parser2, parser3, ...)

Constructs an ExprParser from multiple given parsers. When to match a value, all parsers actually need to
parse correctly, otherwise the `ParseError` from the first non-matching parser is rethrown.
If all parsers match, then the return value from the last parser is returned.

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.AllOf(EP.Isa(Number), EP.SatisfiesPredicate(isodd), 3);

julia> parse_expr(parser, 3)
3
julia> parse_expr(parser, "something")
ERROR: ParseError: Expected type `Number`, got `something` of type `String`.
julia> parse_expr(parser, 4)
ERROR: ParseError: Predicate `isodd` returned false on expr `4`.
julia> parse_expr(parser, 5)
ERROR: ParseError: Using default `==` comparison, but parser `3` ≠ value `5`.
```
"""
struct AllOf{T} <: ExprParser
  several::T
  AllOf(several...) = new{typeof(several)}(several)
end

function parse_expr(parser::AllOf, expr)
  for subparser in parser.several[1:end-1]
    parse_expr(subparser, expr)
  end
  # return last match, if all else matched
  parse_expr(parser.several[end], expr)
end


text:
module Cartesian

export @ngenerate, @nsplat, @nloops, @nref, @ncall, @nexprs, @nextract, @nall, @ntuple, @nif, ngenerate

const CARTESIAN_DIMS = 4

### @ngenerate, for auto-generation of separate versions of functions for different dimensionalities
# Examples (deliberately trivial):
#     @ngenerate N returntype myndims{T,N}(A::Array{T,N}) = N
# or alternatively
#     function gen_body(N::Int)
#         quote
#             return $N
#         end
#     end
#     eval(ngenerate(:N, returntypeexpr, :(myndims{T,N}(A::Array{T,N})), gen_body))
# The latter allows you to use a single gen_body function for both ngenerate and
# when your function maintains its own method cache (e.g., reduction or broadcasting).
#
# Special syntax for function prototypes:
#   @ngenerate N returntype function myfunction(A::AbstractArray, I::NTuple{N, Int}...)
# for N = 3 translates to
#   function myfunction(A::AbstractArray, I_1::Int, I_2::Int, I_3::Int)
# and for the generic (cached) case as
#   function myfunction(A::AbstractArray, I::Int...)
#     @nextract N I I
# with N = length(I). N should _not_ be listed as a parameter of the function unless
# earlier arguments use it that way.
# To avoid ambiguity, it would be preferable to have some specific syntax for this, such as
#   myfunction(A::AbstractArray, I::Int...N)
# where N can be an integer or symbol. Currently T...N generates a parser error.
macro ngenerate(itersym, returntypeexpr, funcexpr)
    isfuncexpr(funcexpr) || error("Requires a function expression")
    esc(ngenerate(itersym, returntypeexpr, funcexpr.args[1], N->sreplace!(copy(funcexpr.args[2]), itersym, N)))
end

# @nsplat takes an expression like
#    @nsplat N 2:3 myfunction(A, I::NTuple{N,Real}...) = getindex(A, I...)
# and generates
#    myfunction(A, I_1::Real, I_2::Real) = getindex(A, I_1, I_2)
#    myfunction(A, I_1::Real, I_2::Real, I_3::Real) = getindex(A, I_1, I_2, I_3)
#    myfunction(A, I::Real...) = getindex(A, I...)
# An @nsplat function _cannot_ have any other Cartesian macros in it.
# If you omit the range, it uses 1:CARTESIAN_DIMS.
macro nsplat(itersym, args...)
    local rng
    if length(args) == 1
        rng = 1:CARTESIAN_DIMS
        funcexpr = args[1]
    elseif length(args) == 2
        rangeexpr = args[1]
        funcexpr = args[2]
        if !isa(rangeexpr, Expr) || rangeexpr.head != :(:) || length(rangeexpr.args) != 2
            error("First argument must be a from:to expression")
        end
        rng = rangeexpr.args[1]:rangeexpr.args[2]
    else
        error("Wrong number of arguments")
    end
    isfuncexpr(funcexpr) || error("Second argument must be a function expression")
    prototype = funcexpr.args[1]
    body = funcexpr.args[2]
    varname, T = get_splatinfo(prototype, itersym)
    isempty(varname) && error("Last argument must be a splat")
    explicit = [Expr(:function, resolvesplat!(copy(prototype), varname, T, N),
                     resolvesplats!(copy(body), varname, N)) for N in rng]
    protosplat = resolvesplat!(copy(prototype), varname, T, 0)
    protosplat.args[end] = Expr(:..., protosplat.args[end])
    splat = Expr(:function, protosplat, body)
    esc(Expr(:block, explicit..., splat))
end

generate1(itersym, prototype, bodyfunc, N::Int, varname, T) =
    Expr(:function, spliceint!(sreplace!(resolvesplat!(copy(prototype), varname, T, N), itersym, N)),
         resolvesplats!(bodyfunc(N), varname, N))

function ngenerate(itersym, returntypeexpr, prototype, bodyfunc, dims=1:CARTESIAN_DIMS, makecached::Bool = true)
    varname, T = get_splatinfo(prototype, itersym)
    # Generate versions for specific dimensions
    fdim = [generate1(itersym, prototype, bodyfunc, N, varname, T) for N in dims]
    if !makecached
        return Expr(:block, fdim...)
    end
    # Generate the generic cache-based version
    if isempty(varname)
        setitersym, extractvarargs = :(), N -> nothing
    else
        s = symbol(varname)
        setitersym = hasparameter(prototype, itersym) ? (:(@assert $itersym == length($s))) : (:($itersym = length($s)))
        extractvarargs = N -> Expr(:block, map(popescape, _nextract(N, s, s).args)...)
    end
    fsym = funcsym(prototype)
    dictname = symbol(string(fsym)*"_cache")
    fargs = funcargs(prototype)
    if !isempty(varname)
        fargs[end] = Expr(:..., fargs[end].args[1])
    end
    flocal = funcrename(copy(prototype), :_F_)
    F = Expr(:function, resolvesplat!(prototype, varname, T), quote
             $setitersym
             if !haskey($dictname, $itersym)
                 gen1 = Base.Cartesian.generate1($(symbol(itersym)), $(Expr(:quote, flocal)), $bodyfunc, $itersym, $varname, $T)
                 $(dictname)[$itersym] = eval(quote
                     local _F_
                     $gen1
                     _F_
                 end)
             end
             ($(dictname)[$itersym]($(fargs...)))::$returntypeexpr
         end)
    Expr(:block, fdim..., quote
            let $dictname = Dict{Int,Function}()
            $F
            end
        end)
end

isfuncexpr(ex::Expr) =
    ex.head == :function || (ex.head == :(=) && typeof(ex.args[1]) == Expr && ex.args[1].head == :call)
isfuncexpr(arg) = false

sreplace!(arg, sym, val) = arg
function sreplace!(ex::Expr, sym, val)
    for i = 1:length(ex.args)
        ex.args[i] = sreplace!(ex.args[i], sym, val)
    end
    ex
end
sreplace!(s::Symbol, sym, val) = s == sym ? val : s

# If using the syntax that will need "desplatting",
#     myfunction(A::AbstractArray, I::NTuple{N, Int}...)
# return the variable name (as a string) and type
function get_splatinfo(ex::Expr, itersym::Symbol)
    if ex.head == :call
        a = ex.args[end]
        if  isa(a, Expr) && a.head == :... && length(a.args) == 1
            b = a.args[1]
            if isa(b, Expr) && b.head == :(::)
                varname = string(b.args[1])
                c = b.args[2]
                if isa(c, Expr) && c.head == :curly && c.args[1] == :NTuple && c.args[2] == itersym
                    T = c.args[3]
                    return varname, T
                end
            end
        end
    end
    "", Nothing
end

# Replace splatted with desplatted for a specific number of arguments
function resolvesplat!(prototype, varname, T::Union(Type,Symbol,Expr), N::Int)
    if !isempty(varname)
        prototype.args[end] = N > 0 ? Expr(:(::), symbol(string(varname, "_1")), T) :
                                      Expr(:(::), symbol(varname), T)
        for i = 2:N
            push!(prototype.args, Expr(:(::), symbol(string(varname, "_", i)), T))
        end
    end
    prototype
end

# Return the generic splatting form, e.g.,
#     myfunction(A::AbstractArray, I::Int...)
function resolvesplat!(prototype, varname, T::Union(Type,Symbol,Expr))
    if !isempty(varname)
        svarname = symbol(varname)
        prototype.args[end] = Expr(:..., :($svarname::$T))
    end
    prototype
end

# Desplatting function calls: replace func(a, b, I...) with func(a, b, I_1, I_2, I_3)
resolvesplats!(arg, varname, N) = arg
function resolvesplats!(ex::Expr, varname, N::Int)
    if ex.head == :call
        for i = 2:length(ex.args)-1
            resolvesplats!(ex.args[i], varname, N)
        end
        a = ex.args[end]
        if isa(a, Expr) && a.head == :... && a.args[1] == symbol(varname)
            ex.args[end] = symbol(string(varname, "_1"))
            for i = 2:N
                push!(ex.args, symbol(string(varname, "_", i)))
            end
        else
            resolvesplats!(a, varname, N)
        end
    else
        for i = 1:length(ex.args)
            resolvesplats!(ex.args[i], varname, N)
        end
    end
    ex
end

# Remove any function parameters that are integers
function spliceint!(ex::Expr)
    if ex.head == :escape
        return esc(spliceint!(ex.args[1]))
    end
    ex.head == :call || error(string(ex, " must be a call"))
    if isa(ex.args[1], Expr) && ex.args[1].head == :curly
        args = ex.args[1].args
        for i = length(args):-1:1
            if isa(args[i], Int)
                splice!(args, i)
            end
        end
    end
    ex
end

function popescape(ex::Expr)
    while ex.head == :escape
        ex = ex.args[1]
    end
    ex
end

# Extract the "function name"
function funcsym(prototype::Expr)
    prototype = popescape(prototype)
    prototype.head == :call || error(string(prototype, " must be a call"))
    tmp = prototype.args[1]
    if isa(tmp, Expr) && tmp.head == :curly
        tmp = tmp.args[1]
    end
    return tmp
end

function funcrename(prototype::Expr, name::Symbol)
    prototype = popescape(prototype)
    prototype.head == :call || error(string(prototype, " must be a call"))
    tmp = prototype.args[1]
    if isa(tmp, Expr) && tmp.head == :curly
        tmp.args[1] = name
    else
        prototype.args[1] = name
    end
    return prototype
end

function hasparameter(prototype::Expr, sym::Symbol)
    prototype = popescape(prototype)
    prototype.head == :call || error(string(prototype, " must be a call"))
    tmp = prototype.args[1]
    if isa(tmp, Expr) && tmp.head == :curly
        for i = 2:length(tmp.args)
            if tmp.args[i] == sym
                return true
            end
        end
    end
    false
end

# Extract the symbols of the function arguments
funcarg(s::Symbol) = s
funcarg(ex::Expr) = ex.args[1]
function funcargs(prototype::Expr)
    prototype = popescape(prototype)
    prototype.head == :call || error(string(prototype, " must be a call"))
    map(a->funcarg(a), prototype.args[2:end])
end

### Cartesian-specific macros

# Generate nested loops
macro nloops(N, itersym, rangeexpr, args...)
    _nloops(N, itersym, rangeexpr, args...)
end

_nloops(N::Int, itersym::Symbol, arraysym::Symbol, args::Expr...) = _nloops(N, itersym, :(d->1:size($arraysym,d)), args...)

function _nloops(N::Int, itersym::Symbol, rangeexpr::Expr, args::Expr...)
    if rangeexpr.head != :->
        error("Second argument must be an anonymous function expression to compute the range")
    end
    if !(1 <= length(args) <= 3)
        error("Too many arguments")
    end
    body = args[end]
    ex = Expr(:escape, body)
    for dim = 1:N
        itervar = inlineanonymous(itersym, dim)
        rng = inlineanonymous(rangeexpr, dim)
        preexpr = length(args) > 1 ? inlineanonymous(args[1], dim) : (:(nothing))
        postexpr = length(args) > 2 ? inlineanonymous(args[2], dim) : (:(nothing))
        ex = quote
            for $(esc(itervar)) = $(esc(rng))
                $(esc(preexpr))
                $ex
                $(esc(postexpr))
            end
        end
    end
    ex
end

# Generate expression A[i1, i2, ...]
macro nref(N, A, sym)
    _nref(N, A, sym)
end

function _nref(N::Int, A::Symbol, ex)
    vars = [ inlineanonymous(ex,i) for i = 1:N ]
    Expr(:escape, Expr(:ref, A, vars...))
end

# Generate f(arg1, arg2, ...)
macro ncall(N, f, sym...)
    _ncall(N, f, sym...)
end

function _ncall(N::Int, f, args...)
    pre = args[1:end-1]
    ex = args[end]
    vars = [ inlineanonymous(ex,i) for i = 1:N ]
    Expr(:escape, Expr(:call, f, pre..., vars...))
end

# Generate N expressions
macro nexprs(N, ex)
    _nexprs(N, ex)
end

function _nexprs(N::Int, ex::Expr)
    exs = [ inlineanonymous(ex,i) for i = 1:N ]
    Expr(:escape, Expr(:block, exs...))
end

# Make variables esym1, esym2, ... = isym
macro nextract(N, esym, isym)
    _nextract(N, esym, isym)
end

function _nextract(N::Int, esym::Symbol, isym::Symbol)
    aexprs = [Expr(:escape, Expr(:(=), inlineanonymous(esym, i), :(($isym)[$i]))) for i = 1:N]
    Expr(:block, aexprs...)
end

function _nextract(N::Int, esym::Symbol, ex::Expr)
    aexprs = [Expr(:escape, Expr(:(=), inlineanonymous(esym, i), inlineanonymous(ex,i))) for i = 1:N]
    Expr(:block, aexprs...)
end

# Check whether variables i1, i2, ... all satisfy criterion
macro nall(N, criterion)
    _nall(N, criterion)
end

function _nall(N::Int, criterion::Expr)
    if criterion.head != :->
        error("Second argument must be an anonymous function expression yielding the criterion")
    end
    conds = [Expr(:escape, inlineanonymous(criterion, i)) for i = 1:N]
    Expr(:&&, conds...)
end

macro ntuple(N, ex)
    _ntuple(N, ex)
end

function _ntuple(N::Int, ex)
    vars = [ inlineanonymous(ex,i) for i = 1:N ]
    Expr(:escape, Expr(:tuple, vars...))
end

# if condition1; operation1; elseif condition2; operation2; else operation3
# You can pass one or two operations; the second, if present, is used in the final "else"
macro nif(N, condition, operation...)
    # Handle the final "else"
    ex = esc(inlineanonymous(length(operation) > 1 ? operation[2] : operation[1], N))
    # Make the nested if statements
    for i = N-1:-1:1
        ex = Expr(:if, esc(inlineanonymous(condition,i)), esc(inlineanonymous(operation[1],i)), ex)
    end
    ex
end

## Utilities

# Simplify expressions like :(d->3:size(A,d)-3) given an explicit value for d
function inlineanonymous(ex::Expr, val)
    if ex.head != :->
        error("Not an anonymous function")
    end
    if !isa(ex.args[1], Symbol)
        error("Not a single-argument anonymous function")
    end
    sym = ex.args[1]
    ex = ex.args[2]
    exout = lreplace(ex, sym, val)
    exout = poplinenum(exout)
    exprresolve(exout)
end

# Given :i and 3, this generates :i_3
inlineanonymous(base::Symbol, ext) = symbol(string(base)*"_"*string(ext))

# Replace a symbol by a value or a "coded" symbol
# E.g., for d = 3,
#    lreplace(:d, :d, 3) -> 3
#    lreplace(:i_d, :d, 3) -> :i_3
#    lreplace(:i_{d-1}, :d, 3) -> :i_2
# This follows LaTeX notation.
lreplace(ex, sym::Symbol, val) = lreplace!(copy(ex), sym, val, Regex("_"*string(sym)*"(\$|(?=_))"))
lreplace!(arg, sym::Symbol, val, r) = arg
function lreplace!(s::Symbol, sym::Symbol, val, r::Regex)
    if (s == sym)
        return val
    end
    symbol(replace(string(s), r, "_"*string(val)))
end
function lreplace!(ex::Expr, sym::Symbol, val, r)
    # Curly-brace notation, which acts like parentheses
    if ex.head == :curly && length(ex.args) == 2 && isa(ex.args[1], Symbol) && endswith(string(ex.args[1]), "_")
        excurly = lreplace!(ex.args[2], sym, val, r)
        return symbol(string(ex.args[1])*string(exprresolve(excurly)))
    end
    for i in 1:length(ex.args)
        ex.args[i] = lreplace!(ex.args[i], sym, val, r)
    end
    ex
end

poplinenum(arg) = arg
function poplinenum(ex::Expr)
    if ex.head == :block
        if length(ex.args) == 1
            return ex.args[1]
        elseif length(ex.args) == 2 && ex.args[1].head == :line
            return ex.args[2]
        end
    end
    ex
end

exprresolve(arg) = arg
function exprresolve(ex::Expr)
    for i = 1:length(ex.args)
        ex.args[i] = exprresolve(ex.args[i])
    end
    # Handle simple arithmetic
    if ex.head == :call && in(ex.args[1], (:+, :-, :*, :/)) && all([isa(ex.args[i], Number) for i = 2:length(ex.args)])
        return eval(ex)
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-) && length(ex.args) == 3 && ex.args[3] == 0
        # simplify x+0 and x-0
        return ex.args[2]
    end
    # Resolve array references
    if ex.head == :ref && isa(ex.args[1], Array)
        for i = 2:length(ex.args)
            if !isa(ex.args[i], Real)
                return ex
            end
        end
        return ex.args[1][ex.args[2:end]...]
    end
    # Resolve conditionals
    if ex.head == :if
        try
            tf = eval(ex.args[1])
            ex = tf?ex.args[2]:ex.args[3]
        catch
        end
    end
    ex
end

end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module Docs

using ..Pages

export servedocs

function make(modulepath=@__FILE__;build="")
    makefile = normpath(joinpath(modulepath,"..","..","docs","make.jl"))
    if isfile(makefile) && isequal(basename(makefile),"make.jl")
        buildfolder = joinpath(dirname(makefile),build)
        if isfile(joinpath(buildfolder,"index.html"))
            Base.Filesystem.rm(joinpath(buildfolder,"index.html"))
            Base.Filesystem.rm(joinpath(buildfolder,"search_index.js"))
            Base.Filesystem.rm(joinpath(buildfolder,"assets"),recursive=true)
            Base.Filesystem.rm(joinpath(buildfolder,"search"),recursive=true)
        end

        include(makefile)
    else
        @warn "$(makefile) is not a valid makefile."
    end
end

function servedocs(modulepath=@__FILE__;root="/",port=8000,build="")
    makefile = normpath(joinpath(modulepath,"..","..","docs","make.jl"))
    if isequal(basename(makefile),"make.jl")
        folder = dirname(makefile)
        if isdir(joinpath(folder,build))
            buildfolder = joinpath(folder,build)
            Pages.servefile(joinpath(buildfolder,"index.html"),root)
            Pages.servefile(joinpath(buildfolder,"search_index.js"))
            Pages.servefolder(joinpath(buildfolder,"assets"),"/assets")
            Pages.servefolder(joinpath(buildfolder,"search"),"/search")
        end
        @async Pages.start(port)
    else
        @warn "$(makefile) does not point to a valid make.jl"
    end
end

end # module


text:
@testset "1.2.1.1 (a+b x+c x^2)^p" begin
    (a, b, c, k, n, p, x, ) = @variables a b c k n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p*when*a=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x+c*x^2)^(p/2)=#

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

    @test_int [(b*x + c*x^2)^(7/2), x, 6, -((35*b^6*(b + 2*c*x)*sqrt(b*x + c*x^2))/(16384*c^4)) + (35*b^4*(b + 2*c*x)*(b*x + c*x^2)^(3/2))/(6144*c^3) - (7*b^2*(b + 2*c*x)*(b*x + c*x^2)^(5/2))/(384*c^2) + ((b + 2*c*x)*(b*x + c*x^2)^(7/2))/(16*c) + (35*b^8*atanh((sqrt(c)*x)/sqrt(b*x + c*x^2)))/(16384*c^(9/2))]

    @test_int [(3*im*x + 4*x^2)^(7/2), x, 6, (25515*(3*im + 8*x)*sqrt(3*im*x + 4*x^2))/4194304 + (945*(3*im + 8*x)*(3*im*x + 4*x^2)^(3/2))/131072 + (21*(3*im + 8*x)*(3*im*x + 4*x^2)^(5/2))/2048 + (1/64)*(3*im + 8*x)*(3*im*x + 4*x^2)^(7/2) + (229635*im*asin(1 - (8*im*x)/3))/16777216]
    @test_int [(3*im*x + 4*x^2)^(5/2), x, 5, (405*(3*im + 8*x)*sqrt(3*im*x + 4*x^2))/32768 + (15*(3*im + 8*x)*(3*im*x + 4*x^2)^(3/2))/1024 + (1/48)*(3*im + 8*x)*(3*im*x + 4*x^2)^(5/2) + (3645*im*asin(1 - (8*im*x)/3))/131072]
    @test_int [(3*im*x + 4*x^2)^(3/2), x, 4, (27*(3*im + 8*x)*sqrt(3*im*x + 4*x^2))/1024 + (1/32)*(3*im + 8*x)*(3*im*x + 4*x^2)^(3/2) + (243*im*asin(1 - (8*im*x)/3))/4096]
    @test_int [(3*im*x + 4*x^2)^(1/2), x, 3, (1/16)*(3*im + 8*x)*sqrt(3*im*x + 4*x^2) + (9/64)*im*asin(1 - (8*im*x)/3)]

    @test_int [(3*x - 4*x^2)^(7/2), x, 6, -((25515*(3 - 8*x)*sqrt(3*x - 4*x^2))/4194304) - (945*(3 - 8*x)*(3*x - 4*x^2)^(3/2))/131072 - (21*(3 - 8*x)*(3*x - 4*x^2)^(5/2))/2048 - (1/64)*(3 - 8*x)*(3*x - 4*x^2)^(7/2) - (229635*asin(1 - (8*x)/3))/16777216]
    @test_int [(3*x - 4*x^2)^(5/2), x, 5, -((405*(3 - 8*x)*sqrt(3*x - 4*x^2))/32768) - (15*(3 - 8*x)*(3*x - 4*x^2)^(3/2))/1024 - (1/48)*(3 - 8*x)*(3*x - 4*x^2)^(5/2) - (3645*asin(1 - (8*x)/3))/131072]
    @test_int [(3*x - 4*x^2)^(3/2), x, 4, -((27*(3 - 8*x)*sqrt(3*x - 4*x^2))/1024) - (1/32)*(3 - 8*x)*(3*x - 4*x^2)^(3/2) - (243*asin(1 - (8*x)/3))/4096]
    @test_int [(3*x - 4*x^2)^(1/2), x, 3, (-(1/16))*(3 - 8*x)*sqrt(3*x - 4*x^2) - (9/64)*asin(1 - (8*x)/3)]

    @test_int [sqrt(6*x - x^2), x, 3, (-(1/2))*(3 - x)*sqrt(6*x - x^2) - (9/2)*asin(1 - x/3)]
    @test_int [sqrt(5*x - 9*x^2), x, 3, (-(1/36))*(5 - 18*x)*sqrt(5*x - 9*x^2) - (25/216)*asin(1 - (18*x)/5)]
    @test_int [(x - x^2)^(3/2), x, 4, (-(3/64))*(1 - 2*x)*sqrt(x - x^2) - (1/8)*(1 - 2*x)*(x - x^2)^(3/2) - (3/128)*asin(1 - 2*x)]

    @test_int [sqrt(4*x + x^2), x, 3, (1/2)*(2 + x)*sqrt(4*x + x^2) - 4*atanh(x/sqrt(4*x + x^2))]
    @test_int [sqrt(-8*x + x^2), x, 3, (-(1/2))*(4 - x)*sqrt(-8*x + x^2) - 16*atanh(x/sqrt(-8*x + x^2))]
    @test_int [sqrt(-x + x^2), x, 3, (-(1/4))*(1 - 2*x)*sqrt(-x + x^2) - (1/4)*atanh(x/sqrt(-x + x^2))]

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

    @test_int [1/(b*x + c*x^2)^(7/2), x, 3, -((2*(b + 2*c*x))/(5*b^2*(b*x + c*x^2)^(5/2))) + (32*c*(b + 2*c*x))/(15*b^4*(b*x + c*x^2)^(3/2)) - (256*c^2*(b + 2*c*x))/(15*b^6*sqrt(b*x + c*x^2))]

    @test_int [1/(3*im*x + 4*x^2)^(1/2), x, 2, (1/2)*im*asin(1 - (8*im*x)/3)]
    @test_int [1/(3*im*x + 4*x^2)^(3/2), x, 1, (2*(3*im + 8*x))/(9*sqrt(3*im*x + 4*x^2))]
    @test_int [1/(3*im*x + 4*x^2)^(5/2), x, 2, (2*(3*im + 8*x))/(27*(3*im*x + 4*x^2)^(3/2)) + (64*(3*im + 8*x))/(243*sqrt(3*im*x + 4*x^2))]
    @test_int [1/(3*im*x + 4*x^2)^(7/2), x, 3, (2*(3*im + 8*x))/(45*(3*im*x + 4*x^2)^(5/2)) + (128*(3*im + 8*x))/(1215*(3*im*x + 4*x^2)^(3/2)) + (4096*(3*im + 8*x))/(10935*sqrt(3*im*x + 4*x^2))]

    @test_int [1/(3*x - 4*x^2)^(1/2), x, 2, (-(1/2))*asin(1 - (8*x)/3)]
    @test_int [1/(3*x - 4*x^2)^(3/2), x, 1, -((2*(3 - 8*x))/(9*sqrt(3*x - 4*x^2)))]
    @test_int [1/(3*x - 4*x^2)^(5/2), x, 2, -((2*(3 - 8*x))/(27*(3*x - 4*x^2)^(3/2))) - (64*(3 - 8*x))/(243*sqrt(3*x - 4*x^2))]
    @test_int [1/(3*x - 4*x^2)^(7/2), x, 3, -((2*(3 - 8*x))/(45*(3*x - 4*x^2)^(5/2))) - (128*(3 - 8*x))/(1215*(3*x - 4*x^2)^(3/2)) - (4096*(3 - 8*x))/(10935*sqrt(3*x - 4*x^2))]

    @test_int [1/sqrt(b*x - b^2*x^2), x, 2, -(asin(1 - 2*b*x)/b)]
    @test_int [1/sqrt(b*x + b^2*x^2), x, 2, (2*atanh((b*x)/sqrt(b*x + b^2*x^2)))/b]

    @test_int [1/sqrt(6*x - x^2), x, 2, -asin(1 - x/3)]
    @test_int [1/sqrt(4*x + x^2), x, 2, 2*atanh(x/sqrt(4*x + x^2))]
    @test_int [1/sqrt(-2*x + x^2), x, 2, 2*atanh(x/sqrt(-2*x + x^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x+c*x^2)^(p/3)=#

    @test_int [(b*x + c*x^2)^(4/3), x, 6, (3*(-((c*x*(b + c*x))/b^2))^(1/3)*(b + 2*c*x)*(b*x + c*x^2)^(4/3))/(55*c*(-((c*(b*x + c*x^2))/b^2))^(4/3)) + (3*(-((c*x*(b + c*x))/b^2))^(4/3)*(b + 2*c*x)*(b*x + c*x^2)^(4/3))/(22*c*(-((c*(b*x + c*x^2))/b^2))^(4/3)) + (2^(1/3)*3^(3/4)*sqrt(2 - sqrt(3))*b^2*(b*x + c*x^2)^(4/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(55*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(4/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [(b*x + c*x^2)^(1/3), x, 5, (3*(-((c*x*(b + c*x))/b^2))^(1/3)*(b + 2*c*x)*(b*x + c*x^2)^(1/3))/(10*c*(-((c*(b*x + c*x^2))/b^2))^(1/3)) + (3^(3/4)*sqrt(2 - sqrt(3))*b^2*(b*x + c*x^2)^(1/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(5*2^(2/3)*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(1/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(2/3), x, 4, (2^(1/3)*3^(3/4)*sqrt(2 - sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(2/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(c*(b + 2*c*x)*(b*x + c*x^2)^(2/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(5/3), x, 5, (3*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(5/3))/(2*c*(-((c*x*(b + c*x))/b^2))^(2/3)*(b*x + c*x^2)^(5/3)) + (2^(1/3)*3^(3/4)*sqrt(2 - sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(5/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(c*(b + 2*c*x)*(b*x + c*x^2)^(5/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(8/3), x, 6, (3*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(8/3))/(5*c*(-((c*x*(b + c*x))/b^2))^(5/3)*(b*x + c*x^2)^(8/3)) + (21*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(8/3))/(5*c*(-((c*x*(b + c*x))/b^2))^(2/3)*(b*x + c*x^2)^(8/3)) + (14*2^(1/3)*3^(3/4)*sqrt(2 - sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(8/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(5*c*(b + 2*c*x)*(b*x + c*x^2)^(8/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]

    @test_int [(b*x + c*x^2)^(5/3), x, 8, (15*(-((c*x*(b + c*x))/b^2))^(2/3)*(b + 2*c*x)*(b*x + c*x^2)^(5/3))/(364*c*(-((c*(b*x + c*x^2))/b^2))^(5/3)) + (3*(-((c*x*(b + c*x))/b^2))^(5/3)*(b + 2*c*x)*(b*x + c*x^2)^(5/3))/(26*c*(-((c*(b*x + c*x^2))/b^2))^(5/3)) - (15*(b + 2*c*x)*(b*x + c*x^2)^(5/3))/(182*2^(1/3)*c*(-((c*(b*x + c*x^2))/b^2))^(5/3)*(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))) - (15*3^(1/4)*sqrt(2 + sqrt(3))*b^2*(b*x + c*x^2)^(5/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.E(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(364*2^(1/3)*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(5/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2))) + (5*3^(3/4)*b^2*(b*x + c*x^2)^(5/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(91*2^(5/6)*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(5/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [(b*x + c*x^2)^(2/3), x, 7, (3*(-((c*x*(b + c*x))/b^2))^(2/3)*(b + 2*c*x)*(b*x + c*x^2)^(2/3))/(14*c*(-((c*(b*x + c*x^2))/b^2))^(2/3)) - (3*(b + 2*c*x)*(b*x + c*x^2)^(2/3))/(7*2^(1/3)*c*(-((c*(b*x + c*x^2))/b^2))^(2/3)*(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))) - (3*3^(1/4)*sqrt(2 + sqrt(3))*b^2*(b*x + c*x^2)^(2/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.E(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(14*2^(1/3)*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(2/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2))) + (2^(1/6)*3^(3/4)*b^2*(b*x + c*x^2)^(2/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(7*c*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(2/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(1/3), x, 6, -((3*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(1/3))/(2^(1/3)*c*(b*x + c*x^2)^(1/3)*(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3)))) - (3*3^(1/4)*sqrt(2 + sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(1/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.E(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(2*2^(1/3)*c*(b + 2*c*x)*(b*x + c*x^2)^(1/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2))) + (2^(1/6)*3^(3/4)*b^2*(-((c*(b*x + c*x^2))/b^2))^(1/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(c*(b + 2*c*x)*(b*x + c*x^2)^(1/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(4/3), x, 7, (3*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(4/3))/(c*(-((c*x*(b + c*x))/b^2))^(1/3)*(b*x + c*x^2)^(4/3)) + (3*2^(2/3)*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(4/3))/(c*(b*x + c*x^2)^(4/3)*(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))) + (3*3^(1/4)*sqrt(2 + sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(4/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.E(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(2^(1/3)*c*(b + 2*c*x)*(b*x + c*x^2)^(4/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2))) - (2*2^(1/6)*3^(3/4)*b^2*(-((c*(b*x + c*x^2))/b^2))^(4/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(c*(b + 2*c*x)*(b*x + c*x^2)^(4/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]
    @test_int [1/(b*x + c*x^2)^(7/3), x, 8, (3*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(7/3))/(4*c*(-((c*x*(b + c*x))/b^2))^(4/3)*(b*x + c*x^2)^(7/3)) + (15*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(7/3))/(2*c*(-((c*x*(b + c*x))/b^2))^(1/3)*(b*x + c*x^2)^(7/3)) + (15*(b + 2*c*x)*(-((c*(b*x + c*x^2))/b^2))^(7/3))/(2^(1/3)*c*(b*x + c*x^2)^(7/3)*(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))) + (15*3^(1/4)*sqrt(2 + sqrt(3))*b^2*(-((c*(b*x + c*x^2))/b^2))^(7/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.E(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(2*2^(1/3)*c*(b + 2*c*x)*(b*x + c*x^2)^(7/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2))) - (5*2^(1/6)*3^(3/4)*b^2*(-((c*(b*x + c*x^2))/b^2))^(7/3)*(1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))*sqrt((1 + 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3) + 2*2^(1/3)*(-((c*x*(b + c*x))/b^2))^(2/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)*Elliptic.F(asin((1 + sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))), -7 + 4*sqrt(3)))/(c*(b + 2*c*x)*(b*x + c*x^2)^(7/3)*sqrt(-((1 - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))/(1 - sqrt(3) - 2^(2/3)*(-((c*x*(b + c*x))/b^2))^(1/3))^2)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x+c*x^2)^(p/4)=#

    @test_int [(b*x + c*x^2)^(5/4), x, 5, -((5*b^2*(b + 2*c*x)*(b*x + c*x^2)^(1/4))/(84*c^2)) + ((b + 2*c*x)*(b*x + c*x^2)^(5/4))/(7*c) + (5*b^5*(-((c*(b*x + c*x^2))/b^2))^(3/4)*Elliptic.F((1/2)*asin(1 + (2*c*x)/b), 2))/(84*sqrt(2)*c^3*(b*x + c*x^2)^(3/4))]
    @test_int [(b*x + c*x^2)^(3/4), x, 4, ((b + 2*c*x)*(b*x + c*x^2)^(3/4))/(5*c) - (3*b^3*(-((c*(b*x + c*x^2))/b^2))^(1/4)*Elliptic.E((1/2)*asin(1 + (2*c*x)/b), 2))/(10*sqrt(2)*c^2*(b*x + c*x^2)^(1/4))]
    @test_int [(b*x + c*x^2)^(1/4), x, 4, ((b + 2*c*x)*(b*x + c*x^2)^(1/4))/(3*c) - (b^3*(-((c*(b*x + c*x^2))/b^2))^(3/4)*Elliptic.F((1/2)*asin(1 + (2*c*x)/b), 2))/(3*sqrt(2)*c^2*(b*x + c*x^2)^(3/4))]
    @test_int [1/(b*x + c*x^2)^(1/4), x, 3, (sqrt(2)*b*(-((c*(b*x + c*x^2))/b^2))^(1/4)*Elliptic.E((1/2)*asin(1 + (2*c*x)/b), 2))/(c*(b*x + c*x^2)^(1/4))]
    @test_int [1/(b*x + c*x^2)^(3/4), x, 3, (2*sqrt(2)*b*(-((c*(b*x + c*x^2))/b^2))^(3/4)*Elliptic.F((1/2)*asin(1 + (2*c*x)/b), 2))/(c*(b*x + c*x^2)^(3/4))]
    @test_int [1/(b*x + c*x^2)^(5/4), x, 4, -((4*(b + 2*c*x))/(b^2*(b*x + c*x^2)^(1/4))) + (4*sqrt(2)*(-((c*(b*x + c*x^2))/b^2))^(1/4)*Elliptic.E((1/2)*asin(1 + (2*c*x)/b), 2))/(b*(b*x + c*x^2)^(1/4))]
    @test_int [1/(b*x + c*x^2)^(9/4), x, 5, -((4*(b + 2*c*x))/(5*b^2*(b*x + c*x^2)^(5/4))) + (48*c*(b + 2*c*x))/(5*b^4*(b*x + c*x^2)^(1/4)) - (48*sqrt(2)*c*(-((c*(b*x + c*x^2))/b^2))^(1/4)*Elliptic.E((1/2)*asin(1 + (2*c*x)/b), 2))/(5*b^3*(b*x + c*x^2)^(1/4))]
    @test_int [1/(b*x + c*x^2)^(13/4), x, 6, -((4*(b + 2*c*x))/(9*b^2*(b*x + c*x^2)^(9/4))) + (112*c*(b + 2*c*x))/(45*b^4*(b*x + c*x^2)^(5/4)) - (448*c^2*(b + 2*c*x))/(15*b^6*(b*x + c*x^2)^(1/4)) + (448*sqrt(2)*c^2*(-((c*(b*x + c*x^2))/b^2))^(1/4)*Elliptic.E((1/2)*asin(1 + (2*c*x)/b), 2))/(15*b^5*(b*x + c*x^2)^(1/4))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x+c*x^2)^p*when*p*symbolic=#

    @test_int [(b*x + c*x^2)^p, x, 1, -(((-((c*x)/b))^(-1 - p)*(b*x + c*x^2)^(1 + p)*HypergeometricFunctions._₂F₁(-p, 1 + p, 2 + p, (b + c*x)/b))/(b*(1 + p)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p*when*b=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+c*x^2)^p=#

    @test_int [(a + c*x^2)^4, x, 2, a^4*x + (4/3)*a^3*c*x^3 + (6/5)*a^2*c^2*x^5 + (4/7)*a*c^3*x^7 + (c^4*x^9)/9]
    @test_int [(a + c*x^2)^3, x, 2, a^3*x + a^2*c*x^3 + (3/5)*a*c^2*x^5 + (c^3*x^7)/7]
    @test_int [(a + c*x^2)^2, x, 2, a^2*x + (2/3)*a*c*x^3 + (c^2*x^5)/5]
    @test_int [(a + c*x^2)^1, x, 1, a*x + (c*x^3)/3]
    @test_int [1/(a + c*x^2)^1, x, 1, atan((sqrt(c)*x)/sqrt(a))/(sqrt(a)*sqrt(c))]
    @test_int [1/(a + c*x^2)^2, x, 2, x/(2*a*(a + c*x^2)) + atan((sqrt(c)*x)/sqrt(a))/(2*a^(3/2)*sqrt(c))]
    @test_int [1/(a + c*x^2)^3, x, 3, x/(4*a*(a + c*x^2)^2) + (3*x)/(8*a^2*(a + c*x^2)) + (3*atan((sqrt(c)*x)/sqrt(a)))/(8*a^(5/2)*sqrt(c))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+c*x^2)^(p/2)=#

    @test_int [(a + c*x^2)^(5/2), x, 5, (5*a^2*x*sqrt(a + c*x^2))/16 + (5*a*x*(a + c*x^2)^(3/2))/24 + (x*(a + c*x^2)^(5/2))/6 + (5*a^3*atanh((sqrt(c)*x)/sqrt(a + c*x^2)))/(16*sqrt(c))]
    @test_int [(a + c*x^2)^(3/2), x, 4, (3*a*x*sqrt(a + c*x^2))/8 + (x*(a + c*x^2)^(3/2))/4 + (3*a^2*atanh((sqrt(c)*x)/sqrt(a + c*x^2)))/(8*sqrt(c))]
    @test_int [(a + c*x^2)^(1/2), x, 3, (1/2)*x*sqrt(a + c*x^2) + (a*atanh((sqrt(c)*x)/sqrt(a + c*x^2)))/(2*sqrt(c))]
    @test_int [1/(a + c*x^2)^(1/2), x, 2, atanh((sqrt(c)*x)/sqrt(a + c*x^2))/sqrt(c)]
    @test_int [1/(a + c*x^2)^(3/2), x, 1, x/(a*sqrt(a + c*x^2))]
    @test_int [1/(a + c*x^2)^(5/2), x, 2, x/(3*a*(a + c*x^2)^(3/2)) + (2*x)/(3*a^2*sqrt(a + c*x^2))]
    @test_int [1/(a + c*x^2)^(7/2), x, 3, x/(5*a*(a + c*x^2)^(5/2)) + (4*x)/(15*a^2*(a + c*x^2)^(3/2)) + (8*x)/(15*a^3*sqrt(a + c*x^2))]
    @test_int [1/(a + c*x^2)^(9/2), x, 4, x/(7*a*(a + c*x^2)^(7/2)) + (6*x)/(35*a^2*(a + c*x^2)^(5/2)) + (8*x)/(35*a^3*(a + c*x^2)^(3/2)) + (16*x)/(35*a^4*sqrt(a + c*x^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p*when*b^2-4*a*c=0=#

    @test_int [(4 + 12*x + 9*x^2)^(3/2), x, 1, (1/12)*(2 + 3*x)*(4 + 12*x + 9*x^2)^(3/2)]
    @test_int [(4 + 12*x + 9*x^2)^(1/2), x, 1, (1/6)*(2 + 3*x)*sqrt(4 + 12*x + 9*x^2)]
    @test_int [1/(4 + 12*x + 9*x^2)^(1/2), x, 2, ((2 + 3*x)*log(2 + 3*x))/(3*sqrt(4 + 12*x + 9*x^2))]
    @test_int [1/(4 + 12*x + 9*x^2)^(3/2), x, 1, -(1/(6*(2 + 3*x)*sqrt(4 + 12*x + 9*x^2)))]

    @test_int [sqrt(4 - 12*x + 9*x^2), x, 1, (-(1/6))*(2 - 3*x)*sqrt(4 - 12*x + 9*x^2)]
    @test_int [1/sqrt(4 - 12*x + 9*x^2), x, 2, -(((2 - 3*x)*log(2 - 3*x))/(3*sqrt(4 - 12*x + 9*x^2)))]

    @test_int [sqrt(-4 + 12*x - 9*x^2), x, 1, (-(1/6))*(2 - 3*x)*sqrt(-4 + 12*x - 9*x^2)]
    @test_int [1/sqrt(-4 + 12*x - 9*x^2), x, 2, -(((2 - 3*x)*log(2 - 3*x))/(3*sqrt(-4 + 12*x - 9*x^2)))]

    @test_int [sqrt(-4 - 12*x - 9*x^2), x, 1, (1/6)*(2 + 3*x)*sqrt(-4 - 12*x - 9*x^2)]
    @test_int [1/sqrt(-4 - 12*x - 9*x^2), x, 2, ((2 + 3*x)*log(2 + 3*x))/(3*sqrt(-4 - 12*x - 9*x^2))]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p=#

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

    @test_int [((-1 + b^2)/(4*c) + b*x + c*x^2)^5, x, 3, (1 - b - 2*c*x)^6/(384*c^6) - (5*(1 - b - 2*c*x)^7)/(896*c^6) + (5*(1 - b - 2*c*x)^8)/(1024*c^6) - (5*(1 - b - 2*c*x)^9)/(2304*c^6) + (1 - b - 2*c*x)^10/(2048*c^6) - (1 - b - 2*c*x)^11/(22528*c^6)]
    @test_int [((-4 + b^2)/(4*c) + b*x + c*x^2)^5, x, 3, (2 - b - 2*c*x)^6/(12*c^6) - (5*(2 - b - 2*c*x)^7)/(56*c^6) + (5*(2 - b - 2*c*x)^8)/(128*c^6) - (5*(2 - b - 2*c*x)^9)/(576*c^6) + (2 - b - 2*c*x)^10/(1024*c^6) - (2 - b - 2*c*x)^11/(22528*c^6)]
    @test_int [((-9 + b^2)/(4*c) + b*x + c*x^2)^5, x, 3, (81*(3 - b - 2*c*x)^6)/(128*c^6) - (405*(3 - b - 2*c*x)^7)/(896*c^6) + (135*(3 - b - 2*c*x)^8)/(1024*c^6) - (5*(3 - b - 2*c*x)^9)/(256*c^6) + (3*(3 - b - 2*c*x)^10)/(2048*c^6) - (3 - b - 2*c*x)^11/(22528*c^6)]
    @test_int [((-16 + b^2)/(4*c) + b*x + c*x^2)^5, x, 3, (8*(4 - b - 2*c*x)^6)/(3*c^6) - (10*(4 - b - 2*c*x)^7)/(7*c^6) + (5*(4 - b - 2*c*x)^8)/(16*c^6) - (5*(4 - b - 2*c*x)^9)/(144*c^6) + (4 - b - 2*c*x)^10/(512*c^6) - (4 - b - 2*c*x)^11/(22528*c^6)]

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

    @test_int [1/(2 + 4*x + 3*x^2), x, 2, atan((2 + 3*x)/sqrt(2))/sqrt(2)]
    @test_int [1/(4 - 2*sqrt(3)*x + x^2), x, 2, -atan(sqrt(3) - x)]
    @test_int [1/(2 + 4*x - 3*x^2), x, 2, -(atanh((2 - 3*x)/sqrt(10))/sqrt(10))]
    @test_int [1/(2 + 5*x + 3*x^2), x, 3, -log(1 + x) + log(2 + 3*x)]
    @test_int [1/(2 + 5*x - 3*x^2), x, 3, (-(1/7))*log(2 - x) + (1/7)*log(1 + 3*x)]
    @test_int [1/(3 + 4*x + x^2), x, 3, -atanh(2 + x), (1/2)*log(1 + x) - (1/2)*log(3 + x)]
    @test_int [1/(1 + pi*x + 2*x^2), x, 2, (-2*atanh((pi + 4*x)/sqrt(-8 + pi^2)))/sqrt(-8 + pi^2)]
    @test_int [1/(1 + pi*x - 2*x^2), x, 2, (-2*atanh((pi - 4*x)/sqrt(8 + pi^2)))/sqrt(8 + pi^2)]
    @test_int [1/(1 + pi*x + 3*x^2), x, 2, (2*atan((pi + 6*x)/sqrt(12 - pi^2)))/sqrt(12 - pi^2)]
    @test_int [1/(1 + pi*x - 3*x^2), x, 2, (-2*atanh((pi - 6*x)/sqrt(12 + pi^2)))/sqrt(12 + pi^2)]
    @test_int [1/(a + c*x + b*x^2), x, 2, (2*atan((c + 2*b*x)/sqrt(4*a*b - c^2)))/sqrt(4*a*b - c^2)]
    @test_int [1/(b + 2*a*x + b*x^2), x, 2, -(atanh((a + b*x)/sqrt(a^2 - b^2))/sqrt(a^2 - b^2))]
    @test_int [1/(b + 2*a*x - b*x^2), x, 2, -(atanh((a - b*x)/sqrt(a^2 + b^2))/sqrt(a^2 + b^2))]

    @test_int [1/(2 + 4*x + 3*x^2)^2, x, 3, (2 + 3*x)/(4*(2 + 4*x + 3*x^2)) + (3*atan((2 + 3*x)/sqrt(2)))/(4*sqrt(2))]
    @test_int [1/(2 + 4*x - 3*x^2)^2, x, 3, -((2 - 3*x)/(20*(2 + 4*x - 3*x^2))) - (3*atanh((2 - 3*x)/sqrt(10)))/(20*sqrt(10))]
    @test_int [1/(2 + 5*x + 3*x^2)^2, x, 4, -((5 + 6*x)/(2 + 5*x + 3*x^2)) + 6*log(1 + x) - 6*log(2 + 3*x)]
    @test_int [1/(2 + 5*x - 3*x^2)^2, x, 4, -((5 - 6*x)/(49*(2 + 5*x - 3*x^2))) - (6/343)*log(2 - x) + (6/343)*log(1 + 3*x)]
    @test_int [1/(a + c*x + b*x^2)^2, x, 3, (c + 2*b*x)/((4*a*b - c^2)*(a + c*x + b*x^2)) + (4*b*atan((c + 2*b*x)/sqrt(4*a*b - c^2)))/(4*a*b - c^2)^(3/2)]
    @test_int [1/(b + 2*a*x + b*x^2)^2, x, 3, -((a + b*x)/(2*(a^2 - b^2)*(b + 2*a*x + b*x^2))) + (b*atanh((a + b*x)/sqrt(a^2 - b^2)))/(2*(a^2 - b^2)^(3/2))]
    @test_int [1/(b + 2*a*x - b*x^2)^2, x, 3, -((a - b*x)/(2*(a^2 + b^2)*(b + 2*a*x - b*x^2))) - (b*atanh((a - b*x)/sqrt(a^2 + b^2)))/(2*(a^2 + b^2)^(3/2))]

    @test_int [1/((a/b)^(2/n) + x^2 - 2*(a/b)^(1/n)*x*cos((pi - 2*k*pi)/n)), x, 2, (-(a/b)^(-n^(-1)))*atan(cot((pi - 2*k*pi)/n) - (x*csc((pi - 2*k*pi)/n))/(a/b)^n^(-1))*csc((pi - 2*k*pi)/n)]

    @test_int [1/(a*b + sqrt(b^2 - 4*a*b^3)*x - b^2*x^2), x, 3, 2*(atanh((-sqrt(b^2 - 4*a*b^3) + 2*b^2*x)/b)/b), -(log(b + sqrt(b^2 - 4*a*b^3) - 2*b^2*x)/b) + log(b - sqrt(b^2 - 4*a*b^3) + 2*b^2*x)/b]
    @test_int [1/(a*b - sqrt(b^2 - 4*a*b^3)*x - b^2*x^2), x, 3, 2*(atanh((sqrt(b^2 - 4*a*b^3) + 2*b^2*x)/b)/b), -(log(b - sqrt(b^2 - 4*a*b^3) - 2*b^2*x)/b) + log(b + sqrt(b^2 - 4*a*b^3) + 2*b^2*x)/b]

    @test_int [1/(1 + x^2 + 2*x*cos(1/7)), x, 2, atan((x + cos(1/7))*csc(1/7))*csc(1/7)]
    @test_int [1/(1 + x^2 + 2*x*cos(pi/7)), x, 2, atan(cot(pi/7) + x*csc(pi/7))*csc(pi/7)]

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

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

    @test_int [sqrt(5 - 6*x + 9*x^2), x, 3, (-(1/6))*(1 - 3*x)*sqrt(5 - 6*x + 9*x^2) + (2/3)*asinh((1/2)*(-1 + 3*x))]
    @test_int [sqrt(3 - 4*x - 4*x^2), x, 3, (1/4)*(1 + 2*x)*sqrt(3 - 4*x - 4*x^2) + asin(1/2 + x)]
    @test_int [sqrt(-8 + 6*x + 9*x^2), x, 3, (1/6)*(1 + 3*x)*sqrt(-8 + 6*x + 9*x^2) - (3/2)*atanh((1 + 3*x)/sqrt(-8 + 6*x + 9*x^2))]
    @test_int [sqrt(2 + 4*x + 3*x^2), x, 3, (1/6)*(2 + 3*x)*sqrt(2 + 4*x + 3*x^2) + asinh((2 + 3*x)/sqrt(2))/(3*sqrt(3))]
    @test_int [sqrt(2 + 4*x - 3*x^2), x, 3, (-(1/6))*(2 - 3*x)*sqrt(2 + 4*x - 3*x^2) - (5*asin((2 - 3*x)/sqrt(10)))/(3*sqrt(3))]
    @test_int [sqrt(2 + 5*x + 3*x^2), x, 3, (1/12)*(5 + 6*x)*sqrt(2 + 5*x + 3*x^2) - atanh((5 + 6*x)/(2*sqrt(3)*sqrt(2 + 5*x + 3*x^2)))/(24*sqrt(3))]
    @test_int [sqrt(2 + 5*x - 3*x^2), x, 3, (-(1/12))*(5 - 6*x)*sqrt(2 + 5*x - 3*x^2) - (49*asin((1/7)*(5 - 6*x)))/(24*sqrt(3))]
    @test_int [sqrt(-2 + 4*x + 3*x^2), x, 3, (1/6)*(2 + 3*x)*sqrt(-2 + 4*x + 3*x^2) - (5*atanh((2 + 3*x)/(sqrt(3)*sqrt(-2 + 4*x + 3*x^2))))/(3*sqrt(3))]
    @test_int [sqrt(-2 + 4*x - 3*x^2), x, 3, (-(1/6))*(2 - 3*x)*sqrt(-2 + 4*x - 3*x^2) + atan((2 - 3*x)/(sqrt(3)*sqrt(-2 + 4*x - 3*x^2)))/(3*sqrt(3))]
    @test_int [sqrt(-2 + 5*x + 3*x^2), x, 3, (1/12)*(5 + 6*x)*sqrt(-2 + 5*x + 3*x^2) - (49*atanh((5 + 6*x)/(2*sqrt(3)*sqrt(-2 + 5*x + 3*x^2))))/(24*sqrt(3))]
    @test_int [sqrt(-2 + 5*x - 3*x^2), x, 3, (-(1/12))*(5 - 6*x)*sqrt(-2 + 5*x - 3*x^2) - asin(5 - 6*x)/(24*sqrt(3))]

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

    @test_int [1/sqrt(5 - 6*x + 9*x^2), x, 2, asinh((-1 + 3*x)/2)/3]
    @test_int [1/sqrt(3 - 4*x - 4*x^2), x, 2, asin(1/2 + x)/2]
    @test_int [1/sqrt(-8 + 6*x + 9*x^2), x, 2, (1/3)*atanh((1 + 3*x)/sqrt(-8 + 6*x + 9*x^2))]
    @test_int [1/sqrt(2 + 4*x + 3*x^2), x, 2, asinh((2 + 3*x)/sqrt(2))/sqrt(3)]
    @test_int [1/sqrt(2 + 4*x - 3*x^2), x, 2, -(asin((2 - 3*x)/sqrt(10))/sqrt(3))]
    @test_int [1/sqrt(2 + 5*x + 3*x^2), x, 2, atanh((5 + 6*x)/(2*sqrt(3)*sqrt(2 + 5*x + 3*x^2)))/sqrt(3)]
    @test_int [1/sqrt(2 + 5*x - 3*x^2), x, 2, -(asin((5 - 6*x)/7)/sqrt(3))]
    @test_int [1/sqrt(-2 + 4*x + 3*x^2), x, 2, atanh((2 + 3*x)/(sqrt(3)*sqrt(-2 + 4*x + 3*x^2)))/sqrt(3)]
    @test_int [1/sqrt(-2 + 4*x - 3*x^2), x, 2, -(atan((2 - 3*x)/(sqrt(3)*sqrt(-2 + 4*x - 3*x^2)))/sqrt(3))]
    @test_int [1/sqrt(-2 + 5*x + 3*x^2), x, 2, atanh((5 + 6*x)/(2*sqrt(3)*sqrt(-2 + 5*x + 3*x^2)))/sqrt(3)]
    @test_int [1/sqrt(-2 + 5*x - 3*x^2), x, 2, -(asin(5 - 6*x)/sqrt(3))]

    @test_int [1/sqrt((4*c + b^2)/(4*c) + b*x + c*x^2), x, 2, asinh((b + 2*c*x)/(2*sqrt(c)))/sqrt(c)]
    @test_int [1/sqrt((4*c - b^2)/(4*c) + b*x - c*x^2), x, 2, -(asin((b - 2*c*x)/(2*sqrt(c)))/sqrt(c))]
    @test_int [1/sqrt((c - b^2)/(4*c) + b*x - c*x^2), x, 2, -(asin((b - 2*c*x)/sqrt(c))/sqrt(c))]

    @test_int [1/(2 + 3*x + x^2)^(3/2), x, 1, (-2*(3 + 2*x))/sqrt(2 + 3*x + x^2)]
    @test_int [1/(27 - 24*x + 4*x^2)^(3/2), x, 1, (3 - x)/(9*sqrt(27 - 24*x + 4*x^2))]

    @test_int [x/(5 - 4*x - x^2)^(3/2), x, 1, (5 - 2*x)/(9*sqrt(5 - 4*x - x^2))]

    @test_int [1/(5 - 4*x - x^2)^(5/2), x, 2, (2 + x)/(27*(5 - 4*x - x^2)^(3/2)) + (2*(2 + x))/(243*sqrt(5 - 4*x - x^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*x+c*x^2)^p*when*p*symbolic=#

    @test_int [(a + b*x + c*x^2)^p, x, 1, -((2^(1 + p)*(-((b - sqrt(b^2 - 4*a*c) + 2*c*x)/sqrt(b^2 - 4*a*c)))^(-1 - p)*(a + b*x + c*x^2)^(1 + p)*HypergeometricFunctions._₂F₁(-p, 1 + p, 2 + p, (b + sqrt(b^2 - 4*a*c) + 2*c*x)/(2*sqrt(b^2 - 4*a*c))))/(sqrt(b^2 - 4*a*c)*(1 + p)))]

    @test_int [(3 + 4*x + 5*x^2)^p, x, 2, 5^(-1 - p)*11^p*(2 + 5*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (-(1/11))*(2 + 5*x)^2)]
    @test_int [(3 + 4*x + 4*x^2)^p, x, 2, 2^(-1 + p)*(1 + 2*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (-(1/2))*(1 + 2*x)^2)]
    @test_int [(3 + 4*x + 3*x^2)^p, x, 2, 3^(-1 - p)*5^p*(2 + 3*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (-(1/5))*(2 + 3*x)^2)]
    @test_int [(3 + 4*x + 2*x^2)^p, x, 2, (1 + x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, -2*(1 + x)^2)]
    @test_int [(3 + 4*x + 1*x^2)^p, x, 1, -((2^(1 + 2*p)*(-2 - 2*x)^(-1 - p)*(3 + 4*x + x^2)^(1 + p)*HypergeometricFunctions._₂F₁(-p, 1 + p, 2 + p, (3 + x)/2))/(1 + p))]
    @test_int [(3 + 4*x + 0*x^2)^p, x, 1, (3 + 4*x)^(1 + p)/(4*(1 + p))]
    @test_int [(3 + 4*x - 1*x^2)^p, x, 2, (-7^p)*(2 - x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (1/7)*(2 - x)^2)]
    @test_int [(3 + 4*x - 2*x^2)^p, x, 2, (-5^p)*(1 - x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (2/5)*(1 - x)^2)]
    @test_int [(3 + 4*x - 3*x^2)^p, x, 2, (-3^(-1 - p))*13^p*(2 - 3*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (1/13)*(2 - 3*x)^2)]
    @test_int [(3 + 4*x - 4*x^2)^p, x, 2, (-2^(-1 + 2*p))*(1 - 2*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (1/4)*(1 - 2*x)^2)]
    @test_int [(3 + 4*x - 5*x^2)^p, x, 2, (-5^(-1 - p))*19^p*(2 - 5*x)*HypergeometricFunctions._₂F₁(1/2, -p, 3/2, (1/19)*(2 - 5*x)^2)]
    end

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.1 (a+b x+c x^2)^p.jl


text:
using NBodySimulator
using StochasticDiffEq

const T = 120.0 # °K
const T0 = 90.0 # °K
const kb = 8.3144598e-3 # kJ/(K*mol)
const ϵ = T * kb
const σ = 0.34 # nm
const ρ = 1374/1.6747# Da/nm^3
const m = 39.95# Da
const N = 8
const L = (m*N/ρ)^(1/3)#10.229σ
const R = 0.5*L
const v_dev = sqrt(kb * T / m)
const bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)

const τ = 0.5e-4 # ps or 1e-12 s
const t1 = 0.0
const t2 = 3τ

const parameters = LennardJonesParameters(ϵ, σ, R)
const lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
#const thermostat = AndersenThermostat(90, 0.01/τ)
#thermostat = BerendsenThermostat(90, 2000τ)
const thermostat = LangevinThermostat(70, 10.0)
#thermostat = NoseHooverThermostat(T0, 200τ)
const pbc = CubicPeriodicBoundaryConditions(L)
const simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat, kb);
#result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
#result = @time run_simulation_sde(simulation, ISSEM(symplectic=true,theta=0.5))
result = @time run_simulation(simulation, EM(), dt=τ)

#(rs, grf) = rdf(result)
#(ts, dr2) = msd(result)

t = t1:τ:result.solution.t[end-1]
temper = @time temperature.(result, t)

#using Plots
#pl=plot(t, temper, ylim=[0,200], xlabel="t, ps", ylabel = "T, °K", label="Temperature, °K", linewidth=2)
#plot!(pl, t, T0*ones(length(t)), label = "90 °K", linewidth=2)
#plot!(pl, title="Andersen thermostat at dt*v=$(thermostat.ν*τ) ")
#plot!(pl, title="Berendsen thermostat at tau=$(thermostat.τ) ps")
#plot!(pl, title="Nose-Hoover thermostat at tau=$(thermostat.τ) ps")

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

#@time save_to_pdb(result, "D:/liquid argon simulation $Nactual molecules and $timesteps steps $time_now.pdb" )

#using JLD
#save("d:/nosehoover thermostat for water liquid argon $T0 and $(thermostat.τ) _$time_now.jld", "t",t, "temper", temper, "T0", T0, "τ", thermostat.τ)

#=
using MAT

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

file = matopen("d:/liquid argon $(length(bodies)) omm units at temperature $T0 andgamma $(thermostat.γ)  $timesteps timesteps $T0 _$time_now.mat", "w")
write(file, "t", collect(t))
write(file, "temper", temper)
write(file, "T0", T0)
if thermostat isa LangevinThermostat
    write(file, "gamma", thermostat.γ)
end
#write(file, "rs", rs)
#write(file, "ts", ts)
#write(file, "grf", grf)
#write(file, "dr2", dr2)
close(file)
=#


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DensityProduct{
    N,
    D<:NTuple{N,AbstractDensity},
    B<:AbstractParamBounds
} <: AbstractDensity
    densities::D
    bounds::B

    # ToDo: Check for equal number of parameters in DensityProduct ctor.

    DensityProduct(densities::D, bounds::B) where {N,D<:NTuple{N,AbstractDensity},B<:AbstractParamBounds} =
        new{N,D,B}(densities, bounds)

end

Base.parent(density::DensityProduct) = density.densities

param_bounds(density::DensityProduct) = density.bounds

nparams(density::DensityProduct) = nparams(density.bounds)

import Base.*
function *(a::AbstractDensity, b::AbstractDensity)
    nparams(a) != nparams(b) && throw(ArgumentError("Can't multiply densities with different number of arguments"))
    new_bounds = param_bounds(a) ∩ param_bounds(b)
    _unsafe_prod(a, b, new_bounds)
end

_unsafe_prod(a::AbstractDensity, b::AbstractDensity, new_bounds::AbstractParamBounds) =
    DensityProduct((a,b), new_bounds)

_unsafe_prod(a::AbstractDensity, b::DensityProduct, new_bounds::AbstractParamBounds) =
    DensityProduct((a, b.densities...), new_bounds)

_unsafe_prod(a::DensityProduct, b::DensityProduct, new_bounds::AbstractParamBounds) =
    DensityProduct((a.densities..., b.densities...), new_bounds)

_unsafe_prod(a::DensityProduct, b::AbstractDensity, new_bounds::AbstractParamBounds) =
    DensityProduct((a.densities...,b), new_bounds)

function density_logval(density::DensityProduct, params::AbstractVector{<:Real})
    ds = density.densities
    isempty(ds) && throw(ArgumentError("Can't evaluate density_logval on empty DensityProduct"))
    sum(map(d -> density_logval(d, params), ds))
end


text:
# A little trick for travis
using PkgBenchmark

pkgpath = dirname(dirname(pathof(SymbolicUtils)))
# move it out of the repository so that you can check out different branches
script = tempname() * ".jl"
benchpath = joinpath(pkgpath, "benchmark", "benchmarks.jl")
cp(benchpath, script)

j = judge(pkgpath, "master", retune=true, script=script)

println("MASTER BRANCH")
println(j.baseline_results)

println("THIS BRANCH")
println(j.target_results)

println("DIFFGROUP")
println(j.benchmarkgroup)

println("MARKDOWN")
export_markdown(stdout, j, export_invariants=true)


text:
using ITensors
using Random
using Test

using ChainRulesCore: rrule_via_ad

include("utils/chainrulestestutils.jl")

using Zygote: ZygoteRuleConfig, gradient

Random.seed!(1234)

@testset "ChainRules rrules: basic ITensor operations" begin
  i = Index(2, "i")
  j = Index(2, "j")
  A = randomITensor(i', dag(i))
  V = randomITensor(i)
  Ac = randomITensor(ComplexF64, i', dag(i))
  B = randomITensor(i', dag(i))
  C = ITensor(3.4)
  D = randomITensor(i', j)

  test_rrule(getindex, ITensor(3.4); check_inferred=false)
  test_rrule(getindex, A, 1, 2; check_inferred=false)
  test_rrule(contract, A', A; check_inferred=false)
  test_rrule(*, 3.2, A; check_inferred=false)
  test_rrule(*, A, 4.3; check_inferred=false)
  test_rrule(+, A, B; check_inferred=false)
  test_rrule(prime, A; check_inferred=false)
  test_rrule(prime, A, 2; check_inferred=false)
  test_rrule(prime, A; fkwargs=(; tags="i"), check_inferred=false)
  test_rrule(prime, A; fkwargs=(; tags="x"), check_inferred=false)
  test_rrule(setprime, D, 2; check_inferred=false)
  test_rrule(noprime, D; check_inferred=false)
  test_rrule(replaceprime, A, 1 => 2; check_inferred=false)
  test_rrule(replaceprime, A, 1, 2; check_inferred=false)
  test_rrule(swapprime, A, 0 => 1; check_inferred=false)
  test_rrule(swapprime, A, 0, 1; check_inferred=false)
  test_rrule(addtags, A, "x"; check_inferred=false)
  test_rrule(addtags, A, "x"; fkwargs=(; plev=1), check_inferred=false)
  test_rrule(removetags, A, "i"; check_inferred=false)
  test_rrule(replacetags, A, "i" => "j"; check_inferred=false)
  test_rrule(replacetags, A, "i", "j"; check_inferred=false)
  test_rrule(settags, A, "x"; check_inferred=false)
  test_rrule(settags, A, "x"; fkwargs=(; plev=1), check_inferred=false)
  test_rrule(
    swaptags, randomITensor(Index(2, "i"), Index(2, "j")), "i" => "j"; check_inferred=false
  )
  test_rrule(
    swaptags, randomITensor(Index(2, "i"), Index(2, "j")), "i", "j"; check_inferred=false
  )
  test_rrule(replaceind, A, i' => sim(i); check_inferred=false)
  test_rrule(replaceind, A, i', sim(i); check_inferred=false)
  test_rrule(replaceinds, A, (i, i') => (sim(i), sim(i)); check_inferred=false)
  test_rrule(replaceinds, A, (i, i'), (sim(i), sim(i)); check_inferred=false)
  test_rrule(swapind, A, i', i; check_inferred=false)
  test_rrule(swapinds, A, (i',), (i,); check_inferred=false)
  test_rrule(itensor, randn(2, 2), i', i; check_inferred=false)
  test_rrule(itensor, randn(2, 2), [i', i]; check_inferred=false)
  test_rrule(itensor, randn(4), i', i; check_inferred=false)
  test_rrule(ITensor, randn(2, 2), i', i; check_inferred=false)
  test_rrule(ITensor, randn(2, 2), [i', i]; check_inferred=false)
  test_rrule(ITensor, randn(4), i', i; check_inferred=false)
  test_rrule(ITensor, 2.3; check_inferred=false)
  test_rrule(dag, A; check_inferred=false)
  test_rrule(permute, A, reverse(inds(A)); check_inferred=false)

  test_rrule(ZygoteRuleConfig(), apply, A, V; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (A, B)
    AT = ITensor(A, i, j)
    BT = ITensor(B, j, i)
    return (BT * AT)[1]
  end
  args = (rand(2, 2), rand(2, 2))
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  args = (rand(4), rand(4))
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  args = (rand(4), rand(2, 2))
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  a, b, k, l, m, n, u, v = Index.([10, 50, 5, 2, 10, 12, 7, 50])
  args = (
    randomITensor(a, b, k),
    randomITensor(a, l, m),
    randomITensor(b, u, n),
    randomITensor(u, v),
    randomITensor(k, v),
    randomITensor(l, m, n),
  )
  f = (args...) -> contract([args...])[] # Left associative
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  seq = ITensors.optimal_contraction_sequence([args...])
  f = (args...) -> contract([args...]; sequence=seq)[] # sequence
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = function (x)
    b = itensor([0, 0, 1, 1], i, j)
    k = itensor([0, 1, 0, 0], i, j)
    T = itensor([0 x x^2 1; 0 0 sin(x) 0; 0 cos(x) 0 exp(x); x 0 0 0], i', j', i, j)
    return x * real((b' * T * k)[])
  end
  args = (0.3,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  #f = function (x)
  #  b = itensor([0, 0, 1, 1], i, j)
  #  k = itensor([0, 1, 0, 0], i, j)
  #  T = itensor([0 x x^2 1; 0 0 sin(x) 0; 0 cos(x) 0 exp(x); x 0 0 0], i, j, i', j')
  #  return x * real((b' * T * k)[])
  #end
  #args = (0.3,)
  #test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  f = x -> sin(scalar(x)^3)
  args = (C,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> sin(x[]^3)
  args = (C,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = adjoint
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x * y)[1, 1]
  args = (A', A)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(x, 2)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> x'[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> addtags(x, "x")[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (prime(x) * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ((x'' * x') * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x'' * (x' * x))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y, z) -> (x * y * z)[1, 1]
  args = (A'', A', A)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x''' * x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x''' * x'' * x' * x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x + y)[1, 1]
  args = (A, B)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (2x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + 2x)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x + 2 * mapprime(x' * x, 2 => 1))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = (x, y) -> (x * y)[]
  args = (A, δ(dag(inds(A))))
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x * x)[]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x * δ(dag(inds(x))))[]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    y = x' * x
    tr = δ(dag(inds(y)))
    return (y * tr)[]
  end
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    y = x'' * x' * x
    tr = δ(dag(inds(y)))
    return (y * tr)[]
  end
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x^2 * δ((i', i)))[1, 1]
  args = (6.2,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> (x^2 * δ(i', i))[1, 1]
  args = (5.2,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> itensor([x^2 x; x^3 x^4], i', i)
  args = (2.54,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ITensor([x^2 x; x^3 x^4], i', i)
  args = (2.1,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> ITensor(x)
  args = (2.12,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
    return real((dag(T) * T)[])
  end
  args = (2.8,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  args = (2.8 + 3.1im,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    v = itensor([exp(-3.2x), cos(2x^2)], j)
    T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
    return real((dag(v') * T * v)[])
  end
  args = (2.8,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #args = (2.8 + 3.1im,)
  #test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = function (x)
    return real((x^3 * ITensor([sin(x) exp(-2x); 3x^3 x+x^2], j', dag(j)))[1, 1])
  end
  args = (3.4 + 2.3im,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(permute(x, reverse(inds(x))))[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  f = x -> prime(x; plev=1)[1, 1]
  args = (A,)
  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)

  W = itensor([1 1] / √2, i)
  f = x -> inner(W', exp(x), W)
  args = (A,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1e-3,
    atol=1e-3,
  )

  f = x -> inner(V', exp(x), V)
  args = (A,)
  test_rrule(
    ZygoteRuleConfig(),
    f,
    args...;
    rrule_f=rrule_via_ad,
    check_inferred=false,
    rtol=1e-4,
    atol=1e-4,
  )
end

@testset "ChainRules rrules: op" begin
  s = siteinds("Qubit", 4)

  # RX
  args = (0.2,)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> op("Rx", s, 1; θ=x)[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # RY
  args = (0.2,)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> op("Ry", s, 1; θ=x)[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # RZ
  args = (0.2,)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> op("Rz", s, 1; ϕ=x)[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # Rn
  args = (0.2, 0.3, 0.4)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> op("Rn", s, 1; θ=x[1], ϕ=x[2], λ=x[3])[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args; rrule_f=rrule_via_ad, check_inferred=false)
  end

  basis = vec(collect(Iterators.product(fill([1, 2], 2)...)))
  # CRx
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("CRx", s, (1, 2); θ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end

  # CRy
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("CRy", s, (1, 2); θ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end

  # CRz
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("CRz", s, (1, 2); ϕ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # Rn
  args = (0.2, 0.3, 0.4)
  for σ in basis, σ′ in basis
    f = x -> op("CRn", s, (1, 2); θ=x[1], ϕ=x[2], λ=x[3])[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args; rrule_f=rrule_via_ad, check_inferred=false)
  end

  # Rxx
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("Rxx", s, (1, 2); ϕ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # Ryy
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("Ryy", s, (1, 2); ϕ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # Rzz
  args = (0.2,)
  for σ in basis, σ′ in basis
    f = x -> op("Rzz", s, (1, 2); ϕ=x)[σ..., σ′...]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end

  # algebra with non-parametric gates
  args = (0.2,)
  # addition
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> x * op("H + Y", s[1])[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  #subtraction
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> x * op("H - Y", s[1])[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # product
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> x * op("H * Y", s[1])[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end
  # composite
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> x * op("H + X * Y", s[1])[σ, σ′]
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  end

  ## algebra with parametric gates
  #args = (0.2,)
  ## addition
  #for σ in [1, 2], σ′ in [1, 2]
  #  f = x -> x * op("H + Rx", s[1]; θ = x)[σ, σ′]
  #  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #end
  ##subtraction
  #for σ in [1, 2], σ′ in [1, 2]
  #  f = x -> x * op("H - Rx", s[1]; θ = x)[σ, σ′]
  #  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #end
  ### product
  #for σ in [1, 2], σ′ in [1, 2]
  #  f = x -> x * op("Rx * Y", s[1]; θ = x)[σ, σ′]
  #  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #end
  ## composite
  #for σ in [1, 2], σ′ in [1, 2]
  #  f = x -> x * op("Rx * Y - Ry", s[1]; θ = x)[σ, σ′]
  #  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #end
  #
  ## two-qubit composite algebra with parametric gate
  #args = (0.2,)
  #for σ in basis, σ′ in basis
  #  f = x -> op("Rxx + CX * CZ - Ryy", s, (1, 2); ϕ = x)[σ..., σ′...]
  #  test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  #end

  # functions
  f = x -> exp(ITensor(Op("Ry", 1; θ=x), q))[1, 1]

  # RX
  args = (0.2,)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> exp(ITensor(Op("Rx", 1; θ=x), s))[σ, σ′]
    test_rrule(
      ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false, atol=1e-6
    )
  end

  # RY
  args = (0.2,)
  for σ in [1, 2], σ′ in [1, 2]
    f = x -> exp(ITensor(Op("Ry", 1; θ=x), s))[σ, σ′]
    test_rrule(
      ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false, atol=1e-6
    )
  end
end

@testset "MPS ($ElType)" for ElType in (Float64, ComplexF64)
  Random.seed!(1234)
  n = 4
  ϵ = 1e-8
  s = siteinds("S=1/2", n; conserve_qns=true)
  function heisenberg(n)
    os = OpSum()
    for j in 1:(n - 1)
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
      os += "Sz", j, "Sz", j + 1
    end
    return os
  end
  H = MPO(heisenberg(n), s)
  ψ = randomMPS(s, n -> isodd(n) ? "Up" : "Dn"; linkdims=2)

  f = x -> inner(x, x)
  args = (ψ,)
  d_args = gradient(f, args...)
  @test norm(d_args[1] - 2 * args[1]) ≈ 0 atol = 1e-13

  f = x -> inner(x', H, x)
  args = (ψ,)
  d_args = gradient(f, args...)
  @test norm(d_args[1]' - 2 * H * args[1]) ≈ 0 atol = 1e-13

  f = x -> inner(x', x)
  args = (ψ,)
  @test_throws ErrorException gradient(f, args...)

  f = x -> inner(x, H, x)
  args = (ψ,)
  @test_throws ErrorException gradient(f, args...)

  # apply on MPS
  s = siteinds("S=1/2", n)
  ϕ = randomMPS(ElType, s)
  ψ = randomMPS(ElType, s)
  f = function (x)
    U = [op("Ry", s[2]; θ=x), op("CX", s[1], s[2]), op("Rx", s[3]; θ=x)]
    ψθ = apply(U, ψ)
    return abs2(inner(ϕ, ψθ))
  end
  θ = 0.5
  ∇f = f'(θ)
  ∇num = (f(θ + ϵ) - f(θ)) / ϵ
  @test ∇f ≈ ∇num atol = 1e-5
end

@testset "MPS rrules" begin
  Random.seed!(1234)
  s = siteinds("S=1/2", 4)
  ψ = randomMPS(s)
  args = (ψ,)
  f = x -> inner(x, x)
  # TODO: Need to make MPS type compatible with FiniteDifferences.
  #test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  d_args = gradient(f, args...)
  @test norm(d_args[1] - 2 * args[1]) ≈ 0 atol = 1e-13

  args = (ψ,)
  f = x -> inner(prime(x), prime(x))
  # TODO: Need to make MPS type compatible with FiniteDifferences.
  #test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
  d_args = gradient(f, args...)
  @test norm(d_args[1] - 2 * args[1]) ≈ 0 atol = 1e-13

  ψ = randomMPS(ComplexF64, s)
  ψtensors = ITensors.data(ψ)
  ϕ = randomMPS(ComplexF64, s)
  f = function (x)
    ψ̃tensors = [x^j * ψtensors[j] for j in 1:length(ψtensors)]
    ψ̃ = MPS(ψ̃tensors)
    return abs2(inner(ϕ, ψ̃))
  end
  x = 0.5
  ϵ = 1e-10
  @test f'(x) ≈ (f(x + ϵ) - f(x)) / ϵ atol = 1e-6

  ρ = randomMPO(s)
  f = function (x)
    ψ̃tensors = [x^j * ψtensors[j] for j in 1:length(ψtensors)]
    ψ̃ = MPS(ψ̃tensors)
    return real(inner(ψ̃', ρ, ψ̃))
  end
  @test f'(x) ≈ (f(x + ϵ) - f(x)) / ϵ atol = 1e-6
end

#@testset "MPO rules" begin
#  Random.seed!(1234)
#  s = siteinds("S=1/2", 2)
#
#  #ρ = randomMPO(s)
#  #ρtensors = ITensors.data(ρ)
#  #ϕ = randomMPS(ComplexF64, s)
#  #f = function (x)
#  #  ρ̃tensors  = [2 * x * ρtensors[1],  log(x) * ρtensors[2]]
#  #  ρ̃ = MPO(ρ̃tensors)
#  #  #@show typeof(ρ̃)
#  #  return real(inner(ϕ', ρ̃, ϕ))
#  #end
#  #x = 3.0
#  #ϵ = 1e-8
#  #@show (f(x+ϵ) - f(x)) / ϵ
#  #@show f'(x)
#  ##@test f'(x) ≈ (f(x+ϵ) - f(x)) / ϵ atol = 1e-6
#  #
#
#  #ϕ = randomMPO(s)
#  #f = function (x)
#  #  ψ̃tensors  = [2 * x * ψtensors[1],  log(x) * ψtensors[2]]
#  #  ψ̃ = MPS(ψ̃tensors)
#  #  return abs2(inner(ϕ, ψ̃))
#  #end
#  #x = 3.0
#  #ϵ = 1e-8
#  #@test f'(x) ≈ (f(x+ϵ) - f(x)) / ϵ atol = 1e-6
#
#  #ρ = randomMPO(s)
#end
@testset "MPO: apply" begin
  Random.seed!(1234)
  ϵ = 1e-8
  n = 3
  s = siteinds("Qubit", n)
  function ising(n, h)
    os = OpSum()
    for j in 1:(n - 1)
      os += -1, "Z", j, "Z", j + 1
      os += -h, "X", j
    end
    os += -h, "X", n
    return os
  end
  H = MPO(ising(n, 1.0), s)

  # apply on MPO with apply_dag=true
  ϕ = randomMPS(ComplexF64, s; linkdims=10)
  f = function (x)
    U = [op("Ry", s[2]; θ=x), op("CX", s[1], s[2]), op("Rx", s[3]; θ=x)]
    Hθ = apply(U, H; apply_dag=true)
    return real(inner(ϕ', Hθ, ϕ))
  end
  θ = 0.5
  ∇f = f'(θ)
  ∇num = (f(θ + ϵ) - f(θ)) / ϵ
  @test ∇f ≈ ∇num atol = 1e-5

  # apply on MPO with apply_dag=false
  f = function (x)
    U = [op("Ry", s[2]; θ=x), op("CX", s[1], s[2]), op("Rx", s[3]; θ=x)]
    Hθ = apply(U, H; apply_dag=false)
    return real(inner(ϕ', Hθ, ϕ))
  end
  θ = 0.5
  ∇f = f'(θ)
  ∇num = (f(θ + ϵ) - f(θ)) / ϵ
  @test ∇f ≈ ∇num atol = 1e-5

  # multiply two MPOs
  V = randomMPO(s)
  f = function (x)
    U = [op("Ry", s[2]; θ=x), op("CX", s[1], s[2]), op("Rx", s[3]; θ=x)]
    Hθ = apply(U, H; apply_dag=false)
    X = replaceprime(V' * Hθ, 2 => 1)
    return real(inner(ϕ', X, ϕ))
  end

  θ = 0.5
  ∇f = f'(θ)
  ∇num = (f(θ + ϵ) - f(θ)) / ϵ
  @test ∇f ≈ ∇num atol = 1e-5

  # trace(MPO)
  V1 = randomMPO(s)
  V2 = randomMPO(s)
  f = function (x)
    U = [op("Ry", s[2]; θ=x), op("CX", s[1], s[2]), op("Rx", s[3]; θ=x)]
    Hθ = apply(U, H; apply_dag=false)
    X = V1''' * Hθ'' * V2' * Hθ
    return real(tr(X; plev=4 => 0))
  end

  θ = 0.5
  ∇f = f'(θ)
  ∇num = (f(θ + ϵ) - f(θ)) / ϵ
  @test ∇f ≈ ∇num atol = 1e-5
end


text:
function pollEvent!()
    #SDL2.Event() = [SDL2.Event(NTuple{56, Uint8}(zeros(56,1)))]
    SDL_Event() = Array{UInt8}(zeros(56))
    e = SDL_Event()
    success = (SDL2.PollEvent(e) != 0)
    return e,success
end
function getEventType(e::Array{UInt8})
    # HAHA This is still pretty janky, but I guess that's all you can do w/ unions.
    bitcat(UInt32, e[4:-1:1])
end
function getEventType(e::SDL2.Event)
    bitcat(UInt32, e[4:-1:1])
end

function bitcat(::Type{T}, arr)::T where T<:Number
    out = zero(T)
    for x in arr
        out = out << T(sizeof(x)*8)
        out |= convert(T, x)  # the `convert` prevents signed T from promoting to Int64.
    end
    out
end


text:
# using CUDA
using CUDA
# using Flux
# using GeometricFlux

nbins = 32
ncol = 100
items = Int(1e6)
hist = zeros(Float32, nbins, ncol)
δ = rand(Float32, items)
idx = rand(1:nbins, items, ncol)
𝑖 = collect(1:items)
𝑗 = collect(1:ncol)

hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)
𝑖_gpu = CuArray(𝑖)
𝑗_gpu = CuArray(𝑗)

# CPU
function hist_cpu!(hist, δ, idx, 𝑖, 𝑗)
    Threads.@threads for j in 𝑗
        @inbounds for i in 𝑖
            hist[idx[i], j] += δ[i]
        end
    end
    return
end

function kernel_1!(h::CuDeviceMatrix{T}, x::CuDeviceVector{T}, id, 𝑖, 𝑗) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k), x[𝑖[i]])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu_1!(h::CuMatrix{T}, x::CuVector{T}, id::CuMatrix{Int}, 𝑖, 𝑗; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_1!(h, x, id, 𝑖, 𝑗)
    return
end

@time hist_cpu!(hist, δ, idx)
CUDA.@time hist_gpu_1!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)

nbins = 32
ncol = 100
items = Int(2e6)
K = 1
hist = zeros(Float32, nbins, 3, ncol)
δ = rand(Float32, items, 3)
idx = rand(1:nbins, items, ncol)
𝑖 = collect(1:items)
𝑗 = collect(1:ncol)

hist_gpu = CuArray(hist)
δ_gpu = CuArray(δ)
idx_gpu = CuArray(idx)
𝑖_gpu = CuArray(𝑖)
𝑗_gpu = CuArray(𝑗)

function kernel_2!(h::CuDeviceArray{T,3}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds k1 = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 1, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k1), x[𝑖[i],1])
        @inbounds k2 = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 2, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k2), x[𝑖[i],2])
        @inbounds k3 = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 3, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, k3), x[𝑖[i],3])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu_2!(h::CuArray{T,3}, x::CuMatrix{T}, id::CuMatrix{Int}, 𝑖, 𝑗; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_2!(h, x, id, 𝑖, 𝑗)
    return
end

CUDA.@time hist_gpu_2!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, MAX_THREADS=1024)

hist_gpu_1 = Array(hist_gpu)
hist_gpu_2 = Array(hist_gpu)
diff1 = hist_gpu_2 - hist_gpu_1

######################################################################################################
# best approach: loop on K indicators
######################################################################################################
function kernel_3!(h::CuDeviceArray{T,3}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗, K) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        for k in 1:K
            @inbounds pt = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], k, 𝑗[j])
            @inbounds CUDA.atomic_add!(pointer(h, pt), x[𝑖[i],k])
        end
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu_3!(h::CuArray{T,3}, x::CuMatrix{T}, id::CuMatrix{Int}, 𝑖, 𝑗, K; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks=blocks threads=threads kernel_3!(h, x, id, 𝑖, 𝑗, K)
    return
end

hist_gpu_1 = Array(hist_gpu)
hist_gpu_2 = Array(hist_gpu)
diff2 = hist_gpu_2 - hist_gpu_1
diff2 - diff1

CUDA.@time hist_gpu_3!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, 3, MAX_THREADS=1024)

######################################################################################################
# 3D kernel - instead of iterating on K - Less efficient than the loop on Ks
######################################################################################################
function kernel_3D!(h::CuDeviceArray{T,3}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗, K) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    k = threadIdx().z + (blockIdx().z - 1) * blockDim().z
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds pt = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], k, 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, pt), x[𝑖[i],k])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu_3D!(h::CuArray{T,3}, x::CuMatrix{T}, id::CuMatrix{Int}, 𝑖, 𝑗, K; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_k = min(MAX_THREADS, K)
    thread_j = min(MAX_THREADS ÷ thread_k, length(𝑗))
    thread_i = min(MAX_THREADS ÷ (thread_k * thread_j), length(𝑖))
    threads = (thread_i, thread_j, thread_k)
    blocks = ceil.(Int, (length(𝑖), length(𝑗), K) ./ threads)
    @cuda blocks=blocks threads=threads kernel_3D!(h, x, id, 𝑖, 𝑗, K)
    return
end

CUDA.@time hist_gpu_3D!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, 3, MAX_THREADS=1024)

hist_gpu_1 = Array(hist_gpu)
hist_gpu_2 = Array(hist_gpu)
diff1 = hist_gpu_2 - hist_gpu_1

######################################################################################################
# 3D kernel - instead of iterating on K - No collision approach - single i thread - bad!
######################################################################################################
function kernel_3D2!(h::CuDeviceArray{T,3}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗, K) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    k = threadIdx().z + (blockIdx().z - 1) * blockDim().z
    if i <= length(𝑖) && j <= length(𝑗)
        # @inbounds pt = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], k, 𝑗[j])
        @inbounds h[id[𝑖[i], 𝑗[j]], k, 𝑗[j]] += x[𝑖[i],k]
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function hist_gpu_3D2!(h::CuArray{T,3}, x::CuMatrix{T}, id::CuMatrix{Int}, 𝑖, 𝑗, K; MAX_THREADS=1024) where {T<:AbstractFloat}
    thread_k = min(MAX_THREADS, K)
    thread_j = min(MAX_THREADS ÷ thread_k, length(𝑗))
    thread_i = 1
    threads = (thread_i, thread_j, thread_k)
    blocks = ceil.(Int, (length(𝑖), length(𝑗), K) ./ threads)
    @cuda blocks=blocks threads=threads kernel_3D2!(h, x, id, 𝑖, 𝑗, K)
    return
end

CUDA.@time hist_gpu_3D2!(hist_gpu, δ_gpu, idx_gpu, 𝑖_gpu, 𝑗_gpu, 3, MAX_THREADS=1024)

hist_gpu_1 = Array(hist_gpu)
hist_gpu_2 = Array(hist_gpu)
diff1 = hist_gpu_2 - hist_gpu_1