text:
struct CellPosition
    row::Int
    column::Int
end

"""
A `CellRef` represents a cell location given by row and column identifiers.

`CellRef("A6")` indicates a cell located at column `1` and row `6`.

Example:

```julia
cn = XLSX.CellRef("AB1")
println( XLSX.row_number(cn) ) # will print 1
println( XLSX.column_number(cn) ) # will print 28
println( string(cn) ) # will print out AB1
```

As a convenience, `@ref_str` macro is provided.

```julia
cn = XLSX.ref"AB1"
println( XLSX.row_number(cn) ) # will print 1
println( XLSX.column_number(cn) ) # will print 28
println( string(cn) ) # will print out AB1
```
"""
struct CellRef
    name::String
    row_number::Int
    column_number::Int
end

abstract type AbstractCell end

mutable struct Cell <: AbstractCell
    ref::CellRef
    datatype::String
    style::String
    value::String
    formula::String
end

struct EmptyCell <: AbstractCell
    ref::CellRef
end

abstract type AbstractCellDataFormat end

struct EmptyCellDataFormat <: AbstractCellDataFormat end

"""
Keeps track of formatting information.
"""
struct CellDataFormat <: AbstractCellDataFormat
    id::UInt
end

const CellValueType = Union{String, Missing, Float64, Int, Bool, Dates.Date, Dates.Time, Dates.DateTime}

"""
CellValue is a Julia type of a value read from a Spreadsheet.
"""
struct CellValue
    value::CellValueType
    styleid::AbstractCellDataFormat
end

"""
A `CellRange` represents a rectangular range of cells in a spreadsheet.

`CellRange("A1:C4")` denotes cells ranging from `A1` (upper left corner) to `C4` (bottom right corner).

As a convenience, `@range_str` macro is provided.

```julia
cr = XLSX.range"A1:C4"
```
"""
struct CellRange
    start::CellRef
    stop::CellRef

    function CellRange(a::CellRef, b::CellRef)

        top = row_number(a)
        bottom = row_number(b)
        left = column_number(a)
        right = column_number(b)

        @assert left <= right && top <= bottom "Invalid CellRange. Start cell should be at the top left corner of the range."

        return new(a, b)
    end
end

struct ColumnRange
    start::Int # column number
    stop::Int  # column number

    function ColumnRange(a::Int, b::Int)
        @assert a <= b "Invalid ColumnRange. Start column must be located before end column."
        return new(a, b)
    end
end

struct SheetCellRef
    sheet::String
    cellref::CellRef
end

struct SheetCellRange
   sheet::String
   rng::CellRange
end

struct SheetColumnRange
    sheet::String
    colrng::ColumnRange
end

abstract type MSOfficePackage end

struct EmptyMSOfficePackage <: MSOfficePackage
end

"""
Relationships are defined in ECMA-376-1 Section 9.2.
This struct matches the `Relationship` tag attribute names.

A `Relashipship` defines relations between the files inside a MSOffice package.
Regarding Spreadsheets, there are two kinds of relationships:

    * package level: defined in `_rels/.rels`.
    * workbook level: defined in `xl/_rels/workbook.xml.rels`.

The function `parse_relationships!(xf::XLSXFile)` is used to parse
package and workbook level relationships.
"""
struct Relationship
    Id::String
    Type::String
    Target::String
end

const CellCache = Dict{Int, Dict{Int, Cell}} # row -> ( column -> cell )

"""
Iterates over Worksheet cells. See `eachrow` method docs.
Each element is a `SheetRow`.

Implementations: SheetRowStreamIterator, WorksheetCache.
"""
abstract type SheetRowIterator end

mutable struct SheetRowStreamIteratorState
    zip_io::ZipFile.Reader
    xml_stream_reader::EzXML.StreamReader
    is_open::Bool # indicated if zip_io and xml_stream_reader are opened
    row::Int # number of current row. It´s set to 0 in the start state.
end

mutable struct WorksheetCache{I<:SheetRowIterator} <: SheetRowIterator
    cells::CellCache # SheetRowNumber -> Dict{column_number, Cell}
    rows_in_cache::Vector{Int} # ordered vector with row numbers that are stored in cache
    row_index::Dict{Int, Int} # maps a row number to the index of the row number in rows_in_cache
    stream_iterator::I
    stream_state::Union{Nothing, SheetRowStreamIteratorState}
end

mutable struct Worksheet
    package::MSOfficePackage # parent XLSXFile
    sheetId::Int
    relationship_id::String # r:id="rId1"
    name::String
    dimension::CellRange
<<<<<<< HEAD
    cache::Union{Missing, WorksheetCache}

    function Worksheet(package::MSOfficePackage, sheetId::Int, relationship_id::String, name::String, dimension::CellRange)
=======
    cache::Union{WorksheetCache, Nothing}

    function Worksheet(package::MSOfficePackage, sheetId::Int, relationship_id::String, name::String, dimension::CellRange)
        return new(package, sheetId, relationship_id, name, dimension, nothing)
>>>>>>> upstream/master
    end
end

struct SheetRowStreamIterator <: SheetRowIterator
    sheet::Worksheet
end

"""
Shared String Table
"""
mutable struct SharedStrings
    unformatted_strings::Vector{String}
    formatted_strings::Vector{String}
    hashmap::Dict{UInt64, Int} # shared string hash -> index in strings vector.
    is_loaded::Bool # for lazy-loading of sst XML file
end

const DefinedNameValueTypes = Union{SheetCellRef, SheetCellRange, Int, Float64, String, Missing}

"""
Workbook is the result of parsing file `xl/workbook.xml`.
"""
mutable struct Workbook
    package::MSOfficePackage # parent XLSXFile
    sheets::Vector{Worksheet} # workbook -> sheets -> <sheet name="Sheet1" r:id="rId1" sheetId="1"/>. sheetId determines the index of the WorkSheet in this vector.
    date1904::Bool              # workbook -> workbookPr -> attribute date1904 = "1" or absent
    relationships::Vector{Relationship} # contains workbook level relationships
    sst::SharedStrings # shared string table
    buffer_styles_is_float::Dict{Int, Bool}      # cell style -> true if is float
    buffer_styles_is_datetime::Dict{Int, Bool}   # cell style -> true if is datetime
    workbook_names::Dict{String, DefinedNameValueTypes} # definedName
    worksheet_names::Dict{Tuple{Int, String}, DefinedNameValueTypes} # definedName. (sheetId, name) -> value.
<<<<<<< HEAD
    styles_xroot::Union{Missing, EzXML.Node}
=======
    styles_xroot::Union{EzXML.Node, Nothing}
>>>>>>> upstream/master
end

"""
`XLSXFile` stores all XML data from an Excel file.

`filepath` is the filepath of the source file for this XLSXFile.
`data` stored the raw XML data. It maps internal XLSX filenames to XMLDocuments.
`workbook` is the result of parsing `xl/workbook.xml`.
"""
mutable struct XLSXFile <: MSOfficePackage
    filepath::AbstractString
    use_cache_for_sheet_data::Bool # indicates wether Worksheet.cache will be fed while reading worksheet cells.
    io::ZipFile.Reader
    io_is_open::Bool
    files::Dict{String, Bool} # maps filename => isread bool
    data::Dict{String, EzXML.Document} # maps filename => XMLDocument
    binary_data::Dict{String, Vector{UInt8}} # maps filename => file content in bytes
    workbook::Workbook
    relationships::Vector{Relationship} # contains package level relationships
    is_writable::Bool # indicates wether this XLSX file can be edited

    function XLSXFile(filepath::AbstractString, use_cache::Bool, is_writable::Bool)
        check_for_xlsx_file_format(filepath)
        io = ZipFile.Reader(filepath)
        xl = new(filepath, use_cache, io, true, Dict{String, Bool}(), Dict{String, EzXML.Document}(), Dict{String, Vector{UInt8}}(), EmptyWorkbook(), Vector{Relationship}(), is_writable)
        xl.workbook.package = xl
        finalizer(close, xl)
        return xl
    end
end

#
# Iterators
#

struct SheetRow
    sheet::Worksheet
    row::Int
    rowcells::Dict{Int, Cell} # column -> value
end

struct Index # based on DataFrames.jl
    lookup::Dict{Symbol, Int} # column label -> table column index
    column_labels::Vector{Symbol}
    column_map::Dict{Int, Int} # table column index (1-based) -> sheet column index (cellref based)

    function Index(column_range::Union{ColumnRange, AbstractString}, column_labels::Vector{Symbol})
        column_range = convert(ColumnRange, column_range)
        @assert length(unique(column_labels)) == length(column_labels) "Column labels must be unique."

        lookup = Dict{Symbol, Int}()
        for (i, n) in enumerate(column_labels)
            lookup[n] = i
        end

        column_map = Dict{Int, Int}()
        for (i, n) in enumerate(column_range)
            column_map[i] = decode_column_number(n)
        end
        return new(lookup, column_labels, column_map)
    end
end

struct TableRowIterator{I<:SheetRowIterator}
    itr::I
    index::Index
    first_data_row::Int
    stop_in_empty_row::Bool
<<<<<<< HEAD
    stop_in_row_function::Union{Function, Nothing}
=======
    stop_in_row_function::Union{Nothing, Function}
>>>>>>> upstream/master
end

struct TableRow
    row::Int # Index of the row in the table. This is not relative to the worksheet cell row.
    index::Index
    cell_values::Vector{CellValueType}
end

struct TableRowIteratorState{S}
    table_row_index::Int
    sheet_row_index::Int
    sheet_row_iterator_state::S
end


text:
@defcomp biodiversity begin
    # Number of species
    nospecies = Variable(index=[time])

    # additive parameter
    bioloss = Parameter(default = 0.003)

    # multiplicative parameter
    biosens = Parameter(default = 0.001)

    # Temperature
    temp = Parameter(index=[time])

    # benchmark temperature change
    dbsta = Parameter()

    # Number of species in the year 2000
    nospecbase = Parameter()

    function run_timestep(p, v, d, t)

        if gettime(t) > 2000
            dt = abs(p.temp[t] - p.temp[t - 1])

            v.nospecies[t] = max(
            p.nospecbase / 100,
            v.nospecies[t - 1] * (1.0 - p.bioloss - p.biosens * dt * dt / p.dbsta / p.dbsta)
            )
        else
            v.nospecies[t] = p.nospecbase
        end
    end
end

ext: jl


text:
module subspace

import Base.==
import ..operators

using ..bases, ..states, ..operators

export SubspaceBasis, projector

"""
A basis describing a subspace in a higher dimensional Hilbert space.

Arguments
---------
superbasis
    Basis of the higher dimensional Hilbert space.
basisstates
    States in respect to the superbasis that span the subspace.
"""
type SubspaceBasis <: Basis
    shape::Vector{Int}
    superbasis::Basis
    basisstates::Vector{Ket}
    basisstates_hash::UInt

    function SubspaceBasis(superbasis::Basis, basisstates::Vector{Ket})
        for state = basisstates
            if state.basis != superbasis
                raise(ArgumentError("The basis of the basisstates has to be the superbasis."))
            end
        end
        basisstates_hash = hash([hash(x.data) for x=basisstates])
        new(Int[length(basisstates)], superbasis, basisstates, basisstates_hash)
    end
end

SubspaceBasis(basisstates::Vector{Ket}) = SubspaceBasis(basisstates[1].basis, basisstates)

==(b1::SubspaceBasis, b2::SubspaceBasis) = b1.superbasis==b2.superbasis && b1.basisstates_hash==b2.basisstates_hash

proj(u::Ket, v::Ket) = dagger(v)*u/(dagger(u)*u)*u

"""
Orthonormalize the given SubspaceBasis using the modified Gram-Schmidt process.
"""
function orthonormalize(b::SubspaceBasis)
    V = b.basisstates
    U = Ket[]
    for (k, v)=enumerate(V)
        u = deepcopy(v)
        for i=1:k-1
            u -= proj(U[i], u)
        end
        normalize!(u)
        push!(U, u)
    end
    return SubspaceBasis(U)
end

"""
Operator projecting states from one subspace to another.
"""
function operators.projector(b1::SubspaceBasis, b2::SubspaceBasis)
    if b1.superbasis != b2.superbasis
        throw(ArgumentError("Both subspace bases have to have the same superbasis."))
    end
    T1 = operators.projector(b1, b1.superbasis)
    T2 = operators.projector(b2.superbasis, b2)
    return T1*T2
end

"""
Operator projecting states from the superspace to the subspace.
"""
function operators.projector(b1::SubspaceBasis, b2::Basis)
    if b1.superbasis != b2
        throw(ArgumentError("Second basis has to be the superbasis of the first one."))
    end
    data = zeros(Complex128, length(b1), length(b2))
    for (i, state) = enumerate(b1.basisstates)
        data[i,:] = state.data
    end
    return DenseOperator(b1, b2, data)
end

"""
Operator projecting states from the subspace to the superspace.
"""
function operators.projector(b1::Basis, b2::SubspaceBasis)
    if b1 != b2.superbasis
        throw(ArgumentError("First basis has to be the superbasis of the second one."))
    end
    data = zeros(Complex128, length(b1), length(b2))
    for (i, state) = enumerate(b2.basisstates)
        data[:,i] = state.data
    end
    return DenseOperator(b1, b2, data)
end

end # module


text:
#!/bin/env julia

println("Hello, World!")

max_stars_repo_name: mstrangfeld/programming-languages

max_issues_repo_name: mstrangfeld/programming-languages

max_forks_repo_name: mstrangfeld/programming-languages


text:
using Optim, LinearAlgebra

# Augmented Lagrange Method code for inequality contrataint
function Augmented_Lagrange_Mult_Inequality_Box(f, g, lower, upper, x, k_max; α=1000., γ=1.01, tol=1e-1000)
    λ, nite = zeros(length(g(x))), 0
    for k in 1 : k_max; nite = k
        ψ(x) = max.(g(x),-λ./(2α))
        F(x) = f(x) + λ ⋅ g(x) + α * sum(ψ(x).^2)
        #x = Optim.minimizer(optimize(F, lower, upper, x, Fminbox(NelderMead()))) # Ordem zero
        x = Optim.minimizer(optimize(F, lower, upper, x, Fminbox(BFGS());autodiff = :forward)) # Ordem 1
        λ .+= 2*α*ψ(x)
        if norm(ψ(x),Inf)<=tol; break; end # Ver Arora pg 481
        α *= γ;
    end
    return x, nite, α, λ
end;

#modulo de elasticidade
E = 3.0E7;
#altura da trelica
L = 10.0;
#angulo do carregamento e carregamento
#cenário 1
#theta = -pi/4;P = 40000.0;
#cenário 2
#theta = -pi/2;P = 30000.0;
#cenário 3
theta = -3*pi/4;P = 20000.0;

#limite inferior das variáveis de projeto
lb=[0.1, 0.1, 0.1, -10.0, -10.0, -10.0];
#limite superior das variáveis de projeto
ub=[6.1, 6, 6, 10.0, 10.0, 10.0];

#chute inicial
#x0 = [A1,A2,A3,x1,x2,x3]
x0 = [6.01,6.0,6.0,-5.0,0.0,5.0];

#Area das barras
A1(x) = x[1];
A2(x) = x[2];
A3(x) = x[3];

#Comprimento das barras
L1(x) = sqrt(L ^ 2 + x[4] ^ 2);
L2(x) = sqrt(L ^ 2 + x[5] ^ 2);
L3(x) = sqrt(L ^ 2 + x[6] ^ 2);

#Matriz de rigidez
k11(x) = E * (A1(x) * x[4] ^ 2 / L1(x) ^ 3 + A2(x) * x[5] ^ 2 / L2(x) ^ 3 + A3(x) * x[6] ^ 2 / L3(x) ^ 3);
k12(x) = E * (A1(x) * L * x[4] / L1(x) ^ 3 + A2(x) * L * x[5] / L2(x) ^ 3 + A3(x) * L * x[6] / L3(x) ^ 3);
k21(x) = k12(x);
k22(x) = E * (A1(x) * L ^ 2 / L1(x) ^ 3 + A2(x) * L ^ 2 / L2(x) ^ 3 + A3(x) * L ^ 2 / L3(x) ^ 3);

K(x) = [k11(x) k12(x); k21(x) k22(x)];

#Vetor de forcas
f1 = P*cos(theta);
f2 = P*sin(theta);
F = [f1,f2];

#Calculo dos deslocamentos
u(x) = K(x)\F;

#Calculo das tensoes
S1(x) = -E / L1(x) ^ 2 * (u(x)[2] * L + u(x)[1] * x[4]);
S2(x) = -E / L2(x) ^ 2 * (u(x)[2] * L + u(x)[1] * x[5]);
S3(x) = -E / L3(x) ^ 2 * (u(x)[2] * L + u(x)[1] * x[6]);

#Restricoes de desigualdades
g1(x) = abs(S1(x)) / 0.5000e4 - 0.1e1;
g2(x) = abs(S2(x)) / 0.20000e5 - 0.1e1;
g3(x) = abs(S3(x)) / 0.50000e4 - 0.1e1;
g(x) = [g1(x), g2(x), g3(x)];

#Funcao Objetivo - Volume das barras
V(x) = A1(x)*L1(x) + A2(x)*L2(x) + A3(x)*L3(x);

#Solve
@time xₑ, k, α, λ = Augmented_Lagrange_Mult_Inequality_Box(V, g, lb,ub, x0, 100)
println("x* : ", xₑ, ", iterations: ", k, ", α : ", α, ", λ : ", λ)


text:
using POMDPs
using FiniteHorizonPOMDPs

struct FHExampleState
    position::Int64
    epoch::Int64
    done::Bool # are we in a terminal state?
end

# initial state constructor
FHExampleState(position::Int64, epoch::Int64)::FHExampleState = FHExampleState(position, epoch, false)
FHExampleState(position::Int64)::FHExampleState = FHExampleState(position, 0, false)
FHExampleState(position::Int64, done::Bool)::FHExampleState = FHExampleState(position, 0, done)

struct FHExample <: MDP{FHExampleState, Symbol} # Note that our MDP is parametarized by the state and the action
    no_states::Int64 # number od statesL
    horizon::Int64
    actions::Vector{Symbol}
    actionCost::Float64
    actionsImpact::Base.ImmutableDict{Symbol, Int64}
    reward_states::Vector{Int64}
    reward::Float64
    discount_factor::Float64 # discount factor
    noise::Float64
end

function POMDPs.isterminal(mdp::FHExample, state::FHExampleState)::Bool
    return state.done
end

function isreward(mdp::FHExample, position::Int64)::Bool
    return position in mdp.reward_states
end

function POMDPs.reward(mdp::FHExample, s::FHExampleState, a::Symbol, sp::FHExampleState)::Float64
    isreward(mdp, sp.position) ? mdp.reward : mdp.actionCost
end

# returns mdp.no_states * 2 long Array with first mdp.no_states corresponding to current epoch's states whose util value is to be updated
# and another mdp.no_states vector of epoch + 1 states as way to store and uset their util value for computations
# first half of vector is true or false depending on terminality of each state, second half is only true because solvers are not iterate through them
function FiniteHorizonPOMDPs.stage_states(mdp::FHExample, epoch::Int64)::Array{FHExampleState}
    mdp_states = FHExampleState[]
    for i=1:mdp.no_states
        push!(mdp_states, FHExampleState(i, epoch, isreward(mdp, i) || epoch==mdp.horizon))
    end

    return mdp_states
end

# Implementation of function that are replacing POMDPs.states, POMDPs.actions and POMDPs.state_index in FiniteHorizonPOMDPs Solver
function FiniteHorizonPOMDPs.stage_stateindex(mdp::FHExample, s::FHExampleState, epoch::Int64)::Int64
    return s.position + (epoch - 1) * mdp.no_states
end

function FiniteHorizonPOMDPs.stage_stateindex(mdp::FHExample, s::FHExampleState)::Int64
    return s.position
end

function POMDPs.actions(mdp::FHExample, s::FHExampleState)
    mdp.actions
end

function POMDPs.actions(mdp::FHExample)
    mdp.actions
end

# Creates (mdp.horizon - 1) * mdp.no_states states to be evaluated and mdp.no_states sink states
function POMDPs.states(mdp::FHExample)::Array{FHExampleState}
    mdp_states = FHExampleState[]
    for e=1:mdp.horizon
        for i=1:mdp.no_states
            push!(mdp_states, FHExampleState(i, e, isreward(mdp, i) || e == mdp.horizon))
        end
    end

    return mdp_states
end

POMDPs.stateindex(mdp::FHExample, s::FHExampleState)::Int64 = s.position + (s.epoch - 1) * mdp.no_states

POMDPs.actions(mdp::FHExample)::Vector{Symbol} = mdp.actions

POMDPs.actionindex(mdp::FHExample, a::Symbol)::Int64 = findall(x->x==a, POMDPs.actions(mdp))[1]

POMDPs.discount(mdp::FHExample)::Number = mdp.discount_factor

# returns transition distributions - works only for 1D Gridworld with possible moves to left and to right
function POMDPs.transition(mdp::FHExample, s::FHExampleState, a::Symbol)::SparseCat{Vector{FHExampleState},Vector{Float64}}
    sp = FHExampleState[]
    prob = Float64[]

    # add original transition target and probability
    position = s.position + mdp.actionsImpact[a]
    push!(sp, FHExampleState(position, s.epoch + 1, isreward(mdp, position)))
    push!(prob, 1. - mdp.noise)

    # add noise transition target and probability
    noise_action = a == :l ? :r : :l
    position = s.position + mdp.actionsImpact[noise_action]
    push!(sp, FHExampleState(position, s.epoch + 1, isreward(mdp, position)))
    push!(prob, mdp.noise)

    return SparseCat(sp, prob)
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m5_opt_rho1_9_low_relnorm(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m5_opt_rho1_9_low_relnorm!(A_copy)
end

@inline function exp_sastre_m5_opt_rho1_9_low_relnorm!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=-0.31551655352795516
    coeff2=0.4801289831900131
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.31551655352795516
    coeff2=0.4801289831900131
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.357355571021064
    coeff2=0.4105933675776958
    coeff3=0.09852026651116912
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.053558315750694814
    coeff2=-0.22014016998509814
    coeff3=0.922091412659715
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.004620908141128493
    coeff2=0.055057556018866645
    coeff3=-0.351290300034291
    coeff4=0.9607057168526784
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.011326704451019963
    coeff2=0.002770498744594938
    coeff3=0.00045550004840354903
    coeff4=0.00012416495659100344
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.11176905543585865
    coeff2=0.5925762886270637
    coeff3=0.4621864853147797
    coeff4=0.038392790243209074
    coeff5=0.24031719666304735
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.1215251552594178
    coeff2=0.2145595643740629
    coeff3=0.04592597180399795
    coeff4=0.021375614386974554
    coeff5=0.512024337751819
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=1.11114318963272
    coeff2=0.4876551310270585
    coeff3=0.3710914185617396
    coeff4=0.08701143872628835
    coeff5=5.180154142825969
    coeff6=0.6495647957423022
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.9070495211990466
    coeff2=0.5443460967617421
    coeff3=0.5323384854781406
    coeff4=0.12934538808055288
    coeff5=5.16322755510199
    coeff6=0.6107117158792567
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.0818796601588853
    coeff2=0.003683007236493971
    coeff3=0.008733812199626457
    coeff4=0.005536448121469381
    coeff5=0.003015592881979151
    coeff6=0.01836716605825683
    coeff7=0.9780133820655359
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
#  breadth first search
function bfs(goal, start)
    reset_timer!(to::TimerOutput)
    node = createtree(start)
    if node.state == goal
        return node
    end
    frontier = [node]
    reached = [node.state]
    while !isempty(frontier)
        @timeit to "Moves" acts = moves(node.state)
        @timeit to "Pop" node = popfirst!(frontier)
        @timeit to "Expand" states = expansion(node.state)

        for (currentstate, currentaction) in zip(states, acts)
            if currentstate == goal
                sovle = addnode(currentaction, node, currentstate)
                displaysolution(sovle)
                return sovle
            end
            if !in(currentstate, reached)

                @timeit to "Node Creation" currentnode = addnode(currentaction, node, currentstate)
                @timeit to "Push" push!(frontier, currentnode)
            end
        end
    end

    println("No solutions found!"^10)
    return nothing
end

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
using DataFrames
using CSV
using FileIO
using Statistics

include("experiments/experimentalutils.jl")

resultsFolder = mainfolder * "experiments/kernelsizetest/"
files = readdir(resultsFolder)

results = []
for f in files
    if isfile(resultsFolder * f)
        push!(results, CSV.read(resultsFolder * f))
    end
end
results = vcat(results...)

aggres = []
for (g, d) in Base.product(unique(results[:γ]), unique(results[:dataset]))
    ddf = results[(results[:dataset] .== d) .& (results[:γ] .== g), :]
    mean_auc_pxv = maximum(ddf[:auc_pxv])
    mean_auc_pxis = maximum(ddf[:auc_pxis])
    mean_auc_pz = maximum(ddf[:auc_pz])
    push!(aggres, DataFrame(dataset = d, γ = g, auc_pxv = mean_auc_pxv, auc_pxis = mean_auc_pxis, auc_pz = mean_auc_pz))
end
aggres = vcat(aggres...)


text:
include("utils.jl")

@testset "Elliptope" begin
    M = Elliptope(4, 2)
    @test repr(M) == "Elliptope(4, 2)"
    @test manifold_dimension(M) == 3
    @test get_embedding(M) == Euclidean(4, 2)
    @test representation_size(M) == (4, 2)
    q = [1.0 0.0; 0.0 1.0; 1/sqrt(2) -1/sqrt(2); 1/sqrt(2) 1/sqrt(2)]
    @test is_point(M, q, true; atol=10^-15)
    @test base_manifold(M) === M
    qN = [2.0 0.0; 0.0 1.0; 1/sqrt(2) -1/sqrt(2); 1/sqrt(2) 1/sqrt(2)]
    @test_throws DomainError is_point(M, qN, true)
    Y = [0.0 1.0; 1.0 0.0; 0.0 0.0; 0.0 0.0]
    @test is_vector(M, q, Y, true)
    YN = [0.1 1.0; 1.0 0.1; 0.0 0.0; 0.0 0.0]
    @test_throws DomainError is_vector(M, q, YN, true)
    qE = similar(q)
    embed!(M, qE, q)
    qE2 = embed(M, q)
    @test qE == q
    @test qE2 == q
    q2 = [4.0/5 3.0/5; 3.0/5.0 -4.0/5.0; 1.0 0.0; 0.0 1.0]
    q3 = [12.0/13.0 5.0/13.0; 1.0 0.0; -12.0/13.0 5.0/13.0; 0.0 1.0]
    @test is_vector(
        M,
        q2,
        vector_transport_to(M, q, Y, q2, ProjectionTransport());
        atol=10^-15,
    )

    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{4,2,Float64,8})

    for T in types
        pts = [convert(T, q), convert(T, q2), convert(T, q3)]
        @testset "Type $T" begin
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_reverse_diff=false,
                test_forward_diff=false,
                test_project_tangent=true,
                test_project_point=true,
                test_exp_log=false,
                default_inverse_retraction_method=nothing,
                default_retraction_method=ProjectionRetraction(),
                is_tangent_atol_multiplier=1,
                test_inplace=true,
            )
        end
    end
end


text:
using POMDPs, POMDPGym, POMDPSimulators, POMDPPolicies, Distributions
using Crux, Flux, BSON, StaticArrays, Random
using FileIO
using Plots

using ParallelTreeSampling

include("utils.jl")

unicodeplots()

# Basic MDP.
tprob = 0.7
Random.seed!(0)
randcosts = Dict(POMDPGym.GWPos(i,j) => rand() for i = 1:10, j=1:10)
# zerocosts = Dict(POMDPGym.GWPos(i,j) => 0.0 for i = 1:10, j=1:10)
mdp = GridWorldMDP(costs=randcosts, cost_penalty=0.1, tprob=tprob)

# Learn a policy that solves it.
# policy = DiscreteNetwork(Chain(x -> (x .- 5f0) ./ 5f0, Dense(2, 32, relu), Dense(32, 4)), [:up, :down, :left, :right])
# policy = solve(DQN(π=policy, S=state_space(mdp), N=100000, ΔN=4, buffer_size=10000, log=(;period=5000)), mdp)
# atable = Dict(s => action(policy, [s...]) for s in states(mdp.g))
# BSON.@save "demo/gridworld_policy_table.bson" atable

atable = BSON.load("demo/policies/gridworld_policy_table.bson")[:atable]

# Define the adversarial mdp.
adv_rewards = deepcopy(randcosts)
# adv_rewards = deepcopy(zerocosts)
for (k,v) in mdp.g.rewards
    if v < 0
        adv_rewards[k] += -10*v
    end
end

amdp = GridWorldMDP(rewards=adv_rewards, tprob=1., discount=1., terminate_from=mdp.g.terminate_from)

# Define action probability for the adv_mdp.
action_probability(mdp, s, a) = (a == atable[s][1]) ? tprob : ((1. - tprob) / (length(actions(mdp)) - 1.))

function distribution(mdp::Union{POMDP,MDP}, s)
    xs = POMDPs.actions(mdp, s)
    ps = [action_probability(mdp, s, x) for x in xs]
    ps ./= sum(ps)
    px = GenericDiscreteNonParametric(xs, ps)
    return px
end

fixed_s = rand(initialstate(amdp))

N = 100_000
base_N = 10_000_000
c = 0.0
α = 0.001; β = 0.1; γ = 0.3;
vloss = 0.0
is_baseline = false

path = "data"

alpha_list = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

function baseline(trial=0)
    baseline_out = run_baseline(amdp, fixed_s, distribution; N=base_N)
    filename = string("gridworld_baseline_$(base_N)",
                      (trial > 0 ? string("_", trial) : "") ,".jld2")
    save(joinpath(path, filename),
         Dict("risks:" => baseline_out[1], "states:" => baseline_out[2]))

    d = Dict([(alpha, Dict()) for alpha in alpha_list])
    for alpha in alpha_list
        m = eval_metrics(baseline_out[1]; alpha)
        d[alpha][:N] = length(baseline_out[1]); d[alpha][:mean] = round(m.mean, digits=3);
        d[alpha][:var] = round(m.var, digits=3); d[alpha][:cvar] = round(m.cvar, digits=3);
        d[alpha][:worst] = round(m.worst, digits=3);
    end
    return d
end

function mcts(trial=0)
    mcts_out, planner = run_mcts(
        amdp, fixed_s, distribution; N=N, c=c, vloss=vloss, α=α, β=β, γ=γ)
    mcts_info = mcts_out[4]
    filename = string("gridworld_mcts_$(N)",
                      (trial > 0 ? string("_", trial) : "") ,".jld2")
    save(joinpath(path, filename),
         Dict("risks:" => mcts_out[1], "states:" => mcts_out[2],
              "weights:" => mcts_out[3], "tree:" => mcts_info[:tree]))
    search_t = round(mcts_info[:search_time_s]; digits=3)

    d = Dict([(alpha, Dict()) for alpha in alpha_list])
    for alpha in alpha_list
        m = eval_metrics(mcts_out[1]; weights=exp.(mcts_out[3]), alpha=alpha)
        d[alpha][:N] = length(mcts_out[1]); d[alpha][:mean] = round(m.mean, digits=3);
        d[alpha][:var] = round(m.var, digits=3); d[alpha][:cvar] = round(m.cvar, digits=3);
        d[alpha][:worst] = round(m.worst, digits=3);
    end
    return d, search_t
end

metrics = []
times = []
num_trials = 1
for trial in 1:num_trials
    Random.seed!(trial)
    if is_baseline
        push!(metrics, baseline(trial))
    else
        metric_d, search_t = mcts(trial)
        push!(metrics, metric_d)
        push!(times, search_t)
    end
end

if is_baseline
    println("Baseline metrics")
else
    println("TIS metrics: N=$(N), c=$(c), α=$(α), β=$(β), γ=$(γ)")
end

for (idx, metric) in enumerate(metrics)
    println("Run $(idx)")
    for alpha in alpha_list
        m_d = metric[alpha]
        println(string("  [Alpha=$(alpha)] N: $(m_d[:N]), Mean: $(m_d[:mean]), ",
                       "VaR: $(m_d[:var]), CVaR: $(m_d[:cvar]), Worst: $(m_d[:worst])"))
    end
end
println("Times: $(times)")


text:
"""
    graph = shortIntGraph(a::SparseMatrixCSC)

Convert the indices in a graph to 32-bit ints.
This takes less storage, but does not speed up much.
"""
shortIntGraph(a::SparseMatrixCSC) = SparseMatrixCSC{Float64,Int32}(convert(Int32,a.m), convert(Int32,a.n), convert(Array{Int32,1},a.colptr), convert(Array{Int32,1},a.rowval), a.nzval)

"""
  graph = floatGraph(a::SparseMatrixCSC)

Convert the nonzero entries in a graph to Float64.
"""
floatGraph(a::SparseMatrixCSC) = SparseMatrixCSC{Float64,Int64}(a.m, a.n, a.colptr, a.rowval, convert(Array{Float64,1},a.nzval))

"""
    l = lap(a)

Create a Laplacian matrix from an adjacency matrix. We might want to do this differently, say by enforcing symmetry
"""
lap(a) = sparse(Diagonal(a*ones(size(a)[1]))) - a

"""
    a,d = adj(sddm)

Create an adjacency matrix and a diagonal vector from an SDD M-matrix.
That is, from a Laplacian with added diagonal weights
"""
function adj(sddm::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    a = sparse(Diagonal(diag(sddm))) - sddm
    d = sddm*ones(size(sddm,1))
    return a,d
end

function adj(sddm::Array{Tv,2}) where Tv
    return adj(sparse(sddm))
end

"""
Add a new vertex to a with weights to the other vertices corresponding to diagonal surplus weight.

This is an efficient way of writing [a d; d' 0]
"""
function extendMatrix(a::SparseMatrixCSC{Tv,Ti}, d::Array{Tv,1}) where {Tv,Ti}

    @assert size(a,1) == length(d)

    if sum(abs.(d)) == 0
        return a
    end

    dpos = d.*(d.>0)

    n = length(d)

    ai,aj,av=findnz(a)

    ai2 = [ai;1:n;(n+1)*ones(Int,n)]
    aj2 = [aj;(n+1)*ones(Int,n);1:n]
    av2 = [av;dpos;dpos]
    a2 = sparse(ai2,aj2,av2)

    return a2

end

"""
    wt1 = unweight(a)

Create a new graph in that is the same as the original, but with all edge weights 1
"""
function unweight(ain::SparseMatrixCSC{Tval,Tind}) where {Tval,Tind}
    a = copy(ain)
    m = length(a.nzval)
    for i in 1:m
        a.nzval[i] = one(Tval)
    end
    return a
end # unweight

"""
    unweight!(a)

Change the weight of every edge in a to 1
"""
function unweight!(a::SparseMatrixCSC{Tval,Tind}) where {Tval,Tind}
    m = length(a.nzval)
    for i in 1:m
        a.nzval[i] = one(Tval)
    end
end # unweight

function unweight!(ijv::IJV)
    ijv.v .= 1.0
end

"""
    b = mapweight(a, x->rand())

Create a new graph that is the same as the original, but with f applied to each nonzero entry of a. For example, to make the weight of every edge uniform in [0,1], we could write.
"""
function mapweight(a::SparseMatrixCSC{Tval,Tind},f) where {Tval,Tind}
    b = triu(a,1)
    for i in 1:length(b.nzval)
        b.nzval[i] = f(b.nzval[i])
    end
    b = b + b'

    return b

end # mapweight

"""
    wted = uniformWeight(unwted)

Put a uniform [0,1] weight on every edge.  This is an example of how to use mapweight."""
uniformWeight_ver(::Type{V06}, a::SparseMatrixCSC)  = mapweight(a,x->rand_ver(V06))
uniformWeight_ver(::Type{Vcur}, a::SparseMatrixCSC)  = mapweight(a,x->rand_ver(Vcur))
uniformWeight(a::SparseMatrixCSC)  = uniformWeight_ver(Vcur, a::SparseMatrixCSC)

"""
    uniformWeight!(a)

Set the weight of every edge to random uniform [0,1]
"""
uniformWeight!(mat::SparseMatrixCSC) = uniformWeight_ver!(Vcur, mat)

function uniformWeight_ver!(ver, mat::SparseMatrixCSC)
    for i in 1:length(mat.nzval)
        mat.nzval[i] = rand_ver(ver)
    end
end

"""
    aprod = productGraph(a0, a1)

The Cartesian product of two graphs.  When applied to two paths, it gives a grid.
"""
function product_graph(a0::SparseMatrixCSC, a1::SparseMatrixCSC)
  n0 = size(a0)[1]
  n1 = size(a1)[1]
  a = kron(sparse(I,n0,n0),a1) + kron(a0,sparse(I, n1, n1));

end # productGraph

function product_graph(b::IJV{Tva,Tia}, a::IJV{Tvb,Tib}) where {Tva, Tvb, Tia, Tib}

    Ti = promote_type(Tia, Tib)

    n = a.n * b.n

    @assert length(a.i) == a.nnz

    a_edge_from = kron(ones(Ti, a.nnz), a.n*collect(0:(b.n-1)))
    ai = a_edge_from + kron(a.i, ones(Ti, b.n))
    aj = a_edge_from + kron(a.j, ones(Ti, b.n))
    av = kron(a.v, ones(b.n))

    b_edge_from = kron(collect(1:a.n), ones(Ti, b.nnz))
    bi = b_edge_from + kron(ones(Ti, a.n), (b.i .- 1) .* a.n)
    bj = b_edge_from + kron(ones(Ti, a.n), (b.j .- 1) .* a.n)
    bv = kron(ones(a.n), b.v)

    return IJV(n, length(av)+length(bv),
        [ai; bi], [aj; bj], [av; bv])
end

"""
    ap = power(a::SparseMatrixCSC, k::Int)

Returns the kth power of a.
"""
function power(a::SparseMatrixCSC, k::Int)
  ap = a^k
  ap = ap - sparse(Diagonal(diag(ap)))
end

"""
    a_new = thicken_once(a)

Creates one edge for every vertex in a of degree > 1
by connecting two of its random neighbors.
To use this to thicken a, return unweight(a + a_new).

```
a = grid2(5)
a2 = unweight(a + thicken_once(a))
(x,y) = grid2coords(5,5);
plotGraph(a2,x,y)
```
"""
function thicken_once(a::SparseMatrixCSC; ver = Vcur)
    n = a.n
    e_new = zeros(Int,n,2)
    ptr = 0
    for i in 1:n
        d = deg(a, i)
        if d > 1
            nbr1 = 0
            nbr2 = 0
            while nbr1 == nbr2
                nbr1 = rand_ver(ver, 1:d)
                nbr2 = rand_ver(ver, 1:d)
            end
            ptr += 1
            e_new[ptr,:] = [nbri(a,i,nbr1) nbri(a,i,nbr2)]
        end
    end

    a_new = sparse(e_new[1:ptr,1], e_new[1:ptr,2], 1, n, n)
    return unweight(a_new + a_new')
end

"""
    a_new = thicken(A,k)

Create a new graph with at least k times as many edges as A
By connecting nodes with common neighbors at random.
When this stops working (not enough new edges),
repeat on the most recently produced graph.
If k is too big, it is decreased so the average degree will not
be pushed much above n/2.

When called without k, it just runs thicken_once.

For example:
```
a = grid2(5)
a2 = thicken(a,3)
(x,y) = grid2coords(5,5);
plotGraph(a2,x,y)
```
"""
function thicken(a::SparseMatrixCSC,k; ver = Vcur)
    n = a.n

    k = min(k, round(Int, n^2/nnz(a)/2 ))

    ne0 = nnz(a)/2

    target = min(ne0*k, n*(n-1)/2)

    a_new = a

    max_tries = 5*k
    ntries = 0

    while nnz(a_new)/2 < target && ntries <= max_tries
        ntries += 1

        before = nnz(a_new)/2

        a_new = unweight(a_new + thicken_once(a, ver=ver))
        if nnz(a_new)/2 - before < n/4
            a = a_new
        end
    end
    return a_new

end

thicken(a; ver=Vcur) = unweight(a + thicken_once(a, ver=ver))

"""
    graph = generalizedNecklace(A, H, k::Int64)

Constructs a generalized necklace graph starting with two graphs A and H. The
resulting new graph will be constructed by expanding each vertex in H to an
instance of A. k random edges will be generated between components. Thus, the
resulting graph may have weighted edges.
"""
function generalizedNecklace(A::SparseMatrixCSC{Tv, Ti}, H::SparseMatrixCSC, k::Int64; ver=Vcur) where {Tv, Ti}
  a = findnz(A)
  h = findnz(H)

  # these are square matrices
  n = A.n
  m = H.n

  newI = Ti[]
  newJ = Ti[]
  newW = Tv[]

  # duplicate the vertices in A so that each vertex in H corresponds to a copy of A
  for i in 1:m
    newI = append!(newI, a[1] .+ n * (i - 1))
    newJ = append!(newJ, a[2] .+ n * (i - 1))
    newW = append!(newW, a[3])
  end

  # for each edge in H, add k random edges between two corresponding components
  # multiedges will be concatenated to a single edge with higher cost
  for i in 1:length(h[1])
    u = h[1][i]
    v = h[2][i]

    if (u < v)
      #component x is from 1 + (x - 1) * n to n + (x - 1) * n
      for edgeToAdd in 1:k
        newU = rand_ver(ver,1:n) + n * (u - 1)
        newV = rand_ver(ver,1:n) + n * (v - 1)
        append!(newI, [newU, newV])
        append!(newJ, [newV, newU])
        append!(newW, [1, 1])
      end
    end
  end

  return sparse(newI, newJ, newW)
end # generalizedNecklace

function generalized_necklace(A::IJV, H::IJV, k::Integer; ver=Vcur)

    if ver == V06
        H = IJV(sparse(H))
    end

    # these are square matrices
    n = A.n
    m = H.n

    # duplicate the vertices in A so that each vertex in H corresponds to a copy of A
    newW = kron(ones(m), A.v)
    newI = kron(ones(Int, m), A.i) + kron(n * collect(0:(m-1)), ones(Int, A.nnz) )
    newJ = kron(ones(Int, m), A.j) + kron(n * collect(0:(m-1)), ones(Int, A.nnz) )

    # for each edge in H, add k random edges between two corresponding components
    # multiedges will be concatenated to a single edge with higher weight

    Hi = Array{Int}(undef, div(k*H.nnz,2))
    Hj = Array{Int}(undef, div(k*H.nnz,2))
    Hv = ones(div(k*H.nnz,2))

    ind = 0
    for i in 1:H.nnz
      u = H.i[i]
      v = H.j[i]

      if (u < v)
        #component x is from 1 + (x - 1) * n to n + (x - 1) * n
        for j in 1:k
          newU = rand_ver(ver, 1:n) + n * (u - 1)
          newV = rand_ver(ver, 1:n) + n * (v - 1)

          ind += 1
          Hi[ind] = newU
          Hj[ind] = newV

        end
      end
    end

    return IJV(n*m, length(newI) + k*H.nnz, [newI; Hi; Hj] , [newJ; Hj; Hi], [newW; Hv; Hv])

  end # generalizedNecklace

"""
    U = edgeVertexMat(a)

The signed edge-vertex adjacency matrix
"""
function edgeVertexMat(mat::SparseMatrixCSC)
    (ai,aj) = findnz(triu(mat,1))
    m = length(ai)
    n = size(mat)[1]
    return sparse(collect(1:m),ai,1.0,m,n) - sparse(collect(1:m),aj,1.0,m,n)
end

"""
    H = line_graph(G::SparseMatrixCSC)

Let G = (V, E) be a graph. The line graph of G is the graph whose vertices are the edges of G in which two are connected if they share an endpoint in G.
That is, (u, v),(w, z) is an edge of the line graph if one of {u, v} is the same as one of {w, z}
"""
function line_graph(G::SparseMatrixCSC)
    u = abs.(edgeVertexMat(G));
    uu = u*u';
    return uu - Diagonal(diag(uu))
end

"""
    U = wtedEdgeVertexMat(a)

The signed and weighted edge-vertex adjacency matrix, so U'*U = L
"""
function wtedEdgeVertexMat(mat::SparseMatrixCSC)
    (ai,aj,av) = findnz(triu(mat,1))
    m = length(ai)
    n = size(mat)[1]
    v = av.^(1/2)
    return sparse(collect(1:m),ai,v,m,n) - sparse(collect(1:m),aj,v,m,n)
end

"""
    graph = subsampleEdges(a::SparseMatrixCSC, p::Float64)

Create a new graph from the old, but keeping edge edge with probability `p`
"""
function subsampleEdges(a::SparseMatrixCSC, p::Float64; ver=Vcur)
  (ai, aj, av) = findnz(triu(a))
  n = size(a)[1]
  mask = rand_ver(ver, length(ai)) .< p
  as = sparse(ai[mask],aj[mask],av[mask],n,n)
  as = as + as'
  return as

end # subsampleEdges

"""
    graph = two_lift(a, flip::AbstractArray{Bool,1})
    graph = two_lift(a)
    graph = two_lift(a, k::Integer)

Creats a 2-lift of a.  `flip` is a boolean indicating which edges cross.
In the third version, k is the number of edges that cross.
"""
function two_lift(a::SparseMatrixCSC, flip::AbstractArray{Bool,1})
    (ai,aj,av) = findnz(triu(a))
    m = length(ai)
    n = size(a)[1]
    a0 = sparse(ai[flip],aj[flip],1,n,n)
    a1 = sparse(ai[.!(flip)],aj[.!(flip)],1,n,n)
    a00 = a0 + a0'
    a11 = a1 + a1'
    return [a00 a11; a11 a00]
end

function two_lift(a::IJV{Tv,Ti}, flip::AbstractArray{Bool,1}) where {Tv,Ti}

    ii = Array{Ti}(undef, a.nnz)
    jj = Array{Ti}(undef, a.nnz)
    vv = Array{Tv}(undef, a.nnz)

    ind = 0
    cnt = 0

    for i in 1:a.nnz
        if (a.i[i] < a.j[i])

            cnt += 1
            if flip[cnt]
                ind += 1
                ii[ind] = a.i[i]
                jj[ind] = a.j[i]
                vv[ind] = a.v[i]

                ind += 1
                ii[ind] = a.i[i]+a.n
                jj[ind] = a.j[i]+a.n
                vv[ind] = a.v[i]
            else
                ind += 1
                ii[ind] = a.i[i]
                jj[ind] = a.j[i]+a.n
                vv[ind] = a.v[i]

                ind += 1
                ii[ind] = a.i[i]+a.n
                jj[ind] = a.j[i]
                vv[ind] = a.v[i]
            end
        end
    end

    return IJV(2*a.n, 2*a.nnz, [ii;jj], [jj;ii], [vv;vv])
end

two_lift(a; ver=Vcur) = two_lift(a,rand_ver(ver, false:true,div(nnz(a),2)))

two_lift(a, k::Integer; ver=Vcur) = two_lift(a,randperm_ver(ver, div(nnz(a),2)) .> k)

"""
    graph = joinGraphs(a, b, k::Integer)

 Create a disjoint union of graphs a and b,
 and then put k random edges between them
"""
function join_graphs(a::SparseMatrixCSC{Tval,Tind}, b::SparseMatrixCSC{Tval,Tind}, k::Integer; ver=Vcur) where {Tval,Tind}
    na = size(a)[1]
    nb = size(b)[1]

    (ai,aj,av) = findnz(a)
    (bi,bj,bv) = findnz(b)
    bi = bi .+ na
    bj = bj .+ na

    ji = rand_ver(ver, 1:na,k)
    jj = rand_ver(ver, 1:nb,k) .+ na

    ab = sparse([ai;bi;ji;jj],[aj;bj;jj;ji],[av;bv;ones(Tval,2*k)],na+nb,na+nb)
end

function join_graphs(a::IJV, b::IJV, k::Integer; ver=Vcur)

    ji = rand_ver(ver, 1:a.n,k)
    jj = rand_ver(ver, 1:b.n,k) .+ a.n

    gi = [a.i; b.i .+ a.n ; ji ; jj]
    gj = [a.j; b.j .+ a.n ; jj; ji]
    gv = [a.v ; b.v ; ones(2*k)]

    return IJV(a.n + b.n, a.nnz + b.nnz + 2*k, gi, gj, gv)
end

"""
    graph = join_graphs!(a::IJV, b::IJV, k::Integer)

 Create a disjoint union of graphs a and b,
 and then put k random edges between them, merging b into a.
"""
function join_graphs!(a::IJV, b::IJV, k::Integer; ver=Vcur)

    ji = rand_ver(ver, 1:a.n,k)
    jj = rand_ver(ver, 1:b.n,k) .+ a.n

    append!(a.i, [b.i .+ a.n ; ji ; jj] )
    append!(a.j, [b.j .+ a.n ; jj; ji])
    append!(a.v, [b.v ; ones(2*k)])
    a.n += b.n
    a.nnz += b.nnz + 2*k

end

"""
    graph = disjoin(a,b)

 Create a disjoint union of graphs a and b,
  with no edges between them.
"""
disjoin(a::SparseMatrixCSC,b::SparseMatrixCSC) = join_graphs(a,b,0)

disjoin(a::IJV, b::IJV) = IJV(a.n+b.n, a.nnz + b.nnz,
        [a.i ; b.i .+ a.n],
        [a.j ; b.j .+ a.n],
        [a.v;b.v])

function disjoin!(a::IJV, b::IJV)
    append!(a.i, b.i .+ a.n )
    append!(a.j, b.j .+ a.n )
    append!(a.v, b.v)
    a.n += b.n
    a.nnz += b.nnz
end

"""
    b = firstn(a::IJV, n::Integer)

Only keep the first n vertices of a.
"""
function firstn(a::IJV, n::Integer)
    mask = (a.i .<= n) .& (a.j .<= n)
    return IJV(n, sum(mask), a.i[mask], a.j[mask], a.v[mask])
end

"""
    plot_graph(gr,x,y;color=[0,0,1],dots=true,setaxis=true,number=false)

Plots graph gr with coordinates (x,y)
"""
function plot_graph(gr,x,y;color=[0,0,1],dots=true,setaxis=true,number=false)

    if isa(color, Vector) && length(color) == 3
        col = RGB(color...)
    else
        col = color
    end

    p = plot(;legend=false, axis=false, xticks=false, yticks=false)

    (ai,aj,av) = findnz(triu(gr))
    for i in 1:length(ai)
        s = [ai[i]; aj[i]]
        plot!(p, x[s], y[s], linecolor=col)
    end

    if dots
        scatter!(p,x, y , markercolor=col, markerstrokecolor=false)
    end

    if number
        annotate!(p, x, y, collect(1:length(x)))
    end

    minx = minimum(x)
    maxx = maximum(x)
    miny = minimum(y)
    maxy = maximum(y)
    delx = maxx - minx
    dely = maxy - miny

    plot!(p; ylims = (miny - dely/20, maxy + dely/20))
    plot!(p; xlims = (minx - delx/20, maxx + delx/20))

    display(p)
    return p
end

"""
    plot_graph(gr,x,y,z;color=[0,0,1],dots=true,setaxis=true,number=false)

Plots graph gr with coordinates (x,y,z)
"""
function plot_graph(gr,x,y,z;color=[0,0,1],dots=true,setaxis=true,number=false)

    if isa(color, Vector) && length(color) == 3
        col = RGB(color...)
    else
        col = color
    end

    p = plot(;legend=false, axis=false, xticks=false, yticks=false)

    (ai,aj,av) = findnz(triu(gr))
    for i in 1:length(ai)
        s = [ai[i]; aj[i]]
        plot!(p, x[s], y[s], z[s], linecolor=col)
    end

    if dots
        scatter!(p,x, y, z,  markercolor=col, markerstrokecolor=false)
    end

    if number
        annotate!(p, x, y, z,  collect(1:length(x)))
    end

    minx = minimum(x)
    maxx = maximum(x)
    miny = minimum(y)
    maxy = maximum(y)
    minz = minimum(z)
    maxz = maximum(z)
    delx = maxx - minx
    dely = maxy - miny
    delz = maxz - minz

    plot!(p; ylims = (miny - dely/20, maxy + dely/20))
    plot!(p; xlims = (minx - delx/20, maxx + delx/20))
    plot!(p; zlims = (minz - delz/20, maxz + delz/20))

    display(p)
    return p
end

"""
    spectral_drawing(a)

Computes spectral coordinates, and then uses plot_graph to draw
"""
function spectral_drawing(a)

    x, y = spectral_coords(a)
    plot_graph(a,x,y)

end # spectral_drawing

"""
    x, y = spectral_coords(a)

Computes the spectral coordinates of a graph.
If more than 2 coords are desired, you can use
~~~
    x, y, z = spectral_coords(a; k = 3)
~~~
"""
function spectral_coords(a; k = 2)

    E = eigs(lap(a), nev = (k+1), which=:SR)
    V = E[2]
    return tuple([V[:,i] for i in 2:(k+1)]...)

end # spectral_coords

"""
    d = diagmat(a)

Returns the diagonal weighted degree matrix(as a sparse matrix) of a graph
"""
function diagmat(a::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}

  return sparse(Diagonal(vec(sum(a,dims=1))))

end # diagmat


text:
mutable struct Tree <: GitHubType
    sha::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    tree::Union{Vector, Nothing}
    truncated::Union{Bool, Nothing}
end

Tree(data::Dict) = json2github(Tree, data)
namefield(tree::Tree) = tree.sha

@api_default function tree(api::GitHubAPI, repo, tree_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/git/trees/$(name(tree_obj))"; options...)
    return Tree(result)
end

@api_default function create_tree(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/git/trees"; options...)
    return Tree(result)
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks

# Benchmark function

function benchmark_time_stepper(Arch, N, TimeStepper)
    grid = RegularRectilinearGrid(size=(N, N, N), extent=(1, 1, 1))
    model = NonhydrostaticModel(architecture=Arch(), grid=grid, timestepper=TimeStepper)

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Ns = [192]
TimeSteppers = [:QuasiAdamsBashforth2, :RungeKutta3]

# Run and summarize benchmarks

print_system_info()
suite = run_benchmarks(benchmark_time_stepper; Architectures, Ns, TimeSteppers)

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :TimeSteppers, :Ns], by=(string, string, identity))
benchmarks_pretty_table(df, title="Time stepping benchmarks")

if GPU in Architectures
    df_Δ = gpu_speedups_suite(suite) |> speedups_dataframe
    sort!(df_Δ, [:TimeSteppers, :Ns], by=(string, identity))
    benchmarks_pretty_table(df_Δ, title="Time stepping CPU to GPU speedup")
end

for Arch in Architectures
    suite_arch = speedups_suite(suite[@tagged Arch], base_case=(Arch, Ns[1], :QuasiAdamsBashforth2))
    df_arch = speedups_dataframe(suite_arch, slowdown=true)
    sort!(df_arch, :TimeSteppers, by=string)
    benchmarks_pretty_table(df_arch, title="Time stepping relative performance ($Arch)")
end


text:
function FDSFactory(arg0::Field, arg1::jint, arg2::jint)
    return FDSFactory((Field, jint, jint), arg0, arg1, arg2)
end

function build(obj::FDSFactory, arg0::Vector{RealFieldElement})
    return jcall(obj, "build", FieldDerivativeStructure, (Vector{RealFieldElement},), arg0)
end

function build(obj::FDSFactory, arg0::Vector{jdouble})
    return jcall(obj, "build", FieldDerivativeStructure, (Vector{jdouble},), arg0)
end

function constant(obj::FDSFactory, arg0::RealFieldElement)
    return jcall(obj, "constant", FieldDerivativeStructure, (RealFieldElement,), arg0)
end

function constant(obj::FDSFactory, arg0::jdouble)
    return jcall(obj, "constant", FieldDerivativeStructure, (jdouble,), arg0)
end

function get_compiler(obj::FDSFactory)
    return jcall(obj, "getCompiler", DSCompiler, ())
end

function get_derivative_field(obj::FDSFactory)
    return jcall(obj, "getDerivativeField", Field, ())
end

function get_value_field(obj::FDSFactory)
    return jcall(obj, "getValueField", Field, ())
end

function variable(obj::FDSFactory, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "variable", FieldDerivativeStructure, (jint, RealFieldElement), arg0, arg1)
end

function variable(obj::FDSFactory, arg0::jint, arg1::jdouble)
    return jcall(obj, "variable", FieldDerivativeStructure, (jint, jdouble), arg0, arg1)
end


text:
module TestSyntax

using FGenerators
using Test

macro expansion_error(ex)
    quote
        err = try
            $Base.@eval $__module__ $ex
            nothing
        catch _err
            _err
        end
        $Test.@test err !== nothing
        err
    end
end

msgof(err) = sprint(showerror, err)

@testset begin
    err = @expansion_error @fgenerator(notype) do
    end
    @test occursin("requires `generator::GeneratorType`", msgof(err))

    err = @expansion_error @fgenerator(x::T) do arg
    end
    @test occursin("cannot be used with arguments after `do`", msgof(err))

    err = @expansion_error @fgenerator(function (; a); end, x::T)
    @test occursin("Got a function with keyword argument(s)", msgof(err))

    err = @expansion_error @fgenerator function ()
        @yield 1 2 3
    end
    @test occursin("`@yield` requires exactly one argument.", msgof(err))
end

end  # module


text:
@testset "1.1.2.6 (g x)^m (a+b x^2)^p (c+d x^2)^q (e+f x^2)^r" begin
    (A, B, a, b, c, d, e, m, p, x, ) = @variables A B a b c d e m p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*x)^m*(a+b*x^2)^p*(c+d*x^2)^q*(e+f*x^2)^r=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*x^2)^p*(c+d*x^2)^q*(A+B*x^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*x^2)^p*(c+d*x^2)^q*(A+B*x^2)=#

    #= ::Subsubsection::Closed:: =#
    #=q>0=#

    @test_int [(e*x)^m*(a + b*x^2)^3*(c + d*x^2)*(A + B*x^2), x, 2, (a^3*A*c*(e*x)^(1 + m))/(e*(1 + m)) + (a^2*(3*A*b*c + a*B*c + a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (a*(3*A*b*(b*c + a*d) + a*B*(3*b*c + a*d))*(e*x)^(5 + m))/(e^5*(5 + m)) + (b*(3*a*B*(b*c + a*d) + A*b*(b*c + 3*a*d))*(e*x)^(7 + m))/(e^7*(7 + m)) + (b^2*(b*B*c + A*b*d + 3*a*B*d)*(e*x)^(9 + m))/(e^9*(9 + m)) + (b^3*B*d*(e*x)^(11 + m))/(e^11*(11 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(c + d*x^2)*(A + B*x^2), x, 2, (a^2*A*c*(e*x)^(1 + m))/(e*(1 + m)) + (a*(2*A*b*c + a*B*c + a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + ((a*B*(2*b*c + a*d) + A*b*(b*c + 2*a*d))*(e*x)^(5 + m))/(e^5*(5 + m)) + (b*(b*B*c + A*b*d + 2*a*B*d)*(e*x)^(7 + m))/(e^7*(7 + m)) + (b^2*B*d*(e*x)^(9 + m))/(e^9*(9 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(c + d*x^2)*(A + B*x^2), x, 2, (a*A*c*(e*x)^(1 + m))/(e*(1 + m)) + ((A*b*c + a*B*c + a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + ((b*B*c + A*b*d + a*B*d)*(e*x)^(5 + m))/(e^5*(5 + m)) + (b*B*d*(e*x)^(7 + m))/(e^7*(7 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(c + d*x^2)*(A + B*x^2), x, 2, (A*c*(e*x)^(1 + m))/(e*(1 + m)) + ((B*c + A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (B*d*(e*x)^(5 + m))/(e^5*(5 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(c + d*x^2)*(A + B*x^2), x, 3, ((b*B*c + A*b*d - a*B*d)*(e*x)^(1 + m))/(b^2*e*(1 + m)) + (B*d*(e*x)^(3 + m))/(b*e^3*(3 + m)) + ((A*b - a*B)*(b*c - a*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*b^2*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(c + d*x^2)*(A + B*x^2), x, 3, -((d*(A*b*(1 + m) - a*B*(3 + m))*(e*x)^(1 + m))/(2*a*b^2*e*(1 + m))) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2))/(2*a*b*e*(a + b*x^2)) + ((a*B*(b*c*(1 + m) - a*d*(3 + m)) + A*b*(a*d*(1 + m) + b*(c - c*m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*b^2*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(c + d*x^2)*(A + B*x^2), x, 3, -(((A*b*(a*d*(1 - m) - b*c*(3 - m)) - a*B*(b*c*(1 + m) - a*d*(3 + m)))*(e*x)^(1 + m))/(8*a^2*b^2*e*(a + b*x^2))) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2))/(4*a*b*e*(a + b*x^2)^2) + ((A*b*(1 - m)*(b*c*(3 - m) + a*d*(1 + m)) + a*B*(1 + m)*(a*d*(3 + m) + b*(c - c*m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*b^2*e*(1 + m))]

    @test_int [(e*x)^m*(a + b*x^2)^3*(c + d*x^2)^2*(A + B*x^2), x, 2, (a^3*A*c^2*(e*x)^(1 + m))/(e*(1 + m)) + (a^2*c*(3*A*b*c + a*B*c + 2*a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (a*(a*B*c*(3*b*c + 2*a*d) + A*(3*b^2*c^2 + 6*a*b*c*d + a^2*d^2))*(e*x)^(5 + m))/(e^5*(5 + m)) + ((a*B*(3*b^2*c^2 + 6*a*b*c*d + a^2*d^2) + A*b*(b^2*c^2 + 6*a*b*c*d + 3*a^2*d^2))*(e*x)^(7 + m))/(e^7*(7 + m)) + (b*(3*a^2*B*d^2 + 3*a*b*d*(2*B*c + A*d) + b^2*c*(B*c + 2*A*d))*(e*x)^(9 + m))/(e^9*(9 + m)) + (b^2*d*(2*b*B*c + A*b*d + 3*a*B*d)*(e*x)^(11 + m))/(e^11*(11 + m)) + (b^3*B*d^2*(e*x)^(13 + m))/(e^13*(13 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(c + d*x^2)^2*(A + B*x^2), x, 2, (a^2*A*c^2*(e*x)^(1 + m))/(e*(1 + m)) + (a*c*(a*B*c + 2*A*(b*c + a*d))*(e*x)^(3 + m))/(e^3*(3 + m)) + ((2*a*B*c*(b*c + a*d) + A*(b^2*c^2 + 4*a*b*c*d + a^2*d^2))*(e*x)^(5 + m))/(e^5*(5 + m)) + ((a^2*B*d^2 + 2*a*b*d*(2*B*c + A*d) + b^2*c*(B*c + 2*A*d))*(e*x)^(7 + m))/(e^7*(7 + m)) + (b*d*(2*b*B*c + A*b*d + 2*a*B*d)*(e*x)^(9 + m))/(e^9*(9 + m)) + (b^2*B*d^2*(e*x)^(11 + m))/(e^11*(11 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(c + d*x^2)^2*(A + B*x^2), x, 2, (a*A*c^2*(e*x)^(1 + m))/(e*(1 + m)) + (c*(A*b*c + a*B*c + 2*a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + ((a*d*(2*B*c + A*d) + b*c*(B*c + 2*A*d))*(e*x)^(5 + m))/(e^5*(5 + m)) + (d*(2*b*B*c + A*b*d + a*B*d)*(e*x)^(7 + m))/(e^7*(7 + m)) + (b*B*d^2*(e*x)^(9 + m))/(e^9*(9 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(c + d*x^2)^2*(A + B*x^2), x, 2, (A*c^2*(e*x)^(1 + m))/(e*(1 + m)) + (c*(B*c + 2*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (d*(2*B*c + A*d)*(e*x)^(5 + m))/(e^5*(5 + m)) + (B*d^2*(e*x)^(7 + m))/(e^7*(7 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(c + d*x^2)^2*(A + B*x^2), x, 3, ((a^2*B*d^2 - a*b*d*(2*B*c + A*d) + b^2*c*(B*c + 2*A*d))*(e*x)^(1 + m))/(b^3*e*(1 + m)) + (d*(2*b*B*c + A*b*d - a*B*d)*(e*x)^(3 + m))/(b^2*e^3*(3 + m)) + (B*d^2*(e*x)^(5 + m))/(b*e^5*(5 + m)) + ((A*b - a*B)*(b*c - a*d)^2*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*b^3*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(c + d*x^2)^2*(A + B*x^2), x, 4, -((d*(A*b*(2*b*c*(1 + m) - a*d*(3 + m)) - a*B*(2*b*c*(3 + m) - a*d*(5 + m)))*(e*x)^(1 + m))/(2*a*b^3*e*(1 + m))) - (d^2*(A*b*(3 + m) - a*B*(5 + m))*(e*x)^(3 + m))/(2*a*b^2*e^3*(3 + m)) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2)^2)/(2*a*b*e*(a + b*x^2)) + ((b*c - a*d)*(a*B*(b*c*(1 + m) - a*d*(5 + m)) + A*b*(a*d*(3 + m) + b*(c - c*m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*b^3*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(c + d*x^2)^2*(A + B*x^2), x, 4, (d*(b*c*(1 + m) - a*d*(3 + m))*(A*b*(1 + m) - a*B*(5 + m))*(e*x)^(1 + m))/(8*a^2*b^3*e*(1 + m)) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2)^2)/(4*a*b*e*(a + b*x^2)^2) + ((b*c - a*d)*(e*x)^(1 + m)*(c*(A*b*(3 - m) + a*B*(1 + m)) - d*(A*b*(1 + m) - a*B*(5 + m))*x^2))/(8*a^2*b^2*e*(a + b*x^2)) - ((a*d*(b*c*(1 + m) - a*d*(3 + m))*(A*b*(1 + m) - a*B*(5 + m)) - b*c*(A*b*(3 - m) + a*B*(1 + m))*(a*d*(1 + m) + b*(c - c*m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*b^3*e*(1 + m))]

    @test_int [(e*x)^m*(a + b*x^2)^3*(c + d*x^2)^3*(A + B*x^2), x, 2, (a^3*A*c^3*(e*x)^(1 + m))/(e*(1 + m)) + (a^2*c^2*(a*B*c + 3*A*(b*c + a*d))*(e*x)^(3 + m))/(e^3*(3 + m)) + (3*a*c*(a*B*c*(b*c + a*d) + A*(b^2*c^2 + 3*a*b*c*d + a^2*d^2))*(e*x)^(5 + m))/(e^5*(5 + m)) + ((3*a*B*c*(b^2*c^2 + 3*a*b*c*d + a^2*d^2) + A*(b^3*c^3 + 9*a*b^2*c^2*d + 9*a^2*b*c*d^2 + a^3*d^3))*(e*x)^(7 + m))/(e^7*(7 + m)) + ((a^3*B*d^3 + 9*a*b^2*c*d*(B*c + A*d) + 3*a^2*b*d^2*(3*B*c + A*d) + b^3*c^2*(B*c + 3*A*d))*(e*x)^(9 + m))/(e^9*(9 + m)) + (3*b*d*(a^2*B*d^2 + b^2*c*(B*c + A*d) + a*b*d*(3*B*c + A*d))*(e*x)^(11 + m))/(e^11*(11 + m)) + (b^2*d^2*(3*b*B*c + A*b*d + 3*a*B*d)*(e*x)^(13 + m))/(e^13*(13 + m)) + (b^3*B*d^3*(e*x)^(15 + m))/(e^15*(15 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(c + d*x^2)^3*(A + B*x^2), x, 2, (a^2*A*c^3*(e*x)^(1 + m))/(e*(1 + m)) + (a*c^2*(2*A*b*c + a*B*c + 3*a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (c*(a*B*c*(2*b*c + 3*a*d) + A*(b^2*c^2 + 6*a*b*c*d + 3*a^2*d^2))*(e*x)^(5 + m))/(e^5*(5 + m)) + ((6*a*b*c*d*(B*c + A*d) + a^2*d^2*(3*B*c + A*d) + b^2*c^2*(B*c + 3*A*d))*(e*x)^(7 + m))/(e^7*(7 + m)) + (d*(a^2*B*d^2 + 3*b^2*c*(B*c + A*d) + 2*a*b*d*(3*B*c + A*d))*(e*x)^(9 + m))/(e^9*(9 + m)) + (b*d^2*(3*b*B*c + A*b*d + 2*a*B*d)*(e*x)^(11 + m))/(e^11*(11 + m)) + (b^2*B*d^3*(e*x)^(13 + m))/(e^13*(13 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(c + d*x^2)^3*(A + B*x^2), x, 2, (a*A*c^3*(e*x)^(1 + m))/(e*(1 + m)) + (c^2*(A*b*c + a*B*c + 3*a*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (c*(3*a*d*(B*c + A*d) + b*c*(B*c + 3*A*d))*(e*x)^(5 + m))/(e^5*(5 + m)) + (d*(3*b*c*(B*c + A*d) + a*d*(3*B*c + A*d))*(e*x)^(7 + m))/(e^7*(7 + m)) + (d^2*(3*b*B*c + A*b*d + a*B*d)*(e*x)^(9 + m))/(e^9*(9 + m)) + (b*B*d^3*(e*x)^(11 + m))/(e^11*(11 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(c + d*x^2)^3*(A + B*x^2), x, 2, (A*c^3*(e*x)^(1 + m))/(e*(1 + m)) + (c^2*(B*c + 3*A*d)*(e*x)^(3 + m))/(e^3*(3 + m)) + (3*c*d*(B*c + A*d)*(e*x)^(5 + m))/(e^5*(5 + m)) + (d^2*(3*B*c + A*d)*(e*x)^(7 + m))/(e^7*(7 + m)) + (B*d^3*(e*x)^(9 + m))/(e^9*(9 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(c + d*x^2)^3*(A + B*x^2), x, 3, -(((a^3*B*d^3 + 3*a*b^2*c*d*(B*c + A*d) - a^2*b*d^2*(3*B*c + A*d) - b^3*c^2*(B*c + 3*A*d))*(e*x)^(1 + m))/(b^4*e*(1 + m))) + (d*(a^2*B*d^2 + 3*b^2*c*(B*c + A*d) - a*b*d*(3*B*c + A*d))*(e*x)^(3 + m))/(b^3*e^3*(3 + m)) + (d^2*(3*b*B*c + A*b*d - a*B*d)*(e*x)^(5 + m))/(b^2*e^5*(5 + m)) + (B*d^3*(e*x)^(7 + m))/(b*e^7*(7 + m)) + ((A*b - a*B)*(b*c - a*d)^3*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*b^4*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(c + d*x^2)^3*(A + B*x^2), x, 4, -((d*(A*b*(3*b^2*c^2*(1 + m) - 3*a*b*c*d*(3 + m) + a^2*d^2*(5 + m)) - a*B*(3*b^2*c^2*(3 + m) - 3*a*b*c*d*(5 + m) + a^2*d^2*(7 + m)))*(e*x)^(1 + m))/(2*a*b^4*e*(1 + m))) - (d^2*(A*b*(3*b*c*(3 + m) - a*d*(5 + m)) - a*B*(3*b*c*(5 + m) - a*d*(7 + m)))*(e*x)^(3 + m))/(2*a*b^3*e^3*(3 + m)) - (d^3*(A*b*(5 + m) - a*B*(7 + m))*(e*x)^(5 + m))/(2*a*b^2*e^5*(5 + m)) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2)^3)/(2*a*b*e*(a + b*x^2)) + ((b*c - a*d)^2*(a*B*(b*c*(1 + m) - a*d*(7 + m)) + A*b*(a*d*(5 + m) + b*(c - c*m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*b^4*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(c + d*x^2)^3*(A + B*x^2), x, 5, -((d*(A*b*(2*b^2*c^2*(3 + 2*m - m^2) + 3*a*b*c*d*(3 + 4*m + m^2) - a^2*d^2*(15 + 8*m + m^2)) + a*B*(2*b^2*c^2*(1 + m)^2 - 3*a*b*c*d*(15 + 8*m + m^2) + a^2*d^2*(35 + 12*m + m^2)))*(e*x)^(1 + m))/(8*a^2*b^4*e*(1 + m))) - (d^2*(A*b*(3 + m)*(b*c*(3 - m) + a*d*(5 + m)) + a*B*(b*c*(3 + 4*m + m^2) - a*d*(35 + 12*m + m^2)))*(e*x)^(3 + m))/(8*a^2*b^3*e^3*(3 + m)) + ((A*b*(b*c*(3 - m) + a*d*(3 + m)) + a*B*(b*c*(1 + m) - a*d*(7 + m)))*(e*x)^(1 + m)*(c + d*x^2)^2)/(8*a^2*b^2*e*(a + b*x^2)) + ((A*b - a*B)*(e*x)^(1 + m)*(c + d*x^2)^3)/(4*a*b*e*(a + b*x^2)^2) + ((b*c - a*d)*(A*b*(2*a*b*c*d*(3 - 2*m - m^2) + b^2*c^2*(3 - 4*m + m^2) + a^2*d^2*(15 + 8*m + m^2)) + a*B*(b^2*c^2*(1 - m^2) + 2*a*b*c*d*(5 + 6*m + m^2) - a^2*d^2*(35 + 12*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*b^4*e*(1 + m))]

    #= ::Subsubsection::Closed:: =#
    #=q<0=#

    @test_int [(e*x)^m*(a + b*x^2)^4*(A + B*x^2)/(c + d*x^2), x, 3, ((a^4*B*d^4 + b^4*c^3*(B*c - A*d) - 4*a*b^3*c^2*d*(B*c - A*d) + 6*a^2*b^2*c*d^2*(B*c - A*d) - 4*a^3*b*d^3*(B*c - A*d))*(e*x)^(1 + m))/(d^5*e*(1 + m)) + (b*(4*a^3*B*d^3 - b^3*c^2*(B*c - A*d) + 4*a*b^2*c*d*(B*c - A*d) - 6*a^2*b*d^2*(B*c - A*d))*(e*x)^(3 + m))/(d^4*e^3*(3 + m)) + (b^2*(6*a^2*B*d^2 + b^2*c*(B*c - A*d) - 4*a*b*d*(B*c - A*d))*(e*x)^(5 + m))/(d^3*e^5*(5 + m)) - (b^3*(b*B*c - A*b*d - 4*a*B*d)*(e*x)^(7 + m))/(d^2*e^7*(7 + m)) + (b^4*B*(e*x)^(9 + m))/(d*e^9*(9 + m)) - ((b*c - a*d)^4*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*d^5*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2), x, 3, ((a^3*B*d^3 - b^3*c^2*(B*c - A*d) + 3*a*b^2*c*d*(B*c - A*d) - 3*a^2*b*d^2*(B*c - A*d))*(e*x)^(1 + m))/(d^4*e*(1 + m)) + (b*(3*a^2*B*d^2 + b^2*c*(B*c - A*d) - 3*a*b*d*(B*c - A*d))*(e*x)^(3 + m))/(d^3*e^3*(3 + m)) - (b^2*(b*B*c - A*b*d - 3*a*B*d)*(e*x)^(5 + m))/(d^2*e^5*(5 + m)) + (b^3*B*(e*x)^(7 + m))/(d*e^7*(7 + m)) + ((b*c - a*d)^3*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*d^4*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2), x, 3, ((a^2*B*d^2 + b^2*c*(B*c - A*d) - 2*a*b*d*(B*c - A*d))*(e*x)^(1 + m))/(d^3*e*(1 + m)) - (b*(b*B*c - A*b*d - 2*a*B*d)*(e*x)^(3 + m))/(d^2*e^3*(3 + m)) + (b^2*B*(e*x)^(5 + m))/(d*e^5*(5 + m)) - ((b*c - a*d)^2*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*d^3*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2), x, 3, -(((b*B*c - A*b*d - a*B*d)*(e*x)^(1 + m))/(d^2*e*(1 + m))) + (b*B*(e*x)^(3 + m))/(d*e^3*(3 + m)) + ((b*c - a*d)*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*d^2*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(A + B*x^2)/(c + d*x^2), x, 2, (B*(e*x)^(1 + m))/(d*e*(1 + m)) - ((B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (1 + m + 2)/2, -((d*x^2)/c)))/(c*d*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2), x, 4, ((A*b - a*B)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (1 + m + 2)/2, -((b*x^2)/a)))/(a*(b*c - a*d)*e*(1 + m)) + ((B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (1 + m + 2)/2, -((d*x^2)/c)))/(c*(b*c - a*d)*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2), x, 5, ((A*b - a*B)*(e*x)^(1 + m))/(2*a*(b*c - a*d)*e*(a + b*x^2)) + ((A*b*(b*c*(1 - m) - a*d*(3 - m)) + a*B*(a*d*(1 - m) + b*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*(b*c - a*d)^2*e*(1 + m)) - (d*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*(b*c - a*d)^2*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2), x, 6, ((A*b - a*B)*(e*x)^(1 + m))/(4*a*(b*c - a*d)*e*(a + b*x^2)^2) + ((A*b*(b*c*(3 - m) - a*d*(7 - m)) + a*B*(a*d*(3 - m) + b*c*(1 + m)))*(e*x)^(1 + m))/(8*a^2*(b*c - a*d)^2*e*(a + b*x^2)) + ((A*b*(a^2*d^2*(15 - 8*m + m^2) - 2*a*b*c*d*(5 - 6*m + m^2) + b^2*c^2*(3 - 4*m + m^2)) + a*B*(b^2*c^2*(1 - m^2) - 2*a*b*c*d*(3 + 2*m - m^2) - a^2*d^2*(3 - 4*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*(b*c - a*d)^3*e*(1 + m)) + (d^2*(B*c - A*d)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(c*(b*c - a*d)^3*e*(1 + m))]

    @test_int [(e*x)^m*(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2)^2, x, 4, -((b*(3*a^2*d^2*(A*d*(1 + m) - B*c*(3 + m)) - 3*a*b*c*d*(A*d*(3 + m) - B*c*(5 + m)) + b^2*c^2*(A*d*(5 + m) - B*c*(7 + m)))*(e*x)^(1 + m))/(2*c*d^4*e*(1 + m))) - (b^2*(3*a*d*(A*d*(3 + m) - B*c*(5 + m)) - b*c*(A*d*(5 + m) - B*c*(7 + m)))*(e*x)^(3 + m))/(2*c*d^3*e^3*(3 + m)) - (b^3*(A*d*(5 + m) - B*c*(7 + m))*(e*x)^(5 + m))/(2*c*d^2*e^5*(5 + m)) - ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^3)/(2*c*d*e*(c + d*x^2)) + ((b*c - a*d)^2*(a*d*(A*d*(1 - m) + B*c*(1 + m)) + b*c*(A*d*(5 + m) - B*c*(7 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*d^4*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2)^2, x, 4, -((b*(2*a*d*(A*d*(1 + m) - B*c*(3 + m)) - b*c*(A*d*(3 + m) - B*c*(5 + m)))*(e*x)^(1 + m))/(2*c*d^3*e*(1 + m))) - (b^2*(A*d*(3 + m) - B*c*(5 + m))*(e*x)^(3 + m))/(2*c*d^2*e^3*(3 + m)) - ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^2)/(2*c*d*e*(c + d*x^2)) - ((b*c - a*d)*(a*d*(A*d*(1 - m) + B*c*(1 + m)) + b*c*(A*d*(3 + m) - B*c*(5 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*d^3*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2)^2, x, 3, -((B*(a*d*(1 + m) - b*c*(3 + m))*(e*x)^(1 + m))/(2*c*d^2*e*(1 + m))) - ((b*c - a*d)*(e*x)^(1 + m)*(A + B*x^2))/(2*c*d*e*(c + d*x^2)) + ((a*d*(A*d*(1 - m) + B*c*(1 + m)) + b*c*(A*d*(1 + m) - B*c*(3 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*d^2*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(A + B*x^2)/(c + d*x^2)^2, x, 2, -(((B*c - A*d)*(e*x)^(1 + m))/(2*c*d*e*(c + d*x^2))) + ((A*d*(1 - m) + B*c*(1 + m))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*d*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2)^2, x, 5, ((B*c - A*d)*(e*x)^(1 + m))/(2*c*(b*c - a*d)*e*(c + d*x^2)) + (b*(A*b - a*B)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*(b*c - a*d)^2*e*(1 + m)) + ((b*c*(B*c*(1 - m) - A*d*(3 - m)) + a*d*(A*d*(1 - m) + B*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*(b*c - a*d)^2*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2)^2, x, 6, (d*(A*b*c - 2*a*B*c + a*A*d)*(e*x)^(1 + m))/(2*a*c*(b*c - a*d)^2*e*(c + d*x^2)) + ((A*b - a*B)*(e*x)^(1 + m))/(2*a*(b*c - a*d)*e*(a + b*x^2)*(c + d*x^2)) + (b*(A*b*(b*c*(1 - m) - a*d*(5 - m)) + a*B*(a*d*(3 - m) + b*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*(b*c - a*d)^3*e*(1 + m)) - (d*(b*c*(B*c*(3 - m) - A*d*(5 - m)) + a*d*(A*d*(1 - m) + B*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*(b*c - a*d)^3*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2)^2, x, 7, -((d*(A*(4*a^2*d^2 - b^2*c^2*(3 - m) + a*b*c*d*(11 - m)) - a*B*c*(a*d*(11 - m) + b*c*(1 + m)))*(e*x)^(1 + m))/(8*a^2*c*(b*c - a*d)^3*e*(c + d*x^2))) + ((A*b - a*B)*(e*x)^(1 + m))/(4*a*(b*c - a*d)*e*(a + b*x^2)^2*(c + d*x^2)) + ((A*b*(b*c*(3 - m) - a*d*(9 - m)) + a*B*(a*d*(5 - m) + b*c*(1 + m)))*(e*x)^(1 + m))/(8*a^2*(b*c - a*d)^2*e*(a + b*x^2)*(c + d*x^2)) + (b*(a*B*(b^2*c^2*(1 - m^2) - 2*a*b*c*d*(5 + 4*m - m^2) - a^2*d^2*(15 - 8*m + m^2)) + A*b*(a^2*d^2*(35 - 12*m + m^2) - 2*a*b*c*d*(7 - 8*m + m^2) + b^2*c^2*(3 - 4*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*(b*c - a*d)^4*e*(1 + m)) + (d^2*(b*c*(B*c*(5 - m) - A*d*(7 - m)) + a*d*(A*d*(1 - m) + B*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(2*c^2*(b*c - a*d)^4*e*(1 + m))]

    @test_int [(e*x)^m*(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2)^3, x, 5, -((b*(2*a^2*d^2*(1 + m)*(A*d*(3 - m) + B*c*(1 + m)) + 3*a*b*c*d*(3 + m)*(A*d*(1 + m) - B*c*(5 + m)) - b^2*c^2*(5 + m)*(A*d*(3 + m) - B*c*(7 + m)))*(e*x)^(1 + m))/(8*c^2*d^4*e*(1 + m))) - (b^2*(a*d*(3 + m)*(A*d*(3 - m) + B*c*(1 + m)) + b*c*(5 + m)*(A*d*(3 + m) - B*c*(7 + m)))*(e*x)^(3 + m))/(8*c^2*d^3*e^3*(3 + m)) - ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^3)/(4*c*d*e*(c + d*x^2)^2) + ((a*d*(A*d*(3 - m) + B*c*(1 + m)) + b*c*(A*d*(3 + m) - B*c*(7 + m)))*(e*x)^(1 + m)*(a + b*x^2)^2)/(8*c^2*d^2*e*(c + d*x^2)) - ((b*c - a*d)*(a^2*d^2*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) + b^2*c^2*(5 + m)*(A*d*(3 + m) - B*c*(7 + m)) + 2*a*b*c*d*(A*d*(3 - 2*m - m^2) + B*c*(5 + 6*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*d^4*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2)^3, x, 4, (b*(a*d*(1 + m) - b*c*(3 + m))*(A*d*(1 + m) - B*c*(5 + m))*(e*x)^(1 + m))/(8*c^2*d^3*e*(1 + m)) - ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^2)/(4*c*d*e*(c + d*x^2)^2) - ((b*c - a*d)*(e*x)^(1 + m)*(a*(A*d*(3 - m) + B*c*(1 + m)) - b*(A*d*(1 + m) - B*c*(5 + m))*x^2))/(8*c^2*d^2*e*(c + d*x^2)) + ((a*d*(a*d*(1 - m) + b*c*(1 + m))*(A*d*(3 - m) + B*c*(1 + m)) - b*c*(a*d*(1 + m) - b*c*(3 + m))*(A*d*(1 + m) - B*c*(5 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*d^3*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2)^3, x, 3, -(((b*c - a*d)*(e*x)^(1 + m)*(A + B*x^2))/(4*c*d*e*(c + d*x^2)^2)) + ((b*c*(A*d*(1 + m) - B*c*(3 + m)) + a*d*(A*d*(3 - m) - B*(c - c*m)))*(e*x)^(1 + m))/(8*c^2*d^2*e*(c + d*x^2)) + ((a*d*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) + b*c*(1 + m)*(A*d*(1 - m) + B*c*(3 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*d^2*e*(1 + m))]
    @test_int [(e*x)^m*(a + b*x^2)^0*(A + B*x^2)/(c + d*x^2)^3, x, 2, -(((B*c - A*d)*(e*x)^(1 + m))/(4*c*d*e*(c + d*x^2)^2)) + ((A*d*(3 - m) + B*c*(1 + m))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(2, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(4*c^3*d*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^1*(A + B*x^2)/(c + d*x^2)^3, x, 6, ((B*c - A*d)*(e*x)^(1 + m))/(4*c*(b*c - a*d)*e*(c + d*x^2)^2) + ((b*c*(B*c*(3 - m) - A*d*(7 - m)) + a*d*(A*d*(3 - m) + B*c*(1 + m)))*(e*x)^(1 + m))/(8*c^2*(b*c - a*d)^2*e*(c + d*x^2)) + (b^2*(A*b - a*B)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*(b*c - a*d)^3*e*(1 + m)) + ((b^2*c^2*(B*c*(1 - m) - A*d*(5 - m))*(3 - m) - a^2*d^2*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) + 2*a*b*c*d*(B*c*(3 + 2*m - m^2) + A*d*(5 - 6*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*(b*c - a*d)^3*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^2*(A + B*x^2)/(c + d*x^2)^3, x, 7, (d*(2*A*b*c - 3*a*B*c + a*A*d)*(e*x)^(1 + m))/(4*a*c*(b*c - a*d)^2*e*(c + d*x^2)^2) + ((A*b - a*B)*(e*x)^(1 + m))/(2*a*(b*c - a*d)*e*(a + b*x^2)*(c + d*x^2)^2) + (d*(A*(4*b^2*c^2 - a^2*d^2*(3 - m) + a*b*c*d*(11 - m)) - a*B*c*(b*c*(11 - m) + a*d*(1 + m)))*(e*x)^(1 + m))/(8*a*c^2*(b*c - a*d)^3*e*(c + d*x^2)) + (b^2*(A*b*(b*c*(1 - m) - a*d*(7 - m)) + a*B*(a*d*(5 - m) + b*c*(1 + m)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(2*a^2*(b*c - a*d)^4*e*(1 + m)) - (d*(b^2*c^2*(B*c*(3 - m) - A*d*(7 - m))*(5 - m) - a^2*d^2*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) + 2*a*b*c*d*(B*c*(5 + 4*m - m^2) + A*d*(7 - 8*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*(b*c - a*d)^4*e*(1 + m))]
    @test_int [(e*x)^m/(a + b*x^2)^3*(A + B*x^2)/(c + d*x^2)^3, x, 8, -((d*(A*(2*a^2*d^2 - b^2*c^2*(3 - m) + a*b*c*d*(13 - m)) - a*B*c*(a*d*(11 - m) + b*c*(1 + m)))*(e*x)^(1 + m))/(8*a^2*c*(b*c - a*d)^3*e*(c + d*x^2)^2)) + ((A*b - a*B)*(e*x)^(1 + m))/(4*a*(b*c - a*d)*e*(a + b*x^2)^2*(c + d*x^2)^2) + ((A*b*(b*c*(3 - m) - a*d*(11 - m)) + a*B*(a*d*(7 - m) + b*c*(1 + m)))*(e*x)^(1 + m))/(8*a^2*(b*c - a*d)^2*e*(a + b*x^2)*(c + d*x^2)^2) + (d*(A*(b*c + a*d)*(b^2*c^2*(3 - m) + a^2*d^2*(3 - m) - 2*a*b*c*d*(9 - m)) + a*B*c*(2*a*b*c*d*(11 - m) + b^2*c^2*(1 + m) + a^2*d^2*(1 + m)))*(e*x)^(1 + m))/(8*a^2*c^2*(b*c - a*d)^4*e*(c + d*x^2)) + (b^2*(a*B*(b^2*c^2*(1 - m^2) - 2*a*b*c*d*(7 + 6*m - m^2) - a^2*d^2*(35 - 12*m + m^2)) + A*b*(a^2*d^2*(63 - 16*m + m^2) - 2*a*b*c*d*(9 - 10*m + m^2) + b^2*c^2*(3 - 4*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(8*a^3*(b*c - a*d)^5*e*(1 + m)) + (d^2*(b^2*c^2*(B*c*(5 - m) - A*d*(9 - m))*(7 - m) - a^2*d^2*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) + 2*a*b*c*d*(B*c*(7 + 6*m - m^2) + A*d*(9 - 10*m + m^2)))*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((d*x^2)/c)))/(8*c^3*(b*c - a*d)^5*e*(1 + m))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*x^2)^p*(c+d*x^2)^q*(A+B*x^2)*with*p*symbolic=#

    @test_int [(e*x)^m*(a + b*x^2)^p*(c + d*x^2)^3*(A + B*x^2), x, 6, -((1/(b^4*e*(3 + m + 2*p)*(5 + m + 2*p)*(7 + m + 2*p)*(9 + m + 2*p)))*((a^3*B*d^3*(105 + 71*m + 15*m^2 + m^3) - a^2*b*d^2*(5 + m)*(A*d*(3 + m)*(9 + m + 2*p) + 2*B*c*(30 + 13*m + m^2 + 2*p + 2*m*p)) + a*b^2*c*d*(2*A*d*(216 + m^3 + 84*p + 8*p^2 + 4*m^2*(5 + p) + m*(123 + 44*p + 4*p^2)) + B*c*(267 + m^3 + 40*p + 4*p^2 + m^2*(21 + 4*p) + m*(143 + 44*p + 4*p^2))) - b^3*c^2*(48*B*c + A*d*(513 + m^3 + 366*p + 92*p^2 + 8*p^3 + m^2*(23 + 6*p) + m*(183 + 92*p + 12*p^2))))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))) + ((a^2*B*d^2*(35 + 12*m + m^2) + b^2*c*(24*B*c + A*d*(99 + m^2 + 40*p + 4*p^2 + 4*m*(5 + p))) - a*b*d*(A*d*(5 + m)*(9 + m + 2*p) + B*c*(65 + m^2 + 2*p + 2*m*(9 + p))))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(c + d*x^2))/(b^3*e*(5 + m + 2*p)*(7 + m + 2*p)*(9 + m + 2*p)) - ((a*B*d*(7 + m) - b*(6*B*c + A*d*(9 + m + 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(c + d*x^2)^2)/(b^2*e*(7 + m + 2*p)*(9 + m + 2*p)) + (B*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(c + d*x^2)^3)/(b*e*(9 + m + 2*p)) + (1/(b^4*e*(1 + m)*(3 + m + 2*p)*(5 + m + 2*p)*(7 + m + 2*p)*(9 + m + 2*p)))*(((a*(1 + m)*(a^3*B*d^3*(105 + 71*m + 15*m^2 + m^3) - a^2*b*d^2*(5 + m)*(A*d*(3 + m)*(9 + m + 2*p) + 2*B*c*(30 + 13*m + m^2 + 2*p + 2*m*p)) + a*b^2*c*d*(2*A*d*(216 + m^3 + 84*p + 8*p^2 + 4*m^2*(5 + p) + m*(123 + 44*p + 4*p^2)) + B*c*(267 + m^3 + 40*p + 4*p^2 + m^2*(21 + 4*p) + m*(143 + 44*p + 4*p^2))) - b^3*c^2*(48*B*c + A*d*(513 + m^3 + 366*p + 92*p^2 + 8*p^3 + m^2*(23 + 6*p) + m*(183 + 92*p + 12*p^2)))) - b*c*(3 + m + 2*p)*(2*b*c*(2 + p)*(2*b*c*(3 + p)*(a*B*(1 + m) - A*b*(9 + m + 2*p)) + (b*c - a*d)*(1 + m)*(a*B*(7 + m) - A*b*(9 + m + 2*p))) + (1 + m)*(b*c*(2*b*c*(3 + p)*(a*B*(1 + m) - A*b*(9 + m + 2*p)) + (b*c - a*d)*(1 + m)*(a*B*(7 + m) - A*b*(9 + m + 2*p))) - a*(2*b*c*d*(3 + p)*(a*B*(1 + m) - A*b*(9 + m + 2*p)) + d*(b*c - a*d)*(1 + m)*(a*B*(7 + m) - A*b*(9 + m + 2*p)) + 4*(b*c - a*d)*(a*B*d*(7 + m) - b*(6*B*c + A*d*(9 + m + 2*p)))))))*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/(1 + (b*x^2)/a)^p)]
    @test_int [(e*x)^m*(a + b*x^2)^p*(c + d*x^2)^2*(A + B*x^2), x, 5, ((a^2*B*d^2*(15 + 8*m + m^2) + b^2*c*(8*B*c + A*d*(7 + m + 2*p)^2) - a*b*d*(A*d*(3 + m)*(7 + m + 2*p) + B*c*(27 + m^2 + 2*p + 2*m*(6 + p))))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))/(b^3*e*(3 + m + 2*p)*(5 + m + 2*p)*(7 + m + 2*p)) - ((a*B*d*(5 + m) - b*(4*B*c + A*d*(7 + m + 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(c + d*x^2))/(b^2*e*(5 + m + 2*p)*(7 + m + 2*p)) + (B*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(c + d*x^2)^2)/(b*e*(7 + m + 2*p)) - (1/(b^3*e*(1 + m)*(3 + m + 2*p)*(5 + m + 2*p)*(7 + m + 2*p)))*(((b*c*(3 + m + 2*p)*(2*b*c*(2 + p)*(a*B*(1 + m) - A*b*(7 + m + 2*p)) + (b*c - a*d)*(1 + m)*(a*B*(5 + m) - A*b*(7 + m + 2*p))) - a*(1 + m)*(2*b*c*d*(2 + p)*(a*B*(1 + m) - A*b*(7 + m + 2*p)) + d*(b*c - a*d)*(1 + m)*(a*B*(5 + m) - A*b*(7 + m + 2*p)) + 2*(b*c - a*d)*(a*B*d*(5 + m) - b*(4*B*c + A*d*(7 + m + 2*p)))))*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/(1 + (b*x^2)/a)^p)]
    @test_int [(e*x)^m*(a + b*x^2)^p*(c + d*x^2)^1*(A + B*x^2), x, 4, -(((a*B*d*(3 + m) - b*(2*A*d + B*c*(5 + m + 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))/(b^2*e*(3 + m + 2*p)*(5 + m + 2*p))) + (d*(e*x)^(1 + m)*(a + b*x^2)^(1 + p)*(A + B*x^2))/(b*e*(5 + m + 2*p)) - ((A*b*(3 + m + 2*p)*(a*d*(1 + m) - b*c*(5 + m + 2*p)) - a*(1 + m)*(a*B*d*(3 + m) - b*(2*A*d + B*c*(5 + m + 2*p))))*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/(1 + (b*x^2)/a)^p/(b^2*e*(1 + m)*(3 + m + 2*p)*(5 + m + 2*p))]
    @test_int [(e*x)^m*(a + b*x^2)^p*(A + B*x^2)/(c + d*x^2)^1, x, 6, -(((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^p*AppellF1((1 + m)/2, -p, 1, (3 + m)/2, -((b*x^2)/a), -((d*x^2)/c)))/((1 + (b*x^2)/a)^p*(c*d*e*(1 + m)))) + (B*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/((1 + (b*x^2)/a)^p*(d*e*(1 + m)))]
    @test_int [(e*x)^m*(a + b*x^2)^p*(A + B*x^2)/(c + d*x^2)^2, x, 7, ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))/(2*c*(b*c - a*d)*e*(c + d*x^2)) - ((a*d*(A*d*(1 - m) + B*c*(1 + m)) - b*c*(A*d*(1 - m - 2*p) + B*c*(1 + m + 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^p*AppellF1((1 + m)/2, -p, 1, (3 + m)/2, -((b*x^2)/a), -((d*x^2)/c)))/((1 + (b*x^2)/a)^p*(2*c^2*d*(b*c - a*d)*e*(1 + m))) - (b*(B*c - A*d)*(1 + m + 2*p)*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/((1 + (b*x^2)/a)^p*(2*c*d*(b*c - a*d)*e*(1 + m)))]
    @test_int [(e*x)^m*(a + b*x^2)^p*(A + B*x^2)/(c + d*x^2)^3, x, 8, ((B*c - A*d)*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))/(4*c*(b*c - a*d)*e*(c + d*x^2)^2) + ((a*d*(A*d*(3 - m) + B*c*(1 + m)) + b*c*(B*c*(1 - m - 2*p) - A*d*(5 - m - 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^(1 + p))/(8*c^2*(b*c - a*d)^2*e*(c + d*x^2)) + (1/(8*c^3*d*(b*c - a*d)^2*e*(1 + m)))*(((a^2*d^2*(1 - m)*(A*d*(3 - m) + B*c*(1 + m)) - 2*a*b*c*d*(B*c*(1 + m)*(1 - m - 2*p) + A*d*(1 - m)*(3 - m - 2*p)) + b^2*c^2*(1 - m - 2*p)*(A*d*(3 - m - 2*p) + B*c*(1 + m + 2*p)))*(e*x)^(1 + m)*(a + b*x^2)^p*AppellF1((1 + m)/2, -p, 1, (3 + m)/2, -((b*x^2)/a), -((d*x^2)/c)))/(1 + (b*x^2)/a)^p) - (b*(a*d*(A*d*(3 - m) + B*c*(1 + m)) + b*c*(B*c*(1 - m - 2*p) - A*d*(5 - m - 2*p)))*(1 + m + 2*p)*(e*x)^(1 + m)*(a + b*x^2)^p*HypergeometricFunctions._₂F₁((1 + m)/2, -p, (3 + m)/2, -((b*x^2)/a)))/((1 + (b*x^2)/a)^p*(8*c^2*d*(b*c - a*d)^2*e*(1 + m)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*x^2)^(p/2)*(c+d*x^2)^q*(A+B*x^2)=#

    @test_int [(sqrt(a + b*x^2)*(A + B*x^2)*(c + d*x^2))/x, x, 5, A*c*sqrt(a + b*x^2) - ((a + b*x^2)^(3/2)*(2*a*B*d - 5*b*(B*c + A*d) - 3*b*B*d*x^2))/(15*b^2) - sqrt(a)*A*c*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [((a + b*x^2)*(A + B*x^2)*sqrt(c + d*x^2))/x, x, 5, a*A*sqrt(c + d*x^2) - ((c + d*x^2)^(3/2)*(2*b*B*c - 5*(A*b + a*B)*d - 3*b*B*d*x^2))/(15*d^2) - a*A*sqrt(c)*atanh(sqrt(c + d*x^2)/sqrt(c))]
    end

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.6 (g x)^m (a+b x^2)^p (c+d x^2)^q (e+f x^2)^r.jl


text:
# # Porting Yao with Zygote
using Zygote, YaoExtensions, Yao
using Random; Random.seed!(2)

# first include this patch file in your working directory.
include("zygote_patch.jl")

# ## setup a circuit
# wrap differentiable nodes with `Diff`, please try `print(c)` to see the marked nodes with `[∂]`.
c = variational_circuit(5) |> autodiff(:BP)
dispatch!(c, :random)

# ## define a loss function
h = rand_hermitian(1<<5)
function loss2(reg::ArrayReg, params::Vector{T}, circuit) where T
    reg = dispatch_diff_apply!(copy(reg), circuit, params)
    v = vec(reg.state)
    v' * h *v |> real
end

# ## obtaining the gradient
reg = rand_state(5)
params = rand!(parameters_of_diff(c))
grad_input, grad_params, _ = Zygote.gradient(loss2, reg, params, c)

# ## Check gradients
using Test

# to check gradients of input wave function, we can introduce a one-hot vector
# perturbation to get numerical finite-difference gradient.
η = 1e-5
perturb_input_real = η*product_state(bit"00111")
perturb_input_imag = im*η*product_state(bit"00111")

num_gradient = (loss2(reg + perturb_input_real, params, c) - loss2(reg - perturb_input_real, params, c)) / 2η +
                im*(loss2(reg + perturb_input_imag, params, c) - loss2(reg - perturb_input_imag, params, c)) / 2η
exact_gradient = select(grad_input, bit"00111").state[]
@test num_gradient ≈ exact_gradient

# check the gradient of circuit parameters
index = 1
perturb_params = (p=zeros(length(params)); p[index]=1; p)
loss_pos = loss2(reg, params+η*perturb_params, c)
loss_neg = loss2(reg, params-η*perturb_params, c)
@test grad_params[index] ≈ (loss_pos - loss_neg) / 2η


text:
module SimpleMultigrid

# dependencies
using IterativeSolvers, LinearAlgebra, StaticArrays, SparseArrays, Primes, Printf

# export statements
export laplace1d, laplace2d, laplace3d, elliptic1d, elliptic2d, elliptic3d
export Injection, FullWeighting, Cubic
export GaussSeidel, Jacobi, RedBlackGaussSeidel
export coarsen
export MultigridMethod, V, W, F, V_cycle, W_cycle, F_cycle, \, size

# import statements
import Base: expand, show, iterate, \, size

# include source files
include("compose_matrix.jl")
include("grid_transfer_operators.jl")
include("iterative_methods.jl")
include("grids.jl")
include("multigrid.jl")

end # module


text:
"""
    ismasked(data, I...)

Check if a cell is masked, using the `mask` array.

Used used internally during simulations to skip masked cells.

If `mask` was not passed to the `Output` constructor or `sim!`
it defaults to `nothing` and `false` is always returned.
"""
ismasked(data::AbstractSimData, I...) = ismasked(mask(data), I...)
ismasked(data::GridData, I...) = ismasked(mask(data), I...)
ismasked(mask::Nothing, I...) = false
ismasked(mask::AbstractArray, I...) = @inbounds !(mask[I...])

"""
    isinferred(output::Output, ruleset::Ruleset)
    isinferred(output::Output, rules::Rule...)

Test if a custom rule is inferred and the return type is correct when
`applyrule` or `applyrule!` is run.

Type-stability can give orders of magnitude improvements in performance.
"""
isinferred(output::Output, rules::Rule...) = isinferred(output, rules)
isinferred(output::Output, rules::Tuple) = isinferred(output, Ruleset(rules...))
function isinferred(output::Output, ruleset::Ruleset)
    simdata = _updaterules(rules(ruleset), SimData(output, ruleset))
    map(rules(simdata)) do rule
        isinferred(simdata, rule)
    end
    return true
end
isinferred(simdata::AbstractSimData, rule::Rule) = _isinferred(simdata, rule)
function isinferred(simdata::AbstractSimData,
    rule::Union{NeighborhoodRule,Chain{<:Any,<:Any,<:Tuple{<:NeighborhoodRule,Vararg}}}
)
    grid = simdata[neighborhoodkey(rule)]
    r = max(1, radius(rule))
    T = eltype(grid)
    S = 2r + 1
    buffer = SArray{Tuple{S,S},T,2,S^2}(Tuple(zero(T) for i in 1:S^2))
    rule = _setbuffer(rule, buffer)
    return _isinferred(simdata, rule)
end
function isinferred(simdata::AbstractSimData, rule::SetCellRule)
    rkeys, rgrids = _getreadgrids(rule, simdata)
    wkeys, wgrids = _getwritegrids(rule, simdata)
    simdata = @set simdata.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
    readval = _readcell(simdata, rkeys, 1, 1)
    @inferred applyrule!(simdata, rule, readval, (1, 1))
    return true
end

function _isinferred(simdata, rule)
    rkeys, rgrids = _getreadgrids(rule, simdata)
    wkeys, wgrids = _getwritegrids(rule, simdata)
    simdata = @set simdata.grids = _combinegrids(rkeys, rgrids, wkeys, wgrids)
    readval = _readcell(simdata, rkeys, 1, 1)

    _example_writeval(grids::Tuple) = map(_example_writeval, grids)
    _example_writeval(grid::WritableGridData) = grid[1, 1]

    ex_writeval = Tuple(_example_writeval(wgrids))
    writeval = @inferred applyrule(simdata, rule, readval, (1, 1))
    typeof(Tuple(writeval)) == typeof(ex_writeval) ||
        error("return type `$(typeof(Tuple(writeval)))` doesn't match grids `$(typeof(ex_writeval))`")
    return true
end

# _zerogrids
# Generate a Vector of zero valued grids
_zerogrids(initgrid::AbstractArray, length) = [zero(initgrid) for f in 1:length]
_zerogrids(initgrids::NamedTuple, length) =
    [map(grid -> zero(grid), initgrids) for f in 1:length]

# _asiterable
# Return some iterable value from a
# Symbol, Tuple or tuple type
@inline _asiterable(x) = (x,)
@inline _asiterable(x::Symbol) = (x,)
@inline _asiterable(x::Type{<:Tuple}) = x.parameters
@inline _asiterable(x::Tuple) = x
@inline _asiterable(x::AbstractArray) = x

# _astuple
# Wrap a value in a tuple if the matching keys are not a tuple
# we cant just dispatch on state, as it may be meant to be a tuple.
@inline _astuple(rule::Rule, state) = _astuple(_readkeys(rule), state)
@inline _astuple(keys::Tuple, state) = state
@inline _astuple(key, state) = (state,)

# _asnamedtuple => NamedTuple
# Returns a NamedTuple given a NamedTuple or an Array.
# the Array will be called _default_.
@inline _asnamedtuple(x::NamedTuple) = x
@inline _asnamedtuple(x::AbstractArray) = (_default_=x,)
@inline function _asnamedtuple(e::Extent)
    @set! e.init = _asnamedtuple(init(e))
    @set e.padval = _samenamedtuple(init(e), padval(e))
end

# _samenamedtuple => NamedTuple
# Returns a NamedTuple with length and keys matching the `init`
# NamedTuple, for another NamedTuple, a Tuple, or a scalar.
@inline _samenamedtuple(init::NamedTuple{K}, x::NamedTuple{K}) where K = x
@noinline _samenamedtuple(init::NamedTuple{K}, x::NamedTuple{J}) where {K,J} =
    error("Keys $K and $J do not match")
@inline _samenamedtuple(init::NamedTuple{K}, x::Tuple) where K = NamedTuple{K}(x)
@inline _samenamedtuple(init::NamedTuple, x) = map(_ -> x, init)

# Unwrap a Val or Val type to its internal value
_unwrap(x) = x
_unwrap(::Val{X}) where X = X
_unwrap(::Type{<:Val{X}}) where X = X


text:
using DataFrames
using LinearAlgebra
using MixedModels
using Random
using StableRNGs
using Tables
using Test

using MixedModels: dataset

include("modelcache.jl")

@testset "simulate!" begin
    @testset "LMM" begin
        ds = dataset(:dyestuff)
        fm = only(models(:dyestuff))
        # # just in case the fit was modified in a previous test
        # refit!(fm, vec(float.(ds.yield)))
        resp₀ = copy(response(fm))
        # type conversion of ints to floats
        simulate!(StableRNG(1234321), fm, β=[1], σ=1)
        refit!(fm,resp₀)
        refit!(simulate!(StableRNG(1234321), fm))
        @test deviance(fm) ≈ 322.6582 atol=0.001
        refit!(fm, float(ds.yield))
        # Global/implicit RNG method
        Random.seed!(1234321)
        refit!(simulate!(fm))
        # just make sure this worked, don't check fit
        # (because the RNG can change between Julia versions)
        @test response(fm) ≠ resp₀
        simulate!(fm, θ = fm.θ)
        @test_throws DimensionMismatch refit!(fm, zeros(29))
        # restore the original state
        refit!(fm, vec(float.(ds.yield)))
    end

    @testset "Poisson" begin
        center(v::AbstractVector) = v .- (sum(v) / length(v))
        grouseticks = DataFrame(dataset(:grouseticks))
        grouseticks.ch = center(grouseticks.height)
        gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(), fast=true)  # fails in pirls! with fast=false
        gm4sim = refit!(simulate!(StableRNG(42), deepcopy(gm4)))
        @test isapprox(gm4.β, gm4sim.β; atol=norm(stderror(gm4)))
    end

    @testset "Binomial" begin
        cbpp = dataset(:cbpp)
        gm2 = fit(MixedModel, first(gfms[:cbpp]), cbpp, Binomial(), wts=float(cbpp.hsz))
        gm2sim = refit!(simulate!(StableRNG(42), deepcopy(gm2)), fast=true)
        @test isapprox(gm2.β, gm2sim.β; atol=norm(stderror(gm2)))
    end

    @testset "_rand with dispersion" begin
        @test_throws ArgumentError MixedModels._rand(StableRNG(42), Normal(), 1, 1, 1)
        @test_throws ArgumentError MixedModels._rand(StableRNG(42), Gamma(), 1, 1, 1)
        @test_throws ArgumentError MixedModels._rand(StableRNG(42), InverseGaussian(), 1, 1, 1)
    end
end

@testset "bootstrap" begin
    fm = only(models(:dyestuff))
    # two implicit tests
    # 1. type conversion of ints to floats
    # 2. test method for default RNG
    parametricbootstrap(1, fm, β=[1], σ=1)

    bsamp = parametricbootstrap(MersenneTwister(1234321), 100, fm, use_threads=false)
    @test isa(propertynames(bsamp), Vector{Symbol})
    @test length(bsamp.objective) == 100
    @test keys(first(bsamp.bstr)) == (:objective, :σ, :β, :se, :θ)
    @test isa(bsamp.σs, Vector{<:NamedTuple})
    @test length(bsamp.σs) == 100
    allpars = DataFrame(bsamp.allpars)
    @test isa(allpars, DataFrame)
    cov = shortestcovint(shuffle(1.:100.))
    # there is no unique shortest coverage interval here, but the left-most one
    # is currently returned, so we take that. If this behavior changes, then
    # we'll have to change the test
    @test first(cov) == 1.
    @test last(cov) == 95.

    coefp = DataFrame(bsamp.coefpvalues)

    @test isa(coefp, DataFrame)
    @test coefp.iter == 1:100
    @test only(unique(coefp.coefname)) == Symbol("(Intercept)")
    @test propertynames(coefp) == [:iter, :coefname, :β, :se, :z, :p]

    @testset "threaded bootstrap" begin
        bsamp_threaded = parametricbootstrap(MersenneTwister(1234321), 100, fm, use_threads=true)
        # even though it's bad practice with floating point, exact equality should
        # be a valid test here -- if everything is working right, then it's the exact
        # same operations occuring within each bootstrap sample, which IEEE 754
        # guarantees will yield the same result
        @test sort(bsamp_threaded.σ) == sort(bsamp.σ)
        @test sort(bsamp_threaded.θ) == sort(bsamp.θ)
        @test sort(columntable(bsamp_threaded.β).β) == sort(columntable(bsamp.β).β)
        @test sum(issingular(bsamp)) == sum(issingular(bsamp_threaded))
    end

    @testset "Bernoulli simulate! and GLMM boostrap" begin
        contra = dataset(:contra)
        gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=true)
        bs = parametricbootstrap(StableRNG(42), 100, gm0)
        bsci = combine(groupby(DataFrame(bs.β), :coefname),
                       :β => shortestcovint => :ci)
        bsci.lower = first.(bsci.ci)
        bsci.upper = last.(bsci.ci)
        select!(bsci, Not(:ci))
        ciwidth = 2 .* stderror(gm0)
        waldci = DataFrame(coef=fixefnames(gm0),
                           lower=fixef(gm0) .- ciwidth,
                           upper=fixef(gm0) .+ ciwidth)

        # coarse tolerances because we're not doing many bootstrap samples
        @test all(isapprox.(bsci.lower, waldci.lower; atol=0.5))
        @test all(isapprox.(bsci.upper, waldci.upper; atol=0.5))

        σbar = mean(MixedModels.tidyσs(bs)) do x; x.σ end
        @test σbar ≈ 0.56 atol=0.1
        apar = filter!(row -> row.type == "σ", DataFrame(MixedModels.allpars(bs)))
        @test !("Residual" in apar.names)
        @test mean(apar.value) ≈ σbar

        # can't specify dispersion for families without that parameter
        @test_throws ArgumentError parametricbootstrap(StableRNG(42), 100, gm0; σ=2)
        @test sum(issingular(bs)) == 0
    end
end


text:
using LaserTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using ContinuousPOMDPTreeSearchExperiments
using Plots
using QMDP
using JLD
using DESPOT

@everywhere begin
    using POMDPs
    using POMDPToolbox
    using ContinuousPOMDPTreeSearchExperiments
    using DiscreteValueIteration
    using ParticleFilters
    using POMCPOW
    using POMCP
    using LaserTag
    using QMDP
    using DESPOT

    N = 1000

    solvers = Dict{String, Union{Policy, Solver}}(

    #=
        "pomcpow" => begin
            # ro = MoveTowards()
            solver = POMCPOWSolver(tree_queries=500_000,
                                   criterion=MaxUCB(20.0),
                                   final_criterion=MaxTries(),
                                   max_depth=100,
                                   enable_action_pw=false,
                                   # k_action=4.0,
                                   # alpha_action=1/8,
                                   k_observation=4.0,
                                   alpha_observation=1/20,
                                   estimate_value=FOValue(ValueIterationSolver()),
                                   check_repeat_act=false,
                                   check_repeat_obs=false,
                                   init_N=InevitableInit(),
                                   init_V=InevitableInit(),
                                   rng=MersenneTwister(13)
                                  )
            solver
        end,
        =#

        # "move_towards_sampled" => MoveTowardsSampled(MersenneTwister(17)),

        "qmdp" => QMDPSolver(max_iterations=1000),

        "despot" => DESPOTSolver{LTState,
                      Int,
                      CMeas,
                      LaserBounds,
                      MersenneStreamArray}(bounds = LaserBounds(),
                                           random_streams=MersenneStreamArray(MersenneTwister(1)),
                                           rng=MersenneTwister(3),
                                           next_state=LTState([1,1], [1,1], false),
                                           curr_obs=CMeas(),
                                           time_per_move=-1.0,
                                           max_trials=500_000
                                          ),

        #=
        "discrete_pomcp" => begin
            ro = translate_policy(ToNextML(mdp(pomdp)), mdp(pomdp), dpomdp, dpomdp)
            solver = POMCPSolver(tree_queries=10_000,
                                   c=40.0,
                                   max_depth=20,
                                   estimate_value=FORollout(ro),
                                   rng=MersenneTwister(13)
                                  )
            dpolicy = solve(solver, dpomdp)
            translate_policy(dpolicy, dpomdp, pomdp, dpomdp)
        end,
        =#
    )
end

@show N

rdict = Dict{String, Any}()
for (k,sol) in solvers
    prog = Progress(N, desc="Simulating...")
    @show k
    rewards = pmap(prog, 1:N) do i
        pomdp = gen_lasertag(rng=MersenneTwister(i+600_000))
        if isa(sol,Solver)
            p = solve(deepcopy(sol), pomdp)
        else
            p = sol
        end
        hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
        up_rng = MersenneTwister(i+100_000)
        up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(100_000), 0.05, up_rng)
        hist = simulate(hr, pomdp, p, up)
        discounted_reward(hist)
    end
    @show mean(rewards)
    rdict[k] = rewards
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "laser_pomcpow_run_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
println("saving to $filename...")
@save(filename, solvers, rdict)
println("done.")