text:
struct klu_l_symbolic
    symmetry::Cdouble
    est_flops::Cdouble
    lnz::Cdouble
    unz::Cdouble
    Lnz::Ptr{Cdouble}
    n::Clong
    nz::Clong
    P::Ptr{Clong}
    Q::Ptr{Clong}
    R::Ptr{Clong}
    nzoff::Clong
    nblocks::Clong
    maxblock::Clong
    ordering::Clong
    do_btf::Clong
    structural_rank::Clong
end

struct klu_l_numeric
    n::Clong
    nblocks::Clong
    lnz::Clong
    unz::Clong
    max_lnz_block::Clong
    max_unz_block::Clong
    Pnum::Ptr{Clong}
    Pinv::Ptr{Clong}
    Lip::Ptr{Clong}
    Uip::Ptr{Clong}
    Llen::Ptr{Clong}
    Ulen::Ptr{Clong}
    LUbx::Ptr{Ptr{Cvoid}}
    LUsize::Ptr{Csize_t}
    Udiag::Ptr{Cvoid}
    Rs::Ptr{Cdouble}
    worksize::Csize_t
    Work::Ptr{Cvoid}
    Xwork::Ptr{Cvoid}
    Iwork::Ptr{Clong}
    Offp::Ptr{Clong}
    Offi::Ptr{Clong}
    Offx::Ptr{Cvoid}
    nzoff::Clong
end

struct klu_l_common_struct
    tol::Cdouble
    memgrow::Cdouble
    initmem_amd::Cdouble
    initmem::Cdouble
    maxwork::Cdouble
    btf::Clong
    ordering::Clong
    scale::Clong
    user_order::Ptr{Cvoid}
    user_data::Ptr{Cvoid}
    halt_if_singular::Clong
    status::Clong
    nrealloc::Clong
    structural_rank::Clong
    numerical_rank::Clong
    singular_col::Clong
    noffdiag::Clong
    flops::Cdouble
    rcond::Cdouble
    condest::Cdouble
    rgrowth::Cdouble
    work::Cdouble
    memusage::Csize_t
    mempeak::Csize_t
end

const klu_l_common = klu_l_common_struct

# typedef int ( * ARKRhsFn ) ( realtype t , N_Vector y , N_Vector ydot , void * user_data )
const ARKRhsFn = Ptr{Cvoid}

# typedef int ( * ARKRootFn ) ( realtype t , N_Vector y , realtype * gout , void * user_data )
const ARKRootFn = Ptr{Cvoid}

# typedef int ( * ARKEwtFn ) ( N_Vector y , N_Vector ewt , void * user_data )
const ARKEwtFn = Ptr{Cvoid}

# typedef int ( * ARKRwtFn ) ( N_Vector y , N_Vector rwt , void * user_data )
const ARKRwtFn = Ptr{Cvoid}

# typedef void ( * ARKErrHandlerFn ) ( int error_code , const char * module , const char * function , char * msg , void * user_data )
const ARKErrHandlerFn = Ptr{Cvoid}

# typedef int ( * ARKAdaptFn ) ( N_Vector y , realtype t , realtype h1 , realtype h2 , realtype h3 , realtype e1 , realtype e2 , realtype e3 , int q , int p , realtype * hnew , void * user_data )
const ARKAdaptFn = Ptr{Cvoid}

# typedef int ( * ARKExpStabFn ) ( N_Vector y , realtype t , realtype * hstab , void * user_data )
const ARKExpStabFn = Ptr{Cvoid}

# typedef int ( * ARKVecResizeFn ) ( N_Vector y , N_Vector ytemplate , void * user_data )
const ARKVecResizeFn = Ptr{Cvoid}

# typedef int ( * ARKPostProcessFn ) ( realtype t , N_Vector y , void * user_data )
const ARKPostProcessFn = Ptr{Cvoid}

# typedef int ( * ARKPostProcessStepFn ) ( realtype t , N_Vector y , void * user_data )
const ARKPostProcessStepFn = Ptr{Cvoid}

# typedef int ( * ARKStepStagePredictFn ) ( realtype t , N_Vector zpred , void * user_data )
const ARKStepStagePredictFn = Ptr{Cvoid}

const realtype = Cdouble

struct _generic_N_Vector_Ops
    nvgetvectorid::Ptr{Cvoid}
    nvclone::Ptr{Cvoid}
    nvcloneempty::Ptr{Cvoid}
    nvdestroy::Ptr{Cvoid}
    nvspace::Ptr{Cvoid}
    nvgetarraypointer::Ptr{Cvoid}
    nvsetarraypointer::Ptr{Cvoid}
    nvgetcommunicator::Ptr{Cvoid}
    nvgetlength::Ptr{Cvoid}
    nvlinearsum::Ptr{Cvoid}
    nvconst::Ptr{Cvoid}
    nvprod::Ptr{Cvoid}
    nvdiv::Ptr{Cvoid}
    nvscale::Ptr{Cvoid}
    nvabs::Ptr{Cvoid}
    nvinv::Ptr{Cvoid}
    nvaddconst::Ptr{Cvoid}
    nvdotprod::Ptr{Cvoid}
    nvmaxnorm::Ptr{Cvoid}
    nvwrmsnorm::Ptr{Cvoid}
    nvwrmsnormmask::Ptr{Cvoid}
    nvmin::Ptr{Cvoid}
    nvwl2norm::Ptr{Cvoid}
    nvl1norm::Ptr{Cvoid}
    nvcompare::Ptr{Cvoid}
    nvinvtest::Ptr{Cvoid}
    nvconstrmask::Ptr{Cvoid}
    nvminquotient::Ptr{Cvoid}
    nvlinearcombination::Ptr{Cvoid}
    nvscaleaddmulti::Ptr{Cvoid}
    nvdotprodmulti::Ptr{Cvoid}
    nvlinearsumvectorarray::Ptr{Cvoid}
    nvscalevectorarray::Ptr{Cvoid}
    nvconstvectorarray::Ptr{Cvoid}
    nvwrmsnormvectorarray::Ptr{Cvoid}
    nvwrmsnormmaskvectorarray::Ptr{Cvoid}
    nvscaleaddmultivectorarray::Ptr{Cvoid}
    nvlinearcombinationvectorarray::Ptr{Cvoid}
    nvdotprodlocal::Ptr{Cvoid}
    nvmaxnormlocal::Ptr{Cvoid}
    nvminlocal::Ptr{Cvoid}
    nvl1normlocal::Ptr{Cvoid}
    nvinvtestlocal::Ptr{Cvoid}
    nvconstrmasklocal::Ptr{Cvoid}
    nvminquotientlocal::Ptr{Cvoid}
    nvwsqrsumlocal::Ptr{Cvoid}
    nvwsqrsummasklocal::Ptr{Cvoid}
end

const N_Vector_Ops = Ptr{_generic_N_Vector_Ops}

struct _generic_N_Vector
    content::Ptr{Cvoid}
    ops::N_Vector_Ops
end

const N_Vector = Ptr{_generic_N_Vector}

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

const SUNLinearSolver_Ops = Ptr{_generic_SUNLinearSolver_Ops}

struct _generic_SUNLinearSolver
    content::Ptr{Cvoid}
    ops::SUNLinearSolver_Ops
end

const SUNLinearSolver = Ptr{_generic_SUNLinearSolver}

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

const SUNMatrix_Ops = Ptr{_generic_SUNMatrix_Ops}

struct _generic_SUNMatrix
    content::Ptr{Cvoid}
    ops::SUNMatrix_Ops
end

const SUNMatrix = Ptr{_generic_SUNMatrix}

struct _generic_SUNNonlinearSolver_Ops
    gettype::Ptr{Cvoid}
    initialize::Ptr{Cvoid}
    setup::Ptr{Cvoid}
    solve::Ptr{Cvoid}
    free::Ptr{Cvoid}
    setsysfn::Ptr{Cvoid}
    setlsetupfn::Ptr{Cvoid}
    setlsolvefn::Ptr{Cvoid}
    setctestfn::Ptr{Cvoid}
    setmaxiters::Ptr{Cvoid}
    getnumiters::Ptr{Cvoid}
    getcuriter::Ptr{Cvoid}
    getnumconvfails::Ptr{Cvoid}
end

const SUNNonlinearSolver_Ops = Ptr{_generic_SUNNonlinearSolver_Ops}

struct _generic_SUNNonlinearSolver
    content::Ptr{Cvoid}
    ops::SUNNonlinearSolver_Ops
end

const SUNNonlinearSolver = Ptr{_generic_SUNNonlinearSolver}

struct ARKodeButcherTableMem
    q::Cint
    p::Cint
    stages::Cint
    A::Ptr{Ptr{realtype}}
    c::Ptr{realtype}
    b::Ptr{realtype}
    d::Ptr{realtype}
end

const ARKodeButcherTable = Ptr{ARKodeButcherTableMem}

# typedef int ( * ARKLsJacFn ) ( realtype t , N_Vector y , N_Vector fy , SUNMatrix Jac , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const ARKLsJacFn = Ptr{Cvoid}

# typedef int ( * ARKLsMassFn ) ( realtype t , SUNMatrix M , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const ARKLsMassFn = Ptr{Cvoid}

# typedef int ( * ARKLsPrecSetupFn ) ( realtype t , N_Vector y , N_Vector fy , booleantype jok , booleantype * jcurPtr , realtype gamma , void * user_data )
const ARKLsPrecSetupFn = Ptr{Cvoid}

# typedef int ( * ARKLsPrecSolveFn ) ( realtype t , N_Vector y , N_Vector fy , N_Vector r , N_Vector z , realtype gamma , realtype delta , int lr , void * user_data )
const ARKLsPrecSolveFn = Ptr{Cvoid}

# typedef int ( * ARKLsMassPrecSetupFn ) ( realtype t , void * user_data )
const ARKLsMassPrecSetupFn = Ptr{Cvoid}

# typedef int ( * ARKLsMassPrecSolveFn ) ( realtype t , N_Vector r , N_Vector z , realtype delta , int lr , void * user_data )
const ARKLsMassPrecSolveFn = Ptr{Cvoid}

# typedef int ( * ARKLsJacTimesSetupFn ) ( realtype t , N_Vector y , N_Vector fy , void * user_data )
const ARKLsJacTimesSetupFn = Ptr{Cvoid}

# typedef int ( * ARKLsJacTimesVecFn ) ( N_Vector v , N_Vector Jv , realtype t , N_Vector y , N_Vector fy , void * user_data , N_Vector tmp )
const ARKLsJacTimesVecFn = Ptr{Cvoid}

# typedef int ( * ARKLsMassTimesSetupFn ) ( realtype t , void * mtimes_data )
const ARKLsMassTimesSetupFn = Ptr{Cvoid}

# typedef int ( * ARKLsMassTimesVecFn ) ( N_Vector v , N_Vector Mv , realtype t , void * mtimes_data )
const ARKLsMassTimesVecFn = Ptr{Cvoid}

# typedef int ( * ARKLsLinSysFn ) ( realtype t , N_Vector y , N_Vector fy , SUNMatrix A , SUNMatrix M , booleantype jok , booleantype * jcur , realtype gamma , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const ARKLsLinSysFn = Ptr{Cvoid}

const sunindextype = Int64

# typedef int ( * ARKLocalFn ) ( sunindextype Nlocal , realtype t , N_Vector y , N_Vector g , void * user_data )
const ARKLocalFn = Ptr{Cvoid}

# typedef int ( * ARKCommFn ) ( sunindextype Nlocal , realtype t , N_Vector y , void * user_data )
const ARKCommFn = Ptr{Cvoid}

@cenum MRISTEP_ID::UInt32 begin
    MRISTEP_ARKSTEP = 0
end

# typedef int ( * MRIStepPreInnerFn ) ( realtype t , N_Vector * f , int nvecs , void * user_data )
const MRIStepPreInnerFn = Ptr{Cvoid}

# typedef int ( * MRIStepPostInnerFn ) ( realtype t , N_Vector y , void * user_data )
const MRIStepPostInnerFn = Ptr{Cvoid}

# typedef int ( * CVRhsFn ) ( realtype t , N_Vector y , N_Vector ydot , void * user_data )
const CVRhsFn = Ptr{Cvoid}

# typedef int ( * CVRootFn ) ( realtype t , N_Vector y , realtype * gout , void * user_data )
const CVRootFn = Ptr{Cvoid}

# typedef int ( * CVEwtFn ) ( N_Vector y , N_Vector ewt , void * user_data )
const CVEwtFn = Ptr{Cvoid}

# typedef void ( * CVErrHandlerFn ) ( int error_code , const char * module , const char * function , char * msg , void * user_data )
const CVErrHandlerFn = Ptr{Cvoid}

# typedef int ( * CVLocalFn ) ( sunindextype Nlocal , realtype t , N_Vector y , N_Vector g , void * user_data )
const CVLocalFn = Ptr{Cvoid}

# typedef int ( * CVCommFn ) ( sunindextype Nlocal , realtype t , N_Vector y , void * user_data )
const CVCommFn = Ptr{Cvoid}

# typedef int ( * CVLsJacFn ) ( realtype t , N_Vector y , N_Vector fy , SUNMatrix Jac , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const CVLsJacFn = Ptr{Cvoid}

const CVDlsJacFn = CVLsJacFn

# typedef int ( * CVLsPrecSetupFn ) ( realtype t , N_Vector y , N_Vector fy , booleantype jok , booleantype * jcurPtr , realtype gamma , void * user_data )
const CVLsPrecSetupFn = Ptr{Cvoid}

# typedef int ( * CVLsPrecSolveFn ) ( realtype t , N_Vector y , N_Vector fy , N_Vector r , N_Vector z , realtype gamma , realtype delta , int lr , void * user_data )
const CVLsPrecSolveFn = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesSetupFn ) ( realtype t , N_Vector y , N_Vector fy , void * user_data )
const CVLsJacTimesSetupFn = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesVecFn ) ( N_Vector v , N_Vector Jv , realtype t , N_Vector y , N_Vector fy , void * user_data , N_Vector tmp )
const CVLsJacTimesVecFn = Ptr{Cvoid}

# typedef int ( * CVLsLinSysFn ) ( realtype t , N_Vector y , N_Vector fy , SUNMatrix A , booleantype jok , booleantype * jcur , realtype gamma , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const CVLsLinSysFn = Ptr{Cvoid}

const CVSpilsPrecSetupFn = CVLsPrecSetupFn

const CVSpilsPrecSolveFn = CVLsPrecSolveFn

const CVSpilsJacTimesSetupFn = CVLsJacTimesSetupFn

const CVSpilsJacTimesVecFn = CVLsJacTimesVecFn

# typedef int ( * CVQuadRhsFn ) ( realtype t , N_Vector y , N_Vector yQdot , void * user_data )
const CVQuadRhsFn = Ptr{Cvoid}

# typedef int ( * CVSensRhsFn ) ( int Ns , realtype t , N_Vector y , N_Vector ydot , N_Vector * yS , N_Vector * ySdot , void * user_data , N_Vector tmp1 , N_Vector tmp2 )
const CVSensRhsFn = Ptr{Cvoid}

# typedef int ( * CVSensRhs1Fn ) ( int Ns , realtype t , N_Vector y , N_Vector ydot , int iS , N_Vector yS , N_Vector ySdot , void * user_data , N_Vector tmp1 , N_Vector tmp2 )
const CVSensRhs1Fn = Ptr{Cvoid}

# typedef int ( * CVQuadSensRhsFn ) ( int Ns , realtype t , N_Vector y , N_Vector * yS , N_Vector yQdot , N_Vector * yQSdot , void * user_data , N_Vector tmp , N_Vector tmpQ )
const CVQuadSensRhsFn = Ptr{Cvoid}

# typedef int ( * CVRhsFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector yBdot , void * user_dataB )
const CVRhsFnB = Ptr{Cvoid}

# typedef int ( * CVRhsFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector yBdot , void * user_dataB )
const CVRhsFnBS = Ptr{Cvoid}

# typedef int ( * CVQuadRhsFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector qBdot , void * user_dataB )
const CVQuadRhsFnB = Ptr{Cvoid}

# typedef int ( * CVQuadRhsFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector qBdot , void * user_dataB )
const CVQuadRhsFnBS = Ptr{Cvoid}

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

# typedef int ( * CVLocalFnB ) ( sunindextype NlocalB , realtype t , N_Vector y , N_Vector yB , N_Vector gB , void * user_dataB )
const CVLocalFnB = Ptr{Cvoid}

# typedef int ( * CVCommFnB ) ( sunindextype NlocalB , realtype t , N_Vector y , N_Vector yB , void * user_dataB )
const CVCommFnB = Ptr{Cvoid}

# typedef int ( * CVLsJacFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector fyB , SUNMatrix JB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const CVLsJacFnB = Ptr{Cvoid}

const CVDlsJacFnB = CVLsJacFnB

# typedef int ( * CVLsJacFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , SUNMatrix JB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const CVLsJacFnBS = Ptr{Cvoid}

const CVDlsJacFnBS = CVLsJacFnBS

# typedef int ( * CVLsPrecSetupFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector fyB , booleantype jokB , booleantype * jcurPtrB , realtype gammaB , void * user_dataB )
const CVLsPrecSetupFnB = Ptr{Cvoid}

# typedef int ( * CVLsPrecSetupFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , booleantype jokB , booleantype * jcurPtrB , realtype gammaB , void * user_dataB )
const CVLsPrecSetupFnBS = Ptr{Cvoid}

# typedef int ( * CVLsPrecSolveFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector fyB , N_Vector rB , N_Vector zB , realtype gammaB , realtype deltaB , int lrB , void * user_dataB )
const CVLsPrecSolveFnB = Ptr{Cvoid}

# typedef int ( * CVLsPrecSolveFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , N_Vector rB , N_Vector zB , realtype gammaB , realtype deltaB , int lrB , void * user_dataB )
const CVLsPrecSolveFnBS = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesSetupFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector fyB , void * jac_dataB )
const CVLsJacTimesSetupFnB = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesSetupFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , void * jac_dataB )
const CVLsJacTimesSetupFnBS = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesVecFnB ) ( N_Vector vB , N_Vector JvB , realtype t , N_Vector y , N_Vector yB , N_Vector fyB , void * jac_dataB , N_Vector tmpB )
const CVLsJacTimesVecFnB = Ptr{Cvoid}

# typedef int ( * CVLsJacTimesVecFnBS ) ( N_Vector vB , N_Vector JvB , realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , void * jac_dataB , N_Vector tmpB )
const CVLsJacTimesVecFnBS = Ptr{Cvoid}

# typedef int ( * CVLsLinSysFnB ) ( realtype t , N_Vector y , N_Vector yB , N_Vector fyB , SUNMatrix AB , booleantype jokB , booleantype * jcurB , realtype gammaB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const CVLsLinSysFnB = Ptr{Cvoid}

# typedef int ( * CVLsLinSysFnBS ) ( realtype t , N_Vector y , N_Vector * yS , N_Vector yB , N_Vector fyB , SUNMatrix AB , booleantype jokB , booleantype * jcurB , realtype gammaB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const CVLsLinSysFnBS = Ptr{Cvoid}

const CVSpilsPrecSetupFnB = CVLsPrecSetupFnB

const CVSpilsPrecSetupFnBS = CVLsPrecSetupFnBS

const CVSpilsPrecSolveFnB = CVLsPrecSolveFnB

const CVSpilsPrecSolveFnBS = CVLsPrecSolveFnBS

const CVSpilsJacTimesSetupFnB = CVLsJacTimesSetupFnB

const CVSpilsJacTimesSetupFnBS = CVLsJacTimesSetupFnBS

const CVSpilsJacTimesVecFnB = CVLsJacTimesVecFnB

const CVSpilsJacTimesVecFnBS = CVLsJacTimesVecFnBS

# typedef int ( * IDAResFn ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector rr , void * user_data )
const IDAResFn = Ptr{Cvoid}

# typedef int ( * IDARootFn ) ( realtype t , N_Vector y , N_Vector yp , realtype * gout , void * user_data )
const IDARootFn = Ptr{Cvoid}

# typedef int ( * IDAEwtFn ) ( N_Vector y , N_Vector ewt , void * user_data )
const IDAEwtFn = Ptr{Cvoid}

# typedef void ( * IDAErrHandlerFn ) ( int error_code , const char * module , const char * function , char * msg , void * user_data )
const IDAErrHandlerFn = Ptr{Cvoid}

# typedef int ( * IDABBDLocalFn ) ( sunindextype Nlocal , realtype tt , N_Vector yy , N_Vector yp , N_Vector gval , void * user_data )
const IDABBDLocalFn = Ptr{Cvoid}

# typedef int ( * IDABBDCommFn ) ( sunindextype Nlocal , realtype tt , N_Vector yy , N_Vector yp , void * user_data )
const IDABBDCommFn = Ptr{Cvoid}

# typedef int ( * IDALsJacFn ) ( realtype t , realtype c_j , N_Vector y , N_Vector yp , N_Vector r , SUNMatrix Jac , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const IDALsJacFn = Ptr{Cvoid}

const IDADlsJacFn = IDALsJacFn

# typedef int ( * IDALsPrecSetupFn ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector rr , realtype c_j , void * user_data )
const IDALsPrecSetupFn = Ptr{Cvoid}

# typedef int ( * IDALsPrecSolveFn ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector rr , N_Vector rvec , N_Vector zvec , realtype c_j , realtype delta , void * user_data )
const IDALsPrecSolveFn = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesSetupFn ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector rr , realtype c_j , void * user_data )
const IDALsJacTimesSetupFn = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesVecFn ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector rr , N_Vector v , N_Vector Jv , realtype c_j , void * user_data , N_Vector tmp1 , N_Vector tmp2 )
const IDALsJacTimesVecFn = Ptr{Cvoid}

const IDASpilsPrecSetupFn = IDALsPrecSetupFn

const IDASpilsPrecSolveFn = IDALsPrecSolveFn

const IDASpilsJacTimesSetupFn = IDALsJacTimesSetupFn

const IDASpilsJacTimesVecFn = IDALsJacTimesVecFn

# typedef int ( * IDAQuadRhsFn ) ( realtype tres , N_Vector yy , N_Vector yp , N_Vector rrQ , void * user_data )
const IDAQuadRhsFn = Ptr{Cvoid}

# typedef int ( * IDASensResFn ) ( int Ns , realtype t , N_Vector yy , N_Vector yp , N_Vector resval , N_Vector * yyS , N_Vector * ypS , N_Vector * resvalS , void * user_data , N_Vector tmp1 , N_Vector tmp2 , N_Vector tmp3 )
const IDASensResFn = Ptr{Cvoid}

# typedef int ( * IDAQuadSensRhsFn ) ( int Ns , realtype t , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector rrQ , N_Vector * rhsvalQS , void * user_data , N_Vector yytmp , N_Vector yptmp , N_Vector tmpQS )
const IDAQuadSensRhsFn = Ptr{Cvoid}

# typedef int ( * IDAResFnB ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , void * user_dataB )
const IDAResFnB = Ptr{Cvoid}

# typedef int ( * IDAResFnBS ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrBS , void * user_dataB )
const IDAResFnBS = Ptr{Cvoid}

# typedef int ( * IDAQuadRhsFnB ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rhsvalBQ , void * user_dataB )
const IDAQuadRhsFnB = Ptr{Cvoid}

# typedef int ( * IDAQuadRhsFnBS ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rhsvalBQS , void * user_dataB )
const IDAQuadRhsFnBS = Ptr{Cvoid}

const N_Vector_S = Ptr{N_Vector}

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

# typedef int ( * IDABBDLocalFnB ) ( sunindextype NlocalB , realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector gvalB , void * user_dataB )
const IDABBDLocalFnB = Ptr{Cvoid}

# typedef int ( * IDABBDCommFnB ) ( sunindextype NlocalB , realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , void * user_dataB )
const IDABBDCommFnB = Ptr{Cvoid}

# typedef int ( * IDALsJacFnB ) ( realtype tt , realtype c_jB , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , SUNMatrix JacB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const IDALsJacFnB = Ptr{Cvoid}

const IDADlsJacFnB = IDALsJacFnB

# typedef int ( * IDALsJacFnBS ) ( realtype tt , realtype c_jB , N_Vector yy , N_Vector yp , N_Vector * yS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrB , SUNMatrix JacB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B , N_Vector tmp3B )
const IDALsJacFnBS = Ptr{Cvoid}

const IDADlsJacFnBS = IDALsJacFnBS

# typedef int ( * IDALsPrecSetupFnB ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , realtype c_jB , void * user_dataB )
const IDALsPrecSetupFnB = Ptr{Cvoid}

# typedef int ( * IDALsPrecSetupFnBS ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrB , realtype c_jB , void * user_dataB )
const IDALsPrecSetupFnBS = Ptr{Cvoid}

# typedef int ( * IDALsPrecSolveFnB ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , N_Vector rvecB , N_Vector zvecB , realtype c_jB , realtype deltaB , void * user_dataB )
const IDALsPrecSolveFnB = Ptr{Cvoid}

# typedef int ( * IDALsPrecSolveFnBS ) ( realtype tt , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrB , N_Vector rvecB , N_Vector zvecB , realtype c_jB , realtype deltaB , void * user_dataB )
const IDALsPrecSolveFnBS = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesSetupFnB ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , realtype c_jB , void * user_dataB )
const IDALsJacTimesSetupFnB = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesSetupFnBS ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrB , realtype c_jB , void * user_dataB )
const IDALsJacTimesSetupFnBS = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesVecFnB ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector yyB , N_Vector ypB , N_Vector rrB , N_Vector vB , N_Vector JvB , realtype c_jB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B )
const IDALsJacTimesVecFnB = Ptr{Cvoid}

# typedef int ( * IDALsJacTimesVecFnBS ) ( realtype t , N_Vector yy , N_Vector yp , N_Vector * yyS , N_Vector * ypS , N_Vector yyB , N_Vector ypB , N_Vector rrB , N_Vector vB , N_Vector JvB , realtype c_jB , void * user_dataB , N_Vector tmp1B , N_Vector tmp2B )
const IDALsJacTimesVecFnBS = Ptr{Cvoid}

const IDASpilsPrecSetupFnB = IDALsPrecSetupFnB

const IDASpilsPrecSetupFnBS = IDALsPrecSetupFnBS

const IDASpilsPrecSolveFnB = IDALsPrecSolveFnB

const IDASpilsPrecSolveFnBS = IDALsPrecSolveFnBS

const IDASpilsJacTimesSetupFnB = IDALsJacTimesSetupFnB

const IDASpilsJacTimesSetupFnBS = IDALsJacTimesSetupFnBS

const IDASpilsJacTimesVecFnB = IDALsJacTimesVecFnB

const IDASpilsJacTimesVecFnBS = IDALsJacTimesVecFnBS

# typedef int ( * KINSysFn ) ( N_Vector uu , N_Vector fval , void * user_data )
const KINSysFn = Ptr{Cvoid}

# typedef void ( * KINErrHandlerFn ) ( int error_code , const char * module , const char * function , char * msg , void * user_data )
const KINErrHandlerFn = Ptr{Cvoid}

# typedef void ( * KINInfoHandlerFn ) ( const char * module , const char * function , char * msg , void * user_data )
const KINInfoHandlerFn = Ptr{Cvoid}

# typedef int ( * KINBBDCommFn ) ( sunindextype Nlocal , N_Vector u , void * user_data )
const KINBBDCommFn = Ptr{Cvoid}

# typedef int ( * KINBBDLocalFn ) ( sunindextype Nlocal , N_Vector uu , N_Vector gval , void * user_data )
const KINBBDLocalFn = Ptr{Cvoid}

# typedef int ( * KINLsJacFn ) ( N_Vector u , N_Vector fu , SUNMatrix J , void * user_data , N_Vector tmp1 , N_Vector tmp2 )
const KINLsJacFn = Ptr{Cvoid}

const KINDlsJacFn = KINLsJacFn

# typedef int ( * KINLsPrecSetupFn ) ( N_Vector uu , N_Vector uscale , N_Vector fval , N_Vector fscale , void * user_data )
const KINLsPrecSetupFn = Ptr{Cvoid}

# typedef int ( * KINLsPrecSolveFn ) ( N_Vector uu , N_Vector uscale , N_Vector fval , N_Vector fscale , N_Vector vv , void * user_data )
const KINLsPrecSolveFn = Ptr{Cvoid}

# typedef int ( * KINLsJacTimesVecFn ) ( N_Vector v , N_Vector Jv , N_Vector uu , booleantype * new_uu , void * J_data )
const KINLsJacTimesVecFn = Ptr{Cvoid}

const KINSpilsPrecSetupFn = KINLsPrecSetupFn

const KINSpilsPrecSolveFn = KINLsPrecSolveFn

const KINSpilsJacTimesVecFn = KINLsJacTimesVecFn

struct _N_VectorContent_ManyVector
    num_subvectors::sunindextype
    global_length::sunindextype
    subvec_array::Ptr{N_Vector}
    own_data::Cint
end

const N_VectorContent_ManyVector = Ptr{_N_VectorContent_ManyVector}

@cenum N_Vector_ID::UInt32 begin
    SUNDIALS_NVEC_SERIAL = 0
    SUNDIALS_NVEC_PARALLEL = 1
    SUNDIALS_NVEC_OPENMP = 2
    SUNDIALS_NVEC_PTHREADS = 3
    SUNDIALS_NVEC_PARHYP = 4
    SUNDIALS_NVEC_PETSC = 5
    SUNDIALS_NVEC_CUDA = 6
    SUNDIALS_NVEC_RAJA = 7
    SUNDIALS_NVEC_OPENMPDEV = 8
    SUNDIALS_NVEC_TRILINOS = 9
    SUNDIALS_NVEC_MANYVECTOR = 10
    SUNDIALS_NVEC_MPIMANYVECTOR = 11
    SUNDIALS_NVEC_MPIPLUSX = 12
    SUNDIALS_NVEC_CUSTOM = 13
end

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

const N_VectorContent_Serial = Ptr{_N_VectorContent_Serial}

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

const DlsMat = Ptr{_DlsMat}

@cenum __JL_Ctag_173::UInt32 begin
    PREC_NONE = 0
    PREC_LEFT = 1
    PREC_RIGHT = 2
    PREC_BOTH = 3
end

@cenum __JL_Ctag_174::UInt32 begin
    MODIFIED_GS = 1
    CLASSICAL_GS = 2
end

# typedef int ( * ATimesFn ) ( void * A_data , N_Vector v , N_Vector z )
const ATimesFn = Ptr{Cvoid}

# typedef int ( * PSetupFn ) ( void * P_data )
const PSetupFn = Ptr{Cvoid}

# typedef int ( * PSolveFn ) ( void * P_data , N_Vector r , N_Vector z , realtype tol , int lr )
const PSolveFn = Ptr{Cvoid}

@cenum SUNLinearSolver_Type::UInt32 begin
    SUNLINEARSOLVER_DIRECT = 0
    SUNLINEARSOLVER_ITERATIVE = 1
    SUNLINEARSOLVER_MATRIX_ITERATIVE = 2
end

@cenum SUNLinearSolver_ID::UInt32 begin
    SUNLINEARSOLVER_BAND = 0
    SUNLINEARSOLVER_DENSE = 1
    SUNLINEARSOLVER_KLU = 2
    SUNLINEARSOLVER_LAPACKBAND = 3
    SUNLINEARSOLVER_LAPACKDENSE = 4
    SUNLINEARSOLVER_PCG = 5
    SUNLINEARSOLVER_SPBCGS = 6
    SUNLINEARSOLVER_SPFGMR = 7
    SUNLINEARSOLVER_SPGMR = 8
    SUNLINEARSOLVER_SPTFQMR = 9
    SUNLINEARSOLVER_SUPERLUDIST = 10
    SUNLINEARSOLVER_SUPERLUMT = 11
    SUNLINEARSOLVER_CUSOLVERSP_BATCHQR = 12
    SUNLINEARSOLVER_CUSTOM = 13
end

@cenum SUNMatrix_ID::UInt32 begin
    SUNMATRIX_DENSE = 0
    SUNMATRIX_BAND = 1
    SUNMATRIX_SPARSE = 2
    SUNMATRIX_SLUNRLOC = 3
    SUNMATRIX_CUSPARSE = 4
    SUNMATRIX_CUSTOM = 5
end

# typedef int ( * SUNNonlinSolSysFn ) ( N_Vector y , N_Vector F , void * mem )
const SUNNonlinSolSysFn = Ptr{Cvoid}

# typedef int ( * SUNNonlinSolLSetupFn ) ( booleantype jbad , booleantype * jcur , void * mem )
const SUNNonlinSolLSetupFn = Ptr{Cvoid}

# typedef int ( * SUNNonlinSolLSolveFn ) ( N_Vector b , void * mem )
const SUNNonlinSolLSolveFn = Ptr{Cvoid}

# typedef int ( * SUNNonlinSolConvTestFn ) ( SUNNonlinearSolver NLS , N_Vector y , N_Vector del , realtype tol , N_Vector ewt , void * mem )
const SUNNonlinSolConvTestFn = Ptr{Cvoid}

@cenum SUNNonlinearSolver_Type::UInt32 begin
    SUNNONLINEARSOLVER_ROOTFIND = 0
    SUNNONLINEARSOLVER_FIXEDPOINT = 1
end

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

const SUNLinearSolverContent_Band = Ptr{_SUNLinearSolverContent_Band}

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

const SUNLinearSolverContent_Dense = Ptr{_SUNLinearSolverContent_Dense}

# typedef sunindextype ( * KLUSolveFn ) ( sun_klu_symbolic * , sun_klu_numeric * , sunindextype , sunindextype , double * , sun_klu_common * )
const KLUSolveFn = Ptr{Cvoid}

struct _SUNLinearSolverContent_KLU
    last_flag::Cint
    first_factorize::Cint
    symbolic::Ptr{klu_l_symbolic}
    numeric::Ptr{klu_l_numeric}
    common::klu_l_common
    klu_solver::KLUSolveFn
end

const SUNLinearSolverContent_KLU = Ptr{_SUNLinearSolverContent_KLU}

struct _SUNLinearSolverContent_LapackBand
    N::sunindextype
    pivots::Ptr{sunindextype}
    last_flag::sunindextype
end

const SUNLinearSolverContent_LapackBand = Ptr{_SUNLinearSolverContent_LapackBand}

struct _SUNLinearSolverContent_LapackDense
    N::sunindextype
    pivots::Ptr{sunindextype}
    last_flag::sunindextype
end

const SUNLinearSolverContent_LapackDense = Ptr{_SUNLinearSolverContent_LapackDense}

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

const SUNLinearSolverContent_PCG = Ptr{_SUNLinearSolverContent_PCG}

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

const SUNLinearSolverContent_SPBCGS = Ptr{_SUNLinearSolverContent_SPBCGS}

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

const SUNLinearSolverContent_SPFGMR = Ptr{_SUNLinearSolverContent_SPFGMR}

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

const SUNLinearSolverContent_SPGMR = Ptr{_SUNLinearSolverContent_SPGMR}

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

const SUNLinearSolverContent_SPTFQMR = Ptr{_SUNLinearSolverContent_SPTFQMR}

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

const SUNMatrixContent_Band = Ptr{_SUNMatrixContent_Band}

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

const SUNMatrixContent_Dense = Ptr{_SUNMatrixContent_Dense}

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

const SUNMatrixContent_Sparse = Ptr{_SUNMatrixContent_Sparse}

struct _SUNNonlinearSolverContent_FixedPoint
    Sys::SUNNonlinSolSysFn
    CTest::SUNNonlinSolConvTestFn
    m::Cint
    imap::Ptr{Cint}
    damping::Cint
    beta::realtype
    R::Ptr{realtype}
    gamma::Ptr{realtype}
    cvals::Ptr{realtype}
    df::Ptr{N_Vector}
    dg::Ptr{N_Vector}
    q::Ptr{N_Vector}
    Xvecs::Ptr{N_Vector}
    yprev::N_Vector
    gy::N_Vector
    fold::N_Vector
    gold::N_Vector
    delta::N_Vector
    curiter::Cint
    maxiters::Cint
    niters::Clong
    nconvfails::Clong
    ctest_data::Ptr{Cvoid}
end

const SUNNonlinearSolverContent_FixedPoint = Ptr{_SUNNonlinearSolverContent_FixedPoint}

struct _SUNNonlinearSolverContent_Newton
    Sys::SUNNonlinSolSysFn
    LSetup::SUNNonlinSolLSetupFn
    LSolve::SUNNonlinSolLSolveFn
    CTest::SUNNonlinSolConvTestFn
    delta::N_Vector
    jcur::Cint
    curiter::Cint
    maxiters::Cint
    niters::Clong
    nconvfails::Clong
    ctest_data::Ptr{Cvoid}
end

const SUNNonlinearSolverContent_Newton = Ptr{_SUNNonlinearSolverContent_Newton}

const ARK_NORMAL = Cint(1)

const ARK_ONE_STEP = Cint(2)

const ARK_INTERP_MAX_DEGREE = Cint(5)

const ARK_INTERP_HERMITE = Cint(0)

const ARK_INTERP_LAGRANGE = Cint(1)

const ARK_SUCCESS = Cint(0)

const ARK_TSTOP_RETURN = Cint(1)

const ARK_ROOT_RETURN = Cint(2)

const ARK_WARNING = Cint(99)

const ARK_TOO_MUCH_WORK = -1

const ARK_TOO_MUCH_ACC = -2

const ARK_ERR_FAILURE = -3

const ARK_CONV_FAILURE = -4

const ARK_LINIT_FAIL = -5

const ARK_LSETUP_FAIL = -6

const ARK_LSOLVE_FAIL = -7

const ARK_RHSFUNC_FAIL = -8

const ARK_FIRST_RHSFUNC_ERR = -9

const ARK_REPTD_RHSFUNC_ERR = -10

const ARK_UNREC_RHSFUNC_ERR = -11

const ARK_RTFUNC_FAIL = -12

const ARK_LFREE_FAIL = -13

const ARK_MASSINIT_FAIL = -14

const ARK_MASSSETUP_FAIL = -15

const ARK_MASSSOLVE_FAIL = -16

const ARK_MASSFREE_FAIL = -17

const ARK_MASSMULT_FAIL = -18

const ARK_CONSTR_FAIL = -19

const ARK_MEM_FAIL = -20

const ARK_MEM_NULL = -21

const ARK_ILL_INPUT = -22

const ARK_NO_MALLOC = -23

const ARK_BAD_K = -24

const ARK_BAD_T = -25

const ARK_BAD_DKY = -26

const ARK_TOO_CLOSE = -27

const ARK_VECTOROP_ERR = -28

const ARK_NLS_INIT_FAIL = -29

const ARK_NLS_SETUP_FAIL = -30

const ARK_NLS_SETUP_RECVR = -31

const ARK_NLS_OP_ERR = -32

const ARK_INNERSTEP_ATTACH_ERR = -33

const ARK_INNERSTEP_FAIL = -34

const ARK_OUTERTOINNER_FAIL = -35

const ARK_INNERTOOUTER_FAIL = -36

const ARK_POSTPROCESS_FAIL = -37

const ARK_POSTPROCESS_STEP_FAIL = -37

const ARK_POSTPROCESS_STAGE_FAIL = -38

const ARK_USER_PREDICT_FAIL = -39

const ARK_INTERP_FAIL = -40

const ARK_UNRECOGNIZED_ERROR = -99

const HEUN_EULER_2_1_2 = Cint(0)

const DEFAULT_ERK_2 = HEUN_EULER_2_1_2

const BOGACKI_SHAMPINE_4_2_3 = Cint(1)

const DEFAULT_ERK_3 = BOGACKI_SHAMPINE_4_2_3

const ZONNEVELD_5_3_4 = Cint(3)

const DEFAULT_ERK_4 = ZONNEVELD_5_3_4

const CASH_KARP_6_4_5 = Cint(6)

const DEFAULT_ERK_5 = CASH_KARP_6_4_5

const VERNER_8_5_6 = Cint(10)

const DEFAULT_ERK_6 = VERNER_8_5_6

const FEHLBERG_13_7_8 = Cint(11)

const DEFAULT_ERK_8 = FEHLBERG_13_7_8

const SDIRK_2_1_2 = Cint(100)

const DEFAULT_DIRK_2 = SDIRK_2_1_2

const ARK324L2SA_DIRK_4_2_3 = Cint(104)

const DEFAULT_DIRK_3 = ARK324L2SA_DIRK_4_2_3

const SDIRK_5_3_4 = Cint(107)

const DEFAULT_DIRK_4 = SDIRK_5_3_4

const ARK548L2SA_DIRK_8_4_5 = Cint(111)

const DEFAULT_DIRK_5 = ARK548L2SA_DIRK_8_4_5

const ARK324L2SA_ERK_4_2_3 = Cint(2)

const DEFAULT_ARK_ETABLE_3 = ARK324L2SA_ERK_4_2_3

const ARK436L2SA_ERK_6_3_4 = Cint(4)

const DEFAULT_ARK_ETABLE_4 = ARK436L2SA_ERK_6_3_4

const ARK548L2SA_ERK_8_4_5 = Cint(9)

const DEFAULT_ARK_ETABLE_5 = ARK548L2SA_ERK_8_4_5

const DEFAULT_ARK_ITABLE_3 = ARK324L2SA_DIRK_4_2_3

const ARK436L2SA_DIRK_6_3_4 = Cint(109)

const DEFAULT_ARK_ITABLE_4 = ARK436L2SA_DIRK_6_3_4

const DEFAULT_ARK_ITABLE_5 = ARK548L2SA_DIRK_8_4_5

const BILLINGTON_3_3_2 = Cint(101)

const TRBDF2_3_3_2 = Cint(102)

const KVAERNO_4_2_3 = Cint(103)

const CASH_5_2_4 = Cint(105)

const CASH_5_3_4 = Cint(106)

const KVAERNO_5_3_4 = Cint(108)

const KVAERNO_7_4_5 = Cint(110)

const ARK437L2SA_DIRK_7_3_4 = Cint(112)

const ARK548L2SAb_DIRK_8_4_5 = Cint(113)

const MIN_DIRK_NUM = Cint(100)

const MAX_DIRK_NUM = Cint(113)

const SAYFY_ABURUB_6_3_4 = Cint(5)

const FEHLBERG_6_4_5 = Cint(7)

const DORMAND_PRINCE_7_4_5 = Cint(8)

const KNOTH_WOLKE_3_3 = Cint(12)

const ARK437L2SA_ERK_7_3_4 = Cint(13)

const ARK548L2SAb_ERK_8_4_5 = Cint(14)

const MIN_ERK_NUM = Cint(0)

const MAX_ERK_NUM = Cint(14)

const ARKLS_SUCCESS = Cint(0)

const ARKLS_MEM_NULL = -1

const ARKLS_LMEM_NULL = -2

const ARKLS_ILL_INPUT = -3

const ARKLS_MEM_FAIL = -4

const ARKLS_PMEM_NULL = -5

const ARKLS_MASSMEM_NULL = -6

const ARKLS_JACFUNC_UNRECVR = -7

const ARKLS_JACFUNC_RECVR = -8

const ARKLS_MASSFUNC_UNRECVR = -9

const ARKLS_MASSFUNC_RECVR = -10

const ARKLS_SUNMAT_FAIL = -11

const ARKLS_SUNLS_FAIL = -12

const DEFAULT_MRI_TABLE_3 = KNOTH_WOLKE_3_3

const CV_ADAMS = Cint(1)

const CV_BDF = Cint(2)

const CV_NORMAL = Cint(1)

const CV_ONE_STEP = Cint(2)

const CV_SUCCESS = Cint(0)

const CV_TSTOP_RETURN = Cint(1)

const CV_ROOT_RETURN = Cint(2)

const CV_WARNING = Cint(99)

const CV_TOO_MUCH_WORK = -1

const CV_TOO_MUCH_ACC = -2

const CV_ERR_FAILURE = -3

const CV_CONV_FAILURE = -4

const CV_LINIT_FAIL = -5

const CV_LSETUP_FAIL = -6

const CV_LSOLVE_FAIL = -7

const CV_RHSFUNC_FAIL = -8

const CV_FIRST_RHSFUNC_ERR = -9

const CV_REPTD_RHSFUNC_ERR = -10

const CV_UNREC_RHSFUNC_ERR = -11

const CV_RTFUNC_FAIL = -12

const CV_NLS_INIT_FAIL = -13

const CV_NLS_SETUP_FAIL = -14

const CV_CONSTR_FAIL = -15

const CV_NLS_FAIL = -16

const CV_MEM_FAIL = -20

const CV_MEM_NULL = -21

const CV_ILL_INPUT = -22

const CV_NO_MALLOC = -23

const CV_BAD_K = -24

const CV_BAD_T = -25

const CV_BAD_DKY = -26

const CV_TOO_CLOSE = -27

const CV_VECTOROP_ERR = -28

const CV_UNRECOGNIZED_ERR = -99

const CVDIAG_SUCCESS = Cint(0)

const CVDIAG_MEM_NULL = -1

const CVDIAG_LMEM_NULL = -2

const CVDIAG_ILL_INPUT = -3

const CVDIAG_MEM_FAIL = -4

const CVDIAG_INV_FAIL = -5

const CVDIAG_RHSFUNC_UNRECVR = -6

const CVDIAG_RHSFUNC_RECVR = -7

const CVLS_SUCCESS = Cint(0)

const CVLS_MEM_NULL = -1

const CVLS_LMEM_NULL = -2

const CVLS_ILL_INPUT = -3

const CVLS_MEM_FAIL = -4

const CVLS_PMEM_NULL = -5

const CVLS_JACFUNC_UNRECVR = -6

const CVLS_JACFUNC_RECVR = -7

const CVLS_SUNMAT_FAIL = -8

const CVLS_SUNLS_FAIL = -9

const CV_SIMULTANEOUS = Cint(1)

const CV_STAGGERED = Cint(2)

const CV_STAGGERED1 = Cint(3)

const CV_CENTERED = Cint(1)

const CV_FORWARD = Cint(2)

const CV_HERMITE = Cint(1)

const CV_POLYNOMIAL = Cint(2)

const CV_NO_QUAD = -30

const CV_QRHSFUNC_FAIL = -31

const CV_FIRST_QRHSFUNC_ERR = -32

const CV_REPTD_QRHSFUNC_ERR = -33

const CV_UNREC_QRHSFUNC_ERR = -34

const CV_NO_SENS = -40

const CV_SRHSFUNC_FAIL = -41

const CV_FIRST_SRHSFUNC_ERR = -42

const CV_REPTD_SRHSFUNC_ERR = -43

const CV_UNREC_SRHSFUNC_ERR = -44

const CV_BAD_IS = -45

const CV_NO_QUADSENS = -50

const CV_QSRHSFUNC_FAIL = -51

const CV_FIRST_QSRHSFUNC_ERR = -52

const CV_REPTD_QSRHSFUNC_ERR = -53

const CV_UNREC_QSRHSFUNC_ERR = -54

const CV_NO_ADJ = -101

const CV_NO_FWD = -102

const CV_NO_BCK = -103

const CV_BAD_TB0 = -104

const CV_REIFWD_FAIL = -105

const CV_FWD_FAIL = -106

const CV_GETY_BADT = -107

const CVDIAG_NO_ADJ = -101

const CVLS_NO_ADJ = -101

const CVLS_LMEMB_NULL = -102

const IDA_NORMAL = Cint(1)

const IDA_ONE_STEP = Cint(2)

const IDA_YA_YDP_INIT = Cint(1)

const IDA_Y_INIT = Cint(2)

const IDA_SUCCESS = Cint(0)

const IDA_TSTOP_RETURN = Cint(1)

const IDA_ROOT_RETURN = Cint(2)

const IDA_WARNING = Cint(99)

const IDA_TOO_MUCH_WORK = -1

const IDA_TOO_MUCH_ACC = -2

const IDA_ERR_FAIL = -3

const IDA_CONV_FAIL = -4

const IDA_LINIT_FAIL = -5

const IDA_LSETUP_FAIL = -6

const IDA_LSOLVE_FAIL = -7

const IDA_RES_FAIL = -8

const IDA_REP_RES_ERR = -9

const IDA_RTFUNC_FAIL = -10

const IDA_CONSTR_FAIL = -11

const IDA_FIRST_RES_FAIL = -12

const IDA_LINESEARCH_FAIL = -13

const IDA_NO_RECOVERY = -14

const IDA_NLS_INIT_FAIL = -15

const IDA_NLS_SETUP_FAIL = -16

const IDA_NLS_FAIL = -17

const IDA_MEM_NULL = -20

const IDA_MEM_FAIL = -21

const IDA_ILL_INPUT = -22

const IDA_NO_MALLOC = -23

const IDA_BAD_EWT = -24

const IDA_BAD_K = -25

const IDA_BAD_T = -26

const IDA_BAD_DKY = -27

const IDA_VECTOROP_ERR = -28

const IDA_UNRECOGNIZED_ERROR = -99

const IDALS_SUCCESS = Cint(0)

const IDALS_MEM_NULL = -1

const IDALS_LMEM_NULL = -2

const IDALS_ILL_INPUT = -3

const IDALS_MEM_FAIL = -4

const IDALS_PMEM_NULL = -5

const IDALS_JACFUNC_UNRECVR = -6

const IDALS_JACFUNC_RECVR = -7

const IDALS_SUNMAT_FAIL = -8

const IDALS_SUNLS_FAIL = -9

const IDA_SIMULTANEOUS = Cint(1)

const IDA_STAGGERED = Cint(2)

const IDA_CENTERED = Cint(1)

const IDA_FORWARD = Cint(2)

const IDA_HERMITE = Cint(1)

const IDA_POLYNOMIAL = Cint(2)

const IDA_NO_QUAD = -30

const IDA_QRHS_FAIL = -31

const IDA_FIRST_QRHS_ERR = -32

const IDA_REP_QRHS_ERR = -33

const IDA_NO_SENS = -40

const IDA_SRES_FAIL = -41

const IDA_REP_SRES_ERR = -42

const IDA_BAD_IS = -43

const IDA_NO_QUADSENS = -50

const IDA_QSRHS_FAIL = -51

const IDA_FIRST_QSRHS_ERR = -52

const IDA_REP_QSRHS_ERR = -53

const IDA_NO_ADJ = -101

const IDA_NO_FWD = -102

const IDA_NO_BCK = -103

const IDA_BAD_TB0 = -104

const IDA_REIFWD_FAIL = -105

const IDA_FWD_FAIL = -106

const IDA_GETY_BADT = -107

const IDALS_NO_ADJ = -101

const IDALS_LMEMB_NULL = -102

const KIN_SUCCESS = Cint(0)

const KIN_INITIAL_GUESS_OK = Cint(1)

const KIN_STEP_LT_STPTOL = Cint(2)

const KIN_WARNING = Cint(99)

const KIN_MEM_NULL = -1

const KIN_ILL_INPUT = -2

const KIN_NO_MALLOC = -3

const KIN_MEM_FAIL = -4

const KIN_LINESEARCH_NONCONV = -5

const KIN_MAXITER_REACHED = -6

const KIN_MXNEWT_5X_EXCEEDED = -7

const KIN_LINESEARCH_BCFAIL = -8

const KIN_LINSOLV_NO_RECOVERY = -9

const KIN_LINIT_FAIL = -10

const KIN_LSETUP_FAIL = -11

const KIN_LSOLVE_FAIL = -12

const KIN_SYSFUNC_FAIL = -13

const KIN_FIRST_SYSFUNC_ERR = -14

const KIN_REPTD_SYSFUNC_ERR = -15

const KIN_VECTOROP_ERR = -16

const KIN_ETACHOICE1 = Cint(1)

const KIN_ETACHOICE2 = Cint(2)

const KIN_ETACONSTANT = Cint(3)

const KIN_NONE = Cint(0)

const KIN_LINESEARCH = Cint(1)

const KIN_PICARD = Cint(2)

const KIN_FP = Cint(3)

const KINBBDPRE_SUCCESS = Cint(0)

const KINBBDPRE_PDATA_NULL = -11

const KINBBDPRE_FUNC_UNRECVR = -12

const KINLS_SUCCESS = Cint(0)

const KINLS_MEM_NULL = -1

const KINLS_LMEM_NULL = -2

const KINLS_ILL_INPUT = -3

const KINLS_MEM_FAIL = -4

const KINLS_PMEM_NULL = -5

const KINLS_JACFUNC_ERR = -6

const KINLS_SUNMAT_FAIL = -7

const KINLS_SUNLS_FAIL = -8

const SUNDIALS_VERSION = "5.2.0"

const SUNDIALS_VERSION_MAJOR = Cint(5)

const SUNDIALS_VERSION_MINOR = Cint(2)

const SUNDIALS_VERSION_PATCH = Cint(0)

const SUNDIALS_VERSION_LABEL = ""

const SUNDIALS_DOUBLE_PRECISION = Cint(1)

const SUNDIALS_INT64_T = Cint(1)

const SUNDIALS_INDEX_TYPE = int64_t

const SUNDIALS_MPI_COMM_F2C = Cint(0)

const SUNDIALS_DENSE = Cint(1)

const SUNDIALS_BAND = Cint(2)

const FCMIX_CVODE = Cint(1)

const FCMIX_IDA = Cint(2)

const FCMIX_KINSOL = Cint(3)

const FCMIX_ARKODE = Cint(4)

const dcopy_f77 = nothing

const dscal_f77 = nothing

const dgemv_f77 = nothing

const dtrsv_f77 = nothing

const dsyrk_f77 = nothing

const dgbtrf_f77 = nothing

const dgbtrs_f77 = nothing

const dgetrf_f77 = nothing

const dgetrs_f77 = nothing

const dgeqp3_f77 = nothing

const dgeqrf_f77 = nothing

const dormqr_f77 = nothing

const dpotrf_f77 = nothing

const dpotrs_f77 = nothing

const scopy_f77 = nothing

const sscal_f77 = nothing

const sgemv_f77 = nothing

const strsv_f77 = nothing

const ssyrk_f77 = nothing

const sgbtrf_f77 = nothing

const sgbtrs_f77 = nothing

const sgetrf_f77 = nothing

const sgetrs_f77 = nothing

const sgeqp3_f77 = nothing

const sgeqrf_f77 = nothing

const sormqr_f77 = nothing

const spotrf_f77 = nothing

const spotrs_f77 = nothing

const SUNLS_SUCCESS = Cint(0)

const SUNLS_MEM_NULL = -801

const SUNLS_ILL_INPUT = -802

const SUNLS_MEM_FAIL = -803

const SUNLS_ATIMES_FAIL_UNREC = -804

const SUNLS_PSET_FAIL_UNREC = -805

const SUNLS_PSOLVE_FAIL_UNREC = -806

const SUNLS_PACKAGE_FAIL_UNREC = -807

const SUNLS_GS_FAIL = -808

const SUNLS_QRSOL_FAIL = -809

const SUNLS_VECTOROP_ERR = -810

const SUNLS_RES_REDUCED = Cint(801)

const SUNLS_CONV_FAIL = Cint(802)

const SUNLS_ATIMES_FAIL_REC = Cint(803)

const SUNLS_PSET_FAIL_REC = Cint(804)

const SUNLS_PSOLVE_FAIL_REC = Cint(805)

const SUNLS_PACKAGE_FAIL_REC = Cint(806)

const SUNLS_QRFACT_FAIL = Cint(807)

const SUNLS_LUFACT_FAIL = Cint(808)

const SUNMAT_SUCCESS = Cint(0)

const SUNMAT_ILL_INPUT = -701

const SUNMAT_MEM_FAIL = -702

const SUNMAT_OPERATION_FAIL = -703

const SUNMAT_MATVEC_SETUP_REQUIRED = -704

const MPI_SUNREALTYPE = MPI_DOUBLE

const MPI_SUNINDEXTYPE = MPI_INT64_T

const PVEC_REAL_MPI_TYPE = MPI_SUNREALTYPE

const PVEC_INTEGER_MPI_TYPE = MPI_SUNINDEXTYPE

const SUN_NLS_SUCCESS = Cint(0)

const SUN_NLS_CONTINUE = +901

const SUN_NLS_CONV_RECVR = +902

const SUN_NLS_MEM_NULL = -901

const SUN_NLS_MEM_FAIL = -902

const SUN_NLS_ILL_INPUT = -903

const SUN_NLS_VECTOROP_ERR = -904

const SUN_NLS_EXT_FAIL = -905

# Skipping MacroDefinition: _SUNDIALS_STRUCT_ struct

const BIG_REAL = DBL_MAX

const SMALL_REAL = DBL_MIN

const UNIT_ROUNDOFF = DBL_EPSILON

const booleantype = Cint

const SUNFALSE = Cint(0)

const SUNTRUE = Cint(1)

const SUNKLU_ORDERING_DEFAULT = Cint(1)

const SUNKLU_REINIT_FULL = Cint(1)

const SUNKLU_REINIT_PARTIAL = Cint(2)

const sun_klu_symbolic = nothing

const sun_klu_numeric = nothing

const sun_klu_common = nothing

const sun_klu_analyze = nothing

const sun_klu_factor = nothing

const sun_klu_refactor = nothing

const sun_klu_rcond = nothing

const sun_klu_condest = nothing

const sun_klu_defaults = nothing

const sun_klu_free_symbolic = nothing

const sun_klu_free_numeric = nothing

const xgbtrf_f77 = dgbtrf_f77

const xgbtrs_f77 = dgbtrs_f77

const xgetrf_f77 = dgetrf_f77

const xgetrs_f77 = dgetrs_f77

const SUNPCG_MAXL_DEFAULT = Cint(5)

const SUNSPBCGS_MAXL_DEFAULT = Cint(5)

const SUNSPFGMR_MAXL_DEFAULT = Cint(5)

const SUNSPFGMR_MAXRS_DEFAULT = Cint(0)

const SUNSPFGMR_GSTYPE_DEFAULT = MODIFIED_GS

const SUNSPGMR_MAXL_DEFAULT = Cint(5)

const SUNSPGMR_MAXRS_DEFAULT = Cint(0)

const SUNSPGMR_GSTYPE_DEFAULT = MODIFIED_GS

const SUNSPTFQMR_MAXL_DEFAULT = Cint(5)

const CSC_MAT = Cint(0)

const CSR_MAT = Cint(1)


text:
using LinearAlgebra

include("Fatoracao_QR.jl")

function algoritmo_QR(E::Matrix, tol = 1e-10)
    # E deve ser uma matriz quadrada
    M = deepcopy(E)
    m, n = size(M)
    if m != n
        error("A entrada deve ser uma matriz quadrada")
    end
    X = I
    iter = 0
    diagonal = false
    while diagonal == false
        if iter == 1000
            break
        end
        # verificamos se M é diagonal
        diagonal = true
        for i in 1:m
            for j in 1:n
                if i != j && abs(M[i,j]) > tol
                    diagonal = false
                end
            end
        end
        # se M não é diagonal, então
        # I) calculamos a fatoração QR de M
        Q, R = fatoracao_QR(M)
        # II) M = R * Q
        M = R * Q
        # multiplicando para obter a matriz X de autovetores no fim do processo
        X = X * Q
        iter = iter + 1
    end
    return E, M, X, iter
end

ext: jl

max_stars_repo_name: EleonoraAvello/GrupoBetaPacote

max_issues_repo_name: EleonoraAvello/GrupoBetaPacote

max_forks_repo_name: EleonoraAvello/GrupoBetaPacote


text:
@testset "Primality check (zerodim subroutine)" begin

    R, (x, y) = Singular.PolynomialRing(Singular.QQ, ["x", "y"])

    @test check_primality_zerodim(Singular.Ideal(R, x^2 - 1, y^2 - 4)) == false

    @test check_primality_zerodim(Singular.Ideal(R, (x + 5) * (x^3 - 7), y - 3)) == false

    @test check_primality_zerodim(Singular.Ideal(R, x^3 - 5, y - 1)) == true

    @test check_primality_zerodim(Singular.Ideal(R, x^2 + 1, y^3 - 3 * x + x + 5)) == true

end


text:
import Base: show

mutable struct SampleModel <: CmdStanModels
    name::AbstractString;              # Name of the Stan program
    model::AbstractString;             # Stan language model program
    num_threads::Int64;                # Number of C++ threads

    check_num_chains::Bool;            # Enforce either C++ or Julia chains
    use_cpp_chains::Bool;              # Enable C++ threads and chains
    num_cpp_chains::Int64;             # Number of C++ chains in each exec process
    num_julia_chains::Int64;           # Number of julia chains ( == processes)
    num_chains::Int64;                 # Actual number of chains

    # Sample fields
    num_samples::Int;                  # Number of draws after warmup
    num_warmups::Int;                  # Number of warmup draws
    save_warmup::Bool;                 # Store warmup_samples
    thin::Int;                         # Thinning of draws
    seed::Int;                         # Seed section of cmd to run cmdstan
    refresh::Int                       # Display progress in output files
    init_bound::Int                    # Bound for initial param values

    # Adapt fields
    engaged::Bool;                     # Adaptation enganged?.
    gamma::Float64;                    # Adaptation regularization scale
    delta::Float64;                    # Adaptation target acceptance statistic
    kappa::Float64;                    # Adaptation relaxation exponent
    t0::Int;                           # Adaptation iteration offset
    init_buffer::Int;                  # Width initial adaptation interval
    term_buffer::Int;                  # Width of final adaptation interval
    window::Int;                       # Initial width slow adaptation interval

    # Algorithm fields
    algorithm::Symbol;                 # :hmc or :fixed_param
    # HMC specific fields
    engine::Symbol;                    # :nuts or :static (default = :nuts)
    # NUTS specific field
    max_depth::Int;                    # Maximum tree depth (> 0, default=10)
    # Static specific field
    int_time::Float64;                 # Static integration time

    # HMC remaining fields
    metric::Symbol;                    # :diag_e, :unit_e, :dense_e
    metric_file::AbstractString;       # Precompiled Euclidean metric
    stepsize::Float64;                 # Stepsize for discrete evolution
    stepsize_jitter::Float64;          # Uniform random jitter of stepsize (%)

    # Output files
    output_base::AbstractString;       # Used for file paths to be created
    # Tmpdir setting
    tmpdir::AbstractString;            # Holds all created files
    # Cmdstan path
    exec_path::AbstractString;         # Path to the cmdstan excutable

    use_json::Bool;                    # Use JSON for data and init files

    # Data and init file paths
    data_file::Vector{AbstractString}; # Array of data files input to cmdstan
    init_file::Vector{AbstractString}; # Array of init files input to cmdstan
    # Generated command line vector
    cmds::Vector{Cmd};                 # Array of cmds to be spawned/pipelined
    # Files created by cmdstan
    sample_file::Vector{String};       # Sample file array (.csv)
    log_file::Vector{String};          # Log file array
    diagnostic_file::Vector{String};   # Diagnostic file array

    # Stansummary settings
    summary::Bool;                     # Store cmdstan's summary as a .csv file
    print_summary::Bool;               # Print the summary

    # CMDSTAN_HOME
    cmdstan_home::AbstractString;      # Directory where cmdstan can be found
end

"""

Create a SampleModel and compile Stan Language Model.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `name::AbstractString`               # Name for the model
* `model::AbstractString`              # Stan model source
```

### Optional positional argument
```julia
* `tmpdir`                             # Directory where output files are stored
                                       # Default: `mktempdir()`
```

Note: On Windows I have seen issues using `tmpdir`.

"""
function SampleModel(name::AbstractString, model::AbstractString,
    tmpdir=mktempdir())

    !isdir(tmpdir) && mkdir(tmpdir)

    update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))

    output_base = joinpath(tmpdir, name)
    exec_path = executable_path(output_base)
    cmdstan_home = CMDSTAN_HOME

    error_output = IOBuffer()
    is_ok = cd(cmdstan_home) do
        success(pipeline(`$(make_command()) -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`;
            stderr = error_output))
    end

    if !is_ok
        throw(StanModelError(name, String(take!(error_output))))
    end

    SampleModel(name, model,
        # num_threads
        4,
        # check_num_chains, use_cpp_chains
        true, false,
        # num_cpp_chains
        1,
        # num_julia_chains
        1,
        # num_chains
        4,
        # num_samples, num_warmups, save_warmups
        1000, 1000, false,
        # thin, seed, refresh, init_bound
        1, -1, 100, 2,
        # Adapt fields
        # engaged, gamma, delta, kappa, t0, init_buffer, term_buffer, window
        true, 0.05, 0.8, 0.75, 10, 75, 50, 25,
        # algorithm fields
        :hmc,                          # or :static
        # engine, max_depth
        :nuts, 10,
        # Static engine specific fields
        2pi,
        # metric, metric_file, stepsize, stepsize_jitter
        :diag_e, "", 1.0, 0.0,

        # Ouput settings
        output_base,                   # Output base
        tmpdir,                        # Tmpdir settings
        exec_path,                     # exec_path

        true,                          # Use JSON for cmdstan input files
        AbstractString[],              # Data files
        AbstractString[],              # Init files

        Cmd[],                         # Command lines
        String[],                      # Sample .csv files
        String[],                      # Log files
        String[],                      # Diagnostic files
        true,                          # Create stansummary result
        false,                         # Display stansummary result
        cmdstan_home
    )
end

function Base.show(io::IO, ::MIME"text/plain", m::SampleModel)
    println("\nModel name:")
    println(io, "  name =                    $(m.name)")
    println("\nC++ threads per forked process:")
    println(io, "  num_threads =             $(m.num_threads)")
    println(io, "  use_cpp_chains =          $(m.use_cpp_chains)")
    println(io, "  check_num_chains =        $(m.check_num_chains)")
    println("\nC++ chains per forked process:")
    println(io, "  num_cpp_chains =          $(m.num_cpp_chains)")
    println("\nNo of forked Julia processes:")
    println(io, "  num_julia_chains =        $(m.num_julia_chains)")
    println("\nActual number of chains:")
    println(io, "  num_chains =              $(m.num_chains)")

    println(io, "\nSample section:")
    println(io, "  num_samples =             ", m.num_samples)
    println(io, "  num_warmups =             ", m.num_warmups)
    println(io, "  save_warmup =             ", m.save_warmup)
    println(io, "  thin =                    ", m.thin)
    println(io, "  seed =                    ", m.seed)
    println(io, "  refresh =                 ", m.refresh)
    println(io, "  init_bound =              ", m.init_bound)
    println(io, "\nAdapt section:")
    println(io, "  engaged =                 ", m.engaged)
    println(io, "  gamma =                   ", m.gamma)
    println(io, "  delta =                   ", m.delta)
    println(io, "  kappa =                   ", m.kappa)
    println(io, "  t0 =                      ", m.t0)
    println(io, "  init_buffer =             ", m.init_buffer)
    println(io, "  term_buffer =             ", m.term_buffer)
    println(io, "  window =                  ", m.window)
    if m.algorithm ==:hmc
        println("\nAlgorithm section:")
        println(io, "\n  algorithm =               $(m.algorithm)")
        if m.engine == :nuts
            println(io, "\n    NUTS section:")
            println(io, "      engine =              $(m.engine)")
            println(io, "      max_depth =           ", m.max_depth)
        elseif m.engine == :static
            println(io, "\n  STATIC section:")
            println(io, "    engine =               :static")
            println(io, "    int_time =             ", m.int_time)
        end
        println(io, "\n  Metric section:")
        println(io, "    metric =                ", m.metric)
        println(io, "    stepsize =              ", m.stepsize)
        println(io, "    stepsize_jitter =       ", m.stepsize_jitter)
    else
        if m.algorithm == :fixed_param
            println(io, "    algorithm =         :fixed_param")
        else
            println(io, "    algorithm =         Unknown")
        end
    end
    println(io, "\nData and init files:")
    println(io, "  use_json =                ", m.use_json)
    println(io, "\nStansummary section:")
    println(io, "  summary                   ", m.summary)
    println(io, "  print_summary             ", m.print_summary)
    println(io, "\nOther:")
    println(io, "  output_base =             ", m.output_base)
    println(io, "  tmpdir =                  ", m.tmpdir)
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
#-----------------------------#
# see page 136 in Sarkka book #
#-----------------------------#

function smoothingrecursion(y; A = A, H = H, Q = Q, R = R, m = m, P = P)

    N = length(y); @assert(N == length(m) == length(P))

    # Store here state distribution

    μstate, Σstate = Array{Vector{Float64}}(undef, N), Array{Matrix{Float64}}(undef, N)

    # Start recursion with prior mean and prior covariance

    mˢₖ₊₁, Pˢₖ₊₁ = m[end], P[end]

    μstate[N], Σstate[N] = mˢₖ₊₁, Pˢₖ₊₁

    for k in N-1:-1:1

        m⁻ₖ₊₁ = A * m[k]

        P⁻ₖ₊₁ = A * P[k] * A' + Q

        Gₖ = (P[k] * A') / P⁻ₖ₊₁

        mˢₖ = m[k] + Gₖ * (mˢₖ₊₁ - m⁻ₖ₊₁)

        Pˢₖ = P[k] + Gₖ * (Pˢₖ₊₁ - P⁻ₖ₊₁) * (Gₖ)'

        # store mean and covariance of Gaussian distribution N(xₖ|mˢₖ, Pˢₖ)
        Pˢₖ = (Pˢₖ + (Pˢₖ)') / 2

        μstate[k], Σstate[k] = mˢₖ, Pˢₖ

        # previous state

        mˢₖ₊₁, Pˢₖ₊₁ = mˢₖ, Pˢₖ

    end

    return μstate, Σstate

end

function testsmoothing(seed=1)

    yclean, y, x, A, H, Q, R, m₀, P₀ = simulatedata(seed=seed)

    μfilter, Σfilter = filteringrecursion(y; A = A, H = H, Q = Q, R = R, m₀ = m₀, P₀ = P₀)

    μsmooth, Σsmooth = smoothingrecursion(y; A = A, H = H, Q = Q, R = R, m = μfilter, P = Σfilter)

    figure()
    subplot(211)
    plot([yᵢ[1]        for yᵢ in yclean], label = "simulated 1")
    plot([(H*xᵢ)[1]    for xᵢ in x], label = "simulated 1")
    legend()

    subplot(212)
    plot([xᵢ[1]    for xᵢ in x], label = "simulated 1")
    plot([xᵢ[1]    for xᵢ in μfilter], label = "filtered 1")
    plot([xᵢ[1]    for xᵢ in μsmooth], "--", label = "smoothed 1")
    legend()

    acc = accfilter = accsmooth = 0.0

    for i in 1:length(y)

        acc += sum((yclean[i] - H*x[i]).^2) / length(y)

        accfilter += sum((yclean[i] - H*μfilter[i]).^2)/ length(y)

        accsmooth += sum((yclean[i] - H*μsmooth[i]).^2)/ length(y)

    end

    @show acc, accfilter, accsmooth

    accx = accfilterx = accsmoothx = 0.0

    for i in 1:length(y)

        accx += sum((x[i] - x[i]).^2)

        accfilterx += sum((x[i] - μfilter[i]).^2)

        accsmoothx += sum((x[i] - μsmooth[i]).^2)

    end

    @show accx, accfilterx, accsmoothx

    x, μfilter, μsmooth

end


text:
export @fermiview

macro fermiview()
    quote
        try
            let mol = Molecule()
                Alchemy.plot(mol.atoms)
            end
        catch UndefVarError
            @error "Fermi.jl could not be found, try running `using Fermi` first."
        end
    end |> esc
end

max_stars_repo_name: gustavojra/Alchemy

max_issues_repo_name: gustavojra/Alchemy

max_forks_repo_name: gustavojra/Alchemy


text:
function basic_routing_model_unitary(m::Model, rd::RoutingData, ::PathFormulation)
    @variable(m, 0 <= routing[d in demands(rd), p in rd.demand_to_path_ids[d]] <= 1)
    @constraint(m, c_path[d in demands(rd)], sum(routing[d, p] for p in rd.demand_to_path_ids[d]) == 1)
    return RoutingModel(rd, m, UnitaryFlows, routing, constraints_convexity=c_path)
end

function basic_routing_model_unscaled(m::Model, rd::RoutingData, dm, ::PathFormulation)
    # Either dm is known (::Dict{Edge{Int}, Float64}) or is it a vector of variables (::DenseAxisArray).
    @variable(m, routing[d in demands(rd), p in rd.demand_to_path_ids[d]] >= 0)
    @constraint(m, c_path[d in demands(rd)], sum(routing[d, p] for p in rd.demand_to_path_ids[d]) == dm[d])
    return RoutingModel(rd, m, UnscaledFlows, routing, constraints_convexity=c_path)
end

function total_flow_in_edge(rm::RoutingModel, e::Edge, dm::Dict{Edge{Int}, Float64}, ::PathFormulation, ::Val{UnitaryFlows})
    flow = AffExpr(0.0)

    paths = find_path_ids_with_edge(rm.data, e)
    for p in paths
        d = rm.data.path_id_to_demand[p]
        if d in keys(dm) && ! iszero(dm[d])
            add_to_expression!(flow, dm[d] * rm.routing[d, p])
        end
    end

    return flow
end

function total_flow_in_edge(rm::RoutingModel, e::Edge, ::PathFormulation, ::Val{UnscaledFlows})
    paths = find_path_ids_with_edge(rm.data, e)
    if isempty(paths) # sum() over an empty array throws an error...
        return AffExpr(0.0)
    end
    return sum(rm.routing[rm.data.path_id_to_demand[p], p] for p in paths)
end


text:
using Compat
using Compat.Test
import Compat.Pkg
using MathProgBase
using CPLEX

@testset "$folder" for folder in ["C_API", "MathProgBase", "MathOptInterface"]
    @testset "$(file)" for file in readdir(folder)
        include(joinpath(folder, file))
    end
end


text:
import NumericalIntegration

"""
    trapz(x, y)

Compute the area under the curve of 2-dimensional points (x, y) using
the trapezoidal method.
"""
function trapz end

trapz(x, y) = trapz(promote(x, y)...)

function trapz(
        x::AbstractVector{T},
        y::AbstractVector{T},
        )::T where T
    if length(x) != length(y)
        error("length(x) != length(y)")
    end
    if length(x) == 0
        error("length(x) == 0")
    end
    if !all(x .== sort(x; rev = false))
        error("x needs to be sorted in ascending order")
    end
    twoI::T = zero(T)
    for k = 2:length(x)
        twoI += ( y[k] + y[k-1] ) * ( x[k] - x[k-1] )
    end
    I_verify::T = NumericalIntegration.integrate(
        x,
        y,
        NumericalIntegration.Trapezoidal(),
        )
    @assert isapprox(twoI/2, I_verify; atol=0.00000001)
    return I_verify
end


text:
using Test

include("../src/RSA.jl")

using .RSA

@testset "Secondary functions" begin

    @test RSA.exp!(2, 47) == 199807
    @test RSA.exp!(3, 51, 678) == 93

    @test RSA.inverse(2, 5) == 3
    @test RSA.inverse(17, 39) == 23
    @test_throws ErrorException RSA.inverse(2, 4)

    @test RSA.seive(1) == []
    @test RSA.seive(13) == [2, 3, 5, 7, 11, 13]

    @test RSA.phi!(10) == 4
    @test RSA.phi!(7) == 6
    @test RSA.phi!(33500000) == 13200000

    @test RSA.get_random_int(3) <= 7
    @test RSA.get_random_int(8) <= 255

    @test RSA.jacobi!(1001, 9907) == -1
    @test RSA.jacobi!(19, 45) == 1
    @test RSA.jacobi!(8,21) == -1
    @test RSA.jacobi!(3, 15) == 0
    @test_throws ErrorException RSA.jacobi!(13, 28)
    @test RSA.jacobi!(28, 13) == -1
end

@testset "Primality functions" begin

    @test RSA.miller_rabin(561) == false
    @test RSA.miller_rabin(29) == true
    @test RSA.miller_rabin(221) == false

    @test RSA.solovay_strassen(561) == false
    @test RSA.solovay_strassen(29) == true
    @test RSA.solovay_strassen(221) == false

    @test RSA.fermat_test(5) == true
    @test RSA.fermat_test(4) == false
    @test RSA.fermat_test(341) == false
end

@testset "Processing" begin

    a =""
    @test RSA.recover_string!(RSA.process_string!(a)) == a
    a ="123"
    @test RSA.recover_string!(RSA.process_string!(a)) == a
    a = "1234567890"
    @test RSA.recover_string!(RSA.process_string!(a)) == a
end

@testset "Encryption" begin
    B = RSASystem()
    a = "ABCDEFG"
    @test decrypt!(B, encrypt!(a, B.public_key)) == a
    a = "AB"
    @test decrypt!(B, encrypt!(a, B.public_key)) == a
    a = "ABABABABABABABABABABABABABABABAB123"
    @test decrypt!(B, encrypt!(a, B.public_key)) == a
end


text:
f(x) = cos(x) * x
sumsq(xs) = sum(xs.^2)

function foo1(as :: Vector{Float64}, b :: Float64)
    p = length(as)
    p = p - 1
    if p > 0
        if as[1] > 0
            y = [sin(a) for a in as] .* f(b)
        else
            y = [1.1*p, foo1(as[2:end], 1.1)]
        end
    else
        y = -as.*f(b)
    end
    f(sumsq(y)) + 5.5555
end

max_stars_repo_name: tomjaguarpaw/knossos-ksc

max_issues_repo_name: tomjaguarpaw/knossos-ksc

max_forks_repo_name: tomjaguarpaw/knossos-ksc


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 2e1d1184-4d32-453f-804b-a0d19d2f3184
include("..\\..\\jl\\commasplit.jl")

# ╔═╡ 03375a9f-a288-4f3b-a6f6-595fba65b042
include("..\\..\\jl\\str2int.jl")

# ╔═╡ 20114b30-2029-4513-9103-5b27cb1885cb
include("..\\..\\jl\\flatten.jl")

# ╔═╡ 87c8c810-657e-11ec-2114-4b4ad4af6ec8
md"""
# Day 6: Lanternfish

**Reference:** [Day 6 - Advent of Code 2021](https://adventofcode.com/2021/day/6)

**Source:** [ljk233/AdventOfCode2021/days/06/d6.jl](https://github.com/ljk233/AdventOfCode2021/blob/main/days/06/d6.jl)
"""

# ╔═╡ b86823bf-1386-449f-8940-17aa45a0d242
md"""
## Imports
"""

# ╔═╡ e05f10b5-f57e-42ca-8ee6-5df06c0baa59
md"""
## Functions
"""

# ╔═╡ 47e4af0d-02ba-4ec8-ae11-87230eff4d32
"""
	scan(line: String) -> Integer[]
"""
scan(line) = line |> commasplit |> flatten .|> str2int

# ╔═╡ 70ceeb78-51ad-4821-8d40-b30883fb8a71
"""
	survey(ocean: Integer[]) -> Integer[]
"""
function survey(ocean)
	state = zeros(Integer, 8)
	for lf ∈ ocean
		state[lf] += 1
	end
	pushfirst!(state, 0)
	return state
end

# ╔═╡ 16929d30-5cbe-4f6e-9034-2f2b0c78065b
"""
	simulate!(ocean: Integer[]) -> NoReturn
"""
function simulate!(ocean)
	push!(ocean, popfirst!(ocean))
	ocean[7] += last(ocean)
	return nothing
end

# ╔═╡ 1a9694c2-923b-4f98-83f6-2b2bf3ff55ab
"""
	main(f: String, days: Integer) -> Integer
"""
function main(f, days)
	ocean::Vector{Integer} = readline(f) |> scan |> survey
	for day ∈ 1:days
		simulate!(ocean)
	end
	return sum(ocean)
end

# ╔═╡ 4c202df3-62b5-4446-b3e1-0afae48210b8
md"""
## Part One

Find a way to simulate lanternfish.

**How many lanternfish would there be after 80 days?**
"""

# ╔═╡ c2a8fd18-c49e-48e0-8780-0591a3c59fa6
md"""
**Test #1.**
Expected result = 5934.
"""

# ╔═╡ a5ac2615-cae0-4a16-901d-687fd46f3145
main("test.in", 80)

# ╔═╡ a957d84e-b1a4-4a95-8539-fcaf07f1f1e5
md"""
**Solution.**
Confirmed result = 388419.
"""

# ╔═╡ d6b93da4-1ef1-48d6-83fe-8424acbcfaf2
main("data.in", 80)

# ╔═╡ cede829f-55f0-4f4d-96b5-a380acd16394
md"""
## Part Two

Suppose the lanternfish live forever and have unlimited food and space.
Would they take over the entire ocean?

**How many lanternfish would there be after 256 days?**
"""

# ╔═╡ 2d234229-dbea-40cc-be6e-0267a04bd1b0
md"""
**Test #2.**
Expected result = 26984457539.
"""

# ╔═╡ e099ec40-a9fe-4446-9db7-9619cb3adc5d
main("test.in", 256)

# ╔═╡ a9688396-74a2-4989-95bf-381b7904f521
md"""
**Solution.**
Confirmed result = 1740449478328..
"""

# ╔═╡ 488943e4-d1ba-4818-b8d9-1f886dc7529f
main("data.in", 256)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─87c8c810-657e-11ec-2114-4b4ad4af6ec8
# ╟─b86823bf-1386-449f-8940-17aa45a0d242
# ╠═2e1d1184-4d32-453f-804b-a0d19d2f3184
# ╠═03375a9f-a288-4f3b-a6f6-595fba65b042
# ╠═20114b30-2029-4513-9103-5b27cb1885cb
# ╟─e05f10b5-f57e-42ca-8ee6-5df06c0baa59
# ╠═47e4af0d-02ba-4ec8-ae11-87230eff4d32
# ╠═70ceeb78-51ad-4821-8d40-b30883fb8a71
# ╠═16929d30-5cbe-4f6e-9034-2f2b0c78065b
# ╠═1a9694c2-923b-4f98-83f6-2b2bf3ff55ab
# ╟─4c202df3-62b5-4446-b3e1-0afae48210b8
# ╟─c2a8fd18-c49e-48e0-8780-0591a3c59fa6
# ╠═a5ac2615-cae0-4a16-901d-687fd46f3145
# ╟─a957d84e-b1a4-4a95-8539-fcaf07f1f1e5
# ╠═d6b93da4-1ef1-48d6-83fe-8424acbcfaf2
# ╟─cede829f-55f0-4f4d-96b5-a380acd16394
# ╟─2d234229-dbea-40cc-be6e-0267a04bd1b0
# ╠═e099ec40-a9fe-4446-9db7-9619cb3adc5d
# ╟─a9688396-74a2-4989-95bf-381b7904f521
# ╠═488943e4-d1ba-4818-b8d9-1f886dc7529f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# generalized ensemble algorithms such as
# * muca / WL
# * population annealing

#Idea: These structs could actually carry information! E.g the current weight factor for WangLandau (then mutable)
#      in principle also the weights ...
struct Multicanonical end
struct WangLandau end #todo (in terms of log_dos instead of log_weight)

###### discuss: maybe move this to importance sampling (because it is importance sampling)
@doc """
    accept(alg::Multicanonical, rng::AbstractRNG, log_weight::Vector{Float64}, arg_new::T, arg_new::T)::Bool where T

Multicanonical acceptance function
```p(x\\to x^') = \\text{min}\\left(1, W(x^')/W(x)\\right)```
"""
function accept(alg::Multicanonical, rng::AbstractRNG, log_weight::Histogram, arg_new::T, arg_old::T)::Bool where T
    difference = log_weight[args_new] - log_weight[args_old]
    if difference > 0
        return true
    elseif rand(rng) < exp(difference)
        return true
    else
        return false
    end
end

"""

mode can be
* :simple
* :extrapolate
* ... ?
"""
function update_weights(alg::Multicanonical, log_weight::Histogram, histogram::Histogram; mode::Symbol=:simple)
end


text:
using Test, MakieCore

# Main tests live in Makie.jl, but we should write some unit tests going forward!
using MakieCore: @recipe, Attributes, Plot
import MakieCore: plot!, convert_arguments, used_attributes, plot
import MakieCore: Observable, PointBased

struct AbstractTimeseriesSolution
    results::Vector{Float32}
end

function plot!(plot::Plot(AbstractTimeseriesSolution))
    # plot contains any keyword arguments that you pass to plot(series; kw...)
    var = get(plot, :var, Observable(5))
    density!(plot, map((v, r)-> v .* r.results, var, plot[1]))
end

struct Test2
    series::Any
end

struct Solution
    data::Any
end

function plot!(plot::Plot(Test2))
    arg1 = plot[1]
    scatter!(plot, arg1[].series)
    ser = AbstractTimeseriesSolution(arg1[].series)
    sol = Solution(arg1[].series)
    plot!(plot, ser, var = 10)
    scatter!(plot, sol, attribute = 3, color=:red)
end

used_attributes(::Any, x::Solution) = (:attribute,)

# Convert for all point based types (lines, scatter)
function convert_arguments(p::MakieCore.PointBased, x::Solution; attribute = 1.0)
    return convert_arguments(p, x.data .* attribute)
end


text:
abstract type AbstractDataContainer end

n_sample(d::AbstractDataContainer) = d.n_sample
n_dim(d::AbstractDataContainer) = d.n_dim
n_output(d::AbstractDataContainer) = 1

input(d::AbstractDataContainer) = d.X
output(d::AbstractDataContainer) = d.y

struct DataContainer{Tx<:Real,TX<:AbstractVector,Ty<:Real,TY<:AbstractArray} <:
       AbstractDataContainer
    X::TX # Feature vectors
    y::TY # Output (-1,1 for classification, real for regression, bitmatrix for multiclass)
    n_sample::Int # Number of samples
    n_dim::Int # Number of features per sample
end

function wrap_data(X::TX, y::TY) where {TX,TY<:AbstractMatrix}
    size(y, 1) == length(X) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(size(y, 1)))",
    )
    Tx = eltype(first(X))
    Ty = eltype(first(y))
    return DataContainer{Tx,TX,Ty,TY}(X, y, length(X), length(first(X)))
end

function wrap_data(X::TX, y::TY) where {TX,TY<:AbstractVector{<:Real}}
    length(y) == length(X) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length(y)))",
    )
    Tx = eltype(first(X))
    Ty = eltype(first(y))
    return DataContainer{Tx,TX,Ty,TY}(X, y, length(X), length(first(X)))
end

function wrap_data(X::TX, y::TY) where {TX,TY<:AbstractVector}
    all(length.(y) .== length(X)) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length.(y))))",
    )
    Tx = eltype(first(X))
    Ty = eltype(first(y))
    return DataContainer{Tx,TX,Ty,TY}(X, y, length(X), length(first(X)))
end

struct MODataContainer{Tx<:Real,TX<:AbstractVector,TY<:AbstractVector} <:
       AbstractDataContainer
    X::TX # Feature vectors
    y::TY # Output ({-1,1} for classification, real for regression, matrix for multiclass)
    n_sample::Int # Number of samples
    n_dim::Int # Number of features per sample
    n_output::Int # Number of outputs
end

function wrap_modata(X::TX, y::TY) where {TX,TY<:AbstractVector}
    all(length.(y) .== length(X)) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length.(y))))",
    )
    Tx = eltype(first(X))
    return MODataContainer{Tx,TX,TY}(X, y, length(X), length(first(X)), length(y))
end

n_output(d::MODataContainer) = d.n_output

function wrap_X(X::AbstractMatrix{T}, obsdim::Int=1) where {T<:Real}
    return KernelFunctions.vec_of_vecs(X; obsdim=obsdim), T
end

function wrap_X(X::AbstractVector{T}, ::Int=1) where {T<:Real}
    return X, T
end

function wrap_X(X::AbstractVector{<:AbstractVector{T}}, ::Int=1) where {T<:Real}
    return X, T
end


text:
using Test

import Petri
using LabelledArrays
using AlgebraicPetri
using AlgebraicPetri.Epidemiology
using Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Programs

@testset "Core" begin
    include("core.jl")
end

@testset "Types" begin
    include("types.jl")
end

@testset "Petri" begin
    include("petri.jl")
end

@testset "Epidemiology" begin
    include("epidemiology.jl")
end

@testset "BilayerNetworks" begin
    include("bilayernetworks.jl")
end


text:
# # [Linear NH mountain waves (Topography)](@id EX-LIN_NH-docs)
#
# ## Description of experiment
# 1) Dry linear Non-hydrostatic Mountain Waves
# This example of a non-linear hydrostatic mountain wave can be classified as an initial value
# problem.
#
# The atmosphere is dry and the flow impinges against a witch of Agnesi mountain of heigh $h_m=1m$
# and base parameter $a=1000$m and centered in $x_c = 72km$ in a 2D domain
# $\Omega = 144km \times 30 km$. The mountain is defined as
#
# ```math
#  z = \frac{h_m}{1 + \frac{x - x_c}{a}}
# ```
# The 2D problem is setup in 3D by using 1 element in the y direction.
# To damp the upward moving gravity waves, a Reyleigh absorbing layer is added at $z = 10,000m$.
#
# The initial atmosphere is defined such that it has a stability frequency $N=0.01 s^{-1}$, where
#
# ```math
#  N^2 = g\frac{\rm d \ln \theta}{ \rm dz}
# ```
# so that
#
# ```math
# \theta = \theta_0 \exp\left(\frac{N^2 z}{g} \right),
# ```
# ```math
# \pi = 1 + \frac{g^2}{c_p \theta_0 N^2}\left(\exp\left(\frac{-N^2 z}{g} \right)\right)
# ```
#
# where $\theta_0 = 280K$.
#
# so that
#
# $ρ = \frac{p_{sfc}}{R_{gas}\theta}pi^{c_v/R_{gas}}$
# and $T = \theta \pi$
#
# 2) Boundaries
#    - `Impenetrable(FreeSlip())` - Top and bottom: no momentum flux, no mass flux through
#      walls.
#    - `Impermeable()` - non-porous walls, i.e. no diffusive fluxes through
#       walls.
#    - Agnesi topography built via meshwarp.
#    - Laterally periodic
# 3) Domain - 144,000 m (horizontal) x 1360 m (horizontal) x 30,000m (vertical) (infinite domain in y)
# 4) Resolution - 340 m X 200 m effective resolution
# 5) Total simulation time - 18,000 s
# 6) Overrides defaults for
#    - CPU Initialisation
#    - Time integrator
#    - Sources
#
#md # !!! note
#md #     This experiment setup assumes that you have installed the
#md #     `ClimateMachine` according to the instructions on the landing page.
#md #     We assume the users' familiarity with the conservative form of the
#md #     equations of motion for a compressible fluid (see the
#md #     [AtmosModel](@ref AtmosModel-docs) page).
#md #
#md #     The following topics are covered in this example
#md #     - Defining the initial conditions
#md #     - Applying source terms
#md #     - Add an idealized topography defined by a warping function
#
# ## Boilerplate (Using Modules)
#
# The setup of this problem is taken from Case 6 of [giraldoRestelli2008a](@cite)
#
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Atmos.SubgridScale: C_smag
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# ## [Initial Conditions](@id init)
#md # !!! note
#md #     The following variables are assigned in the initial condition
#md #     - `state.ρ` = Scalar quantity for initial density profile
#md #     - `state.ρu`= 3-component vector for initial momentum profile
#md #     - `state.energy.ρe`= Scalar quantity for initial total-energy profile
#md #       humidity
function init_agnesi_hs_lin!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    ## Problem float-type
    FT = eltype(state)
    param_set = parameter_set(bl)

    ## Unpack constant parameters
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v

    c::FT = c_v / R_gas
    c2::FT = R_gas / c_p

    ## Define initial thermal field as isothermal
    Tiso::FT = 250.0
    θ0::FT = Tiso

    ## Assign a value to the Brunt-Vaisala frquencey:
    Brunt::FT = 0.01
    Brunt2::FT = Brunt * Brunt
    g2::FT = _grav * _grav

    π_exner::FT = exp(-_grav * z / (c_p * Tiso))
    θ::FT = θ0 * exp(Brunt2 * z / _grav)
    ρ::FT = p0 / (R_gas * θ) * (π_exner)^c

    ## Compute perturbed thermodynamic state:
    T = θ * π_exner
    e_int = internal_energy(param_set, T)
    ts = PhaseDry(param_set, e_int, ρ)

    ## initial velocity
    u = FT(10.0)

    ## State (prognostic) variable assignment
    e_kin = FT(0)                                       # kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)# potential energy
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)         # total energy

    state.ρ = ρ
    state.ρu = SVector{3, FT}(ρ * u, 0, 0)
    state.energy.ρe = ρe_tot
end

function setmax(f, xmax, ymax, zmax)
    function setmaxima(xin, yin, zin)
        return f(xin, yin, zin; xmax = xmax, ymax = ymax, zmax = zmax)
    end
    return setmaxima
end

# Define a warping function to build an analytic topography:
function warp_agnesi(xin, yin, zin; xmax = 1000.0, ymax = 1000.0, zmax = 1000.0)

    FT = eltype(xin)

    ac = FT(1000)
    hm = FT(1)
    xc = FT(0.5) * xmax
    zdiff = hm / (FT(1) + ((xin - xc) / ac)^2)

    ## Linear relaxation towards domain maximum height
    x, y, z = xin, yin, zin + zdiff * (zmax - zin) / zmax
    return x, y, z
end

# ## [Model Configuration](@id config-helper)
# We define a configuration function to assist in prescribing the physical
# model. The purpose of this is to populate the
# `AtmosLESConfiguration` with arguments
# appropriate to the problem being considered.
function config_agnesi_hs_lin(
    ::Type{FT},
    N,
    resolution,
    xmax,
    ymax,
    zmax,
) where {FT}
    ##
    ## Explicit Rayleigh damping:
    ##
    ## ``
    ##   \tau_s = \alpha * \sin\left(0.5\pi \frac{z - z_s}{zmax - z_s} \right)^2,
    ## ``
    ## where
    ## ``sponge_ampz`` is the wave damping coefficient (1/s)
    ## ``z_s`` is the level where the Rayleigh sponge starts
    ## ``zmax`` is the domain top
    ##
    ## Setup the parameters for the gravity wave absorbing layer
    ## at the top of the domain
    ##
    ## u_relaxation(xvelo, vvelo, wvelo) contains the background velocity values to which
    ## the sponge relaxes the vertically moving wave
    u_relaxation = SVector(FT(10), FT(0), FT(0))

    ## Wave damping coefficient (1/s)
    sponge_ampz = FT(0.5)

    ## Vertical level where the absorbing layer starts
    z_s = FT(10000.0)

    ## Pass the sponge parameters to the sponge calculator
    rayleigh_sponge =
        RayleighSponge{FT}(zmax, z_s, sponge_ampz, u_relaxation, 2)

    ## Setup the source terms for this problem:
    source = (Gravity(), rayleigh_sponge)

    ## Define the reference state:
    T_virt = FT(280)
    temp_profile_ref = IsothermalProfile(param_set, T_virt)
    ref_state = HydrostaticState(temp_profile_ref)
    nothing # hide

    _C_smag = FT(0.0)
    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = Vreman(_C_smag),
        moisture = DryModel(),
        tracers = NoTracers(),
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = init_agnesi_hs_lin!,
        source = source,
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "Agnesi_NH_LINEAR",      # Problem title [String]
        N,                       # Polynomial order [Int]
        resolution,              # (Δx, Δy, Δz) effective resolution [m]
        xmax,                    # Domain maximum size [m]
        ymax,                    # Domain maximum size [m]
        zmax,                    # Domain maximum size [m]
        param_set,               # Parameter set.
        init_agnesi_hs_lin!,     # Function specifying initial condition
        model = model,           # Model type
        meshwarp = setmax(warp_agnesi, xmax, ymax, zmax),
    )

    return config
end

# Define a `main` method (entry point)
function main()

    FT = Float64

    ## Define the polynomial order and effective grid spacings:
    N = 4

    ## Define the domain size and spatial resolution
    Nx = 20
    Ny = 20
    Nz = 20
    xmax = FT(144000)
    ymax = FT(4000)
    zmax = FT(30000)
    Δx = xmax / FT(Nx)
    Δy = ymax / FT(Ny)
    Δz = zmax / FT(Nz)
    resolution = (Δx, Δy, Δz)

    t0 = FT(0)
    timeend = FT(100)

    ## Define the max Courant for the time time integrator (ode_solver).
    ## The default value is 1.7 for LSRK144:
    CFL = FT(1.5)

    ## Assign configurations so they can be passed to the `invoke!` function
    driver_config = config_agnesi_hs_lin(FT, N, resolution, xmax, ymax, zmax)

    ## Define the time integrator:
    ## We chose an explicit single-rate LSRK144 for this problem
    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = CFL,
    )

    ## Set up the spectral filter to remove the solutions spurious modes
    ## Define the order of the exponential filter: use 32 or 64 for this problem.
    ## The larger the value, the less dissipation you get:
    filterorder = 64
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end
    ## End exponential filter

    ## Invoke solver (calls `solve!` function for time-integrator),
    ## pass the driver, solver and diagnostic config information.
    result = ClimateMachine.invoke!(
        solver_config;
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )

    ## Check that the solution norm is reasonable.
    @test FT(0.9) < result < FT(1)
    ## Some energy is lost from the sponge, so we cannot
    ## expect perfect energy conservation.
end

# Call `main`
main();


text:
module lattice

include("load_kwant.jl")
import .. AbstractKwantObject

abstract type AbstractLattice <:AbstractKwantObject end

struct Monatomic <: AbstractLattice o::PyObject end
struct Polyatomic <: AbstractLattice o::PyObject end

(lat::Monatomic)(i::Integer,j::Integer) = pycall(lat.o,PyObject,i,j)

square(a=1;kwargs...) = Monatomic(kwant.lattice.square(a;kwargs...))

general(arg;kwargs...) = Monatomic(kwant.lattice.general(arg;kwargs...))
general(arg1,arg2;kwargs...) = Polyatomic(kwant.lattice.general(arg1,arg2;kwargs...))

# function Base.getproperty(k::AbstractLattice, name::Symbol)
    # if name==:o
        # p = getfield(k,:o)
    # else
        # p = getproperty(k.o,name)
    # end
# #     if name==:sublattices
# #         p = Monatomic.(p)
#     end
    # return p
# end

end # module

# function Base.getproperty(k::AbstractLattice, name::Symbol)
#     if name==:o
#         p = getfield(k,:o)
#     else
#         p = getproperty(k.o,name)
#     end
#     if name==:neighbors
#         p = _Neighbors(p)
#     end
#     return p
# end


text:
"""
    lagrange(x, xd, yd)
adalah fungsi yang digunakan untuk mencari nilai interpolasi pada titik/vektor `x`, jika
diketahui suatu himpunan pasangan terurut `(xd,yd)`.

# Example
```jl
julia> xd = [1,2,3,5];

julia> yd = [1.06 1.12 1.34 1.78];

julia> lagrange(4,xd,yd)
1.6000000000000003

julia> lagrange([2.5,4,5.5],xd,yd)
3-element Array{Float64,1}:
 1.2175
 1.6000000000000003
 1.8025000000000002
```
return solusi hampiran interpolasi `y`.
"""
function lagrange(x, xd, yd)
  #% sedikit trik agar bisa menghitung nilai interpolasi pada banyak
  #% titik sekaligus. Intinya, tetap diproses satu titik demi titik.
  m=length(x);
  y=zeros(m)
  if m > 1
    for i=1:m
      # proses satu titik demi satu titik
      y[i]=lagrange(x[i], xd, yd);
    end
    return y
  end
  #% periksa jumlah titik dan tentukan derajat polinom
  ntitik = length(xd);
  n = ntitik-1;  #% derajat maksimum, sesuai jumlah titik yang ada
  #% hitung L_{n,k}(x) untuk k=1:(n+1)
  Ln=ones(1,ntitik);
  for i=1:ntitik
    for j=1:ntitik
      if i!=j
        Ln[i] = Ln[i] * (x-xd[j])/(xd[i]-xd[j]);;
      end
    end
  end
  #% hitung P_{N}(x)
  y = 0;
  for i=1:ntitik
    y = y + yd[i]*Ln[i];
  end
  return y
end


text:
function mean_std_video_gpu(han::Tracker_Handles,total_frame_num)
    mean_std_video_gpu(han.wt.vid_name,total_frame_num)
end

function mean_std_video_gpu(vid_name::String,total_frame_num,max_intensity=255,loading_size=500)

    (w,h,fps)=get_vid_dims(vid_name)
    start_frame = 0

    load_number = div(total_frame_num - start_frame, loading_size)-1

    temp_frames=zeros(UInt8,w,h,loading_size)

    temp_frames2 = convert(KnetArray,zeros(Float32,w,h,loading_size))

    running_mean = convert(KnetArray,zeros(Float32,w,h,1))
    running_mean_i = convert(KnetArray,zeros(Float32,w,h,1))

    running_std = convert(KnetArray,zeros(Float32,w,h,1))
    running_std_i = convert(KnetArray,zeros(Float32,w,h,1))

    @ffmpeg_env run(WhiskerTracking.ffmpeg_cmd(start_frame / fps,vid_name,loading_size,"test5.yuv"))
    read!("test5.yuv",temp_frames)

    temp_frames2[:]=convert(KnetArray{Float32,3},temp_frames)
    running_mean_i[:,:,1] = mean(temp_frames2,dims=3) ./ max_intensity

    for i=1:load_number
        @ffmpeg_env run(WhiskerTracking.ffmpeg_cmd((i*loading_size + start_frame)/ fps,vid_name,loading_size,"test5.yuv"))
        read!("test5.yuv",temp_frames)
        temp_frames2[:]=convert(KnetArray{Float32,3},temp_frames)

        x_k = mean(temp_frames2,dims=3) ./ max_intensity

        running_mean = running_mean_i .+ (x_k .- running_mean_i) ./ (i+1)

        running_std = running_std_i .+ (x_k .- running_mean_i).*(x_k - running_mean)

        running_mean_i = running_mean
        running_std_i = running_std
    end

    rm("test5.yuv")

    (convert(Array,running_mean), convert(Array,running_std))
end


text:
module FHist

export Hist1D, binedges, bincounts, bincenters, binerrors, nbins, integral
export sample, lookup, cumulative, normalize, restrict, rebin
export atomic_push!

export Hist2D, project, profile, transpose

export Hist3D

using StatsBase, RecipesBase, UnicodePlots, Statistics
import LinearAlgebra: normalize, normalize!
using Base.Threads: SpinLock

const _default_overflow = false

for (H, N) in ((:Hist1D, 1), (:Hist2D, 2), (:Hist3D, 3))
    @eval struct $H{T<:Real,E} <: AbstractHistogram{T,$N,E}
        hist::Histogram{T,$N,E}
        sumw2::Array{Float64, $N}
        hlock::SpinLock
        overflow::Bool
        function $H(h::Histogram{T,$N,E}, sw2 = copy(h.weights); overflow=_default_overflow) where {T,E}
            return new{T,E}(h, sw2, SpinLock(), overflow)
        end
    end
end

include("./utils.jl")
include("./hist1d.jl")
include("./hist2d.jl")
include("./hist3d.jl")
include("./displays.jl")
include("./arithmatics.jl")
include("./plot.jl")
end


text:
# IFRK

"""
    IFRK(u,Δt,plan_intfact,r₁;[rk::RKParams=RK31])

Construct an integrator to advance a system of the form

du/dt - Au = r₁(u,t)

The resulting integrator will advance the state `u` by one time step, `Δt`.

# Arguments

- `u` : example of state vector data
- `Δt` : time-step size
- `plan_intfact` : constructor to set up integrating factor operator for `A` that
              will act on type `u` (by left multiplication) and return same type as `u`
- `r₁` : operator acting on type `u` and `t` and returning `u`
"""
struct IFRK{NS,FH,FR1,TU}

  # time step size
  Δt :: Float64

  rk :: RKParams

  # Integrating factors
  H :: Vector{FH}

  r₁ :: FR1  # function of u and t, returns TU

  # scratch space
  qᵢ :: TU
  w :: Vector{TU}

end

function (::Type{IFRK})(u::TU,Δt::Float64,
                          plan_intfact::FI,rhs::FR1;
                          rk::RKParams{NS}=RK31) where {TU,FI,FR1,NS}

    # check for methods for r₁ and r₂
    if hasmethod(rhs,(TU,Float64))
        r₁ = rhs
    else
        error("No valid operator for r₁ supplied")
    end

    # scratch space
    qᵢ = deepcopy(u)
    w = [deepcopy(u) for i = 1:NS] # one extra for last step in tuple form

    dclist = diff([0;rk.c])

    # construct an array of operators for the integrating factor. Each
    # one can act on data of type `u` and return data of the same type.
    # e.g. we can call Hlist[1]*u to get the result.
    if TU <: Tuple
      (FI <: Tuple && length(plan_intfact) == length(u)) ||
                error("plan_intfact argument must be a tuple")
      Hlist = map(dc -> map((plan,ui) -> plan(dc*Δt,ui),plan_intfact,u),unique(dclist))
    else
      Hlist = map(dc -> plan_intfact(dc*Δt,u),unique(dclist))
    end

    H = [Hlist[i] for i in indexin(dclist,unique(dclist))]

    htype,_ = typeof(H).parameters

    # fuse the time step size into the coefficients for some cost savings
    rkdt = deepcopy(rk)
    rkdt.a .*= Δt
    rkdt.c .*= Δt

    ifrksys = IFRK{NS,htype,typeof(r₁),TU}(Δt,rkdt,H,r₁,qᵢ,w)

    #pre-compile
    ifrksys(0.0,u)

    return ifrksys
end

function Base.show(io::IO, scheme::IFRK{NS,FH,FR1,TU}) where {NS,FH,FR1,TU}
    println(io, "Order-$NS IF-RK integator with")
    println(io, "   State of type $TU")
    println(io, "   Time step size $(scheme.Δt)")
end

# Advance the IFRK solution by one time step
# This form works when u is a tuple of state vectors
function (scheme::IFRK{NS,FH,FR1,TU})(t::Float64,u::TU) where {NS,FH,FR1,TU <: Tuple}
  @get scheme (Δt,rk,H,r₁,qᵢ,w)

  # H[i] corresponds to H(i,i+1) = H((cᵢ - cᵢ₋₁)Δt)
  # Each of the coefficients includes the time step size

  i = 1
  tᵢ₊₁ = t
  for I in eachindex(u)
    qᵢ[I] .= u[I]
  end

  if NS > 1
    # first stage, i = 1
    w[i] = r₁(u,tᵢ₊₁)
    for I in eachindex(u)
      w[i][I] .*= rk.a[i,i]  # gᵢ

      # diffuse the scratch vectors
      qᵢ[I] .= H[i][I]*qᵢ[I] # qᵢ₊₁ = H(i,i+1)qᵢ
      w[i][I] .= H[i][I]*w[i][I] # H(i,i+1)gᵢ
      u[I] .= qᵢ[I] .+ w[i][I]
    end
    tᵢ₊₁ = t + rk.c[i]

    # stages 2 through NS-1
    for i = 2:NS-1

      w[i] = r₁(u,tᵢ₊₁)
      for I in eachindex(u)
        w[i-1][I] ./= rk.a[i-1,i-1] # w(i,i-1)

        w[i][I] .*= rk.a[i,i] # gᵢ

        # diffuse the scratch vectors and assemble
        qᵢ[I] .= H[i][I]*qᵢ[I] # qᵢ₊₁ = H(i,i+1)qᵢ
        for j = 1:i
          w[j][I] .= H[i][I]*w[j][I] # for j = i, this sets H(i,i+1)gᵢ
        end
        u[I] .= qᵢ[I] .+ w[i][I] # r₁
        for j = 1:i-1
          u[I] .+= rk.a[i,j]*w[j][I]
        end
      end
      tᵢ₊₁ = t + rk.c[i]

    end
    i = NS
    for I in eachindex(u)
      w[i-1][I] ./= rk.a[i-1,i-1] # w(i,i-1)
    end
  end

  # final stage (assembly)
  w[i] = r₁(u,tᵢ₊₁)
  for I in eachindex(u)
    w[i][I] .*= rk.a[i,i]
    u[I] .= qᵢ[I] .+ w[i][I] # r₁
    for j = 1:i-1
      u[I] .+= rk.a[i,j]*w[j][I] # r₁
    end
    u[I] .= H[i][I]*u[I]
  end
  t = t + rk.c[i]

  return t, u

end

# Advance the IFRK solution by one time step
function (scheme::IFRK{NS,FH,FR1,TU})(t::Float64,u::TU) where {NS,FH,FR1,TU}
  @get scheme (Δt,rk,H,r₁,qᵢ,w)

  # H[i] corresponds to H(i,i+1) = H((cᵢ - cᵢ₋₁)Δt)
  # Each of the coefficients includes the time step size

  i = 1
  tᵢ₊₁ = t
  qᵢ .= u

  if NS > 1
    # first stage, i = 1
    w[i] .= rk.a[i,i].*r₁(u,tᵢ₊₁) # gᵢ

    # diffuse the scratch vectors
    qᵢ .= H[i]*qᵢ # qᵢ₊₁ = H(i,i+1)qᵢ
    w[i] .= H[i]*w[i] # H(i,i+1)gᵢ
    u .= qᵢ .+ w[i]
    tᵢ₊₁ = t + rk.c[i]

    # stages 2 through NS-1
    for i = 2:NS-1
      w[i-1] ./= rk.a[i-1,i-1] # w(i,i-1)
      w[i] .= rk.a[i,i].*r₁(u,tᵢ₊₁) # gᵢ

      # diffuse the scratch vectors and assemble
      qᵢ .= H[i]*qᵢ # qᵢ₊₁ = H(i,i+1)qᵢ
      for j = 1:i
         w[j] .= H[i]*w[j] # for j = i, this sets H(i,i+1)gᵢ
      end
      u .= qᵢ .+ w[i] # r₁
      for j = 1:i-1
        u .+= rk.a[i,j]*w[j]
      end
      tᵢ₊₁ = t + rk.c[i]

    end
    i = NS
    w[i-1] ./= rk.a[i-1,i-1] # w(i,i-1)
  end

  # final stage (assembly)
  u .= qᵢ .+ rk.a[i,i].*r₁(u,tᵢ₊₁) # r₁
  for j = 1:i-1
    u .+= rk.a[i,j]*w[j] # r₁
  end
  u .= H[i]*u
  t = t + rk.c[i]

  return t, u

end


text:
using Test
using LinearAlgebra
using SparseArrays
using CompScienceMeshes

T = Float64
tol = eps(T) * 10^3
const third = one(T)/3

# test creation of a rectangle
m = meshrectangle(1.0,1.0,0.5);
@test numvertices(m) == 9
@test numcells(m) == 8

# test edge detection
edges = skeleton(m, 1)
@test numcells(edges) == 16

# test volume pair construction from a list of edges
vp = cellpairs(m,edges)
internal_edges = count(x->x>0, vp[2,:])
@test internal_edges == 8

# test the relative orientation routine
@test CompScienceMeshes.relorientation(index(1,2,3),index(1,2,3,4)) == +4
@test CompScienceMeshes.relorientation(index(2,3,4),index(1,2,3,4)) == -1
@test CompScienceMeshes.relorientation(index(1,3,2),index(1,2,3,4)) == -4
@test CompScienceMeshes.relorientation(index(2,4,3),index(1,2,3,4)) == +1

o, e1, e2, e3 = euclidianbasis(3)
cell = simplex(o, e1, e2)

A  = volume(cell)
@test abs(A - 0.5) < tol

mp = neighborhood(cell,[third, third]);
r = cartesian(mp)
@test norm(r - point(third, third, zero(T))) < tol

# repeat the test but now on a random cell
randpoint() = point(2*rand(T)-1, 2*rand(T)-1, 2*rand(T)-1)
v1 = randpoint()
v2 = randpoint()
v3 = randpoint()
cell = simplex(v1, v2, v3)

mp = neighborhood(cell,[third, third]);
r = cartesian(mp)

@test norm(third*(v1 + v2 + v3) - r) < tol
@test abs(dot(cross(r-v1, r-v2), r-v3)) < tol

# test that the normal is directed outwards
m = meshcircle(1.0, 2π/51)
#for i in 1:numcells(m)
for cl in cells(m)
    #p = chart(m, cells(m,i))
    v = vertices(m, cl)
    p = simplex(v)
    _c = (p.vertices[1] + p.vertices[2]) / 2
    @test dot(_c, p.normals[1]) > 0
end

## test cellpairs on non-oriented meshes
using CompScienceMeshes
using Test

p1 = point(1,0,0)
p2 = point(0,1,0)
p3 = point(-1,0,0)
p4 = point(0,-1,0)

i1 = index(1,2,3)
i2 = index(1,4,2)

m = Mesh([p1,p2,p3,p4], [i1,i2])
e = skeleton(m,1)

@test numcells(m) == 2
@test numcells(e) == 5

cps = cellpairs(m,e)

#Test unstructed mesh rectangle
r1 = meshrectangle(1.0,1.0,1/1)
# r2 = meshrectangle(1.0,1.0,1/1,structured=false)
r2 = readmesh(joinpath(@__DIR__,"assets/rectangle1.in"))

@test numcells(r1) == 2
@test numcells(r2) == 4
@test numvertices(r2) == 5
@test numcells(boundary(r2)) == 4
@test numcells(CompScienceMeshes.interior(r2)) == 4
@test numcells(boundary(r1)) == numcells(boundary(r2))

#Test mesh cuboid parts
deltaT = 1/1
box = meshcuboid(1.0,1.0,1.0,deltaT)
top = meshcuboid(1.0,1.0,1.0,deltaT,physical="TopPlate")
sides = meshcuboid(1.0,1.0,1.0,deltaT,physical="SidePlates")
bot = meshcuboid(1.0,1.0,1.0,deltaT,physical="BottomPlate")
obox = meshcuboid(1.0,1.0,1.0,deltaT,physical="OpenBox")

@test numcells(boundary(box)) == 0
@test numcells(boundary(top)) == numcells(boundary(bot))  == 4
@test numvertices(box) == 14
@test numcells(box) == 24
@test numcells(skeleton(box,1)) == 36
@test numcells(top) + numcells(bot) + numcells(sides) == numcells(box)
@test numcells(obox) == numcells(box) - numcells(top)

e1 = skeleton(top,1)
E1 = skeleton(box,1)

e1_cells = cells(e1)
E1_cells = cells(E1)

# cols = vec([6 16 17 13 18 19 20 21])
# Z11 = sparse(collect(1:8),cols,ones(8),numcells(e1),numcells(E1))
S11 = CompScienceMeshes.embedding(e1,E1)
@test size(S11) == (length(e1), length(E1))
vals = nonzeros(S11)
rows = rowvals(S11)
for col in axes(S11,2)
    for k in nzrange(S11,col)
        row = rows[k]
        val = vals[k]
        val == 0 && continue
        @test e1_cells[row] == E1_cells[col]
    end
end
# @test S11 == Z11