text:
# Julia wrapper for header: ceed.h
# Automatically generated using Clang.jl
#! format: off

function CeedInit(resource, ceed)
    ccall((:CeedInit, libceed), Cint, (Cstring, Ptr{Ceed}), resource, ceed)
end

function CeedGetResource(ceed, resource)
    ccall((:CeedGetResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource)
end

function CeedIsDeterministic(ceed, isDeterministic)
    ccall((:CeedIsDeterministic, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDeterministic)
end

function CeedView(ceed, stream)
    ccall((:CeedView, libceed), Cint, (Ceed, Ptr{FILE}), ceed, stream)
end

function CeedDestroy(ceed)
    ccall((:CeedDestroy, libceed), Cint, (Ptr{Ceed},), ceed)
end

function CeedErrorReturn(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:CeedErrorReturn, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function CeedErrorStore(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:CeedErrorStore, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function CeedErrorAbort(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:CeedErrorAbort, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function CeedErrorExit(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    ccall((:CeedErrorExit, libceed), Cint, (Ceed, Cstring, Cint, Cstring, Cint, Cstring, Ptr{Cvoid}), arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function CeedSetErrorHandler(ceed, eh)
    ccall((:CeedSetErrorHandler, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, eh)
end

function CeedGetErrorMessage(arg1, errmsg)
    ccall((:CeedGetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg)
end

function CeedResetErrorMessage(arg1, errmsg)
    ccall((:CeedResetErrorMessage, libceed), Cint, (Ceed, Ptr{Cstring}), arg1, errmsg)
end

function CeedGetPreferredMemType(ceed, type)
    ccall((:CeedGetPreferredMemType, libceed), Cint, (Ceed, Ptr{CeedMemType}), ceed, type)
end

function CeedVectorCreate(ceed, len, vec)
    ccall((:CeedVectorCreate, libceed), Cint, (Ceed, CeedInt, Ptr{CeedVector}), ceed, len, vec)
end

function CeedVectorSetArray(vec, mtype, cmode, array)
    ccall((:CeedVectorSetArray, libceed), Cint, (CeedVector, CeedMemType, CeedCopyMode, Ptr{CeedScalar}), vec, mtype, cmode, array)
end

function CeedVectorSetValue(vec, value)
    ccall((:CeedVectorSetValue, libceed), Cint, (CeedVector, CeedScalar), vec, value)
end

function CeedVectorSyncArray(vec, mtype)
    ccall((:CeedVectorSyncArray, libceed), Cint, (CeedVector, CeedMemType), vec, mtype)
end

function CeedVectorTakeArray(vec, mtype, array)
    ccall((:CeedVectorTakeArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
end

function CeedVectorGetArray(vec, mtype, array)
    ccall((:CeedVectorGetArray, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
end

function CeedVectorGetArrayRead(vec, mtype, array)
    ccall((:CeedVectorGetArrayRead, libceed), Cint, (CeedVector, CeedMemType, Ptr{Ptr{CeedScalar}}), vec, mtype, array)
end

function CeedVectorRestoreArray(vec, array)
    ccall((:CeedVectorRestoreArray, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array)
end

function CeedVectorRestoreArrayRead(vec, array)
    ccall((:CeedVectorRestoreArrayRead, libceed), Cint, (CeedVector, Ptr{Ptr{CeedScalar}}), vec, array)
end

function CeedVectorNorm(vec, type, norm)
    ccall((:CeedVectorNorm, libceed), Cint, (CeedVector, CeedNormType, Ptr{CeedScalar}), vec, type, norm)
end

function CeedVectorReciprocal(vec)
    ccall((:CeedVectorReciprocal, libceed), Cint, (CeedVector,), vec)
end

function CeedVectorView(vec, fpfmt, stream)
    ccall((:CeedVectorView, libceed), Cint, (CeedVector, Cstring, Ptr{FILE}), vec, fpfmt, stream)
end

function CeedVectorGetLength(vec, length)
    ccall((:CeedVectorGetLength, libceed), Cint, (CeedVector, Ptr{CeedInt}), vec, length)
end

function CeedVectorDestroy(vec)
    ccall((:CeedVectorDestroy, libceed), Cint, (Ptr{CeedVector},), vec)
end

function CeedRequestWait(req)
    ccall((:CeedRequestWait, libceed), Cint, (Ptr{CeedRequest},), req)
end

function CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
    ccall((:CeedElemRestrictionCreate, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
end

function CeedElemRestrictionCreateStrided(ceed, nelem, elemsize, ncomp, lsize, strides, rstr)
    ccall((:CeedElemRestrictionCreateStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, ncomp, lsize, strides, rstr)
end

function CeedElemRestrictionCreateBlocked(ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
    ccall((:CeedElemRestrictionCreateBlocked, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedMemType, CeedCopyMode, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, mtype, cmode, offsets, rstr)
end

function CeedElemRestrictionCreateBlockedStrided(ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr)
    ccall((:CeedElemRestrictionCreateBlockedStrided, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedInt}, Ptr{CeedElemRestriction}), ceed, nelem, elemsize, blksize, ncomp, lsize, strides, rstr)
end

function CeedElemRestrictionCreateVector(rstr, lvec, evec)
    ccall((:CeedElemRestrictionCreateVector, libceed), Cint, (CeedElemRestriction, Ptr{CeedVector}, Ptr{CeedVector}), rstr, lvec, evec)
end

function CeedElemRestrictionApply(rstr, tmode, u, ru, request)
    ccall((:CeedElemRestrictionApply, libceed), Cint, (CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, tmode, u, ru, request)
end

function CeedElemRestrictionApplyBlock(rstr, block, tmode, u, ru, request)
    ccall((:CeedElemRestrictionApplyBlock, libceed), Cint, (CeedElemRestriction, CeedInt, CeedTransposeMode, CeedVector, CeedVector, Ptr{CeedRequest}), rstr, block, tmode, u, ru, request)
end

function CeedElemRestrictionGetCompStride(rstr, compstride)
    ccall((:CeedElemRestrictionGetCompStride, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, compstride)
end

function CeedElemRestrictionGetNumElements(rstr, numelem)
    ccall((:CeedElemRestrictionGetNumElements, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numelem)
end

function CeedElemRestrictionGetElementSize(rstr, elemsize)
    ccall((:CeedElemRestrictionGetElementSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, elemsize)
end

function CeedElemRestrictionGetLVectorSize(rstr, lsize)
    ccall((:CeedElemRestrictionGetLVectorSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, lsize)
end

function CeedElemRestrictionGetNumComponents(rstr, numcomp)
    ccall((:CeedElemRestrictionGetNumComponents, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numcomp)
end

function CeedElemRestrictionGetNumBlocks(rstr, numblk)
    ccall((:CeedElemRestrictionGetNumBlocks, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, numblk)
end

function CeedElemRestrictionGetBlockSize(rstr, blksize)
    ccall((:CeedElemRestrictionGetBlockSize, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, blksize)
end

function CeedElemRestrictionGetMultiplicity(rstr, mult)
    ccall((:CeedElemRestrictionGetMultiplicity, libceed), Cint, (CeedElemRestriction, CeedVector), rstr, mult)
end

function CeedElemRestrictionView(rstr, stream)
    ccall((:CeedElemRestrictionView, libceed), Cint, (CeedElemRestriction, Ptr{FILE}), rstr, stream)
end

function CeedElemRestrictionDestroy(rstr)
    ccall((:CeedElemRestrictionDestroy, libceed), Cint, (Ptr{CeedElemRestriction},), rstr)
end

function CeedBasisCreateTensorH1Lagrange(ceed, dim, ncomp, P, Q, qmode, basis)
    ccall((:CeedBasisCreateTensorH1Lagrange, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, CeedQuadMode, Ptr{CeedBasis}), ceed, dim, ncomp, P, Q, qmode, basis)
end

function CeedBasisCreateTensorH1(ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis)
    ccall((:CeedBasisCreateTensorH1, libceed), Cint, (Ceed, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, basis)
end

function CeedBasisCreateH1(ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis)
    ccall((:CeedBasisCreateH1, libceed), Cint, (Ceed, CeedElemTopology, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedBasis}), ceed, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, basis)
end

function CeedBasisView(basis, stream)
    ccall((:CeedBasisView, libceed), Cint, (CeedBasis, Ptr{FILE}), basis, stream)
end

function CeedBasisApply(basis, nelem, tmode, emode, u, v)
    ccall((:CeedBasisApply, libceed), Cint, (CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector), basis, nelem, tmode, emode, u, v)
end

function CeedBasisGetDimension(basis, dim)
    ccall((:CeedBasisGetDimension, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, dim)
end

function CeedBasisGetTopology(basis, topo)
    ccall((:CeedBasisGetTopology, libceed), Cint, (CeedBasis, Ptr{CeedElemTopology}), basis, topo)
end

function CeedBasisGetNumComponents(basis, numcomp)
    ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, numcomp)
end

function CeedBasisGetNumNodes(basis, P)
    ccall((:CeedBasisGetNumNodes, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P)
end

function CeedBasisGetNumNodes1D(basis, P1d)
    ccall((:CeedBasisGetNumNodes1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P1d)
end

function CeedBasisGetNumQuadraturePoints(basis, Q)
    ccall((:CeedBasisGetNumQuadraturePoints, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q)
end

function CeedBasisGetNumQuadraturePoints1D(basis, Q1d)
    ccall((:CeedBasisGetNumQuadraturePoints1D, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, Q1d)
end

function CeedBasisGetQRef(basis, qref)
    ccall((:CeedBasisGetQRef, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qref)
end

function CeedBasisGetQWeights(basis, qweight)
    ccall((:CeedBasisGetQWeights, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, qweight)
end

function CeedBasisGetInterp(basis, interp)
    ccall((:CeedBasisGetInterp, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp)
end

function CeedBasisGetInterp1D(basis, interp1d)
    ccall((:CeedBasisGetInterp1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, interp1d)
end

function CeedBasisGetGrad(basis, grad)
    ccall((:CeedBasisGetGrad, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad)
end

function CeedBasisGetGrad1D(basis, grad1d)
    ccall((:CeedBasisGetGrad1D, libceed), Cint, (CeedBasis, Ptr{Ptr{CeedScalar}}), basis, grad1d)
end

function CeedBasisDestroy(basis)
    ccall((:CeedBasisDestroy, libceed), Cint, (Ptr{CeedBasis},), basis)
end

function CeedGaussQuadrature(Q, qref1d, qweight1d)
    ccall((:CeedGaussQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d)
end

function CeedLobattoQuadrature(Q, qref1d, qweight1d)
    ccall((:CeedLobattoQuadrature, libceed), Cint, (CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), Q, qref1d, qweight1d)
end

function CeedQRFactorization(ceed, mat, tau, m, n)
    ccall((:CeedQRFactorization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt), ceed, mat, tau, m, n)
end

function CeedSymmetricSchurDecomposition(ceed, mat, lambda, n)
    ccall((:CeedSymmetricSchurDecomposition, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, mat, lambda, n)
end

function CeedSimultaneousDiagonalization(ceed, matA, matB, x, lambda, n)
    ccall((:CeedSimultaneousDiagonalization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, matA, matB, x, lambda, n)
end

function CeedQFunctionCreateInterior(ceed, vlength, f, source, qf)
    ccall((:CeedQFunctionCreateInterior, libceed), Cint, (Ceed, CeedInt, CeedQFunctionUser, Cstring, Ptr{CeedQFunction}), ceed, vlength, f, source, qf)
end

function CeedQFunctionCreateInteriorByName(ceed, name, qf)
    ccall((:CeedQFunctionCreateInteriorByName, libceed), Cint, (Ceed, Cstring, Ptr{CeedQFunction}), ceed, name, qf)
end

function CeedQFunctionCreateIdentity(ceed, size, inmode, outmode, qf)
    ccall((:CeedQFunctionCreateIdentity, libceed), Cint, (Ceed, CeedInt, CeedEvalMode, CeedEvalMode, Ptr{CeedQFunction}), ceed, size, inmode, outmode, qf)
end

function CeedQFunctionAddInput(qf, fieldname, size, emode)
    ccall((:CeedQFunctionAddInput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode)
end

function CeedQFunctionAddOutput(qf, fieldname, size, emode)
    ccall((:CeedQFunctionAddOutput, libceed), Cint, (CeedQFunction, Cstring, CeedInt, CeedEvalMode), qf, fieldname, size, emode)
end

function CeedQFunctionSetContext(qf, ctx)
    ccall((:CeedQFunctionSetContext, libceed), Cint, (CeedQFunction, CeedQFunctionContext), qf, ctx)
end

function CeedQFunctionView(qf, stream)
    ccall((:CeedQFunctionView, libceed), Cint, (CeedQFunction, Ptr{FILE}), qf, stream)
end

function CeedQFunctionApply(qf, Q, u, v)
    ccall((:CeedQFunctionApply, libceed), Cint, (CeedQFunction, CeedInt, Ptr{CeedVector}, Ptr{CeedVector}), qf, Q, u, v)
end

function CeedQFunctionDestroy(qf)
    ccall((:CeedQFunctionDestroy, libceed), Cint, (Ptr{CeedQFunction},), qf)
end

function CeedQFunctionContextCreate(ceed, ctx)
    ccall((:CeedQFunctionContextCreate, libceed), Cint, (Ceed, Ptr{CeedQFunctionContext}), ceed, ctx)
end

function CeedQFunctionContextSetData(ctx, mtype, cmode, size, data)
    ccall((:CeedQFunctionContextSetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, CeedCopyMode, Csize_t, Ptr{Cvoid}), ctx, mtype, cmode, size, data)
end

function CeedQFunctionContextGetData(ctx, mtype, data)
    ccall((:CeedQFunctionContextGetData, libceed), Cint, (CeedQFunctionContext, CeedMemType, Ptr{Cvoid}), ctx, mtype, data)
end

function CeedQFunctionContextRestoreData(ctx, data)
    ccall((:CeedQFunctionContextRestoreData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
end

function CeedQFunctionContextView(ctx, stream)
    ccall((:CeedQFunctionContextView, libceed), Cint, (CeedQFunctionContext, Ptr{FILE}), ctx, stream)
end

function CeedQFunctionContextDestroy(ctx)
    ccall((:CeedQFunctionContextDestroy, libceed), Cint, (Ptr{CeedQFunctionContext},), ctx)
end

function CeedOperatorCreate(ceed, qf, dqf, dqfT, op)
    ccall((:CeedOperatorCreate, libceed), Cint, (Ceed, CeedQFunction, CeedQFunction, CeedQFunction, Ptr{CeedOperator}), ceed, qf, dqf, dqfT, op)
end

function CeedCompositeOperatorCreate(ceed, op)
    ccall((:CeedCompositeOperatorCreate, libceed), Cint, (Ceed, Ptr{CeedOperator}), ceed, op)
end

function CeedOperatorSetField(op, fieldname, r, b, v)
    ccall((:CeedOperatorSetField, libceed), Cint, (CeedOperator, Cstring, CeedElemRestriction, CeedBasis, CeedVector), op, fieldname, r, b, v)
end

function CeedCompositeOperatorAddSub(compositeop, subop)
    ccall((:CeedCompositeOperatorAddSub, libceed), Cint, (CeedOperator, CeedOperator), compositeop, subop)
end

function CeedOperatorLinearAssembleQFunction(op, assembled, rstr, request)
    ccall((:CeedOperatorLinearAssembleQFunction, libceed), Cint, (CeedOperator, Ptr{CeedVector}, Ptr{CeedElemRestriction}, Ptr{CeedRequest}), op, assembled, rstr, request)
end

function CeedOperatorLinearAssembleDiagonal(op, assembled, request)
    ccall((:CeedOperatorLinearAssembleDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
end

function CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)
    ccall((:CeedOperatorLinearAssembleAddDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
end

function CeedOperatorLinearAssemblePointBlockDiagonal(op, assembled, request)
    ccall((:CeedOperatorLinearAssemblePointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
end

function CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)
    ccall((:CeedOperatorLinearAssembleAddPointBlockDiagonal, libceed), Cint, (CeedOperator, CeedVector, Ptr{CeedRequest}), op, assembled, request)
end

function CeedOperatorMultigridLevelCreate(opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict)
    ccall((:CeedOperatorMultigridLevelCreate, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, opCoarse, opProlong, opRestrict)
end

function CeedOperatorMultigridLevelCreateTensorH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
    ccall((:CeedOperatorMultigridLevelCreateTensorH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
end

function CeedOperatorMultigridLevelCreateH1(opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
    ccall((:CeedOperatorMultigridLevelCreateH1, libceed), Cint, (CeedOperator, CeedVector, CeedElemRestriction, CeedBasis, Ptr{CeedScalar}, Ptr{CeedOperator}, Ptr{CeedOperator}, Ptr{CeedOperator}), opFine, PMultFine, rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict)
end

function CeedOperatorCreateFDMElementInverse(op, fdminv, request)
    ccall((:CeedOperatorCreateFDMElementInverse, libceed), Cint, (CeedOperator, Ptr{CeedOperator}, Ptr{CeedRequest}), op, fdminv, request)
end

function CeedOperatorView(op, stream)
    ccall((:CeedOperatorView, libceed), Cint, (CeedOperator, Ptr{FILE}), op, stream)
end

function CeedOperatorApply(op, in, out, request)
    ccall((:CeedOperatorApply, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request)
end

function CeedOperatorApplyAdd(op, in, out, request)
    ccall((:CeedOperatorApplyAdd, libceed), Cint, (CeedOperator, CeedVector, CeedVector, Ptr{CeedRequest}), op, in, out, request)
end

function CeedOperatorDestroy(op)
    ccall((:CeedOperatorDestroy, libceed), Cint, (Ptr{CeedOperator},), op)
end

function CeedIntPow(base, power)
    ccall((:CeedIntPow, libceed), CeedInt, (CeedInt, CeedInt), base, power)
end

function CeedIntMin(a, b)
    ccall((:CeedIntMin, libceed), CeedInt, (CeedInt, CeedInt), a, b)
end

function CeedIntMax(a, b)
    ccall((:CeedIntMax, libceed), CeedInt, (CeedInt, CeedInt), a, b)
end
# Julia wrapper for header: ceed-cuda.h
# Automatically generated using Clang.jl

function CeedQFunctionSetCUDAUserFunction(qf, f)
    ccall((:CeedQFunctionSetCUDAUserFunction, libceed), Cint, (CeedQFunction, Cint), qf, f)
end
# Julia wrapper for header: ceed-backend.h
# Automatically generated using Clang.jl

function CeedMallocArray(n, unit, p)
    ccall((:CeedMallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
end

function CeedCallocArray(n, unit, p)
    ccall((:CeedCallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
end

function CeedReallocArray(n, unit, p)
    ccall((:CeedReallocArray, libceed), Cint, (Csize_t, Csize_t, Ptr{Cvoid}), n, unit, p)
end

function CeedFree(p)
    ccall((:CeedFree, libceed), Cint, (Ptr{Cvoid},), p)
end

function CeedRegister(prefix, init, priority)
    ccall((:CeedRegister, libceed), Cint, (Cstring, Ptr{Cvoid}, UInt32), prefix, init, priority)
end

function CeedIsDebug(ceed, isDebug)
    ccall((:CeedIsDebug, libceed), Cint, (Ceed, Ptr{Bool}), ceed, isDebug)
end

function CeedGetParent(ceed, parent)
    ccall((:CeedGetParent, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent)
end

function CeedGetDelegate(ceed, delegate)
    ccall((:CeedGetDelegate, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, delegate)
end

function CeedSetDelegate(ceed, delegate)
    ccall((:CeedSetDelegate, libceed), Cint, (Ceed, Ceed), ceed, delegate)
end

function CeedGetObjectDelegate(ceed, delegate, objname)
    ccall((:CeedGetObjectDelegate, libceed), Cint, (Ceed, Ptr{Ceed}, Cstring), ceed, delegate, objname)
end

function CeedSetObjectDelegate(ceed, delegate, objname)
    ccall((:CeedSetObjectDelegate, libceed), Cint, (Ceed, Ceed, Cstring), ceed, delegate, objname)
end

function CeedGetOperatorFallbackResource(ceed, resource)
    ccall((:CeedGetOperatorFallbackResource, libceed), Cint, (Ceed, Ptr{Cstring}), ceed, resource)
end

function CeedSetOperatorFallbackResource(ceed, resource)
    ccall((:CeedSetOperatorFallbackResource, libceed), Cint, (Ceed, Cstring), ceed, resource)
end

function CeedGetOperatorFallbackParentCeed(ceed, parent)
    ccall((:CeedGetOperatorFallbackParentCeed, libceed), Cint, (Ceed, Ptr{Ceed}), ceed, parent)
end

function CeedSetDeterministic(ceed, isDeterministic)
    ccall((:CeedSetDeterministic, libceed), Cint, (Ceed, Bool), ceed, isDeterministic)
end

function CeedSetBackendFunction(ceed, type, object, fname, f)
    ccall((:CeedSetBackendFunction, libceed), Cint, (Ceed, Cstring, Ptr{Cvoid}, Cstring, Ptr{Cvoid}), ceed, type, object, fname, f)
end

function CeedGetData(ceed, data)
    ccall((:CeedGetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data)
end

function CeedSetData(ceed, data)
    ccall((:CeedSetData, libceed), Cint, (Ceed, Ptr{Cvoid}), ceed, data)
end

function CeedVectorGetCeed(vec, ceed)
    ccall((:CeedVectorGetCeed, libceed), Cint, (CeedVector, Ptr{Ceed}), vec, ceed)
end

function CeedVectorGetState(vec, state)
    ccall((:CeedVectorGetState, libceed), Cint, (CeedVector, Ptr{UInt64}), vec, state)
end

function CeedVectorAddReference(vec)
    ccall((:CeedVectorAddReference, libceed), Cint, (CeedVector,), vec)
end

function CeedVectorGetData(vec, data)
    ccall((:CeedVectorGetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data)
end

function CeedVectorSetData(vec, data)
    ccall((:CeedVectorSetData, libceed), Cint, (CeedVector, Ptr{Cvoid}), vec, data)
end

function CeedElemRestrictionGetCeed(rstr, ceed)
    ccall((:CeedElemRestrictionGetCeed, libceed), Cint, (CeedElemRestriction, Ptr{Ceed}), rstr, ceed)
end

function CeedElemRestrictionGetStrides(rstr, strides)
    ccall((:CeedElemRestrictionGetStrides, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, strides)
end

function CeedElemRestrictionGetOffsets(rstr, mtype, offsets)
    ccall((:CeedElemRestrictionGetOffsets, libceed), Cint, (CeedElemRestriction, CeedMemType, Ptr{Ptr{CeedInt}}), rstr, mtype, offsets)
end

function CeedElemRestrictionRestoreOffsets(rstr, offsets)
    ccall((:CeedElemRestrictionRestoreOffsets, libceed), Cint, (CeedElemRestriction, Ptr{Ptr{CeedInt}}), rstr, offsets)
end

function CeedElemRestrictionIsStrided(rstr, isstrided)
    ccall((:CeedElemRestrictionIsStrided, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, isstrided)
end

function CeedElemRestrictionHasBackendStrides(rstr, hasbackendstrides)
    ccall((:CeedElemRestrictionHasBackendStrides, libceed), Cint, (CeedElemRestriction, Ptr{Bool}), rstr, hasbackendstrides)
end

function CeedElemRestrictionGetELayout(rstr, layout)
    ccall((:CeedElemRestrictionGetELayout, libceed), Cint, (CeedElemRestriction, Ptr{NTuple{3, CeedInt}}), rstr, layout)
end

function CeedElemRestrictionSetELayout(rstr, layout)
    ccall((:CeedElemRestrictionSetELayout, libceed), Cint, (CeedElemRestriction, Ptr{CeedInt}), rstr, layout)
end

function CeedElemRestrictionGetData(rstr, data)
    ccall((:CeedElemRestrictionGetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data)
end

function CeedElemRestrictionSetData(rstr, data)
    ccall((:CeedElemRestrictionSetData, libceed), Cint, (CeedElemRestriction, Ptr{Cvoid}), rstr, data)
end

function CeedBasisGetCollocatedGrad(basis, colograd1d)
    ccall((:CeedBasisGetCollocatedGrad, libceed), Cint, (CeedBasis, Ptr{CeedScalar}), basis, colograd1d)
end

function CeedHouseholderApplyQ(A, Q, tau, tmode, m, n, k, row, col)
    ccall((:CeedHouseholderApplyQ, libceed), Cint, (Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt), A, Q, tau, tmode, m, n, k, row, col)
end

function CeedBasisGetCeed(basis, ceed)
    ccall((:CeedBasisGetCeed, libceed), Cint, (CeedBasis, Ptr{Ceed}), basis, ceed)
end

function CeedBasisIsTensor(basis, istensor)
    ccall((:CeedBasisIsTensor, libceed), Cint, (CeedBasis, Ptr{Bool}), basis, istensor)
end

function CeedBasisGetData(basis, data)
    ccall((:CeedBasisGetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data)
end

function CeedBasisSetData(basis, data)
    ccall((:CeedBasisSetData, libceed), Cint, (CeedBasis, Ptr{Cvoid}), basis, data)
end

function CeedBasisGetTopologyDimension(topo, dim)
    ccall((:CeedBasisGetTopologyDimension, libceed), Cint, (CeedElemTopology, Ptr{CeedInt}), topo, dim)
end

function CeedBasisGetTensorContract(basis, contract)
    ccall((:CeedBasisGetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract)
end

function CeedBasisSetTensorContract(basis, contract)
    ccall((:CeedBasisSetTensorContract, libceed), Cint, (CeedBasis, Ptr{CeedTensorContract}), basis, contract)
end

function CeedTensorContractCreate(ceed, basis, contract)
    ccall((:CeedTensorContractCreate, libceed), Cint, (Ceed, CeedBasis, Ptr{CeedTensorContract}), ceed, basis, contract)
end

function CeedTensorContractApply(contract, A, B, C, J, t, tmode, Add, u, v)
    ccall((:CeedTensorContractApply, libceed), Cint, (CeedTensorContract, CeedInt, CeedInt, CeedInt, CeedInt, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, Ptr{CeedScalar}, Ptr{CeedScalar}), contract, A, B, C, J, t, tmode, Add, u, v)
end

function CeedTensorContractGetCeed(contract, ceed)
    ccall((:CeedTensorContractGetCeed, libceed), Cint, (CeedTensorContract, Ptr{Ceed}), contract, ceed)
end

function CeedTensorContractGetData(contract, data)
    ccall((:CeedTensorContractGetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data)
end

function CeedTensorContractSetData(contract, data)
    ccall((:CeedTensorContractSetData, libceed), Cint, (CeedTensorContract, Ptr{Cvoid}), contract, data)
end

function CeedTensorContractDestroy(contract)
    ccall((:CeedTensorContractDestroy, libceed), Cint, (Ptr{CeedTensorContract},), contract)
end

function CeedQFunctionRegister(arg1, arg2, arg3, arg4, init)
    ccall((:CeedQFunctionRegister, libceed), Cint, (Cstring, Cstring, CeedInt, CeedQFunctionUser, Ptr{Cvoid}), arg1, arg2, arg3, arg4, init)
end

function CeedQFunctionSetFortranStatus(qf, status)
    ccall((:CeedQFunctionSetFortranStatus, libceed), Cint, (CeedQFunction, Bool), qf, status)
end

function CeedQFunctionGetCeed(qf, ceed)
    ccall((:CeedQFunctionGetCeed, libceed), Cint, (CeedQFunction, Ptr{Ceed}), qf, ceed)
end

function CeedQFunctionGetVectorLength(qf, vlength)
    ccall((:CeedQFunctionGetVectorLength, libceed), Cint, (CeedQFunction, Ptr{CeedInt}), qf, vlength)
end

function CeedQFunctionGetNumArgs(qf, numinputfields, numoutputfields)
    ccall((:CeedQFunctionGetNumArgs, libceed), Cint, (CeedQFunction, Ptr{CeedInt}, Ptr{CeedInt}), qf, numinputfields, numoutputfields)
end

function CeedQFunctionGetSourcePath(qf, source)
    ccall((:CeedQFunctionGetSourcePath, libceed), Cint, (CeedQFunction, Ptr{Cstring}), qf, source)
end

function CeedQFunctionGetUserFunction(qf, f)
    ccall((:CeedQFunctionGetUserFunction, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionUser}), qf, f)
end

function CeedQFunctionGetContext(qf, ctx)
    ccall((:CeedQFunctionGetContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx)
end

function CeedQFunctionGetInnerContext(qf, ctx)
    ccall((:CeedQFunctionGetInnerContext, libceed), Cint, (CeedQFunction, Ptr{CeedQFunctionContext}), qf, ctx)
end

function CeedQFunctionIsIdentity(qf, isidentity)
    ccall((:CeedQFunctionIsIdentity, libceed), Cint, (CeedQFunction, Ptr{Bool}), qf, isidentity)
end

function CeedQFunctionGetData(qf, data)
    ccall((:CeedQFunctionGetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data)
end

function CeedQFunctionSetData(qf, data)
    ccall((:CeedQFunctionSetData, libceed), Cint, (CeedQFunction, Ptr{Cvoid}), qf, data)
end

function CeedQFunctionGetFields(qf, inputfields, outputfields)
    ccall((:CeedQFunctionGetFields, libceed), Cint, (CeedQFunction, Ptr{Ptr{CeedQFunctionField}}, Ptr{Ptr{CeedQFunctionField}}), qf, inputfields, outputfields)
end

function CeedQFunctionFieldGetName(qffield, fieldname)
    ccall((:CeedQFunctionFieldGetName, libceed), Cint, (CeedQFunctionField, Ptr{Cstring}), qffield, fieldname)
end

function CeedQFunctionFieldGetSize(qffield, size)
    ccall((:CeedQFunctionFieldGetSize, libceed), Cint, (CeedQFunctionField, Ptr{CeedInt}), qffield, size)
end

function CeedQFunctionFieldGetEvalMode(qffield, emode)
    ccall((:CeedQFunctionFieldGetEvalMode, libceed), Cint, (CeedQFunctionField, Ptr{CeedEvalMode}), qffield, emode)
end

function CeedQFunctionContextGetCeed(cxt, ceed)
    ccall((:CeedQFunctionContextGetCeed, libceed), Cint, (CeedQFunctionContext, Ptr{Ceed}), cxt, ceed)
end

function CeedQFunctionContextGetState(ctx, state)
    ccall((:CeedQFunctionContextGetState, libceed), Cint, (CeedQFunctionContext, Ptr{UInt64}), ctx, state)
end

function CeedQFunctionContextGetContextSize(ctx, ctxsize)
    ccall((:CeedQFunctionContextGetContextSize, libceed), Cint, (CeedQFunctionContext, Ptr{Csize_t}), ctx, ctxsize)
end

function CeedQFunctionContextGetBackendData(ctx, data)
    ccall((:CeedQFunctionContextGetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
end

function CeedQFunctionContextSetBackendData(ctx, data)
    ccall((:CeedQFunctionContextSetBackendData, libceed), Cint, (CeedQFunctionContext, Ptr{Cvoid}), ctx, data)
end

function CeedOperatorGetCeed(op, ceed)
    ccall((:CeedOperatorGetCeed, libceed), Cint, (CeedOperator, Ptr{Ceed}), op, ceed)
end

function CeedOperatorGetNumElements(op, numelem)
    ccall((:CeedOperatorGetNumElements, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numelem)
end

function CeedOperatorGetNumQuadraturePoints(op, numqpts)
    ccall((:CeedOperatorGetNumQuadraturePoints, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numqpts)
end

function CeedOperatorGetNumArgs(op, numargs)
    ccall((:CeedOperatorGetNumArgs, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numargs)
end

function CeedOperatorIsSetupDone(op, issetupdone)
    ccall((:CeedOperatorIsSetupDone, libceed), Cint, (CeedOperator, Ptr{Bool}), op, issetupdone)
end

function CeedOperatorGetQFunction(op, qf)
    ccall((:CeedOperatorGetQFunction, libceed), Cint, (CeedOperator, Ptr{CeedQFunction}), op, qf)
end

function CeedOperatorIsComposite(op, iscomposite)
    ccall((:CeedOperatorIsComposite, libceed), Cint, (CeedOperator, Ptr{Bool}), op, iscomposite)
end

function CeedOperatorGetNumSub(op, numsub)
    ccall((:CeedOperatorGetNumSub, libceed), Cint, (CeedOperator, Ptr{CeedInt}), op, numsub)
end

function CeedOperatorGetSubList(op, suboperators)
    ccall((:CeedOperatorGetSubList, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperator}}), op, suboperators)
end

function CeedOperatorGetData(op, data)
    ccall((:CeedOperatorGetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data)
end

function CeedOperatorSetData(op, data)
    ccall((:CeedOperatorSetData, libceed), Cint, (CeedOperator, Ptr{Cvoid}), op, data)
end

function CeedOperatorSetSetupDone(op)
    ccall((:CeedOperatorSetSetupDone, libceed), Cint, (CeedOperator,), op)
end

function CeedOperatorGetFields(op, inputfields, outputfields)
    ccall((:CeedOperatorGetFields, libceed), Cint, (CeedOperator, Ptr{Ptr{CeedOperatorField}}, Ptr{Ptr{CeedOperatorField}}), op, inputfields, outputfields)
end

function CeedOperatorFieldGetElemRestriction(opfield, rstr)
    ccall((:CeedOperatorFieldGetElemRestriction, libceed), Cint, (CeedOperatorField, Ptr{CeedElemRestriction}), opfield, rstr)
end

function CeedOperatorFieldGetBasis(opfield, basis)
    ccall((:CeedOperatorFieldGetBasis, libceed), Cint, (CeedOperatorField, Ptr{CeedBasis}), opfield, basis)
end

function CeedOperatorFieldGetVector(opfield, vec)
    ccall((:CeedOperatorFieldGetVector, libceed), Cint, (CeedOperatorField, Ptr{CeedVector}), opfield, vec)
end

function CeedMatrixMultiply(ceed, matA, matB, matC, m, n, kk)
    ccall((:CeedMatrixMultiply, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt, CeedInt), ceed, matA, matB, matC, m, n, kk)
end


text:
module ParallelScattering
using Distributed, SharedArrays
using LinearAlgebra, Interpolations, Wavelets, FFTW
using ContinuousWavelets
using SpecialFunctions, LinearAlgebra
using HDF5, JLD
using Plots

using ScatteringTransform
import ScatteringTransform: Scattered, ScatteredFull, ScatteredOut, depth, eltypes, scatteringTransform, makeTuple, processArgs
export scatteringTransform, stFlux, stParallel, eltypes, depth, ScatteredOut, ScatteredFull

# extra constructor methods for ScatteredFull and ScatteredOut
function ScatteredFull(layers::scatteringTransform{S,1}, X::Array{T,N}; totalScales = [-1 for i = 1:depth(layers)+1], outputSubsample = (-1, -1)) where {T<:Real,N,S}
    if N == 1
        X = reshape(X, (size(X)..., 1))
    end

    n, q, dataSizes, outputSizes, resultingSize =
        calculateSizes(layers, outputSubsample, size(X),
            totalScales = totalScales)
    numInLayer = getQ(layers, n, totalScales; product = false)
    addedLayers = [numInLayer[min(i, 2):i] for i = 1:depth(layers)+1]
    if 1 == N
        zerr = [zeros(T, n[i], addedLayers[i]...) for i = 1:depth(layers)+1]
        output = [zeros(T, resultingSize[i], addedLayers[i]...) for i = 1:depth(layers)+1]
    else
        zerr = [zeros(T, n[i], q[i], size(X)[2:end]...) for i = 1:depth(layers)+1]
        output = [zeros(T, resultingSize[i], addedLayers[i]..., size(X)[2:end]...)
                  for i = 1:depth(layers)+1]
        @info "" [size(x) for x in output]
    end
    zerr[1][:, 1, axes(X)[2:end]...] = copy(X)
    return ScatteredFull{T,N + 1}(depth(layers), 1, zerr, output)
end

function ScatteredOut(layers::ST, X::Array{T,N};
    totalScales = [-1 for i = 1:depth(layers)+1],
    outputSubsample = (-1, -1)) where {ST<:scatteringTransform,T<:Real,N,S}
    if N == 1
        X = reshape(X, (size(X)..., 1))
    end

    n, q, dataSizes, outputSizes, resultingSize =
        calculateSizes(layers, outputSubsample, size(X), totalScales = totalScales)
    addedLayers = getListOfScaleDims(layers, n, totalScales)
    @info addedLayers
    if 1 == N
        output = [zeros(T, resultingSize[i], addedLayers[i]...) for i = 1:depth(layers)+1]
    else
        output = [zeros(T, resultingSize[i], addedLayers[i]..., size(X)[2:end]...)
                  for i = 1:depth(layers)+1]
    end
    @info [size(x) for x in output]
    return ScatteredOut{T,N + 1}(depth(layers), 1, output)
end

@doc """
        st(X::Array{T, N}, layers::stParallel, nonlinear::nl; fullOr::fullType=fullType(),# subsam::Sub = bspline(), thin::Bool=true, outputSubsample::Tuple{Int, Int}=(-1,-1), subsam::Bool=true, totalScales = [-1 for i=1:depth(layers)+1], percentage = .9, fftPlans = -1) where {T <: Real, S <: Union, N, nl <: Function, Sub <: resamplingMethod}
        stParallel(m, Xlength; CWTType=morl, subsampling = 2, outputSize=[1,1], varargs...) -> layers
1D scattering transform using the stParallel layers. you can switch out the nonlinearity as well as the method of subsampling. Finally, the stType is a string. If it is "full", it will produce all paths. If it is "decreasing", it will only keep paths of increasing scale. If it is "collating", you must also include a vector of matrices.
# Arguments
- `nonlinear` : a type of nonlinearity. Should be a function that acts on Complex numbers
- `thin` : determines whether to wrap the output into a format that can be indexed using paths. `thin` cannot.
- `totalScales`, if positive, gives the number of non-averaging wavelets.
- `outputSubsample` is a tuple, with the first number indicating a rate, and the second number giving the minimum allowed size. If some of the entries are less than 1, it has different behaviour:
    + `(<1, x)` : subsample to x elements for each path.
    + `(<1, <1)` : no ssubsampling
    + `(x, <1)` : subsample at a rate of x, with at least one element kept in each path
- `fullOr::fullType=fullType()` : the structure of the transform either
       `fullType()`, `collatingType()` or `decreasingType()`. At the moment,
       only `fullType()` is functional.
- `fftPlans = false` if not `false`, it should be a 2D array of `Future`s, where the first index is the layer, and the second index the core. See `createFFTPlans` if you want to do this.
"""
struct stParallel{T,Dimension,Depth,subsampType,outType} <: scatteringTransform{Dimension,Depth}
    n::Tuple{Vararg{Int,Dimension}} # the dimensions of a single entry
    shears::Array{T,1} # the array of the transforms; the final of these is
    # used only for averaging, so it has length m+1
    subsampling::subsampType # for each layer, the rate of
    # subsampling. There is one of these for layer zero as well, since the
    # output is subsampled, so it should have length m+1
    outputSize::outType # a list of the size of a single output example
    # dimensions in each layer. The first index is layer, the second is
    # dimension (e.g. a 3 layer shattering transform is 3×2) TODO: currently
    # unused for the 1D case
end
function Base.show(io::IO, l::stParallel{T,D,Depth}) where {T,D,Depth}
    print(io, "stParallel{$T,$D} depth $(Depth), input size $(l.n), subsampling rates $(l.subsampling), outputsizes = $(l.outputSize)")
end
function eltypes(f::stParallel)
    return (eltype(f.shears), length(f.n))
end
# the fully specified version
# TODO: use a layered transform parameter to determine if we're returning a st or thinst instead
@doc """
 """
function stParallel(m, Xlength; CWTType = morl, subsampling = 2, outputSize = [1, 1], varargs...)
    CWTType = makeTuple(m + 1, CWTType)
    subsampling = makeTuple(m + 1, subsampling)
    pairedArgs = processArgs(m + 1, varargs)
    shears = map(x -> ContinuousWavelets.wavelet(x[1]; x[2]...), zip(CWTType, pairedArgs))
    stParallel{typeof(shears[1]),1,m,typeof(subsampling),typeof(outputSize)}((Xlength,), shears, subsampling, outputSize)
end

include("subsampling.jl")
export sizes, bsplineType, bilinearType
include("modifiedTransforms.jl")
include("Utils.jl")
export calcuateSizes, calculateThinStSizes, createFFTPlans, remoteMultiply,
    createRemoteFFTPlan, computeAllWavelets, plotAllWavelets, getQ
include("nonlinearities.jl")
export spInverse, aTanh, Tanh, ReLU, piecewiseLinear, plInverse
include("transform.jl")
export st, transformMidLayer!, transformFinalLayer!
include("pathMethods.jl")
export pathToThinIndex, plotCoordinate, numberSkipped, logabs, maxPooling, numScales,
    incrementKeeper, numInLayer
include("postProcessing.jl")
export MatrixAggrigator, reshapeFlattened, loadSyntheticMatFile, transformFolder, flatten, rollSt
end

ext: jl


text:
# ---
# title: 713. Subarray Product Less Than K
# id: problem713
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/subarray-product-less-than-k/description/>
# hidden: true
# ---
#
# Your are given an array of positive integers `nums`.
#
# Count and print the number of (contiguous) subarrays where the product of all
# the elements in the subarray is less than `k`.
#
# **Example 1:**
#
#
#
#     Input: nums = [10, 5, 2, 6], k = 100
#     Output: 8
#     Explanation: The 8 subarrays that have product less than 100 are: [10], [5], [2], [6], [10, 5], [5, 2], [2, 6], [5, 2, 6].
#     Note that [10, 5, 2] is not included as the product of 100 is not strictly less than k.
#
#
# **Note:**
#
# * `0 < nums.length <= 50000`.
# * `0 < nums[i] < 1000`.
# * `0 <= k < 10^6`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module FirstOrder

using LinearAlgebra

import ..LocalDescent:
search,
LineSearch, StrongBacktracking

struct Termination
    max_iter
    ϵ_abs
    ϵ_rel
    ϵ_grad
    Termination(;
        max_iter = 10_000,
        ϵ_abs = eps(),
        ϵ_rel = eps(),
        ϵ_grad = eps()) =
        new(max_iter, ϵ_abs, ϵ_rel, ϵ_grad)
    Termination(ϵ;
        max_iter = 10_000) =
        new(max_iter, ϵ, ϵ, ϵ)
end

mutable struct TerminationConditions
    abs
    rel
    grad
    TerminationConditions(term::Termination) =
    new((fx, fx_next)->fx - fx_next < term.ϵ_abs,
        (fx, fx_next)->fx - fx_next < term.ϵ_abs * abs(fx),
        ∇fx_next->norm(∇fx_next) < term.ϵ_grad)
end

function descent_until(term::Termination, descent_method, f, ∇f, X, trace)
    term_cond = TerminationConditions(term)
    for _ = 1:term.max_iter
        x = X[:, end]
        x_next = descent_method(f, ∇f, x)
        if trace
            X = hcat(X, x_next)
        else
            X = x_next
        end
        fx, fx_next, ∇fx_next = f(x), f(x_next), ∇f(x_next)
        if term_cond.abs(fx, fx_next) ||
            term_cond.rel(fx, fx_next) ||
            term_cond.grad(∇fx_next)
            return X
        end
    end
    @warn "descent_until: max number of iterations reached"
    X
end

abstract type FirstOrderMethods end

# Maximum Gradient Descent

struct MaximumGradientDescent <: FirstOrderMethods end

function search(::MaximumGradientDescent, f, ∇f, x_0; term = Termination(), trace = false)
    descent_until(term, descent_step_maxgd, f, ∇f, x_0, trace)
end

function descent_step_maxgd(f, ∇f, x)
    """
    Maximum Gradient Descent

    A step of maximum gradient descent produces the value of x(k+1)
    where f(x) in the direction of maximum descent is at its local
    minimum.
    """
    g = ∇f(x)
    d = direction_maxgd(g)
    search(LineSearch(), f, x, d)
end

function direction_maxgd(g)
    -g / norm(g)
end

# Gradient Descent

struct GradientDescent <: FirstOrderMethods
    α
    GradientDescent(; α = 0.001) = new(α)
end

function search(params::GradientDescent, f, ∇f, x_0; term = Termination(), trace = false)
    descent_step = (f, ∇f, x)->descent_step_gd(params.α, f, ∇f, x)
    descent_until(term, descent_step, f, ∇f, x_0, trace)
end

function descent_step_gd(α, f, ∇f, x)
    """
    Gradient Descent
    Algorithm 5.1

    A step of a gradient descent algorithm with fixed learning rate α.
    """
    g = ∇f(x)
    x - α * g
end

# Conjugate Gradient Descent

struct ConjugateGradientDescent <: FirstOrderMethods end

mutable struct CGDProblemState
    g
    d
    CGDProblemState(dim) = new(ones(dim), zeros(dim))
end

function search(params::ConjugateGradientDescent, f, ∇f, x_0; term = Termination(), trace = false)
    state = CGDProblemState(size(x_0))
    descent_step = (f, ∇f, x)->descent_step_cgd!(state, f, ∇f, x)
    descent_until(term, descent_step, f, ∇f, x_0, trace)
end

function descent_step_cgd!(state, f, ∇f, x)
    """
    Conjugate Gradient Descent
    Algorithm 5.2

    https://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
    """
    g = ∇f(x)
    d = direction_cgd(g, state.g, state.d)
    state.g = g
    state.d = d
    search(StrongBacktracking(), f, ∇f, x, d)
end

function direction_cgd(g, gm1, dm1)
    """
    Polak-Ribiere update
    Notation:
        g is for g(k)
        dm1 is for d(k-1)
        gm1 indicates g(k-1)
    """
    β_PR = (g' * (g - gm1)) / (gm1'gm1)
    β = max(β_PR, 0)
    -g + β * dm1
end

end # module


text:
import Term: Theme

@testset "\e[34mtheme" begin
    show(devnull, Theme())  # coverage
    @test true
end

# TODO: write more tests


text:
#################### TransformDistribution ####################

const TransformDistribution{T<:ContinuousUnivariateDistribution} =
  Union{T, Truncated{T}}

function link(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    logit((x - a) / (b - a))
  elseif lowerbounded
    log(x - a)
  elseif upperbounded
    log(b - x)
  else
    x
  end
end

function invlink(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    (b - a) * invlogit(x) + a
  elseif lowerbounded
    exp(x) + a
  elseif upperbounded
    b - exp(x)
  else
    x
  end
end

function logpdf(d::TransformDistribution, x::Real, transform::Bool)
  lp = logpdf(d, x)
  if transform
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      lp += log((x - a) * (b - x) / (b - a))
    elseif lowerbounded
      lp += log(x - a)
    elseif upperbounded
      lp += log(b - x)
    end
  end
  lp
end

#################### RealDistribution ####################

const RealDistribution =
          Union{Cauchy, Gumbel, Laplace, Logistic, NoncentralT, Normal,
                NormalCanon, TDist}

link(d::RealDistribution, x::Real) = x

invlink(d::RealDistribution, x::Real) = x

logpdf(d::RealDistribution, x::Real, transform::Bool) = logpdf(d, x)

#################### PositiveDistribution ####################

const PositiveDistribution =
          Union{BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet,
                Gamma, InverseGamma, InverseGaussian, Kolmogorov, LogNormal,
                NoncentralChisq, NoncentralF, Rayleigh, Weibull}

link(d::PositiveDistribution, x::Real) = log(x)

invlink(d::PositiveDistribution, x::Real) = exp(x)

function  logpdf(d::PositiveDistribution, x::Real, transform::Bool)
  lp = logpdf(d, x)
  transform ? lp + log(x) : lp
end

#################### UnitDistribution ####################

const UnitDistribution =
          Union{Beta, KSOneSided, NoncentralBeta}

link(d::UnitDistribution, x::Real) = logit(x)

invlink(d::UnitDistribution, x::Real) = invlogit(x)

function logpdf(d::UnitDistribution, x::Real, transform::Bool)
  lp = logpdf(d, x)
  transform ? lp + log(x * (1.0 - x)) : lp
end

const PhylogeneticDistribution = Union{PhyloDist, MultiplePhyloDist}


text:
##
#
# This is a simple example for travel time calculation in 3D using the Fast Marching approach for the factored eikonal equation on a regular grid.
# See more details in:
# Eran Treister and Eldad Haber, A fast marching algorithm for the factored eikonal equation, Journal of Computational Physics, 324, 210-225, 2016.
#
# Please cite this paper if you use our package.
#
##
using jInv.Mesh
using FactoredEikonalFastMarching

plotting = false;
if plotting
	using PyPlot;
	using jInvVis;
end

#### SETTING THE PARAMETERS OF THE PROBLEM

dim = 3;
n_nodes = 0;
if dim == 2
	n_nodes = [33;33];
	domain = [0.0,1.0,0.0,1.0]; # Domain is in KMs
	src = div(n_nodes,2);
	src[2] = 1;
else
	n_nodes = [33;33;33];
	domain = [0.0,1.0,0.0,1.0,0.0,1.0]; # Domain is in KMs
	src = div(n_nodes,2);
	src[3] = 1;
end
n_tup = tuple(n_nodes...);
mNodal = ones(Float64,n_tup); # m is 1/v^2 whre v is the velocity in km/sec

#### SETTING THE PARAMETERS OF THE FAST MARCHING SOLVER

M = getRegularMesh(domain,n_nodes-1);
mem = getEikonalTempMemory(n_nodes); # here we allocate temporary memory. Reuse this memory for multiple sources on the same mesh.
eikParam = getEikonalParam(M,mNodal,src,true);
solveFastMarchingUpwindGrad(eikParam, mem);
mem = 0; # deallocate memory.
T = copy(eikParam.T1); # here we get the factor tau1.

## Below we extract the true travel time (multiply it wit tau0)
# This function works on the solution vector itself.
# Use getAnalytic2DeikonalSolutionAll and getAnalytic3DeikonalSolutionAll to get the analytical solution itself for a constant velocity model.
if M.dim==2
	selfMultiplyWithAnalyticSolution2D(n_nodes,M.h,src,T);
else
	selfMultiplyWithAnalyticSolution3D(n_nodes,M.h,src,T)
end
# Here, T is the solution of the eikonal equation, the answer is in seconds.

#### PLOTTING THE RESULTS

if plotting
	plotModel(T);
end


text:
const use_stdcall = startswith(basename(libmpi), "msmpi")

# This is required in addition to __init__() so that we can call library at precompilation time.
if Sys.isunix()
    Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
end

macro mpicall(expr)
    @assert expr isa Expr && expr.head == :call && expr.args[1] == :ccall

    # On unix systems we call the global symbols to allow for LD_PRELOAD interception
    # It can be emulated in Windows (via Libdl.dllist), but this is not fast.
    if Sys.isunix() && expr.args[2].head == :tuple &&
            (VERSION ≥ v"1.5-" || expr.args[2].args[1] ≠ :(:MPI_Get_library_version))
        expr.args[2] = expr.args[2].args[1]
    end

    # Microsoft MPI uses stdcall calling convention
    # this only affects 32-bit Windows
    # unfortunately we need to use ccall to call Get_library_version
    # so check using library name instead
    if use_stdcall
        insert!(expr.args, 3, :stdcall)
    end
    return esc(expr)
end

function Get_library_version()
    # There is no way to query at runtime what the length of the buffer should be.
    # https://github.com/mpi-forum/mpi-issues/issues/159
    # 8192 is the maximum value of MPI_MAX_LIBRARY_VERSION_STRING across known
    # implementations.
    buf = Array{UInt8}(undef, 8192)
    buflen = Ref{Cint}()

    libfilename, = split(basename(libmpi),'.')
    @mpicall ccall((:MPI_Get_library_version, libmpi), Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    resize!(buf, buflen[])
    return String(buf)
end

"""
    MPI_LIBRARY_VERSION_STRING :: String

The full version string provided by the library

# External links
$(_doc_external("MPI_Get_library_version"))
"""
const MPI_LIBRARY_VERSION_STRING = Get_library_version()

"""
    MPIImpl

An enum corresponding to known MPI implementations

- `UnknownMPI`: unable to determine MPI implementation
- `MPICH`: [MPICH](https://www.mpich.org/)
- `OpenMPI`: [Open MPI](https://www.open-mpi.org/)
- `MicrosoftMPI`: [Microsoft MPI](https://docs.microsoft.com/en-us/message-passing-interface/microsoft-mpi)
- `IntelMPI`: [Intel MPI](https://software.intel.com/en-us/mpi-library)
- `SpectrimMPI`: [IBM Spectrum MPI](https://www.ibm.com/us-en/marketplace/spectrum-mpi)
- `MVAPICH`: [MVAPICH](http://mvapich.cse.ohio-state.edu/)
- `CrayMPICH`: Part of the Cray Message Passing Toolkit (MPT)

# See also

- [`MPI_LIBRARY`](@ref)
"""
@enum MPIImpl begin
    UnknownMPI
    MPICH
    OpenMPI
    MicrosoftMPI
    IntelMPI
    IBMSpectrumMPI
    MVAPICH
    CrayMPICH
end

"""
    impl, version = identify_implementation()

Attempt to identify the MPI implementation based on
[`MPI_LIBRARY_VERSION_STRING`](@ref). Returns a triple of values:

- `impl`: a value of type [`MPIImpl`](@ref)
- `version`: a `VersionNumber` of the library, or `nothing` if it cannot be determined.

This function is only intended for internal use. Users should use [`MPI_LIBRARY`](@ref),
[`MPI_LIBRARY_VERSION`](@ref).
"""
function identify_implementation()
    impl = UnknownMPI
    version = v"0"

    if startswith(MPI_LIBRARY_VERSION_STRING, "MPICH")
        impl = MPICH
        # "MPICH Version:\t%s\n" /  "MPICH2 Version:\t%s\n"
        if (m = match(r"^MPICH2? Version:\t(\d+.\d+.\d+\w*)\n", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Open MPI")
        # Open MPI / Spectrum MPI
        impl = occursin("IBM Spectrum MPI", MPI_LIBRARY_VERSION_STRING) ? IBMSpectrumMPI : OpenMPI

        if (m = match(r"^Open MPI v(\d+.\d+.\d+\w*)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Microsoft MPI")
        impl = MicrosoftMPI
        # "Microsoft MPI %u.%u.%u.%u%S"
        # ignore last 2 (build numbers)
        if (m = match(r"^Microsoft MPI (\d+.\d+)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "Intel")
        impl = IntelMPI

        # "Intel(R) MPI Library 2019 Update 4 for Linux* OS"
        if (m = match(r"^Intel\(R\) MPI Library (\d+)(?: Update (\d+))?", MPI_LIBRARY_VERSION_STRING)) !== nothing
            if m.captures[2] === nothing
                version = VersionNumber(m.captures[1])
            else
                version = VersionNumber(m.captures[1]*"."*m.captures[2])
            end
        end

    elseif startswith(MPI_LIBRARY_VERSION_STRING, "MVAPICH2")
        impl = MVAPICH
        # "MVAPICH2 Version      :\t%s\n")
        if (m = match(r"^MVAPICH2? Version\s*:\t(\S*)\n", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif occursin("CRAY MPICH", MPI_LIBRARY_VERSION_STRING)
        impl = CrayMPICH
        # "MPI VERSION    : CRAY MPICH version 7.7.10 (ANL base 3.2)\n"
        if (m = match(r"CRAY MPICH version (\d+.\d+.\d+)", MPI_LIBRARY_VERSION_STRING)) !== nothing
            version = VersionNumber(m.captures[1])
        end
    end

    return impl, version
end

const MPI_LIBRARY, MPI_LIBRARY_VERSION = identify_implementation()

"""
    MPI_LIBRARY :: MPIImpl

The current MPI implementation: this is determined by

# See also
- [`MPIImpl`](@ref)
"""
MPI_LIBRARY

"""
    MPI_LIBRARY_VERSION :: VersionNumber

The version of the MPI library
"""
MPI_LIBRARY_VERSION

function Get_version()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    @mpicall ccall((:MPI_Get_version, libmpi), Cint,
                   (Ptr{Cint}, Ptr{Cint}), major, minor)
    VersionNumber(major[], minor[])
end

"""
    MPI_VERSION :: VersionNumber

The supported version of the MPI standard.

# External links
$(_doc_external("MPI_Get_version"))
"""
const MPI_VERSION = Get_version()


text:
using Test
using MATDaemon
using MATDaemon: JLCallOptions, KeyValueContainer, matlabify

using Distributed
using GarishPrint
using MAT
using Pkg

# This environment variable tells the `jlcall.jl` api script where to find the MATDaemon workspace folder
ENV["MATDAEMON_WORKSPACE"] = mktempdir(; prefix = ".jlcall_", cleanup = true)

# Try loading MATLAB
RUN_MATLAB_TESTS = false
try
    @eval using MATLAB: MEngineError, mxcall
    mxcall(:addpath, 0, realpath(joinpath(@__DIR__, "..", "api")))
    global RUN_MATLAB_TESTS = true
catch e
    @warn "`import MATLAB` failed; skipping MATLAB tests" exception=(e, catch_backtrace())
end

#### Julia tests

include("utils.jl")
include("julia_tests.jl")

#### MATLAB tests

if RUN_MATLAB_TESTS
    include("matlab_tests.jl")
end


text:
module FastPow
export @fastpow

include("addition-chain-pow.jl")
include("exp_n.jl")

fastpow(x) = x
function fastpow(expr::Expr)
    if Meta.isexpr(expr, :call, 3) && expr.args[1] == :^
        if expr.args[2] === 1
            return :($exp_1($(expr.args[3])))
        elseif expr.args[2] === -1
            return :($exp_minus1($(expr.args[3])))
        elseif expr.args[2] === 2
            if expr.args[3] isa Integer && expr.args[3] < 0
                return :(exp2($(expr.args[3])))
            else
                return :($exp_2($(expr.args[3])))
            end
        elseif expr.args[2] === 10
            if expr.args[3] isa Integer && expr.args[3] < 0
                return :(exp10($(expr.args[3])))
            else
                return :($exp_10($(expr.args[3])))
            end
        elseif expr.args[3] isa Integer
            return addition_chain_pow(expr.args[2], expr.args[3])
        end
    end
    return Expr(expr.head, fastpow.(expr.args)...)
end

macro fastpow(expr)
    return esc(fastpow(expr))
end

end # module


text:
# function to calc intensity at given x,y coord.
function line_profile!(mid::T, rot_shift::T, conv_blueshift::T, lambdas::AA{T,1},
                       prof::AA{T,1}, wsp::SynthWorkspace{T}) where T<:AF
    # calculate line center given rot. and conv. doppler shift -> λrest * (1 + z)
    λΔD = mid * (one(T) + rot_shift + conv_blueshift/c_ms)

    # synthesize the line profile given bisector and width input data
    line_from_bis!(λΔD, lambdas, prof, wsp.wavt, wsp.dept, wsp.widt,
                   wsp.lwavgrid, wsp.rwavgrid, wsp.allwavs, wsp.allints)
    return nothing
end

function line_from_bis!(mid::T, lambdas::AA{T,1}, prof::AA{T,1},
                        wavm::AA{T,1}, depm::AA{T,1}, widm::AA{T,1},
                        lwavgrid::AA{T,1}, rwavgrid::AA{T,1},
                        allwavs::AA{T,1}, allints::AA{T,1}) where T<:AF
    # set wavgrids to line center to start
    lwavgrid .= (mid .- (0.5 .* widm .- wavm))
    rwavgrid .= (mid .+ (0.5 .* widm .+ wavm))
    rwavgrid[1] = lwavgrid[1] + 1e-3            # TODO: fix to deal with nodes

    # concatenate into one big array
    len = length(rwavgrid)
    itr = len:-1:1
    allwavs[(len+1):end] .= rwavgrid
    allints[(len+1):end] .= depm
    allwavs[1:len] .= view(lwavgrid, itr)
    allints[1:len] .= view(depm, itr)

    # interpolate onto original lambda grid, extrapolate to continuum
    xs = (allwavs,)
    ys = allints
    it = Gridded(Linear())
    itp1 = extrapolate(interpolate!(T, xs, ys, it), 1.0)
    prof .*= itp1.(lambdas)
    return nothing
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
module DataFramesMeta

importall Base
importall DataFrames
using DataFrames

# Basics:
export @with, @ix, @where, @orderby, @transform, @by, @based_on, @select
export where, orderby, transform, select

include("compositedataframe.jl")
include("linqmacro.jl")
include("byrow.jl")

##############################################################################
##
## @with
##
##############################################################################

replace_syms(x, membernames) = x
function replace_syms(e::Expr, membernames)
    if e.head == :call && length(e.args) == 2 && e.args[1] == :^
        return e.args[2]
    elseif e.head == :.     # special case for :a.b
        return Expr(e.head, replace_syms(e.args[1], membernames),
                            typeof(e.args[2]) == Expr && e.args[2].head == :quote ? e.args[2] : replace_syms(e.args[2], membernames))
    elseif e.head != :quote
        return Expr(e.head, (isempty(e.args) ? e.args : map(x -> replace_syms(x, membernames), e.args))...)
    else
        if haskey(membernames, e.args[1])
            return membernames[e.args[1]]
        else
            a = gensym()
            membernames[e.args[1]] = a
            return a
        end
    end
end

function with_helper(d, body)
    membernames = Dict{Symbol, Symbol}()
    body = replace_syms(body, membernames)
    funargs = map(x -> :( getindex($d, $(Meta.quot(x))) ), collect(keys(membernames)))
    funname = gensym()
    return(:( function $funname($(collect(values(membernames))...)) $body end; $funname($(funargs...)) ))
end

macro with(d, body)
    esc(with_helper(d, body))
end

##############################################################################
##
## @ix - row and row/col selector
##
##############################################################################

ix_helper(d, arg) = :( let d = $d; $d[@with($d, $arg),:]; end )
ix_helper(d, arg, moreargs...) = :( let d = $d; getindex(d, @with(d, $arg), $(moreargs...)); end )

macro ix(d, args...)
    esc(ix_helper(d, args...))
end

##############################################################################
##
## @where - select row subsets
##
##############################################################################

where(d::AbstractDataFrame, arg) = d[arg, :]
where(d::AbstractDataFrame, f::Function) = d[f(d), :]
where(g::GroupedDataFrame, f::Function) = g[Bool[f(x) for x in g]]

where_helper(d, arg) = :( where($d, _DF -> @with(_DF, $arg)) )

macro where(d, arg)
    esc(where_helper(d, arg))
end

##############################################################################
##
## select - select columns
##
##############################################################################

select(d::AbstractDataFrame, arg) = d[ arg]

##############################################################################
##
## @orderby
##
##############################################################################

function orderby(d::AbstractDataFrame, args...)
    D = typeof(d)(args...)
    d[sortperm(D), :]
end
orderby(d::AbstractDataFrame, f::Function) = d[sortperm(f(d)), :]
orderby(g::GroupedDataFrame, f::Function) = g[sortperm([f(x) for x in g])]
orderbyconstructor(d::AbstractDataFrame) = (x...) -> DataFrame(Any[x...])
orderbyconstructor(d) = x -> x

# I don't esc just the input because I want _DF to be visible to the user
macro orderby(d, args...)
    esc(:(let _D = $d;  DataFramesMeta.orderby(_D, _DF -> DataFramesMeta.@with(_DF, DataFramesMeta.orderbyconstructor(_D)($(args...)))); end))
end

##############################################################################
##
## transform & @transform
##
##############################################################################

function transform(d::Union(AbstractDataFrame, Associative); kwargs...)
    result = copy(d)
    for (k, v) in kwargs
        result[k] = isa(v, Function) ? v(d) : v
    end
    return result
end

function transform(g::GroupedDataFrame; kwargs...)
    result = DataFrame(g)
    idx2 = cumsum(Int[size(g[i],1) for i in 1:length(g)])
    idx1 = [1; 1 + idx2[1:end-1]]
    for (k, v) in kwargs
        first = v(g[1])
        result[k] = Array(eltype(first), size(result, 1))
        result[idx1[1]:idx2[1], k] = first
        for i in 2:length(g)
            result[idx1[i]:idx2[i], k] = v(g[i])
        end
    end
    return result
end

function transform_helper(x, args...)
    # convert each kw arg value to: _DF -> @with(_DF, arg)
    newargs = [args...]
    for i in 1:length(args)
        newargs[i].args[2] = :( _DF -> @with(_DF, $(newargs[i].args[2]) ) )
    end
    :( transform($x, $(newargs...)) )
end

macro transform(x, args...)
    esc(transform_helper(x, args...))
end

##############################################################################
##
## @based_on - summarize a grouping operation
##
##############################################################################

macro based_on(x, args...)
    esc(:( DataFrames.based_on($x, _DF -> DataFramesMeta.@with(_DF, DataFrames.DataFrame($(args...)))) ))
end

##############################################################################
##
## @by - grouping
##
##############################################################################

macro by(x, what, args...)
    esc(:( DataFrames.by($x, $what, _DF -> DataFramesMeta.@with(_DF, DataFrames.DataFrame($(args...)))) ))
end

##############################################################################
##
## @select - select and transform columns
##
##############################################################################

expandargs(x) = x

function expandargs(e::Expr)
    if e.head == :quote && length(e.args) == 1
        return Expr(:kw, e.args[1], Expr(:quote, e.args[1]))
    else
        return e
    end
end

function expandargs(e::Tuple)
    res = [e...]
    for i in 1:length(res)
        res[i] = expandargs(e[i])
    end
    return res
end

function select(d::Union(AbstractDataFrame, Associative); kwargs...)
    result = typeof(d)()
    for (k, v) in kwargs
        result[k] = v
    end
    return result
end

macro select(x, args...)
    esc(:(let _DF = $x; DataFramesMeta.@with(_DF, select(_DF, $(DataFramesMeta.expandargs(args)...))); end))
end

##############################################################################
##
## Extras for GroupedDataFrames
##
##############################################################################

combnranges(starts, ends) = [[starts[i]:ends[i] for i in 1:length(starts)]...;]

DataFrame(g::GroupedDataFrame) = g.parent[g.idx[combnranges(g.starts, g.ends)], :]

Base.getindex(gd::GroupedDataFrame, I::AbstractArray{Int}) = GroupedDataFrame(gd.parent,
                                                                              gd.cols,
                                                                              gd.idx,
                                                                              gd.starts[I],
                                                                              gd.ends[I])

end # module


text:
push!(LOAD_PATH,"../src/")

using Documenter
using ReversibleJumpMCMC

makedocs(
    sitename = "ReversibleJumpMCMC",
    format = Documenter.HTML(),
    modules = [ReversibleJumpMCMC]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.

deploydocs(;
    repo="github.com/LidkeLab/ReversibleJumpMCMC.jl.git",
    devbranch = "main"
)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Tiago Os��io
#* $Id$
#*-------------------------------------------------------------------
type FirstOrder
	FirstOrder()=begin
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Time Constant",
				:Unit=>"s",
				:Default=>4
			)),
			DanaReal(Dict{Symbol,Any}(
				:Unit=>"1/s"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Unit=>"1/s"
			)),
			DanaReal(),
			DanaReal(Dict{Symbol,Any}(
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"State"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Input signal",
				:PosX=>0,
				:PosY=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Output signal",
				:PosX=>1,
				:PosY=>0.5
			)),
			[
				:(diff(x) = A*x + B*u),
				:(y = C*x + D*u),
			],
			[
				"","",
			],
			[:tau,:A,:B,:C,:D,],
			[:x,:u,:y,]
		)
	end
	tau::DanaReal
	A::DanaReal
	B::DanaReal
	C::DanaReal
	D::DanaReal
	x::control_signal
	u::control_signal
	y::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export FirstOrder
function setEquationFlow(in::FirstOrder)
	addEquation(1)
	addEquation(2)
end
function atributes(in::FirstOrder,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/PIDIncr"
	drive!(fields,_)
	return fields
end
FirstOrder(_::Dict{Symbol,Any})=begin
	newModel=FirstOrder()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Test, CasADi
import LinearAlgebra: cross, ×, Symmetric
import Suppressor: @capture_out

cas_symbol = [:SX, :MX]

include("test_constructors.jl")
include("test_generic.jl")
include("test_importexport.jl")
include("test_mathfuns.jl")
include("test_mathops.jl")
include("test_numbers.jl")
include("test_types.jl")
include("test_utils.jl")

for i in cas_symbol
    @eval begin
        test_constructors($i)
        test_generic($i)
        test_importexport($i)
        test_mathfuns($i)
        test_mathops($i)
        test_numbers($i)
        test_types($i)
        test_utils($i)
    end
end

## Test examples
@testset "Test first example                                " begin
    x = SX("x")
    y = SX("y")
    α = 1
    b = 100
    f = (α - x)^2 + b*(y - x^2)^2

    nlp = Dict("x" => vcat([x ; y]), "f" => f)
    S = casadi.nlpsol("S", "ipopt", nlp, Dict("ipopt" => Dict("print_level" => 0), "verbose" => false))

    sol = S(x0 = [0, 0]);

    @test sol["x"].toarray()[1] ≈ 0.9999999999999899
    @test sol["x"].toarray()[2] ≈ 0.9999999999999792
end

@testset "Test second example                               " begin
    opti = casadi.Opti();

    x = opti._variable()
    y = opti._variable()

    opti.minimize( (y - x^2)^2 )
    opti._subject_to(x^2 + y^2 == 1)
    opti._subject_to(x + y >= 1)

    opti.solver("ipopt", Dict("verbose" => false), Dict("print_level" => 0));
    sol = opti.solve();

    @test sol.value(x) ≈ 0.7861513776531158
    @test sol.value(y) ≈ 0.6180339888825889
end


text:
"""
    solveModel(; TT=52, ndims=20, thetai=0.38, thetas=0.0044)
    solveModel(nprocs:Int64; TT=52, ndims=20, thetai=0.38, thetas=0.0044)

Solve the model for all possible policy choices and replicate Figures 5, 6, and 7 in Chari, Kirpalani, and Phelan (2021).
The function calls [`nopolicy`](@ref) and [`withpolicy`](@ref) with the relevant keywords, and uses the outputs to generate the figures using
[`createFig5`](@ref), [`createFig6`](@ref), and [`createFig7`](@ref).

# Arguments:
- `nprocs::Int64`: Optional. If specified, uses the module [Distributed](https://docs.julialang.org/en/v1/manual/distributed-computing/) for parallelizing computations among `nprocs` workers.
- `TT::Integer`: Number of periods.
- `ndims::Integer`: Number of grid points in each dimension.
- `thetai::Float`: Probability that infected person sends "infected" signal.
- `thetas::Float`: Probability that non-infected person sends "infected" signal.
"""
function solveModel(; TT=52, ndims=20, thetai=0.38, thetas=0.0044)
    println("Solving the model under no policy, and under policies 'notest', 'untarget', 'target', 'isolate'. This may take a while...")

    # Solve model with standard parameters
    @time noint    = nopolicy(;TT);                                            # no intervention at all
    @time notest   = withpolicy("notest"; TT, ndims, thetai, thetas);         # no testing
    @time untarget = withpolicy("untargettest"; TT, ndims, thetai, thetas);   # untargeted testing
    @time target   = withpolicy("targettest"; TT, ndims, thetai, thetas);     # targeted testing
    @time isolate  = withpolicy("isolate"; TT, ndims, thetai, thetas);        # isolation

    # Figures
    println("Generating figures...")
    createFig5(noint, notest, untarget)
    createFig6(notest, untarget, target)
    createFig7(untarget, target, isolate)
end

##########################################################################################
function solveModel(nprocs::Int64; TT=52, ndims=20, thetai=0.38, thetas=0.0044)
    println("Solving the model under no policy, and under policies 'notest', 'untarget', 'target', 'isolate'. This may take a while...")

    # Tell Julia to use several processes for distributed computing
    addprocs(nprocs)

    @everywhere begin
        # Load packages
        @eval using Pkg; Pkg.activate(".") # Activate environment so that workers "see" packages
        @eval using HammerScalpel

    end

    # Solve model with standard parameters
    @time noint    = nopolicy(;TT);                                                   # no intervention at all
    @time notest   = withpolicy("notest", nprocs; TT, ndims, thetai, thetas);        # no testing
    @time untarget = withpolicy("untargettest", nprocs; TT, ndims, thetai, thetas);  # untargeted testing
    @time target   = withpolicy("targettest", nprocs; TT, ndims, thetai, thetas);    # targeted testing
    @time isolate  = withpolicy("isolate", nprocs; TT, ndims, thetai, thetas);       # isolation

    # Figures
    println("Generating figures...")
    createFig5(noint, notest, untarget)
    createFig6(notest, untarget, target)
    createFig7(untarget, target, isolate)

end


text:
# F = E - TS
# S = -k \sum( .... )
function calc_entropy( Focc::Array{Float64,1}, kT::Float64; is_spinpol=false )
    SMALL = 1.e-10
    Nstates = length(Focc)
    e = 0.0
    if is_spinpol
        for ist = 1:Nstates
            if Focc[ist] > SMALL
                e = e + Focc[ist]*log(Focc[ist])
            else
                e = e + (1.0 - Focc[ist])*log(1.0 - Focc[ist])
            end
        end
    else
        for ist = 1:Nstates
            # spin-degenerate case
            if Focc[ist] > SMALL
                e = e + 0.5*Focc[ist]*log(0.5*Focc[ist])
            else
                e = e + (1.0 - 0.5*Focc[ist])*log(1.0 - 0.5*Focc[ist])
            end
        end
        e = 2*e
    end
    # double negative
    return kT*e
end


text:
"""
    emul!(C::GBArray, A::GBArray, B::GBArray, op = BinaryOps.TIMES; kwargs...)::GBArray

Apply the binary operator `op` elementwise on the set intersection of `A` and `B`. Store or
accumulate the result into C. When `op = BinaryOps.TIMES` this is equivalent to `A .* B`,
however any binary operation may be substituted.

As mentioned the pattern of the result is the set intersection of `A` and `B`. For a set
union equivalent see [`eadd!`](@ref).

# Arguments
- `C::GBArray`: the output vector or matrix.
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `B::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `op::MonoidBinaryOrRig = BinaryOps.TIMES`: the binary operation which is applied such that
    `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in both `A` and `B`.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`
"""
function emul!(
    C::GBVecOrMat,
    A::GBArray,
    B::GBArray,
    op::MonoidBinaryOrRig = BinaryOps.TIMES;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A, in2=B)
    size(C) == size(A) == size(B) || throw(DimensionMismatch())
    op = getoperator(op, optype(A, B))
    accum = getaccum(accum, eltype(C))
    if op isa TypedSemiring
        libgb.GrB_Matrix_eWiseMult_Semiring(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    elseif op isa TypedMonoid
        libgb.GrB_Matrix_eWiseMult_Monoid(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    elseif op isa TypedBinaryOperator
        libgb.GrB_Matrix_eWiseMult_BinaryOp(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    else
        throw(ArgumentError("$op is not a valid monoid binary op or semiring."))
    end
    return C
end

"""
    emul(A::GBArray, B::GBArray, op = BinaryOps.TIMES; kwargs...)::GBMatrix

Apply the binary operator `op` elementwise on the set intersection of `A` and `B`. Store or
accumulate the result into C. When `op = BinaryOps.TIMES` this is equivalent to `A .* B`,
however any binary operation may be substituted.

As mentioned the pattern of the result is the set intersection of `A` and `B`. For a set
union equivalent see [`eadd`](@ref).

# Arguments
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `B::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `op::MonoidBinaryOrRig = BinaryOps.TIMES`: the binary operation which is applied such that
    `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in both `A` and `B`.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`

# Returns
- `GBArray`: Output `GBVector` or `GBMatrix` whose eltype is determined by the `eltype` of
    `A` and `B` or the binary operation if a type specific operation is provided.
"""
function emul(
    A::GBArray,
    B::GBArray,
    op::MonoidBinaryOrRig = BinaryOps.TIMES;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    t = inferoutputtype(A, B, op)
    if A isa GBVector && B isa GBVector
        C = GBVector{t}(size(A))
    else
        C = GBMatrix{t}(size(A))
    end
    return emul!(C, A, B, op; mask, accum, desc)
end

"""
    eadd!(C::GBArray, A::GBArray, B::GBArray, op = BinaryOps.PLUS; kwargs...)::GBArray

Apply the binary operator `op` elementwise on the set union of `A` and `B`. Store or
accumulate the result into C. When `op = BinaryOps.PLUS` this is equivalent to `A .+ B`,
however any binary operation may be substituted.

As mentioned the pattern of the result is the set union of `A` and `B`. For a set
intersection equivalent see [`emul!`](@ref).

# Arguments
- `C::GBArray`: the output vector or matrix.
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `B::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `op::MonoidBinaryOrRig = BinaryOps.PLUS`: the binary operation which is applied such that
    `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` and/or `B`.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`
"""
function eadd!(
    C::GBVecOrMat,
    A::GBArray,
    B::GBArray,
    op::MonoidBinaryOrRig = BinaryOps.PLUS;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A, in2 = B)
    size(C) == size(A) == size(B) || throw(DimensionMismatch())
    op = getoperator(op, optype(A, B))
    accum = getaccum(accum, eltype(C))
    if op isa TypedSemiring
        libgb.GrB_Matrix_eWiseAdd_Semiring(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    elseif op isa TypedMonoid
        libgb.GrB_Matrix_eWiseAdd_Monoid(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    elseif op isa TypedBinaryOperator
        libgb.GrB_Matrix_eWiseAdd_BinaryOp(C, mask, accum, op, parent(A), parent(B), desc)
        return C
    else
        throw(ArgumentError("$op is not a valid monoid binary op or semiring."))
    end
    return C
end

"""
    eadd(A::GBArray, B::GBArray, op = BinaryOps.PLUS; kwargs...)::GBArray

Apply the binary operator `op` elementwise on the set union of `A` and `B`. Store or
accumulate the result into C. When `op = BinaryOps.TIMES` this is equivalent to `A .* B`,
however any binary operation may be substituted.

As mentioned the pattern of the result is the set union of `A` and `B`. For a set
intersection equivalent see [`emul`](@ref).

# Arguments
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `B::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `op::MonoidBinaryOrRig = BinaryOps.PLUS`: the binary operation which is applied such that
    `C[i,j] = op(A[i,j], B[i,j])` for all `i,j` present in either `A` or `B`.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`

# Returns
- `GBArray`: Output `GBVector` or `GBMatrix` whose eltype is determined by the `eltype` of
    `A` and `B` or the binary operation if a type specific operation is provided.
"""
function eadd(
    A::GBArray,
    B::GBArray,
    op::MonoidBinaryOrRig = BinaryOps.PLUS;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    t = inferoutputtype(A, B, op)
    if A isa GBVector && B isa GBVector
        C = GBVector{t}(size(A))
    else
        C = GBMatrix{t}(size(A))
    end
    return eadd!(C, A, B, op; mask, accum, desc)
end

function emul!(C, A, B, op::Function; mask = nothing, accum = nothing, desc = nothing)
    emul!(C, A, B, BinaryOp(op); mask, accum, desc)
end

function emul(A, B, op::Function; mask = nothing, accum = nothing, desc = nothing)
    emul(A, B, BinaryOp(op); mask, accum, desc)
end

function eadd!(C, A, B, op::Function; mask = nothing, accum = nothing, desc = nothing)
    eadd!(C, A, B, BinaryOp(op); mask, accum, desc)
end

function eadd(A, B, op::Function; mask = nothing, accum = nothing, desc = nothing)
    eadd(A, B, BinaryOp(op); mask, accum, desc)
end

function Base.:+(A::GBArray, B::GBArray)
    eadd(A, B, BinaryOps.PLUS)
end

function Base.:-(A::GBArray, B::GBArray)
    eadd(A, B, BinaryOps.MINUS)
end

⊕(A, B, op; mask = nothing, accum = nothing, desc = nothing) =
    eadd(A, B, op; mask, accum, desc)
⊗(A, B, op; mask = nothing, accum = nothing, desc = nothing) =
    emul(A, B, op; mask, accum, desc)

⊕(f::Union{Function, BinaryUnion}) = (A, B; mask = nothing, accum = nothing, desc = nothing) ->
    eadd(A, B, f; mask, accum, desc)

⊗(f::Union{Function, BinaryUnion}) = (A, B; mask = nothing, accum = nothing, desc = nothing) ->
    emul(A, B, f; mask, accum, desc)


text:
abstract type OrdinaryDiffEqInterpolation{cacheType} <: DiffEqBase.AbstractDiffEqInterpolation end

struct InterpolationData{F,uType,tType,kType,cacheType} <: OrdinaryDiffEqInterpolation{cacheType}
  f::F
  timeseries::uType
  ts::tType
  ks::kType
  dense::Bool
  cache::cacheType
end

struct CompositeInterpolationData{F,uType,tType,kType,cacheType} <: OrdinaryDiffEqInterpolation{cacheType}
  f::F
  timeseries::uType
  ts::tType
  ks::kType
  alg_choice::Vector{Int}
  dense::Bool
  cache::cacheType
end

DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where {cacheType<:FunctionMapConstantCache} = "left-endpoint piecewise constant"
DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where {cacheType<:FunctionMapCache} = "left-endpoint piecewise constant"
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{DP5ConstantCache,DP5Cache}
  interp.dense ? "specialized 4th order \"free\" interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Rosenbrock23ConstantCache,Rosenbrock32ConstantCache,Rosenbrock23Cache,Rosenbrock32Cache}
  interp.dense ? "specialized 2nd order \"free\" stiffness-aware interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Rodas4ConstantCache,Rodas4Cache}
  interp.dense ? "specialized 3rd order \"free\" stiffness-aware interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{SSPRK22,SSPRK22ConstantCache,SSPRK33,SSPRK33ConstantCache,SSPRK432,SSPRK432ConstantCache}
  interp.dense ? "2nd order \"free\" SSP interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{OwrenZen3Cache,OwrenZen3ConstantCache}
  interp.dense ? "specialized 3rd order \"free\" interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{OwrenZen4Cache,OwrenZen4ConstantCache}
  interp.dense ? "specialized 4th order \"free\" interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{OwrenZen5Cache,OwrenZen5ConstantCache}
  interp.dense ? "specialized 5th order \"free\" interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Tsit5Cache,Tsit5ConstantCache}
  interp.dense ? "specialized 4th order \"free\" interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{BS5ConstantCache,BS5Cache}
  interp.dense ? "specialized 5th order lazy interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Vern6Cache,Vern6ConstantCache}
  interp.dense ? "specialized 6th order lazy interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Vern7Cache,Vern7ConstantCache}
  interp.dense ? "specialized 7th order lazy interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Vern8Cache,Vern8ConstantCache}
  interp.dense ? "specialized 8th order lazy interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{Vern9Cache,Vern9ConstantCache}
  interp.dense ? "specialized 9th order lazy interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:Union{DP8ConstantCache,DP8Cache}
  interp.dense ? "specialized 7th order interpolation" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType
  interp.dense ? "3rd order Hermite" : "1st order linear"
end
function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{cacheType}) where cacheType<:CompositeCache
  interp.dense ? "Automatic order switching interpolation" : "1st order linear"
end

(interp::InterpolationData)(tvals,idxs,deriv,p,continuity::Symbol=:left) = ode_interpolation(tvals,interp,idxs,deriv,p,continuity)
(interp::CompositeInterpolationData)(tvals,idxs,deriv,p,continuity::Symbol=:left) = ode_interpolation(tvals,interp,idxs,deriv,p,continuity)
(interp::InterpolationData)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = ode_interpolation!(val,tvals,interp,idxs,deriv,p,continuity)
(interp::CompositeInterpolationData)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = ode_interpolation!(val,tvals,interp,idxs,deriv,p,continuity)

function InterpolationData(id::InterpolationData,f)
  InterpolationData(f,id.timeseries,
                      id.ts,
                      id.ks,
                      id.dense,
                      id.cache)
end

function CompositeInterpolationData(id::CompositeInterpolationData,f)
  CompositeInterpolationData(f,id.timeseries,
                               id.ts,
                               id.ks,
                               id.alg_choice,
                               id.dense,
                               id.cache)
end


text:
import Distributions, DataFrames, GLM, StatsModels

function standardize(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})
    std_data = standardize(model_data(model))
    std_model = lm(StatsModels.Formula(model.mf.terms), std_data)
    return std_model
end

function model_data(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})
    return model.mf.df
end

function model_description(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})

    modelname = "linear model"
    outcome = model.mf.terms.eterms[1]
    predictors = model.mf.terms.eterms[2:end]
    formula = StatsModels.Formula(model.mf.terms)

    output = Dict(
        "Model" => modelname,
        "Outcome" => outcome,
        "Predictors" => predictors,
        "Formula" => formula,
        "text_description" => "We fitted a linear regression to predict $outcome with $(join(predictors, ", ", " and ")) ($formula)."
        )
    return Report(text=output["text_description"], values=output)
end

function model_effects_existence(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})
    # p value formula
    # TODO: Damn it's complicated for now!

    coef_table = GLM.coeftable(model)
    p = [x.v for x in coef_table.cols[4]]

    output = Dict(
        "p" => p,
        "p_interpretation" => interpret_p.(p),
        "p_formatted" => format_p.(p))
    return output
end

function model_performance(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})

    R2 = GLM.r2(model)
    R2_adj = GLM.adjr2(model)

    output = Dict(
        "R2" => R2,
        "R2_adj" => R2_adj,
        "text_performance" => "The model's explanatory power (R²) is of $(round(R2, digits=2)) (adj. R² = $(round(R2_adj, digits=2))).")

    return Report(text=output["text_performance"], values=output)
end

function model_std_parameters(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel}; CI=95)

    std_model = standardize(model)

    # Gather values
    coefs = GLM.coef(std_model)
    std_errors = GLM.stderror(std_model)
    ci = GLM.confint(std_model, CI/100)
    ci_lower = ci[:, 1]
    ci_higher = ci[:, 2]

    parameters = Dict{Any, Any}(
        "Std_Coef" => coefs,
        "Std_SE" => std_errors,
        "Std_CI" => ci,
        "Std_CI_lower" => ci_lower,
        "Std_CI_higher" => ci_higher)

    return parameters
end

function model_effect_size(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel}, coefs::Vector{<:Number}; rules="cohen1988", add_std::Bool=true)

    effect_size = interpret_d.(coefs, rules=rules)

    output = Dict{Any, Any}("Effect_Size" => effect_size)

    # Effects
    # --------
    # Generate empty vector
    output["text_effect_size"] = string.(zeros(length(output["Effect_Size"])))
    for i in 1:length(output["text_effect_size"])
        effect = " and can be considered as $(output["Effect_Size"][i])"

        if add_std == true
            effect = effect * " (Std. Coef = $(round(coefs[i], digits=2)))."
        else
            effect = effect * "."
        end

        output["text_effect_size"][i] = effect
    end

    return output
end

function model_initial_parameters(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel})

    intercept = GLM.coef(model)[1]

    output = Dict(
        "Intercept" => intercept,
        "text_initial" => "The model's intercept is at $(round(intercept, digits=2)).")
    return Report(text=output["text_initial"], values=output)
end

function model_parameters(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel};
    CI::Number=95,
    std_coefs::Bool=true,
    effect_size="cohen1988")

    # Gather everything
    # ------------------

    parameters = GLM.coefnames(model)
    coefs = GLM.coef(model)
    std_errors = GLM.stderror(model)
    t = GLM.coeftable(model).cols[3]
    DoF = repeat([GLM.dof_residual(model)], length(parameters))
    ci = GLM.confint(model, CI/100)
    ci_bounds = [(100-CI)/2, 100-(100-CI)/2]
    ci_lower = ci[:, 1]
    ci_higher = ci[:, 2]
    LogLikelihood = GLM.loglikelihood(model)
    deviance = GLM.deviance(model)
    sigma = sqrt(deviance/GLM.dof_residual(model))

    parameters = Dict{Any, Any}(
                "Parameter" => parameters,
                "Coef" => coefs,
                "SE" => std_errors,
                "t" => t,
                "DoF" => DoF,
                "CI_level" => CI,
                "CI" => ci,
                "CI_bounds" => ci_bounds,
                "CI_lower" => ci_lower,
                "CI_higher" => ci_higher,
                "CI_lower_formatted" => "CI_$(ci_bounds[1])",
                "CI_higher_formatted" => "CI_$(ci_bounds[2])",
                "LogLikelihood" => LogLikelihood,
                "Deviance" => deviance)

    parameters = merge(parameters, model_description(model).values)
    parameters = merge(parameters, model_performance(model).values)
    parameters = merge(parameters, model_initial_parameters(model).values)
    parameters = merge(parameters, model_effects_existence(model))
    parameters = merge(parameters, model_std_parameters(model, CI=CI))

    if effect_size != nothing
        if std_coefs == true
            parameters = merge(parameters, model_effect_size(model, parameters["Std_Coef"], rules=effect_size, add_std=true))
        else
            parameters = merge(parameters, model_effect_size(model, coefs, rules=effect_size, add_std=false))
        end
    end

    # Effects
    # --------

    # Generate empty vector
    parameters["text_parameters"] = string.(zeros(length(parameters["Parameter"])))
    for (i, var) in enumerate(parameters["Parameter"])
        effect =
        "$var is $(parameters["p_interpretation"][i]) " *
        "(Coef = $(round(parameters["Coef"][i], digits=2)), " *
        "t($(Int(parameters["DoF"][i]))) = $(round(parameters["t"][i], digits=2)), " *
        "$(parameters["CI_level"])% CI "*
        "[$(round(parameters["CI_lower"][i], digits=2)); " *
        "$(round(parameters["CI_higher"][i], digits=2))]" *
        ", $(parameters["p_formatted"][i]))"
        if effect_size != nothing
            effect = effect * parameters["text_effect_size"][i]
        else
            effect = effect * "."
        end

        parameters["text_parameters"][i] = effect
    end

    return parameters
end

"""
    report(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel};
        CI::Number=95,
        std_coefs::Bool=true,
        effect_size="cohen1988")

Describe a linear model.

# Arguments
- `model`: A `LinearModel`.
- `CI`: Confidence interval level.
- `std_coefs`: Interpret effect sizes of standardized, rather than raw, coefs.
- `effect_size`: Can be 'cohen1988', 'sawilowsky2009', or custom set of [Rules](@ref) (See [cohen_d](@ref)). Set to `nothing` to omit interpretation.

# Examples
```julia
using GLM, DataFrames

model = lm(@formula(y ~ Var1), DataFrame(y=[0, 1, 2, 3], Var1=[2, 3, 3.5, 4]))
report(model)
```
"""
function report(model::StatsModels.DataFrameRegressionModel{<:GLM.LinearModel};
    CI::Number=95,
    std_coefs::Bool=true,
    effect_size="cohen1988")

    # Parameters
    parameters = model_parameters(model,
                                CI=CI,
                                std_coefs=std_coefs,
                                effect_size=effect_size)

    # Text
    description = parameters["text_description"]
    performance = parameters["text_performance"]
    initial = parameters["text_initial"]

    text = "$description $performance $initial Within this model:"
    text = text * join("\n  - " .* parameters["text_parameters"][2:end])

    text = format_text(text)

    # Table
    table = hcat(parameters["Parameter"],
                parameters["Coef"],
                parameters["SE"],
                parameters["t"],
                parameters["DoF"],
                parameters["CI_lower"],
                parameters["CI_higher"],
                parameters["p"])

    table = DataFrames.DataFrame(table,
                [:Parameter,
                :Coef,
                :SE,
                :t,
                :DoF,
                Symbol(parameters["CI_lower_formatted"]),
                Symbol(parameters["CI_higher_formatted"]),
                :p])
    return Report(text=text, values=parameters, table=table)
end


text:
function see(x)
    if x<0
        throw(DomainError("Too small"))
    elseif x>1000
        throw(DomainError("Too big"))
    end
    println(x)
end

try
    see(8)
    see(-9)
catch i
    if isa(i,DomainError)
        print(i.val)
    end
end

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using StatsFuns: logistic
using DelimitedFiles

function run_normal_cv(n :: Int64, p :: Int64, k :: Int64, debias :: Bool)
    #construct snpmatrix, covariate files, and true model b
    x, maf = simulate_random_snparray(n, p)
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1) # non-genetic covariates, just the intercept
    true_b = zeros(p)
    true_b[1:k] = randn(k)
    shuffle!(true_b)
    correct_position = findall(x -> x != 0, true_b)
    noise = rand(Normal(0, 0.1), n) # noise vectors from N(0, s)

    #simulate phenotypes (e.g. vector y) via: y = Xb + noise
    y = xbm * true_b + noise

    #specify path and folds
    path = collect(1:20)
    num_folds = 5
    folds = rand(1:num_folds, size(x, 1))

    #compute cross validation
    result = @timed cv_iht(x, z, y, 1, path, folds, num_folds, use_maf = false, glm = "normal", debias=debias)
    mses = result[1]
    runtime = result[2] #seconds
    memory = result[3] / 1e6  #MB

    return mses, runtime, memory
end

function run_normal_cv(n :: Int64, p :: Int64)
    #run the real code to
    Random.seed!(2019)
    k = 10 # number of true predictors
    repeats = 30
    mses = zeros(20, repeats)
    run_time = zeros(repeats)
    memory = zeros(repeats)
    debias = true
    for i in 1:repeats
        mse, t, m = run_normal_cv(n, p, k, debias)
        mses[:, i] .= mse
        run_time[i] = t
        memory[i] = m
    end

    writedlm("normal_cv_mses", mses)
    writedlm("normal_cv_run_times", run_time)
    writedlm("normal_cv_memory", memory)
end

run_normal_cv(5000, 50000)


text:
"The NodeView provides a read-only comfortable interface for graph traversal."
struct NodeView{G<:Graphs.SequenceDistanceGraph}
    ws::WorkSpace
    id::Graphs.NodeID
    graph::G
end

"""
The LinkView provides a read-only interface to links in the a SequenceDistanceGraph.

A LinkView can only be generated using the `next` and `prev` methods of the
`NodeView` type.
"""
struct LinkView{G<:Graphs.SequenceDistanceGraph}
    node_view::NodeView{G}
    distance::Int32
end

"""
    Get a view of node `nodeid` in the base SequenceDistanceGraph of the
    WorkSpace `ws`.
"""
function node(ws::WorkSpace, nodeid::Graphs.NodeID)
    return NodeView(ws, nodeid, ws.sdg)
end

node(lv::LinkView) = lv.node_view

"""
    Get the sequence of the underlying SequenceDistanceGraph node.
"""
@inline sequence(nv::NodeView) = sequence(graph(nv), id(nv))
@inline sequence_unsafe(nv::NodeView) = sequence_unsafe(graph(nv), id(nv))

Base.length(nv::NodeView) = length(sequence_unsafe(nv))

id(nv::NodeView) = nv.id
graph(nv::NodeView) = nv.graph
workspace(nv::NodeView) = nv.ws

function next(nv::NodeView{G}) where {G}
    g = graph(nv)
    fwl = forward_links(g, id(nv))
    r = Vector{LinkView{G}}(undef, length(fwl))
    i = 1
    @inbounds for l in fwl
        r[i] = LinkView{G}(NodeView{G}(workspace(nv), destination(l), g), distance(l))
        i = i + 1
    end
    return r
end

function Base.summary(io::IO, nv::NodeView)
    print(io, "A view of a graph node (node: ", id(nv), ", graph: ", name(graph(nv)), "):")
end

function Base.show(io::IO, ::MIME"text/plain", nv::NodeView)
    summary(io, nv)
    print(io, "\n  ")
    show(io, sequence(nv))
end

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs