text:
export device, send_to_host, send_to_device

using Flux
using CUDA
using Adapt
using Random
using ElasticArrays

import CUDA: device

send_to_host(x) = send_to_device(Val(:cpu), x)

send_to_device(::Val{:cpu}, m) = fmap(x -> adapt(Array, x), m)
send_to_device(::Val{:gpu}, m) = fmap(CUDA.cu, m)

"""
    device(model)

Detect the suitable running device for the `model`.
Return `Val(:cpu)` by default.
"""
device(x) = device(Flux.trainable(x))
device(x::Function) = nothing
device(::CuArray) = Val(:gpu)
device(::Array) = Val(:cpu)
device(x::Tuple{}) = nothing
device(x::NamedTuple{(),Tuple{}}) = nothing
device(x::ElasticArray) = device(x.data)
device(x::SubArray) = device(parent(x))
device(x::Base.ReshapedArray) = device(parent(x))

function device(x::Random.AbstractRNG)
    if x isa CUDA.CURAND.RNG
        Val(:gpu)
    else
        Val(:cpu)
    end
end

function device(x::Union{Tuple,NamedTuple})
    d1 = device(first(x))
    if isnothing(d1)
        device(Base.tail(x))
    else
        d1
    end
end

# recoganize Torch.jl
# device(x::Tensor) = Val(Symbol(:gpu, x.device))


text:
using BinaryProvider
using CxxWrap

const tarball_url = "https://github.com/mwydmuch/ViZDoom/archive/1.1.6.tar.gz"
const tarball_hash = "c127b97dec3365acbd0ff036b4011fc46c48df49b8fa984bc47acfb9303f7db3"
const deps_dir = @__DIR__
const build_dir = joinpath(@__DIR__, "usr")
const vizdoom_dir = joinpath(build_dir, "ViZDoom-1.1.6")
const verbose = "--verbose" in ARGS

install(tarball_url, tarball_hash; prefix=Prefix(build_dir), force=true, verbose=verbose)

# replace old wrapper
run(`rm -rf $(joinpath(vizdoom_dir, "src", "lib_julia"))`)
run(`cp -r $(joinpath(deps_dir, "lib_julia")) $(joinpath(vizdoom_dir, "src"))`)

const JLCxx_DIR = get(ENV, "JLCXX_DIR", joinpath(dirname(CxxWrap.jlcxx_path), "cmake", "JlCxx"))
const cmake_cmd = `cmake -DJlCxx_DIR=$(JLCxx_DIR) -DBUILD_JULIA=ON -DJulia_EXECUTABLE=$(joinpath(Sys.BINDIR, "julia")) .`
@info cmake_cmd
cd(() -> run(cmake_cmd), vizdoom_dir)
cd(() -> run(`make`), vizdoom_dir)

ext: jl


text:
"Uniform sample of subset of preimage of `Y` under `f` unioned with `X`."
function pre_mc{T}(
    Y::SymbolicRandVar{Bool},
    niters::Integer,
    ::Type{T};
    RandVarType::Type = default_randvar(),
    args...)

  init_box = unit_box(LazyBox{Float64}, dims(Y))
  Y_conv = convert(RandVarType, Y)
  pre_mc(Y_conv, init_box, niters, T; args...)
end


text:
# VTK cell definitions

# Definitions copied from the vtkCellType.h file of the VTK source code:
# https://raw.githubusercontent.com/Kitware/VTK/master/Common/DataModel/vtkCellType.h

__precompile__()

module VTKCellTypes

export VTKCellType

struct VTKCellType
    vtk_name::String
    vtk_id::UInt8
    nodes::Int
end

const VTK_CELL_TYPES = Dict{UInt8, VTKCellType}()

# Construction of VTKCellType objects from the integer vtk id
function VTKCellType(vtkid::Integer)
    id = UInt8(vtkid)
    if !haskey(VTK_CELL_TYPES, id)
        throw(ArgumentError("unknown vtkid = $vtkid"))
    end
    VTK_CELL_TYPES[id]
end

# Define a VTK cell type
macro add(ctype, id, nodes)
    name = string(ctype)
    quote
        const $(esc(ctype)) = VTKCellType($name, $id, $nodes)
        VTK_CELL_TYPES[$id] = $(esc(ctype))
    end
end

# Linear cells
@add VTK_EMPTY_CELL 0 0
@add VTK_VERTEX 1  1
@add VTK_POLY_VERTEX 2 -1
@add VTK_LINE 3  2
@add VTK_POLY_LINE 4 -1
@add VTK_TRIANGLE 5  3
@add VTK_TRIANGLE_STRIP 6 -1
@add VTK_POLYGON 7 -1
@add VTK_PIXEL 8  4
@add VTK_QUAD 9  4
@add VTK_TETRA 10 4
@add VTK_VOXEL 11 8
@add VTK_HEXAHEDRON 12 8
@add VTK_WEDGE 13 6
@add VTK_PYRAMID 14 5
@add VTK_PENTAGONAL_PRISM 15 10
@add VTK_HEXAGONAL_PRISM 16 12

# Quadratic, isoparametric cells
@add VTK_QUADRATIC_EDGE 21 3
@add VTK_QUADRATIC_TRIANGLE 22 6
@add VTK_QUADRATIC_QUAD 23 8
@add VTK_QUADRATIC_POLYGON 36 -1
@add VTK_QUADRATIC_TETRA 24 10
@add VTK_QUADRATIC_HEXAHEDRON 25 20
@add VTK_QUADRATIC_WEDGE 26 15
@add VTK_QUADRATIC_PYRAMID 27 13
@add VTK_BIQUADRATIC_QUAD 28 9
@add VTK_TRIQUADRATIC_HEXAHEDRON 29 27
@add VTK_QUADRATIC_LINEAR_QUAD 30 6
@add VTK_QUADRATIC_LINEAR_WEDGE 31 12
@add VTK_BIQUADRATIC_QUADRATIC_WEDGE 32 18
@add VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON 33 24
@add VTK_BIQUADRATIC_TRIANGLE 34 7

# Cubic, isoparametric cell
@add VTK_CUBIC_LINE 35 4

# Special class of cells formed by convex group of points
@add VTK_CONVEX_POINT_SET 41 -1

# Polyhedron cell (consisting of polygonal faces)
@add VTK_POLYHEDRON 42 -1

# Higher order cells in parametric form
@add VTK_PARAMETRIC_CURVE 51 -1
@add VTK_PARAMETRIC_SURFACE 52 -1
@add VTK_PARAMETRIC_TRI_SURFACE 53 -1
@add VTK_PARAMETRIC_QUAD_SURFACE 54 -1
@add VTK_PARAMETRIC_TETRA_REGION 55 -1
@add VTK_PARAMETRIC_HEX_REGION 56 -1

# Higher order cells
@add VTK_HIGHER_ORDER_EDGE 60 -1
@add VTK_HIGHER_ORDER_TRIANGLE 61 -1
@add VTK_HIGHER_ORDER_QUAD 62 -1
@add VTK_HIGHER_ORDER_POLYGON 63 -1
@add VTK_HIGHER_ORDER_TETRAHEDRON 64 -1
@add VTK_HIGHER_ORDER_WEDGE 65 -1
@add VTK_HIGHER_ORDER_PYRAMID 66 -1
@add VTK_HIGHER_ORDER_HEXAHEDRON 67 -1

end # module


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const SUNDIALS_VERSION = "3.1.0"
const SUNDIALS_VERSION_MAJOR = Cint(3)
const SUNDIALS_VERSION_MINOR = Cint(1)
const SUNDIALS_VERSION_PATCH = Cint(0)
const SUNDIALS_VERSION_LABEL = ""
const SUNDIALS_DOUBLE_PRECISION = Cint(1)
const SUNDIALS_INT64_T = Cint(1)
const SUNDIALS_DENSE = Cint(1)
const SUNDIALS_BAND = Cint(2)

# Skipping MacroDefinition: DENSE_COL ( A , j ) ( ( A -> cols ) [ j ] ) #
# Skipping MacroDefinition: DENSE_ELEM ( A , i , j ) ( ( A -> cols ) [ j ] [ i ] ) /*
# * -----------------------------------------------------------------
# * BAND_COL, BAND_COL_ELEM, and BAND_ELEM
# * -----------------------------------------------------------------
# *
# * BAND_COL(A,j) references the diagonal element of the jth column
# * of the N by N band matrix A, 0 <= j <= N-1. The type of the
# * expression BAND_COL(A,j) is realtype *. The pointer returned by
# * the call BAND_COL(A,j) can be treated as an array which is
# * indexed from -(A->mu) to (A->ml).
# *
# * BAND_COL_ELEM references the (i,j)th entry of the band matrix A
# * when used in conjunction with BAND_COL. The index (i,j) should
# * satisfy j-(A->mu) <= i <= j+(A->ml).
# *
# * BAND_ELEM(A,i,j) references the (i,j)th element of the M-by-N
# * band matrix A, where 0 <= i,j <= N-1. The location (i,j) should
# * further satisfy j-(A->mu) <= i <= j+(A->ml).
# *
# * -----------------------------------------------------------------
# */
# Skipping MacroDefinition: BAND_COL ( A , j ) ( ( ( A -> cols ) [ j ] ) + ( A -> s_mu ) ) #
# Skipping MacroDefinition: BAND_COL_ELEM ( col_j , i , j ) ( col_j [ ( i ) - ( j ) ] ) #
# Skipping MacroDefinition: BAND_ELEM ( A , i , j ) ( ( A -> cols ) [ j ] [ ( i ) - ( j ) + ( A -> s_mu ) ] ) /*
# * ==================================================================
# * Exported function prototypes (functions working on dlsMat)
# * ==================================================================
# */

@static if Int == Int64
    const sunindextype = Int64
else
    const sunindextype = Int32
end

const realtype = Cdouble

mutable struct _DlsMat
    _type::Cint
    M::sunindextype
    N::sunindextype
    ldim::sunindextype
    mu::sunindextype
    ml::sunindextype
    s_mu::sunindextype
    data::Ptr{realtype}
    ldata::sunindextype
    cols::Ptr{Ptr{realtype}}
end

const DlsMat = Ptr{_DlsMat}

# begin enum ANONYMOUS_1
const ANONYMOUS_1 = UInt32
const PREC_NONE = (UInt32)(0)
const PREC_LEFT = (UInt32)(1)
const PREC_RIGHT = (UInt32)(2)
const PREC_BOTH = (UInt32)(3)
# end enum ANONYMOUS_1

# begin enum ANONYMOUS_2
const ANONYMOUS_2 = UInt32
const MODIFIED_GS = (UInt32)(1)
const CLASSICAL_GS = (UInt32)(2)
# end enum ANONYMOUS_2

const ATimesFn = Ptr{Void}
const PSetupFn = Ptr{Void}
const PSolveFn = Ptr{Void}
const SUNLS_SUCCESS = Cint(0)
const SUNLS_MEM_NULL = Cint(-1)
const SUNLS_ILL_INPUT = Cint(-2)
const SUNLS_MEM_FAIL = Cint(-3)
const SUNLS_ATIMES_FAIL_UNREC = Cint(-4)
const SUNLS_PSET_FAIL_UNREC = Cint(-5)
const SUNLS_PSOLVE_FAIL_UNREC = Cint(-6)
const SUNLS_PACKAGE_FAIL_UNREC = Cint(-7)
const SUNLS_GS_FAIL = Cint(-8)
const SUNLS_QRSOL_FAIL = Cint(-9)
const SUNLS_RES_REDUCED = Cint(1)
const SUNLS_CONV_FAIL = Cint(2)
const SUNLS_ATIMES_FAIL_REC = Cint(3)
const SUNLS_PSET_FAIL_REC = Cint(4)
const SUNLS_PSOLVE_FAIL_REC = Cint(5)
const SUNLS_PACKAGE_FAIL_REC = Cint(6)
const SUNLS_QRFACT_FAIL = Cint(7)
const SUNLS_LUFACT_FAIL = Cint(8)

# begin enum ANONYMOUS_3
const ANONYMOUS_3 = UInt32
const SUNLINEARSOLVER_DIRECT = (UInt32)(0)
const SUNLINEARSOLVER_ITERATIVE = (UInt32)(1)
const SUNLINEARSOLVER_CUSTOM = (UInt32)(2)
# end enum ANONYMOUS_3

# begin enum SUNLinearSolver_Type
const SUNLinearSolver_Type = UInt32
const SUNLINEARSOLVER_DIRECT = (UInt32)(0)
const SUNLINEARSOLVER_ITERATIVE = (UInt32)(1)
const SUNLINEARSOLVER_CUSTOM = (UInt32)(2)
# end enum SUNLinearSolver_Type

mutable struct _generic_SUNLinearSolver_Ops
    gettype::Ptr{Void}
    setatimes::Ptr{Void}
    setpreconditioner::Ptr{Void}
    setscalingvectors::Ptr{Void}
    initialize::Ptr{Void}
    setup::Ptr{Void}
    solve::Ptr{Void}
    numiters::Ptr{Void}
    resnorm::Ptr{Void}
    lastflag::Ptr{Void}
    space::Ptr{Void}
    resid::Ptr{Void}
    free::Ptr{Void}
end

const SUNLinearSolver_Ops = Ptr{_generic_SUNLinearSolver_Ops}

mutable struct _generic_SUNLinearSolver
    content::Ptr{Void}
    ops::Ptr{_generic_SUNLinearSolver_Ops}
end

const SUNLinearSolver = Ptr{_generic_SUNLinearSolver}

# begin enum ANONYMOUS_4
const ANONYMOUS_4 = UInt32
const SUNMATRIX_DENSE = (UInt32)(0)
const SUNMATRIX_BAND = (UInt32)(1)
const SUNMATRIX_SPARSE = (UInt32)(2)
const SUNMATRIX_CUSTOM = (UInt32)(3)
# end enum ANONYMOUS_4

# begin enum SUNMatrix_ID
const SUNMatrix_ID = UInt32
const SUNMATRIX_DENSE = (UInt32)(0)
const SUNMATRIX_BAND = (UInt32)(1)
const SUNMATRIX_SPARSE = (UInt32)(2)
const SUNMATRIX_CUSTOM = (UInt32)(3)
# end enum SUNMatrix_ID

mutable struct _generic_SUNMatrix_Ops
    getid::Ptr{Void}
    clone::Ptr{Void}
    destroy::Ptr{Void}
    zero::Ptr{Void}
    copy::Ptr{Void}
    scaleadd::Ptr{Void}
    scaleaddi::Ptr{Void}
    matvec::Ptr{Void}
    space::Ptr{Void}
end

const SUNMatrix_Ops = Ptr{_generic_SUNMatrix_Ops}

mutable struct _generic_SUNMatrix
    content::Ptr{Void}
    ops::Ptr{_generic_SUNMatrix_Ops}
end

const SUNMatrix = Ptr{_generic_SUNMatrix}

# begin enum ANONYMOUS_5
const ANONYMOUS_5 = UInt32
const SUNDIALS_NVEC_SERIAL = (UInt32)(0)
const SUNDIALS_NVEC_PARALLEL = (UInt32)(1)
const SUNDIALS_NVEC_OPENMP = (UInt32)(2)
const SUNDIALS_NVEC_PTHREADS = (UInt32)(3)
const SUNDIALS_NVEC_PARHYP = (UInt32)(4)
const SUNDIALS_NVEC_PETSC = (UInt32)(5)
const SUNDIALS_NVEC_CUDA = (UInt32)(6)
const SUNDIALS_NVEC_RAJA = (UInt32)(7)
const SUNDIALS_NVEC_CUSTOM = (UInt32)(8)
# end enum ANONYMOUS_5

# begin enum N_Vector_ID
const N_Vector_ID = UInt32
const SUNDIALS_NVEC_SERIAL = (UInt32)(0)
const SUNDIALS_NVEC_PARALLEL = (UInt32)(1)
const SUNDIALS_NVEC_OPENMP = (UInt32)(2)
const SUNDIALS_NVEC_PTHREADS = (UInt32)(3)
const SUNDIALS_NVEC_PARHYP = (UInt32)(4)
const SUNDIALS_NVEC_PETSC = (UInt32)(5)
const SUNDIALS_NVEC_CUDA = (UInt32)(6)
const SUNDIALS_NVEC_RAJA = (UInt32)(7)
const SUNDIALS_NVEC_CUSTOM = (UInt32)(8)
# end enum N_Vector_ID

mutable struct _generic_N_Vector_Ops
    nvgetvectorid::Ptr{Void}
    nvclone::Ptr{Void}
    nvcloneempty::Ptr{Void}
    nvdestroy::Ptr{Void}
    nvspace::Ptr{Void}
    nvgetarraypointer::Ptr{Void}
    nvsetarraypointer::Ptr{Void}
    nvlinearsum::Ptr{Void}
    nvconst::Ptr{Void}
    nvprod::Ptr{Void}
    nvdiv::Ptr{Void}
    nvscale::Ptr{Void}
    nvabs::Ptr{Void}
    nvinv::Ptr{Void}
    nvaddconst::Ptr{Void}
    nvdotprod::Ptr{Void}
    nvmaxnorm::Ptr{Void}
    nvwrmsnorm::Ptr{Void}
    nvwrmsnormmask::Ptr{Void}
    nvmin::Ptr{Void}
    nvwl2norm::Ptr{Void}
    nvl1norm::Ptr{Void}
    nvcompare::Ptr{Void}
    nvinvtest::Ptr{Void}
    nvconstrmask::Ptr{Void}
    nvminquotient::Ptr{Void}
end

const N_Vector_Ops = Ptr{_generic_N_Vector_Ops}

mutable struct _generic_N_Vector
    content::Ptr{Void}
    ops::Ptr{_generic_N_Vector_Ops}
end

const N_Vector = Ptr{_generic_N_Vector}
const N_Vector_S = Ptr{N_Vector}
const PCG_SUCCESS = Cint(0)
const PCG_RES_REDUCED = Cint(1)
const PCG_CONV_FAIL = Cint(2)
const PCG_PSOLVE_FAIL_REC = Cint(3)
const PCG_ATIMES_FAIL_REC = Cint(4)
const PCG_PSET_FAIL_REC = Cint(5)
const PCG_MEM_NULL = Cint(-1)
const PCG_ATIMES_FAIL_UNREC = Cint(-2)
const PCG_PSOLVE_FAIL_UNREC = Cint(-3)
const PCG_PSET_FAIL_UNREC = Cint(-4)

mutable struct PcgMemRec
    l_max::Cint
    r::N_Vector
    p::N_Vector
    z::N_Vector
    Ap::N_Vector
end

const PcgMem = Ptr{Void}
const SPBCG_SUCCESS = Cint(0)
const SPBCG_RES_REDUCED = Cint(1)
const SPBCG_CONV_FAIL = Cint(2)
const SPBCG_PSOLVE_FAIL_REC = Cint(3)
const SPBCG_ATIMES_FAIL_REC = Cint(4)
const SPBCG_PSET_FAIL_REC = Cint(5)
const SPBCG_MEM_NULL = Cint(-1)
const SPBCG_ATIMES_FAIL_UNREC = Cint(-2)
const SPBCG_PSOLVE_FAIL_UNREC = Cint(-3)
const SPBCG_PSET_FAIL_UNREC = Cint(-4)

# Skipping MacroDefinition: SPBCG_VTEMP ( mem ) ( mem -> r ) #

mutable struct SpbcgMemRec
    l_max::Cint
    r_star::N_Vector
    r::N_Vector
    p::N_Vector
    q::N_Vector
    u::N_Vector
    Ap::N_Vector
    vtemp::N_Vector
end

const SpbcgMem = Ptr{Void}
const SPFGMR_SUCCESS = Cint(0)
const SPFGMR_RES_REDUCED = Cint(1)
const SPFGMR_CONV_FAIL = Cint(2)
const SPFGMR_QRFACT_FAIL = Cint(3)
const SPFGMR_PSOLVE_FAIL_REC = Cint(4)
const SPFGMR_ATIMES_FAIL_REC = Cint(5)
const SPFGMR_PSET_FAIL_REC = Cint(6)
const SPFGMR_MEM_NULL = Cint(-1)
const SPFGMR_ATIMES_FAIL_UNREC = Cint(-2)
const SPFGMR_PSOLVE_FAIL_UNREC = Cint(-3)
const SPFGMR_GS_FAIL = Cint(-4)
const SPFGMR_QRSOL_FAIL = Cint(-5)
const SPFGMR_PSET_FAIL_UNREC = Cint(-6)

# Skipping MacroDefinition: SPFGMR_VTEMP ( mem ) ( mem -> vtemp ) #

mutable struct _SpfgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Z::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

mutable struct SpfgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Z::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

const SpfgmrMem = Ptr{_SpfgmrMemRec}
const SPGMR_SUCCESS = Cint(0)
const SPGMR_RES_REDUCED = Cint(1)
const SPGMR_CONV_FAIL = Cint(2)
const SPGMR_QRFACT_FAIL = Cint(3)
const SPGMR_PSOLVE_FAIL_REC = Cint(4)
const SPGMR_ATIMES_FAIL_REC = Cint(5)
const SPGMR_PSET_FAIL_REC = Cint(6)
const SPGMR_MEM_NULL = Cint(-1)
const SPGMR_ATIMES_FAIL_UNREC = Cint(-2)
const SPGMR_PSOLVE_FAIL_UNREC = Cint(-3)
const SPGMR_GS_FAIL = Cint(-4)
const SPGMR_QRSOL_FAIL = Cint(-5)
const SPGMR_PSET_FAIL_UNREC = Cint(-6)

# Skipping MacroDefinition: SPGMR_VTEMP ( mem ) ( mem -> vtemp ) #

mutable struct _SpgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

mutable struct SpgmrMemRec
    l_max::Cint
    V::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

const SpgmrMem = Ptr{_SpgmrMemRec}
const SPTFQMR_SUCCESS = Cint(0)
const SPTFQMR_RES_REDUCED = Cint(1)
const SPTFQMR_CONV_FAIL = Cint(2)
const SPTFQMR_PSOLVE_FAIL_REC = Cint(3)
const SPTFQMR_ATIMES_FAIL_REC = Cint(4)
const SPTFQMR_PSET_FAIL_REC = Cint(5)
const SPTFQMR_MEM_NULL = Cint(-1)
const SPTFQMR_ATIMES_FAIL_UNREC = Cint(-2)
const SPTFQMR_PSOLVE_FAIL_UNREC = Cint(-3)
const SPTFQMR_PSET_FAIL_UNREC = Cint(-4)

# Skipping MacroDefinition: SPTFQMR_VTEMP ( mem ) ( mem -> vtemp1 ) #

mutable struct SptfqmrMemRec
    l_max::Cint
    r_star::N_Vector
    q::N_Vector
    d::N_Vector
    v::N_Vector
    p::N_Vector
    r::Ptr{N_Vector}
    u::N_Vector
    vtemp1::N_Vector
    vtemp2::N_Vector
    vtemp3::N_Vector
end

const SptfqmrMem = Ptr{Void}

# Skipping MacroDefinition: RCONST ( x ) x #

const BIG_REAL = DBL_MAX
const SMALL_REAL = DBL_MIN
const UNIT_ROUNDOFF = DBL_EPSILON

# Skipping MacroDefinition: booleantype int #

const SUNFALSE = Cint(0)
const SUNTRUE = Cint(1)
const CV_ADAMS = Cint(1)
const CV_BDF = Cint(2)
const CV_FUNCTIONAL = Cint(1)
const CV_NEWTON = Cint(2)
const CV_NORMAL = Cint(1)
const CV_ONE_STEP = Cint(2)
const CV_SUCCESS = Cint(0)
const CV_TSTOP_RETURN = Cint(1)
const CV_ROOT_RETURN = Cint(2)
const CV_WARNING = Cint(99)
const CV_TOO_MUCH_WORK = Cint(-1)
const CV_TOO_MUCH_ACC = Cint(-2)
const CV_ERR_FAILURE = Cint(-3)
const CV_CONV_FAILURE = Cint(-4)
const CV_LINIT_FAIL = Cint(-5)
const CV_LSETUP_FAIL = Cint(-6)
const CV_LSOLVE_FAIL = Cint(-7)
const CV_RHSFUNC_FAIL = Cint(-8)
const CV_FIRST_RHSFUNC_ERR = Cint(-9)
const CV_REPTD_RHSFUNC_ERR = Cint(-10)
const CV_UNREC_RHSFUNC_ERR = Cint(-11)
const CV_RTFUNC_FAIL = Cint(-12)
const CV_MEM_FAIL = Cint(-20)
const CV_MEM_NULL = Cint(-21)
const CV_ILL_INPUT = Cint(-22)
const CV_NO_MALLOC = Cint(-23)
const CV_BAD_K = Cint(-24)
const CV_BAD_T = Cint(-25)
const CV_BAD_DKY = Cint(-26)
const CV_TOO_CLOSE = Cint(-27)
const CVRhsFn = Ptr{Void}
const CVRootFn = Ptr{Void}
const CVEwtFn = Ptr{Void}
const CVErrHandlerFn = Ptr{Void}
const CV_SIMULTANEOUS = Cint(1)
const CV_STAGGERED = Cint(2)
const CV_STAGGERED1 = Cint(3)
const CV_CENTERED = Cint(1)
const CV_FORWARD = Cint(2)
const CV_HERMITE = Cint(1)
const CV_POLYNOMIAL = Cint(2)
const CV_NO_QUAD = Cint(-30)
const CV_QRHSFUNC_FAIL = Cint(-31)
const CV_FIRST_QRHSFUNC_ERR = Cint(-32)
const CV_REPTD_QRHSFUNC_ERR = Cint(-33)
const CV_UNREC_QRHSFUNC_ERR = Cint(-34)
const CV_NO_SENS = Cint(-40)
const CV_SRHSFUNC_FAIL = Cint(-41)
const CV_FIRST_SRHSFUNC_ERR = Cint(-42)
const CV_REPTD_SRHSFUNC_ERR = Cint(-43)
const CV_UNREC_SRHSFUNC_ERR = Cint(-44)
const CV_BAD_IS = Cint(-45)
const CV_NO_QUADSENS = Cint(-50)
const CV_QSRHSFUNC_FAIL = Cint(-51)
const CV_FIRST_QSRHSFUNC_ERR = Cint(-52)
const CV_REPTD_QSRHSFUNC_ERR = Cint(-53)
const CV_UNREC_QSRHSFUNC_ERR = Cint(-54)
const CV_NO_ADJ = Cint(-101)
const CV_NO_FWD = Cint(-102)
const CV_NO_BCK = Cint(-103)
const CV_BAD_TB0 = Cint(-104)
const CV_REIFWD_FAIL = Cint(-105)
const CV_FWD_FAIL = Cint(-106)
const CV_GETY_BADT = Cint(-107)
const CVQuadRhsFn = Ptr{Void}
const CVSensRhsFn = Ptr{Void}
const CVSensRhs1Fn = Ptr{Void}
const CVQuadSensRhsFn = Ptr{Void}
const CVRhsFnB = Ptr{Void}
const CVRhsFnBS = Ptr{Void}
const CVQuadRhsFnB = Ptr{Void}
const CVQuadRhsFnBS = Ptr{Void}

mutable struct CVadjCheckPointRec
    my_addr::Ptr{Void}
    next_addr::Ptr{Void}
    t0::realtype
    t1::realtype
    nstep::Clong
    order::Cint
    step::realtype
end

const IDA_NORMAL = Cint(1)
const IDA_ONE_STEP = Cint(2)
const IDA_YA_YDP_INIT = Cint(1)
const IDA_Y_INIT = Cint(2)
const IDA_SUCCESS = Cint(0)
const IDA_TSTOP_RETURN = Cint(1)
const IDA_ROOT_RETURN = Cint(2)
const IDA_WARNING = Cint(99)
const IDA_TOO_MUCH_WORK = Cint(-1)
const IDA_TOO_MUCH_ACC = Cint(-2)
const IDA_ERR_FAIL = Cint(-3)
const IDA_CONV_FAIL = Cint(-4)
const IDA_LINIT_FAIL = Cint(-5)
const IDA_LSETUP_FAIL = Cint(-6)
const IDA_LSOLVE_FAIL = Cint(-7)
const IDA_RES_FAIL = Cint(-8)
const IDA_REP_RES_ERR = Cint(-9)
const IDA_RTFUNC_FAIL = Cint(-10)
const IDA_CONSTR_FAIL = Cint(-11)
const IDA_FIRST_RES_FAIL = Cint(-12)
const IDA_LINESEARCH_FAIL = Cint(-13)
const IDA_NO_RECOVERY = Cint(-14)
const IDA_MEM_NULL = Cint(-20)
const IDA_MEM_FAIL = Cint(-21)
const IDA_ILL_INPUT = Cint(-22)
const IDA_NO_MALLOC = Cint(-23)
const IDA_BAD_EWT = Cint(-24)
const IDA_BAD_K = Cint(-25)
const IDA_BAD_T = Cint(-26)
const IDA_BAD_DKY = Cint(-27)
const IDA_UNRECOGNISED_ERROR = Cint(-99)
const IDAResFn = Ptr{Void}
const IDARootFn = Ptr{Void}
const IDAEwtFn = Ptr{Void}
const IDAErrHandlerFn = Ptr{Void}
const IDA_SIMULTANEOUS = Cint(1)
const IDA_STAGGERED = Cint(2)
const IDA_CENTERED = Cint(1)
const IDA_FORWARD = Cint(2)
const IDA_HERMITE = Cint(1)
const IDA_POLYNOMIAL = Cint(2)
const IDA_NO_QUAD = Cint(-30)
const IDA_QRHS_FAIL = Cint(-31)
const IDA_FIRST_QRHS_ERR = Cint(-32)
const IDA_REP_QRHS_ERR = Cint(-33)
const IDA_NO_SENS = Cint(-40)
const IDA_SRES_FAIL = Cint(-41)
const IDA_REP_SRES_ERR = Cint(-42)
const IDA_BAD_IS = Cint(-43)
const IDA_NO_QUADSENS = Cint(-50)
const IDA_QSRHS_FAIL = Cint(-51)
const IDA_FIRST_QSRHS_ERR = Cint(-52)
const IDA_REP_QSRHS_ERR = Cint(-53)
const IDA_UNRECONGISED_ERROR = Cint(-99)
const IDA_NO_ADJ = Cint(-101)
const IDA_NO_FWD = Cint(-102)
const IDA_NO_BCK = Cint(-103)
const IDA_BAD_TB0 = Cint(-104)
const IDA_REIFWD_FAIL = Cint(-105)
const IDA_FWD_FAIL = Cint(-106)
const IDA_GETY_BADT = Cint(-107)
const IDAQuadRhsFn = Ptr{Void}
const IDASensResFn = Ptr{Void}
const IDAQuadSensRhsFn = Ptr{Void}
const IDAResFnB = Ptr{Void}
const IDAResFnBS = Ptr{Void}
const IDAQuadRhsFnB = Ptr{Void}
const IDAQuadRhsFnBS = Ptr{Void}

mutable struct IDAadjCheckPointRec
    my_addr::Ptr{Void}
    next_addr::Ptr{Void}
    t0::realtype
    t1::realtype
    nstep::Clong
    order::Cint
    step::realtype
end

const KIN_SUCCESS = Cint(0)
const KIN_INITIAL_GUESS_OK = Cint(1)
const KIN_STEP_LT_STPTOL = Cint(2)
const KIN_WARNING = Cint(99)
const KIN_MEM_NULL = Cint(-1)
const KIN_ILL_INPUT = Cint(-2)
const KIN_NO_MALLOC = Cint(-3)
const KIN_MEM_FAIL = Cint(-4)
const KIN_LINESEARCH_NONCONV = Cint(-5)
const KIN_MAXITER_REACHED = Cint(-6)
const KIN_MXNEWT_5X_EXCEEDED = Cint(-7)
const KIN_LINESEARCH_BCFAIL = Cint(-8)
const KIN_LINSOLV_NO_RECOVERY = Cint(-9)
const KIN_LINIT_FAIL = Cint(-10)
const KIN_LSETUP_FAIL = Cint(-11)
const KIN_LSOLVE_FAIL = Cint(-12)
const KIN_SYSFUNC_FAIL = Cint(-13)
const KIN_FIRST_SYSFUNC_ERR = Cint(-14)
const KIN_REPTD_SYSFUNC_ERR = Cint(-15)
const KIN_ETACHOICE1 = Cint(1)
const KIN_ETACHOICE2 = Cint(2)
const KIN_ETACONSTANT = Cint(3)
const KIN_NONE = Cint(0)
const KIN_LINESEARCH = Cint(1)
const KIN_PICARD = Cint(2)
const KIN_FP = Cint(3)
const KINSysFn = Ptr{Void}
const KINErrHandlerFn = Ptr{Void}
const KINInfoHandlerFn = Ptr{Void}
const ARK_S_MAX = Cint(15)

# Skipping MacroDefinition: ARK_A ( A , i , j ) ( A [ i * ARK_S_MAX + j ] ) /* itask */

const ARK_NORMAL = Cint(1)
const ARK_ONE_STEP = Cint(2)
const HEUN_EULER_2_1_2 = Cint(0)
const BOGACKI_SHAMPINE_4_2_3 = Cint(1)
const ARK324L2SA_ERK_4_2_3 = Cint(2)
const ZONNEVELD_5_3_4 = Cint(3)
const ARK436L2SA_ERK_6_3_4 = Cint(4)
const SAYFY_ABURUB_6_3_4 = Cint(5)
const CASH_KARP_6_4_5 = Cint(6)
const FEHLBERG_6_4_5 = Cint(7)
const DORMAND_PRINCE_7_4_5 = Cint(8)
const ARK548L2SA_ERK_8_4_5 = Cint(9)
const VERNER_8_5_6 = Cint(10)
const FEHLBERG_13_7_8 = Cint(11)
const DEFAULT_ERK_2 = HEUN_EULER_2_1_2
const DEFAULT_ERK_3 = BOGACKI_SHAMPINE_4_2_3
const DEFAULT_ERK_4 = ZONNEVELD_5_3_4
const DEFAULT_ERK_5 = CASH_KARP_6_4_5
const DEFAULT_ERK_6 = VERNER_8_5_6
const DEFAULT_ERK_8 = FEHLBERG_13_7_8
const MIN_ERK_NUM = Cint(0)
const MAX_ERK_NUM = Cint(11)
const SDIRK_2_1_2 = Cint(12)
const BILLINGTON_3_3_2 = Cint(13)
const TRBDF2_3_3_2 = Cint(14)
const KVAERNO_4_2_3 = Cint(15)
const ARK324L2SA_DIRK_4_2_3 = Cint(16)
const CASH_5_2_4 = Cint(17)
const CASH_5_3_4 = Cint(18)
const SDIRK_5_3_4 = Cint(19)
const KVAERNO_5_3_4 = Cint(20)
const ARK436L2SA_DIRK_6_3_4 = Cint(21)
const KVAERNO_7_4_5 = Cint(22)
const ARK548L2SA_DIRK_8_4_5 = Cint(23)
const DEFAULT_DIRK_2 = SDIRK_2_1_2
const DEFAULT_DIRK_3 = ARK324L2SA_DIRK_4_2_3
const DEFAULT_DIRK_4 = SDIRK_5_3_4
const DEFAULT_DIRK_5 = ARK548L2SA_DIRK_8_4_5
const MIN_DIRK_NUM = Cint(12)
const MAX_DIRK_NUM = Cint(23)
const DEFAULT_ARK_ETABLE_3 = ARK324L2SA_ERK_4_2_3
const DEFAULT_ARK_ETABLE_4 = ARK436L2SA_ERK_6_3_4
const DEFAULT_ARK_ETABLE_5 = ARK548L2SA_ERK_8_4_5
const DEFAULT_ARK_ITABLE_3 = ARK324L2SA_DIRK_4_2_3
const DEFAULT_ARK_ITABLE_4 = ARK436L2SA_DIRK_6_3_4
const DEFAULT_ARK_ITABLE_5 = ARK548L2SA_DIRK_8_4_5
const ARK_SUCCESS = Cint(0)
const ARK_TSTOP_RETURN = Cint(1)
const ARK_ROOT_RETURN = Cint(2)
const ARK_WARNING = Cint(99)
const ARK_TOO_MUCH_WORK = Cint(-1)
const ARK_TOO_MUCH_ACC = Cint(-2)
const ARK_ERR_FAILURE = Cint(-3)
const ARK_CONV_FAILURE = Cint(-4)
const ARK_LINIT_FAIL = Cint(-5)
const ARK_LSETUP_FAIL = Cint(-6)
const ARK_LSOLVE_FAIL = Cint(-7)
const ARK_RHSFUNC_FAIL = Cint(-8)
const ARK_FIRST_RHSFUNC_ERR = Cint(-9)
const ARK_REPTD_RHSFUNC_ERR = Cint(-10)
const ARK_UNREC_RHSFUNC_ERR = Cint(-11)
const ARK_RTFUNC_FAIL = Cint(-12)
const ARK_LFREE_FAIL = Cint(-13)
const ARK_MASSINIT_FAIL = Cint(-14)
const ARK_MASSSETUP_FAIL = Cint(-15)
const ARK_MASSSOLVE_FAIL = Cint(-16)
const ARK_MASSFREE_FAIL = Cint(-17)
const ARK_MASSMULT_FAIL = Cint(-18)
const ARK_MEM_FAIL = Cint(-20)
const ARK_MEM_NULL = Cint(-21)
const ARK_ILL_INPUT = Cint(-22)
const ARK_NO_MALLOC = Cint(-23)
const ARK_BAD_K = Cint(-24)
const ARK_BAD_T = Cint(-25)
const ARK_BAD_DKY = Cint(-26)
const ARK_TOO_CLOSE = Cint(-27)
const ARK_POSTPROCESS_FAIL = Cint(-28)
const ARKRhsFn = Ptr{Void}
const ARKRootFn = Ptr{Void}
const ARKEwtFn = Ptr{Void}
const ARKRwtFn = Ptr{Void}
const ARKErrHandlerFn = Ptr{Void}
const ARKAdaptFn = Ptr{Void}
const ARKExpStabFn = Ptr{Void}
const ARKVecResizeFn = Ptr{Void}
const ARKPostProcessStepFn = Ptr{Void}
const ARKSPILS_SUCCESS = Cint(0)
const ARKSPILS_MEM_NULL = Cint(-1)
const ARKSPILS_LMEM_NULL = Cint(-2)
const ARKSPILS_ILL_INPUT = Cint(-3)
const ARKSPILS_MEM_FAIL = Cint(-4)
const ARKSPILS_PMEM_NULL = Cint(-5)
const ARKSPILS_MASSMEM_NULL = Cint(-6)
const ARKSPILS_SUNLS_FAIL = Cint(-10)
const ARKSPILS_MSBPRE = Cint(50)

# Skipping MacroDefinition: ARKSPILS_DGMAX RCONST ( 0.2 ) /* -- TURN INTO A PARAMETER, MAKE THIS THE DEFAULT */
# Skipping MacroDefinition: ARKSPILS_EPLIN RCONST ( 0.05 ) /*===============================================================#
# ARKSPILS user-supplied function prototypes#
# ===============================================================*/

const ARKSpilsPrecSetupFn = Ptr{Void}
const ARKSpilsPrecSolveFn = Ptr{Void}
const ARKSpilsJacTimesSetupFn = Ptr{Void}
const ARKSpilsJacTimesVecFn = Ptr{Void}
const ARKSpilsMassTimesSetupFn = Ptr{Void}
const ARKSpilsMassTimesVecFn = Ptr{Void}
const ARKSpilsMassPrecSetupFn = Ptr{Void}
const ARKSpilsMassPrecSolveFn = Ptr{Void}

# Skipping MacroDefinition: SM_CONTENT_B ( A ) ( ( SUNMatrixContent_Band ) ( A -> content ) ) #
# Skipping MacroDefinition: SM_ROWS_B ( A ) ( SM_CONTENT_B ( A ) -> M ) #
# Skipping MacroDefinition: SM_COLUMNS_B ( A ) ( SM_CONTENT_B ( A ) -> N ) #
# Skipping MacroDefinition: SM_LDATA_B ( A ) ( SM_CONTENT_B ( A ) -> ldata ) #
# Skipping MacroDefinition: SM_UBAND_B ( A ) ( SM_CONTENT_B ( A ) -> mu ) #
# Skipping MacroDefinition: SM_LBAND_B ( A ) ( SM_CONTENT_B ( A ) -> ml ) #
# Skipping MacroDefinition: SM_SUBAND_B ( A ) ( SM_CONTENT_B ( A ) -> s_mu ) #
# Skipping MacroDefinition: SM_LDIM_B ( A ) ( SM_CONTENT_B ( A ) -> ldim ) #
# Skipping MacroDefinition: SM_DATA_B ( A ) ( SM_CONTENT_B ( A ) -> data ) #
# Skipping MacroDefinition: SM_COLS_B ( A ) ( SM_CONTENT_B ( A ) -> cols ) #
# Skipping MacroDefinition: SM_COLUMN_B ( A , j ) ( ( ( SM_CONTENT_B ( A ) -> cols ) [ j ] ) + SM_SUBAND_B ( A ) ) #
# Skipping MacroDefinition: SM_COLUMN_ELEMENT_B ( col_j , i , j ) ( col_j [ ( i ) - ( j ) ] ) #
# Skipping MacroDefinition: SM_ELEMENT_B ( A , i , j ) ( ( SM_CONTENT_B ( A ) -> cols ) [ j ] [ ( i ) - ( j ) + SM_SUBAND_B ( A ) ] ) /*
# * -----------------------------------------------------------------
# * PART III: functions exported by sunmatrix_band
# *
# * CONSTRUCTORS:
# *    SUNBandMatrix
# * OTHER:
# *    SUNBandMatrix_Print
# *    SUNBandMatrix_Rows
# *    SUNBandMatrix_Columns
# *    SUNBandMatrix_LowerBandwidth
# *    SUNBandMatrix_UpperBandwidth
# *    SUNBandMatrix_StoredUpperBandwidth
# *    SUNBandMatrix_LDim
# *    SUNBandMatrix_Data
# *    SUNBandMatrix_Cols
# *    SUNBandMatrix_Column
# * -----------------------------------------------------------------
# */

mutable struct _SUNMatrixContent_Band
    M::sunindextype
    N::sunindextype
    ldim::sunindextype
    mu::sunindextype
    ml::sunindextype
    s_mu::sunindextype
    data::Ptr{realtype}
    ldata::sunindextype
    cols::Ptr{Ptr{realtype}}
end

const SUNMatrixContent_Band = Ptr{_SUNMatrixContent_Band}

# Skipping MacroDefinition: SM_CONTENT_D ( A ) ( ( SUNMatrixContent_Dense ) ( A -> content ) ) #
# Skipping MacroDefinition: SM_ROWS_D ( A ) ( SM_CONTENT_D ( A ) -> M ) #
# Skipping MacroDefinition: SM_COLUMNS_D ( A ) ( SM_CONTENT_D ( A ) -> N ) #
# Skipping MacroDefinition: SM_LDATA_D ( A ) ( SM_CONTENT_D ( A ) -> ldata ) #
# Skipping MacroDefinition: SM_DATA_D ( A ) ( SM_CONTENT_D ( A ) -> data ) #
# Skipping MacroDefinition: SM_COLS_D ( A ) ( SM_CONTENT_D ( A ) -> cols ) #
# Skipping MacroDefinition: SM_COLUMN_D ( A , j ) ( ( SM_CONTENT_D ( A ) -> cols ) [ j ] ) #
# Skipping MacroDefinition: SM_ELEMENT_D ( A , i , j ) ( ( SM_CONTENT_D ( A ) -> cols ) [ j ] [ i ] ) /*
# * -----------------------------------------------------------------
# * PART III: functions exported by sunmatrix_dense
# *
# * CONSTRUCTORS:
# *    SUNDenseMatrix
# * OTHER:
# *    SUNDenseMatrix_Print
# *    SUNDenseMatrix_Rows
# *    SUNDenseMatrix_Columns
# *    SUNDenseMatrix_LData
# *    SUNDenseMatrix_Data
# *    SUNDenseMatrix_Cols
# *    SUNDenseMatrix_Column
# * -----------------------------------------------------------------
# */

mutable struct _SUNMatrixContent_Dense
    M::sunindextype
    N::sunindextype
    data::Ptr{realtype}
    ldata::sunindextype
    cols::Ptr{Ptr{realtype}}
end

const SUNMatrixContent_Dense = Ptr{_SUNMatrixContent_Dense}

# Skipping MacroDefinition: NV_CONTENT_S ( v ) ( ( N_VectorContent_Serial ) ( v -> content ) ) #
# Skipping MacroDefinition: NV_LENGTH_S ( v ) ( NV_CONTENT_S ( v ) -> length ) #
# Skipping MacroDefinition: NV_OWN_DATA_S ( v ) ( NV_CONTENT_S ( v ) -> own_data ) #
# Skipping MacroDefinition: NV_DATA_S ( v ) ( NV_CONTENT_S ( v ) -> data ) #
# Skipping MacroDefinition: NV_Ith_S ( v , i ) ( NV_DATA_S ( v ) [ i ] ) /*
# * -----------------------------------------------------------------
# * PART III: functions exported by nvector_serial
# *
# * CONSTRUCTORS:
# *    N_VNew_Serial
# *    N_VNewEmpty_Serial
# *    N_VMake_Serial
# *    N_VCloneVectorArray_Serial
# *    N_VCloneVectorArrayEmpty_Serial
# * DESTRUCTORS:
# *    N_VDestroy_Serial
# *    N_VDestroyVectorArray_Serial
# * OTHER:
# *    N_VGetLength_Serial
# *    N_VPrint_Serial
# *    N_VPrintFile_Serial
# * -----------------------------------------------------------------
# */

mutable struct _N_VectorContent_Serial
    length::sunindextype
    own_data::Cint
    data::Ptr{realtype}
end

const N_VectorContent_Serial = Ptr{_N_VectorContent_Serial}
const FCMIX_CVODE = Cint(1)
const FCMIX_IDA = Cint(2)
const FCMIX_KINSOL = Cint(3)
const FCMIX_ARKODE = Cint(4)

# Skipping MacroDefinition: SUNMIN ( A , B ) ( ( A ) < ( B ) ? ( A ) : ( B ) ) #
# Skipping MacroDefinition: SUNMAX ( A , B ) ( ( A ) > ( B ) ? ( A ) : ( B ) ) #
# Skipping MacroDefinition: SUNSQR ( A ) ( ( A ) * ( A ) ) #

const CSC_MAT = Cint(0)
const CSR_MAT = Cint(1)

mutable struct _SlsMat
    M::Cint
    N::Cint
    NNZ::Cint
    NP::Cint
    data::Ptr{realtype}
    sparsetype::Cint
    indexvals::Ptr{Cint}
    indexptrs::Ptr{Cint}
    rowvals::Ptr{Ptr{Cint}}
    colptrs::Ptr{Ptr{Cint}}
    colvals::Ptr{Ptr{Cint}}
    rowptrs::Ptr{Ptr{Cint}}
end

const SlsMat = Ptr{_SlsMat}
const CVLocalFn = Ptr{Void}
const CVCommFn = Ptr{Void}
const CVDIAG_SUCCESS = Cint(0)
const CVDIAG_MEM_NULL = Cint(-1)
const CVDIAG_LMEM_NULL = Cint(-2)
const CVDIAG_ILL_INPUT = Cint(-3)
const CVDIAG_MEM_FAIL = Cint(-4)
const CVDIAG_INV_FAIL = Cint(-5)
const CVDIAG_RHSFUNC_UNRECVR = Cint(-6)
const CVDIAG_RHSFUNC_RECVR = Cint(-7)
const CVDLS_SUCCESS = Cint(0)
const CVDLS_MEM_NULL = Cint(-1)
const CVDLS_LMEM_NULL = Cint(-2)
const CVDLS_ILL_INPUT = Cint(-3)
const CVDLS_MEM_FAIL = Cint(-4)
const CVDLS_JACFUNC_UNRECVR = Cint(-5)
const CVDLS_JACFUNC_RECVR = Cint(-6)
const CVDLS_SUNMAT_FAIL = Cint(-7)
const CVDlsJacFn = Ptr{Void}
const CVSPILS_SUCCESS = Cint(0)
const CVSPILS_MEM_NULL = Cint(-1)
const CVSPILS_LMEM_NULL = Cint(-2)
const CVSPILS_ILL_INPUT = Cint(-3)
const CVSPILS_MEM_FAIL = Cint(-4)
const CVSPILS_PMEM_NULL = Cint(-5)
const CVSPILS_SUNLS_FAIL = Cint(-6)
const CVSPILS_MSBPRE = Cint(50)

# Skipping MacroDefinition: CVSPILS_DGMAX RCONST ( 0.2 ) #
# Skipping MacroDefinition: CVSPILS_EPLIN RCONST ( 0.05 ) /*===============================================================#
#  CVSPILS user-supplied function prototypes#
#  ===============================================================*/

const CVSpilsPrecSetupFn = Ptr{Void}
const CVSpilsPrecSolveFn = Ptr{Void}
const CVSpilsJacTimesSetupFn = Ptr{Void}
const CVSpilsJacTimesVecFn = Ptr{Void}
const CVLocalFnB = Ptr{Void}
const CVCommFnB = Ptr{Void}
const CVDIAG_NO_ADJ = Cint(-101)
const CVDLS_NO_ADJ = Cint(-101)
const CVDLS_LMEMB_NULL = Cint(-102)
const CVDlsJacFnB = Ptr{Void}
const CVDlsJacFnBS = Ptr{Void}
const CVSPILS_NO_ADJ = Cint(-101)
const CVSPILS_LMEMB_NULL = Cint(-102)

# Skipping MacroDefinition: CVSPILS_EPLIN RCONST ( 0.05 ) /*-----------------------------------------------------------------#
#  PART I - forward problems#
#  -----------------------------------------------------------------*/

const CVSpilsPrecSetupFnB = Ptr{Void}
const CVSpilsPrecSetupFnBS = Ptr{Void}
const CVSpilsPrecSolveFnB = Ptr{Void}
const CVSpilsPrecSolveFnBS = Ptr{Void}
const CVSpilsJacTimesSetupFnB = Ptr{Void}
const CVSpilsJacTimesSetupFnBS = Ptr{Void}
const CVSpilsJacTimesVecFnB = Ptr{Void}
const CVSpilsJacTimesVecFnBS = Ptr{Void}
const IDABBDLocalFn = Ptr{Void}
const IDABBDCommFn = Ptr{Void}
const IDADLS_SUCCESS = Cint(0)
const IDADLS_MEM_NULL = Cint(-1)
const IDADLS_LMEM_NULL = Cint(-2)
const IDADLS_ILL_INPUT = Cint(-3)
const IDADLS_MEM_FAIL = Cint(-4)
const IDADLS_JACFUNC_UNRECVR = Cint(-5)
const IDADLS_JACFUNC_RECVR = Cint(-6)
const IDADLS_SUNMAT_FAIL = Cint(-7)
const IDADlsJacFn = Ptr{Void}
const IDASPILS_SUCCESS = Cint(0)
const IDASPILS_MEM_NULL = Cint(-1)
const IDASPILS_LMEM_NULL = Cint(-2)
const IDASPILS_ILL_INPUT = Cint(-3)
const IDASPILS_MEM_FAIL = Cint(-4)
const IDASPILS_PMEM_NULL = Cint(-5)
const IDASPILS_SUNLS_FAIL = Cint(-6)
const IDASpilsPrecSetupFn = Ptr{Void}
const IDASpilsPrecSolveFn = Ptr{Void}
const IDASpilsJacTimesSetupFn = Ptr{Void}
const IDASpilsJacTimesVecFn = Ptr{Void}
const IDABBDLocalFnB = Ptr{Void}
const IDABBDCommFnB = Ptr{Void}
const IDADLS_NO_ADJ = Cint(-101)
const IDADLS_LMEMB_NULL = Cint(-102)
const IDADlsJacFnB = Ptr{Void}
const IDADlsJacFnBS = Ptr{Void}
const IDASPILS_NO_ADJ = Cint(-101)
const IDASPILS_LMEMB_NULL = Cint(-102)
const IDASpilsPrecSetupFnB = Ptr{Void}
const IDASpilsPrecSetupFnBS = Ptr{Void}
const IDASpilsPrecSolveFnB = Ptr{Void}
const IDASpilsPrecSolveFnBS = Ptr{Void}
const IDASpilsJacTimesSetupFnB = Ptr{Void}
const IDASpilsJacTimesSetupFnBS = Ptr{Void}
const IDASpilsJacTimesVecFnB = Ptr{Void}
const IDASpilsJacTimesVecFnBS = Ptr{Void}
const KINBBDPRE_SUCCESS = Cint(0)
const KINBBDPRE_PDATA_NULL = Cint(-11)
const KINBBDPRE_FUNC_UNRECVR = Cint(-12)
const KINBBDCommFn = Ptr{Void}
const KINBBDLocalFn = Ptr{Void}
const KINDLS_SUCCESS = Cint(0)
const KINDLS_MEM_NULL = Cint(-1)
const KINDLS_LMEM_NULL = Cint(-2)
const KINDLS_ILL_INPUT = Cint(-3)
const KINDLS_MEM_FAIL = Cint(-4)
const KINDLS_JACFUNC_ERR = Cint(-5)
const KINDLS_SUNMAT_FAIL = Cint(-6)
const KINDlsJacFn = Ptr{Void}
const KINSPILS_SUCCESS = Cint(0)
const KINSPILS_MEM_NULL = Cint(-1)
const KINSPILS_LMEM_NULL = Cint(-2)
const KINSPILS_ILL_INPUT = Cint(-3)
const KINSPILS_MEM_FAIL = Cint(-4)
const KINSPILS_PMEM_NULL = Cint(-5)
const KINSPILS_SUNLS_FAIL = Cint(-6)
const KINSpilsPrecSetupFn = Ptr{Void}
const KINSpilsPrecSolveFn = Ptr{Void}
const KINSpilsJacTimesVecFn = Ptr{Void}
const ARKLocalFn = Ptr{Void}
const ARKCommFn = Ptr{Void}
const ARKDLS_SUCCESS = Cint(0)
const ARKDLS_MEM_NULL = Cint(-1)
const ARKDLS_LMEM_NULL = Cint(-2)
const ARKDLS_ILL_INPUT = Cint(-3)
const ARKDLS_MEM_FAIL = Cint(-4)
const ARKDLS_MASSMEM_NULL = Cint(-5)
const ARKDLS_JACFUNC_UNRECVR = Cint(-6)
const ARKDLS_JACFUNC_RECVR = Cint(-7)
const ARKDLS_MASSFUNC_UNRECVR = Cint(-8)
const ARKDLS_MASSFUNC_RECVR = Cint(-9)
const ARKDLS_SUNMAT_FAIL = Cint(-10)
const ARKDlsJacFn = Ptr{Void}
const ARKDlsMassFn = Ptr{Void}

mutable struct _SUNLinearSolverContent_Band
    N::sunindextype
    pivots::Ptr{sunindextype}
    last_flag::Clong
end

const SUNLinearSolverContent_Band = Ptr{_SUNLinearSolverContent_Band}

mutable struct _SUNLinearSolverContent_Dense
    N::sunindextype
    pivots::Ptr{sunindextype}
    last_flag::Clong
end

const SUNLinearSolverContent_Dense = Ptr{_SUNLinearSolverContent_Dense}
const SUNPCG_MAXL_DEFAULT = Cint(5)

mutable struct _SUNLinearSolverContent_PCG
    maxl::Cint
    pretype::Cint
    numiters::Cint
    resnorm::realtype
    last_flag::Clong
    ATimes::ATimesFn
    ATData::Ptr{Void}
    Psetup::PSetupFn
    Psolve::PSolveFn
    PData::Ptr{Void}
    s::N_Vector
    r::N_Vector
    p::N_Vector
    z::N_Vector
    Ap::N_Vector
end

const SUNLinearSolverContent_PCG = Ptr{_SUNLinearSolverContent_PCG}
const SUNSPBCGS_MAXL_DEFAULT = Cint(5)

mutable struct _SUNLinearSolverContent_SPBCGS
    maxl::Cint
    pretype::Cint
    numiters::Cint
    resnorm::realtype
    last_flag::Clong
    ATimes::ATimesFn
    ATData::Ptr{Void}
    Psetup::PSetupFn
    Psolve::PSolveFn
    PData::Ptr{Void}
    s1::N_Vector
    s2::N_Vector
    r::N_Vector
    r_star::N_Vector
    p::N_Vector
    q::N_Vector
    u::N_Vector
    Ap::N_Vector
    vtemp::N_Vector
end

const SUNLinearSolverContent_SPBCGS = Ptr{_SUNLinearSolverContent_SPBCGS}
const SUNSPFGMR_MAXL_DEFAULT = Cint(5)
const SUNSPFGMR_MAXRS_DEFAULT = Cint(0)
const SUNSPFGMR_GSTYPE_DEFAULT = MODIFIED_GS

mutable struct _SUNLinearSolverContent_SPFGMR
    maxl::Cint
    pretype::Cint
    gstype::Cint
    max_restarts::Cint
    numiters::Cint
    resnorm::realtype
    last_flag::Clong
    ATimes::ATimesFn
    ATData::Ptr{Void}
    Psetup::PSetupFn
    Psolve::PSolveFn
    PData::Ptr{Void}
    s1::N_Vector
    s2::N_Vector
    V::Ptr{N_Vector}
    Z::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

const SUNLinearSolverContent_SPFGMR = Ptr{_SUNLinearSolverContent_SPFGMR}
const SUNSPGMR_MAXL_DEFAULT = Cint(5)
const SUNSPGMR_MAXRS_DEFAULT = Cint(0)
const SUNSPGMR_GSTYPE_DEFAULT = MODIFIED_GS

mutable struct _SUNLinearSolverContent_SPGMR
    maxl::Cint
    pretype::Cint
    gstype::Cint
    max_restarts::Cint
    numiters::Cint
    resnorm::realtype
    last_flag::Clong
    ATimes::ATimesFn
    ATData::Ptr{Void}
    Psetup::PSetupFn
    Psolve::PSolveFn
    PData::Ptr{Void}
    s1::N_Vector
    s2::N_Vector
    V::Ptr{N_Vector}
    Hes::Ptr{Ptr{realtype}}
    givens::Ptr{realtype}
    xcor::N_Vector
    yg::Ptr{realtype}
    vtemp::N_Vector
end

const SUNLinearSolverContent_SPGMR = Ptr{_SUNLinearSolverContent_SPGMR}
const SUNSPTFQMR_MAXL_DEFAULT = Cint(5)

mutable struct _SUNLinearSolverContent_SPTFQMR
    maxl::Cint
    pretype::Cint
    numiters::Cint
    resnorm::realtype
    last_flag::Clong
    ATimes::ATimesFn
    ATData::Ptr{Void}
    Psetup::PSetupFn
    Psolve::PSolveFn
    PData::Ptr{Void}
    s1::N_Vector
    s2::N_Vector
    r_star::N_Vector
    q::N_Vector
    d::N_Vector
    v::N_Vector
    p::N_Vector
    r::Ptr{N_Vector}
    u::N_Vector
    vtemp1::N_Vector
    vtemp2::N_Vector
    vtemp3::N_Vector
end

const SUNLinearSolverContent_SPTFQMR = Ptr{_SUNLinearSolverContent_SPTFQMR}

# Skipping MacroDefinition: SM_CONTENT_S ( A ) ( ( SUNMatrixContent_Sparse ) ( A -> content ) ) #
# Skipping MacroDefinition: SM_ROWS_S ( A ) ( SM_CONTENT_S ( A ) -> M ) #
# Skipping MacroDefinition: SM_COLUMNS_S ( A ) ( SM_CONTENT_S ( A ) -> N ) #
# Skipping MacroDefinition: SM_NNZ_S ( A ) ( SM_CONTENT_S ( A ) -> NNZ ) #
# Skipping MacroDefinition: SM_NP_S ( A ) ( SM_CONTENT_S ( A ) -> NP ) #
# Skipping MacroDefinition: SM_SPARSETYPE_S ( A ) ( SM_CONTENT_S ( A ) -> sparsetype ) #
# Skipping MacroDefinition: SM_DATA_S ( A ) ( SM_CONTENT_S ( A ) -> data ) #
# Skipping MacroDefinition: SM_INDEXVALS_S ( A ) ( SM_CONTENT_S ( A ) -> indexvals ) #
# Skipping MacroDefinition: SM_INDEXPTRS_S ( A ) ( SM_CONTENT_S ( A ) -> indexptrs ) /*
# * -----------------------------------------------------------------
# * PART III: functions exported by sunmatrix_sparse
# *
# * CONSTRUCTORS:
# *    SUNSparseMatrix
# *    SUNSparseFromDenseMatrix
# *    SUNSparseFromBandMatrix
# * OTHER:
# *    SUNSparseMatrix_Print
# *    SUNSparseMatrix_Realloc
# *    SUNSparseMatrix_Rows
# *    SUNSparseMatrix_Columns
# *    SUNSparseMatrix_NNZ
# *    SUNSparseMatrix_NP
# *    SUNSparseMatrix_SparseType
# *    SUNSparseMatrix_Data
# *    SUNSparseMatrix_IndexValues
# *    SUNSparseMatrix_IndexPointers
# * -----------------------------------------------------------------
# */

mutable struct _SUNMatrixContent_Sparse
    M::sunindextype
    N::sunindextype
    NNZ::sunindextype
    NP::sunindextype
    data::Ptr{realtype}
    sparsetype::Cint
    indexvals::Ptr{sunindextype}
    indexptrs::Ptr{sunindextype}
    rowvals::Ptr{Ptr{sunindextype}}
    colptrs::Ptr{Ptr{sunindextype}}
    colvals::Ptr{Ptr{sunindextype}}
    rowptrs::Ptr{Ptr{sunindextype}}
end

const SUNMatrixContent_Sparse = Ptr{_SUNMatrixContent_Sparse}


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

abstract type AbstractMeshData{NDIMS, ElemType} end

"""
    VertexMappedMesh{NDIMS, ElemType, Nboundaries, Tv, Ti} <: AbstractMeshData{NDIMS, ElemType}

`VertexMappedMesh` describes a mesh which is constructed by an reference-to-physical
mapping which can be constructed using only the vertex positions.

Wraps `MeshData` and `boundary_faces` in a dispatchable mesh type.
"""
struct VertexMappedMesh{NDIMS, ElemType, MeshDataT <: MeshData{NDIMS}, Nboundaries, BoundaryFaceT} <: AbstractMeshData{NDIMS, ElemType}
  md::MeshDataT
  boundary_faces::BoundaryFaceT
end

Base.ndims(::VertexMappedMesh{NDIMS}) where {NDIMS} = NDIMS

function Base.show(io::IO, mesh::VertexMappedMesh{NDIMS, ElemType}) where {NDIMS, ElemType}
  @nospecialize mesh # reduce precompilation time
  print(io, "$ElemType VertexMappedMesh with NDIMS = $NDIMS.")
end

function Base.show(io::IO, ::MIME"text/plain", mesh::VertexMappedMesh{NDIMS, ElemType, MeshDataT, Nb}) where {NDIMS, ElemType, MeshDataT, Nb}
  @nospecialize mesh # reduce precompilation time
  if get(io, :compact, false)
    show(io, mesh)
  else
    summary_header(io, "VertexMappedMesh{$NDIMS, $ElemType, $Nb}, ")
    summary_line(io, "number of elements", mesh.md.num_elements)
    summary_line(io, "number of boundaries", length(mesh.boundary_faces))
    for (boundary_name, faces) in mesh.boundary_faces
      summary_line(increment_indent(io), "nfaces on $boundary_name", length(faces))
    end
    summary_footer(io)
  end
end

"""
    VertexMappedMesh(vertex_coordinates::NTuple{NDIMS, Vector{Tv}}, EToV, rd::RefElemData;
                     is_on_boundary = nothing,
                     is_periodic::NTuple{NDIMS, Bool} = ntuple(_->false, NDIMS)) where {NDIMS, Tv}

- `vertex_coordinates` is a tuple of vectors containing x,y,... components of the vertex coordinates
- `EToV` is a 2D array containing element-to-vertex connectivities for each element
- `rd` is a `RefElemData` from `StartUpDG.jl`, and contains information associated with to the
  reference element (e.g., quadrature, basis evaluation, differentiation, etc).
- `is_on_boundary` specifies boundary using a `Dict{Symbol, <:Function}`
- `is_periodic` is a tuple of booleans specifying periodicity = `true`/`false` in the (x,y,z) direction.
"""
function VertexMappedMesh(vertex_coordinates::NTuple{NDIMS, Vector{Tv}}, EToV::Array{Ti,2}, rd::RefElemData;
                          is_on_boundary = nothing,
                          is_periodic::NTuple{NDIMS, Bool} = ntuple(_->false, NDIMS)) where {NDIMS, Tv, Ti}

  md = MeshData(vertex_coordinates, EToV, rd)
  md = StartUpDG.make_periodic(md, is_periodic)
  boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary)
  return VertexMappedMesh{NDIMS, typeof(rd.elementType), typeof(md), length(boundary_faces), typeof(boundary_faces)}(md, boundary_faces)
end

# specialization for NDIMS = 1
function VertexMappedMesh(vertex_coordinates::NTuple{1, Vector{Tv}}, EToV::Array{Ti,2}, rd::RefElemData;
                          is_on_boundary = nothing,
                          is_periodic = (false, )) where {Tv, Ti}

  md = MeshData(vertex_coordinates, EToV, rd)
  md = StartUpDG.make_periodic(md, is_periodic...)
  boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary)
  return VertexMappedMesh{1, typeof(rd.elementType), typeof(md), length(boundary_faces), typeof(boundary_faces)}(md, boundary_faces)
end

"""
    VertexMappedMesh(triangulateIO, rd::RefElemData{2, Tri}, boundary_dict::Dict{Symbol, Int})

- `triangulateIO` is a `TriangulateIO` mesh representation
- `rd` is a `RefElemData` from `StartUpDG.jl`, and contains information associated with to the
reference element (e.g., quadrature, basis evaluation, differentiation, etc).
- `boundary_dict` is a `Dict{Symbol, Int}` which associates each integer `TriangulateIO` boundary tag with a Symbol
"""
function VertexMappedMesh(triangulateIO, rd::RefElemData{2, Tri}, boundary_dict::Dict{Symbol, Int})

  vertex_coordinates, EToV = StartUpDG.triangulateIO_to_VXYEToV(triangulateIO)
  md = MeshData(vertex_coordinates, EToV, rd)
  boundary_faces = StartUpDG.tag_boundary_faces(triangulateIO, rd, md, boundary_dict)
  return VertexMappedMesh{2, typeof(rd.elementType), typeof(md), length(boundary_faces), typeof(boundary_faces)}(md, boundary_faces)
end

end # @muladd


text:
module Utils

using Random: shuffle
using Statistics: mean
using Knet: Knet, AutoGrad, Param, @diff, value, params, grad, progress, progress!, relu, load, save, minibatch, Data, dropout, conv4, RNN, pool, mat

export
    MLPClassifier,
    RNNClassifier,
    BiRNNClassifier,
    CNNClassifier,
    EmbedLayer,
    Layer,
    Chain,
    readtext,
    readfolder,
    unzip,
    shuffledata,
    tokenize,
    skipwords,
    frequencies,
    vocabulary,
    max_sentence,
    addpadding!,
    truncate!,
    doc2ids,
    accuracy,
    zeroone,
    convert2d,
    nll,
    sgdupdate!,
    sgd

struct EmbedLayer; w; b; end

EmbedLayer(numwords::Int, embed::Int, scale=0.01) = EmbedLayer(
    Param(scale * randn(embed, numwords)))

function (l::EmbedLayer)(input)
    embed = l.input[:, permutedims(hcat(input...))]
    output = sum(embed, dims=3)[:,:,end]
end

(m::EmbedLayer)(x, y) = nll(m(x), y)

struct Layer; w; b; f; end

Layer(i::Int, o::Int, scale=0.01, f=relu) = Layer(
    Param(scale * randn(o, i)), Param(zeros(o)), f)

(m::Layer)(x) = m.f.(m.w * x .+ m.b)
(m::Layer)(x, y) = nll(m(x), y)

struct Chain
    layers
    Chain(layers...) = new(layers)
end

(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x, y) = nll(c(x), y)

struct MLPClassifier; input; hidden; b1; output; b2; f; pdrop; end

MLPClassifier(input::Int, embed::Int, hidden::Int, output::Int, f=relu; pdrop=0, scale=0.01) =
    MLPClassifier(Param(randn(embed, input) * scale),
                  Param(randn(hidden, embed) * scale),
                  Param(zeros(hidden)),
                  Param(randn(output, hidden) * scale),
                  Param(zeros(output)),
                  f,
                  pdrop)

function (c::MLPClassifier)(input)
    embed = c.input[:, permutedims(hcat(input...))]
    hiddeninput = sum(embed, dims=3)[:,:,end]
    hiddeninput = dropout(hiddeninput, c.pdrop)
    hiddenoutput = c.f.(c.hidden * hiddeninput .+ c.b1)
    hiddenoutput = dropout(hiddenoutput, c.pdrop)
    outputlayer = c.output * hiddenoutput .+ c.b2

    return outputlayer
end

(c::MLPClassifier)(input,output) = nll(c(input), output)

struct RNNClassifier; input; rnn; output; b; pdrop; end

RNNClassifier(input::Int, embed::Int, hidden::Int, output::Int; pdrop=0, scale=0.01) =
    RNNClassifier(Param(Float32.(randn(embed, input) * scale)),
                  RNN(embed, hidden, rnnType=:gru),
                  Param(randn(output, hidden) * scale),
                  Param(zeros(output)),
                  pdrop)

function (c::RNNClassifier)(input)
    embed = c.input[:, permutedims(hcat(input...))]
    embed = dropout(embed, c.pdrop)
    hiddenoutput = c.rnn(embed)
    hiddenoutput = dropout(hiddenoutput, c.pdrop)

    return c.output * hiddenoutput[:,:,end] .+ c.b
end

(c::RNNClassifier)(input,output) = nll(c(input), output)

struct BiRNNClassifier; input; rnn; output; b; pdrop; end

BiRNNClassifier(input::Int, embed::Int, hidden::Int, output::Int; pdrop=0, scale=0.01, ) =
    BiRNNClassifier(Param(Float32.(randn(embed, input) * scale)),
                  RNN(embed, hidden, rnnType=:gru, bidirectional=true),
                  Param(randn(output, 2 * hidden) * scale),
                  Param(zeros(output)),
                  pdrop)

function (c::BiRNNClassifier)(input)
    embed = c.input[:, permutedims(hcat(input...))]
    embed = dropout(embed, c.pdrop)
    hiddenoutput = c.rnn(embed)
    hiddenoutput = dropout(hiddenoutput, c.pdrop)

    return c.output * hiddenoutput[:,:,end] .+ c.b
end

(c::BiRNNClassifier)(input,output) = nll(c(input), output)

struct CNNClassifier; input; w; b1; output; b2; f; pdrop; end

CNNClassifier(input::Int, embed::Int, w1::Int, w2::Int, cx::Int, cy::Int, hidden::Int, output::Int, f=relu; pdrop=0, scale=0.01) =
    CNNClassifier(Param(randn(embed, input) * scale),
                   Param(randn(w1, w2, cx, cy) * scale),
                   Param(zeros(1, 1, cy, 1)),
                   Param(randn(output, hidden) * scale),
                   Param(zeros(output)),
                   f,
                   pdrop)

function (c::CNNClassifier)(input)
    embed = c.input[:, permutedims(hcat(input...), [1, 2])]
    embed = dropout(embed, c.pdrop)
    cnnoutput = c.f.(pool(conv4(c.w, reshape(
        embed, (size(embed,1), size(embed, 2), 1, size(embed, 3)))) .+ c.b1))
    return c.output * mat(cnnoutput) .+ c.b2
end

(c::CNNClassifier)(input,output) = nll(c(input), output)

struct Conv; w; b; f; p; end

(c::Conv)(x) = c.f.(pool(conv4(c.w, dropout(x, c.p)) .+ c.b))

Conv(w1::Int, w2::Int, cx::Int, cy::Int, f=relu; pdrop=0, scale=0.01) =
    Conv(Param(randn(w1, w2, cx, cy) * scale), zeros(1, 1, cy, 1), f, pdrop)

struct Dense; w; b; f; p; end

(d::Dense)(x) = d.f.(d.w * mat(dropout(x, d.p)) .+ d.b)

Dense(i::Int, o::Int, f=relu; pdrop=0) =
    Dense(Param(randn(o,i) * scale),
          Param(zeros(o)),
          f,
          pdrop)

function readtext(filename)
    file = open(filename, "r")
    text = read(file, String)
    close(file)
    text
end

function readfolder(path; class=0)
    files = readdir(path)
    docs = [tokenize(readtext(path * file)) for file in files]
    outputs = [class for i in docs]
    docs, outputs
end

unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))

shuffledata(x, y) = unzip(shuffle(collect(zip(x, y))))

function tokenize(text)
    punctuations = ['.', '?', '!', ';', '(', ')']
    text = replace(text, "<br />" => "")

    for p in punctuations
        text = replace(text, "$p" => " $p")
    end

    tokenized = split(text)
end

function skipwords(x)
    stopwords = ["this", "is", "a", "an", "the",
                 '.', '?', '!', ';']

    removed = [[word for word in doc if !in(word, stopwords)]
               for doc in x]
end

function frequencies(x)
    freqs = Dict{String, Integer}()

    for doc in x
        for word in doc
            freqs[word] = 1 + get(freqs, word, 0)
        end
    end
    freqs
end

function vocabulary(freqs, threshold=1)
    vocabulary = Dict{String, Integer}()
    vocabulary["<pad>"] = 1
    pid = 2

    frequency_order = sort([(freqs[word], word)
                           for word in keys(freqs)], rev=true)

    for (f, w) in frequency_order
        if f >= threshold
            vocabulary[w] = pid
            pid += 1
        end
    end

    vocabulary["<unk>"] = pid
    vocabulary
end

function maxsentence(x)
    maximum(length(doc) for doc in x)
end

function addpadding!(x; pos="post")
    max_length = maxsentence(x)

    for doc in x
        while length(doc) < max_length
            if pos == "post"
                push!(doc, 1)
            else
                pushfirst!(doc, 1)
            end
        end
    end
    x
end

function truncate!(x, th=300; pos="post")
    for doc in x
        while length(doc) > th
            if pos == "post"
                pop!(doc)
            else
                popfirst!(doc)
            end
        end
    end
    x
end

function doc2ids(x, voc)
    x = [map(lowercase, doc) for doc in x]
    ids = [map(x -> haskey(voc, x) ? voc[x] : voc["<unk>"], doc)
           for doc in x]
end

function convert2d(x)
    permutedims(reshape(hcat(x...), (length(x[1]), length(x))))
end

function nll(scores, y)
    expscores = exp.(scores)
    probabilities = expscores ./ sum(expscores, dims=1)
    answerprobs = (probabilities[y[i],i] for i in 1:length(y))
    mean(-log.(answerprobs))
end

function sgdupdate!(func, args; lr=0.1)
    fval = @diff func(args...)
    for param in params(fval)
        ∇param = grad(fval, param)
        param .-= lr * ∇param
    end
    return value(fval)
end

sgd(func, data; lr=0.1) =
    (sgdupdate!(func, args; lr=lr) for args in data)

accuracy(model, x, y) = mean(y' .== map(i->i[1], findmax(Array(model(x)),dims=1)[2]))

accuracy(model, data) = mean(accuracy(model,x,y) for (x,y) in data)

zeroone(x...) = 1 - accuracy(x...)

end

max_stars_repo_name: alperkesen/imdb-sentiment-analysis

max_issues_repo_name: alperkesen/imdb-sentiment-analysis

max_forks_repo_name: alperkesen/imdb-sentiment-analysis


text:
#in this file, op_a and op_b are a,b operators in Barry's original paper.
#H,V means horizontal and vertical
#p means '(prime) for output modes

using SymPy
using Printf
using PyCall
struct QRelaySym{T<:Tuple}
    aH::T
    bH::T
    aV::T
    bV::T
    apH::T
    bpH::T
    apV::T
    bpV::T
end

function trans(op, a, b, ap, bp, B)
    p = B * [ap; bp]
    op = subs(op, a, p[1])
    op = subs(op, b, p[2])
    return op
end

#2D rotation matrix
function rotmat(theta)
    c = cos(theta)
    s = sin(theta)
    return [c s; -s c]
end

function qrelay_op(n, phi, alpha, delta)
    #operators before BS
    aH = symbols(@sprintf("a_H1:%d", n+1))
    bH = symbols(@sprintf("b_H1:%d", n+1))
    aV = symbols(@sprintf("a_V1:%d", n+1))
    bV = symbols(@sprintf("b_V1:%d", n+1))

    op = 0
    for i=1:n
        op += phi[i] * (aH[i]*bH[i] + aV[i]*bV[i])
    end

    #beam splitter transformation
    B = 1/sqrt(2)*[1 1;-1 1]

    #operators after beam splitter
    apH = symbols(@sprintf("a'_H1:%d", n+1))
    bpH = symbols(@sprintf("b'_H1:%d", n+1))
    apV = symbols(@sprintf("a'_V1:%d", n+1))
    bpV = symbols(@sprintf("b'_V1:%d", n+1))

    for i=1:n-1
        op = trans(op, bH[i], aH[i+1], bpH[i], apH[i+1], B)
        op = trans(op, bV[i], aV[i+1], bpV[i], apV[i+1], B)
    end

    op = trans(op, aH[1], aV[1], apH[1], apV[1], rotmat(alpha))
    op = trans(op, bH[n], bV[n], bpH[n], bpV[n], rotmat(delta))

    syms = QRelaySym(aH, bH, aV, bV, apH, bpH, apV, bpV)

    return syms, op
end

function op_mat(op)
    op = op[:as_poly](domain="C")
    op_a = op.x[:gens]
    nab = op[:length]()
    op_ab = ones(SymPy.Sym, nab)
    coef = zeros(Complex, nab)
    mat = zeros(Int64, length(op_a), nab)
    for (i, (ps, c)) in enumerate(op[:as_dict]())
        for (j, p) in enumerate(ps)
            mat[j, i] = p
            op_ab[i] = op_a[j]^p * op_ab[i]
        end
        coef[i] = c
    end

    return op_a, op_ab, mat, coef
end


text:
# This file contains algorithms needed for the solution of
# the Risch differential equation from chapter 6 of the book
#
#    Manuel Bronstein, Symbolic Integration I: Transcendental Functions,
#    2nd ed., Springer 2005.
#

"""
    WeakNormalizer(f, D) -> q

Weak normalization.

Given a field `k`, a derivation `D` on `k[t]` and `f` in `k(t)`, return
`q` in `k[t]` such that `f-D(q)/q` is weakly normalized with respect to D.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.1, p. 183.
"""
function WeakNormalizer(f::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) || error("rational function f must be in the domain of derivation D")
    dn, ds = SplitFactor(denominator(f), D)
    g = gcd(dn, derivative(dn))
    dstar = divexact(dn, g)
    d1 = divexact(dstar, gcd(dstar, g))
    a, b = gcdx(divexact(denominator(f), d1), d1, numerator(f))
    kz, z = PolynomialRing(base_ring(a), :ζ)
    kzt, t = PolynomialRing(kz, var(parent(a)))
    dd1 = D(d1)
    r = resultant(d1(t), a(t)-z*dd1(t))
    if iszero(r)
        return one(parent(a))
    end
    D1 = CoefficientLiftingDerivation(kz, BaseDerivation(D)) # dummy derivation compatible with r
    ns = [numerator(rationalize_over_Int(x)) for x in constant_roots(r, D1) if
            isrational(x) && isone(denominator(x)) && x>0]
    if isempty(ns)
        return one(numerator(f))
    end
    prod([gcd(a - n*dd1, d1)^n for n in ns])
end

"""
    RdeNormalDenominator(f, g, D) -> (a, b, c, h, ρ)

Normal part of the denomninator.

Given a field `k`, a derivation `D` on `k[t]`, `f`, `g` in `k(t)`
with `f` weakly normalized with respect to `t`, return
either `ρ=0`, in which case the equation `D(y)+f*y=g` has no solution in `k⟨t⟩`,
or `ρ=1`, `a`, `h` in `k[t]` and `b`, `c` in `k⟨t⟩`, such that for any  solution `y`
in `k(t)` of `D(y)+f*y=g`, `q=y*h` in `k⟨t⟩` satisfies `a*D(q)+b*q=c`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.1, p. 185.
"""
function RdeNormalDenominator(f::F, g::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) && iscompatible(g, D) ||
        error("rational functions f and g must be in the domain of derivation D")
    # Note: f must be weakly normalized which we do not check. It is recommended
    # to use this function only for rational functions f which were computed by WeakNormalizer.
    (dn, ds) = SplitFactor(denominator(f), D)
    (en, es) = SplitFactor(denominator(g), D)
    p = gcd(dn, en)
    h = divexact(gcd(en, derivative(en)), gcd(p, derivative(p)))
    if !iszero(rem(dn*h^2, en))
        return zero(h), zero(f), zero(f), zero(h), 0 # no solution
    end
    dn*h, dn*h*f-dn*D(h), dn*h^2*g, h, 1
end

function one_times_n_solve(A::Vector{T}, B::Vector{T}) where T<:FieldElement
    if length(A)!=length(B) || length(A)==0
        return T[] # no solution
    end
    eqfound = false
    c = 0
    for i=1:length(A)
        a = A[i]
        b = B[i]
        if iszero(b)
            if !iszero(a)
                return T[] # no solution
            end
        elseif eqfound
            if c!=a//b
                return T[] # no solution
            end
        else
            c = a//b
            eqfound = true
        end
    end
    eqfound ? [c] : T[]
end

"""
    ParametricLogarithmicDerivative(f, w, D) -> (n, m, v, ρ)

Parametric logarithmic derivative heuristic.

Given a field `k`, a derivation `D` on `k[t]`, `f` in `k(t)` and `w=D(θ)/θ` in `k` such that `θ`
is a hyperexponential monomial over `k(t)`, return either `ρ=0`, in which case
`n*f=D(v)/v+m*D(θ)/θ` has no solution `v≠0` in `k(t)` and
integers `n≠0` and `m`, or `ρ=1` and a solution `n`, `m`, `v` of that equation.

This function may throw exception `AlgorithmFailedError`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 7.3, p. 253.
"""
function ParametricLogarithmicDerivative(f::F, w::F, D::Derivation) where F<:FieldElement
    # base case f,w \in constant field, D must be the null derivation
    iszero(D) || error("base case only for null derivations")
    v = one(parent(f))
    q = f//w
    if !isrational(q) # could be complex for example
        return 0, 0, zero(parent(f)), 0
    end
    q = rationalize_over_Int(f//w)
    m = numerator(q)
    n = denominator(q)
    n, m, v, 1
end

function ParametricLogarithmicDerivative(f::F, w::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    # See Bronstein's book, Section 7.3, p. 253
    Z = zero(f)
    no_solution = (0, 0, Z, 0)
    if degree(numerator(f))<=0 && degree(denominator(f))<=0 &&
        degree(numerator(w))<=0 && degree(denominator(w))<=0
        #f, w \in k(t) are actually in k, so search solution not in k(t) but in k
        n, m, v, ρ = ParametricLogarithmicDerivative(
            constant_coefficient(numerator(f)),
            constant_coefficient(numerator(w)),
            BaseDerivation(D))
        return n, m, v + Z, ρ # make v \in k an element of k(t)
    end
    d = denominator(f)
    e = denominator(w)
    p, a = divrem(numerator(f), d)
    q, b = divrem(numerator(w), e)
    t = gen(parent(d))
    B = max(0, degree(D(t))-1)
    C = max(degree(q), degree(p))
    if degree(q)>B
        ss = one_times_n_solve([coeff(p, i) for i=B+1:C], [coeff(q, i) for i=B+1:C])
        if isempty(ss)
            return no_solution
        end
        if !isrational(ss[1])
            return no_solution
        end
        s = rationalize_over_Int(ss[1])
        M = numerator(s)  # In Bronstein's book, M and N are switched, this is the first and so far only bug I found there
        N = denominator(s)
        Q, v, ρ = InFieldLogarithmicDerivativeOfRadical(N*f - M*w, D)
        if ρ>=1 && !iszero(Q) && !iszero(v)
            return Q*N, Q*M, v, 1
        else
            return no_solution
        end
    end
    if degree(p)>B
        return no_solution
    end
    l = lcm(d, e)
    ln, ls = SplitFactor(l, D)
    z = ls*gcd(ln, derivative(ln))
    if degree(z)<=0 # z \in k
        throw(AlgorithmFailedError("ParametricLogarithmicDerivative\n@ $(@__FILE__):$(@__LINE__)"))
    end
    u1, r1 = divrem(numerator(l*f), z)
    u2, r2 = divrem(numerator(l*w), z)
    ss = one_times_n_solve([coeff(r1, i) for i=0:(degree(z)-1)], [coeff(r2, i) for i=0:(degree(z)-1)])
    if isempty(ss)
        return no_solution
    end
    if !isrational(ss[1])
        return no_solution
    end
    s = rationalize_over_Int(ss[1])
    M = numerator(s)
    N = denominator(s)
    Q, v, ρ = InFieldLogarithmicDerivativeOfRadical(N*f - M*w, D)
    if ρ>=1 && !iszero(Q) && !iszero(v)
        return Q*N, Q*M, v, 1
    else
        return no_solution
    end
end

"""
    RdeSpecialDenomExp(a, b, c, D) -> (A, B, C, h)

Special part of the denominator - hyperexponential case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)/t` in `k`,
`a≠0` in `k[t]` with `gcd(a,t)=1`, and `b`, `c` in `k⟨t⟩`
return  `A`, `B`, `C`, `h` in `k[t]` such that
for any solution `q` in `k⟨t⟩` of `a*D(q)+b*q=c`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=C`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.2, p. 190.
"""
function RdeSpecialDenomExp(a::P, b::F, c::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    !iszero(a) || error("a must be != 0")
    ishyperexponential(D) ||
        error("monomial of derivation D must be hyperexponential")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomial a and rational functions b and c must be in the domain of derivation D")
    isreduced(b, D) && isreduced(c, D) ||
        error("rational functions b and c must be reduced with respect to derivation D")
    t = gen(parent(a))
    degree(gcd(t, a))==0 || error("gcd(a, t) must be == 1")
    p = t
    nb = iszero(b) ? typemax(Int) : valuation(b, p)
    nc = iszero(c) ? typemax(Int) : valuation(c, p)
    n = min(0, nc - min(0, nb))
    if nb==0
        α = constant_coefficient(Remainder(-b//a, p))
        w = coeff(MonomialDerivative(D), 1)
        n0, m, z, ρ = ParametricLogarithmicDerivative(α, w, BaseDerivation(D))
        if  ρ>0 && n0==1 && !iszero(z)
            n = min(n, m)
        end
    end
    N = max(0, -nb, n-nc)
    p_power_N = p^N
    b1 = (b+n*a*divexact(D(p), p))*p_power_N
    @assert isone(denominator(b1))
    c1 = c*p^(N-n)
    @assert isone(denominator(c1))
    a*p_power_N, numerator(b1), numerator(c1), p^(-n)
end

"""
    RdeSpecialDenomTan(a, b, c, D) -> (A, B, C, h)

Special part of the denominator - hypertangent case.

Given a field `k` not containing `√-1`, a derivation `D` on `k[t]` with `D(t)/(t^2+1)` in `k`,
`a≠0` in `k[t]` with `gcd(a,t^2+1)=1`, and `b`, `c` in `k⟨t⟩`,
return  `A`, `B`, `C`, `h` in `k[t]` such that
for any solution `q` in `k⟨t⟩` of `a*D(q)+b*q=c`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=C`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.2, p. 192.
"""
function RdeSpecialDenomTan(a::P, b::F, c::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    !iszero(a) || error("a must be != 0")
    ishypertangent(D) ||
        error("monomial of derivation D must be hypertangent")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomial a and rational functions b and c must be in the domain of derivation D")
    isreduced(b, D) && isreduced(c, D) ||
        error("rational functions b and c must be reduced with respect to derivation D")
    t = gen(parent(a))
    p = t^2+1
    degree(gcd(a, p))==0 || error("gcd(a, t^2+1) must be == 1")
    nb = iszero(b) ? typemax(Int) : valuation(b, p)
    nc = iszero(c) ? typemax(Int) : valuation(c, p)
    n = min(0, nc - min(0, nb))
    if nb==0
        αI_plus_β = Remainder(-b//a, p)
        α = coeff(αI_plus_β, 1)
        β = coeff(αI_plus_β, 0)
        η = constant_coefficient(divexact(MonomialDerivative(D), p))
        D0 = BaseDerivation(D)
        v, ρ = InFieldDerivative(2*β, D0)
        if ρ>0 && !iszero(v)
            _, I, D0I = Complexify(parent(η), D0)
            n0, m, v, ρ = ParametricLogarithmicDerivative(α*I+β, 2*η*I, D0I)
            if  ρ>0 && n0==1 && !iszero(v)
                n = min(n, m)
            end
        end
    end
    N = max(0, -nb, n-nc)
    p_power_N = p^N
    b1 = (b+n*a*divexact(D(p), p))*p_power_N
    @assert isone(denominator(b1))
    c1 = c*p^(N-n)
    @assert isone(denominator(c1))
    a*p_power_N, numerator(b1), numerator(c1), p^(-n)
end

"""
    RdeSpecialDenomTanI(a, b, c, D) -> (A, B, C, h)

Special part of the denominator - hypertangent case, field contains I=√-1.

Given a field `k` containing `√-1`, a derivation `D` on `k[t]` with `D(t)/(t^2+1)` in `k`,
`a≠0` in `k[t]` with `gcd(a,t^2+1)=1`, and `b`, `c` in `k⟨t⟩`,
return  `A`, `B`, `C`, `h` in `k[t]` such that
for any solution `q` in `k⟨t⟩` of `a*D(q)+b*q=c`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=C`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

The implementation is  analogous to `RdeSpecialDenomPrim` and `RdeSpecialDenomTan`,
but now with `p=t-√-1` instead of `p=t` resp. `p=t^2+1`,
see [Bronstein's book](https://link.springer.com/book/10.1007/b138171), last sentence
before algorithm `RdeSpecialDenomTan` in Section 6.2, p. 192 and Exercise 6.1, p. 216.
"""
function RdeSpecialDenomTanI(a::P, b::F, c::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    !iszero(a) || error("a must be != 0")
    ishypertangent(D) ||
        error("monomial of derivation D must be hypertangent")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomial a and rational functions b and c must be in the domain of derivation D")
    isreduced(b, D) && isreduced(c, D) ||
        error("rational functions b and c must be reduced with respect to derivation D")
    t = gen(parent(a))
    degree(gcd(a, t^2 + 1))==0 || error("gcd(a, t^2+1) must be == 1")
    contains_I(parent(a)) || error("field k must contain I=sqrt(-1)")
    I = get_I(base_ring(parent(a)))
    p = t-I
    nb = iszero(b) ? typemax(Int) : valuation(b, p)
    nc = iszero(c) ? typemax(Int) : valuation(c, p)
    n = min(0, nc - min(0, nb))
    if nb==0
        α = constant_coefficient(Remainder(-b//a, p))
        η = constant_coefficient(divexact(MonomialDerivative(D), t^2 + 1))
        n0, m, v, ρ = ParametricLogarithmicDerivative(α, 2*η*I, BaseDerivation(D))
        if  ρ>0 && n0==1 && !iszero(v)
            n = min(n, m)
        end
    end
    N = max(0, -nb, n-nc)
    pp = t^2 + 1
    pp_power_N = pp^N
    b1 = (b+n*a*divexact(D(pp), pp))*pp_power_N
    @assert isone(denominator(b1))
    c1 = c*pp^(N-n)
    @assert isone(denominator(c1))
    a*pp_power_N, numerator(b1), numerator(c1), pp^(-n)
end

"""
    RdeSpecialDenominator(a, b, c, D) -> (A, B, C, h)

Special part of the denominator.

Given a field `k`, a derivation `D` on `k[t]`, `a≠0` in `k[t]`, and `b`, `c` in `k⟨t⟩`,
return  `A`, `B`, `C`, `h` in `k[t]` such that
for any solution `q` in `k⟨t⟩` of `a*D(q)+b*q=c`, `r=q*h` in `k[t]` satisfies `A*D(r)+B*r=C`.

(Here, `k⟨t⟩` denotes the elements of `k(t)` which are reduced w.r.t. `D`.)

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.2, p. 186.
"""
function RdeSpecialDenominator(a::P, b::F, c::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    !iszero(a) || error("a must be != 0")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomial a and rational functions b and c must be in the domain of derivation D")
    isreduced(b, D) && isreduced(c, D) ||
        error("rational functions b and c must be reduced with respect to derivation D")
    if isprimitive(D)
        @assert isone(denominator(b))
        @assert isone(denominator(c))
        return a, numerator(b), numerator(c), one(b)
    elseif ishyperexponential(D)
        t = gen(parent(a))
        d = gcd(a, t)
        a = divexact(a, d)
        b = b//d
        c = c//d
        return  RdeSpecialDenomExp(a, b, c, D)
    elseif ishypertangent(D)
        t = gen(parent(a))
        p = t^2 + 1
        d = gcd(a, p)
        a = divexact(a, d)
        b = b//d
        c = c//d
        if contains_I(parent(a))
            return RdeSpecialDenomTanI(a, b, c, D) #case √-1 in k
        else
            return RdeSpecialDenomTan(a, b, c, D)
        end
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("RdeSpecialDenominator: monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    RdeBoundDegreePrim(a, b, c, D) -> n

Bound on polynomial solutions - primitive case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)` in `k`, and `a`, `b`, `c` in `k[t]`
with `a≠0`, return integer `n` such that `degree(q)≤n` for any solution `q` in `k[t]` of `a*D(q)+b*q=c`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.3, p. 198.
"""
function RdeBoundDegreePrim(a::P, b::P, c::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    isprimitive(D) ||
        error("monomial of derivation D must be primitive")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials a, b, and c must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    if db>da
        n = max(0, dc-db)
    else
        n = max(0, dc-da+1)
    end
    if db==da-1
        α = -leading_coefficient(b)//leading_coefficient(a)
        z, m0, ρ = LimitedIntegrate(α, leading_coefficient(D), BaseDerivative(D))
        if ρ>0 && is_rational(m0)
            m = rationalize_over_Int(m0)
            if denominator(m)==1
                n = max(n, numerator(m))
            end
        end
    end
    D0 = BaseDerivation(D)
    if db==da
        α = -leading_coefficient(b)//leading_coefficient(a)
        z, ρ = InFieldDerivative(α, D0)
        if ρ>0 && !iszero(z)
            β = -leading_coefficient(a*D0(z)+b*z)//(z*leading_coefficient(a))
            w, m0, ρ = LimitedIntegrate(β, leading_coefficient(D), D0)
            if ρ>0 && isrational(m0)
                m = rationalize_over_Int(m0)
                if denominator(m)==1
                    n = max(n, numerator(m))
                end
            end
        end
    end
    n
end

"""
    RdeBoundDegreeBase(a, b, c) -> n

Bound on polynomial solutions - base case.

Given a field `k` and  `a`, `b`, `c` in `k[t]` with `a≠0`, return integer `n`
such that `degree(q)≤n` for any solution `q` in `k[t]` of `a*(d/dt)(q)+b*q=c`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.3, p. 199.
"""
function RdeBoundDegreeBase(a::P, b::P, c::P) where
    {T<:FieldElement, P<:PolyElem{T}}
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    n = max(0, dc - max(db, da - 1))
    if db==da-1
        m0 = -leading_coefficient(b)//leading_coefficient(a)
        if isrational(m0)
            m = rationalize_over_Int(m0)
            if isone(denominator(m))
                n = max(0, numerator(m), dc - db)
            end
        end
    end
    return n
end

"""
    RdeBoundDegreeExp(a, b, c, D) -> n

Bound on polynomial solutions - hyperexponential case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)/t` in `k`, and `a`, `b`, `c` in `k[t]`
with`a≠0`, return integer `n` such that `degree(q)≤n` for any solution `q` in `k[t]` of `a*D(q)+b*q=c`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.3, p. 200.
"""
function RdeBoundDegreeExp(a::P, b::P, c::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    ishyperexponential(D) ||
        error("monomial of derivation D must be hyperexponential")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials a, b, and c must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    n = max(0, dc - max(db, da))
    if da==db
        α = -leading_coefficient(b)//leading_coefficient(a)
        w = coeff(MonomialDerivative(D), 1)
        n0, m, z, ρ = ParametricLogarithmicDerivative(α, w, BaseDerivation(D))
        if  ρ>0 && n0==1 && !iszero(z)
            n = max(n, m)
        end
    end
    n
end

"""
    RdeBoundDegreeNonLinear(a, b, c, D) -> n

Bound on polynomial solutions - nonlinear case.

Given a field `k`, a derivation `D` on `k[t]` with `degree(D(t))≥2`, and `a`, `b`, `c` in `k[t]`
with `a≠0`, return integer `n` such that `degree(q)≤n` for any solution `q` in `k[t]` of `a*D(q)+b*q=c`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.3, p. 201.
"""
function RdeBoundDegreeNonLinear(a::P, b::P, c::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    isnonlinear(D) ||
        error("monomial of derivation D must be nonlinear")
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials a, b, and c must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    da = degree(a)
    db = degree(b)
    dc = degree(c)
    δ = degree(D)
    λ = leading_coefficient(D)
    n = max(0, dc - max(da + δ - 1, db))
    if db==da+δ-1
        m0 = -leading_coefficient(b)//(λ*leading_coefficient(a))
        if isrational(m0)
            m = rationalize_over_Int(m0)
            if isone(denominator(m))
                n = max(0, numerator(m), dc - db)
            end
        end
    end
    n
end

"""
    RdeBoundDegree(a, b, c, D) -> n

Bound on polynomial solutions.

Given a field `k`, a derivation `D` on `k[t]` and `a`, `b`, `c` in `k[t]`
with `a≠0`, return integer `n` such that `degree(q)≤n`
for any solution `q` in `k[t]` of `a*D(q)+b*q=c`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.3, p. 193.
"""
function RdeBoundDegree(a::P, b::P, c::P, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials a, b, and c must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    if isbasic(D) || (isprimitive(D) && isone(MonomialDerivative(D)))
        return RdeBoundDegreeBase(a, b, c)
    elseif isprimitive(D)
        return RdeBoundDegreePrim(a, b, c, D)
    elseif ishyperexponential(D)
        return RdeBoundDegreeExp(a, b, c, D)
    elseif isnonlinear(D)
        return RdeBoundDegreeNonLinear(a, b, c, D)
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("RdeBoundDegree: monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    SPDE(a, b, c, D, n) -> (B, C, m, α, β, ρ)

Rothstein's Special Polynomial Differential Equation algorithm.

Given a field `k`, a derivation `D` on `k[t]`, an integer `n`
and `a`, `b`, `c` in `k[t]`  with `a≠0`, return either
`ρ=0`, in which case the equation `a*D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1`, `B`, `C`, `α`, `β` in `k[t]` and an integer `m`, such that any solution `q` in `k[t]`
of `a*D(q)+b*q=c` must be of the form `q=α*h+β`, where `h` is in `k[t]`, `degree(h)≤m` and
`D(h)+B*h=C`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.4, p. 204.
"""
function SPDE(a::P, b::P, c::P, D::Derivation, n::Int) where
    {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(a, D) && iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials a, b, and c must be in the domain of derivation D")
    !iszero(a) || error("polynomial a must be nonzero")
    Z = zero(a)
    no_solution = (Z, Z, 0, Z, Z, 0)
    if n<0
        if iszero(c)
            return Z, Z, 0, Z, Z, 1
        else
            return no_solution
        end
    end
    g = gcd(a, b)
    if !iszero(rem(c, g))
        return no_solution
    end
    a = divexact(a, g)
    b = divexact(b, g)
    c = divexact(c, g)
    if degree(a)==0
        return divexact(b, a), divexact(c, a), n, one(a), Z, 1
    end
    r, z = gcdx(b, a, c)
    u = SPDE(a, b + D(a), z - D(r), D, n - degree(a))
    if u[6]==0 # no solution of recursive call of SPDE
        return no_solution
    end
    b1, c1, m, α, β, _ = u
    b1, c1, m, a*α, a*β+r, 1
end

"""
    PolyRischDENoCancel1(b, c, D[, n=∞]) -> (q, ρ)

Polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]`, an integer `n`
and `b`, `c` in `k[t]`  with `b≠0` such that either `D=d/dt`or
`degree(b)>max(0, degree(D(t))-1)`, return either
`ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.5, p. 208.
"""
function PolyRischDENoCancel1(b::P, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}      # here typemax(Int) represents +infinity
    iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials b and c must be in the domain of derivation D")
    !iszero(b) || error("polynomial b must be nonzero")
    isbasic(D) || degree(b)>max(0, degree(D)-1) ||
        error("either derivation D must be basic or degree(b)>max(0, degree(D)-1)")
    Z  = zero(b)
    t = gen(parent(b))
    q = Z
    while !iszero(c)
        m = degree(c)-degree(b)
        if n<0 || m<0 || m>n
            return Z, 0 # no solution
        end
        p = (leading_coefficient(c)//leading_coefficient(b))*t^m
        q += p
        n = m - 1
        c -= D(p)+b*p
    end
    q, 1
end

"""
    PolyRischDENoCancel2(b, c, D[, n=∞]) -> (q, B, C, ρ)

Polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]`, an integer `n`
and `b`, `c` in `k[t]` such that `degree(b)<degree(D(t))-1` and either
`D=d/dt` or `degree(D(t))≥2`, return either
`ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`,
or `ρ=2` and `q`, `A`, `B` in `k[t]` such that for any solution `y` in `k[t]`
of degree at most `n` of `D(y)+b*y=c`, `z=y-q` is a solution of `D(z)+B*z=C`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.5, p. 209.
"""
function PolyRischDENoCancel2(b::P, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}      # here typemax(Int) represents +infinity
    iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials b and c must be in the domain of derivation D")
    δ = degree(D)
    iszero(b) || degree(b)<δ-1 || error("degree(b)<degree(D)-1 must hold")
    # Note: In Bronstein's book the convention degree(0)=-infinity is used,
    # but in AbstractAlgebra we have degree(0)=-1, so we have to explicitly
    # include the clause iszero(b) in the above condition.
    isbasic(D) || δ>=2 ||
        error("either derivation D must be basic or degree(D)>=2")
    t = gen(parent(b))
    λ = leading_coefficient(D)
    Z  = zero(b)
    Z0 = zero(λ)
    no_solution = (Z, Z0, Z0, 0)
    q = Z
    while !iszero(c)
        if n==0
            m=0
        else
            m = degree(c) - δ + 1
        end
        if n<0 || m<0 || m>n
            return no_solution
        end
        if m>0
            p = (leading_coefficient(c)//(m*λ))*t^m
        else
            if degree(b)!=degree(c)
                return no_solution
            end
            if degree(b)==0
                return q, constant_coefficient(b), constant_coefficient(c), 2
            end
            p = (leading_coefficient(c)//leading_coefficient(b)) + Z #make p\in k an element of k(t)
        end
        q += p
        n = m - 1
        c -= D(p)+b*p
    end
    q, Z0, Z0, 1
end

"""
    PolyRischDENoCancel3(b, c, D[, n=∞]) -> (q, m, C, ρ)

Polynomial Risch differential equation - no cancellation.

Given a field `k`, a derivation `D` on `k[t]`, an integer `n`
and `b`, `c` in `k[t]` such that `degree(D(t))≥2` and `degree(b)=degree(D(t))-1`,
return either `ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`,
or `ρ=2`, `q`, `C` in `k[t]` and an integer `m`such that for any solution `y` in `k[t]`
of degree at most `n` of `D(y)+b*y=c`, `z=y-q` is a solution in `k[t]`
of degree at most `m` of `D(z)+b*z=C`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.5, p. 210.
"""
function PolyRischDENoCancel3(b::P, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}      # here typemax(Int) represents +infinity
    iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials b and c must be in the domain of derivation D")
    δ = degree(D)
    δ>=2 || error("degree(D) must be >= 2")
    degree(b)==δ-1 || error("degree(b)==degree(D)-1 must hold")
    Z = zero(b)
    no_solution = (Z, 0, Z, 0)
    t = gen(parent(b))
    λ = leading_coefficient(D)
    M = -1
    h = -leading_coefficient(b)//λ
    if isrational(h)
        h = rationalize_over_Int(h)
        if isone(denominator(h)) &&  h>0
            M = numerator(h)
        end
    end
    q = Z
    while !iszero(c)
        m = max(M, degree(c)-δ+1)
        if n<0 || m<0 || m>n
            return no_solution
        end
        u = m*λ + leading_coefficient(b)
        if iszero(u)
            return q, m, c, 2
        end
        if m>0
            p = (leading_coefficient(c)//u)*t^m
        else
            if degree(c)!=δ-1
                return no_solution
            end
            p = (leading_coefficient(c)//leading_coefficient(b)) + Z #make p\in k an element of k[t]
        end
        q += p
        n = m - 1
        c -= D(p)+b*p
    end
    q, 0, Z, 1
end

"""
    PolyRischDECancelPrim(b, c, D[, n=∞]) -> (q, ρ)

Polynomial Risch differential equation, cancellation - primitive case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)` in `k`, an integer `n`,
`b≠0` in `k` and `c` in `k[t]`, return either
`ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.6, p. 212.
"""
function PolyRischDECancelPrim(b::T, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}       # here typemax(Int) represents +infinity
    isprimitive(D) ||
        error("monomial of derivation D must be primitive")
    D0 = BaseDerivation(D)
    iscompatible(b, D0) ||
        error("coefficient b must be in the domain of the base derivation of D")
    iscompatible(c, D) ||
        error("polynomial c must be in the domain of derivation D")
    Z = zero(c)
    no_solution = (Z, 0)
    if b==0  #Note: case b==0 allowed, TODO: test this and update docstring
        q0, ρ = InFieldDerivative(c//1, D) # make poly c a rational function
        q = numerator(q0)
        if ρ<=0 || !isone(denominator(q0)) || degree(q)>n
            return no_solution
        end
        return q, 1
    end
    t = gen(parent(c))
    z, ρ = InFieldLogarithmicDerivative(b, D0)
    if ρ>0
        p0, ρ = InFieldDerivative(z*c, D)
        p = numerator(p0)
        if ρ<=0 || !isone(denominator(p0)) || degree(p)>n
            return no_solution
        end
        return p//z, 1
    end
    if iszero(c)
        return Z, 1 # zero is solution
    end
    if n<degree(c)
        return no_solution
    end
    q = Z
    while !iszero(c)
        m = degree(c)
        if n<m
            return no_solution
        end
        s, ρ = RischDE(b, leading_coefficient(c), D0)
        if ρ<=0
            return Z, 0
        end
        q += s*t^m
        n = m - 1
        c -= b*s*t^m + D(s*t^m)
    end
    q, 1
end

"""
    PolyRischDECancelExp(b, c, D[, n=∞]) -> (q, ρ)

Polynomial Risch differential equation, cancellation - hyperexponential case.

Given a field `k`, a derivation `D` on `k[t]` with `D(t)/t` in `k`, an integer `n`,
`b≠0` in `k` and `c` in `k[t]`, return either
`ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.6, p. 213.
"""
function PolyRischDECancelExp(b::T, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}      # here typemax(Int) represents +infinity
    ishyperexponential(D) ||
        error("monomial of derivation D must be hyperexponential")
    D0 = BaseDerivation(D)
    iscompatible(b, D0) ||
        error("coefficient b must be in the domain of the base derivation of D")
    iscompatible(c, D) ||
        error("polynomial c must be in the domain of derivation D")
    Z = zero(c)
    no_solution = (Z, 0)
    if b==0
        q0, ρ = InFieldDerivative(c//1, D) # make poly c a rational function
        q = numerator(q0)
        if ρ<=0 || !isone(denominator(q0)) || degree(q)>n
            return no_solution
        end
        return q, 1
    end
    t = gen(parent(c))
    H = MonomialDerivative(D)
    w = coeff(H,1) # = Dt/t for hyperexponentialm case
    n1, m, z, ρ = ParametricLogarithmicDerivative(b, w, D0)
    if  ρ>0 && n1==1
        p, ρ = InFieldDerivative(c*z*(t//1)^m, D) # t//1 in k(t)
        if ρ<=0 || !isreduced(p, D)
            return no_solution
        end
        q0 = p//(z*(t//1)^m) # t((1 in k(t)))
        q = numerator(q0)
        if isone(denominator(q0)) && degree(q)<=n
            return q, 1
        else
            return no_solution
        end
    end
    if iszero(c)
        return Z, 1 # zero is solution
    end
    if n<degree(c)
        return no_solution
    end
    q = Z
    while !iszero(c)
        m = degree(c)
        if n<m
            return no_solution
        end
        s, ρ = RischDE(b+m*w, leading_coefficient(c), D0)
        if ρ<=0
            return no_solution
        end
        q += s*t^m
        n = m - 1
        c -= b*s*t^m + D(s*t^m)
    end
    q, 1
end

"""
    PolyRischDECancelTan(b₀, c, D[, n=∞]) -> (q, ρ)

Polynomial Risch differential equation, cancellation - hypertangent case.

Given a field `k` not containing `√-1`, a derivation `D` on `k[t]` with `D(t)/(t^2+1)` in `k`, an integer `n`,
`b₀` in `k` and `c` in `k[t]`, return either
`ρ=0`, in which case the equation `D(q)+(b₀-n*t*D(t)/(t^2+1))*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Section 6.6, p. 215.
"""
function PolyRischDECancelTan(b0::T, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}}      # here typemax(Int) represents +infinity
    ishypertangent(D) ||
        error("monomial of derivation D must be hypertangent")
    D0 = BaseDerivation(D)
    iscompatible(b0, D0) ||
        error("coefficient b0 must be in the domain of the base derivation of D")
    iscompatible(c, D) ||
        error("polynomial c must be in the domain of derivation D")
    Z = zero(c)
    no_solution = (Z, 0)
    if n==0
        if degree(c)<=0
            if !iszero(b0)
                q, ρ = RischDE(b0, constant_coefficient(c), BaseDerivation(D))
                if ρ<=0
                    return no_solution
                else
                    return q + Z, 1
                end
            end
            q, ρ = InFieldDerivative(constant_coefficient(c), BaseDerivation(D))
            if ρ<=0
                return no_solution
            else
                return q + Z, 1
            end
        end
    end
    t = gen(parent(c))
    p = t^2 + 1
    η = constant_coefficient(divexact(MonomialDerivative(D), p))
    r = rem(c, p)
    c0 = coeff(r, 0)
    c1 = coeff(r, 1)
    u, v, ρ   = CoupledDESystem(b0, -n*η, c0, c1, BaseDerivation(D))
    if ρ<=0
        return no_solution
    end
    r = u + v*t
    if n==1
        return r, 1
    end
    c = divexact(c - D(r) - (b0 - n*η*t)*r, p)
    h, ρ = PolyRischDECancelTan(b0, c, D, n-2)
    if ρ<=0
        return no_solution
    end
    p*h + r, 1
end

"""
    PolyRischDE(b, c, D[, n=∞]) -> (q, ρ)

Polynomial Risch differential equation.

Given a field `k`, a derivation `D` on `k[t]`, an integer `n` and
`b`, `c` in `k[t]`, return either
`ρ=0`, in which case the equation `D(q)+b*q=c` has no solution of degree at most `n` in `k[t]`,
or `ρ=1` and a solution `q` in `k[t]` of this equation with `degree(q)≤n`.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Sections 6.6 and 6.7, p_power_N. 206-216.
"""
function PolyRischDE(b::P, c::P, D::Derivation, n::Int=typemax(Int)) where
    {T<:FieldElement, P<:PolyElem{T}} # here typemax(Int) represents +infinity
    iscompatible(b, D) && iscompatible(c, D) ||
        error("polynomials b and c must be in the domain of derivation D")
    δ = degree(D)
    Z = zero(parent(b))
    if !iszero(b) && isbasic(D) || degree(b)>max(0, δ-1)
        return PolyRischDENoCancel1(b, c, D, n)
    elseif (iszero(b) || degree(b)<δ-1) && (isbasic(D) || δ>=2)
        q, b, c, ρ = PolyRischDENoCancel2(b, c, D, n)
        if ρ==2
            q1, ρ = RischDE(b, c, BaseDerivation(D))
            ρ>=1 || return Z, ρ
            q = q1 - q
        end
        return q, 1
    elseif δ>=2 && degree(b)==δ-1
        q, m, c, ρ = PolyRischDENoCancel3(b, c, D, n)
        if ρ<=0
            return Z, 0
        elseif ρ==1
            return q, 1
        elseif ρ==2
            if ishypertangent(D)
                if contains_I(parent(b))
                    # Seems that PolyRischDE was called by CoupledDESystem,
                    # this case should be handled by CoupledDESystem.
                    return c, 101 + max(-1, m) # return new c, m
                else
                    t = gen(parent(b))
                    η = divexact(MonomialDerivative(D), t^2+1)
                    b0 = b + n*t*η
                    @assert degree(b0)<=0
                    return PolyRischDECancelTan(constant_coefficient(b0), c, D, n)
                end
            else
                H = MonomialDerivative(D)
                throw(NotImplementedError("PolyRischDE: no cancellation, degree(b)==δ-1, monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
            end
        else
            @assert false # never reach this point
        end
    # At this point only the cancellation case δ<=1, D!=d/dt is possible;
    # this is only compatible with primitive or hyperexponential.
    elseif contains_I(parent(b))
        # Seems that PolyRischDE was called by CoupledDESystem,
        # this case should be handled by CoupledDESystem.
        return c, 101 + max(-1, n) # return original c, n
    elseif isprimitive(D)
        return PolyRischDECancelPrim(constant_coefficient(b), c, D, n)
    elseif ishyperexponential(D)
        return PolyRischDECancelExp(constant_coefficient(b), c, D, n)
    else
        H = MonomialDerivative(D)
        throw(NotImplementedError("PolyRischDE: cancellation case, monomial derivative $H\n@ $(@__FILE__):$(@__LINE__)"))
    end
end

"""
    RischDE(f, g, D) -> (y, ρ)

Risch differential equation.

Given a field `K`, a derivation `D` on `K` and  `f`, `g` in `K`, return
either `ρ=0`, in which case the equation `D(y)+f*q=g` has no solution in `K`,
or `ρ=1` and a solution `y` in `K` of this equation.

See [Bronstein's book](https://link.springer.com/book/10.1007/b138171), Chapter 6, p. 181.
"""
function RischDE(f::F, g::F, D::Derivation) where F<:FieldElement
    iszero(D) || error("base case only for null derivations")
    #base case => pure algebra problem ...
    iscompatible(f, D) && iscompatible(g,D) ||
        error("coefficients f and g must be in the domain of the base derivation of D")
    if iszero(f)
        if iszero(g)
            return zero(parent(f)), 1
        else
            return zero(parent(f)), 0 # no solution
        end
    else
        return g//f, 1
    end
 end

function RischDE(f::F, g::F, D::Derivation) where
    {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    iscompatible(f, D) && iscompatible(g, D) ||
        error("rational functions f and g must be in the domain of derivation D")
    if iszero(f)
        return InFieldDerivative(g::F, D::Derivation)
    end

    Z = zero(f)
    h0 = WeakNormalizer(f, D)

    f = f - D(h0)//h0
    g = h0*g
    a, b, c, h1, ρ = RdeNormalDenominator(f, g, D)
    ρ>=1 || return Z, ρ

    a, b, c, h2 = RdeSpecialDenominator(a, b, c, D)

    n = RdeBoundDegree(a, b, c, D)

    b, c, n, α, β, ρ = SPDE(a, b, c, D, n)
    ρ>=1 || return Z, ρ

    z, ρ = PolyRischDE(b, c, D, n)
    ρ>=1 || return Z, ρ

    (α*z+β)//(h0*h1*h2), 1
end


text:
module SpringCollab2020FloatierSpaceBlock

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/floatierSpaceBlock" FloatierSpaceBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, tiletype::String="3", disableSpawnOffset::Bool=false, floatinessMultiplier::Number=1.0, bounceBackMultiplier::Number=1.0, dashOffsetMultiplier::Number=1.0)

const placements = Ahorn.PlacementDict(
    "Floatier Space Block (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FloatierSpaceBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    )
)

Ahorn.editingOptions(entity::FloatierSpaceBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::FloatierSpaceBlock) = 8, 8
Ahorn.resizable(entity::FloatierSpaceBlock) = true, true

Ahorn.selection(entity::FloatierSpaceBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FloatierSpaceBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end


text:
"""
    _D_inv_gain_method(parameters::Parameters1PL, latent::Latent1D)
"""
function _D_inv_gain_method(parameters::Parameters1PL, latent::Latent1D)
    return 1 / _item_expected_information(parameters, latent)
end

"""
    _D_inv_gain_method(parameters::Parameters2PL, latent::Latent1D)
"""
function _D_inv_gain_method(parameters::Parameters2PL, latent::Latent1D)
    old_exp_info = copy(parameters.expected_information)
    return LinearAlgebra.det(LinearAlgebra.inv(old_exp_info)) -
           LinearAlgebra.det(LinearAlgebra.inv(
        old_exp_info + _item_expected_information(parameters, latent),
    ))
end

"""
    D_inv_gain_method(item::AbstractItem, examinee::AbstractExaminee)

# Description
It computes the gain in the inverse of the expected information of the item.

# Arguments
  - **`item::AbstractItem`**: The item.
  - **`examinee::AbstractExaminee`**: The examinee at which computing the information.

# Output
It returns a `Float64` scalar.

# References
__Yinhong He & Ping Chen, 2020. "Optimal Online Calibration Designs for Item Replenishment in Adaptive Testing," Psychometrika, Springer;The Psychometric Society, vol. 85(1), pages 35-55, March.__
"""
function D_inv_gain_method(item::AbstractItem, examinee::AbstractExaminee)
    return _D_inv_gain_method(item.parameters, examinee.latent)
end


text:
"""
    has_standard_linear_indexing(A)

yields whether array `A` implements standard linear indexing (1-based).

"""
has_standard_linear_indexing(A::Array) = true
has_standard_linear_indexing(A::AbstractArray) =
    is_one_based_unit_range(eachindex(A))

"""
    is_one_based_unit_range(itr)

yields whether iterator `itr` is a 1-based unit range.

"""
is_one_based_unit_range(itr) = false
is_one_based_unit_range(itr::Base.OneTo) = true
is_one_based_unit_range(itr::AbstractUnitRange{T}) where {T} =
    first(itr) == oneunit(T)

"""
    to_eltype(A, x)

lazily yields `x` converted to the type of the elements of `A`.

"""
to_eltype(A, x) = to_type(eltype(A), x)

"""
    in_range(i, len::Integer)

yields whether `1 ≤ i ≤ len`.

    in_range(i, A::Array)

yields whether `i` is a valid linear index of array `A`.

    in_range(i, R::AbstractUnitRange{<:Integer})

yields whether `i` is in the range `R`.

"""
in_range(i::Integer, len::Integer) = ((i % UInt) - 1 < (len % UInt))
in_range(i::Integer, A::Array) = in_range(i, length(A))
in_range(i::Integer, R::OneTo) = in_range(i, length(R))
in_range(i::Integer, R::AbstractUnitRange{<:Integer}) = (i ∈ R)

"""
    has_bad_values(A, isbad)

yields whether array `A` has bad values according to predicate `isbad`.  For
arrays with floating-point values, `isbad` default to `isnan` if unspecified.

"""
function has_bad_values(A::AbstractArray, isbad)
    flag = false
    @inbounds @simd for i in eachindex(A)
        flag |= isbad(A[i])
    end
    return flag
end

has_bad_values(A::AbstractArray{<:Integer}) = false
has_bad_values(A::AbstractArray{<:AbstractFloat}) =
    has_bad_values(A, isnan)

"""
    typename(x)

yields a short string describing the type of object `x`.  Argument may also be
the object type.

"""
typename(x::Any) = typename(typeof(x))
typename(T::DataType) = string(nameof(T))

require_one_based_indexing(A...) =
    has_offset_axes(A...) && throw_argument_error(
        "arrays must have 1-based indexing")

throw_argument_error(msg::AbstractString) = throw(ArgumentError(msg))
@noinline throw_argument_error(args...) =
    throw_argument_error(string(args...))

throw_dimension_mismatch(msg::AbstractString) = throw(DimensionMismatch(msg))
@noinline throw_dimension_mismatch(args...) =
    throw_dimension_mismatch(string(args...))